Symmetry Breaking in Stochastic Dynamics and Turbulence

Abstract

Symmetries play paramount roles in dynamics of physical systems. All theories of quantum physics and microworld including the fundamental Standard Model are constructed on the basis of symmetry principles. In classical physics, the importance and weight of these principles are the same as in quantum physics: dynamics of complex nonlinear statistical systems is straightforwardly dictated by their symmetry or its breaking, as we demonstrate on the example of developed (magneto)hydrodynamic turbulence and the related theoretical models. To simplify the problem, unbounded models are commonly used. However, turbulence is a mesoscopic phenomenon and the size of the system must be taken into account. It turns out that influence of outer length of turbulence is significant and can lead to intermittency. More precisely, we analyze the connection of phenomena such as behavior of statistical correlations of observable quantities, anomalous scaling, and generation of magnetic field by hydrodynamic fluctuations with symmetries such as Galilean symmetry, isotropy, spatial parity and their violation and finite size of the system.

1. Introduction

The success of physics is to a large extent dictated by its enormous predictive power describing many natural phenomena. Ranging from microscopic distances probed at colliders facilities up to macroscopic scales observed through sophisticated astronomical devices, physics develops theories and models that describe reality to a very high precision. Once one understands basic principles, one could not stop wonder about the incredible efficiency with which fundamental physical laws are constructed. To a large extent, guiding principles in physics are based on a correct recognition of underlying symmetries.

From a historical point of view, the first nontrivial symmetry found was Galilei’s discovery of equivalence of inertial frames whose direct consequence is momentum conservation. A further development in theoretical physics, most notably utilized by E. Noether in her work, uncovers the fundamental role in classical physics played by symmetries. Later on, it was realized that many physical theories can be based on a proper identification of symmetries. The prototypical and most successful example is how the Lorentz covariance of particle physics, underlying principles of quantum mechanics and local gauge symmetry, restrict the permissible form of theory to such an extent that every such attempt results in a kind of quantum field theory [1]. From a modern perspective, this is related to the observation that most quantum field theory models should be interpreted as kind of effective field models that describe physics sufficiently up to a certain energy scale. An input from experiments is still needed to impose restrictions on theory, for example, in mirror symmetry, time reversal and charge conjugation play a fundamental role in the formulation of the standard model [2]. Indications of whether a given interaction is accompanied by some symmetry are inferred by the experiment and not from theoretical reasoning alone.

Symmetry considerations are not restricted to particle physics only. They are important in other research areas as well. In particular, our aim in this article is to describe them in the context of classical physics related mainly to applications in fluid mechanics. To set the context and provide the basic framework for theoretical considerations and relations with other branches of physics, we discuss the underlying ideas more broadly.

The problems considered in this paper belong to statistical physics, which forms a cornerstone of the modern science. Since its foundation as a scientific discipline in the works of Gibbs and Boltzmann, it has evolved to a great depth both in scope and rigorousness. At present, methods primarily devoted to the study of physical systems are applied in such diverse scientific fields as chemistry, biology, economics, sociology and computer science. Such a success might be explained by the generality of the fundamental laws of statistical physics and genuine appearance of systems with great resemblance to the models studied in physics. In general, they could be characterized by a large number of entities (atoms, spins and so on) that interact with each other. It is important to realize how large this number is. For instance, merely one mole of ordinary matter under normal circumstances contains as many as ${10}^{23}$ atoms or molecules. This is an incredibly large value beyond human ability to imagine. The rigorous treatment of such a system in terms of particle dynamics (via classical Euler–Lagrange equations or Schrödinger equation in the quantum case) is apparently meaningless. Nonetheless, it is an experimental fact that under stationary boundary conditions all closed macroscopic systems tend to evolve to the equilibrium state that is characterized by a constant value of its macroscopic (coarse-grained) characteristics, e.g., temperature, volume, magnetization, etc. At mesoscopic scales of space and time, this tendency remains true, but the thermodynamic parameters are slowly varying functions of position and time. In a more precise sense, this property of nature is formulated in the second law of thermodynamics. In an isolated system, it also identifies an equilibrium with the most disordered (equivalent to the most probable) state under given external conditions.

In many instances encountered in equilibrium physics, it is possible to use an approximation in which interactions or fluctuations between microscopic constituents of the model can be neglected or treated as a small perturbation to the ideal situation of noninteracting particles. The ideal gas and van der Waals model are famous examples. Note that in the latter case attractive interaction between gas molecules are effectively taken into account by the corresponding virial term [3].

When matter is more dense, interactions between neighboring particles cannot be neglected. The fundamental property of additivity of energy and entropy [3] in the equilibrium is maintained by the assumption that the system may be considered comprising of noninteracting subsystems (“mesoscopic” elements of matter). It should be noted that, for the very possibility of considering noninteracting particles or subsystems as the ideal situation, particles are assumed to have short-range interactions. At the microscopic level, this corresponds to electrically neutral molecules or clusters of molecules. Problems related to formation and structure of these quantities [4] are beyond the scope of this article.

However, there also exist situations, in which neglecting fluctuations is not appropriate at all. The theory of critical phenomena [5,6], which deals with the second-order phase transitions in macroscopic systems, is a well-known representative. It is observed, e.g., in liquid–vapor transition, $\lambda$—transition in superfluid helium or various magnetic transitions between paramagnetic and ferromagnetic phases. A characteristic feature of these transitions is the appearance of strong fluctuations and correlations between underlying constituents (atoms or spins, subsystems comprising thereof). A parameter used for a quantitative description of correlations is the correlation length. Broadly speaking, it represents the average distance to which atoms or clusters “feel” each other or behave cooperatively. In equilibrium situations, it is of atomic order, which explains why a dilute atomic system usually can be considered consisting of effectively non-interacting atoms. In dense matter, this is rephrased to clusters of atoms.

The first attempt to tackle the problem of phase transitions was based on the use of mean field theory. Loosely speaking, it presumes that correlations between mesoscopic stochastic variables can be treated perturbatively. Thus, in the sense of the central limit theorem, deviations from the ideal Gaussian behavior of fluctuations can be constructed (in quantum field theory, this procedure is also known as the Wick theorem). In most cases, some equation of state for macroscopic quantities can be directly obtained. In situations when the system is far from the critical region, the correlation length is very small and can be effectively neglected. There is no obvious discrepancy with the experimental data. However, it turns out that this approach exhibits large quantitative differences between the theory and the experiment in the vicinity of a critical point. The problem is that here the correlations are very large. In fact, directly at the transition point, the correlation length is divergent. Therefore, perturbation theory is no longer applicable (the coupling constant is very large and by no means can be treated as a small quantity). However, the divergent correlation length reveals another important effect in the physical picture: scale invariance. In contrast to the aforementioned symmetries, this is a case of a dynamically generated (emergent) symmetry. Mathematically, it is an invariance with respect to the special class of transformations that account for a change in the scale at which a physical system is studied. A divergent correlation length means that there is no special scale and the system looks and behaves at every scale in the same way.

A well-known observation in physics is that symmetry might have enormous influence on the properties of the physical system. In terms of scale invariance, the experimentally observed power laws for the functional dependence of various thermodynamic functions can be explained and also the various relations between critical exponents can be quantitatively estimated. Another important property of the second order phase transition is universality. In a simple formulation, it says that the behavior of the system near its critical point is fully determined by universal quantities—dimension of space, and number of components of order parameter, symmetry constraints—that are not characteristic of one system only, but a whole class of systems and also states that universal quantities do not depend on the model-dependent parameters—coupling constants, etc. Thus, in microscopic details very different physical systems such as strongly anisotropic magnetic material (the Ising model) and liquid might have the same critical properties. One only has to know a few very general pieces of information to classify a given physical system according to its critical behavior into some universality class, wherein all the systems behave in the same way.

It is interesting to realize that there is a great intrinsic similarity between quantum models in particle physics and statistical models. A profound relation between quantum field theory and statistical mechanics is revealed through the language of path integrals [5,6]. A classical random field in this framework is completely analogous to a fluctuating quantum field. Field models are often amenable to perturbation methods. Using Feynman diagrammatic technique, terms in perturbative expansion are expressed via Feynman diagrams that are a graphical representation of certain integrals. As a rule, they contain divergences in the range of large and small scales (wavevectors). In particle physics, there are no natural restrictions on scales, therefore it is necessary to find an effective procedure to eliminate these divergences step by step in each order of a concrete perturbation scheme. In perturbation expansions of classical field theories, natural scales usually exist: at small scales, the continuum description breaks down at atomic scales (nanophysics) and there is no reason to go below this lower limit. On the other hand, real quantities of matter are of finite size. In the theoretical description, however, it is convenient and customary to extrapolate results obtained for a finite quantity of homogeneous matter to the whole space when modeling real systems. Below, we demonstrate renormalization methods in the framework of the stochastic model of developed turbulence and related applications.

The method of renormalization group (RG) was proposed in the framework of the quantum field theory in the 1950s [1,7,8,9,10,11]. From the practical point of view, the RG method represents an effective way to determine non-trivial asymptotic behavior of Green functions (correlation functions) in the range of large (ultraviolet) or small (infrared) wavevectors (scales). The asymptotic behavior is non-trivial if, in a given order of a perturbative calculation, the divergences in a certain range of wavevectors appear (e.g., so-called large logarithms), which compensate for the smallness of the coupling constant g. In such case, summation of all terms of a perturbation series is needed. This summation can be carried out by means of the RG approach. Technically, one obtains linear partial differential RG equations for the Green functions. The coefficient functions (RG-functions) in the differential operator (see below) are calculated at a given order of the perturbation scheme. However, the solution of the RG equation represents the sum of an infinite series. For example, if the RG-functions are calculated at the lowest non-trivial order of the perturbation theory and the corresponding RG-equation is solved, the obtained result is a sum of leading logarithms of the whole perturbation series. Moreover, if the RG-functions are calculated with an improved precision, the solution of the RG equation includes corrections to the leading logarithms.

Notwithstanding the similarity of theoretical description of quantum field theory and classical statistical models, it has to be borne in mind that there is an essential difference in the interpretation and use of the RG in statistical physics on the one hand and in particle physics on the other hand. In particle physics, we are interested in an analysis of scaling in ultraviolet (UV) regime corresponding to large momenta. On the other hand, in statistical physics, asymptotic behavior in the opposite infrared (IR) limit of small momenta is usually studied. In both statistical mechanics and hydrodynamic transport problems, the interest in the IR behavior of statistical models is determined by the property of the basic field-theoretic tool—perturbation theory—to reproduce the observed singular behavior of certain physical quantities only in the limit of an infinite (flat) space. However, this infrared limit is usually rather sensitive to the large-scale structure of the model and care has to exercised when passing to the limit. In the case of equilibrium systems, the analysis is based on the Gibbs distribution, but, in the case of steady-state stochastic systems, there is no generic tool to this end which emphasizes the role of symmetry arguments.

The final aim of the theory (either in stochastic dynamics or developed turbulence) is to find the time-space dependence of statistical correlations—mainly those that can be experimentally measured. It turns out that use of quantum field theory methods (RG included) allows deriving a linear differential equation, which contains stable solutions in the asymptotic region of large macroscopic scales.

Solutions take a form of a product of a power-like term with a nontrivial exponent and scaling function of dimensionless variables (the scaling function is not determined by the RG method). To compute critical exponents in the form of asymptotic series, one has to resort to a certain scheme (we often employ variants of dimensional renormalization). Asymptotic properties of the scaling functions are analyzed by the operator product expansion, which is another theoretical tool developed mainly by Wilson, Wegner and Kadanoff [12,13,14]. In the stochastic theory of fully developed turbulence, scaling functions may be singular functions of dimensionless arguments and this can drastically change the critical exponents. The results demonstrate intermittent (multifractal) behavior of statistical correlations of the random fields of concentration of advected particles. Intermittency is a typical mesoscopic phenomenon, which is quantitatively revealed in singular behavior of correlation functions of velocity fluctuations with respect to an external turbulent spatial scale L.

The RG method not only leads to a quantitative description of the behavior near critical points, it also provides a new framework in which the aforementioned scale invariance and universality are naturally explained. At present, the use of the RG method in equilibrium statistical physics is well established and represents an important theoretical tool.

Contrary to the equilibrium physics, there are only few rigorous results in the case of non-equilibrium systems. Some of their properties are reminiscent of the equilibrium systems and thus it seems natural to apply the RG method to them. However, there exist also fundamental differences between them. Non-equilibrium systems might be divided [15] into two broad classes

• Systems with a Hermitian Hamiltonian, whose stationary states are described by Gibbs–Boltzmann distribution. Note that at the beginning they happen to be in a state far from the stationary (equilibrium) state. The dynamic description of such systems is obtained directly from static formulations. Examples include the Landau–Ginzburg equation for time evolution of local magnetization, kinetic Ising model, and models A-H for various models of critical dynamics [16]. All these equations are specific realizations of a rather general Langevin equation [5].
• Systems without Hermitian Hamiltonian or without Hamiltonian description at all, which in general do not need to have a stationary state. The detailed balance condition is not satisfied for them, which implies that Einstein relation between thermal fluctuations and friction forces cannot be stated. Typical examples of such systems cover: fluid in turbulent state, irreversible chemical processes, surface growing models, etc. Other approaches to such systems have to be used via quite general stochastic differential equation, which can be considered as an extension of a Langevin equation or using a master equation [17]. The former equation is suggested for some macroscopic quantity. Neglect of microscopic degrees of freedom is replaced by an introduction of random force. Then, according to underlying physical observations, properties of random force have to be specified. The latter approach is probably more fundamental, but also more difficult to handle.

In what follows, our main interest concerns specific problems of the second type related to hydrodynamics. In hydrodynamics, dissipation of mechanical energy to heat is an essential part of the physics. Since a Hamiltonian description is not possible in this case, we are dealing with steady states of the latter class.

As we know from everyday experience, fluids can exhibit very different behaviors from very simple, e.g., laminar flow, which is very predictable, to very chaotic, as is realized in turbulent motions. Turbulence is important in the analysis of phenomena in a wide range of scales from particle collisions in accelerators [18] and circulation of human blood [19] to the flow of air and water in the atmosphere and oceans, solar wind [20] and clusters of galaxies [21]. Let us stress that all studied systems belong to open systems that need a continuous input of energy in order to maintain steady state.

Theoretical analysis of turbulence is based on the statistical analysis of solutions of the Navier– Stokes problem. Symmetry and similarity arguments have allowed infering important conclusions about the scaling behavior of velocity correlation functions in the case of very large Reynolds numbers (the famous Kolmogorov theory in the first place) [22,23]. However, a more detailed statistical description of this fully developed turbulence as well as the onset of turbulence in a laminar flow are still lacking. Notwithstanding the rapid development of experimental methods [24], one of the major problems in the study of turbulence is the deficit of high-resolution experimental data. Therefore, numerical methods have become important tools in the investigation of turbulence [24] and provide solid benchmarks for testing of analytic results.

It is well-known that weather forecasting can be done for no more than a few days. This is caused by the intrinsic instability of Navier–Stokes (NS) equations, which are believed to describe motion of viscous (non-relativistic) fluids [22]. The formidable task of finding its solution remains one of the last unsolved classical problems [23]. For classification of various fluid states, the Reynolds number $\mathrm{Re}$ has been introduced. It is defined as $\mathrm{Re}=VL/\nu$, where V is typical average flow velocity, L is an external scale (e.g., a dimension of an obstacle, which causes perturbation to the regular flow) and $\nu$ is kinematic viscosity of the medium. It thus expresses the ratio between inertial and friction (dissipation) forces in a given fluid. In the case of low values, $\mathrm{Re}\ll 1$, regular (laminar) flow is observed. With an increasing value of $\mathrm{Re}$, very different phenomena occur ranging from the periodical ones as Kármán vortices to very chaotic irregular motion for the limit of very high values of $\mathrm{Re}\gg 1$ (in practice value $\mathrm{Re}\ge {10}^6$ is large enough) [23,25]. This state of fluid is known as fully developed turbulence.

At first sight, a very complicated problem turns out to be theoretically tractable because of appearance of new symmetries (again kind of emergent symmetry)—statistical symmetries. Kolmogorov postulated hypotheses [23,26] that could explain turbulence and also predict statistical and scaling properties of various correlation and structure functions. The Kolmogorov theory can be considered a kind of theory for “ideal” turbulence in the sense that it assumes the infinite value for the Reynolds number. These hypotheses are still not proved from the first principles—in this case from the Navier– Stokes equations. It should also be borne in mind that Kolmogorov’s hypotheses are stated for the case of homogeneous isotropic flat space but without any specific indication how this limit is approached when the Reynolds number grows without limit.

In contrast to mathematics, the physicist’s approach to turbulence follows a different path. Instead of considering a difficult mathematical problem related to boundary and initial conditions, their effect is replaced by properly chosen random force. The Navier–Stokes equation is amended by an additive random variable, which also mimics continuous input of mechanical energy into the system. The choice of the structure and statistics of the random force is the most essential point in modeling of the large-scale effects on the scale-invariant behavior predicted by the Kolmogorov theory. To this end, random forces concentrated at large spatial scales are used. In the basic setup of the stochastic problem, rotational symmetry of force correlations is assumed, but variation of the symmetry properties of force correlations (e.g., anisotropy and reflection asymmetry) may be used to probe the effect of large-scale properties on the scaling behavior of velocity correlations. The modeled large-scale induced effects include the appearance of a set of anomalous scaling dimensions in corrections to Kolmogorov scaling (multifractality) and magnetohydrodynamic dynamo. On the other hand, the Galilei invariance and small-scale anisotropy have been shown to be stable against large-scale perturbations.

In 1977, D. Forster, D. R. Nelson and M. J. Stephen applied the RG method to calculate the correlations of velocity field [27] governed by the stochastic Navier–Stokes equation with external random forcing. This work was motivated by Wilson’s momentum shell approach to RG, in which tracing out of fast degrees of freedom is supplemented with scale transformation [12,15]. Later, it was shown by C. De Dominicis and P. C. Martin [28] that in the range of small wave numbers the correlations of the velocity field manifest a scaling behavior with the celebrated Kolmogorov exponents. The stochastic NS equation was proposed to justify the Kolmogorov theory and has to be distinguished from the usual NS equations.

The essential idea of applying RG in the theory of developed turbulence consists in elimination of the direct influence of the modes with high wave numbers on observed quantities. Their influence is included in some effective variables, e.g., to the turbulent viscosity. Such an approach based on momentum shell approach was later developed further [29,30]. Let us note that in this paper we consider a different field-theoretic renormalization group technique [23,26,31,32], whose main advantage is more transparent and easier calculations.

Another interesting problem related to turbulence is the advection of some quantity [33,34] (temperature field, concentration field or tracer) by the turbulent field. In addition to the practical importance of such a problem, it is also very interesting from the theoretical point of view. It is still not clear to what extent turbulence is intermittent [23], i.e., what its fractal nature is. On the other hand, advection of a passive scalar quantity by simpler models (e.g., Kraichnan model that is described in detail below) than turbulence exhibits very strong intermittent behavior.

Naively, basic assumptions of Kraichnan-like models can be considered too crude. A typical approach to stochastic dynamics starts from an analysis of ideal systems—homogeneous in space and time, isotropic, incompressible (in the case of fluids), possessing mirror symmetry, etc. In the present review, the corresponding results for fully developed turbulence are summarized. However, real systems almost always exhibit some form of anisotropy, compressibility or violated mirror symmetry. The effect of such deviations from the ideal system on fluctuating random fields has been an object of intensive research activity, whose arguments and conclusions are described. The results have led to the general conclusion that such effects play a very important role. They can drastically change the large-scale behavior predicted by models of ideal systems.

In the original formulation of the Kraichnan model, the velocity field is assumed to be Gaussian, isotropic, incompressible and uncorrelated in time (white noise). More sophisticated models aim to incorporate effects of anisotropy, compressibility and finite correlation time [35]. Recent studies have pointed out some crucial differences between problems with vanishing and finite correlation time [33,36] and between the compressible and incompressible flows [37,38]. We employ the Kraichnan rapid-change ensemble to model the turbulent mixing [39]. Thus, we assume that velocity field is given by time-decorrelated Gaussian variable with the pair velocity function of the following form $\langle vv\rangle \propto \delta (t-t^{\prime })k^{-d-\epsilon }$, where $k=\left|k\right|$ is the wave number and $0<\epsilon <2$ is a free parameter of the theory. The physically most interesting value $\epsilon =\frac{4}{3}$ corresponds to the realistic (“Kolmogorov”) scaling behavior. This model gained popularity in the past mainly because of the insight it provides into the explanation of intermittency and anomalous scaling in turbulent flows. In the context of our study, it is worth mentioning that the Kraichnan ensemble allows a straightforward incorporation of compressibility, which appears complicated if the velocity is modeled by dynamical equations. The Kraichnan ensemble has been generalized further to the case of finite correlation time (see, e.g., [34,36,40] for the passive scalar and [41] for the passive vector fields). However, such synthetic models with non-vanishing correlation time are plagued by the lack of Galilean symmetry.

In Section 2, a short introduction to field-theoretic approach to stochastic dynamic is given and we briefly discuss the choice of the functional representation of the perturbation expansion for the solution of the Langevin equation. Section 3 is reserved for discussion of stochastic Navier–Stokes equation and basics of Kraichnan model. Section 4 is devoted to basic information of renormalization group approach and the related Section 5 to operator product expansion that is a specific method of the RG technique. In Section 6, Ward identities are used to obtain information about energy and momentum transfer in turbulent media. Section 7 is devoted to use of Ward identities in the RG analysis. Section 8 describes a dynamic restoration of initially broken Galilean symmetry. Section 9 is dedicated to a mechanism of spontaneous symmetry breaking in magnetohydrodynamic problem which is responsible for creation of magnetic dynamo. In Section 10, we discuss the effect of anisotropy and Section 11 is reserved for final remarks and comments.

2. Field-Theoretic Formulation

It is a well-known fact [5] that the failure of Landau theory of the second-order phase transition lies in the assumption of analyticity of the energy functional $F=F\left(\phi \right)$, where $\phi \left(x\right)$ is an order parameter configuration. The equilibrium physics of the phase transition is described by the order parameter at the minimum of the energy functional. The fluctuation theory of phase transitions takes as the fundamental quantity the random field $\phi \left(x\right)$, whose probability density function is defined by the energy functional as the effective Hamilton function of the Gibbs distribution. The difference of the fluctuation theory from the mean field theory of the microcanonical ensemble is to take the Landau functional as the fundamental Hamiltonian of the canonical ensemble instead that of an exact microscopic model. To calculate physical quantities, one has to average over all configurations of $\phi \left(x\right)$. Although a complete mathematical proof of the equivalence between microscopic and fluctuation model is missing, the latter approach has a very important and useful property. In contrast to the microscopic model, it is possible to use the RG method in order to analyze its behavior and to obtain quantitative predictions for critical exponents.

According to rather general arguments, many dynamic phenomena in nature exhibit a clear separation of time scales. For instance, typical time and space scales in (classical) critical systems diverge and this allows describing relevant physical quantities in terms of continuous fields. In fact, the latter corresponds simply to slow modes stemming from conservation laws or broken symmetries. On the other hand, fast degrees of freedom enter theoretical description through random noise fields. Similar reasoning applies to other systems that do not exhibit criticality. Famous examples encompass turbulence, reaction–diffusion problems, driven systems. A general class that covers such dynamical systems is known in the literature as stochastic dynamics. From the theoretical point of view, it is important that large scale properties can be properly taken into account through a formalism of Langevin-like equations.

The Langevin approach can be briefly summarized as follows. The aim is to study a slowly varying field or a set of fields $\phi$. Employing physical insight and symmetry reasoning, it is possible to postulate a stochastic differential equation of the form

$${\partial }_t\phi =V\left(\phi \right)+f,$$

(1)

where the functional $V=V\left(\phi \right)$ is local in the time variable, i.e., V depends only on the field $\phi$ and its spatial derivatives at a given time instant. In rare cases, V can be obtained from a microscopic model through a controlled coarse-grained procedure, but mostly its construction requires nontrivial knowledge about physical properties of a given stochastic system [42]. The random force f mimics the neglected rapid degrees of freedom and it is often modeled by means of Gaussian random variables. A nonzero mean value of f can be easily absorbed in the functional V. Thus, the complete statistical information about f is captured by specification of the first two moments

$$\langle f\left(x\right)\rangle =0,\langle f\left(x\right)f\left(x^{\prime }\right)\rangle =\delta (t-t^{\prime })D(\mathit{x},{\mathit{x}}^{\prime }),$$

(2)

where for brevity we have introduced the following notation $x=(t,\mathit{x})$, where $\mathit{x}$ is a d-dimensional vector. Let us note that, when necessary, we write the space dimension d explicitly. This permits a straightforward check of complicated expressions and also plays an important role in perturbative RG techniques such as dimensional regularization.

