Quantum-Field Multiloop Calculations in Critical Dynamics

Abstract

The quantum-field renormalization group method is one of the most efficient and powerful tools for studying critical and scaling phenomena in interacting many-particle systems. The multiloop Feynman diagrams underpin the specific implementation of the renormalization group program. In recent years, multiloop computation has had a significant breakthrough in both static and dynamic models of critical behavior. In the paper, we focus on the state-of-the-art computational techniques for critical dynamic diagrams and the results obtained with their help. The generic nature of the evaluated physical observables in a wide class of field models is manifested in the asymptotic character of perturbation expansions. Thus, the Borel resummation of series is required to process multiloop results. Such a procedure also enables one to take high-order contributions into consideration properly. The paper outlines the resummation framework in dynamic models and the circumstances in which it can be useful. An important resummation criterion is the properties of the higher-order asymptotics of the perturbation theory. In static theories, these properties are determined by the method of instanton analysis. A similar approach is applicable in critical dynamics models. We describe the calculation of these asymptotics in dynamical models and present the results of the corresponding resummation.

1. Introduction

Continuous phase transitions controlled by an order parameter play an essential role in the description of condensed systems. This applies equally to the well-known phase transitions, which have been included in textbooks for hundreds of years, and to the exotic ones, obtained relatively recently in powerful fields or at superhigh pressures. A wide variety of critical points are known: these points represent a set of thermodynamic parameters near which the behavior of the system is non-analytical. The anomalies demonstrate “universality”, expressed, for example, in the existence of universal critical indices. The latter is determined by rather general properties of the system under consideration, for example, by the tensor structure of the order parameter and space dimension.

Despite a rich history, the theory of critical behavior began to develop only in the 20th century, when L.D. Landau proposed the simplest description of continuous phase transitions based on the analysis of symmetry and dimension [1,2]. However, Landau’s theory was based on an incorrect account of fluctuations in the considered systems and therefore gave incorrect predictions. Significant progress in the theory of critical phenomena was associated with the emergence of the so-called Wilson renormalization group (RG) [3,4,5] and the $\epsilon$-expansion. Here, for the first time, it was possible to take into account the renormalization of critical exponents due to fluctuation effects. Similar techniques have yielded both the static properties of the critical behavior and the dynamic properties of equilibrium fluctuations, which describe the growth of the relaxation times in the vicinity of the critical point. Unfortunately, the consistent application of the Wilsonian RG turned out to be associated with considerable computational difficulties. Therefore, in this computational formalism, it was only possible to obtain significant results within the scope of the first-order perturbation theory. Of course, such accuracy did not correspond to experimental expectations and did not allow making practically important predictions, which the theory of condensed matter so badly needed. Fortunately, it turned out that the equivalent and more advanced technique of quantum field RG [6,7,8] solves the problem much better. On its basis, the development of theoretical studies of critical behavior continued.

In the theory of critical behavior, there are quite a few different renormalization group approaches. Among them, in addition to the Wilsonian RG mentioned above, we will name the RG in the real space dimension [9,10,11] as well as the functional RG [12]. Our article is mainly devoted to the standard quantum field RG approach. It allows the calculation of multiloop diagrams most efficiently with proper use of the arbitrariness of renormalization scheme. The result for critical exponents is constructed as a series in $\epsilon$ parameter, which is the deviation of the space dimension from the logarithmic dimension, for which the Lagrangian coupling constant has specific scaling properties. In presenting the principles of multiloop calculations, we assume that the reader is familiar with quantum field renormalization group methods in statistical physics as well as with the theory of ultraviolet renormalization. For our part, we recommend here the classic monographs [6,7,8]. In the process of time-consuming multiloop calculations, it turned out that progress is influenced not only by technical methods but also by the regularization and renormalization schemes used. For dimensional regularization, the minimal subtraction scheme in the massless theory is the most efficient (see, for example, [8] and references there). With this approach, the best results of multiloop calculations in critical dynamics are achieved.

In this article, we will refer to the results obtained in the study of “critical dynamics”. Here, we have in mind the dynamics of fluctuations near the equilibrium state in the vicinity of the critical point. It should be noted that the presented methods for calculating diagrams and for resummation of perturbation series are also successfully applied to models studied by the renormalization group in the real space dimension and to microscopic quantum models. Furthermore, they are also successfully used in various stochastic models that do not tend to the Gibbs equilibrium in the limit of long relaxation time. However, these models still need to sufficiently implement multiloop calculations.

The article is organized as follows. Section 1.1 outlines the main diagram calculation techniques developed for static field theories. The basic models of critical dynamics are described in Section 2. Section 3 presents the currently developed methods for calculating dynamics diagrams and a simple illustration of some calculations. The fluctuation–dissipation theorem leads to the connections between the static and corresponding dynamics models. This leads to significant simplifications in the multiloop dynamics calculations. This is a basis of diagram reduction for the critical dynamics models described in Section 4, where the result of five-loop calculations for the $O\left(n\right)$ symmetric model A is presented. Section 5 is devoted to methods of resuming divergent series in critical dynamics models. The resummed values of dynamic critical exponents z for the model A are presented.

1.1. Basic Diagram Calculation Methods

The standard ${\varphi }^4$ model, illustrating multiloop calculation methods, is considered as the theory of the scalar field $\varphi \left(x\right)$ over the Euclidean space of dimension $d=4-\epsilon$. The system is described by the action

$$\begin{array}{c}\hfill S=\int {\mathrm{d}}^{d}x\left(\frac{1}{2}{(\nabla \varphi )}^2+\frac{{\tau }_0}{2}{\varphi }^2+\frac{g_0}{4!}{\varphi }^4\right).\end{array}$$

(1)

As is customary in the literature, in similar expressions below, we will omit integration over spatial coordinates and time. Although the investigated functional has the meaning of a Hamiltonian, in the quantum field approach, it is denoted by S, inheriting the notation from the relativistic QFT in the Minkowski space. The model under consideration is multiplicatively renormalizable. The renormalization is performed by stretching the fields $\varphi \to Z_{\varphi }\varphi$ and the parameters $g_0=g{\mu }^{\epsilon }{Z}_g$, ${\tau }_0=\tau Z_{\tau }$, where $\mu$ is the renormalization mass. In the minimal subtraction scheme, the renormalization constants have purely pole contributions generated by divergent diagrams

$$Z_i=1+\frac{c_i^{\left(1\right)}\left(g\right)}{\epsilon }+\frac{c_i^{\left(2\right)}\left(g\right)}{{\epsilon }^2}+\dots ,i\in \{\varphi ,g,\tau \}.$$

(2)

The method of studying critical behavior is based on the RG equation (see, for example, [8] and references therein)

$$[s{\partial }_s-\beta {\partial }_g+{\Delta }_{\tau }\tau {\partial }_{\tau }+{\Delta }_n]G_n^R=0,$$

where $s=p/\mu$ is the scale parameter and $G_n^R$ is the n-point correlation function. The so-called renormalization group functions $\beta$ and ${\Delta }_i=d_i+{\gamma }_i$ (here $d_i$ are the canonical dimensions of i, ${\gamma }_i$ correspond to the anomalous dimensions) are expressed in a known way in terms of the renormalization constants [8]. Small s corresponds to the large-scale—or infra-red (IR)—behavior of the system’s correlators. The IR region is reached at $p\to 0$, for example, when p is the only argument of the correlation function in the momentum representation. Thus, the standard problem in this approach is to calculate the contributions of all diagrams to the renormalization constants for a given order of the perturbation expansion. Further, a fairly routine recalculation of them into the $\epsilon$-expansion coefficients of critical exponents is required. Let us list the main techniques for calculating diagrams that have been developed in statics field theories [8,13].

