Universality Classes of Percolation Processes: Renormalization Group Approach

Abstract

Complex systems of classical physics in certain situations are characterized by intensive fluctuations of the quantities governing their dynamics. These include important phenomena such as (continuous) second-order phase transitions, fully developed turbulence, magnetic hydrodynamics, advective–diffusive processes, the kinetics of chemical reactions, percolation, and processes in financial markets. The theoretical goal of the study of such systems is to determine and predict the temporal and spatial dependence of statistical correlations of fluctuating variables. Modern methods of quantum field theory, originally developed in high-energy physics to describe the properties of elementary particles, allow for quantitative analysis of such correlations. Peculiarities of quantum field theory in solving two paradigmatic statistical problems related to percolation are reviewed, and new results on calculating representative universal parameters such as critical exponents that describe asymptotic behavior are presented.

1. Introduction

During the last decades, the field of fractal growth phenomena has exhibited continuous rapid growth, with the majority of important new results having been published in this period. The physics of far-from-equilibrium growth phenomena constitutes just one of the main fields in which fractal geometry is commonly applied. Much effort has been devoted to the study of experimental, theoretical, and numerical information on the geometrical characteristics of fractal growth. Nontrivial features and underlying mechanisms have been obtained by a thorough analysis of appropriate microscopic models [1,2].

Non-equilibrium systems are common in nature, and most observed phenomena may be characterized as being in some sort of a non-equilibrium state [3,4]. Typical examples encompass turbulence [5,6], pattern formations [7], the Earth’s atmosphere [8], living systems [9], etc. Other relevant examples can be found outside the domain of physics, for instance in chemistry, economics, and sociology.

At the fundamental level, non-equilibrium systems are inherently dynamic, and their theoretical analysis cannot be properly captured by standard equilibrium methods. There are many mechanisms by which systems can fall out of equilibrium. For instance, a continuous supply of energy into a hydrodynamic flow causes the transfer of energy down to microscopic scales, where it is dissipated due to viscosity effects. From a quantitative point of view, we expect equilibrium relations such as the well-known fluctuation–dissipation theorem [10,11] to be violated.

An important class of non-equilibrium systems is formed by reaction–diffusion models, which can be used to model a plethora of biological systems, growth processes, ecology, sociology, etc. These might be interpreted as stochastic processes for which it is possible to identify slowly varying collective variables. On the other hand, (fast) microscopic degrees of freedom are modeled by means of appropriate stochastic variables, and their description follows prescribed probabilistic principles [4,11]. In certain cases, emergent collective behavior in which many microscopic entities participate, might occur, allowing continuum approximation to be employed. The main advantage of such an approach lies in the possibility of using sophisticated field-theoretical methods borrowed from quantum and statistical field theory [2,11].

One of the paradigmatic models is directed percolation (DP), which constitutes arguably the simplest model for spreading processes [12,13]. It finds applications in systems which possess a preferred growth direction. Although directed percolation clusters are known to be self-sustaining, their global scaling behavior displays non-trivial features. Initially, DP was developed [14] to model the propagation of fluid through irregularly structured media (i.e., porous material) in the presence of an external force that brings about spatial anisotropy, e.g., gravity. From a biological point of view, DP can be understood as a simple model of the epidemic spreading of an infectious disease [13]. In this example, it is assumed that there is no immunization of sick individuals. During disease proliferation, infected individuals can spread the disease to their immediate neighbors. Although sick individuals can recover, they may also become infected repeatedly. It turns out [11] that such a simple model shows a critical behavior that occurs between a state without sick individuals (the absorbing state) and a fluctuating (active) state with non-zero sick individuals. Other systems with DP-like behavior have been identified [15] in such diverse research fields as particle physics [16,17], population dynamics [18,19], various reaction–diffusion schemes [20], and recently in describing the passage from a laminar flow to a turbulent one [5,21].

The main characteristics of the DP universality class are concisely encapsulated by famous Janssen–Grassberger hypothesis [22,23], according to which a model exhibits asymptotic DP behavior if the following requirements are satisfied: (i) the presence of a unique absorbing state; (ii) a scalar order parameter with solely positive values; (iii) local interactions; and (iv) no other special features, e.g. symmetries, quenched disorder, coupling to additional slow variables, etc. Although this DP hypothesis has not yet been mathematically proven, there is numerical support in its favor [15,24]. Systems fulfilling the aforementioned requirements are expected to show identical large-scale behavior near the critical point, i.e., the same asymptotic regimes and same values of critical indices [10,25]. Analogously to equilibrium critical phenomena, non-equilibrium phase transitions can be categorized into distinct universality classes [11,24] which contain microscopically different systems that nevertheless exhibit the same large-scale asymptotic behavior. Despite the great effort that has been invested in an analytical investigation of the DP process, the solution of its simplest possible version in $(1+1)$ dimensions is not yet known [15].

Typically, non-equilibrium systems suffer from the absence of exact solutions. This prompts the development of independent analytical, numerical, and computational approaches to acquire a deeper understanding of the general properties of non-equilibrium systems. In fact, the usual approach is to attempt to identify some special property of a system and then take advantage of it. In the case of percolation processes, such a distinguished feature is the existence of non-equilibrium continuous phase transition. In this context, fluctuations on all spatial and temporal scales are simultaneously significant, and universal critical behavior might be expected.

In this work, we are mainly interested in the universal behavior of two classes of percolation processes that can be formulated in terms of a suitable reaction–diffusion system. In particular, the directed percolation process consists of three reaction processes

$$2A\stackrel{K_{21}}{⟶}A,A\stackrel{K_{12}}{⟶}2A,A\stackrel{K_{10}}{⟶}⌀,$$

(1)

which can be interpreted [15] as coagulation, offspring production, and self-destruction, respectively. The corresponding K factors stand for the appropriate rate constants for a given reaction process. In Section 2, we provide a concise derivation of the action functional for such a reaction–diffusion system, which allows us to include effects such as randomness and spatial inhomogeneities in individual reaction processes. Arguably the most basic formulation of reaction systems consists in the specification of a corresponding master equation [26]. We briefly outline the basics of how a generic master equation can be reformulated by means of an appropriate Fock space, leading to an action functional that can be systematically treated by the powerful methods of quantum field theory.

Another well-known variation of the percolation process is isotropic percolation (IP), in the epidemic parlance referred to as the general epidemic process [27,28,29,30]. In the IP process, recovered individuals become immune, making them effectively inactive. In the given local neighborhood, the disease becomes locally depleted and leaves behind clusters of immune individuals, which can be satisfactorily described by percolation theory [12].

To make further progress, we concentrate on scaling regimes in the aforementioned percolation systems and on their corresponding continuum limits. Note that for both of these this limit might be derived in a controllable fashion [31]. We adopt a field-theoretic approach, the main advantage of which is the possibility of employing theoretical procedures such as graphical Feynman methods, perturbative renormalization groups, instantons, operator product expansion, etc. [25,32,33]. In particular, we employ a perturbative renormalization group to determine possible asymptotic regimes and calculate critical exponents perturbatively. Renormalization group methods, while far from the most important approach to non-equilibrium problems, often provide great insights into the origin of universality and the justification of scaling relations. Furthermore, and of more relevance for this paper, they provide a calculational framework which allows for systematic evaluation of critical exponents. In this regard, the well-known $\epsilon$-expansion plays an important role, where $\epsilon$ denotes the deviation from the upper critical dimension $d_c$ of an action functional.

In addition, we are interested in a perturbative approach based on the calculation of Feynman graphs in higher, i.e., more than two, multi-loop approximations. We note that unlike other analogous models, such as critical statics or dynamics, there are practically no related calculations for non-equilibrium systems. In the recent past, we managed to accomplish a complete three-loop analysis for the DP model [34,35,36]. Our next goal is to expand the established techniques and apply them in the context of IP, for which universal quantities have been determined up to two-loop precision in the perturbative calculation [13]. In the paper, we provide a short summary of three-loop calculations for the DP model and discuss specific issues related to two-loop approximation of the IP process that confirm previous results [13].

The rest of this article is organized as follows: in Section 2, we present the major steps of the Doi approach; in Section 3, we provide a brief overview of the Langevin equation and related Fokker–Planck equation; Section 4 contains a detailed construction of microscopic models leading to the percolation process; Section 5 is reserved for a general field-theoretic analysis; an RG analysis of directed percolation is performed in Section 6; and its isotropic variant in Section 7 is analyzed together with technical details of the calculation. Finally, we provide some concluding remarks in Section 8.

2. Doi Formalism

In many physical systems, the number of particles is not conserved. Despite the fundamental quantum nature of their microscopic interactions, in many practical circumstances we may adopt a classical description of their behavior. For a theoretical treatment of such classical systems with a non-conserved number of particles, a convenient way to begin is the Doi approach [37,38]. As the topic is well described in the literature, we provide only a short summary of the main points here.

Let us start with particles A that are located on a d-dimensional regular hypercubic lattice with lattice spacing a. Lattice sites are enumerated by integers ($i=1,2,\dots$). We assume that particles A continuously hop between adjacent lattice sites with the hopping factor $D/a^2$ (the factor $a^2$ is eliminated later on, when the continuum limit is taken). It is customary to assume that the reaction process can occur only for those particles located at the same site, and we denote the corresponding rate as $K_0$. The dynamics of such a stochastic problem can be summarized by specifying governing equations for probability distributions $P(t;\{n\left\}\right)$. Here, $\left\{n\right\}$ is a convenient shorthand notation for a given microstate $\{n_1,n_2,\dots \}$, in which there are $n_1$ particles at site 1, $n_2$ at site 2, etc. The corresponding equations, known in the literature as master equations [11,26], encode information about the gains and losses of outgoing and incoming probabilities for the given state $\left\{n\right\}$. In condensed notation, we can formulate them in the following form:

$$\frac{dP(t;\{n\left\}\right)}{dt}=\sum _{\left\{m\right\}}{R}_{m\to n}P(t;\left\{m\right\})-\sum _{\left\{m\right\}}{R}_{n\to m}P(t;\left\{n\right\}),$$

(2)

where $R_{m\to n}$ denotes the transition probability rate from microstate m to microstate n. Doi formalism [11,37,39] of such coupled differential equations (2) can be represented by a suitable Fock space. The main underlying features that allow this are the linearity and first-in-time derivative nature of (2). If we do not pose any site restrictions, we may introduce bosonic operators for each lattice site i with the following bosonic-like commutation relations:

$$\left[{\widehat{a}}_i,{\widehat{a}}_j^{†}\right]={\delta }_{ij},\left[{\widehat{a}}_i,{\widehat{a}}_j\right]=\left[{\widehat{a}}_i^{†},{\widehat{a}}_j^{†}\right]=0.$$

(3)

The ground (vacuum) state of the theory $|0\rangle$ is defined as ${\widehat{a}}_i|0\rangle =0$, where i is any site, and corresponds to the empty system. The bosonic relations (3) imply important commutation relations:

$${\widehat{a}}_i^n{\widehat{a}}_i^{†}=n{\widehat{a}}_i^{n-1}+{\widehat{a}}_i^{†}{\widehat{a}}_i^n,{\widehat{a}}_i{\widehat{a}}_i^{†n}=n{\widehat{a}}_i^{†(n-1)}+{\widehat{a}}_i^{†n}{\widehat{a}}_i.$$

(4)

For the given lattice configuration $\left\{n\right\}=\{n_1,n_2,\dots \}$, we introduce state the $|\{n\}\rangle$ as follows: $|\left\{n\right\}\rangle ={\widehat{a}}_1^{†n_1}{\widehat{a}}_2^{†n_2}\dots |0\rangle .$ In a straightforward manner, we can show the validity of two relations

$${\widehat{a}}_i|\left\{n\right\}\rangle =n_i|\{n_1,n_2,\dots ,n_i-1,\dots \}\rangle ,{\widehat{a}}_i^{†}|\left\{n\right\}\rangle =|\{n_1,n_2,\dots ,n_i+1,\dots \}\rangle .$$

(5)

Then, we have ${\widehat{a}}_i^{†}{\widehat{a}}_i|\left\{n\right\}\rangle =n_i|\left\{n\right\}\rangle$, which legitimizes the standard identification for the particle number operator located at a lattice site i:

$${\widehat{n}}_i={\widehat{a}}_i^{†}{\widehat{a}}_i.$$

(6)

The scalar product between states $|\{n\}\rangle$ and $|\{m\}\rangle$ differs from the usual quantum mechanical form

$$\langle \left\{n\right\}|\left\{m\right\}\rangle =\prod _{i=1}{n}_i!{\delta }_{n_i,m_i},$$

(7)

where ${\delta }_{i,j}$ is the Kronecker symbol and $n!$ denotes standard factorial quotient $1\times 2\times \dots n$.

All relevant information about the studied stochastic process is encoded in the probability distributions $P(t;\{n\left\}\right)$, which in the Doi mechanism [11,37,40] is transferred to the corresponding state vector $|\mathrm{\Phi }\rangle$. This state can be introduced in the following way:

$$|\mathrm{\Phi }\left(t\right)\rangle \equiv \sum _{\left\{n\right\}}P(t;\left\{n\right\})|\left\{n\right\}\rangle =\sum _{\left\{n\right\}}P(t;\left\{n\right\}){\widehat{a}}_1^{†n_1}{\widehat{a}}_2^{†n_2}\dots |0\rangle .$$

(8)

Here, the sum runs over all possible lattice occupations. The corresponding governing dynamic Equation (2) for the state vector $|\mathrm{\Phi }\rangle$ takes the Schrödinger-like form

$$\frac{d}{dt}|\mathrm{\Phi }\left(t\right)\rangle =-\widehat{H}|\mathrm{\Phi }\left(t\right)\rangle ,$$

(9)

where the pseudo-Hamiltonian $\widehat{H}$ (also called the Liouvillean, or Liouville operator [11,41]) appears; its precise form depends on the considered reaction scheme. We can formally integrate Equation (9) with the following result: $|\mathrm{\Phi }\left(t\right)\rangle =exp\left(-\widehat{H}t\right)|\mathrm{\Phi }\left(0\right)\rangle$. To finalize problem formulation, it is necessary to specify the initial state $|\mathrm{\Phi }(0)\rangle$. With the initial probability distribution $P(0;\{n\left\}\right)$, the initial state follows directly from the definition (8). For practical reasons [31], the most convenient choice is the Poisson distribution describing the system in an external number density field.

In addition, interacting A particles perform unbiased random walks between lattice sites. Here, we a brief derivation of the corresponding term ${\widehat{H}}_D$ of the Hamiltonian H. First, let us start with an elementary two-site system in which $n_1$ particles are located at site 1 and $n_2$ particles at site 2, with the unidirectional hopping process $1\to 2$ at the rate $D_0/a^2$. For such a process, it is straightforward to see that the master Equation (2) reads

$$\frac{dP(n_1,n_2)}{dt}=\frac{D_0}{a^2}(n_1+1)P(n_1+1,n_2-1)-\frac{D_0}{a^2}{n}_{1}P(n_1,n_2),$$

(10)

where the combinatorial factors $n_1+1$ and $n_1$ stem from the assumed jump independence of A particles. By multiplying Equation (10) by the expression $a_1^{†n_1}{a}_2^{†n_2}$ and performing sum ${\sum }_{\{n_1,n_2\}}$ over all permissible occupation numbers, we obtain

$$\frac{d|\mathrm{\Phi }\rangle }{dt}=\frac{D_0}{{\widehat{a}}^2}({\widehat{a}}_2^{†}{\widehat{a}}_1-{\widehat{a}}_1^{†}{\widehat{a}}_1)|\mathrm{\Phi }\rangle .$$

(11)

We can easily generalize this expression to the bidirectional case $1\leftrightarrow 2$:

$$\frac{d|\mathrm{\Phi }\rangle }{dt}=-\frac{D_0}{a^2}({\widehat{a}}_2^{†}-{\widehat{a}}_1^{†})({\widehat{a}}_2-{\widehat{a}}_1)|\mathrm{\Phi }\rangle .$$

(12)

From this, we can infer that the corresponding term responsible for diffusion in the Hamiltonian takes the form

$${\widehat{H}}_D=\frac{D_0}{a^2}({\widehat{a}}_2^{†}-{\widehat{a}}_1^{†})({\widehat{a}}_2-{\widehat{a}}_1).$$

(13)

As has been already pointed out, the DP model is based on three reactions (1). Let us start with the binary annihilation process $2A⟶A$ at a given site and obtain the desired term ${\widehat{H}}_R$ of the reaction Hamiltonian. Because any pair of particles might undergo mutual annihilation, we immediately obtain

$$\frac{dP\left(n\right)}{dt}=K_{21}(n+1)nP(n+1)-K_{21}n(n-1)P\left(n\right),$$

(14)

where, similar to the diffusion term, the combinatorial factors are taken into consideration.

We now recast (14) employing the relations (3) and (5) into the following form:

$$\frac{d|\varphi \rangle }{dt}=K_{21}({\widehat{a}}^{†}-{\widehat{a}}^{†2}){\widehat{a}}^2|\varphi \rangle ,$$

(15)

from which the reaction part of the Hamiltonian is as follows:

$${\widehat{H}}_{R1}=-K_{21}({\widehat{a}}^{†}{\widehat{a}}^2-{\widehat{a}}^{†2}{\widehat{a}}^2).$$

(16)

In a similar vein, we can work out the terms for the remaining two reactions (1), with the final result being

$${\widehat{H}}_R=-K_{21}({\widehat{a}}^{†}{\widehat{a}}^2-{\widehat{a}}^{†2}{\widehat{a}}^2)-K_{12}({\widehat{a}}^{†2}\widehat{a}-{\widehat{a}}^{†}\widehat{a})-K_{10}(\widehat{a}-{\widehat{a}}^{†}\widehat{a}).$$

(17)

Results (13) and (17) are easily generalized for the entire lattice system. The total Hamiltonian for the system that accounts for diffusion and percolation processes on the entire lattice is

$$\begin{array}{cc}\hfill {\widehat{H}}_D+{\widehat{H}}_R& =\frac{D_0}{a^2}\sum _{<ij>}({\widehat{a}}_i^{†}-{\widehat{a}}_j^{†})({\widehat{a}}_i-{\widehat{a}}_j)\hfill \\ \hfill & -\sum _i\left[K_{21}({\widehat{a}}_i^{†}{\widehat{a}}_i^2-{\widehat{a}}_i^{†2}{\widehat{a}}_i^2)+K_{12}({\widehat{a}}_i^{†2}{\widehat{a}}_i-{\widehat{a}}_i^{†}{\widehat{a}}_i)+K_{10}({\widehat{a}}_i-{\widehat{a}}_i^{†}{\widehat{a}}_i)\right],\hfill \end{array}$$

(18)

where the former sum runs over all neighboring lattice sites. Such non-Hermitean Hamiltonians are typical for out-of-equilibrium systems [11,20,24]. Non-hermiticity additionally implies that the reaction rates in the master Equation (2) cannot satisfy the condition of detailed balance. Therefore, the equilibrium state in the system cannot be characterized by the usual Gibbs distribution. There are additional differences between the Doi formalism and usual quantum mechanical expressions for calculating physical observables. It can be easily observed that it is not possible to express observables as bilinear products $\langle \mathrm{\Phi }|A|\mathrm{\Phi }\rangle$ in the Doi approach; based on the definition of the state vector (8), we would end up with expressions that are quadratic in probability $P(t;\{n\left\}\right)$.