The crucial difference between critical and genuine non-equilibrium systems lies in the correlation function $D(\mathit{x},{\mathit{x}}^{\prime })$ in Equation (2). In critical dynamics, we know to what the system should relax. It should end up in the thermal equilibrium described by the Gibbs probability distribution ${\mathrm{e}}^{-\mathcal{H}/{\mathrm{k}}_{\mathrm{B}}T}$, which greatly restricts the form of D. On the other hand, in non-equilibrium systems, such a relation is broken and steady states (obtained in the limit of large time) are of much more complicated dynamical nature.

There exist many theoretical approaches, which can be undertaken for an investigation of the stochastic problem in Equations (1) and (2). Remarkable equivalence of stochastic problems with certain quantum field theory models offers a plethora of possibilities to use. Due to work of H.-K. Janssen [43] and C. De Dominicis [44], a given stochastic model can be cast into a path integral formulation, which is amenable to many theoretical methods such as Feynman diagrammatic technique, the RG method, and others. To provide background for later use of functional methods and sake of notation, we recall now the De Dominicis-0Janssen statement in a succinct manner [5]. First, let us rewrite potential term $V\left(\phi \right)$ as the sum $L\phi +n\left(\phi \right)$, where $L\phi$ represents the linear part in the field $\phi$ and $n\left(\phi \right)$ contains non-linearities. Then, stochastic problem in Equation (1) is tantamount to the quantum-field-theory model with double set of fields $\varphi \equiv \{\phi ,{\phi }^{\prime }\}$ and action functional of form [32]

$$\mathcal{S}\left[\varphi \right]=\frac{1}{2}{\phi }^{\prime }D{\phi }^{\prime }+{\phi }^{\prime }[-{\partial }_t\phi +L\phi +n\left(\phi \right)].$$

(3)

The auxiliary prime fields ${\phi }^{\prime }$ were put forward in [45] and are known as Martin–Siggia–Rose response fields. Hereinafter, we have employed a condensed notation, in which integrals over space-time and summations over repeated internal indices are implied. For instance, the second term in the action functional in Equation (3) is a shorthand for the expression

$${\phi }^{\prime }{\partial }_t\phi =\sum \int \mathrm{d}t\int {\mathrm{d}}^{d}x{\phi }^{\prime }(t,\mathit{x}){\partial }_t\phi (t,\mathit{x})\equiv \int \mathrm{d}x{\phi }_i^{\prime }\left(x\right){\partial }_t{\phi }_i\left(x\right),$$

(4)

where the index i numbers different field components. In this work, we are mainly interested in stochastic models concerning the velocity field $v$, which is a vector quantity and thus the appropriate summation over the vector (internal) index must be taken into account as well. In particular, an analogous expression to Equation (4) for the velocity field would be written as

$${\mathit{v}}^{\prime }{\partial }_t\mathit{v}=\int \mathrm{d}t\int {\mathrm{d}}^{d}x{v}_i^{\prime }(t,\mathit{x}){\partial }_t{v}_i(t,\mathit{x})\equiv \int \mathrm{d}x{v}_i^{\prime }\left(x\right){\partial }_t{v}_i\left(x\right).$$

(5)

The main goal of any statistical theory is to predict behavior of various correlation and response functions. Borrowing terminology from quantum field theory we refer to them as Green functions. These are defined as functional averages over the fields $\varphi$ with the weight $exp\mathcal{S}\left[\varphi \right]$, where $\mathcal{S}$ is action functional in Equation (3). Statistical averaging with the weight ${\mathrm{e}}^{\mathcal{S}}$ are denoted as follows

$$\langle \dots \rangle =\int \mathcal{D}\varphi \dots {\mathrm{e}}^{\mathcal{S}\left[\varphi \right]}.$$

(6)

All Green functions are effectively encoded into generating functional $\mathcal{G}$, which takes the form of the functional integral

$$\mathcal{G}\left[A\right]=\int \mathcal{D}\varphi exp\left[\mathcal{S}\left[\varphi \right]+A\varphi \right],$$

(7)

where $A\equiv (A_{\phi },A_{{\phi }^{\prime }})$ is the formal source and

$$A\varphi \equiv \int dx\left[A_{\phi }\left(x\right)\phi \left(x\right)+A_{{\phi }^{\prime }}\left(x\right){\phi }^{\prime }\left(x\right)\right].$$

(8)

Further, $\mathcal{D}\varphi$ in Equation (7) denotes the functional measure, i.e., $\mathcal{D}\varphi \equiv \mathcal{D}\phi \mathcal{D}{\phi }^{\prime }$, and the expression $\int \mathcal{D}\varphi \dots$ corresponds to a functional integral over the infinite dimensional space of all possible field configurations. Taking sufficiently many derivatives of $\mathcal{G}$ with respect to the formal sources at $A=0$ yields any permissible Green function of the theory. For example, the response function $\langle \phi {\phi }^{\prime }\rangle$ can be represented by the following functional integral

$$\langle \phi \left(x\right){\phi }^{\prime }\left(x^{\prime }\right)\rangle =\frac{{\delta }^2\mathcal{G}\left[A\right]}{\delta A_{\phi }\left(x\right)\delta A_{{\phi }^{\prime }}\left(x^{\prime }\right)}{|}_{A=0}=\int \mathcal{D}\varphi \phi \left(x\right){\phi }^{\prime }\left(x^{\prime }\right){\mathrm{e}}^{\mathcal{S}\left[\varphi \right]}.$$

(9)

The normalization factor, which ensures the equality $\mathcal{G}\left(0\right)=1$, has been included into the functional measure $\mathcal{D}\varphi$. The generating functional $\mathcal{G}$ for field-theoretic models might be interpreted as an analog of the partition function in equilibrium statistical physics [6]. The formal Taylor expansion of $\mathcal{G}$ reads

$$\mathcal{G}\left[A\right]=\sum _{n=0}^{\infty }\frac{1}{n!}\int \mathrm{d}{x}_1\cdots \int \mathrm{d}{x}_n{G}_n(x_1,\dots ,x_n)A\left(x_1\right)\cdots A\left(x_n\right),$$

(10)

where A on the left hand side stands for either $A_{\phi }$ or $A_{{\phi }^{\prime }}$ source field, and coefficient functions correspond to full Green functions

$$G_n(x_1,\dots ,x_n)=\langle \phi \left(x_1\right)\cdots \phi \left(x_n\right)\rangle =\frac{{\delta }^n\mathcal{G}\left[A\right]}{\delta A\left(x_1\right)\cdots \delta A\left(x_n\right)}{|}_{A=0}.$$

(11)

The formulation in Equation (3) of the stochastic problem is advantageous for the use of the powerful machinery of field-theoretic methods such as Feynman diagrammatic technique, RG method, and operator product expansion. The starting point of perturbative techniques is a separation of action $\mathcal{S}$ into a free part ${\mathcal{S}}_0$ and part containing nonlinearities ${\mathcal{S}}_{\mathrm{int}}={\phi }^{\prime }n\left(\phi \right).$ This division is not unique, but the necessary condition is the ability to solve the free part exactly. The free part ${\mathcal{S}}_0$ from Equation (3) can be symmetrized in the following way

$${\mathcal{S}}_0\left[\varphi \right]\equiv -\frac{1}{2}\varphi K\varphi \equiv -\frac{1}{2}\left(\begin{array}{c}\phi \\ {\phi }^{\prime }\end{array}\right)\left(\begin{array}{cc}0& {({\partial }_t-L)}^T\\ {\partial }_t-L& -D\end{array}\right)\left(\begin{array}{c}\phi \\ {\phi }^{\prime }\end{array}\right)$$

(12)

with the symmetric matrix K. Here, T denotes transposing, i.e., $K^T(x,x^{\prime })=K(x^{\prime },x).$ The inverse matrix $\Delta =K^{-1}$ defines the set of bare propagators ${\Delta }_{ik}(x,x^{\prime })={\langle {\varphi }_i\left(x\right){\varphi }_k\left(x^{\prime }\right)\rangle }_0$, which we number as follows

$${\Delta }_{12}={\Delta }_{21}^T={({\partial }_t-L)}^{-1},{\Delta }_{11}={\Delta }_{12}D{\Delta }_{21},{\Delta }_{22}=0,$$

(13)

where ${\varphi }_1\equiv \phi$, ${\varphi }_2\equiv {\phi }^{\prime }$. Generalization to a multicomponent field $\phi$ is obvious. The propagator ${\Delta }_{12}$ is retarded, therefore ${\Delta }_{21}={\Delta }_{12}^T$ is advanced. The symmetric propagator ${\Delta }_{11}={\Delta }_{11}^T$ contains both (retarded and advanced) contributions. The interaction part generates vertices with one field ${\phi }^{\prime }$ and two or more fields $\phi$, which are determined by the concrete form of the nonlinear terms in the action of the model. The aforementioned functional representation in Equation (3) permits construction of standard Feynman graphs for Green functions [5,6,46] by means of Wick’s theorem. The lines (propagators) are derived from the quadratic (free) part ${\mathcal{S}}_0$, whereas the interaction part $S_{\mathrm{int}}$ gives rise to vertices. Wick’s theorem (see, e.g., [5,47] for details) for the functional in Equation (7) may be compactly written in the exponential form

$$\begin{array}{c}\hfill \mathcal{G}\left[A\right]=exp(\frac{1}{2}\frac{\delta }{\delta \varphi }\Delta \frac{\delta }{\delta \varphi })exp({\mathcal{S}}_{\mathrm{int}}\left[\varphi \right]+A\varphi ){|}_{\varphi =0},\end{array}$$

(14)

where $\Delta$ is the matrix of propagators in Equation (13) and

$$\frac{\delta }{\delta \varphi }\Delta \frac{\delta }{\delta \varphi }\equiv \int \mathrm{d}x\int \mathrm{d}{x}^{\prime }\frac{\delta }{\delta {\varphi }_i\left(x\right)}{\Delta }_{ik}(x,x^{\prime })\frac{\delta }{\delta {\varphi }_k\left(x^{\prime }\right)}$$

(15)

is a shorthand notation for the universal differential operation and the indices $i,k$ enumerate all fields (response field included) in the model. Expansion of both exponents in Equation (14) leads to the celebrated Feynman diagrammatic technique, which allows perturbative calculation of all Green functions of the theory.

3. Stochastic Approach to Turbulence

The stochastic approach to the Navier–Stokes (NS) equation is analogous to fluctuation theory for critical phenomena mentioned in Section 1. It can be regarded as a microscopic approach to fully developed turbulence. The crucial difference from critical phenomena is that for turbulence there is no counterpart of Hamiltonian (free-energy) functional. The stochastic NS equation neglects such physical effects as influence of the boundaries or the precise form of system’s geometry (e.g., information about the way turbulence is produced), which are in the experiments responsible for creating turbulent instabilities. In a phenomenological sense, the input of energy is modeled by the proper choice of the stochastic force. The main goal of this theory is to justify Kolmogorov hypotheses [23,25]. A general proof of the equivalence between Kolmogorov hypotheses and the stochastic NS equations is still missing; nevertheless, as various studies show, it provides a nontrivial input to scaling behavior observed in turbulent flows [31,33].

The stochastic Navier–Stokes equation, which governs the dynamics of the velocity fluctuations $v_i=v_i\left(x\right),i=1,\dots ,d$, assumes the following form

$${\partial }_t\mathit{v}+(\mathit{v}·\nabla )\mathit{v}-{\nu }_0{\nabla }^2\mathit{v}+\nabla p=f,$$

(16)

where ${\nu }_0$ is the molecular kinematic viscosity, $p=p\left(x\right)$ stands for pressure fluctuations, ∇ is gradient, ${\nabla }^2=\Delta ={\partial }^2/\partial x_i\partial x_i={\partial }_i{\partial }_i$ is Laplace operator, and $\mathit{f}=\mathit{f}\left(x\right)$ represents an external random force per unit mass. For simplicity, we consider incompressible fluid with the solenoidal velocity $\nabla ·\mathit{v}=0$ and unit density of fluid ($\rho =1$). The incompressibility condition permits elimination of the pressure field from the stochastic Navier–Stokes equation (Equation (16)) and we can consider only its transverse components

$${\partial }_t\mathit{v}+P(\mathit{v}·\nabla )\mathit{v}-{\nu }_0{\nabla }^2\mathit{v}=\mathit{f},$$

(17)

where all fields are transverse, and P denotes the transverse projection operator, which in the momentum representation takes the form

$$P_{ij}\left(\mathit{k}\right)={\delta }_{ij}-k_i{k}_j/k^2$$

(18)

with $k=\left|\mathit{k}\right|$ being the magnitude of the wavevector $\mathit{k}$. In view of universality, it is assumed that the large-scale random force $\mathit{f}$ obeys the Gaussian distribution law. Hence, only the mean value and the second moment have to be postulated. The former takes zero value ($\langle f_i\rangle =0$) and pair correlation function is chosen in a general form

$$\langle f_i(t_1,{\mathit{x}}_1)f_j(t_2,{\mathit{x}}_2)\rangle \equiv D_{ij}(x_1,x_2).$$

(19)

It is convenient to specify the kernel function $D_{ij}$ in frequency–momentum representation

$$D_{ij}(x_1,x_2)\equiv \delta (t_1-t_2)d_{ij}({\mathit{x}}_1,{\mathit{x}}_2)=\delta (t_1-t_2)\int \frac{{\mathrm{d}}^{d}k}{{\left(2\pi \right)}^d}{P}_{ij}\left(\mathit{k}\right)d_f\left(k\right){\mathrm{e}}^{i\mathit{k}·(x_1-x_2)},$$

(20)

where d is a dimension of space. To employ the RG technique [28,31], the energy injection $d_f\left(k\right)$ is usually chosen in the power-law form

$$d_f\left(k\right)=D_0{k}^{4-d-2\epsilon }F\left(kL\right)$$

(21)

where L denotes outer integral scale, $D_0$ is the amplitude, and the scaling function $F\left(kL\right)$ possesses the unit asymptotic behavior in the range of large wave numbers $kL\gg 1$. For our purposes, it is sufficient to consider the “massless” theory for which Equation (21) becomes simply

$$d_f\left(k\right)=g_0{\nu }_0^3{k}^{4-d-2\epsilon },$$

(22)

with the additional feature that the corresponding integral in Equation (20) is IR regularized at $m\sim L^{-1}$. The parameter $D_0$ in Equation (22) is rewritten as $g_0{\nu }_0^3$ for dimensional and calculational reasons. The parameter $g_0$ plays the role of the coupling constant, $\epsilon \ge 0$ is a free parameter of the theory. For completeness, let us note that Equation (19) takes in frequency–momentum representation the following form

$$\langle f_i(\omega ,\mathit{k})f_j({\omega }^{\prime },{\mathit{k}}^{\prime })\rangle ={\left(2\pi \right)}^{d+1}\delta (\omega +{\omega }^{\prime })\delta (\mathit{k}+{\mathit{k}}^{\prime })P_{ij}\left(\mathit{k}\right)d_f\left(k\right).$$

(23)

From the mathematical point of view, Equation (16) represents a stochastic partial differential equation, first order in time variable. This allows us to employ the machinery of Section 2. Let us explain how these formal rules are applied to the theory of developed turbulence. According to the aforementioned De Dominicis–Janssen approach, the stochastic model described by Equation (17) is tantamount to the field-theoretic model with the action

$$\begin{array}{c}\hfill {\mathcal{S}}_{\mathrm{NS}}[\mathit{v},{\mathit{v}}^{\prime }]=\frac{1}{2}{v^{\prime }}_i{D}_{ij}{v^{\prime }}_j+{\mathit{v}}^{\prime }·\left[-{\partial }_t\mathit{v}+{\nu }_0{\nabla }^2\mathit{v}-(\mathit{v}·\nabla )\mathit{v}\right],\end{array}$$

(24)

where $D_{ij}$ is introduced into Equation (19), the auxiliary response vector field ${\mathit{v}}^{\prime }$ is solenoidal $(\nabla ·{\mathit{v}}^{\prime }=0)$ as well as the velocity field $\mathit{v}$, and ${\nu }_0$ is the bare (molecular) viscosity coefficient. To distinguish it from the renormalized (turbulent) viscosity $\nu$, which is generated in the process of the renormalization, we denote it and other similar (bare) parameters by the subscript “0”. We stress that this notation is used in the whole work.

Feynman rules for the perturbation theory are constructed by means of the general operation in Equation (14)). The explicit form of the propagators is determined by the quadratic part of the action in Equation (24) and in the frequency–momentum representation they are

$$\begin{array}{c}\hfill {\Delta }_{ij}^{vv}({\omega }_k,\mathit{k})=\frac{P_{ij}\left(\mathit{k}\right)d_f\left(k\right)}{{\omega }_k^2+{\nu }_0^2{k}^4},{\Delta }_{ij}^{v{v}^{\prime }}({\omega }_k,\mathit{k})={\left({\Delta }_{ij}^{v^{\prime }v}({\omega }_k,\mathit{k})\right)}^{*}=\frac{P_{ij}\left(\mathit{k}\right)}{-i{\omega }_k+{\nu }_0{k}^2},{\Delta }_{ij}^{v^{\prime }{v}^{\prime }}({\omega }_k,\mathit{k})=0,\end{array}$$

(25)

where * denotes complex conjugation and the transverse projector appears due to incompressibility condition. In the time–momentum representation, the corresponding expressions are

$$\begin{array}{cc}\hfill {\Delta }_{ij}^{vv}(\mathit{k},t^{\prime }-t)& =\frac{P_{ij}\left(\mathit{k}\right)d_f\left(k\right)}{2{\nu }_0{k}^2}{\mathrm{e}}^{-{\nu }_0{k}^2\left|t^{\prime }-t\right|},\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill {\Delta }_{ij}^{v{v}^{\prime }}(\mathit{k},t^{\prime }-t)& =\theta \left(t^{\prime }-t\right)P_{ij}\left(\mathit{k}\right){\mathrm{e}}^{-{\nu }_0{k}^2\left(t^{\prime }-t\right)},\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill {\Delta }_{ij}^{v^{\prime }v}(\mathit{k},t^{\prime }-t)& =\theta \left(t-t^{\prime }\right)P_{ij}\left(\mathit{k}\right){\mathrm{e}}^{-{\nu }_0{k}^2\left(t^{\prime }-t\right)},\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill {\Delta }_{ij}^{v^{\prime }v}(\mathit{k},t^{\prime }-t)& =0.\hfill \end{array}$$

(29)

Here, the step function $\theta \left(t\right)$ displays an important physical feature of the propagator ${\Delta }^{v{v}^{\prime }}$—its retardation. In fact, ${\Delta }^{v{v}^{\prime }}$ is the leading order contribution to the response function $\langle v{v}^{\prime }\rangle$ of the original model in Equations (16)–(21). The propagator ${\Delta }^{vv}$ represents the leading contribution to the pair correlation function of the velocity field $W_{ij}=\langle v_i{v}_j\rangle$. With coinciding time arguments, the latter is proportional to the kinetic energy spectrum $E\left(k\right)$ in the wavevector representation. This function enters the equation of energy balance describing the transfer of the kinetic energy from the largest spatial scales to the smallest ones, where it dissipates to heat [23]. The vertex factor

$$V_m(x_1,x_2,\dots ,x_m;\varphi )=\frac{{\delta }^{m}V\left[\varphi \right]}{\delta \varphi \left(x_1\right)\delta \varphi \left(x_2\right)\dots \delta \varphi \left(x_m\right)}$$

(30)

is associated to each interaction vertex of a Feynman graph. In Equation (30), the dummy field $\varphi$ is one from the set of all fields $\{v^{\prime },v\}$. The interaction vertex in Equation (24) is cast in a more convenient form

$$-\int \mathrm{d}t\int {\mathrm{d}}^{d}x{\mathit{v}}^{\prime }(\mathit{v}·\nabla )\mathit{v}=-\int \mathrm{d}t\int {\mathrm{d}}^{d}x{v^{\prime }}_i{v}_k{\partial }_k{v}_i=\left({\partial }_k{v^{\prime }}_i\right)v_k{v}_i,$$

(31)

where the incompressibility condition ${\partial }_i{v}_i=0$ and integration by parts have been used. The latter step requires the standard assumption of rapid enough vanishing of velocity in the limit $\left|\mathit{x}\right|\to \infty$. Furthermore, the last expression in Equation (31) corresponds to the shorthand of Equation (4). Rewriting functional in Equation (31) in the symmetric form $v_i{V}_{ijl}{v}_j{v}_l/2$, we derive the explicit form for the corresponding vertex factor in the Fourier representation

$$V_{ijl}=i(k_j{\delta }_{il}+k_l{\delta }_{ij}).$$

(32)

Here, the wavevector $k$ is flowing in the vertex through the field $v^{\prime }$ and is denoted by slash in Figure 2. The propagators (lines) $\Delta$ and vertices V are graphically depicted in Figure 1 and Figure 2.

The theoretical description of the fluid turbulence on the basis of “first principles”, i.e., starting from the stochastic Navier–Stokes (NS) equation [25] remains an open problem. However, considerable progress has been made in understanding simplified model systems sharing certain essential properties with the real problem: stochastic Burgers equation [48], shell models [49] and advection by random “synthetic” velocity fields [33].

A paradigmatic model of a scalar quantity advected passively by a Gaussian random velocity field, uncorrelated in time and self-similar in space, the so-called Kraichnan’s rapid-change model [39], is a famous example. The standard notation for advection problem using the Kraichnan model slightly differs from the one using stochastic Navier–Stokes ensemble. Therefore, in what follows, we give a brief overview of basic physical ideas behind the Kraichnan model and introduce the corresponding notation.

The governing equation for diffusion–advection for field $\theta$ is

$${\partial }_t\theta +(\mathit{v}·\nabla )\theta -D_0{\nabla }^2\theta =f^{\theta },$$

(33)

where $D_0$ is the coefficient of molecular diffusivity and $f^{\theta }\equiv f^{\theta }\left(x\right)$ is a zero-mean Gaussian random noise with the correlation function

$$\langle f^{\theta }\left(x\right)f^{\theta }\left(x^{\prime }\right)\rangle =\delta (t-t^{\prime })C(\mathit{r}/L),\mathit{r}=\mathit{x}-{\mathit{x}}^{\prime }.$$

(34)

The noise $f^{\theta }$ in Equation (33) maintains the steady-state of the system. The particular form of the correlator is not relevant, however. The sole condition which must be satisfied by the function $C(\mathit{r}/L)$ is that it must fall off rapidly at $r\equiv \left|\mathit{r}\right|\gg L$. Here, L is an integral scale related to the stirring.

In accordance with the generalized Kraichnan model [34,50] with finite correlation time taken into account, we assume velocity field generated by a simple linear stochastic equation

$${\partial }_t{v}_i+R{v}_i=f_i^v,$$

(35)

where $R\equiv R\left(x\right)$ is a linear operation to be specified below and $f_i^v\equiv f_i^v\left(x\right)$ is a zero-mean random stirring force with the correlator

$$\langle f_i^v\left(x\right)f_j^v\left(x^{\prime }\right)\rangle \equiv D_{ij}^f(x;x^{\prime })=\frac{1}{{\left(2\pi \right)}^{d+1}}\int \mathrm{d}\omega \int {\mathrm{d}}^{d}k{P}_{ij}\left(\mathit{k}\right)D^f({\omega }_k,\mathit{k}){\mathrm{e}}^{-i(t-t^{\prime })+i\mathit{k}·(x-x^{\prime })}.$$

(36)

It should be noted that, in the SDE in Equation (33), the multiplicative noise due to random velocity is not a white noise in time as in the original Kraichnan model. Therefore, there is no need to specify the interpretation of the SDE. However, in the analysis, the white-noise limit is considered and it should recalled that in this limit the results correspond to the Stratonovich interpretation of the SDE in Equation (33).

The correlator $D^f$ is chosen [34,35,36] in the following form

$$D^f({\omega }_k,\mathit{k})=g_0{\nu }_0^3{(k^2+m^2)}^{2-d/2-\epsilon /2-\eta /2},$$

(37)

with the wavenumber representation of the function $R\left(x\right)$:

$$R\left(\mathit{k}\right)=u_0{\nu }_0{(k^2+m^2)}^{1-\eta /2}.$$

(38)

The positive amplitude factors $g_0$ and $u_0$ are the coupling constants of the model. Furthermore, $g_0$ can be regarded as a formally small parameter of the perturbation theory. The positive exponents $\epsilon$ and $\eta$ ($\epsilon =\mathcal{O}\left(\eta \right)$) are RG expansion parameters. They are analogous to expansion parameter $\epsilon =4-d$ in the ${\phi }^4-$ theory. Now, the expansion is carried out in the $(\epsilon ,\eta )$-plane around the origin $\epsilon =\eta =0$.

Note the presence of two scales in the problem—integral scale L introduced in Equation (34) and momentum scale m, which has appeared in Equation (38). Clearly, they have different physical origins. However, these two scales can be related to each other and for technical purposes [35] it is reasonable to choose $L=1/m$. When not explicitly stated, this relation is always assumed.

In the limit $k\gg m$ the functions in Equations (37) and (38) take on a simple powerlike form

$$D^f({\omega }_k,\mathit{k})=g_0{\nu }_0^3{k}^{4-d-\epsilon -\eta },R\left(k\right)=u_0{\nu }_0{k}^{2-\eta },$$