Alpha-representation. It consists of the integral representation of a power factor in the denominator of the diagram in the form

$$\frac{1}{a^{\alpha }}=\frac{1}{\Gamma \left(\alpha \right)}\underset{0}{\overset{\infty }{\int }}\mathrm{d}t{t}^{\alpha -1}{e}^{-at},$$

(3)

and the subsequent replacement of the orders of integration in the diagram. As a result, the d-dimensional integrals over the momenta are calculated first, and the parametric integrals are processed at the end.

Feynman representation.

$$\frac{1}{a_1^{{\alpha }_1}\dots a_n^{{\alpha }_n}}=\frac{\Gamma ({\alpha }_1+\cdots +{\alpha }_n)}{\Gamma \left({\alpha }_1\right)\dots \Gamma \left({\alpha }_n\right)}\underset{0}{\overset{1}{\int }}\cdots \underset{0}{\overset{1}{\int }}\prod _{j=1}^n\mathrm{d}{z}_j\frac{z_j^{{\alpha }_j-1}\delta ({\sum }_{i=1}^n{z}_i-1)}{{\left({\sum }_{i=1}^n{z}_i{a}_i\right)}^{{\alpha }_1+\cdots +{\alpha }_n}}.$$

(4)

Let us consider an example of applying the mentioned techniques to calculate the pole singularity of the following diagram in the dimensional regularization $d=4-\epsilon$

The calculation of the last expression is simplified if we recall the independence of the calculated pole contribution in $\epsilon$ from the external momentum p. Setting $p=0$ and using the $\alpha$-representation, we obtain

The pole contribution in $\epsilon$ in this expression does not depend on $\tau$, which is the basis for using the “massless” calculation scheme with $\tau =0$. In more complex diagrams, after calculating the d-dimensional integrals over the arguments of the propagators, a nontrivial problem is the evaluation of the integrals over the remaining set of parameters. Let us designate this set of parameters $z_i$ and $t_i$ and list the main techniques employed to compute them.

Sector decomposition method. Within this machinery [14] the parametric multiple integral is divided into sectors $z_1\ge \{z_2,\dots ,z_n\}$, $z_2\ge \{z_1,z_3,\dots ,z_n\}$, etc., and then the appropriate scaling of variables $z_i$ is performed in each sector. As a result, all integrals again take the form $\underset{0}{\overset{1}{\int }}\cdots \underset{0}{\overset{1}{\int }}\prod \mathrm{d}{z}_i$. The integral over the large parameter z is easy to calculate. It produces poles in $\epsilon$. The sector decomposition procedure can be repeated for diagrams with divergent subgraphs to obtain a convergent integral.

Kompaniets–Panzer method. The method [15] is based on the fact that integrals over parameters $t_i$ of alpha-representations in the diagrams without divergent subgraphs or in the subgraphs can be calculated with a certain order of integration using analytical formulas for polylogarithms.

Alternating transitions to the coordinate and momentum representations. A significant simplification in analytical calculations was given by the method presented in [8], using formulas for the Fourier transform $\mathcal{F}$ of arbitrary power functions in an arbitrary space dimension

$$\frac{1}{{(x_1-x_2)}^{2\alpha }}\stackrel{\mathcal{F}}{\Rightarrow }\frac{c\left(\alpha \right)}{k^{2(d/2-\alpha )}}$$

(5)

and formulas for the inverse Fourier transform ${\mathcal{F}}^{-1}$

$$\frac{1}{k^{2\alpha }}\stackrel{{\mathcal{F}}^{-1}}{\Rightarrow }\frac{\overline{c}\left(\alpha \right)}{{(x_1-x_2)}^{2(d/2-\alpha )}}.$$

(6)

Here, we introduce magnitudes

$$c\left(\alpha \right)={\pi }^{d/2}{4}^{d/2-\alpha }\Gamma (d/2-\alpha )/\Gamma \left(\alpha \right),\overline{c}\left(\alpha \right)=c\left(\alpha \right)/{\left(2\pi \right)}^d.$$

For example, calculating a two-loop diagram using this technique looks like this

This expression explicitly contains the required $\epsilon$ pole contributions. For calculating more complex diagrams, the convolution of lines formula is also useful:

$$\int {\mathrm{d}}^{d}x\frac{1}{{(x_1-x)}^{2{\alpha }_1}}\frac{1}{{(x_2-x)}^{2{\alpha }_2}}=\frac{c\left({\alpha }_1\right)c\left({\alpha }_2\right)\overline{c}(d-{\alpha }_1-{\alpha }_2)}{{(x_1-x_2)}^{2({\alpha }_1+{\alpha }_2-d)}}.$$

(7)

Further development of this method was based on the use of the vertex uniqueness method [16] and the $R^{*}$-operation [17,18,19,20]. The described methods made it possible to evaluate the critical exponents in the ${\varphi }^4$ static model up to six-loop diagrams [15] and present the critical indices in the form of the $\epsilon$-expansion up to the corresponding order. As a rule, the $\epsilon$-series obtained in the static models turns out to be asymptotic, so the resulting segment of the series required a resummation procedure.

2. Models of Critical Dynamics

The study of static properties of correlated equilibrium systems in the vicinity of the critical point or close to the continuous phase transitions is a major area of statistical field theory. However, time-dependent phenomena, e.g., relaxation dynamics, are more diverse and no less interesting because they capture inter-mode interaction effects, which die out in equilibrium. A key experimental manifestation of dynamic effects is critical slowing down in the relaxation time, which increases abnormally. The most investigated and known models of critical dynamics are listed in [21] and in more detail discussed in [8,22,23].

Critical dynamics is described by a set of physically meaningful fields ${\varphi }_a$, which obey the system of stochastic differential equations. For a wide class of systems, these equations can be written in the general Langevin form

$$\begin{array}{c}\hfill {\partial }_t{\varphi }_a=-({\alpha }_{ab}+{\beta }_{ab})\frac{\delta S}{\delta {\varphi }_b}+{\eta }_a,\end{array}$$

(8)

where $\alpha$ is a symmetric operator ${\alpha }_{ab}={\alpha }_{ba}$, $\beta$ is a skew-symmetric one ${\beta }_{ab}=-{\beta }_{ba}$ and satisfies the condition $\delta {\beta }_{ab}/\delta {\varphi }_a=0$, where summation over the index a is assumed [8]. It may be verified that this condition is sufficient to ensure relaxation towards an equilibrium state with the Gibbs distribution functional $exp(-S)$. The additive noise ${\eta }_a$ is the Gaussian random field and is fully characterized by the first two moments

$$\begin{array}{c}\hfill \langle {\eta }_a(x,t)\rangle =0,\langle {\eta }_a(x,t){\eta }_b(x^{\prime },t^{\prime })\rangle =2{\alpha }_{ab}{\delta }^d(x-x^{\prime })\delta (t-t^{\prime }).\end{array}$$

(9)

The explicit form of the static action S, which encodes the types of underlined symmetry and strengths of the assorted interactions, depends on the phase transitions or critical points we explore.

Moving on, let us review the paradigmatic models which are at the heart of a tremendous number of research papers on critical dynamical behavior. The simplest one is the model A that encompasses the critical dynamic of a single nonconserved order parameter. The corresponding action is given by