Let us now derive the ensemble mean value of an observable A within the Doi formalism. Here, we further assume that A depends solely on particle numbers, i.e., $A=A\left(\right\{n\left\}\right)$, which is quite natural from the physics point of view. Typical examples of such quantities for chemically interacting systems include the average number of particles (concentration) $n⟷{\sum }_i{\widehat{a}}_i^{†}{\widehat{a}}_i$ and the pair correlation function $C(i,j)⟷{\widehat{a}}_i^{†}{\widehat{a}}_i{\widehat{a}}_j^{†}{\widehat{a}}_j$ between sites i and j.

The ensemble mean value of a quantity A is defined by the following formula:

$$\langle A\left(t\right)\rangle =\sum _{\left\{n\right\}}A\left(\left\{n\right\}\right)P(t;\left\{n\right\}).$$

(19)

For practical purposes, it is helpful to have a state $\langle \mathcal{P}|$ at our disposal, for which the following identity holds:

$$\langle A\left(t\right)\rangle =\sum _{\left\{n\right\}}P(t;\left\{n\right\})\langle \mathcal{P}|\widehat{A}\left(\left\{{\widehat{a}}^{†}\widehat{a}\right\}\right){\widehat{a}}_1^{†n_1}{\widehat{a}}_2^{†n_2}\dots |0\rangle =\langle \mathcal{P}|\widehat{A}|\mathrm{\Phi }\left(t\right)\rangle .$$

(20)

Replacing the formal solution of Equation (9) yields an alternative form of the right-hand side of Equation (20):

$$\langle A\left(t\right)\rangle =\langle \mathcal{P}|\widehat{A}exp(-\widehat{H}t)|\mathrm{\Phi }\left(0\right)\rangle .$$

(21)

In the mean value, the operator expression $\widehat{A}=\widehat{A}\left(\left\{{\widehat{a}}^{†}\widehat{a}\right\}\right)$ corresponds to the standard function $A\left(\right\{n\left\}\right)$ substituting $n_i\to {\widehat{a}}_i^{†}{\widehat{a}}_i$ at each lattice point i, which is enabled by Equation (6). From direct comparison of the relations (19) and (20), we anticipate that the projection state $\langle \mathcal{P}|$ ought to fulfill the following relation:

$$\langle \mathcal{P}|\widehat{A}\left(\left\{{\widehat{a}}^{†}\widehat{a}\right\}\right){\widehat{a}}_1^{†n_1}{\widehat{a}}_2^{†n_2}\dots |0\rangle =A\left(\left\{n\right\}\right)\langle \mathcal{P}|{\widehat{a}}_1^{†n_1}{\widehat{a}}_2^{†n_2}\dots |0\rangle =A\left(\left\{n\right\}\right).$$

(22)

This, in return, leads to two requirements

$$\langle \mathcal{P}|{\widehat{a}}_i^{†}=\langle \mathcal{P}|\langle \mathcal{P}|0\rangle =1,$$

(23)

that should hold for the projection state $\mathcal{P}$ and the operator ${\widehat{a}}^{†}$, respectively, and be satisfied at each site i. Using the commutation relations (3), we further deduce that the projection state of the following form

$$\langle \mathcal{P}|=\langle 0|exp(\sum _i{\widehat{a}}_i)$$

(24)

serves the desired purpose. Note that if the operator $\widehat{A}$ is written in the normal-ordered form [42] then it can be constructed solely in terms of the operators $a_i$ utilizing only the aforementioned properties (23) of the state $\langle \mathcal{P}|$, i.e.,

$$\langle \mathcal{P}|\widehat{A}\left(\left\{{\widehat{a}}^{†}\widehat{a}\right\}\right)=\langle \mathcal{P}|N\left[{\widehat{A}}_N\left(\{{\widehat{a}}^{†},\widehat{a}\}\right)\right]=\langle \mathcal{P}|{\widehat{A}}_N\left(\{1,\widehat{a}\}\right).$$

(25)

We define the normal form ${\widehat{A}}_N$ of the operator $\widehat{A}$ in a standard manner using the relation $\widehat{A}=N\left[{\widehat{A}}_N\left(\{{\widehat{a}}^{†},\widehat{a}\}\right)\right]$, in which N stands for the normal product (i.e., all creation operators are placed to the left of all annihilation operators). The $\widehat{a}$-dependent operator ${\widehat{A}}_N\left(\{1,\widehat{a}\}\right)$ equivalent to the average number of particles n is then ${\widehat{a}}_i$, whereas the correlation function $C(i,j)$ corresponds to the operator expression ${\widehat{a}}_i{\widehat{a}}_j+{\delta }_{ij}{\widehat{a}}_j$. Furthermore, it is convenient to commute the exponential term ${\mathrm{e}}^{{\sum }_i{\widehat{a}}_i}$ to the right in the relation (20). To this end, we use the formula

$${\mathrm{e}}^{{\widehat{a}}_i}{\widehat{a}}_i^{†}=({\widehat{a}}_i^{†}+1){\mathrm{e}}^{{\widehat{a}}_i},$$

(26)

which can be obtained from the commutation relations (4). Using (26), we recast Equation (21) into the form

$$\langle A\left(t\right)\rangle =〈0\left|{\widehat{A}}_N\left(\{1,\widehat{a}\}\right){\mathrm{e}}^{-\widehat{H}(\{{\widehat{a}}^{†}+1,\widehat{a}\}t}{\mathrm{e}}^{{\sum }_i{\widehat{a}}_i}\right|\mathrm{\Phi }\left(0\right)〉,$$

(27)

where we have shifted the creation operator ${\widehat{a}}_i^{†}\to {\widehat{a}}_i^{†}+1$ (Doi shift) in the expression for the operator $\widehat{A}$ (followed by a replacement ${\widehat{a}}_i^{†}\to 0$) and similarly in the original Hamiltonian $\widehat{H}\left(\{a^{†},a\}\right)$.

In many irreversible reaction processes [31], the Poisson distribution represents the most suitable choice for the initial state. Accordingly, at a given site i, the probability of having $n_i$ particles is

$$p\left(n_i\right)={\mathrm{e}}^{-n_0}\frac{n_0^{n_i}}{n_i!}.$$

(28)

Here, a parameter $n_0$ corresponds to the average number of particles; using the definition of the state vector (8), we obtain

$$|\mathrm{\Phi }\left(0\right)\rangle =\sum _{\left\{n_1\right\}}{\mathrm{e}}^{-n_0}\frac{n_0^{n_1}}{n_1!}{\widehat{a}}_1^{†n_1}\sum _{\left\{n_2\right\}}{\mathrm{e}}^{-n_0}\frac{n_0^{n_2}}{n_2!}{\widehat{a}}_2^{†n_2}\dots |0\rangle =\prod _i{\mathrm{e}}^{-n_0}{\mathrm{e}}^{n_0{\widehat{a}}_i^{†}}|0\rangle .$$

(29)

Observe that when the shift ${\widehat{a}}^{†}\to {\widehat{a}}^{†}+1$ is performed, the exponential part ${\mathrm{e}}^{-n_0}$ is effectively removed. The mean particle number of the Poisson distribution can be different for different sites; in this case, the effective initial Poisson state vector ${\mathrm{e}}^{{\sum }_i{\widehat{a}}_i}|\mathrm{\Phi }\left(0\right)\rangle ={\prod }_i{\mathrm{e}}^{n_{0i}{\widehat{a}}_i^{†}}|0\rangle$ describes the system in an external density field $n_{0i}$. This allows us to describe initial localized seeds of given particle numbers by taking the derivatives with respect to components $n_{0i}$ and putting $n_{0i}=0$ thereafter, in the same manner that response functions are calculated in statistical mechanics. This is necessary for the description of epidemic processes as well as for the percolation problem, in which the localized seed at the origin denotes the onset of a spreading process.

A most striking feature of the percolation model is the appearance of a phase transition between the absorbing phase and the active phase, which exhibits emergent universal behavior. To analyze such a transition quantitatively, it is advantageous to obtain a continuous limit of the model, which allows us to employ powerful renormalization group methods. We summarize the major points of such derivation here, and the topic has been covered thoroughly elsewhere [11,31,39,43].

The crucial step includes a calculation of matrix elements for the time evolution operator $exp(-\widehat{H}t)$. To this end, we use the well-known Trotter formula [43]. Accordingly, the exponential factor ${\mathrm{e}}^{-\widehat{H}t}$ can be rewritten through a limiting procedure:

$$exp(-\widehat{H}t)=\underset{\mathrm{\Delta }t\to 0}{lim}{\left(1-\widehat{H}\mathrm{\Delta }t\right)}^{t/\mathrm{\Delta }t}.$$

(30)

We take $N\mathrm{\Delta }t=t$, where N is the number of time slices. At the end of our derivation, we let $N\to \infty$. At each time instant, we insert a whole set of coherent states, which is in fact tantamount to a representation of operators ${\widehat{a}}_i^{†},{\widehat{a}}_i$ by complex numbers $\psi$. Corresponding coherent states can be introduced [43] in the following fashion:

$$|\psi \rangle =exp(-\frac{1}{2}{|\psi |}^2+\psi {\widehat{a}}^{†})|0\rangle .$$

(31)

It can be shown that such states form an eigenstate basis of the annihilation operator, i.e., they satisfy

$$\widehat{a}|\psi \rangle =\psi |\psi \rangle ,\langle \psi |{\widehat{a}}^{†}={\psi }^{\ast }\langle \psi |,$$

(32)

where the star denotes the corresponding complex-conjugate quantity. Further, a distinguished feature of coherent states is that they are not orthogonal and form an over-complete set of states. The overlap between different states is actually provided by the formula

$$\langle {\psi }_1|{\psi }_2\rangle =exp(-\frac{1}{2}|{\psi }_1{|}^2-\frac{1}{2}{|{\psi }_2|}^2+{\psi }_1^{\ast }{\psi }_2).$$

(33)

For a single site, the identity resolution has the following expression:

$$\begin{array}{c}\hfill 1=\sum _n\frac{1}{n!}|n\rangle \langle n|=\sum _{m,n}\frac{1}{n!}|n\rangle \langle m|{\delta }_{mn}=\int \frac{d{\psi }^{\ast }d\psi }{\pi }|\psi \rangle \langle \psi |,\end{array}$$

(34)

where the orthogonality relation

$${\delta }_{mn}=\frac{1}{\pi m!}\int d{\psi }^{\ast }\int d\psi exp(-|\psi {|}^2){\psi }^{\ast m}{\psi }^n$$

(35)

has been employed, ${\mathrm{e}}^{-{|\psi |}^2}$ is the weight function, and we have adopted the measure in the form $d{\psi }^{\ast }d\psi =\mathrm{Re}\psi \mathrm{Im}\psi$. Generalization of the relation (34) to the entire lattice system is straightforward:

$$1=\int \prod _i\frac{d{\psi }_{(i,j)}^{\ast }d{\psi }_{(i,j)}}{\pi }|{\left\{\psi \right\}}_j\rangle \langle {\{\psi \}}_j|,$$

(36)

where ${\left\{\psi \right\}}_j=({\psi }_{(1,j)},{\psi }_{(2,j)},\dots )$ is the set of all eigenvalues corresponding to the annihilation operators ${\widehat{a}}_i$ at each lattice point at the given instant of time $j\mathrm{\Delta }t(j=0,\dots ,N)$. Inserting the relation (36) in each appropriate time instant into the Trotter Formula (30), we arrive at

$$\begin{array}{c}exp(-\widehat{H}t)=\frac{1}{\mathcal{N}}\underset{\mathrm{\Delta }t\to 0}{lim}\int \left[d{\psi }^{\ast }\right]\left[d\psi \right]|{\left\{\psi \right\}}_N\rangle \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hfill \\ \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\times \prod _{j=1}^N\langle {\left\{\psi \right\}}_j|exp(-\widehat{H}\left(\{{\widehat{a}}^{†}+1,\widehat{a}\}\right)\mathrm{\Delta }t)|{\left\{\psi \right\}}_{j-1}\rangle \langle {\left\{\psi \right\}}_0|.\end{array}$$

(37)

Here, $\mathcal{N}$ stands for the normalization constant. We can abbreviate the functional measure as follows:

$$\left[d{\psi }^{\ast }\right]\left[d\psi \right]\equiv \prod _{i,j}d{\psi }_{(i,j)}^{\ast }d{\psi }_{(i,j)}.$$

(38)

Here in after, we assume that Hamiltonian is normal-ordered (this is explicitly the case for the pseudo-Hamiltonian (18) in the DP process). Using relations (32), we obtain

$$\begin{array}{c}\langle {\left\{\psi \right\}}_j|exp(-H\left(\{a^{†}+1,a\}\right)\mathrm{\Delta }t)|{\left\{\psi \right\}}_{j-1}\rangle \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hfill \\ \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=\langle {\left\{\psi \right\}}_j|{\left\{\psi \right\}}_{j-1}\rangle exp(-H({\left\{{\psi }^{\ast }\right\}}_j+1,{\left\{\psi \right\}}_{j-1})\mathrm{\Delta }t),\end{array}$$

(39)

in which the indicated transformation $H({\left\{{\psi }^{\ast }\right\}}_j+1,{\left\{\psi \right\}}_{j-1})$ corresponds to the replacement of operators by their corresponding eigenvalues ${\widehat{a}}_i\to {\psi }_i,{\widehat{a}}_i^{†}\to {\psi }_i^{\ast }$. The remaining expression in Equation (39) is a condensed form of the following expression:

$$\langle {\left\{\psi \right\}}_j|{\left\{\psi \right\}}_{j-1}\rangle =\prod _i\langle {\psi }_{(i,j)}|{\psi }_{(i,j-1)}\rangle$$

(40)

and we can rewrite it with the help of the overlap relation (33) as follows:

$$\langle {\psi }_{(i,j)}|{\psi }_{(i,j-1)}\rangle =exp(-{\psi }_{(i,j)}^{\ast }[{\psi }_{(i,j)}-{\psi }_{(i,j-1)}])exp(-\frac{1}{2}|{\psi }_{(i,j-1)}{|}^2+\frac{1}{2}{|{\psi }_{(i,j)}|}^2).$$

(41)

For the first order in $\mathrm{\Delta }t$, we can express the scalar product in Equation (37) in the following way:

$$\prod _j\langle {\psi }_{(i,j)}|{\psi }_{(i,j-1)}\rangle =exp(-\sum _j{\psi }_{(i,j)}^{\ast }\frac{d{\psi }_{(i,j)}}{dt}\mathrm{\Delta }t+O\left(\mathrm{\Delta }t\right))exp(-\frac{1}{2}|{\psi }_{(i,0)}{|}^2+\frac{1}{2}{|{\psi }_{(i,N)}|}^2).$$

(42)

Then, in the continuum time limit ($\mathrm{\Delta }t\to 0$) we obtain the evolution operator in the functional integral representation

$$\begin{array}{cc}\hfill exp(-\widehat{H}t)& =\int \frac{\left[d\psi \right]\left[d{\psi }^{\ast }\right]}{\mathcal{N}}|{\left\{\psi \right\}}_N\rangle \hfill \\ \hfill & \times exp(\frac{1}{2}|{\psi }_i{\left(t\right)|}^2-\frac{1}{2}{|{\psi }_i\left(0\right)|}^2-{\int }_0^{t}dt[{\psi }_i^{\ast }{\partial }_t{\psi }_i+H\left(\{{\psi }^{\ast }+1,\psi \}\right)])\langle {\left\{\psi \right\}}_0|.\hfill \end{array}$$

(43)

According to expression (27) for the average value of a quantity A, for expectation values such as the mean particle number n we must act on the operator (43) from the left by the vacuum state $\langle 0|$ and from the right by the effective initial state $|exp\left[{\sum }_i{\widehat{a}}_i\right]\mathrm{\Phi }\left(0\right)\rangle$. For the initial state provided by the Poisson distribution (29), the relation

$$\langle {\left\{\psi \right\}}_0|exp[\sum _i{\widehat{a}}_i]\mathrm{\Phi }\left(0\right)\rangle \prod _{i}exp(-\frac{1}{2}{|{\psi }_i\left(0\right)|}^2)=exp(\sum _i\left[n_0{\psi }_i^{\ast }\left(0\right)-|{\psi }_i{\left(0\right)}^2|\right])$$

(44)

is satisfied, where we have used relations (26), (31), and (32). On the other hand, we have

$$\begin{array}{c}\langle 0|{\left\{\psi \right\}}_N\rangle \prod _{i}exp(\frac{1}{2}{|{\psi }_i\left(t\right)|}^2)=\langle 0|exp\left\{\sum _i\left[-\frac{1}{2}{|{\psi }_i\left(t\right)|}^2+{\psi }_i\left(t\right){\widehat{a}}_i^{†}\right]\right\}|0\rangle \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hfill \\ \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} \times \prod _{i}exp\left(\frac{1}{2}{|{\psi }_i\left(t\right)|}^2\right)=\langle 0|\prod _{i}exp\left({\psi }_i\left(t\right){\widehat{a}}_i^{†}\right)|0\rangle =1.\end{array}$$

(45)

By combining relations (43)–(45) and inserting them into expression (21) for the expectation value of quantity A, we obtain the final result in the form of a functional integral:

$$\langle A\left(t\right)\rangle ={\mathcal{N}}^{-1}\int \left[d\psi \right]\left[d{\psi }^{\ast }\right]A_N\left(\{1,\psi \}\right)exp\left[\mathcal{S}\left(\{{\psi }^{\ast },\psi \}\right)\right],$$

(46)

where the functional $A_N\left(\{1,\psi \}\right)$ is determined from the normal form ${\widehat{A}}_N\left(\{1,\widehat{a}\}\right)$ of the operator $\widehat{A}$ by replacing the operators ${\widehat{a}}_i$ with their corresponding eigenvalues ${\psi }_i$, and the field-theoretic functional $\mathcal{S}$ takes a similar form to that employed previously (see e.g., [11]):

$$\mathcal{S}\left(\{{\psi }^{\ast },\psi \}\right)=\sum _i(n_0{\psi }_i^{\ast }\left(0\right)-{|{\psi }_i\left(0\right)|}^2-{\int }_0^{t}dt[{\psi }_i^{\ast }{\partial }_t{\psi }_i+H(\{{\psi }^{\ast }+1\},\left\{\psi \right\})]).$$

(47)

The condition

$$\mathcal{N}=\int \left[d\psi \right]\left[d{\psi }^{\ast }\right]exp\left[\mathcal{S}\left(\{{\psi }^{\ast },\psi \}\right)\right]$$

(48)

fixes the normalization constant $\mathcal{N}$. The initial term $-|{\psi }_i{\left(0\right)|}^2$ in (47) is a consequence of the inherent ambiguity of the interpolation definition of the functional integral [42,44]. “Correctly defining the path integral is equivalent to constructing a renormalized perturbation theory” [44]; thus, we briefly recall the construction of perturbation theory in the present approach. Together with the function $A_N$ and the initial probability density, the Hamiltonian may be extracted from the functional integral using the functional derivatives:

$$\begin{array}{c}exp\left[{\int }_0^{t}dtH(\{{\psi }^{\ast }+1\},\left\{\psi \right\})\right]=exp\left[{\int }_0^{t}dtH\left(\left\{\frac{\delta }{\delta A}+1\right\},\left\{\frac{\delta }{\delta A^{\ast }}\right\}\right)\right]\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hfill \\ \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\times exp\left[{\int }_0^{t}dt\left(A{\psi }^{\ast }+A^{\ast }\psi \right)\right]{|}_{A=A^{\ast }=0}.\end{array}$$

(49)

We are left with the Gaussian integral

$$G_0(A,A^{\ast })=\int \left[d{\psi }^{\ast }\right]\left[d\psi \right]exp\left[-{|\psi \left(0\right)|}^2+{\int }_0^{t}dt\left(-{\psi }^{\ast }{\partial }_t\psi +A{\psi }^{\ast }+A^{\ast }\psi \right)\right]$$

(50)

where, for economy of notation, the spatial subscript has been omitted. We now calculate the interpolation approximation for this integral (in which the subscript refers to consecutive time instants):

$$\begin{array}{c}{G}_0(A,A^{\ast })\approx \int \left[d{\psi }^{\ast }\right]\left[d\psi \right]exp\left[-|{\psi }_0{|}^2-\sum _{j=1}^N\left({\psi }_j^{\ast }{\psi }_j-{\psi }_j^{\ast }{\psi }_{j-1}\right)+\sum _{j=0}^N\left(A_j{\psi }_j^{\ast }+A_j^{\ast }{\psi }_j\right)\mathrm{\Delta }t\right]\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hfill \\ \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=exp\left[\sum _{j=0}^N{A}_j^{\ast }\mathrm{\Delta }t\sum _{l=0}^j{A}_l\mathrm{\Delta }t\right]\underset{\mathrm{\Delta }t\to 0}{\overset{}{\to }}exp{\int }_0^{t}du{\int }_0^{t}d{u}^{\prime }{A}^{\ast }\left(u\right)\theta (u-u^{\prime })A\left(u^{\prime }\right).\end{array}$$

(51)

Perturbation theory is based on given propagators and interaction functionals; therefore, the continuum limit should be taken upon calculation of the Gaussian integral, as in (51). Because the step function $\theta (t-t^{\prime })$ is the fundamental result of the differential equation ${\partial }_{t}f=0$, by reversing the order of the steps we arrive at the Gaussian functional integral in the form

$$\begin{array}{c}{G}_0(A,A^{\ast })=exp{\int }_0^{t}du{\int }_0^{t}d{u}^{\prime }{A}^{\ast }\left(u\right)\theta (u-u^{\prime })A\left(u^{\prime }\right)\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hfill \\ \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=\int \left[d{\psi }^{\ast }\right]\left[d\psi \right]exp{\int }_0^{t}dt\left(-{\psi }^{\ast }{\partial }_t\psi +A{\psi }^{\ast }+A^{\ast }\psi \right).\end{array}$$

(52)

The formal space of integration can then be defined to allow necessary operations with Gaussian integrals [42,44]. Therefore, in perturbation theory, the field-theoretic action functional $\mathcal{S}$ reads

$$\mathcal{S}\left(\{{\psi }^{\ast },\psi \}\right)=\sum _i(n_0{\psi }_i^{\ast }\left(0\right)-{\int }_0^{t}dt\left[{\psi }_i^{\ast }{\partial }_t{\psi }_i+H(\{{\psi }^{\ast }+1\},\left\{\psi \right\})\right]).$$

(53)

We choose the corresponding retarded propagator vanishing at coinciding time arguments. In the standard language of quantum field theories, this is equivalent to simply augmenting the definition of the time-ordered product at coinciding time arguments as the normal-ordered product [42].

The continuum spatial limit can then be obtained in a standard way using the replacement:

$$\sum _i\to \frac{\int d^{d}x}{a^d},{\psi }_i\to \psi (t,x)a^d,{\psi }_i^{\ast }\left(t\right)\to {\psi }^{†}(t,x),n_0\to n_0{a}^d,$$

(54)

which yields the field-theoretic action $\mathcal{S}$ for the directed percolation universality class

$$\begin{array}{cc}\hfill \mathcal{S}& =-{\int }_0^{t}dt\int d^{d}x[{\psi }^{†}{\partial }_t\psi -D_0{\psi }^{†}{\nabla }^2\psi +K_{21}{\psi }^{†}(1+{\psi }^{†}){\psi }^2-K_{12}{\psi }^{†}(1+{\psi }^{†})\psi \hfill \\ \hfill & -K_{10}{\psi }^{†}\psi ]+n_0\int d^{d}x{\psi }^{†}(0,x),\hfill \end{array}$$

(55)

which coincides with the lattice Hamiltonian (18). The derived action functional (55) can be employed in the evaluation of important physical quantities such as the average (27), which can be particularly represented in the following way:

$$\langle A\left(t\right)\rangle =\int \mathcal{D}{\psi }^{†}\int \mathcal{D}\psi A_N\left\{1,\psi \left(t\right)\right\}{\mathrm{e}}^{\mathcal{S}}.$$

(56)

3. Langevin Equation and Fokker–Planck Equation

A number of equations governing the time evolution of biological, chemical, physical, and social systems can be modeled by mean-field-like equations for averages of variables. However, they intrinsically correspond to random systems to a certain extent. In order to properly take fluctuations around the averages into account, a direct way is to insert a source of randomness straight into the mean-field equation. In this way, the quantities obtained as solutions from the mean-field equations become stochastic processes depending on spatial coordinates, i.e., stochastic fields.

The prototypical example of such an approach is the Langevin equation for a random walk. It serves to describe the location r of a test particle subject to random impulses, which can be accounted for by the noise term $\eta$ in the following way:

$$\frac{dr}{dt}=\eta .$$

Here, it is assumed that the random noise $\eta$ is of zero mean and decorrelated in time (white noise), i.e., its correlator takes the form

$$\langle {\eta }_i\left(t\right){\eta }_j\left(t^{\prime }\right)\rangle =D{\delta }_{ij}\delta (t-t^{\prime }).$$

In general, this typical physical formulation is mathematically inconsistent, as for the multiplicative noise (i.e., when the noise variable is multiplied by some function of the random location) it leads to inevitable ambiguities [45].

Let us deal with the mean distance between two positions on the trajectory of the random walk, an indication of the characteristic that the Brownian trajectory, although continuous, is nowhere differentiable [46]. The increments of the Brownian trajectory

$$W\left(t\right)-W\left(t_0\right)={\int }_{t_0}^{t}dt\eta \left(t\right)$$

produce a well-known random process, the so-called Wiener process, the conditional probability density of which is governed by Gaussian-like dependence:

$$p\left(W,t|W_0,t_0\right)=\frac{1}{{\left[4\pi (t-t_0)\right]}^{d/2}}{e}^{-{(W-W_0)}^2/2(t-t_0)}$$

where d is the dimension of space. In particular, the Wiener process is crucial for a rigorously consistent formulation of the Langevin equation via the stochastic differential equation (SDE).

According to Landau theory, close to the critical point, the dynamic behavior of the order parameter $\phi$ near equilibrium is governed by a kinetic equation (the time-dependent Ginzburg–Landau (TDGL) equation) with the following form:

$$\frac{\partial \phi }{\partial t}=-\mathrm{\Gamma }\left(-{\nabla }^2\phi +a\phi +\frac{\lambda }{6}{\phi }^3\right).$$

(57)

In linear response theory, the dynamics of fluctuations close to equilibrium [47] are generally defined using kinetic equations analogous to (57) with an additional random force variable added to the right-hand side.

In general, typical models of critical dynamics are formulated by means of nonlinear Langevin equations:

$$\frac{\partial \phi }{\partial t}=-\mathrm{\Gamma }\frac{\delta H}{\delta \phi }+f:=V\left(\phi \right)+f,$$

(58)

where H is the so-called effective equilibrium Hamiltonian. Gaussian distribution is postulated for a random source field f and the corellator $\langle ff\rangle$ is determined using the corresponding relation to the static equilibrium (fluctuation–dissipation theorem) [11,25]. An appropriate Gaussian distribution is required; for example, the dynamic model A with the non-conserved order parameter [10] is provided by the SDE, which simply follows from Equation (57) by the addition of a white-noise field to the right-hand side.

In reaction–diffusion systems, the most straightforward kinetic formulation of the dynamics of the average numbers of particles is based on the rate equation approach, which consists of postulating a deterministic differential equation for average numbers of particles in a homogeneous system; hence, it does not take into account effects such as boundary conditions, spatial inhomogeneities and randomness in the individual reaction events. Spatial dependence can be accounted for by a diffusion term, which then leads to models of diffusion-limited reactions (DLR) [11,19].

As an elementary example, let us consider the coagulation reaction $A+A\to A$. The corresponding diffusion-limited rate equation for the concentration $\phi =\phi (t,x)$ of the reacting compound A is

$$\frac{\partial \phi }{\partial t}=D{\nabla }^2\phi -k{\phi }^2,$$

where k is the appropriate rate constant and D is a diffusion parameter.

The most direct way to take into account the diverse effects of randomness is to include a random source and sink term in the rate equation:

$$\frac{\partial \phi }{\partial t}=D{\nabla }^2\phi -k{\phi }^2+f.$$

(59)

This is nothing else than a nonlinear Langevin equation for the stochastic field $\phi$. Physically, in the case of concentration we require the field $\phi$ to attain only non-negative values.

There is an essential distinction between the critical dynamics and the reaction–diffusion systems. In the former, deviations of the fluctuating order parameter from the (usually zero) mean can physically assign an arbitrary sign (or direction). In particular, deviations from the equilibrium value are always permissible. On the other hand, reaction–diffusion systems often exhibit a distinguished absorbing (inactive) state that does not permit order parameter fluctuations. Thus, should the system fall into the inactive state, it will remain there permanently. Above all, if the empty state is an inactive state of the reaction model, the random noise term should be introduced multiplied by a coefficient vanishing in the limit $\phi \to 0$ to ensure that the model cannot return from the inactive state purely by means of noise fluctuations. Then, arguably the simplest choice corresponds to

$$\frac{\partial \phi }{\partial t}=D{\nabla }^2\phi -k{\phi }^2+f\phi$$

(60)

as an alternative to Equation (59). This corresponds to an SDE with multiplicative noise.

Further, let us discuss the Langevin equation with multiplicative noise of the general form

$$\frac{\partial \phi }{\partial t}=V\left(\phi \right)+fb\left(\phi \right):=-K\phi +U\left(\phi \right)+fb\left(\phi \right),$$

(61)

where f is assumed to be a Gaussian noise variable with zero mean and the pair correlation function

$$\langle f(t^{\prime },x^{\prime })f(t,x)\rangle =\overline{D}(x^{\prime }-x)=\delta (t^{\prime }-t)D(x^{\prime }-x),$$

(62)

where the condensed notation $x=(t,x)$ has been employed. In Equation (61), the expression $b\left(\phi \right)$ stands for a functional of the stochastic field $\phi$. In what follows, we are interested in applying perturbation theory, for which it is useful to split b into a linear part K and nonlinear part $U\left(\phi \right)$. Functionals b and U are both time-local, i.e., they depend only on the instant of the SDE.

As has been already been pointed out, the Langevin equation with white-in-time noise f is mathematically inconsistent. It is possible to approach this issue differently by starting with the set of correlation functions that consist of a suitable delta sequence in time, i.e.,

$$\langle f(t^{\prime },x^{\prime })f(t,x)\rangle =\overline{D}(t^{\prime },x^{\prime };t,x)\underset{t\to t^{\prime }}{\overset{}{\to }}\delta (t^{\prime }-t)D(x^{\prime },x),$$

(63)

and moving to the white noise limit at a later phase. Mathematically, however, such treatment leads to the solution of the SDE (61) in the Stratonovich interpretation [46]. From a physics point of view, this is a common method to deal with the case of white noise. On the other hand, the Stratonovich formulation brings about a relatively complex technical treatment; therefore, the SDE is usually interpreted according to the Itô prescription [26,46].

Recall that the main reason for introducing the SDE with white noise is to evade explicitly solving the aforementioned limit (63). Accordingly, instead of applying the mathematically precarious (though physically quite admissible) Langevin equation, the stochastic problem provided by Equation (61) and the correlator (62) can be equivalently formulated in terms of the Fokker–Planck equation (FPE). This is an equation for both the conditional probability density $p\left(\phi ,t|{\phi }_0,t_0\right)$ and the probability density $p\left(\phi ,t\right)$ of the stochastic variable $\phi$. Recall that the Fokker–Planck equation and the master equation are particular examples of the universal (forward) Kolmogorov equation; therefore, the two cases discussed here are intimately connected. A smooth approach to prove this connection utilizes the method of Itô calculus; an appropriate proof can be found in the literature [46]. As a result, only the analogy between the quantities specifying the stochastic problem in both formulations is provided here. The main advantage of the FPE is that the equation itself is a completely well-defined partial differential (or functional differential for field variables) equation. The inherent ambiguity of the Langevin problem reveals itself in that the FPE takes different forms for different interpretations of the underlying SDE.

For simplicity, let us consider zero-dimensional field theory. The FPE for the conditional probability density $p\left(\phi ,t|{\phi }_0,t_0\right)$ in the case of the Itô interpretation takes the form

$$\begin{array}{cc}\hfill \frac{\partial }{\partial t}p\left(\phi ,t|{\phi }_0,t_0\right)& =-\frac{\partial }{\partial \phi }\left\{\left[-K\phi +U\left(\phi \right)\right]p\left(\phi ,t|{\phi }_0,t_0\right)\right\}\hfill \\ & +\frac{1}{2}\frac{{\partial }^2}{\partial {\phi }^2}\left[b\left(\phi \right)Db\left(\phi \right)p\left(\phi ,t|{\phi }_0,t_0\right)\right].\hfill \end{array}$$

(64)

On the other hand, if the SDE (61) is formulated in the Stratonovich sense, the FPE takes the form

$$\begin{array}{cc}\hfill \frac{\partial }{\partial t}p\left(\phi ,t|{\phi }_0,t_0\right)& =-\frac{\partial }{\partial \phi }\left\{\left[-K\phi +U\left(\phi \right)\right]p\left(\phi ,t|{\phi }_0,t_0\right)\right\}\hfill \\ & +\frac{1}{2}\frac{\partial }{\partial \phi }\left\{b\left(\phi \right)\frac{\partial }{\partial \phi }\left[Db\left(\phi \right)p\left(\phi ,t|{\phi }_0,t_0\right)\right]\right\}.\hfill \end{array}$$

(65)

The conditional probability density $p\left(\phi ,t|{\phi }_0,t_0\right)$ is the basic solution of the FPE (64) or (65), respectively, with the initial condition taking the following form:

$$p\left(\phi ,t_0|{\phi }_0,t_0\right)=\delta \left(\phi -{\phi }_0\right).$$

Contractions are not entirely visible when the random variable contains several components; for example, the Fokker–Planck equation in the Itô interpretation becomes

$$\begin{array}{cc}\hfill \frac{\partial }{\partial t}p\left(\phi ,t|{\phi }_0,t_0\right)& =-\frac{\partial }{\partial {\phi }_i}\left\{\left[-K_{ij}{\phi }_j+U_i\left(\phi \right)\right]p\left(\phi ,t|{\phi }_0,t_0\right)\right\}\hfill \end{array}$$

(66)

$$\begin{array}{cc}& \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\frac{1}{2}\frac{{\partial }^2}{\partial {\phi }_i\partial {\phi }_j}\left[b_{ik}\left(\phi \right)D_{kl}{b}_{jl}\left(\phi \right)p\left(\phi ,t|{\phi }_0,t_0\right)\right]\hfill \end{array}$$

(67)

for a multi-component stochastic field ${\phi }_i$.

The FPE can be viewed as the Schrödinger equation with imaginary time. Therefore, the Fokker–Planck solution together with the calculation of expectation values can be expressed in a manner similar to quantum field theory [48]. The construction with FPE as a starting point leads to the celebrated Martin–Siggia–Rose solution of the SDE [49] while at the same time avoiding the inherent lack of clarity in the SDE thanks to correction based on the selection of the FPE.

For completeness, we now discuss the FPE (64) corresponding to the Itô interpretation of the Langevin Equation (61). In close analogy with Dirac’s notation in quantum mechanics, we introduce the state vector $|p_t\rangle$ according to the following representation of the PDF:

$$p(\phi ,t)=\langle \phi |p_t\rangle ,$$

which is the solution of the FPE (64) with the initial condition $p(\phi ,0)=p_0\left(\phi \right)$. Further, in order to construct the evolution operator for the state vector, we establish coordinate and momentum operators, similarly to the case in quantum mechanics, using the following relations:

$$\widehat{\pi }f\left(\phi \right)=-\frac{\partial }{\partial \phi }f\left(\phi \right),\widehat{\phi }f\left(\phi \right)=\phi f\left(\phi \right),\left[\widehat{\phi },\widehat{\pi }\right]=1.$$

Note that the latter relation follows directly from the first two defining relations for $\widehat{\pi }$ and $\widehat{\phi }$. In these variables, the FPE for the PDF yields the evolution equation for the state vector in the form

$$\frac{\partial }{\partial t}|p_t\rangle =-\widehat{H}|p_t\rangle ,$$

where the “Hamilton”-operator for the FPE corresponding to the Itô formulation of the SDE, in accordance with (64), takes the form

$$\widehat{H}=\widehat{\pi }K\widehat{\phi }-\widehat{\pi }U\left(\widehat{\phi }\right)-\frac{1}{2}{\widehat{\pi }}^{2}b\left(\widehat{\phi }\right)Db\left(\widehat{\phi }\right).$$

(68)

Note that, in contrast to quantum mechanics, there is no issue with ordering uncertainty in the definition of the Hamilton operator. In this formulation, we may express the conditional PDF as the matrix element

$$p\left(\phi ,t|{\phi }_0,t_0\right)=〈\phi \left|e^{-\widehat{H}(t-t_0)}\right|{\phi }_0〉,$$

(69)

and the functional-integral construction can be derived in an identical way to that in quantum mechanics [44].

4. Functional Representation of General Reaction–Diffusion Scheme

In this section, we analyze the annihilation reaction $A+A\to ⌀$ (A being a particle, agent, etc., depending on a given model [11]) in a more general layout than usually presented [31]. To this end, we include random sinks and sources of the interacting particles to sustain a stationary state. In the majority of papers, this is performed by means of an additive shot noise in the Langevin equation of the given stochastic problem. Because the investigation here starts from the master equation approach, it is not quite justified for our purposes. To the best of our knowledge, however, we are not aware of any unique way to introduce random sources into the master equation corresponding to the random noise term of the mean-field (Langevin) description. We employ the simplest permissible choice, described in detail in the literature [26], which is tantamount to adding reactions $A\to X$ and $Y\to A$ to the entire scheme. Here, Y and X denote the particle baths of the source and the sink, respectively. In a homogeneous system, these reactions give rise to a master equation of the following form:

$$\frac{\mathrm{d}P(t,n)}{\mathrm{d}t}={\mu }_{+}V\left[P(t,n-1)-P(t,n)\right]+{\mu }_{-}\left[(n+1)P(t,n+1)-nP(t,n)\right]\dots$$

(70)

where $P(t,n)$ is the probability at time t of finding n reacting particles in the medium. The ellipsis stands for terms accounting for other processes, such as the annihilation reaction, diffusion, and advection processes. Further, the parameter ${\mu }_{-}$ (${\mu }_{+}$) denotes the reaction constant of the annihilation (creation) process. We choose the transition rate to be proportional to the number of particles n, which can be interpreted as a consequence of independent reactions $A\to X$. Moreover, this option retains the same empty state as an inactive state. The parameter V in Equation (70) is (at least for current purposes) the volume of the homogeneous system, and plays an important role in the transition to the continuum limit of the inhomogeneous system. Furthermore, note that the master Equation (70) leads to the kinetic rate equation

$$\frac{\mathrm{d}\langle n\rangle }{\mathrm{d}t}={\mu }_{+}V-{\mu }_{-}\langle n\rangle +\dots ,$$

(71)

where $\langle n\rangle$ denotes the average number of particles.

Let us recall here that the crucial conception of the Doi formalism [37,40] is to recast the family of master equations for probability distributions in a single equation for a state vector which contains all relevant physical information. The kinetic equation itself is formulated by means of the Liouville operator acting in the appropriate Fock space and produced by the group of master equations. Even though the main steps of this technique have been comprehensively discussed previously [11,50], the insertion of random sink and source terms in the master equation displays unique properties which should be elaborated thoroughly. To this end, let us briefly recall the fundamental quantities and equations of the Doi formalism. To keep the discussion concise, let us examine the probabilities $P(t,n)$ of finding n particles on a given lattice site at time t. Then, the spatial dependency can be defined by marking the particle number by the positions of the lattice and introducing products and sums over the lattice sites. The needed Fock space was introduced previously in Section 2, specifically in Equations (3)–(8). The family of master equations for a death–birth reaction can be recast in a single equation with no explicit dependence on the occupation number:

$$\frac{\mathrm{d}|\mathrm{\Phi }\rangle }{\mathrm{d}t}=-\widehat{H}({\widehat{a}}^{†},\widehat{a})|\mathrm{\Phi }\rangle .$$

(72)

Equation (70) leads to the expressions below in the corresponding “Hamiltonian”

$${\widehat{H}}_g({\widehat{a}}^{†},\widehat{a})={\mu }_{+}V\left(I-{\widehat{a}}^{†}\right)+{\mu }_{-}\left({\widehat{a}}^{†}-I\right)\widehat{a},$$

(73)

where I denotes the identity operator. The dynamic action $S_1$ corresponding to (72) and (73) reads

$${\mathcal{S}}_1({\psi }^{\ast },\psi )=\underset{0}{\overset{\infty }{\int }}dt{\psi }^{\ast }\left(t\right)\left[-{\partial }_t\psi \left(t\right)+{\mu }_{+}V-{\mu }_{-}\psi \left(t\right)\right]+\dots .$$

(74)

Only the expressions corresponding to the random source process and general time-derivative expression are explicitly mentioned here, while the ellipsis denotes those terms stemming from initial conditions and other included reaction processes.

Suppose next that the transition rates ${\mu }_{\pm }$ are arbitrary functions independent in time, with a probability distribution specified through the moments $\langle {\mu }_{\pm }^n\rangle =E_{\pm ,n}$. In order to maintain translation invariance with respect to the time variable, we further consider stationary stochastic processes defined on the entire time axis. Hereinafter, all time integrals in the functionals, and consequently in the Feynman diagrams arising in the perturbation theory, are taken over the whole time axis. Let us now generalize the overall analysis to the spatially inhomogeneous system. To this end, we replace a lattice subscript with the spatial argument, i.e., $\psi \left(t\right)\to {\psi }_i\left(t\right)$. In this way, the volume V turns into the volume element associated with each lattice site. To keep the notation simple, the time integral is replaced by the sum ${\int }_{-\infty }^{\infty }dt\to {\sum }_{\alpha }\mathrm{\Delta }t$ and we make a natural assumption that the transition rates at each time instant and lattice site ${\mu }_{\pm ,\alpha ,i}$ are independent random variables. As a direct consequence, the average over the distribution of random sources reduces to the calculation of the following expectation value

$$\prod _{\alpha ,i}\langle exp({\mu }_{+,\alpha ,i}V{\psi }_{\alpha ,i}^{\ast }\mathrm{\Delta }t-{\mu }_{-,\alpha ,i}{\psi }_{\alpha ,i}^{\ast }{\psi }_{\alpha ,i}\mathrm{\Delta }t)\rangle .$$

(75)

For each specific time instant and lattice site (noting that we assume the moments of ${\mu }_{\pm }$ to be equivalent for all $\alpha$ and i and omit the labels for brevity), this leads to the standard cumulant expansion:

$$\begin{array}{ccc}\hfill \langle {\mathrm{e}}^{\mu b\mathrm{\Delta }t}\rangle & =& 1+b\mathrm{\Delta }t{E}_1+\frac{1}{2}{E}_2{\left(b\mathrm{\Delta }t\right)}^2+\frac{1}{6}{E}_3{\left(b\mathrm{\Delta }t\right)}^3+\cdots \hfill \\ & =& exp(b\mathrm{\Delta }t{E}_1+\frac{E_2-E_1^2}{2}{\left(b\mathrm{\Delta }t\right)}^2+\frac{E_3-3E_1{E}_2+E_1^3}{6}{\left(b\mathrm{\Delta }t\right)}^3+\cdots ).\hfill \end{array}$$

(76)

Here, b denotes either $V{\psi }^{\ast }$ or $-{\psi }^{\ast }\psi$. Thus, for example, the average over the moment ${\mu }_{+}$ takes the form

$$\begin{array}{ccc}\hfill \prod _{\alpha ,i}\langle exp\left({\mu }_{+,\alpha ,i}V{\psi }_{\alpha ,i}^{\ast }\mathrm{\Delta }t\right)\rangle & =& exp\left(\sum _{\alpha }\sum _i\left[\mathrm{\Delta }t{E}_{+1}V{\psi }_{\alpha ,i}^{\ast }+\frac{1}{2}\left(E_{+2}-E_{+1}^2\right){\left(V{\psi }_{\alpha ,i}^{\ast }\mathrm{\Delta }t\right)}^2\right]\right)\hfill \\ & \times & exp\left(\sum _{\alpha }\sum _i\left[\frac{E_{+3}-3E_{+1}{E}_{+2}+E_{+1}^3}{6}{\left(V{\psi }_{\alpha ,i}^{\ast }\mathrm{\Delta }t\right)}^3+\cdots \right]\right).\hfill \end{array}$$

(77)

In the continuum limit, we replace the expression ${\psi }_{\alpha ,i}^{\ast }$ with the stochastic field variable ${\psi }^{†}(t,x)$, while in the limit $V\to 0$ the ratio ${\psi }_{\alpha ,i}/V$ leads to the variable $\psi (t,x)$, respectively. The sum over $\alpha$, together with a small time increment $\mathrm{\Delta }t$, brings about the time integral and the sum over lattice index i, while including the volume element V produces the spatial integral ${\sum }_{i}V\to \int dx$. Then, for the first expression in the exponential in Equation (77) we obtain

$$\sum _{\alpha }\sum _i\mathrm{\Delta }t{E}_{+1}V{\psi }_{\alpha ,i}^{\ast }\to E_{+1}\int dt\int dx{\psi }^{†}(t,x).$$

(78)

However, the continuum limit for cumulants of the second order and higher is not as self-evident. We consider the simplest permissible distribution for ${\mu }_{\pm }$, with the proviso that only the variance term has a finite limit when $V\to 0$ and $\mathrm{\Delta }t\to 0$. On the other hand, the contributions of higher-order cumulants disappear; for example we have

$$\begin{array}{c}\hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\left(E_{+2}-E_{+1}^2\right)V\mathrm{\Delta }t\to {\sigma }_{+},V\to 0,\mathrm{\Delta }t\to 0,\end{array}$$

(79)

$$\begin{array}{c}\hfill \left(E_{+3}-3E_{+1}{E}_{+2}+E_{+1}^3\right){\left(V\mathrm{\Delta }t\right)}^2\to 0,V\to 0,\mathrm{\Delta }t\to 0.\end{array}$$

(80)

Thus, the contribution of the average over ${\mu }_{+}$ to the action functional takes the form

$$S_{+}=\int dt\int dx\left\{E_{+1}{\psi }^{†}(t,x)+\frac{1}{2}{\sigma }_{+}{\left[{\psi }^{†}(t,x)\right]}^2\right\}.$$

(81)

For the average over ${\mu }_{-}$, an analogous argument yields

$${\mathcal{S}}_{-}=\int dt\int dx\left\{-E_{-1}{\psi }^{†}(t,x)\psi (t,x)+\frac{1}{2}{\sigma }_{-}{\left[{\psi }^{†}(t,x)\psi (t,x)\right]}^2\right\}.$$

(82)

These contributions to the action functional can, of course, be produced by appropriately chosen normal distributions of moments ${\mu }_{\pm }$.

The presented method of introducing random sources and sinks has the rather unpleasant property that it does not conserve the total number of particles in the system. To compare with the analysis of this problem using the Langevin approach, the random sources and sinks should be introduced in such a way that the particle number is preserved. The simplest way to approach this problem is to add an expression proportional to the particle number to the random source, i.e., to use the “reaction constant” of the form ${\mu }_{+}V+{\mu }_{1+}n$ instead of just ${\mu }_{+}V$ in the master equation. Then, the source terms on the right-hand side of Equation (70) take the form

$$\begin{array}{ccc}\hfill \frac{\mathrm{d}P(t,n)}{\mathrm{d}t}& =& {\mu }_{+}V\left[P(t,n-1)-P(t,n)\right]\hfill \\ & +& {\mu }_{1+}\left[(n-1)P(t,n-1)-nP(t,n)\right]\dots \hfill \end{array}$$

(83)

It can be observed that the new term in the master equation represents the branching process [26].

The added expression leads to the following contribution to the Hamilton operator:

$${\widehat{H}}_{g2}({\widehat{a}}^{†},\widehat{a})=-{\mu }_{1+}\left({\widehat{a}}^{†}-I\right){\widehat{a}}^{†}\widehat{a}.$$

(84)

By carrying out the steps described previously, we obtain the contribution to the dynamic action in the following form:

$${\mathcal{S}}_{1+}=\int dt\int dx\left\{E_{1+1}{\psi }^{†}\left({\psi }^{†}+1\right)\psi +\frac{1}{2}{\sigma }_{1+}{{\psi }^{†}}^2{\left({\psi }^{†}+1\right)}^2{\psi }^2\right\}.$$

(85)

Now it is straightforward to see that if we eliminate the plain source (i.e., let $E_{+1}={\sigma }_{+}=0$) and choose $E_{1+1}=E_{-1}$, the empty state remains an absorbing one and the mass-like term $\propto {\psi }^{†}\psi$ vanishes in the effective functional; we end up with the effective functional of the random sources and sinks:

$${\mathcal{S}}_{gc}=\int dt\int dx\{E_{1+1}{{\psi }^{†}}^2\psi +\frac{1}{2}{\sigma }_{-}{\left({\psi }^{†}\psi \right)}^2+\frac{1}{2}{\sigma }_{1+}{{\psi }^{†}}^2{\left(1+{\psi }^{†}\right)}^2{\psi }^2\},$$

(86)

which preserves the average number of particles.

The influence of the higher-order parts is considerably distinct in the two cases tractable by a dimensional analysis using the RG method. The term containing the derivative ${\partial }_t$ in the dynamic functional

$$\mathcal{S}=-\int dt\int dx{\psi }^{†}(t,x){\partial }_t\psi (t,x)+\dots$$

(87)

has to be dimensionless in order to provide interesting dynamics. Consequently, the total scaling dimension of the operator ${\psi }^{†}(t,x)\psi (t,x)$ equals the dimension of space, and as such is positive.

First, in the case where the scaling dimension of the field ${\psi }^{†}$ equals zero, i.e., $d_{{\psi }^{†}}=0$, the dimension of the field $\psi$ is positive (i.e., $d_{\psi }=d$) and the operator monomials in the second and third terms in the functional (86) attain the identical scaling dimension. Because they carry the quadratic term ${\psi }^2$, their scaling dimension is larger than that of the cubic term ${{\psi }^{†}}^2\psi$. As a result, the second and third terms in the action (86) are infrared-irrelevant and can be neglected in the ensuing asymptotic analysis.

Second, if both fields have positive scaling dimensions, then we have at minimum one additional field factor in the operator monomials in the second and third terms in (86) compared to the first term, making them irrelevant. Therefore, the corresponding infrared relevant effective functional of the random sinks and sources is reduced to the monomial

$${\mathcal{S}}_{gc}^{\prime }=\int dt\int dx{E}_{1+1}{{\psi }^{†}}^2\psi ,d_{{\psi }^{†}}=0\vee d_{{\psi }^{†}}>0,d_{\psi }>0.$$

(88)

Third, if the field $\psi$ has a vanishing scaling dimension, then the response field ${\psi }^{†}$ has a positive scaling dimension and the monomials with additional powers of response field ${\psi }^{†}$ are necessarily infrared-irrelevant. Thus, the initial point for the following renormalization analysis is the sink and source functional

$${\mathcal{S}}_{gc}^{″}=\int dt\int dx\left\{E_{1+1}{{\psi }^{†}}^2\psi +\frac{1}{2}\left({\sigma }_{-}+{\sigma }_{1+}\right){\left({\psi }^{†}\psi \right)}^2\right\},d_{\psi }=0.$$

(89)

Let us next apply canonical dimension analysis on the functional of the binary annihilation process $A+A\to ⌀$ with diffusion-limited spreading

$$\begin{array}{cc}\hfill {\mathcal{S}}_1& =-\int dt\int dx\{{\psi }^{†}{\partial }_t\psi -D_0{\psi }^{†}{\nabla }^2\psi +{\lambda }_0{D}_0\left[2{\psi }^{†}+{\left({\psi }^{†}\right)}^2\right]{\psi }^2\}+n_0\int dx{\psi }^{†}(0,x).\hfill \end{array}$$

(90)

In the first case with $d_{{\psi }^{†}}=0$, the terms with the quadratic factor ${\psi }^2$ in the action (90) have the same canonical dimensions. However, the monomial responsible for the sink–source term (88) is linear with respect to $\psi$, with a positive canonical dimension, unlike the nonlinear $\psi$ terms of the action (90). Thus, the infrared-relevant interaction part for more than two dimensions is provided by (88), and the corresponding action functional reads

$$\begin{array}{cc}\hfill {\mathcal{S}}_{IR1}& =-\int dt\int dx\left\{{\psi }^{†}{\partial }_t\psi -D_0{\psi }^{†}{\nabla }^2\psi -E_{1+1}{{\psi }^{†}}^2\psi \right\}+n_0\int dx{\psi }^{†}(0,x).\hfill \end{array}$$

(91)

It is easy to show that this dynamic action does not generate Feynman diagrams with closed loops of the density propagator. In turn, this means that fluctuation effects are effectively suppressed. On the other hand, the canonical dimension of the interaction part is negative, and can counterbalance the positive dimensions of the irrelevant interaction monomials. Consequently, the other monomials of the interaction part are actually dangerous irrelevant operators [51]. In this case, it is not possible to draw decisive conclusions about the IR-relevant action functional based on investigation of the canonical dimensions.

In the second case with positive dimensions $d_{\psi }$ and $d_{{\psi }^{†}}$, the quartic term in the functional (90) becomes IR-irrelevant. The remaining third-order terms alone do not generate loop diagrams. As a result, the effects of density fluctuations are generated only when both fields attain the same canonical dimension $d_{\psi }=d_{{\psi }^{†}}=d/2$. In this situation, the IR-relevant action functional reads

$$\begin{array}{cc}\hfill {\mathcal{S}}_{IR2}& =-\int dt\int dx\{{\psi }^{†}{\partial }_t\psi -D_0{\psi }^{†}{\nabla }^2\psi +2{\lambda }_0{D}_0{\psi }^{+}{\psi }^2-E_{1+1}{{\psi }^{†}}^2\psi \}+n_0\int dx{\psi }^{†}(0,x).\hfill \end{array}$$

(92)

In the functional, the canonical dimension of two interaction terms equals $(d/2)-2$ and disappears at the upper critical dimension $d_c=4$, at which remaining interaction monomials have positive dimensions and as such are clearly irrelevant. The functional (92) corresponds to the dynamic action of Reggeon field theory, known in the literature as the Gribov process [52]. Recently [53], the influence of random stirring for this model were studied using the Obukhov–Kraichnan compressible velocity field. For the incompressible Kraichnan model, even two-loop results are available [54].

In the third possible case with $d_{\psi }=0$, the quartic monomial in the functional (90) becomes irrelevant because the field ${\psi }^{†}$ attains a positive dimension. In addition, both monomials of the sink–source functional (89) are irrelevant, and we obtain the infrared relevant effective functional

$$\begin{array}{cc}\hfill {\mathcal{S}}_{IR3}& =-\int dt\int dx\{{\psi }^{†}{\partial }_t\psi -D_0{\psi }^{†}{\nabla }^2\psi +2{\lambda }_0{D}_0{\psi }^{†}{\psi }^2\}+n_0\int dx{\psi }^{†}(0,x).\hfill \end{array}$$

(93)

Reasoning analogous to that employed for action (91) demonstrates that the dimensional analysis with this choice of field dimensions does not guarantee an unambiguous resolution of the relevance of the interaction terms. Recall that the canonical dimensions of auxiliary quantities and the asymptotic behavior of Feynman diagrams are independent of the choice of the values of the field dimensions. Hence, the functional (92) with the explicit classification of irrelevant and relevant interaction parts characterizes the critical scaling behavior amenable to the renormalization analysis.

In short, if the sink and sources are selected in such a way that they preserve the average number of particles in the model, the anomalous scaling behavior in the model corresponds to that of Reggeon field theory.

A distinct case occurs if the primitive source term is included in the considerations. Then, it is possible that the model evolves towards an active phase in which the average number of particles is finite, instead of towards the absorbing (inactive) phase. In this case, the initial point of RG analysis is the effective functional which includes all of the aforementioned terms, i.e.,

$$\begin{array}{cc}\hfill \mathcal{S}& =\int dt\int dx\{-{\psi }^{†}{\partial }_t\psi +D_0{\psi }^{†}{\nabla }^2\psi -{\lambda }_0{D}_0\left[2{\psi }^{†}+{\left({\psi }^{†}\right)}^2\right]{\psi }^2+E_{+1}{\psi }^{†}\hfill \\ \hfill & +\frac{1}{2}{\sigma }_{+}{\left({\psi }^{†}\right)}^2+E_{1+1}{\psi }^{†}\left({\psi }^{†}+1\right)\psi +\frac{1}{2}{\sigma }_{1+}{{\psi }^{†}}^2{\left({\psi }^{†}+1\right)}^2{\psi }^2-E_{-1}{\psi }^{†}\psi \hfill \\ \hfill & +\frac{1}{2}{\sigma }_{-}{\left({\psi }^{†}\psi \right)}^2\}+n_0\int dx{\psi }^{†}(0,x).\hfill \end{array}$$