(39)

which is convenient for actual calculations. The needed IR regularization will be given by restrictions on the region of integrations.

From Equations (35), (36) and (39), the statistics of the velocity field $v$ can be determined. It obeys Gaussian distribution with zero mean and correlator

$$\langle v_i(t,\mathit{x})v_j(0,\mathbf{0})\rangle =\int \frac{\mathrm{d}{\omega }_k}{2\pi }\int \frac{{\mathrm{d}}^{d}k}{{\left(2\pi \right)}^d}{D}_v({\omega }_k,\mathit{k}){\mathrm{e}}^{-i{\omega }_{k}t+i\mathit{k}·x},$$

(40)

where the kernel function $D_v({\omega }_k,\mathit{k})$ is assumed in the form

$$D_v({\omega }_k,\mathit{k})=P_{ij}\left(\mathit{k}\right)\frac{g_{10}{u}_{10}{D}_0^3{k}^{4-d-\epsilon -\eta }}{{\omega }_k^2+u_{10}^2{D}_0^2{\left(k^{2-\eta }\right)}^2}.$$

(41)

The correlator in Equation (41) is directly connected to the energy spectrum via the frequency integral [34,51,52,53,54,55]

$$E\left(k\right)\simeq k^{d-1}\int \mathrm{d}\omega D^v({\omega }_k,k)\simeq \frac{g_0{\nu }_0^2}{u_0}{k}^{1-\epsilon }.$$

(42)

Hence, the coupling constant $g_0$ and the exponent $\epsilon$ characterize the equal-time velocity correlator or, similarly, energy spectrum. Further, the parameter $u_0$ and the exponent $\eta$ are related to the frequency ${\omega }_k\simeq u_0{\nu }_0{k}^{2-\eta }$ (or to the function $R\left(k\right)$, the reciprocal of the correlation time at the wave number k) which describes the mode with wave number k [34,51,52,53,54,55,56,57]. Let us note that in the chosen notation the value $\epsilon =8/3$ corresponds to the well-known Kolmogorov “five-thirds law” for the spatial scaling behavior of the velocity field, and the value $\eta =4/3$ corresponds to the Kolmogorov frequency. A straightforward dimensional analysis reveals that the parameters (charges) $g_0$ and $u_0$ are connected to the ultraviolet (UV) momentum scale $\Lambda$ (of the same order of magnitude as the inverse Kolmogorov length) by the relations

$$g_0\simeq {\Lambda }^{\epsilon +\eta },u_0\simeq {\Lambda }^{\eta }.$$

(43)

In Ref. [50], it was demonstrated that the linear model in Equation (35) (and consequently the Gaussian model in Equation (40) as well) is not invariant under Galilean transformation and, therefore, it effectively neglects important effect of the self-advection of turbulent eddies. As a result of these so-called “sweeping effects” the different time correlations of the velocity are not self-similar and exhibit strong dependence on the integral scale [58,59,60,60]. However, the results presented in Ref. [50] lead to the conclusion that the Gaussian model describe the passive advection reasonably well in the appropriate frame of reference, in which the mean velocity field vanishes. An additional argument to support the model in Equation (40) is that we are mainly interested in the equal-time, Galilean invariant quantities (e.g., structure functions), which are not affected by the sweeping. Therefore, their absence in the Gaussian model in Equation (40) is not relevant [34,36,40].

The kernel function in Equation (41) is written in a very general form and allows studying various special limits, in which the numerical analysis of the resulting equations is simplified and which provide a deeper physical insight. Possible limiting cases are

• The rapid-change model corresponding to the limit $u_{10}\to \infty ,g_{10}^{\prime }\equiv g_{10}/u_{10}=const$. Then, the kernel function becomes

  $$D_v({\omega }_k,\mathit{k})\propto g_{10}^{\prime }{D}_0{k}^{-d-\epsilon +\eta }.$$

  (44)
  The velocity correlator is obviously $\delta$—correlated in the time variable.
• The frozen velocity field arising in the limit $u_{10}\to 0$, in which the kernel function corresponds to

  $$D_v({\omega }_k,\mathit{k})\propto g_0{D}_0^2\pi \delta \left({\omega }_k\right)k^{2-d-\epsilon }.$$

  (45)
• The purely potential velocity field obtained in the limit $\alpha \to \infty$ with $\alpha g_{10}=$ constant. This case is similar to the model of random walks in a random environment with long-range correlations [61,62].
• The turbulent advection, for which $\epsilon =2\eta =8/3$. This choice mimics properties of the fully developed turbulence and yields well-known Kolmogorov scaling [23].

Using Equation (3), the stochastic problem in Equations (33)–(36) can be recast into the equivalent field theoretic model of the doubled set of fields $\varphi \equiv \{\theta ,{\theta }^{\prime },\mathit{v},{\mathit{v}}^{\prime }\}$ with the action functional

$$\mathcal{S}\left[\varphi \right]=\frac{1}{2}{\mathit{v}}^{\prime }{D}^f{\mathit{v}}^{\prime }+{\theta }^{\prime }\left[-{\partial }_t\theta -(\mathit{v}·\nabla )\theta +{\nu }_0{\nabla }^2\theta \right]+{\mathit{v}}^{\prime }·\left[-{\partial }_t\mathit{v}-R\mathit{v}\right],$$

(46)

where $D_{ij}^f$ is defined in Equation (36), and as usual ${\theta }^{\prime }$ and ${\mathit{v}}^{\prime }$ are auxiliary response fields.

Generating functional of full Green functions $\mathcal{G}\left[A\right]$ is defined by Equation (7), where now a linear form $A\phi$ is defined as

$$A\varphi =A_{\theta }\theta +A_{{\theta }^{\prime }}{\theta }^{\prime }+A{v}_i{v}_i+A{v}_i^{\prime }{v}_i^{\prime }.$$

(47)

Following the argument in [34], we set $A_i^{v^{\prime }}=0$ in Equation (47) and carry out the explicit Gaussian integration over the auxiliary vector field, ${\mathit{v}}^{\prime }$ because we are not interested in the Green functions containing the auxiliary field ${\mathit{v}}^{\prime }$. After the integration, we are left with the field-theoretic model with the action functional

$$\mathcal{S}\left[\varphi \right]=-\frac{1}{2}\mathit{v}{\left(D^{\mathit{v}}\right)}^{-1}\mathit{v}+{\theta }^{\prime }[-{\partial }_t\theta -(v·\nabla )\theta +{\nu }_0{\nabla }^2\theta ],$$

(48)

where the second term represents the De Dominicis–Janssen action for the stochastic problem in Equation (33) at fixed velocity field $\mathit{v}$. The first term describes the Gaussian averaging over $\mathit{v}$ specified by the correlator $D^v$. The latter explicitly reads

$${\mathcal{S}}_{\mathrm{vel}}\left[\mathit{v}\right]=\frac{1}{2}\int \mathrm{d}{x}_1\int \mathrm{d}{x}_2{\mathit{v}}_i\left(x_1\right)D_{ij}^{-1}(x_1-x_2){\mathit{v}}_j\left(x_2\right).$$

(49)

The action in Equation (48) is written in a form that is suitable for a straightforward application of the field-theoretic perturbative analysis with the use of the standard Feynman diagrammatic technique. From the quadratic part of the action, we derive the matrix of bare propagators. The wavenumber frequency representations of relevant propagators are: (a) the bare propagator ${\langle \theta {\theta }^{\prime }\rangle }_0$ defined as follows

$${\langle \theta {\theta }^{\prime }\rangle }_0={\langle {\theta }^{\prime }\theta \rangle }_0^{*}=\frac{1}{-i\omega +{\nu }_0{k}^2};$$

(50)

and (b) the bare propagator for the velocity field ${\langle vv\rangle }_0$ that reads

$${\langle v_i{v}_j\rangle }_0=P_{ij}\left(\mathit{k}\right)D_v(\omega ,\mathit{k}).$$

(51)

The triple (interaction) vertex $-{\theta }^{\prime }{v}_j{\partial }_j$ can be rewritten in $\theta ={\theta }^{\prime }{v}_j{V}_j\theta$, where momentum $\mathit{k}$ is flowing into the vertex via the response field ${\theta }^{\prime }$. A graphical representation of the perturbation elements for a Kraichnan-like model is schematically depicted in Figure 3.

Taking as a example the Kraichnan model, let us briefly describe what kind of quantities might be studied by functional techniques. From experimental and theoretical point of view, the main focus is in the behavior of the equal-time structure functions

$$S_N\left(r\right)\equiv \langle {[\theta (t,\mathit{x}+\mathit{r})-\theta (t,\mathit{x})]}^N\rangle$$

(52)

in the inertial range, specified by the inequalities $l\sim 1/\Lambda \ll r\equiv \left|\mathit{r}\right|\ll L=1/m$ (l is an internal length). Brackets $\langle \cdots \rangle$ denote the functional average over fields $\varphi =\{\theta ,{\theta }^{\prime },v\}$ with the weight functional $exp\mathcal{S}\left[\mathsf{\Phi }\right]$ from Equation (48). In the isotropic case, the odd functions $S_{2N+1}$ vanish identically, while for even functions $S_{2N}$ a simple dimensional argument dictates the following form

$$S_{2N}\left(r\right)={\nu }_0^{-N}{r}^{2N}{R}_{2N}(r/l,r/L,g_0,u_0),$$

(53)

where $R_{2N}$ are scaling functions of purely dimensionless variables. In principle, functions $R_{2N}$ can be calculated by means of the usual perturbation theory (i.e., as series in $g_0$). However, this is not a reasonable way to study the inertial-range behavior: the reason is that the coefficients are singular in the limits $r/l\to \infty$ and/or $r/L\to 0$ and compensate for the smallness of $g_0$. To obtain correct IR behavior the entire series have to be summed. Such a summation procedure can be effectively done by the use of the field theoretic RG and OPE (see Section 4 and Section 5).

The RG analysis can be divided into two stages. During the first stage, the multiplicative renormalizability of the model is proved and the differential RG equations for its correlation (structure) functions are derived. The asymptotic behavior of functions similar to the one in Equation (52) for $r/l\gg 1$ and any fixed $r/L$ is governed by IR stable fixed points (see next section) of the RG equations and assumes the form

$$S_{2N}\left(r\right)={\nu }_0^{-N}{r}^{2n}{(r/l)}^{-{\gamma }_N}{R}_{2N}(r/L),r/l\gg 1$$

(54)

with so far unknown scaling functions $R_{2N}(r/L)$. Whenever ${\gamma }_N$ is a nonlinear function of N, we refer to such case as anomalous scaling or multiscaling.

Let us remind that the quantity $\Delta \left[S_{2N}\right]\equiv -2N+{\gamma }_N$ is called the critical dimension. The exponent ${\gamma }_N$, the difference between the critical dimension $\Delta \left[S_{2N}\right]$ and the canonical dimension $-2N$, is known as the anomalous dimension. In the present case, the latter takes a simple form: ${\gamma }_N=nϵ$. For any function $R_N(r/L)$, the representation in Equation (54) implies scaling behavior in the IR region ($r/l\gg 1$, $r/L$ fixed) with definite critical dimensions of all IR relevant parameters, $\Delta \left[S_{2N}\right]=-2N+Nϵ$, ${\Delta }_r=-1$, ${\Delta }_{L^{-1}}=1$ and fixed irrelevant parameters ${\nu }_0$ and l.

In the second stage, the small $r/L$ behavior of the functions $R_{2N}(r/L)$ is analyzed in the general representation in Equation (54) employing the OPE technique (Section 5). It predicts that, in the limit $r/L\to 0$, the functions $R_{2N}(r/L)$ have the asymptotic forms

$$R_{2N}(r/L)=\sum _F{C}_F(r/L){(r/L)}^{{\Delta }_N},$$

(55)

where $C_F$ are coefficients regular in the variable $r/L$. In general, the summation is performed over specific renormalized composite operators F with critical dimensions ${\Delta }_n$. Kraichnan model exhibits nontrivial scaling behavior as some of anomalous exponents ${\Delta }_N$ are negative and singular behavior on L is present. Such situation never occurs in critical phenomena [5,6] where corresponding exponents are positive and lead only to subleading corrections.

More elaborated discussion on anomalous scaling can be found in Section 10, which is devoted to generalization of Kraichnan model. Namely, assumption of isotropy is abandoned and effect of anisotropy is taken into account.

4. Renormalization Group Analysis

Let us briefly summarize main ideas of the quantum-field theory of renormalization and RG technique; a detailed account can be found in monographs [5,6,15,46].

Feynman graphs of Green functions are a convenient graphical representation of perturbation theory. Quantum field theory models are well-known for appearance of UV divergences in loop diagrams. This results from an integration at large momenta. Therefore, it is necessary to find an effective procedure to eliminate these divergences step by step in a controlled manner. Finite diagrams (free of UV divergences) are brought by an iterative renormalization procedure. The inherent ambiguity of this removal of divergences may be used to establish connection between values of Green functions at different scales without having explicit solutions for them. This is the idea of the method of renormalization group (RG) proposed in the framework of high energy physics long time ago [1,7,8,9,10,11].

In statistical physics, RG allows one to extract relevant information about large-scale behavior from the mutual correlation between IR behavior of Green functions and UV divergences at critical dimension [5,15]. Thus, in statistical theories, the RG method can be understood as an effective way to determine non-trivial asymptotic behavior of Green functions in the range of small (infrared) wavevectors (scales).

There is simple criterion how to determine the true asymptotic range in the framework of the RG. One of the RG-functions is the $\beta$-function, which is the coefficient of the operation ${\partial }_g$ in the RG equation. The $\beta$-function is calculated perturbatively as infinite series of powers of the coupling constant g. Non-trivial asymptotic behavior is governed by RG fixed points $g^{*}$, which are roots of the $\beta$ function (solutions of equation $\beta \left(g\right)=0$). A fixed point can be IR or UV stable depending on the behavior of the $\beta$-function in the vicinity of $g^{*}$. Of course, physical theories contain usually many charges and these considerations have to be properly generalized [5,46].

The field theoretic RG is based on non-trivial techniques of UV renormalization. The basic procedure lies in a perturbative calculation of the RG-functions in the framework of a prescribed scheme of regularization [5,6]. To find and analyze UV divergences in a specific field-theoretic model counting of canonical scaling dimensions of fields and parameters of the model is used. The essence of such a power counting is closely connected with the existence of a scale invariance in the model.

For models considered in this paper, it is advantageous to calculate Feynman diagrams in a formal scheme [5] without UV-cut-off $\Lambda$. Then, UV-divergences manifest themselves as poles in a dimensionless parameter $\epsilon$ that measures deviation from a logarithmic theory, i.e., a theory in which all coupling constants become dimensionless. The procedure of multiplicative renormalization removing UV-divergences (in the present case, poles in a parameter $\epsilon$) is the following: the original action $\mathcal{S}[\varphi ,e_0]$ is declared to be unrenormalized; its parameters $e_0$ (the letter $e_0$ stands for the whole set of parameters; for instance, coupling constants, deviation from criticality, viscosity, etc.) are the bare parameters, and they are assumed to be functions of the new renormalized parameters e. The new renormalized action is the functional ${\mathcal{S}}_R\left[\varphi \right]=\mathcal{S}[\varphi Z_{\varphi },e{Z}_e]$ with certain (to be determined perturbatively such that the Green functions generated by the renormalized action are UV finite, i.e., regular in $\epsilon$) renormalization constants of fields $Z_{\varphi }$ (one per each independent component of the field) and parameters $Z_e$. In unrenormalized full Green functions $G_N=\langle \varphi \cdots \varphi \rangle$, the functional averaging $\langle \cdots \rangle$ is performed with the weight functional $exp\mathcal{S}\left[\varphi \right]$, whereas, in renormalized functions $G_N^R$, with the renormalized weight functional $exp{\mathcal{S}}_R\left[\varphi \right]$. The relation between the functionals $\mathcal{S}\left[\varphi \right]$ and ${\mathcal{S}}_R\left[\varphi \right]$ leads to the relation between the corresponding Green functions $G_N^R=Z_{\varphi }^{-N}{G}_N$, where by definition $G_N=G_N(e_0,\epsilon \cdots )$ (ellipsis denotes other arguments such as coordinates or wavenumbers), and, by convention, the quantities $G_N^R$ and $Z_{\varphi }$ are expressed in terms of the parameters e. The correspondence $e_0\leftrightarrow e$ within perturbation theory is assumed to be one-to-one, therefore either of the sets $e_0,e$ can be taken as the independent variables.

For translationally invariant theories, it is much more convenient to deal not with the full Green functions $G_N$, but with their connected parts $W_N$. Their generating functional being through the relation

$$\mathcal{W}\left[A\right]\equiv ln\mathcal{G}\left[A\right].$$

(56)

A further simplification is possible through 1-irreducible functions ${\mathsf{\Gamma }}_N$ (also called one particle irreducible functions or vertex functions). The generating functional for the latter is defined by the functional Legendre transform [47]

$$\mathsf{\Gamma }\left[\alpha \right]\equiv \mathcal{W}\left[A\right]-A\alpha ,$$

(57)

where

$$A\alpha =\int \mathrm{d}x(A_{\phi }\left(x\right){\alpha }_{\phi }\left(x\right)+A_{{\phi }^{\prime }}\left(x\right){\alpha }_{{\phi }^{\prime }}\left(x\right)),{\alpha }_{\phi }=\frac{\delta \mathcal{W}\left[A\right]}{\delta A_{\phi }\left(x\right)},{\alpha }_{{\phi }^{\prime }}=\frac{\delta \mathcal{W}\left[A\right]}{\delta A_{{\phi }^{\prime }}\left(x\right)}.$$

(58)

To simplify notation in practical calculations, it is convenient to relabel $\alpha$-variables back to the original fields $\varphi$. This allows us to rewrite the first relation in Equation (57) compactly as

$$\mathsf{\Gamma }\left[\varphi \right]=\mathcal{S}\left[\varphi \right]+\tilde{\mathsf{\Gamma }}\left[\varphi \right],$$

(59)

where $\tilde{\mathsf{\Gamma }}\left[\varphi \right]$ is the sum of all one particle irreducible (1PI) loop diagrams [5].

Statements of RG theory are readily summarized at the level of corresponding Green functions. For connected and 1PI Green functions, they read

$$\begin{array}{ccc}\hfill W_N^R(e,\epsilon ,\cdots )& =& Z_{\varphi }^{-N}(e,\epsilon )W_N(e_0(e,\epsilon ),\epsilon ,\cdots ),\hfill \end{array}$$

(60)

$$\begin{array}{ccc}\hfill {\mathsf{\Gamma }}_N^R(e,\epsilon ,\cdots )& =& Z_{\varphi }^N(e,\epsilon ){\mathsf{\Gamma }}_N(e_0(e,\epsilon ),\epsilon ,\cdots ),\hfill \end{array}$$

(61)

where the functions $e_0(e,\epsilon ),$$Z_{\varphi }^N(e,\epsilon )$ can be chosen arbitrarily, which implies an arbitrary choice of normalization of the fields and parameters e at given $e_0$. In the present text, we also interchangeably use the following notation for the connected Green functions

$$W_{{\varphi }_1\dots {\varphi }_N}\equiv {\langle {\varphi }_1\dots {\varphi }_N\rangle }_{\mathrm{conn}}.,$$

(62)

and for the 1PI Green functions according to the aforementioned relabeling $\alpha \to \varphi$

$${\mathsf{\Gamma }}_{{\varphi }_1\dots {\varphi }_N}\equiv {\langle {\varphi }_1\dots {\varphi }_N\rangle }_{1-\mathrm{ir}}.$$

(63)

The crucial statement of the theory of renormalization is that for the multiplicatively renormalizable models these functions can be chosen to provide UV finiteness of Green functions as $\epsilon \to 0.$ With this choice, all UV divergences (poles in $\epsilon$) contained in the functions $e_0(e,\epsilon ),$$Z_{\varphi }^N(e,\epsilon )$ are absent in renormalized Green functions $W_N^R(e,\epsilon )$. We note that the UV finiteness in this sense of any one set of Green functions (full, connected, and 1-irreducible) automatically leads to the UV-finiteness of any other. The RG equations are written for the renormalized functions $W_N^R$ which differ from the original unrenormalized functions $W_N$ only by normalization, and, therefore, can be used equally well to analyze the critical scaling. Let us recall an elementary derivation of the RG equations [5,46]. The requirement of elimination of divergences does not uniquely determine the functions $e_0(e,\epsilon )$ and $Z_{\varphi }(e,\epsilon )$. An arbitrariness remains which allows introducing an additional dimensional parameter scale setting parameter (renormalization mass) $\mu$ in these functions (and via them also into $W_N^R$)

$$W_N^R(e,\mu ,\epsilon ,\cdots )=Z_{\varphi }^{-N}(e,\mu ,\epsilon )W_N(e_0(e,\mu ,\epsilon ),\epsilon ,\cdots ).$$

(64)

A change of $\mu$ at fixed $e_0$ leads to a change of $e,Z_{\varphi }$ and $W_R$ for unchanged $W_N(e_0,\epsilon ,\cdots )$. We denote by ${\tilde{\mathcal{D}}}_{\mu }$ the differential operator $\mu {\partial }_{\mu }$ for fixed $e_0$ and apply it to both sides of the equation $Z_{\varphi }{W}_N^R=W_N$ with it. This yields the basic RG differential equation

$$\left[\mu {\partial }_{\mu }+\sum _e{\tilde{\mathcal{D}}}_{\mu }e{\partial }_e+N{\gamma }_{\varphi }\right]W_N^R(e,\mu ,\epsilon ,\cdots )=0,{\gamma }_{\varphi }\equiv {\tilde{\mathcal{D}}}_{\mu }ln{Z}_{\varphi }$$

(65)

where the operator ${\tilde{\mathcal{D}}}_{\mu }$ is expressed in the variables $\mu ,e.$ The coefficients ${\tilde{\mathcal{D}}}_{\mu }e$ and ${\gamma }_{\varphi }$ are called the RG functions and are calculated in terms of various renormalization constants Z. Coupling constants (charges) g are those parameters e, on which the renormalization constants $Z=Z\left(g\right)$ depend. Logarithmic derivatives of charges in Equation (65) are $\beta$ functions

$${\beta }_g={\tilde{\mathcal{D}}}_{\mu }g.$$

(66)

All the RG-functions are UV-finite, i.e., have no poles in $\epsilon$, which is a consequence of the functions $W_N^R$ being UV-finite in Equation (65).

For models considered in the present work, the analysis of divergences should be augmented by the following considerations:

• For any dynamic model in Equation (1), all 1PI Green functions containing only the original fields $\varphi$ are proportional to the closed loops of step functions, hence they vanish, and thus do not generate counterterms.
• If for some reason several external momenta or frequencies occur as an overall factor in all the Feynman diagrams of a particular Green function, the real degree of divergence ${\delta }^{\prime }$ is less than $\delta \equiv d_{\mathsf{\Gamma }}(\epsilon =0)$ by the corresponding number of units.
• Sometimes the divergences formally allowed by dimensionality are absent due to symmetry restrictions, for instance, the Galilean invariance of the fully developed turbulence [31] restricts the form of possible counterterms.
• Nonlocal terms of the model are not renormalized.

In principle, these general considerations permit determining all superficially divergent functions and to explicitly obtain the corresponding counter-terms for any dynamic model.

The most convenient scheme for analytic calculations is the scheme of minimal subtractions (MS) proposed in [63], in which all the renormalization constants Z in the perturbation theory are of the form

$$Z^{\mathrm{MS}}(g,\epsilon )=1+\sum _{n=1}^{\infty }{g}^n\sum _{k=1}^n{\epsilon }^{-k}{c}_{n,k}.$$

(67)

In the dimensional renormalization the contribution to the coefficient of $g^n$ in Equation (67) may be expressed as a Laurent series in $\epsilon$. In the MS scheme, only the singular part of the Laurent expansion of each coefficient is retained. In any other renormalization scheme, the renormalization constant is of the form

$$Z(g,\epsilon )=1+\sum _{n=1}^{\infty }{g}^n\sum _{k=-n}^{\infty }{\epsilon }^k{c}_{n,-k},$$

(68)

where the regular part of each coefficient ${\sum }_{k=0}^{\infty }{\epsilon }^k{c}_{n,-k}$ is, by and large, an arbitrary regular function of $\epsilon$ at the origin.

5. Composite Operators and Operator Product Expansion

In this section, we recall the basic information about renormalization and critical exponents (dimensions) of composite operators, i.e., local products of the basic fields of the model and their derivatives. In the models we are interested in, they are constructed from the velocity field $v,$ scalar field $\theta$ or magnetic field $\mathit{b}$ at the single space-time point $x\equiv (t,\mathit{x})$. Examples are ${\mathit{v}}^n,$ ${\mathit{b}}^n,$ ${\theta }^n,$ ${\partial }_t{\mathit{v}}^n,$$\mathit{v}\Delta \mathit{v}$, ${(\nabla \theta ·\nabla \theta )}^n$ and so on.