$$S=\frac{1}{2}{(\nabla \varphi )}^2+\frac{\tau }{2}{\varphi }^2+\frac{g}{4!}{\varphi }^4,\alpha =const>0.$$

(10)

The model B is determined by the same static action, Equation (10), and describes relaxation of the conserved field, resulting in the differential operator $\alpha =-\lambda {\nabla }^2$ with a positive kinetic coefficient $\lambda >0$.

The rich collective dynamics of a nonconserved field $\varphi$ and conserved scalar one m are captured by the effective action of the model C

$$S=\frac{1}{2}{(\nabla \varphi )}^2+\frac{g_{\varphi }}{4!}{\varphi }^4+\frac{m^2}{2}+\frac{1}{2}{g}_{m\varphi }m{\varphi }^2,\alpha =\left(\begin{array}{cc}{\lambda }_{\varphi }& 0\\ 0& -{\lambda }_m{\nabla }^2\end{array}\right),$$

(11)

with a positive set of coefficients $\{{\lambda }_{\varphi },{\lambda }_m\}$ and two independent coupling constants $\{g_{\varphi },g_{m\varphi }\}$.

The study of fully conserved dynamic of two fields $\varphi$ and m with the equilibrium action, Equation (11), gives rise to the model D, where the operator $\alpha$ has the form

$$\alpha =\left(\begin{array}{cc}-{\lambda }_{\varphi }{\nabla }^2& 0\\ 0& -{\lambda }_m{\nabla }^2\end{array}\right).$$

(12)

Until now, the above definitions of the stochastic models have not incorporated all the possibilities of the entire equation, Equation (8). One can specify the operator $\beta$ responsible for the inter-mode interaction during the relaxation process. On the one hand, additional couplings $\beta$ lead to a more diverse physical picture. On the other hand, this sufficiently complicates the theoretical analysis, which is reflected in the absence of multi-loop renormalization group calculations, and the available results that some authors led are currently controversial.

The static action of the model F of the complex valued field $\psi$ and the real scalar one m is

$$S=\nabla \psi \nabla {\psi }^{+}+\frac{g_{\psi }}{3!}{\left(\psi {\psi }^{+}\right)}^2+\frac{m^2}{2}+g_{m\psi }m\psi {\psi }^{+},$$

(13)

while

$$\begin{array}{ccc}& & \alpha =\left(\begin{array}{ccc}0& {\lambda }_{\psi }& 0\\ {\lambda }_{\psi }& 0& 0\\ 0& 0& -{\lambda }_m{\nabla }^2\end{array}\right),\beta =\left(\begin{array}{ccc}0& i{v}_1& i{v}_2\psi \\ -i{v}_1& 0& -i{v}_2{\psi }^{+}\\ -i{v}_2{\psi }^{+}& i{v}_2\psi & 0\end{array}\right),\hfill \end{array}$$

(14)

where $v_1,v_2$ are additional real coupling constant. We also can arrive at the model E by setting $g_{m\psi }=0$ and $v_1=0$.

The model G describes relaxation of the nonconserved ${\varphi }_a$ and conserved $m_a$ three-component ($a=1,2,3$) real fields with the action

$$S=\frac{1}{2}{(\nabla \varphi )}^2+\frac{g}{4!}{\left({\varphi }^2\right)}^2+\frac{m^2}{2},$$

(15)

and

$$\alpha =\left(\begin{array}{cc}-{\lambda }_{\varphi }& 0\\ 0& -{\lambda }_m{\nabla }^2\end{array}\right).$$

(16)

For notational clarity, one can concatenate these fields and introduce the six-component “doublet” $\phi$, such that ${\phi }_a={\varphi }_a$ and ${\phi }_{a+3}=m_a$. Then, the $6\times 6$ matrix $\beta$ is determined by its elements

$${\beta }_{ab}=0,{\beta }_{a,b+3}=v{\epsilon }_{abc}{\varphi }_c,{\beta }_{a+3,b+3}=v{\epsilon }_{abc}{m}_c,a,b,c=1,2,3,$$

(17)

with a coupling constant v.

Model H is based on the action

$$S_{\mathrm{st}}=\frac{1}{2}{(\nabla \varphi )}^2+\frac{g}{4!}{\varphi }^4+\frac{{\mathrm{v}}_{\perp }^2}{2},$$

(18)

where ${\mathrm{v}}_{\perp }$ is a three-component divergence-free velocity field and

$$\alpha =\left(\begin{array}{cc}-{\lambda }_{\varphi }{\nabla }^2& 0\\ 0& -{\lambda }_{\mathrm{v}}{\nabla }^2\end{array}\right),\beta =\left(\begin{array}{cc}0& v\nabla \varphi \\ -{[v\nabla \varphi ]}^T& 0\end{array}\right).$$

(19)

with a real coupling constant v, and the superscript T denotes the transpose of an operator.

The variety of non-equilibrium physical phenomena and complex systems and, as a consequence, a set of the respective stochastic field models are beyond the scope of this classification. Among them, one can refer to the most notable and key topics: fully developed turbulence [24], advection–diffusion problems [25], share flows [26], chemical reactions [27,28], self-organized criticality [29], random walks [30], and Kardar–Parisi–Zhang dynamics. Overall, the diagrammatic formulation of stochastic field theory is a bit more unwieldy than that of static models. In addition to the emergence of the integration over internal frequencies, there appear new fields, various types of interaction vertices, and propagators, which bare structure may be more complicated. Thus, the implementation of the renormalization group program is a tricky business. For a long time, it was believed that the calculations in critical dynamics lag behind the gains made in static theories by at least two loops of the perturbation expansion.

One of the first who employed diagrammatic methods in the stochastic field theory was Wyld [31]. To analyze Navier–Stokes turbulence, he put forward a systematic perturbation approach analogous to Feynman diagrams. After a while, through the work of Martin, Siggia, and Rose [32], Janssen [33], and De Dominicis [34] it became clear that the study of a certain class of stochastic equations, in particular Equation (8) is tantamount to that of the effective field models with an action $S_{\mathrm{dyn}}$. Thus, dynamical correlation functions were put in one-to-one correspondence to the Green functions determined by the generating functional

$$G(A,A^{\prime })=\int \prod _{a}D{\varphi }_{a}D{\varphi }_a^{\prime }{e}^{S_{\mathrm{dyn}}+A_a{\varphi }_a+A_a^{\prime }{\varphi }_a^{\prime }},$$

(20)

with the action

$$S_{\mathrm{dyn}}={\varphi }_a^{\prime }{\alpha }_{ab}{\varphi }_b^{\prime }-{\varphi }_a^{\prime }\left({\partial }_t{\varphi }_a-({\alpha }_{ab}+{\beta }_{ab})\frac{\delta S_{\mathrm{st}}}{{\delta }_{{\varphi }_b}}\right).$$

(21)

Eventually, this field theory can be renormalized, similar to the usual quantum field theory. Moreover, the renormalization constants of fields and parameters presented in both the static S and the dynamic $S_{\mathrm{dyn}}$ actions are the same. This is the first simplification of calculations in dynamic models.

For the description of the critical dynamics near the superfluid phase transition, in the vicinity of the $\lambda$-point, the phenomenological F and E models have been proposed [21], but the characteristics of the critical dynamics of this phase transition have not been determined. Two-loop calculations did not lead to an unambiguous choice of a stable fixed point of the corresponding RG equation. Recently, the authors of [35,36] made theoretical progress by approaching this problem in an ab initio manner. The system has been investigated within the formalism of time-dependent Green functions at the finite temperature on the basis of the microscopic action