(94)

The stationarity equation corresponding to this functional for the field $\psi$ is

$${\partial }_t\psi -D_0{\nabla }^2\psi =-2{\lambda }_0{D}_0{\psi }^2+E_{+1}+E_{1+1}\psi -E_{-1}\psi .$$

(95)

Note that, as usual, we take ${\psi }^{†}=0$ as the stationary value for the response field ${\psi }^{†}$. On the other hand, the dynamic functional expanded around the stationary value becomes rather involved. To keep the model analysis manageable, let us proceed with a special case corresponding to the choice of $E_{1+1}=E_{-1}$. Then, the re-expansion of the action functional yields

$$\begin{array}{cc}\hfill \mathcal{S}& =\int dt\int dx\{-{\psi }^{†}{\partial }_t\psi +D_0{\psi }^{†}{\nabla }^2\psi -\sqrt{8}\sqrt{E_{+1}{\lambda }_0{D}_0}{\psi }^{†}\psi \hfill \\ \hfill & +\left(\frac{-E_{+1}}{2}+\frac{E_{-1}\sqrt{E_{+1}{\lambda }_0{D}_0}}{\sqrt{2}{\lambda }_0{D}_0}+\frac{E_{+1}{\sigma }_{1+}}{4{\lambda }_0{D}_0}+\frac{E_{+1}{\sigma }_{-}}{4{\lambda }_0{D}_0}+\frac{{\sigma }_{+}}{2}\right){{\psi }^{†}}^2\hfill \\ \hfill & +\frac{E_{+1}{\sigma }_{1+}{{\psi }^{†}}^3}{2{\lambda }_0{D}_0}+\frac{E_{+1}{\sigma }_{1+}{{\psi }^{†}}^4}{4{\lambda }_0{D}_0}+\frac{\sqrt{2}\sqrt{E_{+1}{\lambda }_0{D}_0}{\sigma }_{1+}{{\psi }^{†}}^3\psi }{{\lambda }_0{D}_0}\hfill \\ \hfill & +\left(E_{-1}-\sqrt{2}\sqrt{E_{+1}{\lambda }_0{D}_0}+\frac{\sqrt{E_{+1}{\lambda }_0{D}_0}{\sigma }_{1+}}{\sqrt{2}{\lambda }_0{D}_0}+\frac{\sqrt{E_{+1}{\lambda }_0{D}_0}{\sigma }_{-}}{\sqrt{2}{\lambda }_0{D}_0}\right){{\psi }^{†}}^2\psi \hfill \\ \hfill & +\frac{\sqrt{E_{+1}{\lambda }_0{D}_0}{\sigma }_{1+}{{\psi }^{†}}^4\psi }{\sqrt{2}{\lambda }_0{D}_0}-2{\lambda }_0{D}_0{\psi }^{†}{\psi }^2+\left(-{\lambda }_0{D}_0+\frac{{\sigma }_{1+}}{2}+\frac{{\sigma }_{-}}{2}\right){{\psi }^{†}}^2{\psi }^2\hfill \\ \hfill & +{\sigma }_{1+}{{\psi }^{†}}^3{\psi }^2+\frac{{\sigma }_{1+}{{\psi }^{†}}^4{\psi }^2}{2}\}+n_0\int dx{\psi }^{†}(0,x).\hfill \end{array}$$

(96)

Criticality corresponds to the limit $E_{+1}\to 0$. The variance ${\sigma }_{+}$ vanishes as well, as it is the expectation value of a non-negative random quantity ${\mu }_{+}$. In the critical region, we solely retain the leading $E_{+1}$ and ${\sigma }_{+}$, neglecting them in those terms where they are subleading. This renders the final form of functional somewhat simpler:

$$\begin{array}{cc}\hfill \mathcal{S}& =\int dt\int dx\{-{\psi }^{†}{\partial }_t\psi +D_0{\psi }^{†}{\nabla }^2\psi -\sqrt{8E_{+1}{\lambda }_0{D}_0}{\psi }^{†}\psi +\left(\frac{E_{-1}\sqrt{E_{+1}}}{\sqrt{2{\lambda }_0{D}_0}}+\frac{{\sigma }_{+}}{2}\right){{\psi }^{†}}^2\hfill \\ \hfill & +\frac{E_{+1}{\sigma }_{1+}{{\psi }^{†}}^3}{2{\lambda }_0{D}_0}+\frac{E_{+1}{\sigma }_{1+}{{\psi }^{†}}^4}{4{\lambda }_0{D}_0}+E_{-1}{{\psi }^{†}}^2\psi +\frac{\sqrt{2E_{+1}{\lambda }_0{D}_0}{\sigma }_{1+}{{\psi }^{†}}^3\psi }{{\lambda }_0{D}_0}\hfill \\ \hfill & +\frac{\sqrt{E_{+1}{\lambda }_0{D}_0}{\sigma }_{1+}{{\psi }^{†}}^4\psi }{\sqrt{2}{\lambda }_0{D}_0}-2{\lambda }_0{D}_0{\psi }^{†}{\psi }^2+\left(-{\lambda }_0{D}_0+\frac{{\sigma }_{1+}}{2}+\frac{{\sigma }_{-}}{2}\right){{\psi }^{†}}^2{\psi }^2\hfill \\ \hfill & +{\sigma }_{1+}{{\psi }^{†}}^3{\psi }^2+\frac{{\sigma }_{1+}{{\psi }^{†}}^4{\psi }^2}{2}\}+n_0\int dx{\psi }^{†}(0,x).\hfill \end{array}$$

(97)

Analysis of the scaling dimensions then leads to the following possibilities. In the interaction terms without the critical parameters $E_{+1}$ and ${\sigma }_{+}$, the earlier argumentation is valid; however, in parts endowed with powers of these parameters as coefficients, their positive scaling dimensions have to be reflected. The quadratic part of the dynamic functional (97) indicates that the canonical dimension of the parameter $E_{+1}$ is equal to four. Actually, the canonical dimension of ${\sigma }_{+}$ still needs to be regarded a free parameter of the model.

If we proceed in the same way as pointed out previously, we obtain the following dynamic functional for the infrared scaling limit. For the first case with a vanishing dimension of the response field ${\psi }^{†}$, the third and fourth powers of ${\psi }^{†}$ and independent of $\psi$ or the first power in $\psi$ can be neglected (because of the terms proportional to $\sqrt{E_{+1}}$ or $E_{+1}$) when comparing to the monomials $\propto {{\psi }^{†}}^2$ in the functional (97). Nonlinear monomials of the field $\psi$ are irrelevant versus the linear monomials because of the positive dimension of field $\psi$. In this case, the dynamic functional becomes

$$\begin{array}{cc}\hfill \mathcal{S}& =\int dt\int dx\{-{\psi }^{†}{\partial }_t\psi +D_0{\psi }^{†}{\nabla }^2\psi -2\sqrt{2}\sqrt{E_{+1}{\lambda }_0{D}_0}{\psi }^{†}\psi \hfill \\ \hfill & +\left(\frac{\sqrt{E_{+1}}{E}_{-1}}{\sqrt{2{\lambda }_0{D}_0}}+\frac{{\sigma }_{+}}{2}\right){{\psi }^{†}}^2+E_{-1}{{\psi }^{†}}^2\psi \}+n_0\int dx{\psi }^{†}(0,x).\hfill \end{array}$$

(98)

Again, the interaction part that remains after the dimensional analysis does not generate loop diagrams, even though the field $\psi$ displays nontrivial correlations; however, the canonical dimension of this expression attains negative value, which makes irrelevant terms dangerous and prevents any decisive statements about the relevance of the interaction terms.

In the second case with both fields attaining positive dimensions, i.e., $d_{{\psi }^{†}}>0$ and $d_{\psi }>0$, higher powers than the leading corrections to the quadratic functional of both fields are irrelevant. This reasoning yields the effective functional

$$\begin{array}{cc}\hfill \mathcal{S}& =\int dt\int dx\{-{\psi }^{†}{\partial }_t\psi +D_0{\psi }^{†}{\nabla }^2\psi -2\sqrt{2}\sqrt{E_{+1}{\lambda }_0{D}_0}{\psi }^{†}\psi \hfill \\ \hfill & +\left(\frac{\sqrt{E_{+1}}{E}_{-1}}{\sqrt{2{\lambda }_0{D}_0}}+\frac{{\sigma }_{+}}{2}\right){{\psi }^{†}}^2+E_{-1}{{\psi }^{†}}^2\psi -2{\lambda }_0{D}_0{\psi }^{†}{\psi }^2\}+n_0\int dx{\psi }^{†}(0,x).\hfill \end{array}$$

(99)

Unlike the situation discussed previously, here the cubic monomial $-2{\lambda }_0{D}_0{\psi }^{†}{\psi }^2$ creates loops itself because of the appearance of the correlation function of the field $\psi$. Hence, pairs of action functionals with nontrivial fluctuation contributions are permissible, which we now briefly discuss.

• The first effective action, for which the following inequality $d_{{\psi }^{†}}>d_{\psi }$ is satisfied. In order to preserve the correlation function of the field $\psi$ for the loop diagrams, the variance ${\sigma }_{+}$ has to possess a dimension lower than that of the parameter $\sqrt{E_{+1}}$. This leads to the dynamic functional

  $$\begin{array}{cc}\hfill \mathcal{S}& =\int dt\int dx\{-{\psi }^{†}{\partial }_t\psi +D_0{\psi }^{†}{\nabla }^2\psi -\sqrt{8E_{+1}{\lambda }_0{D}_0}{\psi }^{†}\psi \hfill \\ \hfill & +\frac{{\sigma }_{+}}{2}{{\psi }^{†}}^2-2{\lambda }_0{D}_0{\psi }^{†}{\psi }^2\}+n_0\int dx{\psi }^{†}(0,x)\hfill \end{array}$$

  (100)

  with the critical dimension depending on the canonical dimension of the parameter ${\sigma }_{+}$ similar to a description of tricritical scaling behavior [55]. However, the model is logarithmic at six dimensions, as in addition to the coefficient of the $\propto {{\psi }^{†}}^2$ the functional corresponds to the dynamic ${\phi }^3$ model. The upper critical dimension can be specified by the scaling behavior of the parameter ${\sigma }_{+}$ in the critical limit ${\sigma }_{+}\to 0$, $E_{+1}\to 0$.
• The second effective action, for which the fields have the same dimension, i.e., $d_{{\psi }^{†}}=d_{\psi }=d/2$. Both third-order monomials are infrared-relevant, and the dynamic functional basically becomes (99). In this situation, the canonical dimension of the parameter ${\sigma }_{+}$ is greater than the dimension of $\sqrt{E_{+1}}$. We omit ${\sigma }_{+}$ for simplicity. Therefore, the dynamic functional takes the form

  $$\begin{array}{cc}\hfill \mathcal{S}& =\int dt\int dx\{-{\psi }^{†}{\partial }_t\psi +D_0{\psi }^{†}{\nabla }^2\psi -\sqrt{8E_{+1}{\lambda }_0{D}_0}{\psi }^{†}\psi \hfill \\ \hfill & +\frac{\sqrt{E_{+1}}{E}_{-1}}{\sqrt{2{\lambda }_0{D}_0}}{{\psi }^{†}}^2+E_{-1}{{\psi }^{†}}^2\psi -2{\lambda }_0{D}_0{\psi }^{†}{\psi }^2\}+n_0\int dx{\psi }^{†}(0,x).\hfill \end{array}$$

  (101)
  Recall that this is the effective functional characterizing the Reggeon field theory (Gribov process) with a random source independent of the particle density. Assuming that the rate of change of the density owing to the random sink is proportional to the power of density is rather natural; on the other hand, assuming that the rate of change of the density as the random source is proportional to some power of density is not natural. Consequently, the effective functional (101) probably predicts a critical behavior of the Reggeon field theory (Gribov process) different from that analyzed in the literature.
• The third effective functional, for which the dimensions of fields satisfy conditions $d_{{\psi }^{†}}>d_{\psi }=0$, is the same as the first discussed functional (100).

The investigation of scaling dimensions reveals that we can in fact drop most of the limitations on the probability distribution of the transition rates of types (79) and (80). In fact, even if the higher-order cumulants are finite, the canonical dimensions of corresponding terms in the effective functional grow with the order of the cumulant, with the exception of the situation in which the transition rate does not depend on the particle density.

The dimensional analysis of the relevance of the interaction monomials described above seems to be rather formal. Above all, the arbitrariness of the canonical dimensions of the fields is a rather annoying property. For this reason, it is illuminating to repeat the investigation using the more traditional power counting method. Let us start with the case corresponding to the particle number-conserving sinks and sources. The effective functional is provided by Equation (94) at the criticality, i.e., with ${\sigma }_{+}=E_{+1}=0$ and $E_{-1}=E_{1+1}$,

$$\begin{array}{cc}\hfill \mathcal{S}& =\int dt\int dx\{-{\psi }^{†}{\partial }_t\psi +D_0{\psi }^{†}{\nabla }^2\psi -{\lambda }_0{D}_0\left[2{\psi }^{†}+{\left({\psi }^{†}\right)}^2\right]{\psi }^2\hfill \\ \hfill & +E_{1+1}{{\psi }^{†}}^2\psi +\frac{1}{2}{\sigma }_{1+}{{\psi }^{†}}^2{\left(1+{\psi }^{†}\right)}^2{\psi }^2+\frac{1}{2}{\sigma }_{-}{\left({\psi }^{†}\psi \right)}^2\}+n_0\int dx{\psi }^{†}(0,x).\hfill \end{array}$$

(102)

The divergence index [25,32] of a 1PI Feynman diagram reads

$$\delta =(2+d)L-2I,$$

(103)

where d is the spatial dimension, L corresponds to the number of loops, and I represents the number of internal lines. The standard conditions yield relations between the parameters L and I, the external field arguments $E_{\psi }$ and $E_{{\psi }^{†}}$, and the number of internal vertices $V_{ij}$. Here, the first index i represents the number of fields $\psi$ and the second index j stands for the number of response fields ${\psi }^{†}$ in the effective functional (102):

$$\begin{array}{cc}\hfill L& =1+I-\left(V_{12}+V_{21}+V_{22}+V_{23}+V_{24}\right),\hfill \\ \hfill I& =V_{12}+2V_{21}+2V_{22}+2V_{23}+2V_{24}-E_{\psi },\hfill \\ \hfill I& =2V_{12}+V_{21}+2V_{22}+3V_{23}+4V_{24}-E_{{\psi }^{†}}.\hfill \end{array}$$

(104)

From these equations we obtain, in particular,

$$E_{\psi }-E_{{\psi }^{†}}=-V_{12}+V_{21}-V_{23}-2V_{24}.$$

(105)

By simple algebraic manipulations, we further arrive at

$$\delta =2+d+(d-4)V_{21}+(d-2)V_{22}+d{V}_{23}+(d+2)V_{24}-(d-2)E_{\psi }-2E_{{\psi }^{†}}.$$

(106)

Multiplying the previous relation (105) by an unspecific parameter a and putting this together with Formula (106), we obtain the following relation for the divergence index:

$$\begin{array}{cc}\hfill \delta & =2+d-a{V}_{12}+(d+a-4)V_{21}+(d-2)V_{22}+(d-a)V_{23}\hfill \\ \hfill & +(d-2a+2)V_{24}-(d+a-2)E_{\psi }-(2-a)E_{{\psi }^{†}}.\hfill \end{array}$$

(107)

Here, we can immediately observe that the value of the coefficient a is actually equivalent to prescribing the values of the canonical dimensions of the fields. Because the index $\delta$ does not depend on the parameter a, it might be concluded that the choice of the canonical dimension of the fields does not have any effect on the IR or UV behavior of a 1PI Feynman diagram. We do not have any mass term in the model, and as a result the index $\delta$ with at least one negative coefficient of a vertex number indicates a possible infrared divergence, whereas a positive $\delta$ represents the standard superficial UV divergence.

For the following analysis, we use the notation $v_{nm}{\psi }^n{{\psi }^{†}}^m$ for a general interaction term and specify any particular interaction term simply by using the coefficient factor $v_{nm}$. It follows from the relation (106) that interactions $v_{23}$ and $v_{24}$ always contribute positively to $\delta$ irrespective of the space dimension; hence, they are certainly infrared-irrelevant. From the perspective of the dimensional value, the interaction $v_{21}$ generates infrared divergences below four dimensions, while the interaction $v_{22}$ is irrelevant above two dimensions. Because there is no available method to deal with them directly, we must turn to the usual relation between the ultraviolet and infrared divergences in logarithmic theory [25] and extrapolate its results below the critical dimension using $\epsilon$ expansion. Hence, if the interaction $v_{21}$ is present, the best that can be done is to perform a renormalization analysis around the critical dimension $d_c=4$, in which the interaction $v_{22}$ is irrelevant as well; in this case, we are left with the Reggeon field theory (Gribov process). From the viewpoint of the preceding investigation of the canonical dimension, Formula (106) corresponds to the dynamic functional (92) at the logarithmic dimension $d=4$ with the field dimensions $d_{\psi }=d-2$, $d_{{\psi }^{†}}=2$ and dimensionless coupling constant of the vertex $v_{12}$.

However, using the auxiliary relation (105) we can remove any one vertex number except $V_{22}$, which then induces corrections in the coefficients of the other vertex numbers (apart from $V_{22}$) and an alternative classification of the relevance of a given interaction, although the index of any Feynman diagram is not affected by the change, e.g., choice $a=4-d$ in (107) to remove interaction $V_{21}$ leads to

$$\delta =d+2+(d-4)V_{12}+(d-2)V_{22}+(2d-4)V_{23}+(3d-6)V_{24}-2E_{\psi }-(d-2)E_{{\psi }^{†}}.$$

(108)

Here, $v_{22}$, $v_{23}$, and $v_{24}$ are all irrelevant at spatial dimensions satisfying $d>2$, while $v_{12}$ is relevant at $d<4$. Therefore, if $V_{12}>0$, the sole critical regime accessible to renormalization analysis is that of the Reggeon field theory (Gribov process). From the previous examination, it is evident that removing a vertex number $V_{ij}$ from the formula for the index $\delta$ is equivalent to setting the canonical dimension of this vertex to zero and correspondingly fixing the canonical dimensions of the fields. The goal of power counting is to eliminate those Feynman diagrams of a perturbation calculation of a Green function with a canonical dimension (i.e., the divergence index) greater than that of the leading-order term. The index $\delta$ does not depend on the selection of the canonical dimension of the fields. Indeed, the sole tunable parameter it might depend on is the space dimension. If the canonical dimension of a vertex is zero, this essentially entails that we are seeking the critical behavior of the model at a space dimension in which we do not have any constraints whatsoever on the number of these vertices in the marginal and relevant diagrams.

Because our main interest is to analyze the influence of random sinks and sources on the binary annihilation process, we shall preserve the $v_{21}$ interaction in the effective infrared model in any case; accordingly, we set the canonical dimension of this vertex to zero for the purposes of the physical model. This means that $d_{\psi }=2$ and $d_{{\psi }^{†}}=d-2$, while that the critical regime of the model corresponds to the behavior of the Reggeon field theory (Gribov process), which is obviously part of the universality class distinct from the pure annihilation process [11].