A theoretical analysis of composite operators and their renormalization is important at least for two reasons. First, their critical dimensions and correlation functions can be measured experimentally and for some operators such data are available [64,65]. For instance, in the fully developed turbulence, the mean of the energy dissipation is proportional to the statistical average of the composite operator $v\Delta v$. This quantity enters the equation of energy balance and contributes to the redistribution of the energy of the turbulent motion and its dissipation. Moreover, strong statistical fluctuations of the operator of energy dissipation seem to account for deviations from Kolmogorov’s exponents and lead to the intermittency (multifractality) of the turbulent processes [23]. Second, the general solution of the RG equation (Equation (65)) contains an unknown scaling function depending on dimensionless effective variables (coupling constants, viscosity, etc.). This function can be calculated in the framework of usual perturbation scheme in an expansion parameter but, as mentioned above, in certain asymptotic ranges of scales, this calculation fails. Both experimental and theoretical reasons in theory of turbulence motivate us to study behavior of correlation functions with respect to outer (integral) scale L. Let us elucidate this issue in some detail. As an example we consider pair correlation function for velocity fluctuations $W_2=\langle vv\rangle$ for field-theoretic model in Equation (24). There is no field renormalization in this model [32], therefore the Green function $W_{2R}$ coincide with the unrenormalized function $W_2.$ The only difference lies in a choice of variable and perturbation theory (expansion either in charge g or $g_0$, respectively). In renormalized variable correlation, function $W_2$ depends on $k,g,\nu ,\mu$ and L. From a dimensional consideration, we directly see that $W_2$ can be represented in the form

$$W_2={\nu }^2{k}^{2-d}R(s,g,u),s=k/\mu ,u=kL,$$

(69)

where R is a function of dimensionless parameters and for brevity we have not explicitly written the transverse projection operator. The correlation function $W_2$ satisfies a general RG equation with ${\gamma }_{\varphi }=0$, which is a direct consequence of absence of renormalization of velocity field $v$, and reads

$$D_{RG}{W}_2=0,D_{RG}=\mu {\partial }_{\mu }+\beta \left(g\right){\partial }_g-{\gamma }_{\nu }\left(g\right)\nu {\partial }_{\nu }.$$

(70)

The solution of this equation can be found using the method of characteristics and presented in the form

$$W_2={\overline{\nu }}^2{k}^{2-d}R(1,\overline{g},\overline{u}),\overline{u}=u,$$

(71)

where $\overline{g},$$\overline{\nu }$ are invariant variables, i.e., the first integrals of Equation (70). Using standard RG considerations, the invariant viscosity [31] takes the following form

$$\overline{\nu }=\nu exp\left[{\int }_{\overline{g}}^g\mathrm{d}x\frac{{\gamma }_{\nu }\left(x\right)}{\beta \left(x\right)}\right]={\left(\frac{g{\nu }^3}{\overline{g}{s}^{2\epsilon }}\right)}^{1/3}={\left(\frac{D_0}{\overline{g}{k}^{2\epsilon }}\right)}^{1/3}.$$

(72)

As the parameter s approaches zero, invariant charge $\overline{g}$ approaches IR fixed point $g_{*}$ and $\overline{\nu }\to {\nu }_{*}={\left(D_0/g_{*}\right)}^{1/3}{k}^{-2\epsilon /3}$. Hence, at fixed point $g_{*}$ (far from dissipation scales $k\ll \mu \sim l^{-1}$), the single-time correlation function of velocity field takes the scaling form

$$W_2={\left(D_0/g_{*}\right)}^{2/3}{k}^{2-d-4\epsilon /3}R(1,g_{*},kL).$$

(73)

Setting $\epsilon =2$ gives kinetic energy spectrum $E\left(k\right)=W_2{k}^{d-1}$ that behaves as a power-law function of wavevector k. This coincides with Kolmogorov’s prediction $-5/3$ for the exponent. The remaining scaling function R is not determined yet and in general it is possible to employ perturbation theory and obtain infinite series in parameter $\epsilon$. In particular, in the theory of turbulence, the main interest is in the scaling function $R(1,g^{*},kL)$ in the inertial interval $kL\gg 1$. In the theory of critical phenomena, the asymptotic form of scaling functions for $kL\gg 1$ (formally, $L\to \infty$) is studied using Wilson’s operator product expansion (OPE) (see, e.g., [6,66]). The analog of L in turbulence is played by the correlation length $r_c$ in critical phenomena. It turns out that this technique can be used also in the theory of turbulence and in simplified (toy) models associated with the genuine turbulence (see, e.g., [5,31,33,67,68]).

The generating functional of the correlation functions of the field $\varphi$ with one insertion of the composite operator $F\left(\varphi \right)$ has the form (compare with the generating functional in Equation (7) for the usual correlation functions of $\varphi$)

$$\mathcal{G}[A,F]=\int \mathcal{D}\phi F\left(\phi \right)exp\left[\mathcal{S}\left[\varphi \right]+A\varphi \right].$$

(74)

Since the arguments of the fields in the operator F coincide (giving rise to new closed loops in the Feynman diagrams), correlation functions with these operators contain new UV divergences, which have to removed by an additional renormalization procedure (see, e.g., [5,6,66]). The standard RG equations yield the IR scaling of the renormalized correlation functions with definite critical dimensions ${\Delta }_F\equiv \Delta \left[F\right]$ of a set of basis operators F. Due to the renormalization, $\Delta \left[F\right]$ is not the sum of critical dimensions of the fields and derivatives in F. A detailed analysis of the renormalization procedure of composite operators for the stochastic NS problem can be found in the review [67], and below we restrict ourselves to the necessary information only.

As a rule, composite operators are mixed during the renormalization procedure, i.e., an UV finite renormalized operator $F^R$ (correlation functions with one insertion of $F_R$ do not possess UV divergences) takes the form $F^R=F+$ counterterms, in which “counterterms” stands for a linear combination of the operator F itself and other unrenormalized operators mixing the the operator F. Let $F\equiv \left\{F_{\alpha }\right\}$ denote a closed set of operators mixing only with each other under renormalization. For this set, the matrix of renormalization constants $Z_F\equiv \left\{Z_{\alpha \beta }\right\}$ and the matrix of anomalous dimensions ${\gamma }_F\equiv \left\{{\gamma }_{\alpha \beta }\right\}$ are defined by

$$F_{\alpha }=\sum _{\beta }{Z}_{\alpha \beta }{F}_{\beta }^R,{\gamma }_F=Z_F^{-1}{\tilde{\mathcal{D}}}_{\mu }{Z}_F.$$

(75)

The subsequent matrix of critical dimensions ${\Delta }_F\equiv \left\{{\Delta }_{\alpha \beta }\right\}$ reads

$$\Delta \left[F\right]\equiv {\Delta }_F=d_F^k+{\Delta }_{\omega }{d}_F^{\omega }+{\gamma }_F^{*},$$

(76)

in which $d_F^k$, $d_F^{\omega }$, and $d_F$ denote diagonal matrices of canonical dimensions of the operators of the closed set (the diagonal element corresponding to a particular operator F is equal to the sum of canonical dimensions of all fields, their derivatives and renormalized parameters in F) and ${\gamma }_F^{*}\equiv {\gamma }_F\left(g^{*}\right)$ is the matrix in Equation (75) at the fixed point.

Critical dimensions of the set $F\equiv \left\{F_{\alpha }\right\}$ correspond to the eigenvalues of the matrix ${\Delta }_F$. The basis operators possessing definite critical dimensions are linear combinations of the renormalized operators

$${\overline{F}}_{\alpha }^R=\sum _{\beta }{U}_{\alpha \beta }{F}_{\beta }^R,$$

(77)

where the matrix $U_F=\left\{U_{\alpha \beta }\right\}$ is such that the matrix ${\Delta }_F^{\prime }=U_F{\Delta }_F{U}_F^{-1}$ is diagonal.

Counterterms generated by a given operator F are determined by all possible 1PI Green functions with one insertion of operator F and an arbitrary number of primary fields $\varphi$,

$${\mathsf{\Gamma }}_{N;F}=\langle F(t,\mathit{x})\varphi (t,{\mathit{x}}_1)\cdots \varphi (t,{\mathit{x}}_N)\rangle .$$

(78)

The total canonical dimension (the formal degree of divergence) for these functions is given by

$$d_{\mathsf{\Gamma }}=d_F-N_{\mathsf{\Phi }}{d}_{\mathsf{\Phi }},$$

(79)

where the sum is taken over all types of field arguments. For $d_{\mathsf{\Gamma }}$ is a nonnegative integer.

According to the OPE, the single-time product $F_1(t,{\mathit{x}}_1)F_2(t,{\mathit{x}}_2)$ of two renormalized operators at $\mathit{x}\equiv ({\mathit{x}}_1+{\mathit{x}}_2)/2=const.$, and $\mathit{r}\equiv {\mathit{x}}_1-{\mathit{x}}_2\to 0$ can be represented as follows

$$F_1(t,{\mathit{x}}_1)F_2(t,{\mathit{x}}_2)=\sum _{\alpha }{A}_{\alpha }\left(\mathit{r}\right){\overline{F}}_{\alpha }^R(t,\mathit{x}).$$

(80)

Here, the functions $A_{\alpha }$ are the Wilson coefficients regular in L, whereas ${\overline{F}}_{\alpha }^R$ are all possible renormalized local composite operators of the type in Equation (77) allowed by symmetry arguments, with specific critical dimensions ${\Delta }_{{\overline{F}}_{\alpha }^R}$.

The renormalized correlator $\langle F_1(t,{\mathit{x}}_1)F_2(t,{\mathit{x}}_2)\rangle$ is obtained by averaging Equation (80) with the weight $exp{S}_R$, quantities $\langle {\overline{F}}_{\alpha }^R\rangle \propto L^{-d_{\alpha }}{f}_{\alpha }(g,L\mu )$ involving dimensionless (scaling) functions $f_{\alpha }(g,L\mu )$ appear on the right hand side. Their asymptotic behavior for $L\mu \to 0$ is found from the corresponding RG equations (see [34] for the case of Kraichnan model) and has the form

$$\langle {\overline{F}}_{\alpha }^R\rangle \propto L^{{\Delta }_{{\overline{F}}_{\alpha }^R}}.$$

(81)

From the operator product expansion in Equation (80), we therefore get

$$\langle F_1(t,{\mathit{x}}_1)F_2(t,{\mathit{x}}_2)\rangle =\sum _{{\overline{F}}^R}{C}_{{\overline{F}}^R}{(r/L)}^{{\Delta }_{{\overline{F}}^R}},r/L\to 0,$$

(82)

where the quantities $C_{{\overline{F}}^R}$ generated by the Wilson coefficients $A_{\alpha }$ in Equation (80) are regular in L, the summation is carried out over all possible composite renormalized basis operators ${\overline{F}}^R$ allowed by the symmetry of the left side, and ${\Delta }_{{\overline{F}}^R}$ are their critical dimensions. The leading contributions for $r/L\to 0$ are those with the least dimension ${\Delta }_{{\overline{F}}^R}.$ In the theory of critical phenomena, it is observed that all the nontrivial composite operators have positive critical dimensions ${\Delta }_{{\overline{F}}^R}>0$ for small $\epsilon$ and the most important term in Equation (82) corresponds to the simplest operator ${\overline{F}}^R=1$ with ${\Delta }_{{\overline{F}}^R}=0$, i.e., the function $R(r/L)$ is finite as $L\equiv r_c\to 0$ (see [6]). However, as has been noted in [68] in the model of developed turbulence composite operators with negative critical dimensions exist and are responsible for possible singular behavior of the scaling functions such as N-point correlation functions $W_N=\langle \phi \dots \phi \rangle$ as $r/L\to 0$. We call operators with ${\Delta }_{{\overline{F}}^R}<0$—if they exist—dangerous [67]. This is motivated by the fact that they correspond to contributions to Equation (82) which diverge for $r/L\to 0$. The scaling functions in Equation (82) decomposed in dangerous operators exhibit anomalous scaling behavior which is a manifestation of a nontrivial multifractal (intermittent) nature of the statistical fluctuations of the random fields under consideration and globally all the physical system.

Dangerous composite operators in the stochastic model of turbulence occur only for finite values of the RG expansion parameter $\epsilon$. Let us note that within the $\epsilon$ expansion it is not possible to determine whether or not a given operator is dangerous, if only its critical dimension is not found exactly employing the Schwinger equations, etc. or the Galilean symmetry (see [67,69]). Furthermore, dangerous operators appear in the operator product expansion in the form of infinite families with the spectrum of critical dimensions unbounded from below. Therefore, for a proper analysis of the large L behavior, a summation of their contributions is called for.

6. Schwinger Equations and Conservation Laws

Useful information about composite operators can be gained even without an actual calculation of Feynman diagrams. Exploiting invariance properties of functional integrals provides nontrivial relations known as Schwinger equations [47]. One of the simplest symmetries is the translation invariance of action in Equation (3). It is invariance with respect to a shift $\varphi \to \varphi +\omega$, where $\omega =\omega \left(x\right)$ is a suitably chosen function that vanishes sufficiently fast, i.e., $\omega (\infty )=0$. Such translations do not change the integration measure and as a result the quantity $\int \mathcal{D}\phi F[\varphi +\omega ]$ does not depend on $\omega$ for any functional $F\left[\varphi \right]$. We then easily derive that the first variation with respect to $\omega$ yields a formal relation

$$\int \mathcal{D}\varphi \frac{\delta F\left[\varphi \right]}{\delta \varphi }=0$$

(83)

written in the notation of Equation (7). The following relation is of particular importance

$$\int \mathcal{D}\varphi \frac{\delta }{\delta \varphi }exp[\mathcal{S}\left[\varphi \right]+A\varphi ]=0.$$

(84)

Performing variational derivatives gives us

$$\int \mathcal{D}\varphi \left[\frac{\delta \mathcal{S}\left[\varphi \right]}{\delta \varphi }+A\left(x\right)\right]{\mathrm{e}}^{\mathcal{S}\left[\varphi \right]+A\varphi }=0.$$

(85)

Multiplication by field ${\varphi }_{\alpha }$ inside the functional integral is tantamount to a differentiation with respect to the corresponding source field A. This observation allows us to rewrite Equation (85)

$$\left[\frac{\delta S\left[\varphi \right]}{\delta \varphi }{|}_{\varphi \to \delta /\delta A}+A\left(x\right)\right]\mathcal{G}\left(A\right)=0.$$

(86)

Substituting $\mathcal{G}={\mathrm{e}}^{\mathcal{W}}$, we obtain the corresponding Schwinger equation for $\mathcal{W}\left[A\right]$ where from we can derive the equation for $\mathsf{\Gamma }\left[\alpha \right].$ All these equations are of finite order (for polynomial action) in functional derivatives, and each of them is tantamount to an infinite chain of connected equations for the Green functions—the expansion coefficients of the corresponding functionals [47].

In the following discussion, we need one additional relation that corresponds to the Schwinger equation

$$\int \mathcal{D}\varphi \frac{\delta }{\delta {\phi }^{\prime }}(\phi \left(x\right)exp[\mathcal{S}\left[\varphi \right]+A\varphi ])=0,$$

(87)

where $\varphi$ stands for either the fluctuating field or the corresponding response field.

As discussed in Section 5, composite operators are related to experimentally measurable physical quantities. We illustrate this claim on an example of stochastic hydrodynamics summarized in the field-theoretic action in Equation (24). Our aim is to elucidate transfer of energy in a stationary of turbulent state. The latter condition ensures that time derivatives of averaged quantities are identically zero.

Let us derive equations describing energy and momentum transfer in turbulent flows. To obtain an equation expressing momentum conservation, we employ the first equation in Equation (87), where $\varphi$ consists of altogether two d-dimensional vector fields $\{v,v^{\prime }\}$. First, for $\varphi$, we choose a response field $v^{\prime }$, and we get

$$\int \mathcal{D}\varphi \frac{\delta }{\delta v_i^{\prime }\left(x\right)}exp[{\mathcal{S}}_{\mathrm{NS}}\left[\varphi \right]+A\varphi ]=\int \mathcal{D}\varphi \left[\frac{\delta {\mathcal{S}}_{\mathrm{NS}}\left[\varphi \right]}{\delta v_i^{\prime }\left(x\right)}+A_{v_i^{\prime }}\right]exp[{\mathcal{S}}_{\mathrm{NS}}\left[\varphi \right]+A\varphi ]=0.$$

(88)

Performing indicated derivative, we obtain differential equation

$$\langle A_{v_i^{\prime }}+D_{ij}{v}_j^{\prime }-{\partial }_t{v}_i+{\nu }_0\Delta v_i-(\mathit{v}·\nabla )v_i-{\partial }_{i}p\rangle =0.$$

(89)

Due to transversality of response field ${\mathit{v}}^{\prime }$, a nonlocal term has appeared, which corresponds to the pressure fluctuations

$$p=-\frac{{\partial }_l{\partial }_s}{\Delta }\left(v_l{v}_s\right).$$

(90)

To derive an equation describing energy transfer, we utilize the second Schwinger equation (Equation (87)). In particular, letting ${\phi }^{\prime }\to {\mathit{v}}^{\prime }$ and $\phi \to \mathit{v}$ in Equation (87) yields

$$\int \mathcal{D}\varphi \frac{\delta }{\delta v_i^{\prime }\left(x\right)}\left(v_i\left(x\right){\mathrm{e}}^{{\mathcal{S}}_{\mathrm{NS}}\left[\varphi \right]+A\varphi }\right)=0.$$

(91)

In an analogous manner to Equation (89), we get

$$\langle v_i{A}_{v_i^{\prime }}+v_i{D}_{ij}{v}_j^{\prime }-v_i{\partial }_t{v}_i+{\nu }_0{v}_i{\nabla }^2{v}_i-v_i\left(v_j{\partial }_j\right)v_i-\left(v_j{\partial }_j\right)p\rangle =0$$

(92)

written in a component notation. Note that all quantities have been normalized to unit mass, i.e., density has been set to unity ($\rho =1$). Equations (89) and (92) represent conservation laws for momentum and energy. They can be further rewritten into a physically more transparent form

$$\begin{array}{cc}\hfill \langle {\partial }_t{v}_i+{\partial }_j{\Pi }_{ij}\rangle & =\langle D_{ij}{v}_j^{\prime }\rangle +A_{v_i^{\prime }},\hfill \end{array}$$

(93)

$$\begin{array}{cc}\hfill \langle {\partial }_{t}E+{\partial }_i{S}_i\rangle & =\langle -\mathcal{E}+v_i{D}_{ij}{v}_j^{\prime }+v_i{A}_{v_i^{\prime }}\rangle ,\hfill \end{array}$$

(94)

where $v_i$ might be interpreted as momentum density, $E\equiv v^2/2$ is energy density, ${\Pi }_{ik}$ is tensor of momentum transfer, $S_i$ is vector of energy flow, and $\mathcal{E}$ is a rate of energy dissipation. Direct comparison of Equations (89) and (92) with Equations (93) and (94) yields explicit expressions

$$\begin{array}{cc}\hfill {\Pi }_{ij}& =p{\delta }_{ij}+v_i{v}_j-{\nu }_0({\partial }_i{v}_j+{\partial }_j{v}_i),\hfill \end{array}$$

(95)

$$\begin{array}{cc}\hfill S_i& =p{v}_i-{\nu }_0{v}_j({\partial }_i{v}_j+{\partial }_j{v}_i)+\frac{1}{2}{v}_i{\mathit{v}}^2,\hfill \end{array}$$

(96)

$$\begin{array}{cc}\hfill \mathcal{E}& =\frac{1}{2}{\nu }_0{({\partial }_i{v}_j+{\partial }_j{v}_i)}^2.\hfill \end{array}$$

(97)

We recognize Equation (93) as a stochastic Navier–Stokes equation stirred by random force $D_{ij}{v}_j^{\prime }$ and regular force $A_{v_i^{\prime }}$.

Functional averaging of Equations (93) and (94) according to the prescription in Equation (6) with weigh functional $exp{\mathcal{S}}_{\mathrm{NS}}\left[\varphi \right]$ leads to the balance equation for energy and momentum. Assuming vanishing external force $A_{v_i^{\prime }}$, we obtain the following equation for time derivative of energy

$${\partial }_t〈E〉=-〈\mathcal{E}〉-{\partial }_i〈S_i〉+〈v_i{D}_{ij}{v}_j^{\prime }〉.$$

(98)

It is clear that for homogeneous and isotropic flows at zero external force the mean value $〈F\left(x\right)〉$ of an arbitrary composite operator $F\left(x\right)=F(t,\mathit{x})$ could not depend on the position $\mathit{x}$. Hence, it is constant and consequently all spatial derivatives are identically zero. From Equation (98), we then get for steady state

$$\overline{\mathcal{E}}=\int \mathrm{d}x{D}_{ij}(x,x^{\prime })〈v_i\left(x\right)v_j^{\prime }\left(x^{\prime }\right)〉,$$

(99)

where we have introduced the following abbreviation for mean energy dissipation

$$\overline{\mathcal{E}}\equiv 〈\mathcal{E}〉.$$

(100)

Let us recall (see Equation (19)) that pair correlation for random force can be written as $D_{ij}(x,x^{\prime })=\delta (t-t^{\prime })d_{ij}(\mathit{x},{\mathit{x}}^{\prime }).$ Insertion of this relation into Equation (100) and integrating over time variable t yields

$$\overline{\mathcal{E}}=\int {\mathrm{d}}^{d}x{d}_{ij}(\mathit{x},{\mathit{x}}^{\prime })〈v_i(t,\mathit{x})v_j^{\prime }(t,{\mathit{x}}^{\prime })〉{|}_{t=t^{\prime }}.$$

(101)

To the lowest order in perturbation theory, the retarded response function ${〈v_i(t,\mathit{x})v_j^{\prime }(t,{\mathit{x}}^{\prime })〉}_0$ is $\delta$-correlated in spatial variable. This property holds also for the full response function $〈v_i\left(x\right)v_j^{\prime }\left(x^{\prime }\right)〉$ (what follows from a straightforward analysis of Feynman graphs) and therefore

$$〈v_i\left(x\right)v_j^{\prime }\left(x^{\prime }\right)〉{|}_{t=t^{\prime }}=\frac{1}{2}\delta (\mathit{x}-{\mathit{x}}^{\prime })\left({\delta }_{ij}-\frac{{\partial }_i{\partial }_j}{\Delta }\right).$$

(102)

Insertion of this relation into Equation (101), recalling Equation (20) and integrating through spatial variable ${\mathit{x}}^{\prime }$, we derive

$$\overline{\mathcal{E}}=\frac{1}{2}\int \frac{{\mathrm{d}}^{d}k}{{\left(2\pi \right)}^d}{d}_f\left(k\right)P_{ij}\left(\mathit{k}\right)P_{ij}\left(\mathit{k}\right).$$

(103)

Summation over internal indices i and j corresponds to a calculation for a trace of transverse projection operator $P_{ii}\left(\mathit{k}\right)$, which equals $d-1$. Thus, we finally arrive at the expression

$$\overline{\mathcal{E}}=\frac{d-1}{2{\left(2\pi \right)}^d}\int {\mathrm{d}}^{d}k{d}_f\left(k\right).$$

(104)

This relation reflects an important property of stationary homogeneous turbulence. It expresses the expected fact that, to achieve a stationary state, it is necessary to inject energy into a system in a continuous fashion. In the stochastic approach, this is done through a random force, which compensates energetic losses due to friction processes. These losses are expressed through mean energy dissipation $\overline{\mathcal{E}}$.

7. Ward Identities and Galilean Invariance

We say that a given theory possesses a symmetry, if the corresponding action functional of the theory is zero under the action generating this symmetry. Ward–Takahashi identities mathematically express inherent symmetry of a given field-theoretic action. In stochastic models of turbulence, they correspond to the well-known Galilean invariance. As a rule, these identities provide nontrivial relations between various Green functions of theory and consequently between renormalization constants. Moreover, they are also relevant for an analysis of composite operators.

As pointed out in Section 4, certain divergences present in Feynman diagrams of the perturbative expansion of a Green function might cancel each other, so that the given Green function is in fact UV finite. This compensation might be caused by the inherent fact that the underlying stochastic model describing developed turbulence is invariant with respect to the Galilean transformations. Of course, such and similar mechanisms are quite general in physics. As a further example, we can mention absence of potential UV divergences in quantum electrodynamics, or even quantum chromodynamics describing strong interactions between quarks and gluons. The underlying symmetry in these cases is gauge symmetry.

Now, we show that Galilean invariance in the model in Equation (24) ensures that UV divergences in triple 1PI function

$${\langle {v_i}^{\prime }\left(x_1\right)v_j\left(x_2\right)v_k\left(x_3\right)\rangle }_{1\mathrm{PI}}\equiv {\mathsf{\Gamma }}_{i,jk}(x_1;x_2,x_3)$$

(105)

are actually absent. The essential idea consists in a derivation of certain Ward identity, which takes form of an differential equation relating the triple 1PI Green function ${\mathsf{\Gamma }}_{i;sl}$ with the pair 1PI Green function ${\mathsf{\Gamma }}_{i;j}$, which is an abbreviation for