$$S(\psi ,{\psi }^{+})=\int {\mathrm{d}}^{d}x\underset{C}{\int }\mathrm{d}t\left[{\psi }^{+}\left(i{\partial }_t+\frac{\Delta }{2m_0}+\mu \right)\psi -\frac{g}{2}{\left({\psi }^{+}\psi \right)}^2\right],$$

(22)

for boson particles with the local repulsive interaction of the strength $g>0$. Here $\psi (t,x)$, ${\psi }^{+}(t,x)$ are complex conjugated fields, $m_0$ denotes the particle mass and $\mu$ is the chemical potential in the unit system with $\hslash =1$, $\Delta$ stands for the Laplace operator. The evolution parameter t belongs to the Schwinger–Keldysh contour C, Figure 1, in the complex t-plane [37,38]. The reference points of the contour usually were sent to infinity $t_0\to -\infty$ and $t_f\to \infty$ to enable the use of the Fourier transforms with respect to time.

To construct the perturbation expansion, the fields of the action, Equation (22), were replicated on different parts of the contour by the introduction of new fields ${\psi }_R,{\psi }_R^{+}$,${\psi }_A,$${\psi }_A^{+}$,${\psi }_T,{\psi }_T^{+}$, where the subscripts R, A, and T refer to the corresponding parts of the contour with retarded, advanced, and temperature propagators, respectively. Then, the diagram technique in the model is similar to that in the stochastic dynamics. Three-loop calculation [39] in this model leads to the conclusion that the dynamics near $\lambda$ point is described by the simple two-component A model. As was shown by the analysis of stochastic equations, it is because E and F models for compressible fluids are unstable.

3. Calculation of Dynamic Diagrams

Let us consider for definiteness some basic techniques of computing the analytical expressions corresponding to ultra-violet divergent parts of dynamic diagrams at the example of ones in the simplest model A. Using the the explicit form of the static action, Equation (10), together with the general formulae, Equation (21), one obtains the respective dynamic action

$$S_{\mathrm{dyn}}={\varphi }^{\prime }\alpha {\varphi }^{\prime }-{\varphi }^{\prime }\left({\partial }_t\varphi -\alpha (-\Delta +\tau )\varphi -\frac{\alpha g}{6}{\varphi }^3\right).$$

(23)

The free propagators of the theory are readily established in the $(p,\omega )$ representation and take the form

(24)

The one-loop diagram shown in Figure 2 gives rise to a renormalization of the coupling constant g at the leading order and corresponds (here, for simplicity, we fix the incoming momentum and frequency to be zero) to the following integral

$$\int \frac{{\mathrm{d}}^{d}p\mathrm{d}\omega }{{\left(2\pi \right)}^{d+1}}\frac{1}{-i\omega +\alpha (p^2+\tau )}\frac{2\alpha }{{\omega }^2+{\left[\alpha (p^2+\tau )\right]}^2},$$

(25)

which involves divergences and thus must be suitably regularized, e.g., by applying dimensional regularization $d=4-\epsilon$. As a first step, we need to calculate the frequency integral using the residue theorem. The remaining momentum integral

$$\frac{1}{2\alpha }\int \frac{{\mathrm{d}}^{d}p}{{\left(2\pi \right)}^d}\frac{1}{{(p^2+\tau )}^2},$$

is then amenable to further analytical evaluation based on the above mentioned methods in Section 1. An alternative way of computing loop diagrams, which may be more appropriate in the multi-loop case, is based on the $(p,t)$ representation, where the bare propagators cast in the form

$$\langle \varphi {\varphi }^{\prime }\rangle =\Theta (t-t^{\prime })e^{-\alpha (p^2+\tau )(t-t^{\prime })},\langle \varphi \varphi \rangle =\frac{1}{p^2+\tau }{e}^{-\alpha (p^2+\tau )|t-t^{\prime }|}.$$

(26)

In $\langle \varphi {\varphi }^{\prime }\rangle$ the temporal variable $t^{\prime }$ is the argument of the response field ${\varphi }^{\prime }$. Now, the loop integral of the same diagram, Equation (25), reads as

$$\begin{array}{c}\hfill \int \frac{{\mathrm{d}}^{d}p}{{\left(2\pi \right)}^d}\underset{-\infty }{\overset{\infty }{\int }}\mathrm{d}t\Theta \left(t\right)\frac{1}{(p^2+\tau )}{e}^{-2\alpha (p^2+\tau )t}=\frac{1}{2\alpha }\int \frac{{\mathrm{d}}^{d}p}{{\left(2\pi \right)}^d}\frac{1}{{(p^2+\tau )}^2},\end{array}$$

(27)

and can again be straightforwardly evaluated in a standard manner.

However, the maximal progress in the three-loops pure analytic calculation in the dynamical A model was achieved in [40]. There, an approach with Fourier and inverse Fourier transformations between frequency—momentum and time—coordinate representations was used.

3.1. Method of Calculation

Let us describe the calculation method [40] in more detail. The origin of the approach lies in the possibility of performing in a closed form the Fourier transform of the function

$$\frac{\theta (t-t^{\prime })}{{(t-t^{\prime })}^a}exp\left[-p^2(t-t^{\prime })b\right],$$

(28)

with respect to either momentum p or temporal variable t to arrive at the function of $x-x^{\prime }$ or frequency $\omega$, respectively. Comparison of these results leads to the following explicit expression for the Fourier transform from the $(x,t)$ representation to the $(p,\omega )$ one. In case $Re\left(b\right)>0$, the result for arbitrary values of a is

(29)

where d is the dimension of the (coordinate) space.

The main advantage of the approach considered is that the loop diagrams with lines in the form of r.h.s. Equation (29) can be simply calculated, and we need not calculate corespondent integrals in momenta and frequency. After the inverse Fourier transform of all lines into the $(x,t)$ representation, expressions in the form of l.h.s. Equation (29) must be multiplied.

3.2. Convolution of Lines

The loop-like diagrams can be quite easily calculated by the approach considered, but the convolution of subgraphs similar to Equation (7) is a rather tricky problem. Nevertheless, it should be exploited if the result leads to some new loop diagram or subgraph. Using the profitable properties, Equation (29), the propagator lines are cast into the $(p,\omega )$ representation in which the convolution of two arbitrary lines is turned into a simple product

or

depending on the temporal arguments of the $\theta$-function in the propagators. The Feynman identity, Equation (4),

$$\frac{1}{B_1^{a_1}{B}_2^{a_2}}=\frac{\Gamma (a_1+a_2)}{\Gamma \left(a_1\right)\Gamma \left(a_2\right)}\underset{0}{\overset{1}{\int }}\mathrm{d}{z}_1\underset{0}{\overset{1}{\int }}\mathrm{d}{z}_2\frac{\delta \left(z_1+z_2-1\right)z_1^{a_1-1}{z}_2^{a_2-1}}{{\left[B_1{z}_1+B_2{z}_2\right]}^{a_1+a_2}}$$

produces

In the second expression, the integral over $z_1$ can be split into a sum of two. The first one is over $0<z_1<1/2$, and the second is over $1/2<z_1<1$. These parts contribute to different $\theta$-functions in the $(x,t)$ representation. It is convenient to rescale the $z_1$ variables to obtain integrals in new variables ${\overline{z}}_1$ in the limits $0<{\overline{z}}_1<1$. After the inverse Fourier transform, the convolution integrals over $x_1,t_1$ may be expressed as