From the viewpoint of the classification of terms in power counting, the two dynamic functionals (91) and (93) in the earlier investigation of scaling-dimensions correspond to the removal of the terms $E_{{\psi }^{†}}$ and $E_{\psi }$, respectively, from the formula for the UV index $\delta$, where we obtain

$$\begin{array}{cc}\hfill \delta & =d+2-2V_{12}+(d-2)(V_{21}+V_{22}+V_{23}+V_{24})-d{E}_{\psi }\hfill \end{array}$$

(109)

$$\begin{array}{cc}& \hspace{1em}\hspace{1em} =d+2-2V_{21}+(d-2)(V_{12}+V_{22})+(2d-2)V_{23}+(3d-2)V_{24}-d{E}_{{\psi }^{†}}.\hfill \end{array}$$

(110)

From the former relation (109) we can infer that all other interactions except $v_{12}$ are irrelevant above two spatial dimensions, yielding the dynamic functional (91). However, these are all dangerous irrelevant operators [51], as the negative dimension of the interaction $v_{12}$ can cause the dimension of a Feynman diagram with irrelevant operators to be negative. Because all dangerous operators have the same dimension, they cannot be organized according by their degree of “dangerousness”. Formula (110) gives rise to an analogous situation with respect to interactions other than $v_{21}$, which possess a negative dimension there. Nonetheless, the dimensions of irrelevant operators here are different, and this can provide a basis for classification of certain operators as more dangerous than others.

We remark that all of the derived Formulas (106)–(110) for the divergence index are identical, and demonstrate only particular choices of the ambiguous field dimensions. We observe that in a multi-charge case it is somewhat problematic to obtain a self-consistent conclusion about the importance of various interactions based only on the investigation of scaling dimensions of fields and the operators corresponding to the interactions. In particular, it seems rather misleading to set any of the field dimensions equal to zero at the beginning. Assuming that both field dimensions are non-zero renders them equal to each other, leading to

$$\begin{array}{cc}\hfill \delta & =2+d+\left(\frac{d}{2}-2\right)(V_{12}+V_{21})+(d-2)V_{22}+\left(\frac{3d}{2}-2\right)V_{23}\hfill \\ \hfill & +(2d-2)V_{24}-\frac{d}{2}{E}_{\psi }-\frac{d}{2}{E}_{{\psi }^{†}}.\hfill \end{array}$$

(111)

In this formula, each coefficient of vertex numbers is a function of the space dimension d. We can immediately conclude that the critical scaling behavior amenable to renormalization group analysis is that of the Reggeon field theory (Gribov process) [13,22].

Next, let us consider the functional (96) for the active state. In this case, the index $\delta$ does not reveal the whole truth due to the presence of the mass-like parameter $E_{+1}$; nevertheless, it demonstrates the contribution of the wave vector and frequency integral. While additional types of interactions are generated, the set of constraints limiting their numbers remains the same. Because the mass parameter is not involved in the power counting of ultraviolet divergences in any other way now that new interactions have appeared, it is advantageous to analyze the UV and IR behavior independently. Therefore, power counting is performed in the same way as earlier. Labeling new interactions with a tilde symbol, it is straightforward to obtain

$$\begin{array}{cc}\hfill L& =1+I-\left({\tilde{V}}_{02}+{\tilde{V}}_{03}+{\tilde{V}}_{04}+V_{12}+{\tilde{V}}_{13}+{\tilde{V}}_{14}+V_{21}+V_{22}+V_{23}+V_{24}\right),\hfill \end{array}$$

(112)

$$\begin{array}{cc}\hfill I& =V_{12}+{\tilde{V}}_{13}+{\tilde{V}}_{14}+2V_{21}+2V_{22}+2V_{23}+2V_{24}-E_{\psi },\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hfill \end{array}$$

(113)

$$\begin{array}{cc}\hfill I& =2{\tilde{V}}_{02}+3{\tilde{V}}_{03}+4{\tilde{V}}_{04}+2V_{12}+3{\tilde{V}}_{13}+4{\tilde{V}}_{14}+V_{21}+2V_{22}+3V_{23}+4V_{24}-E_{{\psi }^{†}}.\hfill \end{array}$$

(114)

We can immediately derive the relation

$$E_{\psi }-E_{{\psi }^{†}}=-2{\tilde{V}}_{02}-3{\tilde{V}}_{03}-4{\tilde{V}}_{04}-V_{12}-2{\tilde{V}}_{13}-3{\tilde{V}}_{14}+V_{21}-V_{23}-2V_{24}.$$

(115)

Substituting relations (112) and (113) in the Formula (103) directly follows

$$\begin{array}{cc}\hfill \delta & =\left(d+2\right)(1-{\tilde{V}}_{02}-{\tilde{V}}_{03}-{\tilde{V}}_{04})-2(V_{12}+V_{13}+V_{14})\hfill \\ \hfill & +\left(d-2\right)(V_{21}+V_{22}+V_{23}+V_{24})-d{E}_{\psi }.\hfill \end{array}$$

(116)

Adding the relation (115) multiplied by the parameter a to Formula (116) yields

$$\begin{array}{cc}\hfill \delta & =d+2+\left(2a-d-2\right){\tilde{V}}_{02}+\left(3a-d-2\right){\tilde{V}}_{03}+\left(4a-d-2\right){\tilde{V}}_{04}\hfill \\ \hfill & +\left(a-2\right)V_{12}+\left(2a-2\right){\tilde{V}}_{13}+\left(3a-2\right){\tilde{V}}_{14}+\left(d-a-2\right)V_{21}\hfill \\ \hfill & +\left(d-2\right)V_{22}+\left(d+a-2\right)V_{23}+\left(d+2a-2\right)V_{24}-a{E}_{\psi }-\left(d-a\right)E_{{\psi }^{†}}.\hfill \end{array}$$

(117)

Here, through a proper selection of the value of the coefficient a we can try to find a formulation suitable for the classification of the interaction vertices. Further, the interaction vertex with the negative dimension ${\psi }^{†}{\psi }^{†}$ might be absorbed in the pair correlation function of the field $\psi$ contained in the perturbative components of the graphical representation. When calculating the ultraviolet index, the field dimensions are completely determined by the requirement that the coupling constant of the term $\propto {\psi }^{†}{\psi }^{†}$ be a dimensionless quantity. There is no ambiguity in the formula for the UV index, and by marking the number of correlation functions in a 1PI diagram $\tilde{I}$ we obtain

$$\begin{array}{cc}\hfill L& =1+I+\tilde{I}-\left({\tilde{V}}_{03}+{\tilde{V}}_{04}+V_{12}+{\tilde{V}}_{13}+{\tilde{V}}_{14}+V_{21}+V_{22}+V_{23}+V_{24}\right),\hfill \\ \hfill I+2\tilde{I}& =V_{12}+{\tilde{V}}_{13}+{\tilde{V}}_{14}+2V_{21}+2V_{22}+2V_{23}+2V_{24}-E_{\psi },\hfill \\ \hfill I& =3{\tilde{V}}_{03}+4{\tilde{V}}_{04}+2V_{12}+3{\tilde{V}}_{13}+4{\tilde{V}}_{14}+V_{21}+2V_{22}+3V_{23}+4V_{24}-E_{{\psi }^{†}}.\hfill \end{array}$$

(118)

The formula for the divergence index $\delta$ follows directly from these relations:

$$\begin{array}{cc}\hfill \delta & =(2+d)L-4\tilde{I}-2I=d+2+\left(\frac{d}{2}+1\right)({\tilde{V}}_{03}-E_{{\psi }^{†}})+\left(d+2\right){\tilde{V}}_{04}\hfill \\ \hfill & +\left(\frac{d}{2}-1\right)(V_{12}-E_{\psi })+d{\tilde{V}}_{13}+\left(\frac{3d}{2}+1\right){\tilde{V}}_{14}+\left(\frac{d}{2}-3\right)V_{21}\hfill \\ \hfill & +(d-2)V_{22}+\left(\frac{3d}{2}-1\right)V_{23}+2d{V}_{24}.\hfill \end{array}$$

(119)

The identical relation follows from Equation (117) when the parameter a is selected in a such way that the coefficient in front of ${\tilde{V}}_{02}$ vanishes. From the relation (119), we can infer that when the vertex $v_{21}$ is marginal, everything else is irrelevant and the model corresponding to the field-theoretic functional (100) is recovered.

However, this is not the complete story, as in the case of infrared behavior we are mainly interested in the limit of the vanishing mass parameter. Because the coupling constants of new interactions depend on the parameter $E_{+1}$, they necessarily affect the limit of small $E_{+1}$. The dependence of a Feynman diagram on the mass parameter $E_{+1}$ in the critical region (${\mathit{k}}_i\to 0$, ${\omega }_i\to 0$, $E_{+1}\to 0$, ${\omega }_i=O\left(\sqrt{E_{+1}}\right)$, ${\mathit{k}}_i^2=O\left(\sqrt{E_{+1}}\right)$) can be easily predicted using a suitable scaling of the variables in the loop integrals if these integrals possess negative dimensions. In such a case, the upper limits of a renormalized diagram in terms of integral might be sent to infinity. The outcome of the scaling of integration variables is a power of the parameter $E_{+1}$ multiplied by some function of the reduced frequencies and wave vectors ${\omega }_i/\sqrt{E_{+1}}$, ${\mathit{k}}_i{E}_{+1}^{-1/4}$. The wave number integrals defining this function are IR- and UV-finite for all values of the reduced frequencies and wave vectors; thus, they possess a finite limit when these vanish. Therefore, any infrared singular behavior is signaled by a negative total power of the mass parameter $E_{+1}$, which takes into consideration both the scaling of integration variables (this provides the index $\delta$) and the additional powers of $E_{+1}$ at the vertex factors of the Feynman diagram. As a result, we can employ the following infrared divergence index brought about by the effective functional (96) (here, the parameter ${\sigma }_{+}$ is assumed to be subleading for simplicity):

$$\begin{array}{cc}\hfill {\delta }_{\mathrm{IR}}& =\delta +2{\tilde{V}}_{02}+4{\tilde{V}}_{03}+4{\tilde{V}}_{04}+2{\tilde{V}}_{13}+2{\tilde{V}}_{14}=\left(d+2\right)(1+{\tilde{V}}_{04})\hfill \\ \hfill & +\left(\frac{d}{2}+2\right){\tilde{V}}_{03}+\left(\frac{d}{2}-2\right)(V_{12}+V_{21})+d{\tilde{V}}_{13}+\frac{3d}{2}{\tilde{V}}_{14}\hfill \\ \hfill & +(d-2)V_{22}+\left(\frac{3d}{2}-2\right)V_{23}+(2d-2)V_{24}-\frac{d}{2}(E_{\psi }-E_{{\psi }^{†}}).\hfill \end{array}$$

(120)

to guide the classification of the vertex terms. The same form can be derived from the generic expression (117) by choosing the parameter a in such a way that the coefficient connected with the term ${\tilde{V}}_{02}$ in index ${\delta }_{\mathrm{IR}}$ vanishes. To keep the infrared behavior tractable, all coefficients of the vertex numbers in (120) have to be non-negative. The lowest spatial dimension fulfilling this criterion is the upper critical dimension of the dynamic functional. From relation (120), it is easy to see that $d_c=4$ for action (96) and that the dynamic functional is indeed provided by Equation (101).

5. Field-Theoretic Analysis of General Quantum Field Model

A typical stochastic model can be formulated in terms of the SDE

$${\partial }_t\phi \left(x\right)=U\left(\phi \left(x\right)\right)+f\left(x\right),〈f\left(x\right)f\left(x^{\prime }\right)〉=D(x,x^{\prime }),$$

(121)

where $\phi =\phi \left(x\right)$ is a fluctuating field, the term $U=U\left(\phi \right(x\left)\right)$ is a prescribed t-local functional (although explicit x dependence in U is permissible [25], and we do consider such models in this paper) which contains relevant physical information about the studied system, $f=f\left(x\right)$ is a random force with a Gaussian distribution with zero mean, and $D=D(x,x^{\prime })$ is a pair correlation function reflecting the underlying physics. It can be shown [11,25,56,57,58] that the statistical properties of the solution of the stochastic problem (121) are identical to those of the field-theoretic model described by the following (effective) action functional:

$$\begin{array}{cc}\hfill \mathcal{S}& =\frac{1}{2}\int \mathrm{d}x\int \mathrm{d}{x}^{\prime }{\phi }^{\prime }\left(x\right)D(x,x^{\prime }){\phi }^{\prime }\left(x^{\prime }\right)+\int \mathrm{d}x{\phi }^{\prime }\left(x\right)\left[-{\partial }_t\phi \left(x\right)+U\left(\phi \left(x\right)\right)\right]\hfill \\ \hfill & +\int \mathrm{d}x\left[A_{\phi }\left(x\right)\phi \left(x\right)+A_{{\phi }^{\prime }}\left(x\right){\phi }^{\prime }\left(x\right)\right],\hfill \end{array}$$

(122)

where ${\phi }^{\prime }={\phi }^{\prime }\left(x\right)$ is the Martin-Siggia-Rose response field [49] and $A_{\phi },A_{{\phi }^{\prime }}$ are source terms [11,25].

A crucial parameter in the RG approach is played by the upper critical dimension $d_c$. Above $d_c$, the mean-field theory, which neglects fluctuations of the order parameter, yields accurate values for the universal quantities [25,32]. Below $d_c$, strong fluctuations of the order parameter significantly affect the behavior of the critical system. Thus, the mean-field approach is not appropriate for this regime, and more sophisticated methods are needed. Exactly at the critical dimension $d_c$, the RG procedure leads to mean-field expressions amended with logarithmic corrections [25,32].

In statistical physics, we aim for multiplicative renormalization, which makes sense of otherwise divergent quantities while yielding additional information on large-scale behavior [25]. First, the original action $\mathcal{S}\left[\mathrm{\Phi }\right]$ together with its parameters $e_0$ are declared to be unrenormalized quantities; parameters $e_0$ (the whole set of parameters) are the bare parameters, and are assumed to be some functions of the new renormalized parameters e, while a new renormalized action is considered to be the functional ${\mathcal{S}}_R\left[\mathrm{\Phi }\right]=\mathcal{S}\left[\mathrm{\Phi }{Z}_{\mathrm{\Phi }}\right]$ with certain renormalization constants of fields $Z_{\mathrm{\Phi }}$. In unrenormalized full n-point Green functions $G_n=\langle \mathrm{\Phi }\dots \mathrm{\Phi }\rangle$, the averaging procedure $\langle \dots \rangle$ is carried out with the weight functional $exp\mathcal{S}\left[\mathrm{\Phi }\right]$, while in renormalized functions, $G_{nR}$ i used with weight $exp{\mathcal{S}}_R\left[\mathrm{\Phi }\right]$, respectively. From the multiplicative relationship between the action functionals $\mathcal{S}\left[\mathrm{\Phi }\right]$ and ${\mathcal{S}}_R\left[\mathrm{\Phi }\right]$, we infer a relation between the corresponding Green functions $G_{nR}=Z_{\mathrm{\Phi }}^{-n}{G}_n$, where by definition $G_n=G_n(e_0,\epsilon \dots )$ (the ellipsis denotes other relevant arguments such as wave numbers or coordinates) and by convention the quantities $G_{nR}$ and $Z_{\mathrm{\Phi }}$ are expressed using the renormalized parameters e. It is assumed that the correspondence $e_0\leftrightarrow e$ in perturbation theory is bijective, i.e., one-to-one. Second, the typical problem to be dealt with is the presence of various ultraviolet (UV) divergences in perturbative expansions. In practice, it is customary to use dimensional regularization with the minimal subtraction (MS) scheme [25,32]. In dimensional regularization, the space dimension d becomes a complex parameter and deviation from the upper critical dimension $\epsilon$ serves as a formally small parameter of the theory. It is possible to regularize UV divergent graphs with this method, and MS represents a useful method for calculation of various universal quantities [25,32]. For instance, we have $d_c=4$ for the DP process, whereas $d_c=2$ for binary annihilation process $A+A\to ⌀$. The crucial quantitative predictions of the RG methods for universal quantities such as critical exponents are asymptotic series. In practice, a resummation procedure is needed to obtain reliable numerical estimations [32].

Finally, we have to employ the renormalization procedure itself, which consists of the reorganization of perturbative expressions such that the limit $\epsilon \to 0$ is well-defined and finite [25,32]. The actual UV divergences are hidden in various renormalization constants [11,25,32]. With this choice, all UV-divergences (poles in parameter $\epsilon$) are absorbed in the parameters $e_0(e,\epsilon )$ and RG constants $Z_{\mathrm{\Phi }}^n(e,\epsilon )$. On the other hand, UV divergences are absent in connected renormalized Green functions $W_{nR}(e,\epsilon )$. Note that the condition of elimination of divergences does not precisely fix the functions $e_0(e,\epsilon )$ and RG constants $Z_{\mathrm{\Phi }}(e,\epsilon )$. A considerable arbitrariness remains, which allows us to introduce an additional dimensional parameter, namely, the renormalization mass (scale setting parameter) $\mu$ in these functions (and via them into functions $W_{nR}$):

$$W_{nR}(e,\mu ,\epsilon ,\dots )=Z_{\mathrm{\Phi }}^{-n}(e,\mu ,\epsilon )W_n(e_0(e,\mu ,\epsilon ),\epsilon ,\dots ).$$

(123)

A change in the renormalization mass $\mu$ at fixed bare parameters $e_0$ causes a change of $e,Z_{\mathrm{\Phi }}$ and $W_R$ for unchanged $W_n(e_0,\epsilon ,\dots )$. We label ${\tilde{\mathcal{D}}}_{\mu }$ for the differential operator $\mu {\partial }_{\mu }$ at fixed bare parameters $e_0$ and operate on both sides of the equation $Z_{\mathrm{\Phi }}{W}_{nR}=W_n$ with it. Thus, we obtain the fundamental RG equation

$$\left[\mu {\partial }_{\mu }+\sum _e{\tilde{\mathcal{D}}}_{\mu }e{\partial }_e+N{\gamma }_{\phi }\right]W_{nR}(e,\mu ,\epsilon ,\dots )=0,$$

(124)

where ${\gamma }_{\phi }\equiv {\tilde{\mathcal{D}}}_{\mu }ln{Z}_{\phi }$ and the operator ${\tilde{\mathcal{D}}}_{\mu }$ is expressed in the parameters $\mu$ and e. The coefficients ${\tilde{\mathcal{D}}}_{\mu }e$ and anomalous dimension ${\gamma }_{\phi }$ are called the RG functions, and are expressed in terms of different renormalization constants Z. Charges (coupling constants) g are those quantities e of the model on which the renormalization constants $Z=Z\left(g\right)$ might in principle depend. The logarithmic derivatives of charges in (124) are beta functions: ${\beta }_g={\tilde{\mathcal{D}}}_{\mu }g$. All the RG functions are UV-finite [25], i.e., they do not possess poles in $\epsilon$. This is a direct consequence of the functions $W_{nR}$ being UV-finite in Equation (124). The crucial importance of Equation (124) lies in the fact that it can be mathematically analyzed, for instance, using the method of characteristics. Its solutions can exhibit scaling behavior that can be directly related to various experimentally measurable universal quantities.