$${\langle v_i^{\prime }\left(x_1\right)v_j\left(x_2\right)\rangle }_{1\mathrm{PI}}\equiv {\mathsf{\Gamma }}_{i,j}(x_1;x_2).$$

(106)

These relations stem from invariance property of generating functional with respect to the Galilean transformations. Then, an absence of certain types of UV divergences in ${\mathsf{\Gamma }}_{i,j}$ leads to absence of divergences in ${\mathsf{\Gamma }}_{i,jk}$.

Let us consider a generalized Galilean transformation of fields $\varphi \equiv \{{\mathit{v}}^{\prime },\mathit{v}\}$ defined as follows

$$\varphi \to {\varphi }_w:{\mathit{v}}_w\left(x\right)=\mathit{v}\left(x_w\right)-\mathit{w}\left(t\right),{\mathit{v}}_w^{\prime }\left(x\right)={\mathit{v}}^{\prime }\left(x_w\right),$$

(107)

where

$${\mathit{x}}_w\equiv \mathit{x}+\mathit{u}\left(t\right);t=t^{\prime };\mathit{u}={\int }_{-\infty }^t\mathrm{d}{t}^{\prime }\mathit{w}\left(t^{\prime }\right)={\int }_{-\infty }^{\infty }\mathrm{d}{t}^{\prime }\theta (t-t^{\prime })\mathit{w}\left(t^{\prime }\right).$$

(108)

Here, $\mathit{w}\left(t\right)$ is some velocity vector describing the Galilean transformation. The spatial vector $\mathit{u}\left(t\right)$ is responsible for a shift of spatial coordinate $\mathit{x}$. The transformations in Equation (107) are a generalization of the standard Galilean transformations [70], in which velocity $w$ is constant in time variable. This is brought about by functional integration in which functional space has to be restricted by appropriate conditions. Here, it is required that velocity and response fields $\mathit{v}$ and ${\mathit{v}}^{\prime }$ vanish sufficiently quickly in the limit $\left|t\right|\to \infty$. Of course, arbitrary symmetry transformation of the model must comply with this property.

Insertion of Equation (107) into the action in Equation (24) yields the following relation for the transformed action

$${\mathcal{S}}_{\mathrm{NS}}\left[{\varphi }_w\right]={\mathcal{S}}_{\mathrm{NS}}\left[\varphi \right]+{\mathit{v}}^{\prime }·{\partial }_t\mathit{w}={\mathcal{S}}_{\mathrm{NS}}\left[\varphi \right]-\left({\partial }_t{\mathit{v}}^{\prime }\right)·\mathit{w},$$

(109)

where in the last equation we have transformed time derivative using partial integration. In this derivation, the following relations for variational derivatives have been utilized

$$\begin{array}{cc}\hfill {\delta }_w\mathit{v}\left(x\right)& =(\mathit{u}·\nabla )\mathit{v}\left(x\right)-\mathit{w},\hfill \end{array}$$

(110)

$$\begin{array}{cc}\hfill {\delta }_w{\mathit{v}}^{\prime }\left(x\right)& =(\mathit{u}·\nabla ){\mathit{v}}^{\prime }\left(x\right),\hfill \end{array}$$

(111)

$$\begin{array}{cc}\hfill {\delta }_w{\partial }_t\mathit{v}\left(x\right)& =(\mathit{u}·\nabla ){\partial }_t\mathit{v}\left(x\right)+\left(\mathit{w}·\nabla \right)\mathit{v}\left(x\right)-{\partial }_t\mathit{w},\hfill \end{array}$$

(112)

which can be directly obtained from Equation (107). In infinitesimal form, Equation (109) takes form

$${\delta }_w{\mathcal{S}}_{\mathrm{NS}}\left[\varphi \right]=-\left({\partial }_t{\mathit{v}}^{\prime }\right)·\mathit{w},$$

(113)

where

$${\delta }_w{\mathcal{S}}_{\mathrm{NS}}\left[\varphi \right]\equiv {\mathcal{S}}_{\mathrm{NS}}\left[{\varphi }_w\right]-{\mathcal{S}}_{\mathrm{NS}}\left[\varphi \right]=\mathit{w}·\frac{\delta {\mathcal{S}}_{\mathrm{NS}}\left[{\varphi }_w\right]}{\delta \mathit{w}}{|}_{w=0}.$$

(114)

The implicit assumption in derivation of Equation (113) is smallness of velocity $\mathit{w}$.

In a compact form, the requirement of Galilean invariance for the model in Equation (24) is equivalent to the condition

$$\mathcal{G}\left[A\right]=\mathcal{G}\left[A_w\right]$$

(115)

or in an infinitesimal form

$${\delta }_w\mathcal{G}\left[A\right]=0,\mathit{w}·\frac{\delta \mathcal{G}\left[A_w\right]}{\delta \mathit{w}}{|}_{w=0}=0.$$

(116)

The Ward identities are useful not only for Green functions of basic fields $\mathit{v},{\mathit{v}}^{\prime }$, but for composite operators as well (see Section 5). This motivates an introduction of generalized generating functional that includes composite operators F. It can be presented in the following form

$$\mathcal{G}[A,a]=\int \mathcal{D}\varphi exp\left[{\mathcal{S}}_{\mathrm{NS}}\left[\varphi \right]+A\varphi +aF\left(\varphi \right)\right],$$

(117)

where $aF\left(\varphi \right)$ stands for $aF\left(\varphi \right)\equiv {\sum }_i^N\int \mathrm{d}x{a}_i\left(x\right)F_i(\varphi ,x)$. Sources $a_i\left(x\right)$ correspond to N composite operators $F_i(\varphi ,x)$. In principle, the functional in Equation (117) generates not only all possible Green functions of basic fields, but also full Green functions consisting of arbitrary inclusion of composite operators and fields. In this regard, it is useful to compare this generating functional in Equation (117) with the functional in Equation (74) that generates Green functions containing one inclusion of composite operator F.

The functional measure is obviously translationally invariant. Hence, equality $\mathcal{D}\varphi =\mathcal{D}{\varphi }_v$ is satisfied and Equation (116) can be rewritten in the following way

$$\int \mathcal{D}\varphi {\delta }_w{\mathrm{e}}^{{\mathcal{S}}_{\mathrm{NS}}\left({\varphi }_w\right)+A{\varphi }_w+a{F}_w}=0,$$

(118)

or

$$\int \mathcal{D}\varphi \left\{{\delta }_w{\mathcal{S}}_{\mathrm{NS}}\left[\varphi \right]+A{\delta }_w\varphi +a{\delta }_{w}F\right\}{\mathrm{e}}^{{\mathcal{S}}_{\mathrm{NS}}\left[\varphi \right]+A\varphi +aF}=0.$$

(119)

As source fields A and a are independent, the choice $a=0$ yields Ward identities for Green functions of basic fields, whereas $a\ne 0$ leads to additional Ward identities for Green functions containing contributions from composite operators. Further, we concentrate on derivation of Ward identity for Green function containing solely basic fields $\varphi =\{\mathit{v},{\mathit{v}}^{\prime }\}$, which we derive from Equation (119) inserting $a=0$. We have

$$〈〈-\mathit{w}·{\partial }_t{\mathit{v}}^{\prime }+A{\delta }_w\varphi 〉〉=0,$$

(120)

where double brackets $〈〈\cdots 〉〉$ correspond to the functional averaging with respect to the weight functional $exp\left[\mathcal{S}\right[\varphi ]+A\varphi ]$

$$〈〈\cdots 〉〉=\frac{\int \mathcal{D}\varphi \cdots exp\left[{\mathcal{S}}_{\mathrm{NS}}\left[\varphi \right]+A\varphi \right]}{\int \mathcal{D}\varphi exp\left[{\mathcal{S}}_{\mathrm{NS}}\left[\varphi \right]+A\varphi \right]}.$$

(121)

Using the relations in Equations (108)–(111), we rewrite Equation (120) in the component notation

$$\int \mathrm{d}x{w}_i\left(t\right)〈〈-{\partial }_t{v}_i^{\prime }\left(x\right)+\underset{-\infty }{\overset{t}{\int }}\mathrm{d}{t}^{\prime }\left[A_j\left(x^{\prime }\right){\partial }_i{v}_j\left(x^{\prime }\right)+A_j^{\prime }\left(x^{\prime }\right){\partial }_i{v}_j^{\prime }\left(x^{\prime }\right)\right]-A_i\left(x\right)〉〉=0,$$

(122)

where for brevity we have denoted $x^{\prime }\equiv (t^{\prime },\mathit{x})$, ${\partial }_i=\partial /\partial x_i$, $A_j\left(x\right)\equiv A_{v_j}\left(x\right)$ and $A_j^{\prime }\left(x\right)\equiv A_{v_j^{\prime }}\left(x\right)$. As Galilean velocity $w$ is arbitrary, Equation (122) can be further simplified to

$$\int {\mathrm{d}}^{d}x〈〈-{\partial }_t{v}_i^{\prime }\left(x\right)+\underset{-\infty }{\overset{\infty }{\int }}\mathrm{d}{t}^{\prime }\theta (t-t^{\prime })\left[A_j\left(x^{\prime }\right){\partial }_i{v}_j\left(x^{\prime }\right)+A_j^{\prime }\left(x^{\prime }\right){\partial }_i{v}_j^{\prime }\left(x^{\prime }\right)\right]-A_i\left(x\right)〉〉=0.$$

(123)

Every term in Equation (123) containing fields $\varphi$ might be obtained by taking an appropriate number of derivatives of $exp\left[A\varphi \right]$ with respect to source A. This means that field $\varphi$ can be effectively replaced by variational derivative $\delta /\delta A$

$$\varphi in〈〈\dots 〉〉\iff \frac{\delta }{\delta A}.$$

(124)

Hence, we get

$$\int {\mathrm{d}}^{d}x〈〈-{\partial }_t\frac{\delta }{\delta A_i^{\prime }\left(x\right)}+{\int }_{-\infty }^{\infty }\mathrm{d}{t}^{\prime }\theta (t-t^{\prime })\left[A_j\left(x^{\prime }\right)\frac{\partial }{\partial x_i}\frac{\delta }{\delta A_j\left(x^{\prime }\right)}+A_j^{\prime }\left(x^{\prime }\right)\frac{\partial }{\partial x_i}\frac{\delta }{\delta A_j^{\prime }\left(x^{\prime }\right)}\right]-A_i\left(x\right)〉〉=0.$$

(125)

In this formulation, the whole expression in brackets does not depend on integration over fields $\varphi$ and the equation can be further rewritten in terms of generating functional

$$\int {\mathrm{d}}^{d}x\left\{-{\partial }_t\frac{\delta }{\delta A_i^{\prime }\left(x\right)}+{\int }_{-\infty }^{\infty }\mathrm{d}{t}^{\prime }\theta (t-t^{\prime })\left[A_j\left(x^{\prime }\right)\frac{\partial }{\partial x_i}\frac{\delta }{\delta A_j\left(x^{\prime }\right)}+A_j^{\prime }\left(x^{\prime }\right)\frac{\partial }{\partial x_i}\frac{\delta }{\delta A_j^{\prime }\left(x^{\prime }\right)}\right]-A_i\left(x\right)\right\}\mathcal{G}\left(A\right)=0.$$

(126)

Substituting $\mathcal{G}={\mathrm{e}}^{\mathcal{W}}$, we rewrite Equation (126) into equation for generating functional $\mathcal{W}$ for connected Green functions

$$\int {\mathrm{d}}^{d}x\left\{-{\partial }_t\frac{\delta \mathcal{W}}{\delta A_i^{\prime }\left(x\right)}+{\int }_{-\infty }^{\infty }\mathrm{d}{t}^{\prime }\theta (t-t^{\prime })\left[A_j\left(x^{\prime }\right)\frac{\partial }{\partial x_i}\frac{\delta \mathcal{W}}{\delta A_j\left(x^{\prime }\right)}+A_j^{\prime }\left(x^{\prime }\right)\frac{\partial }{\partial x_i}\frac{\delta \mathcal{W}}{\delta A_j^{\prime }\left(x^{\prime }\right)}\right]-A_i\left(x\right)\right\}=0.$$

(127)

Finally, by means of the functional Legendre transformations in Equation (57), we rewrite Equation (127) in terms of generating functional for 1PI (vertex) functions $\mathsf{\Gamma }$

$$\int {\mathrm{d}}^{d}x\left\{-{\partial }_t{\alpha }_i^{\prime }\left(x\right)+{\int }_{-\infty }^{\infty }\mathrm{d}{t}^{\prime }\theta (t-t^{\prime })\left[\frac{\delta \mathsf{\Gamma }\left(\alpha \right)}{\delta {\alpha }_j\left(x^{\prime }\right)}\frac{\partial {\alpha }_{v_j}\left(x^{\prime }\right)}{\partial x_i}+\frac{\delta \mathsf{\Gamma }\left(\alpha \right)}{\delta {\alpha }_j^{\prime }\left(x^{\prime }\right)}\frac{\partial {\alpha }_{v_j^{\prime }}\left(x^{\prime }\right)}{\partial x_i}\right]-\frac{\delta \mathsf{\Gamma }\left(\alpha \right)}{\delta {\alpha }_i\left(x\right)}\right\}=0,$$

(128)

where

$${\alpha }_i\equiv {\alpha }_{v_i},{\alpha }_i^{\prime }\equiv {\alpha }_{v_i^{\prime }}.$$

(129)

Generating functional $\mathsf{\Gamma }$ can be represented in form of formal Taylor series with respect to sources $\alpha$, in which coefficients by powers ${\alpha }^n$ ($n=0,1,2,\cdots ,\infty$) are 1PI Green functions. The first few terms are

$$\mathsf{\Gamma }\left(\alpha \right)={\alpha }_{v^{\prime }}{\mathsf{\Gamma }}_{v^{\prime }v}{\alpha }_v+\frac{1}{2}{\alpha }_{v^{\prime }}{\mathsf{\Gamma }}_{v^{\prime }vv}{\alpha }_v^2+\cdots ,$$

(130)

where, for instance, the second term actually stands for an expression

$${\alpha }_{v^{\prime }}{\mathsf{\Gamma }}_{v^{\prime }vv}{\alpha }_v^2\equiv \int \mathrm{d}{x}_1\int \mathrm{d}{x}_2\int \mathrm{d}{x}_3{\alpha }_{v_i^{\prime }}\left(x_1\right){\mathsf{\Gamma }}_{i,jk}(x_1;x_2,x_3){\alpha }_{v_j}\left(x_2\right){\alpha }_{v_k}\left(x_3\right).$$

(131)

Insertion of this expansion into Equation (128), we get an infinite system of the Ward identities that relates various 1PI Green functions. For the present discussion, the most relevant is the Ward identity between pair (two-point) and triple (three-point) Green functions. The substitution of the expansion in Equation (130) into Equation (128) yields a formal polynomial in variables ${\alpha }_i$ and ${\alpha }_j^{\prime }$. A comparison of terms proportional to term ${\alpha }_i{\alpha }_j^{\prime }$ gives the following relation

$$\int {\mathrm{d}}^{d}x{\mathsf{\Gamma }}_{i,jk}(x_1;x_2,x)+\left[\theta (t-t_1)\frac{\partial }{\partial x_{1k}}+\theta (t-t_2)\frac{\partial }{\partial x_{2k}}\right]{\mathsf{\Gamma }}_{i,s}(x_1;x_2)=0.$$

(132)

Due to translation, invariance, ${\mathsf{\Gamma }}_{i,j}(x_1;x_2)={\mathsf{\Gamma }}_{i,j}(x_1-x_2;0)$, therefore

$$\frac{\partial }{\partial x_{2k}}=-\frac{\partial }{\partial x_{1k}},k=1,\dots ,d.$$

(133)

Using this relation and integrating over time, the variable t (132) yields

$$\int \int {\mathrm{d}}^{d}x\mathrm{d}t{\mathsf{\Gamma }}_{i,jk}(x_1;x_2,x)+\int \mathrm{d}t\left[\theta (t-t_1)-\theta (t-t_2)\right]\frac{\partial }{\partial x_{1k}}{\mathsf{\Gamma }}_{i,j}(x_1;x_2)=0.$$

(134)

Let us integrate the second term in this equation over time t (${\mathsf{\Gamma }}_{i,j}(x_1;x_2)$ does not depend on time variable t, only on $t_1$ and $t_2$). It is straightforward to show that

$${\int }_{-\infty }^{\infty }\mathrm{d}t\left[\theta (t-t_1)-\theta (t-t_2)\right]={\int }_{t_1}^{\infty }\mathrm{d}t-{\int }_{t_2}^{\infty }\mathrm{d}t={\int }_{t_1}^{t_2}\mathrm{d}t=t_2-t_1.$$

(135)

Employing this relation in Equation (134), we finally obtain the required Ward identity in the coordinate representation

$$\int \mathrm{d}x{\mathsf{\Gamma }}_{i,jk}(x_1;x_2,x)+(t_2-t_1)\frac{\partial {\mathsf{\Gamma }}_{i,j}(x_1;x_2)}{\partial x_{1k}}=0.$$

(136)

For an analysis of UV divergences, it is useful to rewrite the Ward identity into a frequency– momentum (Fourier) representation. Let us determine Fourier transforms of pair and triple Green functions. They read

$${\mathsf{\Gamma }}_{i,j}(x_1;x_2)=\frac{1}{{\left(2\pi \right)}^{2d+2}}\int \mathrm{d}{p}_1\int \mathrm{d}{p}_2{\mathsf{\Gamma }}_{i,j}(p_1;p_2){\mathrm{e}}^{i({\mathit{p}}_1·{\mathit{x}}_1+{\mathit{p}}_2·{\mathit{x}}_2)-i({\omega }_1{t}_1+{\omega }_2{t}_2)},$$

(137)

where $p_i=({\omega }_i,{\mathit{p}}_i)$. Translation invariance in time and space variables dictates

$${\mathsf{\Gamma }}_{i;j}(p_1,p_2)={\left(2\pi \right)}^{d+1}\delta ({\omega }_1+{\omega }_2){\delta }^{\left(d\right)}({\mathit{p}}_1+{\mathit{p}}_2){\mathsf{\Gamma }}_{i,j}\left(p_1\right).$$

(138)

This allows us to simplify Equation (137) into

$${\mathsf{\Gamma }}_{i,j}(x_1;x_2)=\frac{1}{{\left(2\pi \right)}^{d+1}}\int \mathrm{d}p{\mathsf{\Gamma }}_{i;j}\left(p\right){\mathrm{e}}^{i\mathit{p}·({\mathit{x}}_1-{\mathit{p}}_2·{\mathit{x}}_2)-i{\omega }_p(t_1-t_2)}.$$

(139)

A similar consideration applies also for the triple Green function ${\mathsf{\Gamma }}_{i,jk}$ and leads to the Fourier representation

$$\begin{array}{cc}\hfill {\mathsf{\Gamma }}_{i,jk}(x_1;x_2,x_3)& =\int \int \int \frac{\mathrm{d}{p}_1\mathrm{d}{p}_2\mathrm{d}{p}_3}{{\left(2\pi \right)}^{2(d+1)}}\delta \left(\sum _i{\omega }_i\right){\delta }^{\left(d\right)}\left(\sum _i{\mathit{p}}_i\right){\mathsf{\Gamma }}_{i,jk}(p_1;p_2,p_3){\mathrm{e}}^{i{\sum }_{i=1}^3({\mathit{p}}_i{\mathit{x}}_i-{\omega }_i{t}_i)}.\hfill \end{array}$$

(140)

Hence, we derive

$$\int \mathrm{d}x{\mathsf{\Gamma }}_{i,jk}(x_1;x_2,x)=\frac{1}{{\left(2\pi \right)}^{d+1}}\int \mathrm{d}p{\mathsf{\Gamma }}_{i,jk}(p;-p,0){\mathrm{e}}^{i\mathit{p}({\mathit{x}}_1-{\mathit{x}}_2)-i{\omega }_p(t_1-t_2)}.$$

(141)

Insertion of the relations in Equations (139) and (141) into Equation (136) leads to the Ward identity in form

$${\mathsf{\Gamma }}_{i,jk}(p;-p,0)=p_l\frac{\partial }{\partial {\omega }_p}{\mathsf{\Gamma }}_{i,j}\left(p\right).$$

(142)

This relation can be conveniently represented in a graphical form as follows

(143)

From the Ward identity, it directly follows that there are no UV divergences in triple Green functions. In fact, on the right hand side of Equation (143), we have the Green function ${\mathsf{\Gamma }}_{i,j},$ which contains only divergent structures that are proportional to ${\mathit{p}}^2$ (divergences proportional to frequency ${\omega }_p$). Therefore, taking derivative with respect to frequency ${\omega }_p$ ensures that the divergent part on right hand side of Equation (143) vanishes. Hence, the right hand side of Equation (143) is UV finite and therefore also the left hand side of the given equation is UV finite. This concludes the proof of UV finiteness of the triple one-time 1PI Green function ${\mathsf{\Gamma }}_{i,jk}$.

8. Symmetry Restoration

In the previous Section 7, we show that the Galilean symmetry in Equations (107) and (108) restricts appearances of UV divergences in the perturbation theory. In this regard, a natural question arises: What can violation of the Galilean symmetry lead to? We address this issue in the context of the stochastic Navier–Stokes model in Equation (24). Let us note that the following exposition closely follows that in [71], where all necessary technical details can be found. Violation can be achieved by various means. In particular, it can be achieved by modification of time behavior of the force correlator in Equation (20). Let us imagine that there exists microscopic finite correlation time behaving according to Ornstein-Uhlenbeck process [17,72].

Without loss of generality, let us consider a generalization of Equation (23) as

$$〈f_i(\omega ,\mathit{k})f_j({\omega }^{\prime },{\mathit{k}}^{\prime })〉={\left(2\pi \right)}^{d+1}\delta (\omega +{\omega }^{\prime })\delta (\mathit{k}+{\mathit{k}}^{\prime })P_{ij}\left(\mathit{k}\right)d_f(\omega ,\mathit{k}),$$

(144)

where the kernel function $d_f$ now assumes the following form

$$\begin{array}{c}\hfill d_f({\omega }_k,\mathit{k})=D_0{P}_{ij}\left(\mathit{k}\right)\frac{k^{8-d-(y+2\eta )}}{{\omega }_k^2+{\nu }_0^2{u}_0^2{k}^{4-2\eta }}.\end{array}$$

(145)

Exponent $\eta$ is related to the dispersion law ${\omega }_k\propto k^{-2+\eta }$. From dimensional considerations, parameter $D_0$ can be represented as

$$D_0=g_0{\nu }_0^5{u}_0^2,$$

(146)

which can be interpreted as a defining relation for charge $g_0$. Quite general form of Equation (145) allows studying two special cases: the limit $u_0\to 0$ fully corresponds to the time-independent correlation of random force. On the other hand, the limit $u_0\to \infty$ yields the previously analyzed Galilean model in Equation (20).

The propagator in the frequency–momentum representation takes the form

$${〈v_i{v}_j〉}_0=D_0\frac{k^{8-d-(y+2\eta )}}{{\omega }^2+{\nu }_0^2{u}_0^2{k}^{4-2\eta }}\frac{P_{ij}\left(\mathit{k}\right)}{{\omega }^2+{\nu }_0^2{k}^4},$$

(147)

According to the theoretical consideration discussed in Section 2 and Section 4, it can be shown [71] that the model corresponding to Equation (24) with Equation (145) is multiplicatively renormalizable in which two renormalization constants $Z_1$ and $Z_2$ have to be added. The renormalized action functional is

$${\mathcal{S}}_{\mathrm{NS}}\left[\varphi \right]=\frac{1}{2}{v}_i^{\prime }{D}_{ik}{v}_k^{\prime }+v_k^{\prime }\left[-{\partial }_t-Z_1\left(v_i{\partial }_i\right)+Z_2\nu {\partial }^2\right]v_k,$$

(148)

where g, $\nu$, and u are the renormalized counterparts of the bare (original) parameters; the correlation function $D_v$ is written in terms of the renormalized parameters with the use of the relation $g_0{\nu }_0^5{u}_0^2=g{\mu }^{y+2\eta }{\nu }^5{u}^2$; and the reference scale $\mu$ is a new free parameter of the renormalized model. The action of the renormalized model in Equation (148) can be constructed from the original action in Equation (24) by renormalization of fields $v\to Z_{v}v$, $v^{\prime }\to Z_{v^{\prime }}{v}^{\prime }$ and parameters

$$g_0=g{\mu }^y{Z}_g,u_0=u{\mu }^{\eta }{Z}_u,{\nu }_0=\nu Z_{\nu }.$$

(149)

The renormalization constants in Equations (148) and (149) are related as follows

$$Z_{\nu }=Z_2,Z_u=Z_2^{-1},Z_g=Z_1^2{Z}_2^{-3},Z_v=Z_{v^{\prime }}^{-1}=Z_1.$$

(150)

The leading one-loop calculation [71] of Feynman diagrams is relatively straightforward and ensuing anomalous dimensions are

$$\begin{array}{c}\hfill {\gamma }_1=g\frac{1}{d(d+2)}\frac{u}{{(u+1)}^3},{\gamma }_2=g\frac{1}{d(d+2)}\frac{u^{3}d(d-1)+3u^{2}d(d-1)+2u(d^2-d+2)}{4{(u+1)}^3},\end{array}$$

