1. Introduction
Numerous important natural phenomena, spanning a wide range of spatial scales from microscopic to cosmological, are associated with hydrodynamic flows [
1,
2,
3,
4]. At all these scales, fluids can exist in a wide variety of states, each presenting an important challenge from both an experimental and a theoretical point of view. Especially intriguing behavior is observed for turbulent flows [
5,
6], which are generally considered rather a rule than an exception in nature [
3]. Their study represents an exceptionally important task both from a theoretical perspective and from the standpoint of applications. It is noteworthy to mention several important areas where turbulence plays a pivotal role, including diverse challenges in fluid dynamics (such as turbulence in boundary layers [
7], atmospheric turbulence [
8,
9], etc.), issues related to turbulent mixing and homogenization [
10], the study of fusion plasma [
11], and various astrophysical phenomena [
12,
13], among others. However, despite the enormous efforts that have been made to study turbulence for more than a century, the problem remains unsolved.
The modern theoretical approach to turbulence is grounded on the statistical analysis of the solutions to the Navier–Stokes (NS) equation, which is generally accepted as describing the motion of viscous (nonrelativistic) media [
1]. As already noted, the medium during its movement can exist in various states depending on the flow parameters, for which the Reynolds number
is typically introduced for classification. Defined as
, where
V is the typical average flow velocity,
L is a typical macroscopic scale, and
is the kinematic viscosity of the medium,
represents the ratio between inertial and dissipative forces. At low values,
, a regular (laminar) flow is observed. With an increase in
, phenomena ranging from periodic, such as von Kármán vortices, to highly chaotic and irregular motion at very high values of
(practically,
is considered large enough) occur. The latter state of fluid is known as fully developed turbulence. In this state, the so-called inertial interval exhibits a transfer of kinetic energy from macroscopic
L to microscopic (dissipative)
ld scales. For fully developed turbulence, the symmetry and dimension arguments have led to important conclusions about the scaling behavior of velocity correlation functions, primarily through the famous Kolmogorov’s K41 theory. However, more detailed statistical descriptions of both the fully developed turbulence regime and its onset in laminar flow are still needed. Mathematically, this research gap stems from the well-known blow-up problem: the existence of a strong instability inherent in solutions of the NS equations at high Reynolds numbers. Generally speaking, for this equation, even the existence of a global smooth solution is in question at present [
14]. Moreover, both mentioned problems are still unsolved not only for the Navier–Stokes equation but even for the simpler Euler equation (inviscid fluid).
Traditionally, the physical theory of developed turbulence models this instability effect by the addition of a suitably chosen random force to the NS equation. This addition also simulates the continuous input of mechanical energy necessary to sustain the turbulent cascade. The selection of the structure and statistics of this random force is the most critical aspect of the theory. Typically, random forces concentrated at large spatial scales are used to model large-scale effects on the scale-invariant behavior predicted by Kolmogorov’s theory. The standard stochastic problem formulation also assumes rotational symmetry of force correlations but can be expanded by varying the symmetry properties of force correlations (e.g., anisotropy and reflection asymmetry), to explore the impact of large-scale properties on velocity correlation scaling behavior [
15].
In the context of developed turbulence studies, investigating the advection of specific quantities (such as temperature fields, concentration fields, or tracers) [
4,
16] by the turbulent field is paramount. The solution to this problem should shed light on the extent of turbulence intermittency [
5], i.e., its fractal nature. In this case, the Prandtl number
is often used to succinctly describe the flow’s quantitative characteristics [
8,
17]. For all types of admixtures,
is defined as the dimensionless ratio of the kinematic viscosity coefficient
to the corresponding diffusion coefficient
D of a given admixture. Since both the kinematic viscosity and the diffusion coefficient are material- and flow-specific, resulting Prandtl numbers have to be specified under conditions characterizing the flow, hence are often found in property tables along with other material-specific properties [
18,
19]. Yet, in the limit of high Reynolds numbers, fully developed turbulence manifests through effective values of kinematic viscosity and the corresponding diffusion coefficient, rendering the turbulent Prandtl number independent of material and flow specifics. In what follows, we call such
a turbulent Prandtl number [
6,
8].
Among general turbulence theory issues, describing fully developed turbulence in an electrically conducting fluid occupies a special place. The system of equations of motion for an electrically conducting fluid is known as the magnetohydrodynamics (MHD) equation system (see, e.g., [
1,
10,
17,
20]). For the first time, this system (in the case of an inviscid fluid) was investigated by Hannes Alfvén. For the subsequent fruitful application of it to problems in plasma physics, Alfvén was awarded the Nobel Prize in 1970. Nowadays, there are many theoretical methods for studying this system of equations, which constitute the golden fund of theoretical and mathematical physics. As is well known, in electrically conducting fluids, turbulent motions are accompanied by magnetic field fluctuations. However, such conducting fluids, exemplified by ionized gas (plasma) or magma, are rare under typical laboratory conditions. Conversely, astrophysical or geophysical settings, such as the solar wind or turbulent dynamo mechanisms generating large-scale magnetic fields, underscore the enduring interest in MHD turbulence studies [
12,
21,
22].
From a chronological perspective, the Kazantsev–Kraichnan kinematic model holds paramount importance in magnetohydrodynamics [
23,
24,
25]. The main idea of this model is the assumption that while the magnetic field is passively advected by the hydrodynamic flow, the reciprocal effect on the main flow from the magnetic field is negligible. However, this assumption falls short in describing realistic flows. A genuine model of magnetohydrodynamics must account for the interplay between the magnetic and velocity fields. Incorporating this interaction, along with spatial parity nonconservation, suggests the possibility of generating a uniform large-scale magnetic field within the system—a mechanism known as the turbulent dynamo [
20]. Many studies have focused on this phenomenon, and here, we provide references only to those that are close in language to this paper [
26,
27,
28].
As well as Kraichnan-like models, the MHD model should reveal scaling behavior for accompanying quantitative parameters such as correlation or structure functions. A quantitative description of these characteristics poses a natural challenge for the theory to explain and construct reliable approximation schemes for their calculation. As demonstrated by Kolmogorov’s theory, scaling invariance naturally emerges in turbulent flows [
5]. However, exploring behaviors such as anomalous scaling demands extensive, meticulous analysis. This phenomenon is characterized by singular (arguably, powerlike) behavior of some statistical quantities (correlation functions, structure functions, etc.) in the inertial–convective range in the fully developed turbulence regime [
3,
5]. Moreover, it is worth noting that scaling properties manifest themselves more clearly at higher values of Reynolds number
[
29], making the
limit critically important. Thus, employing the well-established renormalization group (RG) theory tools for analytical results on turbulent flows is prudent. Unlike numerical approaches, where achieving the
limit is difficult, the RG method allows for analytical calculations, ultimately allowing for the examination of universal turbulent phenomena, i.e., independent of material microstructure and macroscopic flow conditions [
15]. Additionally, the RG framework readily accommodates various scenarios, including the advection of different admixtures [
29].
The renormalization group approach was previously applied to the kinematic MHD at the level of one-loop approximation in quantum field perturbation theory [
27]. The authors also considered the case of so-called chiral or gyrotropic turbulent media, when the corresponding magnetic fields and velocity fields lack defined spatial parity. Technically, the presence of chirality means that all correlation functions’ transverse structures (assuming the medium is incompressible and magnetic field transversality is ensured by Maxwell’s equations) become a mixture of two quantities: tensor and pseudotensor (with respect to the three-dimensional rotation group
). It has been shown that spatial parity nonconservation, even without magnetic fluctuations, leads to spontaneous symmetry breaking tied to the instability of linearized MHD and the emergence of a large-scale homogeneous mean magnetic field
. This field potentially stabilizes the system in its new “vacuum state”, acting through the turbulent dynamo mechanism. From a technical standpoint, spontaneous symmetry breaking generates new counterterms of mass renormalization type (here, a homogeneous magnetic field gradient plays the mass role) in a form
, nontrivially arising during the ultraviolet (UV) renormalization procedure, starting from one-loop diagrams. Here,
ld represents the naturally available cutoff scale in the theory. However, as demonstrated in [
30,
31], the appearance of these counterterms does not compromise the theory’s multiplicative renormalizability, in a sense that their contribution can always (at least at the one-loop level) be mitigated by an appropriate choice of
, presumably stabilizing the system. Beyond merely positing a stabilizing magnetic field, the paper [
30] predicts a new physical effect: the emergence of disturbances in Alfvén waves perpendicular to the spontaneous field
, leading to their linear temporal growth. This results in long-lived pulses
,
, akin to the massless bosons in the Goldstone model.
In essence, numerous vector impurity models bear similarities to kinematic MHD to varying degrees. As repeatedly noted, the tensor nature of considered impurities significantly influences diffusion and advection processes [
31,
32,
33,
34]. For instance, as has been shown in [
33], introducing spatial parity violation (helicity) into the turbulent flow not only offers a highly realistic physical scenario but also uniquely highlights the studied model’s specific tensorial and interactive properties. Consequently, the authors infer that interaction structures impact diffusion and advection processes more than the tensor nature of the advected field itself. To extend these findings beyond a few selected models, examining the general
A model of passive vector advection, initially defined in [
35], proves useful. To prevent confusion, it should be emphasized that the current general
A model exclusively pertains to passive vector advection and should not be confused with the model
A according to the classification of Halperin and Hohenberg [
36] or with the
A model of the forced Navier–Stokes equation presented in [
37,
38]. The model’s name stems from the widely accepted notation for the
A parameter [
28,
31,
39], which appears in its definition (refer to
Section 2 for details). The
A parameter plays a central role in the appropriate definition of a unified description of a wide range of vector admixture models. Given its requirement to be real without explicit restrictions, it can adopt values
,
, and
, correlating to the kinematic MHD model, the model for passively advected vector fields, and the scenario of the so-called linearized Navier–Stokes equation [
39], respectively. Thus, the general
A model serves as a tool to unite several distinct but physically important cases into a single model. Overall, for this model, a comprehensive generalization to all physically permissible values of
A (including beyond the boundaries of the interval
) was conducted by the authors of [
28]. However, the impact of helical effects in a fully turbulent environment was confined to what was necessary for calculating the turbulent Prandtl number. Therefore, to obtain a comprehensive picture of MHD, it seems highly desirable to calculate the full set of renormalization constants for the general parametric
A model considering the effects of spatial parity violation up to the two-loop order in perturbation theory.
To conduct the above-discussed studies, this endeavor employs quantum field theory methods alongside the established renormalization group technique tools (for instance, see [
15,
29,
40]). The RG method has found extensive application in fully developed turbulence research without admixtures [
30,
41,
42,
43], as well as in studying advection–diffusion processes involving various admixtures, including passive scalar admixtures [
35,
43,
44,
45,
46,
47], magnetic admixtures [
27,
45,
48], and also vector admixtures [
33,
39,
45]. Highlighting the nuances of applying standard quantum field techniques to stochastic hydrodynamics issues is also pertinent. Among various regularization techniques in field theory, a distinguished role is played by dimensional regularization accompanied with
-expansion [
29,
40]. This approach directly yields regularized counterparts for divergent Feynman diagrams, with the associated
-expansion serving as a valuable perturbative method for calculating universal quantities in critical phenomena and stochastic dynamics theory [
29,
40]. It is worth mentioning that the
parameter’s interpretation within stochastic turbulence theory differs from its role in critical phenomena models (e.g., well-known
-theory), where the small
-expansion parameter is usually determined by a difference
from the upper critical dimension [
29,
40]. In the stochastic theory of developed turbulence, its definition has no relation to the dimension of space. Most studies integrate
into the scaling behavior of random forces, strictly considering a three-dimensional context [
5,
15]. The
parameter thus derives from the logarithmic theory where
, aligning with the ideal scenario of (delta-shaped) energy infusion into the system from infinitely vast spatial scales [
15,
41]. Once a model undergoes correct renormalization, perturbative results for universal quantities, such as critical (scaling) exponents, become attainable. These invariably take the form of asymptotic series. Despite the existence of impressive multiloop outcomes in critical and stochastic dynamics [
49,
50,
51], the developed turbulence RG theory predominantly limits analytical calculations to two-loop approximations [
29,
52,
53]. In fact, it might be argued [
42] that the complexity of such endeavors is arguably an order of magnitude greater than first-order analyses. Not only is the number of Feynman diagrams much higher but there are also intrinsic numerical problems with the correct extraction of divergent parts of Feynman diagrams. Therefore, any attempt that might help to make two-loop calculations feasible seems worthwhile to undertake.
A comprehensive two-loop calculation for the dynamo effect in stochastic MHD has not yet been finished. Our goal herein is to make progress in this direction and present particular results of two-loop approximation. In particular, we present details on the two-loop calculation of a three-point interaction vertex, which, by itself, presents quite demanding technical problems. To make the problem tractable, it is necessary to rely on powerful algorithmic and numerical approaches [
54,
55,
56].
This paper is structured as follows.
Section 2 defines the
A model of passive advection of a vector admixture as a generalization of the stochastic MHD and introduces its field theoretical counterpart.
Section 3 is devoted to a detailed renormalization group analysis. The final
Section 4 is devoted to concluding remarks and plans for the future.
2. Formulation of the Quantum Field Model
It is well-known [
1] that the hydrodynamics of a nonrelativistic incompressible viscous fluid is described by the Navier–Stokes equation
Here, represents the velocity field, stands for the (transverse) random force exerted per unit density, and scalar function denotes the pressure per unit density. Throughout the text, the argument of all field functions is defined as , where is a d-dimensional spatial vector and t is the time. The parameter represents the kinematic viscosity, the symbol ∇ is the spatial gradient, is the time derivative, is the Laplace operator, and parentheses indicate the standard dot product of three-dimensional vectors. Let us also stress that for possible verification and transparency in intermediate calculations, we retain general dimension d in actual expressions, although in the end, we are mainly interested in arguably the most relevant three-dimensional case . In other space dimensions, the used approach, strictly speaking, will not work.
In the approximation of an incompressible fluid (low Mach number), the explicit dependence on the pressure in (
1) can be simply eliminated by taking the divergence of (
1). The result is the following elliptic equation,
from which it follows that the pressure can be formally represented in terms of the fundamental solution of Equation (
2). The symbol ⊗ here represents the outer product. Substituting this solution into (
1), it is easy to see that the purely longitudinal contribution of pressure
reduces the longitudinal part of the value
. Therefore, the pressure term in (
1) can be eliminated by discarding the contribution
and at the same time placing a transverse projector in front of
.
In the conventional field theory setting, the problem in (
1) is analyzed throughout the entire space and along the whole time axis, which requires zero divergence-free boundary conditions and the implementation of a retardation condition, which entails that the velocity field
vanishes as
. Note that, generally speaking, these assumptions do not fix either the existence or the uniqueness of the solution of (
1). A comprehensive overview of the current developments in the mathematical theory of the Navier–Stokes equations can be found, e.g., in [
57].
Particularly, developed turbulence across a sufficiently broad range of parameters is described within the approximation of an incompressible fluid. This is adequate as long as the flow velocity is much less than the speed of sound [
58,
59]. As already noted in the introductory
Section 1, investigating the turbulent regime at the level of numerical solutions of Equation (
1) is extremely complicated by the fact that in this regime, all scales (mixing lengths) are interconnected in a cascade-like manner. Consequently, the flow in the boundary layer (the scale at which the pressure gradient becomes comparable to the viscous term) affects everything up to the overall scale of the system. This “entanglement of scales” implies that inevitable tiny errors at the smallest scales of the approximation grid will ultimately have a huge impact on the overall structure of the solution. For example, in the simplest idealized model of the atmosphere (considered as a two-dimensional ideal fluid on the surface of a torus), deviations increase by
times in two months. Clearly, this fact alone renders dynamic weather forecasting for such a duration practically impossible, regardless of the computational power or the density of the data network employed. Therefore, in practical calculations, people often entirely forego solving the NS equation, describing the fluid flow not hydrodynamically but kinetically, i.e., through an ensemble of interacting particles with the same statistics as the real medium (see, for example, the lattice Boltzmann method [
60,
61]). By averaging such a lattice system over all particles in the continuum limit, the Navier–Stokes equations are obtained. However, the ensemble itself can be modeled without any approximations, which yields the correct fluid flow without solving any equations.
Taking these difficulties into account, the physical theory of developed turbulence follows a different route based on the stochastic version of Equation (
1) (see, e.g., [
3,
5,
8,
15]). More specifically, in this approach,
in Equation (
1) is generalized to the random variable with prescribed stochastic properties. This modification aims to simulate physical processes, such as the injection of energy into the system and the influence of boundaries. Consequently, the theory’s main objective shifts to identifying the statistical characteristics of the (stochastic) velocity field
.
Let us turn our attention to magnetohydrodynamics [
12,
21,
22], where, as is known, any movement of a conducting fluid in a magnetic field generates an electromotive force, inducing electric currents. Due to the presence of the magnetic field, these currents create mechanical forces that alter the state of motion of the fluid. Accordingly, when examining the hydrodynamics of an electrically conducting fluid, Equation (
1) has to be suitably modified to incorporate the interaction between the velocity field and the magnetic field and complemented by an appropriate equation of motion for the magnetic field. More specifically, a hydrodynamic description here implies the adoption of three fundamental assumptions: firstly, the fluid should be considered highly conductive; secondly, its motion is nonrelativistic; and finally, the bulk variables must change slowly over time compared to all relevant characteristic time scales. Summarizing the above, magnetohydrodynamics is the simplest sufficient approximation to describe many large-scale low-frequency phenomena. Experience shows that MHD is applicable across a vast range of spatial scales, from nanometer (
m) scales in, e.g., physics of semiconductors, to interstellar and galactic (
m) scales, e.g., galactic arms.
Let us clarify here the theoretical setup needed to write the system of MHD equations in the case of a viscous liquid with finite conductivity. In general, magnetohydrodynamics studies the interaction of the electromagnetic field with a liquid or gaseous moving conductor, treated as a continuous medium. The MHD equations are a combination of Maxwell’s equations for the electromagnetic field and the conventional hydrodynamic equations that describe the motion of the continuous medium. The hydrodynamic equations of an electrically conductive fluid are valid under the conditions of the so-called MHD limit. Broadly speaking, this corresponds to the limit in which the charge and volume densities are obtained from the fluid equations rather than from the Boltzmann equation [
12]. Assuming this limit to be fulfilled, let us derive the equation of motion for the magnetic field in the MHD system. When the fluid is moving at the velocity
, the corresponding electric current density is connected with the electromagnetic fields through generalized Ohm’s law [
22]. Following [
27], we use its simplest form here, omitting the Hall term and others:
Here,
and
are the electric and magnetic fields in the “laboratory” system of coordinates, where we measure the velocity
, and
is the electrical conductivity, which is assumed to be finite but large enough, i.e., the electromagnetic processes are not very fast. The last assumption regarding the value of
allows us to omit displacement currents (and convection currents) compared to conduction currents; therefore,
in Formula (
3) denotes precisely the conduction current [
12,
22]. Then, taking the curl of Equation (
3) using the Maxwell–Faraday equation
, one can obtain
Utilizing another of Maxwell’s equations,
, with
c representing the speed of light and
denoting the vacuum magnetic permeability, and noting that in the case of nonrelativistic motion, the time derivative of the electric field is much smaller than the spatial derivative of the magnetic field, we can eliminate
in favor of
and finally derive the governing dynamic equation for the magnetic field in the form
Here, the parameter
is introduced, which has the meaning of the magnetic diffusion (or magnetic viscosity) coefficient. By elementary manipulations, Equation (
5) can be reduced to the form
which presents a starting point for a majority of studies in the MHD area. For a more detailed exposition, the interested reader is referred to the literature [
12,
20].
Let us now go back to the dynamics of the velocity field in the presence of a magnetic field. To describe it, we will use the Navier–Stokes equation (
1), where the Lorentz force contribution (per unit density)
is explicitly allocated in the force term
. Applying further reasoning similar to that used to obtain Equation (
5), we obtain
Rewriting the first term on the right side of (
7), we arrive at the final desired generalization of the Navier–Stokes equation:
where the second term on the right side is the total pressure, consisting of two terms: gas (or thermal) pressure
p and the so-called magnetic pressure
.
A stochastic version of the MHD model is obtained by generalizing the external force term in the Navier–Stokes equation (
8) to a random variable
, as well as adding the corresponding random force term
modeling a magnetic noise to Equation (
6). Stochastic random forces
and
serve to mimic microscopic effects such as the presence of boundaries or interactions between the mean velocity flow and its fluctuating part [
29]. Notwithstanding their obvious physical relevance, it is widely believed that their actual form is unimportant for interesting universal quantities. Note that from now on, variables
and
in Equations (
6) and (
8) should already be understood as the fluctuation parts of the corresponding total magnetic field and velocity field [
26,
44]. Furthermore, as has been pointed out in introductory
Section 1, here we will consider a generalization of Equation (
6) by introducing a new parameter
A into the nonlinear term
. Such a generalization corresponds to the transition to the general model A [
28,
31]. Note that the conventional formulation of MHD corresponds to the value
[
26,
27]. Eventually, the stochastic MHD can be formulated by two intercoupled stochastic differential equations [
28,
31]
with material (covariant) derivative
, complemented by transversality conditions for all vector fields
In Equations (
9)–(
11), fields
and
are the stochastic magnetic and velocity fields,
, where
is a reciprocal magnetic Prandtl number (by definition, magnetic Prandtl number
). The scalar functions
and
have the meaning of total pressure (see Equation (
8)) and have no significance for further analysis for reasons similar to the text after Formula (
2) and will be omitted. In general, within the framework of homogeneous and isotropic developed turbulence, the random velocity field in Equation (
8) is represented as
, where
is the constant mean flow velocity and
is a relatively small stochastic (fluctuating) component. The statistical characteristics of the latter are the subject of our study, since
can always be eliminated by transitioning to a moving reference frame through an appropriate Galilean transformation
. Given the smallness of
compared to
, any explicit dependence on
t in
is neglected in comparison to
. In terms of the correlators
, this means that the average value of the corresponding dynamic correlator over time
is nothing other than the value of the static (simultaneous) correlator for the distance
. For brevity, we will henceforth denote
as
, always understanding it to refer specifically to the fluctuating component. In contrast to
, the average value of the random field
is usually assumed to be zero, and the corresponding uniform magnetic field required for mass applications is introduced as an additive to the source
.
The influence of Galilean symmetry in the model governed by Equations (
9) and (10) extends beyond the aspects previously discussed. This symmetry encapsulates the foundational physics of the model and has significant implications for the RG analysis (see
Section 3 and [
31]). Moving forward, we interpret the Galilean invariance of Models (
9) and (10) in a broader context, specifically, as symmetry under what we term the generalized Galileo transformation (see, e.g., [
15,
31]). This concept goes beyond a mere constant speed transformation
to include transformations with an arbitrarily variable speed
that rapidly diminishes as
. Being essentially similar to gauge transformations, the generalized Galilean transformation is structured as follows:
Hereinafter,
denotes the Heaviside step function. Transformations (
12) are considered here as an independent Lie group and not as a representation of the corresponding group of coordinate transformations. They generalize the usual Galilean transformations [
1], in which the velocity
is constant and
.
Let us briefly discuss the physical relevance of the parameter
A [
28,
31]. First, we note that the aforementioned Galilean symmetry (
12) demands
A only to be real, yet the interval
already includes quite diverse, physically significant models. As already mentioned in
Section 1, choosing
leads us to the kinematic MHD model [
27]. On the other hand, the choice
leads to passive advection of a vector field in turbulent environments [
35,
39], and finally
yields the model of the linearized Navier–Stokes equations [
39]. Thus, the parameter
A stands in front of the stretching term, and, due to its continuous nature, represents a measure of specific interactions allowed by Galilean symmetry. Varying
A allows for the investigation of many passively advected vector admixtures with different interaction properties. Although
A can take any real value, it is often discussed only within the smallest possible continuous interval that encompasses these three mentioned distinguished cases. However, as shown in [
28], in cases where spatial parity conservation is violated in the model, in general, there are conditions under which it can be considered for all real
A, without restrictions. The general conclusion is that adopting
A as a general parameter of the theory allows for the simultaneous investigation of a multitude of diverse diffusion–advection processes. Particularly intriguing behavior is observed in processes occurring in MHD (
) that possess additional instability in helical environments, which can be stabilized on large scales by the emergence of a macroscopic field
[
28,
30].
The cornerstone of the stochastic theory of developed turbulence is the choice of the random force correlator, along with the corresponding dynamic equations that define the model [
29]. In the general case, in the system under consideration, all elements of the correlator matrix
,
, and
are different from zero. In particular, this situation was considered in [
27], applicable to magnetohydrodynamics. Let us start with the correlation function
of a magnetic noise
, which is assumed to be a transverse Gaussian random vector with zero mean
where
is an integral scale related to the stirring of the magnetic field
and function
is finite in the asymptotic limit
, while in the limit
, it should rapidly decrease. As noted above, the requirement of time decorrelation in
and other correlators is necessary to ensure the Galilean invariance of the model. A detailed form of
is not relevant for the actual RG analysis [
15]. In a physically more realistic formulation, the noise term
could be substituted by something like
, where
represents a strong (≈
–
G), large-scale magnetic field [
45,
62].
Let us stress that, in contrast to the original work [
27], our study explicitly omits both the correlators
and
. Here, we refer to a more general formulation of the problem considered previously in [
27], which shows that both
and
can be discarded as they are zero at the kinetic (Kolmogorov’s) fixed point which we are interested in.
Our primary interest lies in the correlator
. The transverse random force
describes the injection of kinetic energy into the turbulent system on a large spatial scale. Therefore, the actual form of its correlator should match the expected behavior of energy injection, which is supposed to be concentrated in the infrared (IR) region. It turns out that this choice can be reconciled with the RG approach [
29]. It is important to note that the universality of fully developed turbulence is not confined to the specific statistical properties of the random force, including the case with broken spatial parity. Hence, without loss of generality, we prescribe the pair correlation function of zero-mean Gaussian statistics
where the function
of
is the so-called pumping function and represents the pumping rate, and
ld (
ld is the dissipative length) is a typical ultraviolet turbulence scale. The tensor quantity
will be specified later. Let us also note that in stochastic turbulence theory, we also require that
correctly models the physical injection of energy (IR injection). Since for the problems of the theory, the exact macroscopic structure of a turbulent flow is negligible, the simplest cutoff can be taken as the corresponding IR regularization. Technically, this can be implemented by introducing the step function
into the pumping function
, where the parameter
, and
L is the integral scale of turbulence (size of the system). Thus, by adopting such a sharp IR cutoff, we effectively regularize the model in the IR region. It is also worth noting that regularization by cutoff in our case does not break the global symmetry group of Model (
12), unlike in relativistic field theory models.
To apply the standard RG technique effectively, it is crucial for the function
to exhibit a power-law behavior as
k becomes large [
29]. In particular, this requirement is fulfilled when
acquires the specific power-law dependence
Here,
is a coupling constant of the theory,
m represents IR cutoff, and
is the model parameter characterizing the type of pumping. For generality, we write expressions everywhere for an arbitrary dimension of space
d, but we are mainly interested in the case
. It is also worth mentioning that, strictly speaking, such a representation is valid only in the limit
(see, e.g., [
15]). Let us also stress that Formula (
14) represents the correlator in a kind of natural condition: there are two direct physical cutoff scales
L and
ld. Such a form is convenient for the correct setting of the problem (so-called
-renormalization of the model) but not for practical purposes. In practice, when calculating all necessary quantities, the limit of
is always taken, transferring the regularization function to the parameter
.
In a normal fluid (i.e., fluid invariant under Galilean transformations), to each quantity we can assign a specific spatial parity:
and
are vectors, whereas
and
are pseudovectors, and correlator
is a tensor. Due to the additional requirement of transversality (
11), any tensor is necessarily proportional to the transverse projector
, and any pseudotensor is proportional to the structure
, where
is the Levi-Civita tensor. The transition to a helical fluid corresponds to the abandonment of spatial parity conservation and technically implies that the correlator (
13) should be given by a sum of a true tensor and a pseudotensor, i.e., the tensor quantity
in correlator (
13) corresponds to the sum of the transverse projection operator
and a so-called helical term
(summation over repeated index
l is implied) [
30]. Hence, we choose the quantity
in the following form:
Here, the helicity parameter
is constrained within the range
, expressing the positive definiteness of the correlator (
13). It can be naturally interpreted as the amount of reflection symmetry breaking: the case
corresponds to the nonhelical case, and
is a state with maximal broken spatial parity.
Nonlinear stochastic differential equations (
9) and (10) are notorious for their mathematical complexity. However, their behavior in specific circumstances can be effectively analyzed by methods of quantum field theory. In particular, the RG method provides a natural framework for the analysis of macroscopic behavior, i.e., the behavior of various correlation functions in the infrared limit. In order to apply the RG method, it is advantageous to recast the Langevin-like formulation of the model in terms of path integrals. Following the well-known Janssen–De Dominicis procedure [
29,
53], stochastic equations (
9) and (10) are fully equivalent to the field-theoretic model of a double set of (unrenormalized) fields
, where
and
are Martin–Siggia–Rose (MSR) response fields [
63]. For the ensuing application of the RG method, it is necessary to distinguish between bare (unrenormalized) and renormalized parameters [
29,
40]. Therefore, we denote the bare fields with a circle on top “˚”, and later on, we write renormalized parameters and fields without any symbol whatsoever.
The field-theoretic model is then given by the De Dominicis–Janssen–Martin–Siggia–Rose (unrenormalized) action functional [
28], which takes the compact form:
Hereinafter, we employed a condensed notation, in which all integrals over space-time are not explicitly written, and summations over repeated Latin indices
are implicitly assumed. Thus, for instance, the expression
in the action (
16) actually stands for
where
denotes the
i-th component of the gradient.
The starting point of perturbation methods in quantum field theory is to express Green functions in terms of a corresponding sum of Feynman diagrams, whose relevance is controlled by interaction parameters (coupling constants). Throughout this text, by Green functions, we mean both correlation and response functions of the fields. These functions are derived from the functional averages of specific field products, weighted by the exponential of the action
. In graphical terms, these objects can be represented by Feynman diagrams constructed according to the standard diagrammatic rules of quantum field theory. The Feynman diagram technique serves as a convenient graphical tool for organizing various algebraic structures that arise from the basic perturbation theory elements (propagators and interaction vertices) [
29]. The graphical representation of these elements can be straightforwardly obtained. For completeness, we also present the necessary expressions here. In the frequency–momentum representation, the nonvanishing propagators of the model are expressed as
and are illustrated in
Figure 1.
The interaction part in (
16) gives rise to three interaction vertices:
,
and
, depicted in
Figure 2, where
and
stands for the momentum entering the vertex via the auxiliary field (
for
and
for
and
).
In actual calculations, it is often more suitable to work preferably with so-called connected Feynman diagrams, i.e., those diagrams that remain connected even if one of the internal lines is cut [
29,
40]. They can be obtained from the generating functional
defined as
by taking the corresponding number of functional derivatives with respect to external sources
J and then substituting
. For dynamic models, the symbol
, where
is the complete set of fields of the model and
are their MSR-conjugated analogs.
For translationally invariant theories, additional simplification consists of introducing functional
for one-particle irreducible (1PI or vertex) functions [
29,
40]. This is an even more restricted class than that of connected Feynman diagrams since they are obtained from all possible diagrams by throwing away those Feynman graphs which fall into disjoint parts when one internal line is cut. Then, for each external leg in a diagram, we divide out one factor of the propagator. In this way, we arrive at so-called proper vertices [
40,
64] that are convenient for the RG analysis. It can be shown that [
29] the corresponding functional for 1PI diagrams can be obtained by means of functional Legendre transform with respect to external sources
where
x denotes the set of all discrete and continuous arguments of the field. Within the framework of the approach used, the set of independent variables
for the functional
is conveniently divided as follows:
. Then, the 1PI functions
are obtained by differentiation with respect to fields
Let us also note that for reasons of causality, all 1PI Green functions of the form
vanish [
29]. In the next section, we are primarily interested in the quantities of
(for velocity fields and magnetic fields, respectively), as well as
for the vertex
U in
Figure 2, since they are the only ones responsible for the renormalization of the model.
3. Renormalization and Results
To obtain relevant information about the macroscopic behavior of solutions to the stochastic problem obeying Equations (
9) and (10), we adopt the field-theoretic RG method [
29]. This approach utilizes powerful quantum field theory techniques such as Feynman diagrams and RG resummation, among others. The starting point of RG analysis involves determining canonical dimensions and identifying all permissible structures (conveniently represented through Feynman diagrams) that exhibit superficial ultraviolet divergences [
40], as determined by their respective canonical dimensions. Typically, dynamical models like (
16) exhibit two independent scaling dimensions [
29,
53], implying that each field quantity
Q is assigned two independent canonical dimensions: the momentum dimension
and the frequency dimension
. These dimensions are individually determined based on standard normalization conditions and the requirement that each term in the action must be dimensionless (momentum and frequency separately) [
29,
40]. As in most models of critical dynamics, in the model under consideration, one can introduce the total canonical dimension
for a given quantity
Q, expressing invariance with respect to some “cumulative” transformation of coordinates and time. Since in linearized dynamic Equations (
9) and (10), the operations
and
∇ are included in the combination
common to all fields, the total canonical dimension
in our case is determined by the relation
By direct calculation from the action (
16), the canonical dimensions of all relevant parameters are computed. In fact, the obtained values reproduce the result found in [
27]. For the reader’s convenience, we also list it here in
Table 1. From this table, we infer that Model (
16) is logarithmic at
and there is just one coupling constant
.
A detailed analysis of all relevant UV-divergent parts and related technical details can be found in previous papers [
27,
28]. Here, it is only noted that in our case, in addition to the 1PI functions
and
of the general
A-model considered in [
28], an additional vertex
is added. The counterterms to it have the structure of the bare vertex
and define the renormalization constant
. Furthermore, as a direct consequence of Galilean invariance (
12) (the derivative
enters (
16) only in the form of
), the terms
and
do not renormalize [
27]. However, it is important to understand that we have no renormalization of the vertex
only for the trivial case
and for
. For the case
, this statement relies on the antisymmetry of the corresponding vertex
with
[
27], which obviously does not exist, e.g., for
. Thus, further, we present expressions for all quantities in the form of expansions by the parameter
A, not explicitly fixing its value but always keeping in mind the two cases
and
, where the latter corresponding to MHD is of primary interest to us in the context of RG analysis.
Thus, taking into account the comments made above, the model is multiplicatively renormalizable, and its action (
16) can be written in the form
where only renormalized parameters and fields appear. The parameters of action
are associated with the corresponding parameters of
by multiplicative UV renormalization transformations of the following form:
where the whole set of the renormalization constants of the fields and parameters is determined from the relations
At this step, we have introduced the so-called renormalization mass parameter
[
29,
40].
The renormalization constants
,
are calculated using the ultraviolet renormalization procedure, which can be vividly described as follows. To make sense of the diverging Feynman integrals, some intermediate regularization is introduced (in our case, it is the natural ultraviolet cutoff
). Simultaneously, the bare parameters of the model are replaced with their renormalized counterparts by introducing
-dependent counterterms into the action. Renormalizability, then, means that with a special choice of the dependence of
on
(the requirement to cancel the UV-divergent parts of the diagrams), the limit
exists, at least for all integrals necessary for calculating all observables. This limit is defined up to a finite number of constants that can be fixed by setting the values of observables and the normalization of fields. In statistical field theory, the so-called minimal subtraction scheme (MS) is considered the most convenient, where only the divergent parts of the diagrams are subtracted [
29,
40].
As noted in
Section 2, from a technical perspective, the natural cutoff regularization is not very convenient for multiloop calculations, so usually, another step is taken: transitioning to the limit
along with the introduction of a new type of regularization like dimensional, in which singularities are presented as poles by the regulator and the corresponding Green functions are understood as corresponding meromorphic functions. In our case, such a regularization parameter is
, previously introduced in the correlator (
13), which parameterizes the energy pumping into the system. In this setting, the renormalization constants of Model (
24) are simply the principal parts of the Laurent series in the parameter
, whose coefficients, being polynomials in
g, can also parametrically depend on the space dimension
d (however, in this work, we assume
), the gyrotropy parameter
, the parameter
A, and the renormalized magnetic Prandtl number
u. Thus, a unified representation for these constants takes the form
For convenience, we assign specific notation to the sets of coefficients for each constant: for , for , and for . In what follows, to simplify the notation, we do not explicitly indicate the dependence on these parameters in the arguments of the corresponding renormalization constants.
The constants
, and
are determined from the condition of UV finiteness for the corresponding vertex functions, each obviously associated with its own set of diagrams. The last fact proves useful when verifying results with existing special cases. We are already familiar with the one-loop calculation results for specific cases
[
27],
[
44], and
[
65]. Our further goal is to complete the full two-loop calculation for all renormalization constants. In this regard, we present here details of the corresponding calculation for the renormalization constant
.
Let us begin the description of the renormalization constants’ calculation process from the constant
. Here, the previously noted issue regarding each renormalization constant corresponding to its own subclass of Feynman diagrams immediately becomes apparent. As can be directly inferred from the action functional (
24), the constant
corresponds to the pure Navier–Stokes equation (without the inclusion of the magnetic field) and can be calculated from Feynman diagrams constructed solely by propagators and vertices containing the fields
and
[
27]. As a consequence, it cannot depend on parameters
A and
u, but it may depend on the gyrotropy parameter
. Directly, the constant
can be determined through the corresponding Dyson Equation [
40] of the following form
where
represents the renormalized self-energy operator defined as a sum containing all permissible 1PI Feynman diagrams. Hereinafter, all quantities where all variables are replaced by renormalized ones are marked with a subscript “
R”.
The coefficient
representing the one-loop contribution to
was calculated for the first time in the work [
41]. Subsequently, the corresponding two-loop contributions for both the nonhelical (
) and helical cases were detailed in [
42] and [
34], respectively. In addition, in ref. [
42], it was shown that the coefficient
is exactly related to the coefficient
through the relation
. The explicit expression for the coefficient
is rather long and we shall not rewrite it here explicitly (see [
34] for the details). Its value for the most physically important case
, along with the exact value of the coefficient
, is given in
Here, we introduce a typical geometric factor
, where
is the surface of the
d-dimensional unit sphere, and
is the Euler gamma function. At this point, it is important to note that the presence of a spiral contribution in the turbulent system not only does not destroy the stability of the Kolmogorov scaling regime determined by the coefficients indicated in (
29) but also has no impact on its properties at all [
34]. Additionally, it is worth mentioning that this parameter naturally arises when the magnetic field is included. From RG analysis of the three-dimensional MHD problem with a general correlator matrix, it is clear that
in our case plays a role similar to the gauge parameter
in QED [
27].
The RG constant
is calculated in the same vein as
. Direct inspection of Feynman diagrams reveals that only those for which the vertex
is not present contribute to
. This omitting in action (
24) is equivalent to moving to the passive admixture model [
28,
65]. The corresponding Dyson equation in the frequency–momentum representation for the renormalized 1PI Green function
takes a similar (
28) form
Note that from the last relation (
26), it follows that the product of fields in
is not renormalized:
. The coefficient
representing the one-loop contribution to
, is given by
Corresponding two-loop expressions for
coefficients can be found in [
28,
65] for the nonhelical term and the helical term, respectively. Since they are generally rather technical and not very illuminating, we shall not reproduce them here.
The renormalization constant
is up to now known in the perturbation theory only for the MHD problem in the leading one-loop approximation [
27]. For the two-loop calculation, the procedure follows similar steps to the previous cases. However, from a technical point of view, it poses a more challenging task. We determine RG constants
from the three-point 1PI Green function
. The needed perturbative relation can be formally expressed in the form
The subclass of Feynman diagrams that contribute to vertex
can be defined simply as an absolute complement to the two subclasses already discussed, defining
and
. Using (
32), the renormalization constant
is determined from the UV finiteness condition for
, which is more preferable to formulate in terms of a scalar quantity
Here, we take advantage of the fact that
has the bare vertex index structure. The limit
in Expression (
33) does exist, provided the IR regularization of the graphs has been properly taken care of. In terms of quantity
, Equation (
32) can be reorganized as follows:
where
is the loopless (tree) approximation,
denotes scalarized contributions from the one-loop diagrams,
represents the scalarized total contribution from the two-loop diagrams, etc. In diagram language, this equation then takes the form
where the diagrams depicted represent the one-loop contribution
, diagrams
are shown in
Figure 3, and
are their symmetry factors. The symmetry factors
,
, and other two-loop diagrams have symmetry factor equal to 2. Then, the UV finiteness condition for
may be cast in the form
where only the pole contributions to
are meant.
By direct calculation, it is not difficult to obtain the pole contribution to
in the following form:
from which we immediately obtain a one-loop contribution to
in the form
Formula (
38) generalizes the expression obtained in [
27] to the case
. From (
37), one can see that the contribution singular in
of the correction
vanishes at
. This is a general property due to the fact that the model with a passively advected vector field (case
) does not contain divergence in
[
33]. From this observation, it directly follows that
for the case
.
Calculating of the two-loop correction
turns out to be a significantly more complex task. Compared to
and
, not only does the number of Feynman diagrams increase (twenty-nine versus eight diagrams) but also their structure is considerably more complicated. The diagrams are computed in the frequency–momentum
representation, and, after analytic integration over frequencies, the remaining expressions can only be examined numerically in the same manner as in [
43]. The integrals remaining after frequency integration contain parameters such as the reciprocal Prandtl number
u, whose value in [
28] was determined to the two-loop precision. Note that within the framework of the
-expansion, for the second-order calculation, it suffices to compute the two-loop Feynman integrals at the value of
determined at the one-loop fixed point. This is so because the integrals themselves are already proportional to
, which is of order
. Finally, the general pole part of
can be formally expressed through the following structures:
where the coefficients
C,
, and
are functions of
d,
u, and
A and are to be further computed. The nonzero value of the coefficient
C arises specifically in cases where the diagram contains one of the two divergent subgraphs represented in Equation (
35). Its independence from the chirality parameter
stems from the model’s multiplicative renormalizability and the fact that the one-loop contribution does not depend on
. The values of the coefficients
C,
, and
for each specific diagram in
Figure 3 can be found in
Appendix A. In addition, to explain intermediate steps in the aforementioned described calculational procedure, we briefly elaborate in
Appendix B on the used methodology in the analysis of a particular two-loop Feynman diagram.
From Expression (
39), it follows in a straightforward manner that the two-loop contributions
and
to constant
have the form
The coefficient
C at the second-order pole is not of particular interest to us because it does not contribute to the RG analysis. Its value is only needed to verify calculations through cancellation with corresponding one-loop contributions. The corresponding expressions are given in
Appendix A. Therefore, we limit ourselves to presenting the value of the coefficient
for the most physically interesting case of three-dimensional MHD
The RG analysis for all 1PI Green functions
is based on the fundamental differential equation [
29,
40] expressing the independence of the renormalized functionals
from the finite renormalization transformation parameter
at the fixed model’s bare parameters
. The equation can be presented in the following form:
Here, we have slightly abused the notation and returned to the terminology of
Section 2 for 1PI functions. The gamma functions, namely,
,
, and
are commonly defined through the following relations:
,
, where
, i.e., the logarithmic derivative is taken at the constant bare parameters. Concurrently, the beta functions
and
expressing the change in the charges
g and
u under the renormalization group flow are defined as follows
The RG functions
and
are associated with the calculated two-loop constant
. The direct calculation yields two relations,
where the former one is a consequence of the equation
given in (
26).
In MHD, as in the theory of homogeneous and isotropic turbulence, the nontrivial asymptotic large-scale behavior of Green functions is governed by anomalous dimensions (exponents)
, which are determined by corresponding gamma functions
taken at an IR stable fixed point of the corresponding RG equation (
42). In its turn, coordinates of fixed points
in the plane of renormalized couplings are determined from the requirement that all beta functions simultaneously vanish:
and
. The stability of the found fixed point is determined by the positive definiteness of the real part of the matrix
,
in the vicinity of
.
In our model, there is one nontrivial stable fixed point – the so-called kinetic point [
29]. With two-loop accuracy, its coordinates are determined by the following relations:
and
, where the corresponding coefficients of
-expansions were calculated in general form in [
28]. In the case of the three-dimensional MHD, they take the following form:
An important difference between MHD and ordinary hydrodynamic turbulence is the presence of a nontrivial anomalous dimension (exponent)
of the magnetic field, calculated perturbatively. Let us recall that due to the absence of renormalization of the velocity field, the corresponding anomalous exponent of the velocity field is exactly equal to zero:
. The two-loop anomalous exponent
at the kinetic point has the form
which, along with the well-known exact (lacking corrections of order
or higher) values of the anomalous dimensions
and
[
29], completes the set of two-loop anomalous exponents in the model under consideration.
Along with the canonical dimensions
given by (
23), the large-scale behavior of Green functions (critical scaling in the language of RG) is also characterized by the corresponding critical dimensions
where
and
are the corresponding canonical dimensions listed in
Table 1. Therefore, the critical dimensions in MHD are determined by the following relations:
Note that the absence of
entails the fact that the second condition in (
49) includes only exact
. This precisely determines the Kolmogorov spectrum “
” for the statistics of the velocity field correlations at the “physical”
. In its turn, the presence of a nontrivial anomalous exponent
of the magnetic field can have a significant influence on the scaling behavior of the statistical correlations of the magnetic field.
In what follows, we finally restrict ourselves to considering only the case of three-dimensional MHD. Thus, using Expression (
47), one obtains the critical dimension of the magnetic field
For the physical setting where , the critical dimension manifests as . This dimension peaks at and and reaches a minimum of for . This observation underscores the significant influence of helicity on the scaling behavior within at least the regime governed by a kinetic fixed point.