Further simplification is possible, and consists of an introduction of functional $\mathrm{\Gamma }\left(\alpha \right)$ for one-particle irreducible (1PI) Green functions [25,32]. These are graphs that remain connected even when a single internal line is cut [32,42]. It can be easily demonstrated [25] that the functional for 1PI diagrams is related to the functional W of connected Feynman graphs by suitable Legendre transformation:

$$\mathrm{\Gamma }\left(\alpha \right)=W\left(A\right)-\alpha A,\alpha \left(x\right)=\frac{\delta W\left(A\right)}{\delta A\left(x\right)}.$$

(125)

The arguments for the functional $\mathrm{\Gamma }$ are a pair of fields $\alpha =\{{\alpha }_{\psi },{\alpha }_{\tilde{\psi }}\}$, while $A=A\left(\alpha \right)$ is determined by the second relation in Equation (125). To obtain a given 1PI Green function ${\mathrm{\Gamma }}_{(n,m)}$, we only need to make appropriate derivatives with respect to source fields $\alpha$:

$${\mathrm{\Gamma }}_{(n,m)}\equiv \frac{{\delta }^{n+m}}{\delta {\alpha }_{\psi }^n\delta {\alpha }_{\tilde{\psi }}^m}\mathrm{\Gamma }\left(\alpha \right){|}_{\alpha =0}.$$

(126)

By virtue of causality, all 1PI Green functions of the form ${\mathrm{\Gamma }}_{(0,m)}$ identically vanish.

To simplify the notation, it is convenient to rename variables $\alpha$ according to the following rule $\alpha \to \phi$ and obtain a crucial equation [25,32], which relates the bare action functional $\mathcal{S}$ and generating functional $\mathrm{\Gamma }$:

$$\mathrm{\Gamma }\left(\phi \right)=\mathcal{S}\left(\phi \right)+{\mathrm{\Gamma }}_{\mathrm{loop}}\left(\phi \right).$$

(127)

The initial step of the perturbation procedure in any field-theoretic model consists of a suitable representation of Green functions in terms of Feynman graphs, the relevance of which is governed by the charges of the model. Feynman diagrams are a convenient graphical tool for the classification of different algebraic structures that are constructed from elementary perturbation structures (propagators and interaction vertices).

6. Directed Percolation

There are two well-established ways for a field-theoretic construction of the DP process. One of them is based on microscopic considerations and may be regarded as more fundamental, whereas the other follows from phenomenological reasoning. For our purposes here, the knowledge of the most important aspects of the ensuing model is sufficient, as its construction has been covered in detail in the past [13,22,24] and crucial details have been elaborated in Section 2. Thus, we consider the DP model in the epidemical terminology [13]; its guiding principles can be stated as follows:

• The only slow variable of the theory is the density of sick individuals $\psi (t,x)$.
• Time evolution can be characterized by means of a Markov process.
• The individual (e.g., person, animal, plant) is infected, and the corresponding rate is a function of the local density of sick individuals $\psi$. After a certain time amount of time, sick individuals spontaneously recover.
• The state with no sick individuals ($\psi =0$) is considered absorbing and cannot be left. Effectively, the spread of percolation terminates.
• The disease proliferates by means of a diffusion process.
• The effect of fast degrees of freedom can be captured by a proper choice of the stochastic variable that reflects the absorbing condition.

The basic governing quantity for the DP process is the density of active objects $\psi =\psi \left(x\right)$ [11,24]. Universal behavior (critical properties) in the vicinity of the phase transition is then completely described by the following stochastic differential equation:

$${\partial }_t\psi =D_0{\nabla }^2\psi +D_0{\tau }_0\psi -\frac{D_0{u}_0}{2}{\psi }^2+\sqrt{\psi }\zeta ,$$

(128)

where $\zeta =\zeta \left(x\right)$ is the source of fluctuation (noise term), ${\partial }_t=\partial /\partial t$ is the time derivative, ${\nabla }^2={\sum }_{i=1}^d\partial /\partial x_i^2$ is the Laplace operator in d spatial dimensions, $D_0$ is the diffusion parameter, $u_0$ plays the role of small coupling constant, and ${\tau }_0$ represents the deviation from the criticality. For instance, in directed bond (site) percolation on the lattice [24], the parameter $\tau =p-p_c$ corresponds to the deviation of the probability of open bonds p from its critical value $p_c$. In critical static theory, the temperature variable plays an analogous role in the standard ${\phi }^4$-model [13,32]. As has already been mentioned, in actual calculations we employ dimensional regularization; therefore, in all relevant expressions we explicitly write d as a general spatial dimension.

From a mathematical point of view, SDE (128) with multiplicative noise is ambiguous and should be augmented with proper interpretation. The usual choice is the Itô representation [46,57]. The noise term $\zeta \left(x\right)$ is assumed to be a Gaussian random field with zero mean $\langle \zeta \rangle =0$ and a pair correlation function of the following form:

$$\langle \zeta \left(x\right)\zeta \left(x^{\prime }\right)\rangle =u_0{D}_0\delta (x-x^{\prime }),$$

(129)

where $\delta \left(x\right)=\delta \left(t\right){\delta }^{\left(d\right)}\left(x\right)$ is the $(d+1)$-dimensional Dirac delta function. The parenthesis $\langle \cdots \rangle$ represents an averaging procedure over all possible realizations of the given stochastic process (Gaussian noise distribution). The multiplicative property of the noise variable in the governing equation (Equation (128)) is a direct consequence of the presence of a unique absorbing state in which all fluctuations must stop when this state is reached.

In the application of the RG method, it is advantageous to rewrite the Langevin formulation (128)–(129) in terms of functional integrals. The subsequent Janssen–De Dominicis action functional of the DP process [11,13,22,25] takes the following form:

$${\mathcal{S}}_{\mathrm{DP}}=\tilde{\psi }(-{\partial }_t+D_0{\nabla }^2-D_0{\tau }_0)\psi +\frac{D_0{u}_0}{2}({\tilde{\psi }}^2\psi -\tilde{\psi }{\psi }^2),$$

(130)

where $\tilde{\psi }=\tilde{\psi }\left(x\right)$ is an auxiliary response field and integrals over time and space variables are implicitly assumed (compare with the general expression in Equation (122)).

As has been pointed out in the introduction, non-equilibrium systems might exhibit scale-free behavior. Then, power-law dependencies are expected in scaling functions with exponents that depend only on gross features of the system, e.g., space dimension and symmetry properties [24,25]. There are certain differences between dynamic (as DP and BAP) and static models. In contrast to the equilibrium system, there is an important extension of scaling analysis. Dynamical models are described by two independent correlation lengths; in addition to the related spatial correlations, there is an additional length scale ${\xi }_{\Vert }$ describing temporal correlations. Typically, in statistical physics the models are not isotropic in the sense of models in relativistically invariant theories. On the contrary, they exhibit so-called strong anisotropic scaling, with the time variable playing a distinguished role [11,15].

Typically, we expect system correlation lengths in the vicinity of the critical point to diverge satisfying some power-law dependence ${\xi }_{\perp }\sim |{\tau }_0{|}^{-{\nu }_{\perp }},{\xi }_{\Vert }\sim {|{\tau }_0|}^{-{\nu }_{\Vert }}$, where ${\tau }_0$ is as defined in Equation (128). The exponents ${\nu }_{\perp }$ and ${\nu }_{\Vert }$ usually acquire different values. For a given universality class of the dynamic model, it is reasonable to define the dynamic exponent z as a fraction $z={\nu }_{\Vert }/{\nu }_{\perp }$. Then, in the critical region, the correlation lengths in the time and space direction are simply related as ${\xi }_{\Vert }\sim {\xi }_{\perp }^z$. From this, we can infer that the time and space scales are related by a similar relation $t\sim {|x|}^z$. In physics, two specific values are commonly encountered: $z=1$, corresponding to Lorentz-covariant theory, and $z=2$, corresponding to diffusive spreading. Subdiffusive propagation is characterized by values of the exponent $z>2$ and superdiffusive propagation by values $1<z<2$. As with any other static critical exponent, z is a universal parameter characterizing a given universality class.

One additional exponent is needed for complete knowledge of the scaling behavior of the DP process at criticality. Usually, this is the exponent $\beta$ that quantitatively characterizes the scaling of the order parameter $\rho$, i.e., the density of particles in the active phase of a macroscopic system in the critical region. The asymptotic form of this relation reads $\rho \sim {(p_c-p)}^{\beta }\propto {(-{\tau }_0)}^{\beta }$. It is possible to directly apply scaling analysis [24], with the resulting critical exponents being $\beta =1,{\nu }_{\perp }=\frac{1}{2},{\nu }_{\Vert }=1$. These values correspond to the mean-field theory, which is valid above the upper critical dimension $d_c=4$ [13].

The natural task is to calculate these exponents in physically relevant space dimensions. In this context, various methods have been used in the past. In particular, analytical results for the critical exponents of the DP process are currently known up to the two-loop approximation [11,13]. In general, perturbative computations are feasible only for a few first terms in loop expansion, as the complexity of the calculation increases dramatically with each loop order. Practically all non-equilibrium models are analyzed to the two-loop order. Therefore, one of the practical aims we are pursuing is to develop numerical and algorithmic techniques to obtain reliable higher-loop corrections. Recently, we were able to improve existing results up to the three-loop order [35,59,60]. The final results can be found elsewhere [36], together with all relevant technical details. The third-order perturbation theory presents a non-trivial task. It can be argued that its difficulty is an order of magnitude larger than the difficulty of the previous order. First of all, the number of Feynman graphs is much greater, and there are technical issues related to proper extraction of the divergent parts and numerical evaluation of the remaining parts of these Feynman diagrams. For illustration purposes, let us consider the renormalized two-point response function

$$W_{\psi \tilde{\psi }}=\langle \psi \left(x\right)\tilde{\psi }\left(0\right)\rangle .$$

(131)

This quantity satisfies the RG equation, which takes the form of the first-order partial differential equation

$$\left(\mu {\partial }_{\mu }+{\beta }_g{\partial }_g-\tau {\gamma }_{\tau }{\partial }_{\tau }-D{\gamma }_D{\partial }_D-{\gamma }_{\psi }-{\gamma }_{\tilde{\psi }}\right)W_{\psi \tilde{\psi }}=0,$$

(132)

where ${\beta }_g$ is the beta function corresponding to the charge $g=u^2$ and ${\gamma }_Q,Q\in \{\tau ,D,\psi ,\tilde{\psi }\}$ are anomalous dimensions of parameters and fields.

RG functions can be determined by perturbation theory, though such a procedure is usually accompanied by technical difficulties. The critical region function (132) takes the asymptotic form

$$W_{\psi \tilde{\psi }}(t,x,\tau )=t^{-2{\mathrm{\Delta }}_{\psi }/{\mathrm{\Delta }}_{\omega }}f\left(\frac{x}{t^{1/{\mathrm{\Delta }}_{\omega }}},\frac{\tau }{t^{-{\mathrm{\Delta }}_{\tau /{\mathrm{\Delta }}_{\omega }}}}\right).$$

(133)

We stress that the critical exponents ${\nu }_{\perp }$, ${\nu }_{\Vert }$, and z defined previously can be expressed through the critical dimensions ${\mathrm{\Delta }}_{\omega }$, ${\mathrm{\Delta }}_{\tau }$, and their ratios, as follows: ${\nu }_{\perp }=1/{\mathrm{\Delta }}_{\tau }$, ${\nu }_{\Vert }={\mathrm{\Delta }}_{\omega }/{\mathrm{\Delta }}_{\tau }$ and $z={\mathrm{\Delta }}_{\omega }$.

In the case of the DP process on a lattice [24], the response function is the primary building block for a definition of the number of active sites $N(t,\tau )$ propagating from the original seed

$$N(t,\tau )=\int {\mathrm{d}}^{d}x{W}_{\psi {\psi }^{\prime }}(t,x,\tau )\sim t^{\theta }F\left(\tau t^{{\mathrm{\Delta }}_{\tau }/{\mathrm{\Delta }}_{\omega }}\right).$$

(134)

Note here that this quantity is central for seed simulations [24]. The latter relation in Equation (134) stands for the corresponding asymptotic behavior controlled by the critical exponent $\theta$ and valid at criticality. The critical exponent $\theta$ can be expressed using the critical dimensions of the field $\psi$ and frequency $\omega$ through the following formula:

$$\theta \equiv \frac{d-2{\mathrm{\Delta }}_{\psi }}{{\mathrm{\Delta }}_{\omega }}.$$

(135)

Another dynamic quantity is the mean square radius $R^2\left(t\right)$ of propagating particles from the origin at time $t=0$, which can be introduced by the following formula

$$R^2(t,\tau )=\frac{\int {\mathrm{d}}^{d}x{x}^2{W}_{\psi {\psi }^{\prime }}(t,x,\tau )}{\int {\mathrm{d}}^{d}x{W}_{\psi {\psi }^{\prime }}(t,x,\tau )}\sim t^{2/z},$$

(136)

where $z={\mathrm{\Delta }}_{\omega }$ and the latter expression holds in the vicinity of critical point, which corresponds to the limit $\tau \to 0$.

We additionally present the results of our three-loop calculations [34,35,60] up to ${\epsilon }^3$ precision:

$$\begin{array}{cc}\hfill z& =2-\frac{\epsilon }{12}-0.02920905\left(2\right){\epsilon }^2+0.029206\left(4\right){\epsilon }^3,\hfill \end{array}$$

(137)

$$\begin{array}{cc}\hfill {\nu }_{\Vert }& =1+\frac{\epsilon }{12}+0.02238280\left(2\right){\epsilon }^2-0.008168\left(3\right){\epsilon }^3,\hfill \end{array}$$

(138)

$$\begin{array}{cc}\hfill {\nu }_{\perp }& =\frac{1}{2}+\frac{\epsilon }{16}+0.021097830\left(11\right){\epsilon }^2-0.009594\left(2\right){\epsilon }^3,\hfill \end{array}$$

(139)

$$\begin{array}{cc}\hfill \theta & =\frac{\epsilon }{12}+0.037509726\left(12\right){\epsilon }^2-0.032978\left(3\right){\epsilon }^3,\hfill \end{array}$$

(140)

where we have explicitly provided the standard error stemming from the numerical calculations. Note that it is expressed in terms of the least significant digit(s). Details of these calculations will be published elsewhere shortly [36]; here, we only note that in the two-loop approximation our calculations coincide with the known results obtained in [13,22].

7. Isotropic Percolation

The construction of effective action for IP follows analogous steps to the DP process described in the previous Section 6. The major difference lies in the necessity of introducing additional slow variables and related dynamic considerations. The governing principles are now as follows [13,28,29]:

• The relevant dynamic variables of the model are fields $n(t,x)$ and $m(t,x)$, which correspond to the density of active sites and the density of the residual debris, respectively. The actual relation between fields n and m reads $m(t,x)\propto {\int }_{-\infty }^t\mathrm{d}{t}^{\prime }n(t^{\prime },x)$.
• The dynamics of the model can be represented by means of a Markov process.
• In the susceptible medium, the activity is localized and the needed reaction rate depends on n and m, respectively. After a certain period of time, the active site becomes inactive debris.
• The system without active sites, i.e., states for which the condition $n=0$ is satisfied with an arbitrary configuration of m, are called absorbing. The active population is effectively extinct for these states, and the percolation process terminates.
• The disease proliferates between neighboring sites through a diffusion process.
• Neglected microscopic degrees of freedom can be effectively modeled through stochastic variable, which reflects the existence of absorbing states.

Using this heuristic reasoning and immediately taking into account the RG properties of the ensuing model, we can formulate a stochastic model for IP [13,29]. The role of the governing dynamic variable is played by the coarse-grained field $\psi =\psi (t,x)$, which corresponds to the density of active sites. Alternatively, the field $\psi$ might be understood as an order parameter of the model. Positive values of $\psi$ corresponds to an active state; in contrast, in the absorbing state we have $\psi =0$ everywhere.

The dynamic equation for $\psi$ takes the form

$${\partial }_t\psi =D_0{\nabla }^2\psi -D_0{\tau }_0\psi -D_0{u}_0\psi \varphi +\sqrt{\psi }\zeta ,$$

(141)

where $\zeta =\zeta (t,x)$ denotes the source of fluctuation (noise term), ${\partial }_t=\partial /\partial t$ is the time derivative, ${\nabla }^2$ is the Laplace operator in d spatial dimensions, $D_0$ is the diffusion parameter, $u_0$ plays the role of the small coupling constant, and ${\tau }_0$ represents a deviation from the critical value for the injected probability and plays an analogous role to the case of a DP process or ${\phi }^4-$ model [13,32]. Here, it is convenient to write the interaction term in the action (141) as a product $D_0{u}_0$, as $u_0$ plays the role of the charge in perturbation theory. In the RG calculation, it is necessary to distinguish between unrenormalized parameters (with the subscript “0”) and renormalized quantities (without the subscript “0”), as we do here. The debris field $\varphi$ is related to the field $\psi$ by the time integral

$$\varphi (t,x)=D_0\underset{-\infty }{\overset{t}{\int }}d{t}^{\prime }\psi (t^{\prime },x).$$

(142)

We again consider the Gaussian noise variable $\zeta$ with zero mean $\langle \zeta \rangle =0$ and a pair correlation function selected in the form

$$〈\zeta (t^{\prime },x^{\prime })\zeta (t,x)〉=u_0{D}_0\delta (t^{\prime }-t)\delta (x^{\prime }-x),$$

(143)

where the brackets stand for averaging over stochastic realizations. The form of the noise term in Equation (141) is evidently consistent with the existence of absorbing states; thus, the SDE (141) corresponds to a model with multiplicative noise.

In order to take advantage of the RG method, let us recast the stochastic problem (141)–(143) into the corresponding field-theoretic formulation using the approach described in Section 5. The ensuing Janssen–De Dominicis action functional [56,57] takes the form

$$\mathcal{S}=\tilde{\psi }\left(-{\partial }_t+D_0{\nabla }^2-D_0{\tau }_0\right)\psi +\frac{u_0{D}_0}{2}\left({\tilde{\psi }}^2\psi -2\tilde{\psi }\psi \varphi \right),$$

(144)

where $\tilde{\psi }=\tilde{\psi }(t,x)$ corresponds to a response field. For brevity, the integrations in the action (144) over all temporal and spatial coordinates are implicitly assumed, e.g., the first term on the right-hand side coincides with $\tilde{\psi }{\partial }_t\psi =\int \mathrm{d}t\int {\mathrm{d}}^{d}x\tilde{\psi }(t,x){\partial }_t\psi (t,x)$. We remark here that the response functional (144) uniquely contains only information on universal quantities, e.g., critical exponents, and not on model-dependent parameters such as the precise value of ${\tau }_c$. Furthermore, the action (144) possesses a useful asymptotic property, namely, invariance with respect to the time-reversal symmetry (dual transformation)

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

(145)

which is useful in actual perturbative calculations.