(151)

where for brevity the following redefinition of the charge g

$$g\frac{S_d}{{\left(2\pi \right)}^d}\equiv g\overline{S_d}\to g.$$

(152)

has been used. Further, from the relations in Equation (150), we get $\beta$-functions

$${\beta }_g=g(-y-2{\gamma }_1+3{\gamma }_2),{\beta }_u=u(-\eta +{\gamma }_2).$$

(153)

A straightforward analysis of the $\beta$ functions reveals three IR fixed points: the trivial fixed point (zero values of charges’ coordinates) and two nontrivial fixed points. The trivial (Gaussian) fixed point, at which all nonlinearities are IR irrelevant, is

$$g^{*}=0,u^{*}=0$$

(154)

and is IR attractive when both y and $\eta$ are negative.

A further fixed point $\left\{g^{*},u^{*}\right\}$ with the coordinates

$$u^{*}=\frac{-3+\sqrt{1-\frac{16(\alpha -1)}{d(d-1)(3\alpha -1)}}}{2},g^{*}=d(d+2)\frac{{(u^{*}+1)}^3}{u^{*}}\frac{3\alpha -1}{2}y,$$

(155)

is actually of a saddle-point type for, which one eigenvalue of the matrix $\Omega$ is positive and the other negative.

The most interesting fixed point is actually obtained in the limit $u^{*}\to \infty$, from Equations (145) and (147). This case corresponds to an earlier studied model with random force uncorrelated in time [32]. Therefore, we might expect to find the known fixed point of this model. Introducing a new variable $x=1/u$, one obtains a new $\beta$ function of the following form

$${\beta }_x={\tilde{\mathcal{D}}}_{\mu }x=-\frac{1}{u^2}{\beta }_u=x\left[\eta -g\frac{d-1}{4(d+2)}\right].$$

(156)

Actual calculations reveal the existence of the fixed point with the coordinates

$$x^{*}=0,g^{*}=\frac{4(d+2)}{3(d-1)}y,$$

(157)

which coincides with that in Ref. [32] and is IR attractive in the region restricted by inequalities $y>0$ and $\eta >y/3$.

The main conclusion from these considerations is that at the only nontrivial IR attractive fixed point correlation time vanishes. In other words, the Galilean symmetry, broken by the colored random force in Equation (145), is reestablished in the IR limit.

9. Parity Breaking in Magnetohydrodynamic Turbulence

In addition to the passive scalar problems introduced in Section 3, there exists a broad class of problems related to vector admixtures [23,35]. Due to the presence of more degrees of freedom, it is natural to expect that observed behavior can be richer than in scalar case. Among many, much attraction has been gained by a model of magnetohydrodynamics (MHD) [73,74,75,76]. The interplay between the velocity field and the magnetic field is crucial in explaining many phenomena—magnetic dynamo, convective processes, galaxy formation, etc. [77,78,79,80,81]—therefore it is clear that it has to be taken into account properly.

The full stochastic MHD problem is rather complicated, therefore for many purposes for a theoretical description it has been proposed to use a simplified model—the kinematic Kazantsev– Kraichnan model [35]. Its basic assumption is that the vector field of magnetic induction $\mathit{B}$ (hereinafter, referred to as vector) is passively advected by the turbulent flow, but reaction on the velocity field by the magnetic field is neglected. The most notable point of criticism to the Kazantsev–Kraichnan model is the assumption about the velocity field using a simple Gaussian statistical ensemble. From a more fundamental point of view, the velocity field would be generated dynamically. Therefore, we assume here that the velocity field is brought about by the stochastic Navier–Stokes equation expounded in Section 3.

The introduction of magnetic field in the Kazantsev–Kraichnan model comes from a physical approximation called magnetohydrodynamic limit [74,75]. It is assumed that the matter is electrically neutral at large scales and that the mean free path is much shorter than the corresponding Debye length. This assures that the electric displacement field accounting for the overall motion of charged particles can be entirely neglected. The system is then completely described in terms of the variables of density, pressure, and mean velocity field. From the RG point of view the displacement field is IR irrelevant. Hence, there is no need to retain it. Faraday’s law ${\partial }_t\mathit{B}=-\nabla \times \mathit{E}$ together with the generalized Ohm’s law $\mathit{J}=\sigma (\mathit{E}+\mathit{v}\times \mathit{B})$ give rise to the advection–diffusion equation ${\partial }_t\mathit{B}-\nabla \times (\mathit{v}\times \mathit{B})={\kappa }_0{\nabla }^2\mathit{B}$. According to the functional formulation from Section 2, the corresponding stochastic problem assumes the form

$${\partial }_t{b}_i+{\partial }_k(v_k{b}_i-v_i{b}_k)={\kappa }_0{\partial }^2{b}_i+f_i^{\theta },$$

(158)

where $b_i$ is just the fluctuating component of the magnetic induction, ${\kappa }_0$ is the magnetic diffusion coefficient, and a stochastic source term $f_i^{\theta }$ has been added to the right side. This term is the random part of the curl of the current and stems from the intrinsic randomness of the magnetic field [31]. A detailed account of the MHD problem can be found in [73,74,75].

The random source $f_i^{\theta }$ in Equation (158) is assumed to be zero-mean Gaussian with the correlation function

$$\langle f_i^b(t,\mathit{x})f_j^b(t^{\prime },{\mathit{x}}^{\prime })\rangle =\delta (t-t^{\prime })C_{ij}(\mathit{r}/L_{\theta }),\mathit{r}=\mathit{x}-{\mathit{x}}^{\prime },$$

(159)

where $C_{ij}(\mathbf{r}/L_{\theta })$ is a function, whose exact form is irrelevant. It must have a finite limit at $(\mathit{r}/L_{\theta })\to 0$ and vanish at $(\mathit{r}/L_{\theta })\to \infty$. The magnetic field $b_i$ is solenoidal, therefore the terms ${\partial }_k\left(v_i{b}_k\right)$ and $\left(b_k{\partial }_k\right)v_i$ are equal.

In realistic scenarios, the stochastic NS equation (Equation (16)) has to be amended by a Lorentz term responsible for the backward influence of the magnetic field on the velocity field. This is brought about by a familiar force term in conducting fluid of form $\mathit{v}\times \mathit{B}\sim \mathit{J}\sim (\nabla \times \mathit{B})\times \mathit{B}$. The most important consequence of this argument is an important physical effect known as the turbulent dynamo: generation of magnetic field at large scales by the turbulent motion. Let us make an important remark, which should clarify the fundamental difference between equilibrium statistical models and non-equilibrium ones. The turbulent dynamo might be explained by the mechanism of spontaneous symmetry breaking most spectacularly associated with the Higgs mechanism in particle physics. Here, “the ground state” of the turbulent gyrotropic fluid with vanishing mean $\mathit{b}=0$ of the magnetic field is unstable. It is stabilized by the spontaneous generation of a spatially uniform anomalous mean $\mathit{b}\ne 0$. This happens in full analogy with the situation in a ferromagnetic material below $T_c$, in which magnetic order is stabilized by the appearance of spontaneous magnetization. The condition for total magnetization follows from an extremal condition imposed on free energy. On the other hand, for the dynamo problem, the resulting field is determined by the properties of the response function. In particular, it is required that all perturbations from a given stable state have to be damped out sufficiently quickly. Let us note that that the instability manifests itself not at the level of the action function but only at the following one-loop approximation.

Let us discuss the full MHD problem in more detail and introduce one specific generalization. From now on, we assume that the MHD problem is governed by two coupled equations [22]

$$\begin{array}{cc}\hfill {\nabla }_t\mathit{v}& ={\nu }_0{\nabla }^2\mathit{v}+(\mathit{b}·\nabla )\mathit{b}-\nabla Q+{\mathit{f}}^v,\hfill \end{array}$$

(160)

$$\begin{array}{cc}\hfill {\nabla }_t\mathit{b}& ={\nu }_0{u}_0{\nabla }^2\mathit{b}+A(\mathit{b}·\nabla )\mathit{v}-\nabla P+{\mathit{f}}^b,\hfill \end{array}$$

(161)

where $\mathit{b}\equiv \mathit{b}\left(x\right)$ denotes a solenoidal vector field advected by a helical turbulent flow of an incompressible velocity field $\mathit{v}\equiv \mathit{v}\left(x\right)$. Both fields $\mathit{v}$ and $\mathit{b}$ are divergenceless. In other words, they satisfy the incompressibility condition

$$\nabla ·\mathit{b}=\nabla ·\mathit{v}=0.$$

(162)

Further, $u_0$ is the bare inverse Prandtl number [23,26]. Functions $P\equiv P\left(x\right)$ and $Q\equiv Q\left(x\right)$ in Equations (160) and (161) representing pressure fields are not relevant in the following analysis. In fact, due to the solenoidal property in Equation (162), functions P and Q can be expressed in terms of a formal Biot–Savart law [23,26]. A is a dimensionless parameter of the model. Three particular values of the parameter A have been analyzed in detail [35]. First, the value $A=1$ yields the Kazantsev–Kraichnan kinematic dynamo model, where the pressure term in Equation (161) drops out. Second, the choice $A=0$ leads to the model of passively advected vector admixture. Third, the value $A=-1$ corresponds to the linearized NS equation in a background field with given statistics. It is therefore convenient to retain A as a free parameter and analyze all three models together. The model in Equation (161) is often called the generic A model in the literature [35,82,83]. The complete model thus comprises two interconnected stochastic equations (Equations (160) and (161)).

Stochastic random sources $f^v$ and $f^b$ must be included in the coarse-grained setup of turbulence and MHD problems. As usual, it is assumed that the random variable $f^b$ is Gaussian with zero mean and correlation function

$$D_{ij}^b(x;0)\equiv \langle f_i^b\left(x\right)f_j^b\left(0\right)\rangle =\delta \left(t\right)C_{ij}\left(\right|\mathit{x}|/L).$$

(163)

Here, L is the integral scale of stirring of the magnetic field $\mathit{b}$, $C_{ij}$ is finite at $L\to \infty$ and at $\left|x\right|\gg L$ it should rapidly fall off. Apart from these restrictions, $C_{ij}$ needs not be specified [35]. Further, $f^v$ mimics the pumping of kinetic energy in the system at the largest spatial scales and is subject to the condition of real IR energy injection implemented through its specific power-like form suited to the RG approach [5]. However, the results obtained do not depend on specific forcing statistics because of the universality of fully developed turbulence. Hence, we only take advantage of the most suitable forcing in the RG analysis. Following this standard argument [5], we choose the pair correlation function of zero-mean Gaussian statistics

$$D_{ij}^v(x;0)\equiv \langle f_i^v\left(x\right)f_j^v\left(0\right)\rangle =\delta \left(t\right)\int \frac{{\mathrm{d}}^{d}k}{{\left(2\pi \right)}^d}{D}_0{k}^{4-d-2\epsilon }{R}_{ij}\left(\mathit{k}\right){\mathrm{e}}^{i\mathit{k}·\mathit{x}},$$

(164)

where $d=3$ denotes the space dimension; $\mathit{k}$ is the wavevector with the magnitude $k=\left|\mathit{k}\right|$; and $D_0\equiv g_0{\nu }_0^3>0$ is the positive quantity, where $g_0$ is a coupling constant connected to the typical UV momentum scale $\Lambda$ through $g_0\simeq {\Lambda }^{2\epsilon }$. The parameter $\epsilon$ specifies the power-like behavior of the energy pumping at large scales and assumes the value of 2 in the physically relevant IR energy pumping. In the RG analysis of the fully developed turbulence, $\epsilon$ is assumed small in calculations. Its physical value of 2 is inserted into perturbative expansions only as the last step [5]. The tensor quantity $R_{ij}$ determines the spatial parity violation in the model at hand. In symmetric isotropic incompressible turbulent fluid, such as analyzed here, the tensor $R_{ij}\left(k\right)$ corresponds to the sum of the ordinary transverse projection operator $P_{ij}\left(\mathit{k}\right)={\delta }_{ij}-k_i{k}_j/k^2$ and a helical term $H_{ij}\left(\mathit{k}\right)=i\rho {\epsilon }_{ijl}{k}_l/k$, i.e.,

$$R_{ij}\left(\mathit{k}\right)={\delta }_{ij}\left(\mathit{k}\right)-\frac{k_i{k}_i}{{\mathit{k}}^2}+i\rho {\epsilon }_{ijl}\frac{k_l}{\left|\mathit{k}\right|},$$

(165)

where ${\epsilon }_{ijl}$ is the third-rank completely antisymmetric tensor. The real-valued helicity parameter $\rho$ obeys the inequality $\left|\rho \right|\le 1$ dictated by the condition that the correlation function is positive definite in Equation (164). Evidently, the value $\rho =0$ corresponds to the isotropic (non-helical) case, while the value $\rho =1$ corresponds to fully broken spatial parity.

Following the De Dominicis–Janssen approach (Section 2), the stochastic problem of Equations (160) and (161) is equivalent to a field-theoretic model with the set of fields $\mathsf{\Phi }=\{\mathit{v},\mathit{b},{\mathit{v}}^{\prime },{\mathit{b}}^{\prime }\}$, where the fields with primes once more denote the response fields [5,45]. Thus, the field-theoretic model is given by the Dominicis–Janssen action

$$\mathcal{S}\left[\mathsf{\Phi }\right]=\frac{1}{2}\left[v_i^{\prime }{D}_{ij}^v{v}_j^{\prime }+b_i^{\prime }{D}_{ij}^b{b}_j^{\prime }\right]+{\mathit{v}}^{\prime }[-{\nabla }_t\mathit{v}+{\nu }_0{\nabla }^2\mathit{v}+(\mathit{b}·\nabla )\mathit{b}]+{\mathit{b}}^{\prime }[-{\nabla }_t\mathit{b}+{\nu }_0{u}_0{\nabla }^2\mathit{b}+A(\mathit{b}·\nabla )\mathit{v}],$$

(166)

where $D_{ij}^b$ and $D_{ij}^v$ are defined in Equations (163) and (164), respectively. As usual, summations over repeated indices $i,j\in 1,2,3$ are implied. As the original fields $\mathit{v}$ and $\mathit{b}$, the response fields are solenoidal, i.e., $\nabla ·{\mathit{v}}^{\prime }=\nabla ·\mathit{b}=0$. In the frequency–momentum representation, free propagators of the model in Equation (166) are

$$\begin{array}{cccc}\hfill {\langle b_i^{\prime }{b}_j\rangle }_0& ={\langle b_i{b}_j^{\prime }\rangle }_0^{*}=\frac{P_{ij}\left(\mathbf{k}\right)}{i\omega +{\nu }_0{u}_0{k}^2},\hfill & \hfill {\langle b_i{b}_j\rangle }_0& =\frac{C_{ij}\left(\mathbf{k}\right)}{|-i\omega +{\nu }_0{u}_0{k}^2{|}^2},\hfill \end{array}$$

(167)

$$\begin{array}{cccc}\hfill {\langle v_i^{\prime }{v}_j\rangle }_0& ={\langle v_i{v}_j^{\prime }\rangle }_0^{*}=\frac{P_{ij}\left(\mathbf{k}\right)}{i\omega +{\nu }_0{k}^2},\hfill & \hfill {\langle v_i{v}_j\rangle }_0& =\frac{g_0{\nu }_0^3{k}^{4-d-2\epsilon }{R}_{ij}\left(\mathbf{k}\right)}{|-i\omega +{\nu }_0{k}^2{|}^2}.\hfill \end{array}$$

(168)

The function $C_{ij}\left(k\right)$ is the Fourier transform of the function $C_{ij}(r/L)$ introduced in Equation (163). Moreover, the theory includes three interaction vertices

• ${\mathcal{S}}_{b^{\prime }bv}$: $b_i^{\prime }(-v_j{\partial }_j{b}_i+A{b}_j{\partial }_j{v}_i)=b_i^{\prime }{v}_j{V}_{ijl}{b}_l$,
• ${\mathcal{S}}_{v^{\prime }vv}$:$-v_i^{\prime }{v}_j{\partial }_j{v}_i=v_i^{\prime }{v}_j{W}_{ijl}{v}_l/2,$
• ${\mathcal{S}}_{v^{\prime }bb}$:$v_i^{\prime }{b}_j{\partial }_j{b}_i=v_i^{\prime }{v}_j{U}_{ijl}{v}_l/2,$

In the momentum–frequency representation, they are associated with the vertex factors [5]

$$V_{ijl}=i(k_j{\delta }_{il}-A{k}_l{\delta }_{ij}),W_{ijl}=i(k_l{\delta }_{ij}+k_j{\delta }_{il}),U_{ijl}=i(k_l{\delta }_{ij}+k_j{\delta }_{il}).$$

(169)

In all these vertices, momentum $\mathit{k}$ is flowing into through the corresponding response field, i.e., in $V_{ijl}$ it is the response field ${\mathit{b}}^{\prime }$ and in $W_{ijl},U_{ijl}$ the field ${\mathit{v}}^{\prime }$. A graphical representation of interaction vertices is displayed in Figure 4.

The field theoretic model defined by the action in Equation (166) is intrinsically unstable due to the fact that 1PI graphs ${\langle v^{\prime }v\rangle }_{1-\mathrm{IR}}$ and ${\langle b^{\prime }b\rangle }_{1-\mathrm{IR}}$ are not UV finite. To ensure the stabilization of the advection–diffusion system, we assume—inspired by the argument in [84]—that the vector field $\mathit{b}$ fluctuates around a spontaneously generated non-vanishing mean $\mathit{B}=\langle \mathit{b}\rangle \ne 0$ with the magnitude depending on the parameter A. The auxiliary field ${\mathit{b}}^{\prime }$ is still supposed to have vanishing mean. Technically, the spontaneous symmetry breaking is achieved by the substitution of the sum $\mathit{B}+\mathit{b}$ instead of the vector field $\mathit{b}$ in the action in Equation (166). Such a change of variables gives rise to a new action functional with the free part now containing additional terms, whereas the interacting part remains intact. The action functional with broken symmetry is

$$\begin{array}{cc}\hfill \mathcal{S}\left[\varphi \right]& =\frac{1}{2}\left[v_i^{\prime }{D}_{ij}^v{v}_j^{\prime }+b_i^{\prime }{D}_{ij}^b{b}_j^{\prime }\right]+{\mathit{v}}^{\prime }[-{\nabla }_t\mathit{v}+{\nu }_0{\nabla }^2\mathit{v}+(\mathit{b}·\nabla )\mathit{b}]+{\mathit{b}}^{\prime }[-{\nabla }_t\mathit{b}+{\nu }_0{u}_0{\nabla }^2\mathit{b}\hfill \\ \hfill & +A(\mathit{b}·\nabla )\mathit{v}]+{\mathit{v}}^{\prime }(\mathit{B}·\nabla )\mathit{b}+A{\mathit{b}}^{\prime }(\mathit{B}·\nabla )\mathit{v}.\hfill \end{array}$$

(170)

The quadratic part of this action determines propagators of the theory and now it has more involved structure than Equation (166). It is seen that the symmetry breaking brings about cross propagators $\langle v{b}^{\prime }\rangle$, $\langle b{v}^{\prime }\rangle$, $\langle bv\rangle$, and $\langle bb\rangle$. Furthermore, all propagators are more complicated and depend on the uniform magnetic field $B$ explicitly. In the frequency–momentum representation, they are

$$\begin{array}{cccc}\hfill \langle v_i{v}_j\rangle & =\frac{\beta \left(\mathit{k}\right){\beta }^{*}\left(\mathit{k}\right)}{\xi \left(\mathit{k}\right){\xi }^{*}\left(\mathit{k}\right)}{D}^v\left(\mathit{k}\right)R_{ij}\left(\mathit{k}\right),\hfill & \hfill \langle b_i{b}_j\rangle & =A^2\frac{{(\mathit{B}·\mathit{k})}^2}{\xi \left(\mathit{k}\right){\xi }^{*}\left(\mathit{k}\right)}{D}^v\left(\mathit{k}\right)R_{ij}\left(\mathit{k}\right)\hfill \end{array}$$

$$\begin{array}{cccc}\hfill \langle v_i{v}_j^{\prime }\rangle & =\frac{{\beta }^{*}\left(\mathit{k}\right)}{{\xi }^{*}\left(\mathit{k}\right)}{P}_{ij}\left(\mathit{k}\right),\hfill & \hfill \langle b_i{b}_j^{\prime }\rangle & =\frac{{\alpha }^{*}\left(\mathit{k}\right)}{{\xi }^{*}\left(\mathit{k}\right)}{P}_{ij}\left(\mathit{k}\right),\hfill \end{array}$$

(171)

$$\begin{array}{cccc}\hfill \langle b_i{v}_j^{\prime }\rangle & =iA\frac{(\mathit{B}·\mathit{k})}{{\xi }^{*}\left(\mathit{k}\right)}{P}_{ij}\left(\mathit{k}\right),\hfill & \hfill \langle v_i{b}_j^{\prime }\rangle & =i\frac{(\mathit{B}·\mathit{k})}{{\xi }^{*}\left(\mathit{k}\right)}{P}_{ij}\left(\mathit{k}\right),\hfill \end{array}$$

(172)

$$\begin{array}{cc}\hfill \langle b_i{v}_j\rangle & =iA\frac{\beta \left(\mathit{k}\right)(\mathit{B}.\mathbf{k})}{\xi \left(\mathbf{k}\right){\xi }^{*}\left(\mathbf{k}\right)}{D}_v\left(\mathbf{k}\right)R_{ij}\left(\mathbf{k}\right).\hfill \end{array}$$

(173)

Here, following abbreviations have been introduced

$$\alpha \left(\mathit{k}\right)=i\omega +\nu k^2,\beta \left(\mathit{k}\right)=i\omega +u\nu k^2,\xi \left(\mathit{k}\right)=A{(\mathit{B}·\mathit{k})}^2+\alpha \left(\mathit{k}\right)\beta \left(\mathit{k}\right)$$

(174)

for brevity. Propagators are depicted in Figure 5.

Identification of all relevant UV divergences is carried out by the standard power counting. We omit the details, which are analogous to those in Ref. [82]. The model at hand is logarithmic at $\epsilon =0$. In the framework of the minimal subtraction (MS) scheme, this means that all possible UV divergences are poles in $\epsilon$ [5]. Following the analysis of Hnatič and Zalom [82], we arrive at the conclusion that all UV divergences are absorbed in counterterms ${\mathit{v}}^{\prime }{\nabla }^2\mathit{v}$ and ${\mathit{b}}^{\prime }{\nabla }^2\mathit{b}$. Thus, the model is multiplicatively renormalizable with renormalized parameters $g_0$, $u_0$, and ${\nu }_0$

$${\nu }_0=\nu Z_{\nu },g_0=g{\mu }^{2\epsilon }{Z}_g,u_0=u{Z}_u,$$

(175)

where $g,u$, and $\nu$ are the corresponding renormalized parameters. The renormalization mass $\mu$ is introduced as a consequence of the analytic regularization used in calculations. The quantities $Z_i=Z_i(g,u;d,\rho ;\epsilon )$ contain poles in the exponent $\epsilon$. The renormalized action functional can be written as follows

$$\begin{array}{cc}\hfill {\mathcal{S}}_R\left[\varphi \right]& =\frac{1}{2}\left[v_i^{\prime }{D}_{ij}^v{v}_j^{\prime }+b_i^{\prime }{D}_{ij}^b{b}_j^{\prime }\right]+{\mathit{v}}^{\prime }[-{\nabla }_t\mathit{v}+\nu Z_1{\nabla }^2\mathit{v}+Z_3(\mathit{b}·\nabla )\mathit{b}]+{\mathit{b}}^{\prime }[-{\nabla }_t\mathit{b}+\nu u{Z}_2{\nabla }^2\mathit{b}\hfill \\ \hfill & +A(\mathit{b}·\nabla )\mathit{v}]+Z_3{\mathit{v}}^{\prime }(\mathit{B}·\nabla )\mathit{b}+A{\mathit{b}}^{\prime }(\mathit{B}·\nabla )\mathit{v},\hfill \end{array}$$

(176)

where $Z_1$ and $Z_2$ are renormalization constants defined by relations

$$Z_{\nu }=Z_1,Z_g=Z_1^{-3},Z_u=Z_2{Z}_1^{-1}.$$

(177)

Both renormalization constants, $Z_1$ and $Z_2$, correspond to a different class of Feynman diagrams (details below). However, they reveal a common structure in the MS scheme: the nth order of perturbation theory is related to the nth power of g with the corresponding coefficient containing poles in $\epsilon$ of order n and less [5].

The relevant one-loop-order Feynman diagrams are shown in Figure 6 and Figure 7.

One-loop calculation [84] of self-energy graphs ${\Sigma }_{ij}$ for the 1PI function ${\mathsf{\Gamma }}_{b_i^{\prime }{b}_j}$ yields the following expression

$${\Sigma }_{ij}\sim {\epsilon }_{isl}{p}_s{T}_{lj}$$

(178)

where

