Renormalization Approach to the Gribov Process: Numerical Evaluation of Critical Exponents in Two Subtraction Schemes

We study universal quantities characterizing the second order phase transition in the Gribov process. To this end, we use numerical methods for the calculation of the renormalization group functions up to two-loop order in perturbation theory in the famous ε-expansion. Within this procedure the anomalous dimensions are evaluated using two different subtraction schemes: the minimal subtraction scheme and the null-momentum scheme. Numerical calculation of integrals was done on the HybriLIT cluster using the Vegas algorithm from the CUBA library. The comparison with existing analytic calculations shows that the minimal subtraction scheme yields more precise results.


Introduction
The Gribov process is also known in the literature [1][2][3][4][5] as an epidemic with recovery or directed (bond) percolation process. Equivalent processes were studied several times as chemical reactions system known in the literature as Schlögl's first equation [6]. The Gribov process belongs to a class of systems that exhibits a nonequilibrium phase transition between the active and absorbing states. The most prominent feature of such transitions is the nontrivial scaling behavior, which can be fully quantitatively characterized by several independent critical exponents. These exponents can be estimated in a controllable manner by various theoretical techniques.
An important procedure for the calculation of universal quantities is the renormalization group (RG) approach [7]. Using this technique, one can systematically build numerical methods for approximate calculation of the critical exponents. In this paper, the renormalization procedure is carried out in the framework of the ε-expansion, where ε stands for the deviation from the upper critical dimension d c = 4. Numerical calculations are performed in two different subtraction schemes. In the first case, the minimal subtraction (MS) scheme is employed and contributions to the renormalization constants for each Feynman diagram are determined. The second scheme is the Normalization Point (NP) scheme [8][9][10][11][12][13] which is a variant of the zero momentum subtraction (ZM) scheme [7], where the e-mail: lukas.mizisin@gmail.com computation of the RG constants is skipped and contributions of diagrams to the anomalous dimension are calculated. The choice of the scheme is arbitrary since the universal quantities in the form of the ε-expansion are independent of the choice of the renormalization scheme. The main goal of this paper is to develop numerical methods of evaluation of the critical exponents for different subtraction schemes and to determine the most efficient one.

Model and Renormalization
In the terminology of reaction-diffusion problems [1][2][3], the Gribov process describes the time evolution of classical objects R that diffuse in a d-dimensional continuous space. At the same time, they are subject to the following chemical reactions: where the number of objects M and N are kept constant, A and B are chemically inert, while k − 1 , k + 1 , k − 2 , k = 2 denote rate constants. Using the classical master equation, the systems can be easily reformulated in terms of bosonic ladder creation and annihilation operators. Via coherent-state path integrals, one finally arrives at the following field theoretic action functional [4]: where ∂ t = ∂/∂t, ∇ 2 is the Laplace operator, D 0 is the diffusion constant, and ψ andψ correspond to the coherent states of the bosonic creation and annihilation operators. Integrations over time and spatial variables in (2) are implied.
For the later use of ultraviolet renormalization [7], it is convenient to rewrite the action functional as follows: where λ 0 is a positive coupling constant and τ 0 measures the deviation from the threshold value of the critical probability (an analog of the critical temperature in static models). The last term in the action (2) (four point interaction) was neglected due to its infrared irrelevance. The model is studied by field theoretical RG [7] near its critical dimension d c = 4 in the vicinity of the second order phase transition by means of ε-expansion, which is for convenience defined as follows 2ε = 4 − d. The renormalized action functional takes the form where µ is the renormalization mass and Z i , i = 1, 2, 3, 4 are the renormalization constants. The model possesses symmetry with respect to the so-called rapidity symmetry [5], and as a direct consequence two triple vertices are renormalized by the same renormalization constant Z 4 . On the other hand, the renormalized action could be obtained by multiplicative renormalization of the fields and parameters The role of the coupling constant in perturbation theory is played by λ 2 . To keep the notation simple, we define a new charge g ≡ 2λ 2 /((4π) d/2 Γ(d/2)), where Γ(x) is the gamma function. An important role is played by the choice of the subtraction scheme for the calculation of the RG functions. We have considered two different schemes for the evaluation of universal quantities. In the MS scheme [7], the renormalization constants can depend only on the poles in ε and the coupling constant g, i.e. Z i = 1 + ∞ n=1 g n n k=1 c nk ε −k , where c nk is a pure number (not depending on other parameters). The anomalous dimensions in the MS scheme can be determined through the renormalization constants as follows [7]: To calculate the anomalous dimensions in the MS scheme, the residues at poles in ε in renormalization constants were calculated with Sector Decomposition method [14]. The second subtraction procedure is the NP scheme [8][9][10][11]. Within this scheme, the anomalous dimensions γ i are directly expressed in terms of some Feynman diagrams which are free of divergences and thus the calculation of the renormalization constants can be entirely skipped (for details see [10,12,13]).
Integrals were calculated by the modified Monte Carlo methods, which were developed for calculations of integrals in multidimensional space. The integration is the same as in the original Monte Carlo method. The only difference lies in the low-discrepancy sequence of random numbers. The numerical calculation was carried out by the Vegas algorithm within a new implementation. A detailed description of the Vegas algorithm can be found in [15].

Results
The unrenormalized Green functions are independent of the momentum scale µ. This fact implies (for both the MS and NP schemes) the basic RG differential equations for the renormalized one-irreducible Green functions [7] µ∂ µ + β∂ g − τγ τ ∂ τ − Dγ D ∂ D − n ψ γ ψ − nψγψ Γ R n ψ nψ = 0, where n ψ , nψ are the numbers of the corresponding fields entering the Green function under consideration. The anomalous dimensions up to the second order of the perturbation theory take the following form: For the MS scheme, the coefficients C (k) i are numbers. For the NP scheme, they can depend on ε, in the two loop case C (2) i is a number, while C (1) i contains a term proportional to ε. The RG method predicts [4,7] that the large scale behavior is determined by the infrared attractive fixed point, the coordinate of which is found from the condition β g (g * ) = 0. The main objects of interest in our study are the critical exponents defined as: η= 2γ ψ (g * ).
The dynamic exponent z is related to the mean square radius, whereas the second exponent η to the survival probability of the percolating cluster. These quantities are universal. This means that their ε-expansion does not depend on the subtraction scheme (contrary to the anomalous dimensions which are non-universal as well as the coordinate of the fixed point). The expansions of these quantities take the following forms: where C z and C η are displayed in Table 1.

Conclusion
Numerical calculations were carried out for the Gribov process up to the second order perturbation theory by two different subtraction schemes. In both cases the results (anomalous dimensions, critical exponents) are in good agreement with the analytical calculations. However, the values are more precise for the MS scheme and this method will be more suitable for numerical calculation in the higher order of the perturbation theory.