The model (144) is amenable to RG analysis; for our purposes, we aim to determine universal quantities in the form of asymptotic series in a formally small expansion parameter $\epsilon$ defined by the relation

$$d=6-2\epsilon ,$$

(146)

where $d_c=6$ is an upper critical dimension of the model.

From a practical point of view, it is sufficient to investigate only those 1PI Feynman diagrams ${\mathrm{\Gamma }}_{(n_{\tilde{\psi }},n_{\psi })}$ [25,32]. Moreover, dynamic models [11,25], in contrast to static models, require two diverse scales: a frequency scale $d_{\omega }\left[Q\right]$ and a momentum scale $d_k\left[Q\right]$ for arbitrary quantity Q in the action (response functional). The dimensions of quantity Q are determined using the standard normalization conditions

$$\begin{array}{cccc}\hfill & d_k\left[k\right]=-d_k\left[x\right]=1,\hfill & \hfill d_k\left[\omega \right]& =d_k\left[t\right]=0,\hfill \\ \hfill & d_{\omega }\left[k\right]=d_{\omega }\left[x\right]=0,\hfill & \hfill d_{\omega }\left[\omega \right]& =-d_{\omega }\left[t\right]=1\hfill \end{array}$$

(147)

together with a requirement for the response functional to be a dimensionless quantity. Using the information about $d_{\omega }\left[Q\right]$ and $d_k\left[Q\right]$, the total canonical dimension $d\left[Q\right]$ can be expressed as

$$d\left[Q\right]=d_k\left[Q\right]+2d_{\omega }\left[Q\right],$$

(148)

the exact form of which is received by comparing the IR of the most important terms (${\partial }_t\propto {\nabla }^2$) in the response functional (144). The total canonical dimension $d\left[Q\right]$ for the dynamical models represents the same quantities as the momentum dimension in static systems. The dimensions of all bare parameters and fields for the isotropic action are listed in

Table 1. Canonical dimensions of the bare parameters and fields for the IP response functional (144).

Q	$\mathit{k},{\mathit{\partial }}_{\mathit{x}}$	$\mathit{\omega },{\mathit{\partial }}_{\mathit{t}}$	$\tilde{\mathit{\psi }},\mathit{\varphi }$	$\mathit{\psi }$	${\mathit{D}}_0$	${\mathit{\tau }}_0$	${\mathit{u}}_0$
$d_k\left[Q\right]$	1	0	$\frac{d-2}{2}$	$\frac{d+2}{2}$	$-2$	2	$\frac{6-d}{2}$
$d_{\omega }\left[Q\right]$	0	1	0	0	1	0	0
$d\left[Q\right]$	1	2	$\frac{d-2}{2}$	$\frac{d+2}{2}$	0	2	$\frac{6-d}{2}$

. Dimensional analysis shows that the IP action is logarithmic (when the charge $u_0$ becomes dimensionless) at $\epsilon =0$, and manifests as pole terms with respect to $\epsilon$.

The action (144) is amenable to standard perturbative analysis, which yields expansion in Feynman diagrams. The fundamental perturbative components are a propagator and two interaction vertices. For practical calculations, it is convenient to consider a propagator ${\langle \psi \tilde{\psi }\rangle }_0$ in time–momentum representation, in which it takes the following form:

$${〈\psi \tilde{\psi }〉}_0=\theta (t-t^{\prime })exp\left[-D_0(p^2+{\tau }_0)(t-t^{\prime })\right].$$

(149)

The rules for Feynman graphs are supplemented by momentum conservation at each interaction vertex and integration over all time coordinates in a diagram. Note that another acceptable way of defining perturbation elements is possible [30].

The Dyson equation for the self-energy Green function in the frequency–momentum representation has the general form

$${\mathrm{\Gamma }}_{(1,1)}=i\omega Z_1+D(p^2{Z}_2+\tau Z_3)-\mathrm{loop}\mathrm{diagrams}.$$

(150)

Furthermore, the perturbation representation for the vertex function ${\mathrm{\Gamma }}_{(2,1)}$ takes the general form

$${\mathrm{\Gamma }}_{(2,1)}=Du{\mu }^{\epsilon }{Z}_4+\mathrm{loop}\mathrm{diagrams}.$$

(151)

Thus, as a direct consequence of the time-reversal symmetry (145), only two 1PI Green functions ${\mathrm{\Gamma }}_{(1,1)}$ and ${\mathrm{\Gamma }}_{(2,1)}$ have to be considered in order to correctly renormalize the model.

Without mentioning unnecessary technical details, for which we refer the interested readers to the literature [11,13,29], we note here only that the crucial RG property of the action functional (144) is fulfilled, i.e., the model is multiplicatively renormalizable. The perturbation theory is most conveniently analyzed by means of dimensional regularization augmented with the famous $\epsilon$-expansion [25,32]. As usual, the response functional is investigated near the critical point, where the mass-like parameter ${\tau }_0$ attains its critical value.

Direct inspection of the Feynman diagrams reveals that the coupling constant of the perturbation theory is $u_0^2$ rather than $u_0$. Consequently, it is appropriate to introduce an expansion parameter $g_0=u_0^2$ within which to carry out the perturbative analysis.

The renormalized action corresponding to the bare action (144) takes the following form:

$${\mathcal{S}}_R={\tilde{\psi }}_{{}_R}(-Z_1{\partial }_t+Z_{2}D{\nabla }^2-Z_{3}D\tau ){\psi }_{{}_R}+\frac{Du{\mu }^{\epsilon }}{2}[Z_4{\tilde{\psi }}_{{}_R}^2{\psi }_{{}_R}-2Z_5{\tilde{\psi }}_{{}_R}{\psi }_{{}_R}{\varphi }_{{}_R}].$$

(152)

On the other hand, the response functional can be obtained by multiplicative renormalization of all parameters and fields:

$$\begin{array}{cccccc}\hfill \psi & ={\psi }_{{}_R}{Z}_{\psi },\hfill & \hfill \tilde{\psi }& ={\tilde{\psi }}_{{}_R}{Z}_{\tilde{\psi }},\hfill & \hfill D_0& =D{Z}_D,\hfill \\ \hfill u_0& =u{\mu }^{ϵ}{Z}_u,\hfill & \hfill {\tau }_0& =\tau Z_{\tau }+{\tau }_c,\hfill \end{array}$$

(153)

amended with analogous relations between RG constants for terms appearing in the action (152):

$$\begin{array}{cccc}\hfill Z_1& =Z_{\psi }{Z}_{\tilde{\psi }},\hfill & \hfill \hspace{1em}\hspace{1em}\hspace{1em}{Z}_2& =Z_D{Z}_{\psi }{Z}_{\tilde{\psi }}\hfill \\ \hfill Z_3& =Z_{\tau }{Z}_D{Z}_{\psi }{Z}_{\tilde{\psi }},\hfill & \hfill \hspace{1em}\hspace{1em}\hspace{1em}{Z}_4& =Z_D{Z}_u{Z}_{\psi }{Z}_{\tilde{\psi }}^2,\hfill & \hfill \hspace{1em}\hspace{1em}\hspace{1em}{Z}_5& =Z_D^2{Z}_u{Z}_{\psi }^2{Z}_{\tilde{\psi }}.\hfill \end{array}$$

(154)

By straightforward algebraic operations, we can readily invert these relations and obtain RG constants for the fields and parameters in terms of RG constants $Z_1-Z_5$

$$\begin{array}{cccc}\hfill Z_{\psi }& =Z_1{Z}_2^{-1/2},\hfill & \hfill \hspace{1em}\hspace{1em}{Z}_{\tilde{\psi }}& =Z_2^{1/2},\hfill \\ \hfill Z_D& =Z_2{Z}_1^{-1},\hfill & \hfill \hspace{1em}\hspace{1em}{Z}_{\tau }& =Z_3{Z}_2^{-1},\hfill & \hfill \hspace{1em}\hspace{1em}{Z}_u& =Z_4{Z}_2^{-3/2}.\hfill \end{array}$$

(155)

Time-reversal symmetry (145) implies that the RG constants $Z_4$ and $Z_5$ are in fact equal. From the general theory of critical phenomena, we know that a set of independent critical exponents are associated with each universality class that completely determine the universal macroscopic behavior of the model. For the IP model, we have to determine three independent exponents. To obtain these, we can begin with the dynamic response functional (144) and follow the usual steps of perturbative RG analysis [25]. However, it can be shown that this is actually not the case [13], as two exponents can be calculated from a less complex action that corresponds to the quasi-static limit [31]. In this way, the time-dependence of the fields is eliminated and perturbative analysis is considerably simpler. In fact, IP in this limit can be related to the widely recognized Potts model [61]. The main result of this remark is that the renormalization functions ${\beta }_u$, ${\gamma }_{\psi }$ and ${\gamma }_{\tau }$ can be obtained by easier means. We note here that results in this regard are already known up to four-loop order [62].

As mentioned at the beginning of this section, to completely accomplish the three-loop calculation of critical exponents it is necessary to use the action (144). Our main task is to find the divergent parts of the 1PI Green function ${\mathrm{\Gamma }}_{\psi \tilde{\psi }}(\omega ,p)$ that are proportional to the external frequency $\omega$. To this end, we need to determine the partial derivative $\partial {\mathrm{\Gamma }}_{(1,1)}(\omega ,p)/\partial \left(i\omega \right)$, which would allow us to calculate the dynamic exponent z, which cannot be determined by means of the aforementioned quasi-static limit. This exponent describes the difference between scaling in the spatial and temporal directions [11,13,24].

In this part, we present results from one- and two-loop approximations. For future reference, we note that these graphs are conveniently described using the Nickel index for later algorithmic manipulation [63,64]. It is possible to precisely determine the internal structure of each Feynman diagram using Nickel notation, which plays a fundamental role in symbolic operations. The Nickel notation for relevant two-loop Feynman diagrams is provided in

Table 2. Representation of two-loop Feynman graphs for self-energy employing a formalism based on Nickel indices. In the first column, the number in brackets corresponds to the Feynman diagram from Equations (157)–(162).

Diagram	Nickel Index	Symmetry Coefficient
(157)	$e12|e3|33|:0\tilde{\psi }\_\psi \tilde{\psi }\_\varphi \tilde{\psi }|0\psi \_\tilde{\psi }\psi |\varphi \tilde{\psi }\_\psi \tilde{\psi }|$	1
(158)	$e12|e3|33|:0\tilde{\psi }\_\varphi \tilde{\psi }\_\psi \tilde{\psi }|0\psi \_\tilde{\psi }\psi |\varphi \tilde{\psi }\_\psi \tilde{\psi }|$	1
(159)	$e12|23|3|e|:0\tilde{\psi }\_\varphi \tilde{\psi }\_\psi \tilde{\psi }|\tilde{\psi }\varphi \_\psi \tilde{\psi }|\psi \tilde{\psi }|0\psi |$	1
(160)	$e12|23|3|e|:0\tilde{\psi }\_\varphi \tilde{\psi }\_\psi \tilde{\psi }|\tilde{\psi }\psi \_\psi \tilde{\psi }|\varphi \tilde{\psi }|0\psi |$	1
(161)	$e12|23|3|e|:0\tilde{\psi }\_\varphi \tilde{\psi }\_\psi \tilde{\psi }|\psi \tilde{\psi }\_\varphi \tilde{\psi }|\psi \tilde{\psi }|0\psi |$	1
(162)	$e12|23|3|e|:0\tilde{\psi }\_\varphi \tilde{\psi }\_\psi \tilde{\psi }|\varphi \tilde{\psi }\_\psi \tilde{\psi }|\psi \tilde{\psi }|0\psi |$	1

.

The standard evaluation of a single divergent one-loop diagram [29] yields the following contribution:

(156)

where $\mathrm{\Gamma }\left(x\right)$ is Euler’s gamma function. Note that the mass term ${\tau }_0$ has been used as an infrared regulator.

Obtaining the two-loop approximation is more demanding; in this section, there is no particular reason to provide all of thedetailed technical information. Thus, we focus only on the calculation of the RG constant and present the contributions of the Feynman diagrams to $Z_1$. These read:

(157)

(158)

(159)

(160)

(161)

(162)

The total two-loop contribution reads

$$\sum _{2-\mathrm{loop}}=\frac{\mathrm{\Gamma }(1+2\epsilon )}{{\left(4\pi \right)}^{6-2\epsilon }}\frac{u_0^4{\tau }_0^{-2\epsilon }}{64\epsilon }\left[\frac{51}{2\epsilon }+\frac{107}{12}-5ln2+\frac{9}{2}ln3\right].$$

(163)

To further simplify the analysis, it is advantageous to introduce a new charge $g_0$ in the following way:

$$g_0=u_0^2\frac{\mathrm{\Gamma }(1+\epsilon )}{{\left(4\pi \right)}^{3-\epsilon }}.$$

(164)

Using the approximate relation, where $\mathrm{\Gamma }(1+2\epsilon )={\mathrm{\Gamma }}^2(1+\epsilon )+O\left({\epsilon }^2\right)$, we finally arrive at the desired results:

$$\sum _{1-\mathrm{loop}}=-\frac{3g_0{\tau }_0^{-\epsilon }}{8\epsilon },\sum _{2-\mathrm{loop}}=\frac{g_0^2{\tau }_0^{-2\epsilon }}{64\epsilon }\left[\frac{51}{2\epsilon }+\frac{107}{12}-5ln2+\frac{9}{2}ln3\right].$$

(165)

To the two-loop precision 1PI function ${\mathrm{\Gamma }}_{(1,1)}$ can be formally written as

$${\mathrm{\Gamma }}_{(1,1)}=i\omega \left(1+g_0{\tau }_0^{-\epsilon }{\mathrm{\Gamma }}^{\left(1\right)}+g_0^2{\tau }_0^{-2\epsilon }{\mathrm{\Gamma }}^{\left(2\right)}\right).$$

(166)

The exact form of the RG constant $Z_1$ follows from standard algebraic operations [25]:

$$\begin{array}{cc}\hfill {\mathrm{\Gamma }}_{1,1}^R& =i\omega \left[Z_1+g{\left(\frac{{\mu }^2}{\tau }\right)}^{\epsilon }{\mathrm{\Gamma }}^{\left(1\right)}{Z}_1{Z}_g{Z}_{\tau }^{-\epsilon }+g^2{\left(\frac{{\mu }^2}{\tau }\right)}^{2\epsilon }{\mathrm{\Gamma }}^{\left(2\right)}{Z}_1{Z}_g^2{Z}_{\tau }^{-2\epsilon }\right]\hfill \end{array}$$

(167)

$$\begin{array}{cc}\hfill & =i\omega \left[Z_1+g{\left(\frac{{\mu }^2}{\tau }\right)}^{\epsilon }{\mathrm{\Gamma }}^{\left(1\right)}{Z}_1{Z}_g{Z}_{\tau }^{-\epsilon }+g^2{\left(\frac{{\mu }^2}{\tau }\right)}^{2\epsilon }{\mathrm{\Gamma }}^{\left(2\right)}\right]\hspace{1em}\hspace{1em}\hspace{1em}\hfill \end{array}$$

(168)

$$\begin{array}{cc}& \hspace{1em}\hspace{1em} =i\omega \left[Z_1-\frac{3}{4}\frac{1}{2\epsilon }g{\left(\frac{{\mu }^2}{\tau }\right)}^{\epsilon }\left(1+g\left(\frac{17}{8\epsilon }-\frac{5}{12}\right)\right)+g^2{\left(\frac{{\mu }^2}{\tau }\right)}^{2\epsilon }{\mathrm{\Gamma }}^{\left(2\right)}\right],\hfill \end{array}$$

(169)

where the renormalized coupling g enters via $g_0=g{\mu }^{2\epsilon }{Z}_g$ and we require that the RG constant cancel all pole terms up to the two-loop order. From this requirement, we eventually obtain the final form of the RG constant $Z_1$:

$$Z_1=1+\frac{3}{4}\frac{g}{2\epsilon }+\frac{1}{64}\frac{g^2}{2\epsilon }\left[\frac{51}{\epsilon }-\frac{227}{6}+10ln2-9ln3\right],$$

(170)

which is consistent with [13].

From the point of view of renormalization, the unacceptable divergent contribution

$$\frac{ln\left(\tau \right)}{\epsilon }$$

(171)

is eliminated in the final expression, which provides a necessary condition for the correct use of the renormalization procedure [25]. Note that terms such as (171) might appear in every order of the perturbation theory. Thus, the final cancellation of such terms (171) provides a nontrivial check of the employed calculation techniques.

For the IP model, evaluations of the dynamic exponent z have been carried out to two-loop precision. Our aim is to develop a calculation scheme that would allow us to obtain a higher three-loop approximation. To this end, we employ a recently obtained algorithmic method [65,66,67]. The method itself consists of the following steps:

• Construction of all needed three-loop 1PI diagrams, with all permissible connections among internal vertices complying with the Feynman rules (Figure 1), with two external legs for ${\mathrm{\Gamma }}_{\psi \tilde{\psi }}$ and three legs for ${\mathrm{\Gamma }}_{\psi \tilde{\psi }\tilde{\psi }}$.
• Symbolic generation of the integrand using the procedure in Python $2.7$ (GraphState and Graphine packages).
• Calculation of the corresponding part of the Feynman diagrams using the Vegas algorithm.

Compared to directed percolation [34,35] (discussed in Section 6), it is notable that the isotropic percolation model represents a much more demanding task, and would probably present a more difficult problem from a technical point of view as well. In this regard, we have performed the initial steps of the three-loop calculation of the IP model. In particular, existing algorithms for the DP were successfully generalized, and all Feynman diagrams required for a complete three-loop approximation were listed in systematic fashion. The number of three-loop diagrams can be found in

Table 3. The number of three-loop 1PI graphs for the IP model.

IP	1-Loop	2-Loop	3-Loop
$\langle \tilde{\psi }\psi \rangle$	1	6	93
$\langle \tilde{\psi }\tilde{\psi }\psi \rangle$	1	36	957

, and their topological graphical representation is shown in Figure 2 for the 1PI Green function ${\mathrm{\Gamma }}_{(1,1)}$ and in Figure 3 for the 1PI Green function ${\mathrm{\Gamma }}_{(2,1)}$. For the sake of brevity, we do not explicitly write out all permissible diagrams, and restrict ourselves to different topology structures. The total number of independent possibilities for different configurations of internal vertices is indicated in the bracket after each topological structure. In other words, it represents the number of possible diagrams for a given topological structure.

8. Conclusions

In this paper, we have studied the universal properties of two different percolation processes within the field-theoretic approach and performed the main steps of the UV renormalization analysis. We present the results for critical exponents of the directed percolation universality class up to the third-loop approximation. For both models, we employed the MS scheme, in which a small expansion parameter is represented by $\epsilon$, which controls the deviation from the upper-critical space dimension $d_c=4$ (DP) or $d_c=6$ for isotropic processes. The obtained results agree very well with the two-loop calculations obtained previously in [13,22]. For the realistic dimension $d=3(\epsilon =1)$, the resummation technique using Padé approximants was performed. The achieved results are in excellent agreement with the Monte Carlo prediction [36]. We believe that this paper represents a step towards the enhancement of multi-loop calculations in the area of nonequilibrium physics. We hope that it will help with future investigations of other more difficult models.