$$T_{lj}=a\Lambda {\delta }_{lj}-b\left|\mathit{B}\right|({\delta }_{lj}+e_l{e}_j),\mathit{e}\equiv \frac{\mathit{B}}{\left|\mathit{B}\right|}.$$

(179)

Here, $e_i$ are components of the unit vector pointing in the direction of the spontaneous magnetic field. Terms proportional to ${\delta }_{ij}$ might in principle generate instabilities. In the literature, the term proportional to $e_l{e}_j$ is known as exotic. As detailed analysis [84] shows, ${\delta }_{ij}$ terms can be eliminated by imposing the following relation

$$\left|\mathit{B}\right|=\frac{a\Lambda }{b}=\sqrt{\frac{1}{\pi \left|A\right|}}\frac{\mathsf{\Gamma }(d/2+3/2)}{\mathsf{\Gamma }(d/2+1)}{u}_{*}\nu \Lambda ,$$

(180)

where the parameter A should not be equal to 0 or $-1$. At real space dimension $d=3$, we get

$$\left|\mathit{B}\right|=\frac{8}{3\pi \sqrt{\left|A\right|}}{u}_{*}\nu \Lambda .$$

(181)

To get insight into relevance of the last term on the right hand side of Equation (179), let us consider the approximate case of linearized MHD equations in polarized medium

$${\partial }_t\mathit{v}=\nu k^2\mathit{v}+i\gamma \mathit{b},{\partial }_t\mathit{b}=u\nu k^2\mathit{v}+i\gamma \mathit{v}-i\mu [\mathit{k}\times \mathit{e}](\mathit{e}·\mathit{b}),$$

(182)

where for brevity we have denoted $\gamma =\mathit{B}·\mathit{k}$, introduced $\mu$ as notation for certain combination of model parameters and employed the time–momentum representation [84]. Solution is sought in the form of plane waves

$$\mathit{v}\equiv \mathit{v}(t,\mathit{x})=\mathit{v}\left(t\right){\mathrm{e}}^{i\mathit{k}·\mathit{x}},\mathit{b}\equiv \mathit{b}(t,\mathit{x})=\mathit{b}\left(t\right){\mathrm{e}}^{i\mathit{k}·\mathit{x}}.$$

(183)

For an inviscid medium ($\nu =0$) without the exotic term ($\mu =0$), we get a solution known as Alfvén waves for which

$$\mathit{v}\left(t\right)\sim \mathit{b}\left(t\right)\sim {\mathrm{e}}^{-i\omega t},\omega \equiv \gamma .$$

(184)

The solution with the exotic term is more complicated. To find solutions with the exotic term ($\mu \ne 0$), let us first find a convenient orthonormal basis of three vectors ${\mathit{e}}_1,{\mathit{e}}_2,{\mathit{e}}_3$. The following choice

$${\mathit{e}}_1\equiv \frac{\mathit{k}}{\left|\mathit{k}\right|},{\mathit{e}}_2\equiv \frac{\mathit{e}-{\mathit{e}}_{1}cos\theta }{sin\theta },{\mathit{e}}_3\equiv [{\mathit{e}}_1\times {\mathit{e}}_2]=\frac{{\mathit{e}}_1\times \mathit{e}}{sin\theta }$$

(185)

turns out to be case and qualitatively is depicted in Figure 8.

The transverse velocity and magnetic vector fields are perpendicular to the wavevector $k$ or, equivalently, to the basis vector $e_1$. Let us decompose fields into perpendicular components $e_2$ and $e_3$ as follows

$$\mathit{v}\left(t\right)=v_2\left(t\right){\mathit{e}}_2+v_3\left(t\right){\mathit{e}}_3,\mathit{b}\left(t\right)=b_2\left(t\right){\mathit{e}}_2+b_3\left(t\right){\mathit{e}}_3.$$

(186)

It is a straightforward exercise to check that modes $v_i,b_i;i=2,3$ satisfy following equations

$$\begin{array}{cccc}\hfill {\partial }_t{v}_2\left(t\right)& =i\gamma b_2\left(t\right),\hfill & \hfill {\partial }_t{b}_2\left(t\right)& =i\gamma v_2\left(t\right),\hfill \end{array}$$

(187)

$$\begin{array}{cccc}\hfill {\partial }_t{v}_3\left(t\right)& =i\gamma b_3\left(t\right),\hfill & \hfill {\partial }_t{b}_3\left(t\right)& =i\gamma v_3\left(t\right)+2i\lambda b_2\left(t\right).\hfill \end{array}$$

(188)

Solution of these equations can be represented in the form

$$\begin{array}{c}\hfill c{b}_2\left(t\right)=-v_2\left(t\right)=b_2{\mathrm{e}}^{-i\gamma t},\\ \hfill b_3\left(t\right)=[b_3+i\lambda b_{2}t]{\mathrm{e}}^{-i\gamma t},\\ \hfill v_3\left(t\right)=[-b_3-\lambda /\gamma b_2-i\lambda b_{2}t]{\mathrm{e}}^{-i\gamma t}.\end{array}$$

(189)

Thus, field-theoretic methods not only explain the emergence of turbulent dynamo, but also lead to an accompanying physical effect: appearance of disturbances in Alfvén waves perpendicular to the spontaneous field $\mathit{B}$, which leads to their linear growth in time. This generates long-lived pulses $texp(-u\nu k^{2}t)$—similar to Goldstone bosons.

An intriguing unsolved problem remains, however. Although gyrotropy is the cause of the dynamo effect, the value in Equation (180) calculated for the spontaneous magnetic field is independent of the parameter $\rho$. Thus, any however small gyrotropy causes emergence of a finite spontaneous field. It is thus not clear whether or not the dynamo will also occur in a normal fluid in which $\rho =0$.

10. Effect of Strong Anisotropy

Another related problem that can be addressed is the effect of large-scale anisotropy on statistics of the velocity field, and the passively advected fields in the inertial range [34,36,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99]. The classical Kolmogorov–Obukhov theory predicts that the anisotropy generated at large spatial scales by the forcing (boundary conditions, geometry of the mesh, etc.) fades away when the energy is transferred down to the smaller scales by means of the cascade mechanism [23,25]. This picture is corroborated in recent analyses for even correlation functions. Thus, they provide us with some backup to the hypothesis of the restored local isotropy of the turbulence for the velocity field and passively advected field in the inertial range [34,36,90,91,92,93,94,95,97,98,99]. More specifically, the exponents describing the scaling in the inertial range exhibit universality and hierarchy related to the state of anisotropy. In particular, the main contribution to an even function comes from the exponent from the isotropic shell [34,36,89,91,92,93,94,97,98,99]. However, the anisotropy survives in the inertial range being explicit in odd correlation functions. This is inconsistent with the behavior anticipated on grounds of the cascade picture. Further, the skewness factor decreases with the length scale much more slowly than expected [50,85,86,87,88,100,101,102,103], whereas the odd dimensionless ratios of structure functions of higher order (hyperskewness, etc.) increase. This hints at persistent anisotropy at small scales [34,36,90,92]. This appears a universal effect as it is observed for both the scalar [34,36] and vector [92] fields, advected by the Gaussian rapid-change velocity, as well as for the scalar advected by the two-dimensional velocity field generated by the stochastic Navier–Stokes equation [90].

Here, we demonstrate the anomalous scaling behavior of a passive scalar advected by the velocity field due to the Kraichnan model with strong anisotropy. In contrast with the studies in [34,36,50,85,86,87,88,103], where the velocity field was assumed isotropic and the anisotropy was generated at large scales by the imposed linear mean gradient, the uniaxial anisotropy in considered model is present at all scales. The anomalous exponents turn out to depend on the anisotropy parameters and are thus not universal.

The aim of this section is twofold. First, expressions for the structure functions and correlation functions of the scalar gradients are obtained and then the corresponding anomalous exponents are computed in the leading order of the $\epsilon$ expansion. Due to the anisotropy of the velocity flow, the composite operators of different ranks mix strongly with each other under renormalization. As a direct consequence, the corresponding anomalous exponents are given by the eigenvalues of matrices of generic structure (this is in contrast with the case of large-scale anisotropy, in which the matrices are diagonal or triangular). In terms of the zero-mode approach [93,94], this basically means that the $SO\left(d\right)$ decompositions of the correlation functions do not diagonalize differential operators in the corresponding exact equations.

As mentioned in Section 3, in the rapid-change model, the passive advection of a scalar field $\theta \left(x\right)\equiv \theta (t,x)$ is described by the stochastic equation (Equation (33)). The velocity $v\left(x\right)$ correlator now instead of Equation (40) is assumed in the form

$$\langle v_i\left(x\right)v_j\left(x^{\prime }\right)\rangle =D_0\frac{\delta (t-t^{\prime })}{{\left(2\pi \right)}^d}\int {\mathrm{d}}^{d}k{T}_{ij}\left(\mathit{k}\right){(k^2+m^2)}^{-d/2-\epsilon /2}exp[\mathrm{i}\mathit{k}·(\mathit{x}-{\mathit{x}}^{\prime })],$$

(190)

where

$$\frac{D_0}{{\nu }_0}\equiv g_0\equiv {\Lambda }^{\epsilon }$$

(191)

and the relation $m=1/L$ holds. In the isotropic case, the tensor quantity $T_{ij}\left(\mathit{k}\right)$ in Equation (190) is the usual transverse projector $T_{ij}\left(\mathit{k}\right)=P_{ij}\left(\mathit{k}\right)$ (see Equation (18)). The velocity statistics is now assumed to be anisotropic in entire region of scales. The ordinary transverse projector is then replaced by the general transverse structure

$$T_{ij}\left(\mathit{k}\right)=a\left(\psi \right)P_{ij}\left(\mathit{k}\right)+b\left(\psi \right){\tilde{n}}_i\left(\mathit{k}\right){\tilde{n}}_j\left(\mathit{k}\right),$$

(192)

where the unit vector $\mathit{n}$ denotes the singled out direction (${\mathit{n}}^2=1$),

$${\tilde{n}}_i\left(\mathit{k}\right)\equiv P_{ij}\left(\mathit{k}\right)n_j,$$

(193)

and $\psi$ is the angle between the vectors $\mathit{k}$ and $\mathit{n}$. In other words, $(\mathit{n}·\mathit{k})=kcos\psi$ [note that $(\tilde{\mathit{n}}·\mathit{k})=0$].

The positivity of the correlator in Equation (190) leads to the inequalities

$$a\left(\psi \right)>0,a\left(\psi \right)+b\left(\psi \right){sin}^2\psi >0.$$

(194)

In practical calculations, one works with the special case

$$T_{ij}\left(\mathit{k}\right)=[1+{\alpha }_1\frac{{(\mathit{n}·\mathit{k})}^2}{k^2}]P_{ij}\left(\mathit{k}\right)+{\alpha }_2{\tilde{n}}_i\left(\mathit{k}\right){\tilde{n}}_j\left(\mathit{k}\right).$$

(195)

Then, the inequalities in Equation (194) reduce to ${\alpha }_{1,2}>-1$.

Let us note that the quantities in Equations (192) and (195) are invariant with respect to transformation $\mathit{n}\to -\mathit{n}$. The anisotropy permits an introduction of the mixed correlator $\langle \mathit{v}f\rangle \propto \mathit{n}\delta (t-t^{\prime })C^{\prime }(r/L)$ with a function $C^{\prime }(r/L)$ similar to $C(r/L)$ in Equation (34). This violates the evenness in $n$ and gives rise to non-vanishing odd functions $S_{2n+1}$. However, this does not alter RG analysis [104].

The stochastic problem to be analyzed is tantamount to the field theoretic model of the set of fields $\varphi \equiv \{{\theta }^{\prime },\theta ,\mathit{v}\}$ with action

$$\mathcal{S}\left[\varphi \right]=\frac{1}{2}{\theta }^{\prime }{D}_{\theta }{\theta }^{\prime }+{\theta }^{\prime }\left[-{\partial }_t-(\mathit{v}·\nabla )+{\nu }_0{\nabla }^2+\frac{1}{2}{D}_{vij}\left(0\right){\partial }_i{\partial }_j\right]\theta -\frac{1}{2}\mathit{v}{D}_v^{-1}\mathit{v}.$$

(196)

Here, $D_{\theta }$ and $D_v$ are the correlators in Equations (34) and (190), respectively, and

$$D_{vij}\left(0\right)=D_0\int \frac{{\mathrm{d}}^{d}q}{{\left(2\pi \right)}^d}\frac{T_{ij}\left(\mathit{q}\right)}{{(q^2+m^2)}^{d/2+\epsilon /2}}$$

(197)

is the diagonal term (in spatial variables) of the coefficient of the temporal $\delta$ function in the velocity pair correlation function in Equation (190)).

The model in Equatio (196) corresponds to the usual Feynman diagrammatic technique with the triple vertex in Equation (51), propagators in Equation (50) and

$$\begin{array}{c}\hfill {\langle \theta \theta \rangle }_0=C\left(k\right){({\omega }^2+{\nu }_0^2{k}^4)}^{-1},{\langle {\theta }^{\prime }{\theta }^{\prime }\rangle }_0=0,\end{array}$$

(198)

where $C\left(k\right)$ is the Fourier transform of the function $C(r/L)$ in Equation (34). The bare propagator ${\langle vv\rangle }_0\equiv \langle vv\rangle$ is defined by Equation (190) with the transverse projection operator from Equation (192) or Equation (195).

The pair correlation functions $\langle \varphi \varphi \rangle$ of the multicomponent field $\varphi \equiv \{{\theta }^{\prime },\theta ,v\}$ fulfill the Dyson equation. In terms of component fields, we arrive at the system of two equations (cf. [25])

$$\begin{array}{cc}\hfill G^{-1}(\omega ,\mathit{k})& =-\mathrm{i}\omega +{\nu }_0{k}^2-{\Sigma }_{{\theta }^{\prime }\theta }(\omega ,\mathit{k}),\hfill \end{array}$$

(199)

$$\begin{array}{cc}\hfill D(\omega ,\mathit{k})& ={\left|G(\omega ,\mathit{k})\right|}^2[C\left(k\right)+{\Sigma }_{{\theta }^{\prime }{\theta }^{\prime }}(\omega ,\mathit{k})],\hfill \end{array}$$

(200)

where $G(\omega ,\mathit{k})\equiv \langle \theta {\theta }^{\prime }\rangle$ and $D(\omega ,\mathit{k})\equiv \langle \theta \theta \rangle$ are the exact response function and pair correlation function, respectively, and ${\Sigma }_{{\theta }^{\prime }\theta }$ and ${\Sigma }_{{\theta }^{\prime }{\theta }^{\prime }}$ are self-energy operators represented by 1PI graphs. The rest of the self-energy functions ${\Sigma }_{\varphi \varphi }$ in the model in Equation (196) vanish identically.

It is a characteristic feature of models such as that in Equation (196) that all the skeleton multi-loop diagrams of self-energy functions ${\Sigma }_{{\theta }^{\prime }\theta }$, ${\Sigma }_{{\theta }^{\prime }{\theta }^{\prime }}$ contain closed circuits of retarded propagators $\langle \theta {\theta }^{\prime }\rangle$ (it is important that the propagator ${\langle \mathit{v}\mathit{v}\rangle }_0$ in Equation (190) is $\delta$ correlated in time) and hence vanish.

In the presence of anisotropy, a new counterterm of the form ${\theta }^{\prime }{(\mathit{n}·\nabla )}^2\theta$ has to be introduced. There is no such term in the unrenormalized action functional in Equation (196). Thus, the model in Equation (196) in its original form is not multiplicatively renormalizable, and to employ RG techniques we have to extend the model by adding a contribution corresponding to the counterterm to the bare action

$$\mathcal{S}\left[\mathsf{\Phi }\right]={\theta }^{\prime }{D}_{\theta }{\theta }^{\prime }/2+{\theta }^{\prime }\left[-{\partial }_t-(\mathit{v}·\nabla )+{\nu }_0{\nabla }^2+{\chi }_0{\nu }_0{(\mathit{n}·\nabla )}^2\right]\theta -\mathit{v}{D}_v^{-1}\mathit{v}/2.$$

(201)

Here, ${\chi }_0$ is a new dimensionless parameter. For stability of the system, positivity of the viscous contribution ${\nu }_0{k}^2+{\chi }_0{\nu }_0{(\mathit{n}·\mathit{k})}^2$ is required, which leads to the inequality ${\chi }_0>-1$. The real physical value of the new parameter ${\chi }_0$ is zero. However, this does not hinder the application of the RG techniques. It is first assumed to be arbitrary, and the equality ${\chi }_0=0$ is imposed as the initial condition in the solution of equations for invariant variables. The vanishing ${\chi }_0$ corresponds to some nonzero value of its renormalized counterpart, which can be explicitly calculated.

For the action in Equation (201), the nontrivial bare propagators in Equation (198) are replaced by

$$\begin{array}{cc}\hfill {\langle \theta {\theta }^{\prime }\rangle }_0={\langle {\theta }^{\prime }\theta \rangle }_0^{*}& =\frac{1}{-\mathrm{i}\omega +{\nu }_0{k}^2+{\chi }_0{\nu }_0{(\mathit{n}·\mathit{k})}^2},\hfill \end{array}$$

(202)

$$\begin{array}{cc}\hfill {\langle \theta \theta \rangle }_0& =\frac{C\left(k\right)}{|-\mathrm{i}\omega +{\nu }_0{k}^2+{\chi }_0{\nu }_0{(\mathit{n}·\mathit{k})}^2{|}^2}.\hfill \end{array}$$

(203)

Once properly extended, the model becomes multiplicatively renormalizable: due to generation of counterterms, two independent renormalization constants $Z_{1,2}$ are introduced as coefficients of the counterterms. This yields the renormalized action in the form

$${\mathcal{S}}_R\left[\mathsf{\Phi }\right]={\theta }^{\prime }{D}_{\theta }{\theta }^{\prime }/2+{\theta }^{\prime }\left[-{\partial }_t-(\mathit{v}·\nabla )+\nu Z_1{\nabla }^2+\chi \nu Z_2{(\mathit{n}·\nabla )}^2\right]\theta -\mathit{v}{D}_v^{-1}\mathit{v}/2,$$

(204)

or to the multiplicative renormalization of all the parameters ${\nu }_0$, $g_0$ and ${\chi }_0$ of the action in Equation (201):

$${\nu }_0=\nu Z_{\nu },g_0=g{\mu }^{\epsilon }{Z}_g,{\chi }_0=\chi Z_{\chi }.$$

(205)

The correlator (Equation (190)) in Equation (204) is expressed in terms of renormalized variables using Equations (205). From direct comparison of Equations (201), (204) and (205), we obtain the relations

$$Z_1=Z_{\nu },Z_2=Z_{\chi }{Z}_{\nu },Z_g=Z_{\nu }^{-1}.$$

(206)

The beta functions are given by

$${\beta }_g(g,\alpha )\equiv {\tilde{\mathcal{D}}}_{\mu }g=g(-\epsilon -{\gamma }_g)=g(-\epsilon +{\gamma }_{\nu })=g(-\epsilon +{\gamma }_1),$$

(207)

$${\beta }_{\chi }(g,\chi )\equiv {\tilde{\mathcal{D}}}_{\mu }\chi =-\chi {\gamma }_{\chi }=\chi ({\gamma }_1-{\gamma }_2).$$

(208)

The relation between ${\beta }_g$ and ${\gamma }_{\nu }$ in Equation (207) is a consequence of the definitions and the last equality in Equation (206).

One-loop calculation [104] leads to the following expressions for the anomalous dimension ${\gamma }_1\left(g\right)$

$${\gamma }_1\left(g\right)=\frac{g{\overline{S}}_d}{2d(d+2)}\left[(d-1)(d+2)+{\alpha }_1(d+1)+{\alpha }_2\right],$$

(209)

and for ${\gamma }_2(g,\alpha )$

$${\gamma }_2(g,\alpha )=\frac{g{\overline{S}}_d}{2d(d+2)\chi }\left[-2{\alpha }_1+{\alpha }_2(d^2-2)\right],$$

(210)

respectively. Let us note that Equations (209) and (210) are exact [104].

From explicit Equations (209) and (210), we see that the RG equations have just one IR stable fixed point

$$g^{*}{\overline{S}}_d=\frac{2d(d+2)\epsilon }{(d-1)(d+2)+{\alpha }_1(d+1)+{\alpha }_2},{\chi }^{*}=\frac{-2{\alpha }_1+{\alpha }_2(d^2-2)}{(d-1)(d+2)+{\alpha }_1(d+1)+{\alpha }_2}.$$

(211)

At this point, both eigenvalues of the matrix $\Omega$ are equal to $\epsilon$; the values ${\gamma }_1^{*}={\gamma }_2^{*}={\gamma }_{\nu }^{*}=\epsilon$ are exact as a consequence of Equations (207) and (208) [here and below, ${\gamma }_F^{*}\equiv {\gamma }_F(g^{*},{\chi }^{*})$]. The fixed point in Equation (211) is degenerate. By this we mean that its coordinates depend explicitly on the anisotropy parameters ${\alpha }_{1,2}$.

Let us consider the solution of the RG equation in the example of the even two-time structure functions

$$S_{2n}(\mathit{r},\tau )\equiv \langle {[\theta (t,\mathit{x})-\theta (t^{\prime },{\mathit{x}}^{\prime })]}^{2n}\rangle ,\mathit{r}\equiv \mathit{x}-{\mathit{x}}^{\prime },\tau \equiv t-t^{\prime }.$$

(212)

They obey the RG equation ${\mathcal{D}}_{RG}{S}_{2n}=0$ with the differential operator ${\mathcal{D}}_{RG}={\mathcal{D}}_{\mu }+{\beta }_g{\partial }_g+{\beta }_{\chi }{\partial }_{\chi }-{\gamma }_{nu}{\mathcal{D}}_{\nu }$.

In renormalized variables, a dimensional argument leads to

$$S_{2n}(\mathit{r},\tau )={\nu }^{-n}{r}^{2n}{\tilde{R}}_{2n}(\mu r,\tau \nu /r^2,r/L,g,\chi ).$$

(213)

Here, ${\tilde{R}}_{2n}$ is a scaling function of dimensionless variables (dependence on d, $\epsilon$, ${\alpha }_{1,2}$ and the angle between the vectors $\mathit{r}$ and $\mathit{n}$ is implied). From the RG equation, the representation

$$S_{2n}(\mathit{r},\tau )={\left(\overline{\nu }\right)}^{-n}{r}^{2n}{\tilde{R}}_{2n}(1,\tau \overline{\nu }/r^2,r/L,\overline{g},\overline{\chi }),$$

(214)

is easily derived in terms of the invariant variables $\overline{e}=\overline{e}(\mu r,e)$. The identity $\overline{L}\equiv L$ is a consequence of the fact that L is not renormalized. The relation between the bare and invariant variables takes the form

$${\nu }_0=\overline{\nu }{Z}_{\nu }\left(\overline{g}\right),g_0=\overline{g}{r}^{-\epsilon }{Z}_g\left(\overline{g}\right),{\chi }_0=\overline{\chi }{Z}_{\chi }(\overline{g},\overline{\chi }).$$

(215)

Equation (215) defines the invariant variables as functions of the bare parameters in an implicit form. It is valid because expressions on both sides satisfy the same RG equation, and at $\mu r=1$ Equation (215) coincides with Equation (205).

The asymptotic behavior at large $\mu r$ of the invariant variables is governed by the IR stable fixed point: $\overline{g}\to g^{*}$, $\overline{\chi }\to {\chi }^{*}$ for $\mu r\to \infty$. However, in multi-charge problems, the possibility must be considered that, even when the stable IR point exists, not every phase trajectory approaches it in the limit $\mu r\to \infty$. The phase trajectory may first pass outside the natural region of stability (physical region is given by the inequalities $g>0$, $\chi >-1$) or go to infinity within this region. It can be easily verified that the RG flow reaches the fixed point in Equation (211) for arbitrary initial conditions $g_0>0$, ${\chi }_0>-1$, including the physical case ${\chi }_0=0$. Moreover, the large $\mu r$ behavior of the invariant viscosity $\overline{\nu }$ can be extracted explicitly from Equation (215) together with the last relation in Equation (206): $\overline{\nu }=D_0{r}^{\epsilon }/\overline{g}\to D_0{r}^{\epsilon }/g^{*}$ (we remind the reader that $D_0=g_0{\nu }_0$). Then, for $\mu r\to \infty$ and any fixed $mr$, we get

$$S_{2n}(\mathit{r},\tau )=D_0^{-n}{r}^{n(2-\epsilon )}{g^{*}}^n{R}_{2n}(\tau D_0{r}^{{\Delta }_t},r/L),$$

(216)

where the scaling function

$$R_{2n}(D_0\tau r^{{\Delta }_t},r/L)\equiv {\tilde{R}}_{2n}(1,D_0\tau r^{{\Delta }_t},r/L,g_{*},{\alpha }_{*}),$$

(217)

has appeared and ${\Delta }_t\equiv -2+{\gamma }_{\nu }^{*}=-2+\epsilon$ is the critical dimension of time. For the equal-time structure function, we have