(30)

The derived expressions are fairly voluminous, but, despite this, they can be evaluated numerically. Note that integrals in some z type variables can be calculated analytically.

3.3. Example

Here, we demonstrate the calculation method by evaluating the two loops diagram, Figure 3a, in the massless model A. The first step is to represent all propagators, Equation (26), in the form Equation (28). Note the arrow lines already have the necessary form. While the plain line representing the bare propagator $\langle \varphi \varphi \rangle$ can be decomposed into a retarded and an advanced component

$$\langle \varphi (p,t)\varphi (-p,t^{\prime })\rangle =\frac{e^{-\alpha p^2|t-t^{\prime }|}}{p^2}=\frac{e^{-\alpha p^2(t-t^{\prime })}}{p^2}\theta (t-t^{\prime })+\frac{e^{+\alpha p^2(t-t^{\prime })}}{p^2}\theta (t^{\prime }-t).$$

(31)

The next step is to transform both components appropriately. Since they have a common structure, we will describe transformation only for the first one, while the second can be treated analogously

$$\begin{array}{cc}\hfill \frac{e^{-\alpha p^2(t-t^{\prime })}}{p^2}\theta (t-t^{\prime })=\alpha \theta (t-t^{\prime })(t-t^{\prime })& \underset{1}{\overset{\infty }{\int }}\mathrm{d}s{e}^{-\alpha p^2(t-t^{\prime })s}\stackrel{\mathcal{F}}{\Rightarrow }\hfill \\ \hfill & \frac{{\alpha }^{1-d/2}}{{\left(4\pi \right)}^{d/2}}\underset{0}{\overset{1}{\int }}\frac{\mathrm{d}s}{s^{2-d/2}}\frac{\theta (t-t^{\prime })}{{(t-t^{\prime })}^{d/2-1}}exp\left[-\frac{{(x-x^{\prime })}^2}{t-t^{\prime }}\frac{s}{4\alpha }\right].\hfill \end{array}$$

(32)

Here, we got rid of the denominator $p^2$ by introducing an integration over an auxiliary variable s. Now, being Gaussian in the momentum p, the propagator can be brought to the $(x,t)$ representation by employing the momentum Fourier transform. Under the s-integral we recognize the function as that of the form Equation (28) with the parameters $a=d/2-1$ and $b=s/\left(4\alpha \right)$. Performing the same steps as before, one can reduce all propagators in the diagram to a set of propagators of the form Equation (28) in the $(x,t)$ representation.

The diagram, Figure 3a, produces a logarithmic ultra-violet divergence. To calculate this contribution to renormalization constants, we put the momentum and frequency of the external line in the upper right corner equal to zero. Then we have to calculate the amputated graph shown in Figure 3b. To do this, we multiply the lines in the $(x,t)$ presentation in the right loop labeled by 2 and 3. Then using Equation (30) we calculate its convolution with the upper $\langle \varphi {\varphi }^{\prime }\rangle$ line marked by 1. The result of the convolution is a sum of new lines obtained. Then every one of these lines is multiplied by the lower line $\langle \varphi {\varphi }^{\prime }\rangle$ marked by 4. After the transform to the $(p,\omega )$ representation with the use of Equation (29), we obtain the result in the form of integrals over two variables similar to that over s in Equation (32) and over z in Equation (30). The $\Gamma$ function, due to the Fourier transform, contains the leading pole in $\epsilon$. The second pole in $\epsilon$ can be simply extracted from the integral over z. The remaining integrals are readily calculated analytically or numerically. In some cases, this approach can be used to obtain closed-form expressions for the diagrams.

Note that the four loop results in the model A have been achieved using another scheme [41]. Unlike other approaches, a massive theory was considered. This scheme was targeted to pure numerical computations of diagrams at zero momenta and frequencies after the so-called $R^{\prime }$-operation [8]. Nevertheless, the five loops results were achieved in [42] using the method of sector decomposition [14] in the massless theory. To illustrate this approach, let us overwrite the calculation of the same diagram, Figure 3b. In the $(p,t)$ representation, one obtains two denominators after integration over temporal variables t of the vertexes between lines numbered $(2,1,3)$ and $(2,4,3)$. These denominators are quadratic in momenta. The other two denominators one has due to lines 2 and 3. Using Equations (4) and (3), one can evaluate all integrals in momenta and obtain the integrals over four Feynman parameters $z_1,\dots ,z_4$. Then we split the integration domain into sectors: $z_1\ge \{z_2,z_3,z_4\}$, $z_2\ge \{z_1,z_3,z_4\}$, etc. In each sector, the smaller variables are scaled, e.g., $\{z_2,z_3,z_4\}=\{{\overline{z}}_2{z}_1,{\overline{z}}_3{z}_1,{\overline{z}}_4{z}_1\}$, $\{z_1,z_3,z_4\}=\{{\overline{z}}_1{z}_2,{\overline{z}}_3{z}_2,{\overline{z}}_4{z}_2\}$, etc., to restore the unit upper limits of integration. Then the integral in the largest $z_i$ can be calculated analytically, and the procedure is repeated until obtaining the convergent integral in the remaining parameters. Numerical calculations are needed for the remaining integrals obtained.

4. Simplification for Critical Dynamic Models: Diagram Reduction

As one can see, the calculation of the dynamics models faces a significant amount of problems in comparison to static models. However, several approaches can help to reduce the complexity of the considered dynamics model. First, the renormalization constants of the static fields and parameters in the chosen renormalization scheme are combined with those calculated in the static theory, where the calculations are much simpler [8]. However, the main interest of the dynamics model is in renormalization constants that correspond to the dynamics field and parameters.

The complexity of the Feynman diagrams and integrands corresponding to them increases the computational time drastically compared to the same order static models. First of all, the existence of different kinds of fields increases the number of possible propagators and, as a result, the number of diagrams in each order. For example, the statical ${\varphi }^4$ model in the fifth order contains only 11 two-point diagrams. However, the simplest dynamical model A, which has only one more extra field, already has 95 dynamics diagrams. Moreover, another difference from the static case, as we showed in Section 3, is additional integration by time (or frequency). For using the same calculation technique as in a static model (in particular, the sector decomposition method), it is necessary to take the time variable integration in a diagram that generates the so-called time versions. The time version has a nontrivial structure in comparison to the static version of the diagram, which requires a slight modification of the sector decomposition. To go back to the model A, the 95 dynamics diagrams mentioned lead to 1025 corresponding time versions. This fact shows the necessity of advanced techniques for reducing the dynamics diagrams.

The problem with a nontrivial denominator was solved by the adaption of the sector decomposition method to the diagrams of critical dynamics. And the problem with the growing number of time versions was largely overcome due to the developed diagram reduction scheme. The proposed approach allows one to automate the calculation process as much as possible and use most of the tools from the static case. Below, we illustrate these procedures based on the model A.