$$S_{2n}\left(\mathit{r}\right)=D_0^{-n}{r}^{n(2-\epsilon )}{g^{*}}^n{R}_{2n}(r/L).$$

(218)

Here, the definition of $R_{2n}$ is evident from Equation (217). It should be noted that Equations (216)–(218) comprise a proof of the independence of the structure functions in the IR range (large $\mu r$ and any $r/L$) of the viscosity coefficient or, in other words, of the UV scale: the parameters $g_0$ and ${\nu }_0$ appear in Equation (216) only in the combination $D_0=g_0{\nu }_0$. A similar property has been found for the stochastic Navier–Stokes equation [105].

In contrast to the previously mentioned models, the scaling function $\tilde{R}$ in Equation (214) contains two different scales—corresponding to spatial and time differences, respectively. The information about its behavior can again be obtained using the OPE method from Section 5. Therefore, let $F(r,\tau )$ stands for some multiplicatively renormalized quantity. Dimensionality considerations imply

$$F_R(r,\tau )={\nu }^{d_F^{\omega }}{r}^{-d_F}{\tilde{R}}_F(\mu r,\tau \nu /r^2,r/L,g,\alpha ),$$

(219)

where $d_F^{\omega }$ and $d_F$ are the canonical dimension in frequency and the total canonical dimension of F, respectively (see Section 4 or [5]), and $R_F$ is a function of completely dimensionless variables. The analog of Equation (214) has the form

$$F(r,\tau )=Z_F(g,\alpha )F_R=Z_F(\overline{g},\overline{\alpha }){\left(\overline{\nu }\right)}^{d_F^{\omega }}{r}^{-d_F}{\tilde{R}}_F(1,\tau \overline{\nu }/r^2,r/L,\overline{g},\overline{\alpha }).$$

(220)

In the large $\mu r$ limit, one has $Z_F(\overline{g},\overline{\alpha })\simeq \mathrm{const}{(\Lambda r)}^{-{\gamma }_F^{*}}$ [106]. The UV scale comes into this relation from Equation (191). In the IR range ($\Lambda r\sim \mu r$ large, $r/L$ arbitrary), Equation (220) takes on the form

$$F(r,\tau )\simeq \mathrm{const}{\Lambda }^{-{\gamma }_F^{*}}{D}_0^{d_F^{\omega }}{r}^{-\Delta \left[F\right]}{R}_F(D_0\tau r^{{\Delta }_t},r/L).$$

(221)

Here,

$$\Delta \left[F\right]\equiv {\Delta }_F=d_F^k-{\Delta }_t{d}_F^{\omega }+{\gamma }_F^{*},{\Delta }_t=-2+\epsilon$$

(222)

is the critical dimension of the function F. The scaling function $R_F$ is related to ${\tilde{R}}_F$ as in Equation (217). For nonvanishing ${\gamma }_F^{*}$, the function F in the IR range exhibits dependence on $\Lambda$ or, equivalently, on ${\nu }_0$.

As a detailed analysis reveals [104], the operator ${\theta }^N$ requires no counterterms at all, i.e., it is in fact UV finite, ${\theta }^N=Z{\left[{\theta }^N\right]}^R$ with $Z=1$. As a straightforward consequence, the critical dimension of ${\theta }^N\left(x\right)$ is given by Equation (222) with no correction from ${\gamma }_F^{*}$ and reduces to the sum of the critical dimensions of the factors:

$$\Delta \left[{\theta }^N\right]=N\Delta \left[\theta \right]=N(-1+\epsilon /2).$$

(223)

The structure functions in Equation (212) are linear combinations of pair correlators involving the operators ${\theta }^N$, therefore Equation (223) shows that they indeed satisfy the RG equation of the form in Equation (65), as discussed in Section 4.

Tensor composite operators ${\partial }_{i_1}\theta \cdots {\partial }_{i_p}\theta {\left({\partial }_i\theta {\partial }_i\theta \right)}^n$ consisting solely of scalar gradients play an an important role. It is technically convenient to work with the scalar operators obtained by contracting the tensors with the vectors $\mathit{n}$ in the fashion

$$F[N,p]\equiv {\left[(\mathit{n}·\nabla )\theta \right]}^p{\left({\partial }_i\theta {\partial }_i\theta \right)}^n,N\equiv 2n+p.$$

(224)

Canonical dimensions of these operators depend on the total number of the fields $\theta$ and have the following form $d_F=0$, $d_F^{\omega }=-N$. These operators (Equation (224)) mix with each other only in the renormalization procedure. As a consequence, the corresponding infinite renormalization matrix

$$F[N,p]=\sum _{N^{\prime },p^{\prime }}{Z}_{[N,p][N^{\prime },p^{\prime }]}{F}^R[N^{\prime },p^{\prime }]$$

(225)

is block-triangular, i.e., $Z_{[N,p][N^{\prime },p^{\prime }]}=0$ for $N^{\prime }>N$. It is then clear that the critical dimensions associated with the operators $F[N,p]$ are fully determined by the eigenvalues of the finite subblocks with $N^{\prime }=N$.

In the isotropic case, as well as in the presence of large-scale anisotropy, the elements $Z_{[N,p][N,p^{\prime }]}$ vanish for $p<p^{\prime }$, and the block $Z_{[N,p][N,p^{\prime }]}$ is triangular together with the related blocks of the matrices $U_F$ and ${\Delta }_F$ in Equations (77) and (222). Moreover, it can be diagonalized in terms of to irreducible operators (scalars, vectors, and traceless tensors), but even for nonvanishing imposed gradient its eigenvalues coincide with those of the isotropic case. Hence, the introduction of large-scale anisotropy has no effect on critical dimensions of operators in Equation (224) (see [34,36]). In the case of small-scale anisotropy, the operators with different values of p mix in renormalization resulting in the matrix $Z_{[N,p][N,p^{\prime }]}$ is of generic form (neither diagonal nor triangula).

The calculation of the renormalization constants $Z_{[N,p][N,p^{\prime }]}$ can be illustrated within the one-loop approximation. Let $\mathsf{\Gamma }(x;\theta )$ be the generating functional of the 1PI Green functions with the insertion of just one composite operator $F[N,p]$ of Equation (224) containing any number of fields $\theta$. Here, x is the argument of the composite operator and $\theta$ is the functional argument, the “classical counterpart” of the random field $\theta$. The general interest is in the Nth term of the expansion of $\mathsf{\Gamma }(x;\theta )$ in $\theta$, which is denoted as ${\mathsf{\Gamma }}_N(x;\theta )$; it has the form

$${\mathsf{\Gamma }}_N(x;\theta )=\frac{1}{N!}\int \mathrm{d}{x}_1\cdots \int \mathrm{d}{x}_N\theta \left(x_1\right)\cdots \theta \left(x_N\right){\langle F[N,p]\left(x\right)\theta \left(x_1\right)\cdots \theta \left(x_N\right)\rangle }_{1-\mathrm{ir}}.$$

(226)

The matrix of critical dimensions in Equation (222) is given in the one-loop approximation by the expression

$${\Delta }_{[N,p][N,p^{\prime }]}=N\epsilon /2+{\gamma }_{[N,p][N,p^{\prime }]}^{*},$$

(227)

where the asterisk implies the substitution in Equation (211). Details of calculation of ${\gamma }_{[N,p][N,p^{\prime }]}$ can be found in [104].

As already mentioned, the critical dimensions are determined by the eigenvalues of the matrix in Equation (227). It is readily checked that, in the isotropic case (${\alpha }_{1,2}=0$), elements of the matrix with $p^{\prime }>p$ vanish, it becomes triangular, and its eigenvalues are the diagonal elements $\Delta [N,p]\equiv {\Delta }_{[N,p][N,p]}$:

$$\Delta [N,p]=N\epsilon /2+\frac{2p(p-1)-(d-1)(N-p)(d+N+p)}{2(d-1)(d+2)}\epsilon +O\left({\epsilon }^2\right).$$

(228)

From this equation, we observe that, for fixed N and arbitrary $d\ge 2$, the dimension $\Delta [N,p]$ monotonically decreases with p. It reaches a minimum at the minimal allowed value of $p=p_N$, i.e., $p_N=0$ if N is even and $p_N=1$ if N is odd:

$$\Delta [N,p]>\Delta [N,p^{\prime }]\mathrm{if}p>p^{\prime }.$$

(229)

This minimal value $\Delta [N,p_N]$ decreases monotonically as N increases separately for even and odd values of N, viz.

$$0\ge \Delta [2n,0]>\Delta [2n+2,0],\Delta [2n+1,1]>\Delta [2n+3,1].$$

(230)

A similar hierarchy is demonstrated by the critical dimensions of certain tensor operators in the stirred Navier–Stokes turbulence (see Ref. [107] and Section 2.3 of [31]). However, no clear hierarchy is demonstrated by neighboring even and odd dimensions: relations

$$\Delta [2n+1,1]-\Delta [2n,0]=\frac{\epsilon (d+2-4n)}{2(d+2)},\Delta [2n+2,0]-\Delta [2n+1,1]=\frac{\epsilon (2-d)}{2(d+2)}$$

(231)

lead to inequality $\Delta [2n+1,1]>\Delta [2n+2,0]$ for any $d>2$. On the other hand, the relation $\Delta [2n,0]>\Delta [2n+1,1]$ is satisfied only if n is large enough, $n>(d+2)/4$.

Let us denote $\Delta [N,p]$ the eigenvalue of the matrix in Equation (227), which coincides with Equation (228) for ${\alpha }_{1,2}=0$. Because the eigenvalues depend continuously on parameters ${\alpha }_{1,2}$, at least for small enough values of ${\alpha }_{1,2}$, this notation is unambiguous.

The dimension $\Delta [2,0]$ vanishes identically for any ${\alpha }_{1,2}$ in all orders in $\epsilon$. As in the isotropic model, this is proved employing the Schwinger equation

$$\int \mathcal{D}\mathsf{\Phi }\frac{\delta }{\delta {\theta }^{\prime }\left(x\right)}\left[\theta \left(x\right){\mathrm{e}}^{{\mathcal{S}}_R\left[\mathsf{\Phi }\right]+A\mathsf{\Phi }}\right]=0.$$

(232)

It can be easily checked [104] that $\Delta [2,0]\equiv 0$, while the to the leading order in $\epsilon$ one obtains

$$\begin{array}{cc}\hfill \Delta [2,2]/\epsilon & =2+\{-(d-2)d(d+2)(d+4)F_0^{*}-(d+2)(d+4)(2+(d-2){\alpha }_1\hfill \\ \hfill & +d{\alpha }_2)F_1^{*}++3(d+4)(d-2{\alpha }_1+2d{\alpha }_2)F_2^{*}+15d({\alpha }_1-{\alpha }_2)F_3^{*}\}/\{(d-1)\hfill \\ \hfill & \times (d+4)[(d-1)(d+2)+(d+1){\alpha }_1+{\alpha }_2]\},\hfill \end{array}$$

(233)

where $F_n^{*}\equiv F(1,1/2+n;d/2+n;-{\alpha }_{*})$ with ${\alpha }_{*}$ from Equation (211).

For $N>2$, the eigenvalues may be found in closed form only within an expansion in ${\alpha }_{1,2}$. The explicit expressions can be found in [104]. They illustrate two features which appear to hold for all N:

• The most relevant anisotropy correction is of order $\left({\alpha }_{1,2}\right)$ for $p\ne 0$ and $\left({\alpha }_{1,2}^2\right)$ for $p=0$. This means that ${\gamma }^{*}[N,0]$ are anisotropy independent in the linear approximation.
• This leading contribution depends on ${\alpha }_{1,2}$ only through the combination ${\alpha }_3\equiv 2{\alpha }_1+d{\alpha }_2$.

This conjecture is supported by the following expressions for $N=6,8$ and $p=0$:

$$\begin{array}{c}\hfill {\gamma }^{*}[6,0]/\epsilon =\frac{-2(d+6)}{(d+2)}-\frac{12{(d-2)}^2(d+1)(d^2+14d+48){\alpha }_3^2}{{(d-1)}^{2}d{(d+2)}^4{(d+4)}^2},\end{array}$$

(234)

$$\begin{array}{c}\hfill {\gamma }^{*}[8,0]/\epsilon =\frac{-4(d+8)}{(d+2)}-\frac{24{(d-2)}^2(d+1)(d^2+18d+80){\alpha }_3^2}{{(d-1)}^{2}d{(d+2)}^4{(d+4)}^2}.\end{array}$$

(235)

Using the operator product expansion (Section 5), we infer for the scaling function $R(r/L)$ in the representation in Equation (218) for the correlator $\langle F_1\left(x\right)F_2\left(x^{\prime }\right)\rangle$: the expression

$$R(r/L)=\sum _F{A}_F(\frac{r}{L}{)}^{{\Delta }_F},r\ll L,$$

(236)

where the coefficients $A_F$ are regular in ${(r/L)}^2$.

Now, let us discuss the equal-time structure functions $S_N$ in Equation (52). From this, it is assumed that the mixed correlator $\langle \mathit{v}f\rangle$ differs from zero; this does not influence the critical dimensions. However, it gives rise to non-vanishing odd structure functions. As a rule, the operators entering into the OPE are those generated by Taylor expansions together with all admissible operators which admix to them in the renormalization procedure [5,6]. The leading term of the Taylor expansion of $S_N$ is the operator $F[N,N]$ from Equation (224). Renormalization brings about all the operators $F[N^{\prime },p]$ with $N^{\prime }\le N$ and all allowed values of p. The operators with $N^{\prime }>N$ (whose contributions would be more relevant) are absent in Equation (236), because they do not appear in the Taylor expansion of $S_N$ and do not admix to the terms therein. Hence, the RG representation in Equation (216) together with the OPE representation in Equation (236) leads in the inertial range to the asymptotic expression for the structure function:

$$S_N\left(\mathit{r}\right)=D_0^{-N/2}{r}^{N(1-\epsilon /2)}\sum _{N^{\prime }\le N}\sum _p\left\{C_{N^{\prime },p}{(r/L)}^{\Delta [N^{\prime },p]}+\cdots \right\}.$$

(237)

The second sum is taken over all values of p, allowed for a given $N^{\prime }$, numerical coefficients $C_{N^{\prime },p}$ depend on $\epsilon$, d, ${\alpha }_{1,2}$ and the angle $\vartheta$ between $r$ and $n$. The ellipsis denotes the contributions of the operators other than $F[N,p]$, e.g., ${\partial }^2\theta {\partial }^2\theta$. They generate terms of order ${(r/L)}^{2+\left(\epsilon \right)}$ and higher and are discarded in the following.

A few general remarks are in order:

• If the mixed correlator $\langle vf\rangle$ is missing, the odd structure functions vanish. At the same time, the contributions to even functions are given only by the operators with even values of $N^{\prime }$. Only contributions with $p=0$ remain in the isotropic case (${\alpha }_{1,2}=0$) [104,108]. When anisotropy is present, ${\alpha }_{1,2}\ne 0$, the operators with $p\ne 0$ assume nonvanishing means, and their dimensions $\Delta [N^{\prime },p]$ enter the right hand side of Equation (237).
• The most relevant term of the small $r/L$ behavior is obviously related to the contribution with the smallest value of the exponent $\Delta [N^{\prime },p]$. Now, we employ the hierarchy relations in Equations (229) and (230), which hold for ${\alpha }_{1,2}=0$ and consequently remain valid at least for ${\alpha }_{1,2}\ll 1$. We conclude that, for sufficiently weak anisotropy, the leading term in (237) is given by the dimension $\Delta [N,0]$ for any $S_N$. In all particular cases considered, this hierarchy remains for finite values of the anisotropy parameters as well. The contribution with $\Delta [N,0]$ stays leading for such N and ${\alpha }_{1,2}$.
• It is possible that the inequalities in Equations (229) and (230) may be violated for some values of n, d and ${\alpha }_{1,2}$, so that the leading contribution to Equation (237) is given by a dimension with $N^{\prime }\ne N$ and/or $p>0$.
• The introduction of the mixed correlator $\langle vf\rangle \propto n\delta (t-t^{\prime })C^{\prime }(r/L)$ explicitly violates the evenness in $n$ and generates non-vanishing odd functions $S_{2n+1}$ and leads to to the contributions with odd $N^{\prime }$ to the expansion in Equation (237) for even functions. If the relations in Equations (229) and (230) are satisfied, in even functions, the leading term is still given by the contribution with $\Delta [N,0]$. If the relations in Equation (231) are fulfilled, in the odd function $S_{2n+1}$, the leading term is given by the dimension $\Delta [2n,0]$ for $n<(d+2)/4$ and by $\Delta [2n+1,1]$ for $n>(d+2)/4$. Let us note that, for the model with an imposed gradient, for all n, the leading terms of $S_{2n+1}$ are given by the dimensions $\Delta [2n+1,1]$ [34,36].

Furthermore, it is permissible that the matrix in Equation (227) for some ${\alpha }_{1,2}$ has a pair of complex conjugate eigenvalues, $\Delta$ and ${\Delta }^{*}$. In this case, the small $r/L$ behavior of the scaling function $\xi (r/L)$ in Equation (237) would contain oscillating terms of the form

$${(r/L)}^{\mathrm{Re}\Delta }\left\{C_{1}cos\left[\mathrm{Im}\Delta (r/L)\right]+C_{2}sin\left[\mathrm{Im}\Delta (r/L)\right]\right\},$$

with constant factors $C_i$.

Representations similar to Equations (218) and (237) can be readily written for arbitrary equal-time pair correlation functions, provided their canonical and critical dimensions are known. In particular, for the operators $F[N,p]$ in the IR region ($\Lambda r\to \infty$, $r/L$ fixed), we arrive at

$$\langle F[N_1,p_1]F[N_2,p_2]\rangle ={\nu }_0^{-(N_1+N_2)/2}\sum _{N,p}\sum _{N^{\prime },p^{\prime }}{(\Lambda r)}^{-{\Delta }_{[N,p]}-{\Delta }_{[N^{\prime },p^{\prime }]}}{R}_{N,p;N^{\prime },p^{\prime }}(r/L).$$

(238)

Here, the indices N and $N^{\prime }$ satisfy the inequalities $N\le N_1$ and $N^{\prime }\le N_2$. The indices p and $p^{\prime }$ assume all allowed values for given N and $N^{\prime }$. The small $r/L$ behavior of the functions $R_{N,p;N^{\prime },p^{\prime }}(r/L)$ is of the form

$${\xi }_{N,p;N^{\prime },p^{\prime }}\left(mr\right)=\sum _{N^{″},p^{″}}{C}_{N^{″},p^{″}}{(r/L)}^{\Delta [N^{″},p^{″}]},$$

(239)

with the constraint $N^{″}\le N+N^{\prime }$ and corresponding values of $p^{″}$; $C_{N^{″},p^{″}}$ are numerical coefficients.

Thus far, we have considered the particular case of the velocity correlator given by Equations (190) and (195). Explicit calculations [104] show that contributions to the renormalization constants (as a direct consequence, to the coordinates of the fixed point and the anomalous dimensions as well) come solely from the even polynomials in the expansion in Equation (192). This is why odd polynomials were absent in Equation (192) from the very beginning. Further, from the explicit form [104] of self-energy graph ${\Sigma }_{{\theta }^{\prime }\theta }$, it can be shown that only the coefficients $a_l$ with $l=0,1$ and $b_l$ with $l=0,1,2$ contribute to the constants $Z_{1,2}$ in Equation (204) and, consequently, to the basic RG functions in Equations (207) as well as Equation (208) to the coordinates of the fixed point in Equation (211). Hence, the fixed point in the generic model in Equation (192) is fully parameterized by these five coefficients. The higher coefficients appear only through the positivity conditions in Equation (194).

Further, for $\chi =0$, only coefficients $a_l$ with $l\le 2$ and $b_l$ with $l\le 3$ might contribute to the integrals $H_n$ and, therefore, to the one-loop critical dimensions in Equation (227). As a consequence, calculation of the latter is significantly simplified in the special case $a_0=1$, $a_1=0$ and $b_l=0$ for $l\le 2$ in Equation (192). The coordinates of the fixed point in Equation (211) coincide with those in the isotropic model. In particular, ${\alpha }_{*}=0$, and the anomalous exponents will depend on the two parameters $a_2$ and $b_3$ solely. In all analyzed cases, the general picture is akin to that outlined above for the case in Equation (195). For instance, the hierarchy of the critical dimensions, stated in the inequalities in Equations (229)–(231), persists also to this case. We conclude that the special case in Equation (195) represents adequately the main features of the generic model in Equation (192).

The exponents are related to the critical dimensions of composite operators in Equation (224) constructed from the scalar gradients. In contrast to the isotropic flow, these operators in the model under consideration mix in renormalization in such a way that the matrices of their critical dimensions are neither of diagonal nor triangular form. These matrices are calculated explicitly to the order $\left(\epsilon \right)$ However, their eigenvalues (anomalous exponents) can be calculated only as perturbative series in ${\alpha }_{1,2}$ (Equations (234) and (235)) or in a numerical fashion.

In the limit of zero anisotropy, the exponents can be related to definite tensor composite operators constructed from the scalar gradients, and exhibit certain hierarchy connected to the degree of anisotropy. It can be summarized as follows: the lower is the rank, the lower is the dimension and, therefore, the contribution to the inertial-range behavior is more important.

Leading terms of the even (odd) structure functions are caused by the scalar (vector) operators. For the case of finite anisotropy, the exponents cannot be related to individual operators (which are mixed in renormalization procedure), but the aforementioned hierarchy is present for all the studied cases.

A short comment about the second-order structure function $S_2\left(r\right)$ is appropriate here. It can be studied employing the RG and zero-mode techniques [104]; as in the isotropic case [39,109,110,111]. Its leading term is $S_2\propto r^{2-\epsilon }$, but the amplitude now depends on ${\alpha }_{1,2}$ and the angle between the vectors $r$ and $n$ from Equation (195). The first correction due to anisotropy is ${(r/L)}^{\Delta [2,2]}$ with the exponent $\Delta [2,2]=O\left(\epsilon \right)$ from Equation (233).

In isotropic velocity field, the anisotropy inserted at large scales by external forcing or imposed mean gradient, persists in the inertial range and shows in odd correlation functions: the skewness factor $S_3/S_2^{3/2}$ decreases at $r/L\to 0$ but slowly (see [50,85,86,87,88,100,101,102,103]), while the higher-order ratios $S_{2n+1}/S_2^{n+1/2}$ increase (see, e.g., [34,36,90,92]).

In the studied model, the inertial-range behavior of the skewness factor is described by the expression $S_3/S_2^{3/2}\propto {(r/L)}^{\Delta [3,1]}$. For ${\alpha }_{1,2}\to 0$, the exponent $\Delta [3,1]$ is given by Equation (228) with $n=3$ and $p=1$. It is positive and coincides with the previous result [85,86]. By numerical means, it can be shown that, if the anisotropy becomes strong enough, the exponent $\Delta [3,1]$ becomes negative and the skewness factor increases towards the depth of the inertial interval. Let us note that the higher-order odd ratios increase already when the anisotropy is weak.

11. Conclusions

The crucial problem in many phenomena encountered in physics is the proper identification of underlying symmetry and the mechanism that leads to its violation. The solution to such a task is by no means obvious and requires deep understanding of the physical problem. In non-equilibrium statistical physics, which deals with systems containing many interacting degrees of freedom, symmetries are present at different levels of theoretical description.

In this article, our main aim is to summarize and elucidate approaches suitable for theoretical analysis of symmetries and their violations in problems related to various aspects of fluid dynamics. In particular, we concentrate on an overall analysis of symmetries emerging in fully developed turbulence and related problems. All studied problems are of classical nature and share the common property that they cannot be described within the standard equilibrium statistical physics and thus constitute systems far from the thermal equilibrium. Nevertheless, it is possible to study them with the use of powerful and versatile methods borrowed from high energy physics. We demonstrate how to construct a functional-integral representation of classical stochastic problems based on the introduction of an action functional analogous to that of particle physics. Use of the action functional is especially illuminating in revealing symmetries of models and their consequences. The functional representation is analyzed using Feynman diagrammatic technique and renormalization group approach. The latter is especially convenient for an analysis of large-scale behavior. To this end, it is necessary to find R stable fixed point(s) of the RG $\beta$ functions of the model at hand. The fundamental difference between the RG in high-energy physics and statistical systems is that only the opposite spatial scales are of interest. In high-energy physics, small spatial scales are studied and the asymptotic behavior of models at large momenta is vital. In classical statistical systems, there is a maximal spatial scale L. The goal of the theory is then to determine how much L affects relevant physical quantities. We apply field-theoretic methods to classical non-equilibrium systems and as a paradigmatic example present the stochastic approach to fully developed turbulence. We review the basic setup, and briefly describe the RG method and the related topic of operator product expansion. Examples are given of the use of the latter as an important tool in the analysis of experimentally measurable quantities. The important Galilei symmetry is described in the functional representation and its consequences derived. Consequences of breaking of Galilei symmetry in the stochastic Navier–Stokes model with colored noise are discussed with the use of the RG. Important effects of parity violation and strong anisotropy in turbulent flows are analyzed in detail.