Time versions. As an example, we consider three-loop ${\varphi }^{\prime }{\varphi }^{\prime }$-diagram. Due to the translational invariance of the propagators (see Equation (26)) over time, the time at one of the vertices (for example, $t_0$ for (33) can be considered to be fixed. The area of integration over the remaining times $t_1$ and $t_2$ is limited by the presence of the function $\theta (t_0-t_1)$ in the propagator $\langle \varphi \left(t_0\right){\varphi }^{\prime }\left(t_1\right)\rangle$. This area can be divided into the following regions (time versions): $(t_0>t_1>t_2)$, $(t_0>t_2>t_1)$, and $(t_2>t_0>t_1)$. In accordance with these inequalities, the relative position of the vertices in the diagram (in decreasing order from left to right) was fixed. We represent this diagram as a sum of three contributions:

(33)

Integration over relative time in each of the time intervals is easily performed. The result can be represented in the form of diagrams containing only integration over momentum:

The lines without any hatch with momentum k correspond to the “static propagator” $1/E_k$, where $E_k=k^2+\tau$. The lines with a hatch correspond to 1. Dashed lines correspond to a factor in the denominator, equal to the sum of the “energies” $E_k$ lines that cross dashed lines. In particular, for the third time version, we have

Diagram Reduction Scheme The integration over time significantly decreases the number of considering integrands. However, the diagram reduction scheme that was proposed in Refs. [41,42], can solve this issue. The idea is based on grouping the time versions according to certain rules, which leads to the construction of the modified diagrams, the number of which is much smaller than the original.

Using diagram reduction, it is easy to show that the combination of the dynamics diagram is equal to some of one static diagrams with the same topology. It corresponds to the correlation between dynamic and static renormalization constants. To illustrate it, let us consider the sum of two-time versions for one dynamic diagram; see Equation (34). The time version here implies that the times in the vertexes in the picture are ordered due to the vertical cross-sections. The coefficients in front of them correspond to the symmetric factor.

(34)

Achieved progress. Since the 1970s, several works [43,44,45,46,47] used a grouping of the diagrams to simplify the answer. However, none of this procedure was generalized. Thus, the automatization of the procedure was not even remotely possible. For the dynamical model A, in recent years, the generic rules of diagram grouping—the diagram reduction scheme—were proposed in [41,42]. That can help reduce the number of considered diagrams and simplify corresponding integrands. As a result, it leads to a significant reduction in computational time. For model A, the amount of reduced diagrams in the fifth order becomes 201, which is five times less than without reduction. The generalization of the diagram–reduction method allowed calculating the critical index z until th fifth order:

$$\begin{array}{cc}\hfill & z\left(n\right)=2+{\epsilon }^2\frac{(2+n)}{2{(8+n)}^2}(-1+6ln(4/3))+{\epsilon }^3\frac{(2+n)}{8{(8+n)}^4}(162.461908017+31.9023677877n\hfill \\ \hfill & -1.27351796946n^2)+{\epsilon }^4\frac{(2+n)}{32{(8+n)}^6}(-23752.4\left(16\right)-6929.0\left(8\right)n-770.28\left(13\right)n^2\hfill \\ \hfill & -170.470\left(6\right)n^3-4.24423766241n^4)+{\epsilon }^5\frac{(2+n)}{128{(8+n)}^8}(1.986\left(25\right){10}^7+1.038\left(23\right){10}^{7}n\hfill \\ \hfill & +2.15\left(8\right){10}^6{n}^2+1.92\left(13\right){10}^5{n}^3-3.8\left(10\right){10}^3{n}^4-626\left(29\right)n^5+7.500527164701581n^6).\hfill \end{array}$$

(35)

The achieved progress was expended on more complicated models. For example, in the E-model, the leading order turned out to be insufficient for determining the IR of a stable fixed point; therefore, in the papers [45,48] (due to disagreements between the proposed results), the second-order $\epsilon$-expansion was calculated. Later [49] recalculate the second order. However, it also turned out to need to be more sufficient to answer the question about the type of critical dynamic behavior. The third order of perturbation theory can change the IR stability of the fixed points. In work [50], the recipe for a reduction diagram was proposed, which can be easily implemented for higher orders too.

5. Resummations of the Divergent Series

Quantum-field perturbation series usually have zero radii of convergence. The $\epsilon$ series considered in this article are not included. They cannot be used for physical $\epsilon =1$.

5.1. Resummations in Static Models

Using the instanton approach, the technique of the resummations of the convergent quantum-field perturbation series was developed in the base of $\epsilon$-expansion in static equilibrium statistical models. Let us revisit the main features of the instanton approach in the example of the static theory with the action

$$S=\frac{1}{2}{(\nabla \varphi )}^2+\frac{g}{4!}{\varphi }^4,$$

(36)

where $\varphi \left(x\right)$ is the basic field, g is the coupling constant, and all necessary integrations in coordinates (and times in dynamic models) and summation in indices of fields and partial derivatives are implied. We will use the notation $X^{\left[N\right]}$ for the N-th order contribution to the perturbation expansion in a small parameter of an arbitrary quantity X. The large-order asymptotics for the expansion in g of the k-point correlation function in a model, Equation (36), can be determined using the expression

$$G_k^{\left[N\right]}(x_1\dots x_k)=\int \mathcal{D}\varphi \varphi \left(x_1\right)\dots \varphi \left(x_k\right)\frac{1}{2\pi i}\oint \frac{\mathrm{d}g}{g^{N+1}}{e}^{-S}$$

(37)

based on the Cauchy residue theorem. The instanton approach is, in fact, the steepest-descent method for approximating the path integral over $\varphi$ and g for large order $N\to \infty$. The respective saddle points satisfy the stationarity equations $\delta S/\delta \varphi =0$ and $\partial S/\partial g=-N/g$, which are explicitly rewritten for the action, Equation (36), in the form of the partial differential equation with the normalization condition

$$-\Delta {\varphi }_s+\frac{g_{\mathrm{st}}}{6}{\varphi }_s^3=0,\frac{g_{\mathrm{st}}}{4!}\int {\mathrm{d}}^{d}x{\varphi }_s^4=-N.$$

(38)

The set of spherically symmetric solutions—instantons—has been found by the author [51] in space dimension $d=4$

$${\varphi }_s\left(x\right)=\frac{\sqrt{3N}}{\pi }\frac{y}{{(x-x_0)}^2+y^2},g_{\mathrm{st}}=-\frac{16{\pi }^2}{N},$$

(39)

where an arbitrary vector $x_0$ is called the center of the instanton, and a scalar magnitude y refers to the instanton radius. The onset of arbitrariness in classical solutions is rooted in the translation and scale invariance of the original action, Equation (36). All saddle points with different values of $x_0$ and y yield the same contribution to the action in the path integral, Equation (37), and, therefore, this arbitrariness has to be properly taken into consideration as in the case of gauge field theories, which are successively treated via the celebrated Faddeev–Popov trick. To this end, the following integral representation of unity

$$\begin{array}{c}1=\int {\mathrm{d}}^d{x}_0\underset{-\infty }{\overset{+\infty }{\int }}\mathrm{d}ln{y}^2\delta \left[-\frac{g}{4!}\int {\mathrm{d}}^{d}x{\varphi }^4\left(x\right)ln{\left(\frac{x-x_0}{y}\right)}^2\right]\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\times {\delta }^d\left[-\frac{g}{4!}\int {\mathrm{d}}^{d}x{\varphi }^4\left(x\right)(x-x_0)\right]{\left[-\frac{g}{4!}\int {\mathrm{d}}^{d}x{\varphi }^4\left(x\right)\right]}^{d+1},\hfill \end{array}$$

(40)

is inserted into the path integral, Equation (37). The detailed computations of the large-order asymptotic for diverse Green functions or the renormalization group ones are technically cumbersome and, inter alia, associated with the appropriate elimination of ultra-violet divergences [51,52,53], but the generic form of the N-th expansion coefficient for an arbitrary quantity F can be presented in the form

$$F^{\left[N\right]}=c{a}^N{N}^{b}N!\left(1+\mathcal{O}\left(N^{-1}\right)\right),$$

(41)

where the constant a is solely determined by the value of the field action at the saddle point; the constant b and the amplitude c depend on the certain quantity F wherein c may be a spatial function.

As the series is divergent, Equation (41), to sum a large number of terms for $\epsilon$-expansion, Equation (35), is senseless. Then one has to use Borel resummation [6,7,54]. Let us outline the main features of Borel resummation. First, let Q be a function determined as a series in the parameter $\epsilon$:

$$Q\left(\epsilon \right)=\sum _{k\ge 0}{Q}_k{\epsilon }^k,$$

(42)

and its high-order asymptotics (HOA) behavior specified by Equation (41). The Borel transform of the series, Equation (42), is defined by the formula

$$Q\left(\epsilon \right)={\int }_0^{+\infty }\mathrm{d}t{t}^b{e}^{-t}B\left(\epsilon t\right),B\left(x\right)=\sum _{k\ge 0}{B}_k{x}^k,B_k=\frac{Q_k}{\Gamma (k+b+1)},$$

(43)

where b is an arbitrary parameter. The knowledge of HOA, Equation (41), and some suppositions concerning analytical properties of $B\left(x\right)$ function (see [6,7,54]) allow us to resum the series Equation (42) in compliance with expressions Equation (43) and to obtain more accurate value of $Q\left(\epsilon \right)$ value. HOA Equation (41) implies the convergence of the series Equation (43) in the circle $\left|x\right|<1/a$. The nearest singularities are situated on the real negative semi-axis in the point $(-1/a,0)$. Note the integration in Equation (43) is performed over the whole positive axis $t\in [0,+\infty )$. The integration contour crosses the convergence circle boundary of Equation (43) in the point $1/a$. The analytic continuation of Equation (43) outside the convergence circle can be constructed by conformal mapping of the complex plane or by Padé-approximation [7,54]. Conformal maps usually have the form [6,54]

$$u\left(x\right)=\frac{\sqrt{1+ax}-1}{\sqrt{1+ax}+1}⟺x\left(u\right)=\frac{4u}{a{(u-1)}^2}.$$

(44)

The integration contour here belongs to a convergence circle; the point $x=-1/a$ is mapped to $u=-1$, and the infinity point goes to $u=1$. A small value of x in Equation (44) corresponds to small u. Then, the series, Equation (43), can be rewritten in terms of u variable

$$B\left(x\right)=\sum _{n\ge 0}{x}^n{B}_n⟹B\left(u\right)=\sum _{n\ge 0}{u}^n{U}_n,$$

where

$$U_0=B_0,U_n=\sum _{m=1}^n{B}_m{(4/a)}^m{C}_{n+m-1}^{n-m},n\ge 1.$$

(45)

Correspondingly, the conformal Borel transform of the function Q has the form

$$Q\left(\epsilon \right)=\sum _{k\ge 0}{U}_k\underset{0}{\overset{+\infty }{\int }}\mathrm{d}t{t}^b{e}^{-t}{\left(u\left(\epsilon t\right)\right)}^k.$$

(46)

Additionally, the use of the Padé-Borel-Leroy (PBL) resummation for Equation (43) are based on the Padé approximants

$$P_M^L\left(x\right)=\frac{{\sum }_{i=0}^L{p}_i{x}^i}{1+{\sum }_{j=1}^M{q}_j{x}^j},M\ge 1L+M=N.$$

(47)

Coefficients $\left\{p,q\right\}$ are fixed by the matching condition for the coefficients of the series Equation (43) and coefficients for the Taylor series of the function Equation (47) at $x\to 0$:

$$\begin{array}{c}\hfill B_0=p_0,B_n=\frac{1}{n!}{\left.\frac{{\mathrm{d}}^n}{\mathrm{d}{x}^n}{P}_M^L\left(x\right)\right|}_{x=0},n=1..N.\end{array}$$

(48)

Solving this system with respect to $\left\{p,q\right\}$ variables, one obtains

$$\begin{array}{ccc}\hfill p_s& =& p_s(B_0,B_1,\dots ,B_N),s=1,\dots ,L,\hfill \\ \hfill q_{s^{\prime }}& =& q_{s^{\prime }}(B_0,B_1,\dots ,B_N),s^{\prime }=1,\dots ,M.\hfill \end{array}$$

(49)

The Borel-transformed expression is defined now as

$$Q\left(\epsilon \right)=\underset{0}{\overset{+\infty }{\int }}\mathrm{d}t{t}^b{e}^{-t}{P}_M^L\left(\epsilon t\right).$$

(50)

This technique must be considered as less accurate because it does not use the whole information known about HOA, Equation (41).

5.2. Resummation in Dynamics

As mentioned above, the dynamic action for the model A has the form Equation (23) and in the massless case $\tau =0$ reads as

$$S_{\mathrm{dyn}}={\varphi }^{\prime }\alpha {\varphi }^{\prime }-{\varphi }^{\prime }\left({\partial }_t\varphi -\alpha \left(-\Delta \varphi +\frac{g}{3!}{\varphi }^3\right)\right).$$

(51)

Hence, the stationarity equations have the form

$$\begin{array}{cc}\hfill & \frac{\delta S_{\mathrm{dyn}}}{\delta {\varphi }^{\prime }}=0\Rightarrow 2\alpha {\varphi }^{\prime }-{\partial }_t\varphi +\alpha \left(-\Delta \varphi +\frac{g}{6}{\varphi }^3\right)=0,\hfill \\ \hfill & \frac{\delta S_{\mathrm{dyn}}}{\delta \varphi }=0\Rightarrow {\partial }_t{\varphi }^{\prime }+\alpha \left(-\Delta {\varphi }^{\prime }+\frac{g}{2}{\varphi }^2{\varphi }^{\prime }\right)=0.\hfill \end{array}$$

This is a system of partial nonlinear differential equations. The principal question that should be at the back of our minds “Does the nonzero solution exist?”. The key to the answer to this question lies in a scaling transformation of the derived equations. Let us consider the auxiliary identities

$$\int {\mathrm{d}}^{d}x\mathrm{d}t{\varphi }^{\prime }\frac{\delta S_{\mathrm{dyn}}}{\delta {\varphi }^{\prime }}=0,\int {\mathrm{d}}^{d}x\mathrm{d}t\varphi \frac{\delta S_{\mathrm{dyn}}}{\delta \varphi }=0,$$

(52)

and change variables according to $\varphi (t,x)\to \varphi ({\lambda }^{\prime }t,\lambda x)$, where the scale parameters ${\lambda }^{\prime },\lambda$ take arbitrary real values. Then $S_{\mathrm{dyn}}=S_{\mathrm{dyn}}({\lambda }^{\prime },\lambda )$ and we have two additional scaling equations $\lambda {\partial }_{\lambda }{S}_{\mathrm{dyn}}{|}_{\lambda ,{\lambda }^{\prime }=1}=0$, ${\lambda }^{\prime }{\partial }_{{\lambda }^{\prime }}{S}_{\mathrm{dyn}}{|}_{\lambda ,{\lambda }^{\prime }=1}=0$, which follow from Equation (52). For vanishing at infinity, namely $\left|t\right|\to \infty ,\left|x\right|\to \infty$, fields, we obtain a system of linear equations for the functionals

$$F_1\equiv \int {\mathrm{d}}^{d}x\mathrm{d}t{\left({\varphi }^{\prime }\right)}^2,F_2\equiv \int {\mathrm{d}}^{d}x\mathrm{d}t{\varphi }^{\prime }{\partial }_t\varphi ,F_3\equiv \int {\mathrm{d}}^{d}x\mathrm{d}t{\varphi }^{\prime }\Delta \varphi ,F_4\equiv \int {\mathrm{d}}^{d}x\mathrm{d}t{\varphi }^{\prime }{\varphi }^3,$$

contributing to the dynamic action in the form

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & d\alpha F_1-d{F}_2+(d-2)\alpha F_3-\frac{d\alpha g}{3!}{F}_4=0,\hfill \\ \hfill & 2\alpha F_1-F_2+\alpha F_3-\frac{\alpha g}{3!}{F}_4=0,\hfill \\ \hfill & \alpha F_1+\alpha F_3-\frac{\alpha g}{3!}{F}_4=0,\hfill \\ \hfill & F_2-\alpha F_3+\frac{\alpha g}{2}{F}_4=0.\hfill \end{array}\end{array}$$

(53)

following from the stationarity and scaling equations. It can be proved that these equations have no nonzero solution. Thus, a nontrivial solution of the stationarity equations in dynamic models does not exist in the class of functions decreasing at $t\to \pm \infty$.

5.3. Instanton Analysis of Models A-H

The boundary contribution becomes essential when there is no stationarity point in the steepest-descent calculation of a numeric integral. It was stated in [55,56] that within the instanton analysis for the near-equilibrium dynamic models, the necessary boundary contribution is produced by functions that do not decrease in the limit $t\to +\infty$. Let us illustrate this statement by considering the steepest-descent calculation of the parametric and path integral

$$\begin{array}{c}\hfill {\displaystyle \int \mathcal{D}\varphi \mathcal{D}{\varphi }^{\prime }\varphi \dots {\varphi }^{\prime }\frac{1}{2\pi i}\oint \frac{\mathrm{d}g}{g}{e}^{{\displaystyle -S_{\mathrm{dyn}}-Nlng}}}.\end{array}$$

The stationarity equations are

$$\begin{array}{c}\hfill -\frac{\partial {\varphi }^{\prime }}{\partial t}+{\varphi }^{\prime }\left(\alpha +\beta \right)\frac{{\delta }^2{S}_{\mathrm{dyn}}}{\delta \varphi \delta \varphi }+{\varphi }^{\prime }\frac{\delta \beta }{\delta \varphi }\frac{\delta S_{\mathrm{dyn}}}{\delta \varphi }=0,\end{array}$$

(54)

$$\begin{array}{c}\hfill -2\alpha {\varphi }^{\prime }+\frac{\partial \varphi }{\partial t}+\left(\alpha +\beta \right)\frac{\delta S_{\mathrm{dyn}}}{\delta \varphi }=0.\end{array}$$

(55)

At the finite time interval $t\in [t_0,T]$ an additional condition appears:

$${\varphi }^{\prime }(T,x)=0.$$

(56)

It can be noticed that the dynamic instanton ${\varphi }_d$ is, in particular, a nontrivial solution of the equation

$$\begin{array}{c}\hfill -\frac{\partial \varphi }{\partial t}+\left(\alpha -\beta \right)\frac{\delta S}{\delta \varphi }=0.\end{array}$$

(57)

This equation can be considered instead of Equations (54) and (55). The stationarity equation Equation (57) has two rather obvious time-independent solutions: $\varphi =0$ and the static instanton ${\varphi }_{\mathrm{st}}$. It was shown that the dynamic instanton behaves as

$$\underset{t\to -\infty }{lim}{\varphi }_d=0,\underset{t\to \infty }{lim}{\varphi }_d={\varphi }_{\mathrm{st}},$$

consistent with the Gibbsian limit of the dynamic model. It was proved [55,56] that the dynamic action on the dynamic instanton solution asymptotically coincides with the static action on the static instanton:

$$S_{\mathrm{dyn}}\left(\varphi ,{\varphi }^{\prime }\right)=S\left({\varphi }_{\mathrm{st}}\right),$$

and equation ${\partial }_g{S}_{\mathrm{dyn}}=N$ leads to $g_{\mathrm{dyn}}=g_{\mathrm{st}}$. Thus, it was shown that the asymptotic properties of the dynamic model at the leading order in N are determined by the static instanton solution, which leads to the factorial growth of the large-order contributions as in the static case. Therefore, the large-order behavior of an arbitrary quantity F (correlation or response function or critical exponent) may be expressed by the general structure, Equation (41). The most essential for resummation schemes constants a in Equation (41) have been determined in [55,56] for all models of critical dynamics listed in Section 2. The exponent b in the $\epsilon$-expansion contribution to the dynamic index z in the $O\left(n\right)$ symmetric dynamic theories with Gibbsian static limits was determined as

$$b=3+\frac{n}{2}.$$

Properties of the response and dynamic parts of the correlation functions were also discussed in [55,56], and findings obtained there were used in [57] to compare results of different resummation schemes for A model z index in four orders of $\epsilon$-expansion. It was obtained that the results of different approaches are inconsistent. Namely, the intervals of estimated possible errors are not intersected. The only inclusions are the results of conform—the Borel transform and, in a sense, the Padé–Borel one based on the known form of the Equation (41) Borel image [55,56]. It proves the usefulness of all additional information in the resummation procedure. Using HOA information, one obtains the most reliable resummation results concerning the value of critical indexes and errors estimated.

The results of four-loop z index resummations for model A obtained in [57] are presented in

Table 1. Critical exponent z for $O\left(n\right)$ symmetric model A.

$\mathit{d}=3$	$\mathit{n}=1$	$\mathit{n}=2$	$\mathit{n}=3$
$P_{1,-1/a}^4$	$2.023\pm 0.006$	$2.024\pm 0.005$	$2.024\pm 0.005$
CB	$2.013\begin{array}{c}+0.011\hfill \\ -0.0\hfill \end{array}$	$2.014\begin{array}{c}+0.011\hfill \\ -0.0\hfill \end{array}$	$2.014\begin{array}{c}+0.011\hfill \\ -0.0\hfill \end{array}$

.

Here, $P_{1,-1/a}^4$ implies Padé–Borel transform with poles of Borel image fixed in $-1/a$ Equation (41) point (CB—the conformal Borel transform). Five-loop calculation in the model A of critical dynamics leads to modified CB resummation result $z=2.02368$ for $d=3$, $n=1$ [42]. The estimations of index z for different n and $d=2,3$ can be found there also.

6. Conclusions

This review was devoted to modern methods for calculating critical dynamic diagrams and the results obtained with their help. For a long period of time, no noticeable progress has been observed in this area, since the calculations in dynamics are rather cumbersome compared to the static model. However, the developed modern methods and the automation of their processes, which we have described in this article, help to study significant advances in this area.

Since the calculation of the dynamic model is based on the statistical model, we started with well-known approaches for static models, such as the Alpha representation, Feynman representation, sector decomposition method, Kompaniets–Panzer method, and alternating transitions to coordinate and momentum representations. Then, we extended them to a dynamic case and showed modern techniques for simplifying calculations of dynamic diagrams based on their connection with static. We also explain the Borel resummation procedure for divergent perturbation series in dynamical models and described the necessary basis for resummation.

The advantage of these approaches is that they can be automated and easily implemented. Furthermore, these methods have already been successfully used to promote higher orders of perturbation theory. We hope that the given overview will be useful for further research on the critical dynamics and behavior of stochastic scaling.