Renormalization Group Analysis of Turbulent Hydrodynamics

Turbulent hydrodynamics is characterised by universal scaling properties of its structure functions. The basic framework for investigations of these functions has been set by Kolmogorov in 1941. His predictions for the scaling exponents, however, deviate from the numbers found in experiments and numerical simulations. It is a challenge for theoretical physics to derive these deviations on the basis of the Navier-Stokes equations. The renormalisation group is believed to be a very promising tool for the analysis of turbulent systems, but a derivation of the scaling properties of the structure functions has so far not been achieved. In this work, we recall the problems involved, present an approach in the framework of the exact renormalisation group to overcome them, and present first numerical results.


Motivation
Since the first work on statistical hydrodynamics by Kolmogorov [20], turbulence basically remained an unsolved problem of classical physics, even though the fundamentals seem to be simple enough -the Navier-Stokes equations just express the conservation of momentum of an incompressible fluid. Nevertheless, experimental results on the moments of velocity differences (the structure functions to be defined below) have yet to be understood on a theoretical basis. While Kolmogorov's assumptions of scale-independence fail, the intrinsic scale-dependence that is responsible for the intermittent exponents have not been deduced from first principles so far. A promising approach seems to be the Renormalization Group (RG), which aims to describe the dependence of the correlation functions of a given field theory on the scale on which the system is observed. Beginning with the works of Forster, Nelson and Stephen [11], numerous attempts have been made to apply the various formulations of the RG to turbulent hydrodynamics, but until today the observables proposed by Kolmogorov could not be deduced in accordance with experiment. For some work on the RG approach to turbulence, see [8,1]. The formulation of the RG most suitable for the study of turbulence is, in our opinion, the so-called Exact Renormalization Group (ERG), due to Wilson [32,33]. Although the distinction between "field theoretic" and "exact" RG is merely artificial, the latter explicitly involves the generating functional of the theory to be studied, which can be simplified and re-formulated to suit the analytic methods involved. It is especially helpful, as we shall see, for the analysis of a theory with constraints, like in our case the incompressibility condition. We will show how to apply the well-known Faddeev-Popov gauge-fixing method to work out the functional integral measure. In this step we differ from previous work, which in one way or another omitted the incompressibility condition and/or the pressure term of the Navier-Stokes equations. The RG-flow, as we shall discuss in detail, can be understood as the continuous way of calculating all Feynman graphs of the theory. Keeping this in mind, we established a numerical algorithm that simulates the RG-flow by integrating out the corresponding graphs. We therefore take advantage of the freedom of choice of a cutoff-function for the propagator. We end up with lengthy, but straightforward rate equations, that can be iterated quickly and up to a high number of involved field operators. We will also show that the predictions of Kolmogorov can be found as the trivial scaling solutions of this theory, and that non-trivial structures in coupling space exist. We have tested the algorithm on well-known theories; the identification of intermittent exponents in turbulence, however, has not yet been accomplished.

Turbulent Preliminaries
We shall introduce the basics that are needed in this work, and sum up Kolmogorov's predictions from 1941 (K41) as far as we aim to falsify them, rather than to give a full account of Kolmogorov's theory. Reviews can be found for example in [25,13,29,24].

Navier-Stokes Equations
We shall mainly work with the full Navier-Stokes Equations (NSE) given by where v denotes the D-dimensional velocity field (in Navier-Stokes turbulence, D is either 2 or 3), ν the kinematic viscosity, p the scalar pressure field and ρ the density of the fluid. Eq. (1) expresses the conservation of momentum in an infinitesimal volume element of the fluid. As these are three equations for four degrees of freedom, it is clear that we need an additional constraint, determining the pressure field. In incompressible turbulence, this constraint simply states that the velocity field is divergencefree, Starting from any given initial condition, the velocity field obeying Eq. (1) has to develop into a constant field, as the dissipative term converts more and more energy into heat.
If we aim to understand fully developed turbulence, statistically homogeneous in space and time, and statistically isotropic in space, we need a mechanism to insert energy into the system, so that an equilibrium flow can develop. The standard way of providing this is to add a stochastic force (stirring force) f α to Eq. (1) that is long-range correlated: The idea is to bring energy into the flow on large scales, let large structures decay freely into smaller ones until the energy is finally dissipated into heat (Richardson-cascade). We model the stochastic force to be Gaussian distributed, with δ-correlation in time, and a long-range correlation function in space: where ǫ is the local energy dissipation rate (not to be confused with the parameter of the ǫ-expansion). ∇ −2 de-notes the fundamental solution of the Laplacian, e.g. in 3 dimensions: Different forms have been tried for F , though in the context of the NSE it is widely believed that the form of the stochastic force does not influence the statistical characteristics of turbulence. On the other hand, it should not be forgotten that e.g. in Burgers turbulence, the intermittent exponents (to be defined below) clearly depend on the choice of the stirring. Eq. (2) is sufficient to eliminate the pressure term, as can be seen when deriving the solenoidal NSE, as we shall do in the following. In a first step, let us operate onto Eq. (3) with a divergence operator. Then the first and the third terms drop out, as the field is divergence free: We then invert the Laplacian, arriving at which is the above mentioned condition for the pressure field. One might ask whether the inversion of the Laplacian leads to a unique solution for p. Now, two solutions might at best differ by a harmonic function, which is either constant or growing without limits. The second option is not physical, the first one not relevant as we are only working with pressure differences.
To obtain the solenoidal NSE, replace the "solved" pressure field into Eq. (3): Observe that the operator P αβ = δ αβ − ∂α∂ β ∇ 2 is precisely the transverse projector known from electrodynamics; a feature that we will use later on. From now on we shall investigate the solenoidal NSE. It is important to keep in mind that these are only equivalent to the full NSE as long as incompressibility is ensured. Also observe that Eq. (10) is non-local, as it involves the inverse of the Laplacian operator, the integral kernel of which is of the form (6).

Gauge-Invariance
As a preparation for a Faddeev-Popov method, that we will need in the following section, we shall now rewrite the formula in a gauge-invariant way. As the transverse operator P projects a field onto its incompressible parts, and knowing that the fields we are interested in are transverse a priori, it is easy to see that so that we are free to replace v by P v in Eq. (3): Now, even if the resulting equation looks more complicated, it is invariant under the same gauge transformation known in electrodynamics, Constraint (2) is still required, but it now acts as a gauge fixing term.

Structure Functions and Intermittent Exponents
In 1941, Kolmogorov published his famous article, introducing a statistical framework for turbulent hydrodynamics [20]. As the theory is Galilean invariant, he proposed that the observables should be functions of velocity differences, and more specific, he considered the so-called velocity increment the difference of the velocities at two points separated by a vector x, projected onto the unit vector in x-direction.
Suitable observables appear to be the structure functions of order p. They are defined as the p-th moment of the distribution of the absolute value of the velocity increment: The average is taken over all points r of a realization of the turbulent flow. In the case of homogeneous turbulence, we suppose this to be equivalent to an average over all histories v(x, t), which we will perform in the functional integral.
Kolmogorov proposed the existence of a smallest length scale λ, also called the "dissipation scale", below which physics is no longer dominated by turbulence, but by dissipation. Dimensional analysis leads to where ǫ is the (constant) dissipation rate. Assuming that the turbulent cascade of decaying vortices happens on scales much larger than λ, we might hope that observables don't depend on it and are thus self-similar, which means power-law functions of the scale: Again, Kolmogorov deduced by dimensional analysis that It has long been stressed that the fundamental assumption, namely the independence of the smallest scale λ, is by no means natural, and is in general not fulfilled in critical systems. This could lead to a scale dependent dissipation rate (or viscosity) and the breakdown of scaling law (17). Even though general agreement on this point seems to be common, the scale dependence could not yet be deduced. In case that a typical (macroscopic) length scale L can be identified in the system, the Reynolds number is defined as L might be the radius of an obstacle of the flow, or, more useful for our purposes, the correlation length of a choice of the stochastic force. Eq. (19) coincides with the more common definition if the typical velocity U is defined to be

Burgers Equation
As we shall refer to it later, we mention Burgers' equation as a simpler model for turbulence, which can be seen as fully compressible hydrodynamics. Further constraints are not needed as a pressure term is missing. Burgers turbulence, or Burgulence, is investigated in one to three space dimensions (see e.g. [3]). It is local, and its solutions typically contain structures that resemble the fundamental solution of the Hopf-eq.
as shown in Fig. (1). These are regularized by the viscosity term to a finite step width (similar to the behaviour of the KPZ-equation, see [21]). v(x; t = const.) x Shock structure in Burgers' con guration Monte-Carlo simulation [10]. The graph shows v(x, t) as a function of x in periodic boundary conditions at constant t.
As published in a different article [10], we investigated a generating functional of Burgers equation by means of Monte-Carlo methods. These investigations have shown that the generating functional formalism is suitable for the study of turbulent systems, even if almost-singular structures are involved.

Exponents for Burgers Turbulence
In general, the intermittent exponents of burgulence strongly depend on the form of the stirring force; for an overview see e.g. [14]. In the canonical case, that of freely decaying burgulence in 1+1 dimensions, analytic calculations lead to ξ p = min (p, 1) .
This scaling can be interpreted as the signature of a bifractal decay channel.

Generating Functional
The basis of any analysis using the Exact Renormalization Group (ERG), first introduced by Wilson, or many other field theoretical methods, is the generating functional. For turbulence several attempts have been made to define it as a functional integral, all in one way or another omitting the condition of incompressibility. We shall begin with the derivation of the Martin-Siggia-Rose functional [23] for the solenoidal NSE, and than show how to respect condition (2) by means of a Faddeev-Popov procedure.

Fine-Grained Distribution
Our starting point is the so-called fine-grained probability distribution for the field v, obtained by counting all possible solutions to the NSE.
where N is defined in Eq. (12) 1 , and we defined the abbreviation A few remarks are in order: -Of course, N −1 is not to be understood as the (not necessarily existing) inverse of an operator, but as a multi-valued operator counting any solution v for a given realization of the random force f . Observe that we are working with a functional δ-function, meaning that we are searching for histories v(x, t) that solve the NSE for all x and t, rather than a realization v(x, t 1 ) at a given time t 1 , depending on some initial condition. -The average · is to be understood as an average over all realizations of f , replacing the spatial average in (15). It is the basic assumption of this work that for homogenous and isotropic turbulence, these averages are interchangeable; we are supported by our results for Burgulence published in [10].
Making the average over all realizations of the stochastic force explicit yields Now the argument of the δ-function is multiplied by N, which leads to a functional determinant that we shall discuss in the next paragraph: The δ-function is now written by means of a functional Fourier-transformation, introducing a new formal (bookkeeping) field u: If we formally define the action S 1 by where the determinant still has to be evaluated, we can identify the elements for our Feynman-rules. For the time being, we want to focus the attention on two parts, which together lead to the famous θ(0)-problem: Notice that the corresponding bare two-point-function, also called response function, is proportional to the Green's function of the diffusion equation, applied to transverse fields: We have to choose the retarded Green's function to ensure causality of the theory, so uvv-vertex: In a graphical expansion, it is required to calculate the following loop: As the retarded Green's function of (not only) the diffusion equation is proportional to θ(t 2 − t 1 ), this loop is proportional to the seemingly ambiguous quantity θ(0). We will show in a following paragraph how this is solved by a choice of discretization of the functional integral.
To illustrate what the above action actually implies, we integrate out the non-physical fields f and u for a moment.
Observing that first the integration over f , and then the one over u are Gaussian, we can simplify Z to: with We see that field configurations not solving the NSE are admitted, but suppressed by a Gaussian weight that we can control numerically, and that can be related to a physical constant by Eq. (5). This expression is already a generating functional, as we are able to extract all correlation functions using functional derivatives. Before we can start to simulate RGflows, we still have to bring the determinant into a suitable form.

Functional Determinant
A straightforward way of writing the functional determinant is by using ghost fields, where we dropped a field-independent term and defined I to be the non-linear part of N , We thus end up with the functional where It is important to keep in mind that the ghost fields, though anticommuting with each other, don't pose any particular problems in the numerical procedure. In principle, it is very much feasible to take them into account, and we performed several RG-runs while keeping track of the ghost terms. The determinant has a simple graphical interpretation, as we see directly that it exactly cancels out the u-v-loop shown in Fig. (2): From eqs. (31) and (36), we see that the ψ * ψ-and the uv-propagator are identical. The reader can easily check that ψ * , δN δv ψ leads to two terms similar to uN[v], but with u replaced by ψ * and one v-fields replaced by ψ. When this vertex is closed to a loop by means of a ψ * ψ-propagator, this is numerically identical to graph (2). This greatly simplifies matters for us, as our numerical program sorts and calculates contributions to the RGflow according to their graphical representation. Rather than simulating two additional fields, and calculating all the graphs, we are thus allowed to drop a certain class of graphs. The cancellation of certain averages can be proven even non-perturbatively, using the BRS-invariance of the action.

BRS-Invariance
If the functional determinant is expressed in terms of ghost fields, this yields an extra symmetry, also called BRS-invariance. We observe that the action (39) is indeed invariant under the following infinitesimal simultaneous changes: that for obvious reason have also been called "half a supersymmetry". From the Ward identities of this symmetry, the desired result follows on a non-perturbative level: Both sides of this equation can be interpreted as the sum of the corresponding graphs as explained in the previous section. For details we refer the reader to the explict proof in [27].

The θ(0)-Problem
We have seen in (3.1) the typical appearance of the seemingly arbitrary number θ(0) in the retarded Green's function of the diffusion equation. It is of course linked to the Itō-Stratanovich-dilemma; Zinn-Justin [35] shows explicitly how it is connected to the problem of operator ordering. To clear the mist, we want to show in this section how θ(0) is fixed by the choice of discretization of the functional integral, namely the time derivative.
In any case the θ-function is required to obey We further require the fundamental theorem of calculus to be valid, so We now have to define a differentiation. We distinguish different choices by a parameter α, suggesting with −1 ≤ α ≤ 1. Notice that α = 0 corresponds to the symmetric, and α = −1 to the pure backward derivative. We conclude (49) The term involving the second θ-function does not contribute since t ≤ 0 in the integral. The first θ-function only gives a contribution for We thus find Thus θ(0) is fixed by the choice of the time derivative. Notice that for the symmetric (Stratanovich) derivative, α = 0, we find θ(0) = 1 2 , while in the pure backward (Itō) case we get θ(0) = 0, as expected.

Incompressibility Condition
Before we demonstrate how to derive a generating functional respecting condition (2) for the full Navier-Stokes equations, let us briefly discuss three common positions concerning incompressibility: -Incompressibility doesn't matter. Though not stated explicitly, this appears to be the point of view taken in those works that neglect the additional condition (2) without any further comment. This is more than questionable, because there is a whole class of equations that only differ by the compressibility condition, and give different statistics. An obvious example is Burgers' equation (22), which models fully compressible fluids and shows bifractal scaling of the structure functions. -Incompressibility is enforced by a choice of boundary conditions. This argument goes as follows: In the solenoidal form, it can be clearly seen that an incompressible flow stays incompressible even without enforcing it by a particular equation, as the random force term on the one hand, and the former pressure term on the other hand lead to incompressible contributions to the flow. On the other hand, all compressible parts of a given flow tend to die out due to the dissipation term. This leads to the observation that condition (2) is always ensured if we impose an incompressible flow for t → −∞. This observation is indeed useful for some calculations, e.g. direct numerical simulations, as they have a definite boundary at t = 0. In our case, we aim for statistics over whole histories of the fields, for all times -and in negative t-direction, as the argument itself states, any compressible part of the flow is going to be amplified. Now the slightly compressible flows and the incompressible ones lie dense to each other in functional space, so that in any numerical application we would lose control of the boundary conditions completely. -Incompressibility is ensured by some standard integration measure. What is meant here is that incompressibility is considered to be included in an integration measure Dv inc (we are going to proceed in the same way in the next paragraph), but then the measure is treated as a standard Lebesgue-like measure. This is certainly not legitimate; e.g. an integral like looks like a Gaussian integral, but will in general not be so.
We thus conclude that incompressibility should well be taken care of. It might be surprising that the way to do this is rather straightforward and well known for gauge invariant field theories.

"Strong" attempt using δ-function
We implement condition (2) through a second functional Fourier transformation of a δ-function, The advantage of this is that the generating functional is exact, apart from the definitions of some constants. The major drawback is that we will have to apply the Renormalization Group in a discretized way, which requires an approximation of the generating functional. Now a constraint like (53) can, after approximation, lead to a field theory that has no solutions at all, as it is driven from the manifold of solutions further away by every iteration of the RG-process. From a technical point of view, the additional field makes simulation times longer and most calculations more elaborate, but does not represent qualitative problem.

Faddeev-Popov Method
A proper way to treat the compressibility condition is by means of the Faddeev-Popov method. Let us first generalize our considerations to the case ∂ α v α = c. We multiply the functional by 1: where U denotes the gauge transformation as defined in (13). The determinant can be written explicitly as and can be dropped as it is independent of any field. We now multiply the functional by another (field-independent) factor, namely Now our original action is gauge-invariant, meaning that the generating functional does not depend on the value of c. We might choose a particular one, or as well integrate over all possible c: The reason to integrate over all c is to evaluate the δfunction, finally arriving at The gauge fixing term appears as a Gaussian distribution, and only in the limit κ → 0 incompressibility is enforced strictly in the integrand. This method overcomes the disadvantages of the δ-function method, as deviations from solutions aren't completely forbidden, but only suppressed by a Gaussian weight. The limit can be taken when results have been obtained. We shall from now on work with the functional (62), but formulated in wavenumber space, which is: We discarded convention (26) because in the next section we shall directly manipulate some of the terms that would otherwise be hidden by the notation.

Galilean Invariance
A remark concerning Galilean invariance, as analyzed by Hochberg and Berera [5], is here in order. The path inte- In any attempt to calculate averages or correlation functions of the field itself, this might be a problem that could be overcome by another application of the Faddeev-Popovmethod. In our case, however, it is of no practical con-sequence, as we are aiming for the averages of velocity differences, which are gauge invariant.

Non-local Interactions
The action we derived still contains interactions of the type which are non-local in physical space, but of a very simple form in wavenumber space. Independent of the space where the theory is formulated, we need to rewrite (67) in a local way: In physical space, non-local interactions are at best cumbersome; in wavenumber space, we are going to sort the terms of the action according to their power of momenta, so we try to avoid 1 p 2 -interactions. We shall proceed in two steps: we will first re-define nonphysical fields, and then introduce new fields to transport the non-locality of interactions. We will end up with a lengthy, but local action that suits our needs for further analysis.
To make the non-local nature of the interactions more obvious, the functions within this paragraph are functions on physical space.

Field Redefinition
If we redefine the non-physical fields in the following way: and introduce we can rewrite the action as The functional determinant of these transformations is field independent and can thus be omitted.

Introducing new fields
Some non-local terms of the general form for some fields K and L still survive. We are able to solve this problem by introducing new fields. These can be understood as transporting fields that "carry" the non-local interaction from one place to another, thus replacing it by two local interactions and a propagator. Formally this means inserting a Gaussian integral of the type: where we defineM to be a new field. This leads to a new kinetic term 1 2 (M , ∇ 2M ) in the action, which is independent of all physical fields. We now shift the variables aŝ The point is that So we get rid of the original, non-local interaction (75) by replacing it by a new kinetic term for M and new interactions. Two of them are still non-local, but of the much simpler type for example. We can deal with them by the same method to get a local action finally. We add again a Gaussian integral for a new field, sayN , and definê from which we end up with The constant λ is needed so that the new fields get a definite dimension. This procedure might appear unusual, but is in fact quite common. The reader may think of electrodynamics: When a light is shining and is seen from a distance, this might be called a non-local interaction. In theory, a local electromagnetic field is introduced, that propagates from the source to your eye. Now, the equations of motion for the electromagnetic field turn out to be fairly simple, so that we are used to think of the field as being fundamental. In the present case, think of the fields N and M as merely helpful auxiliary quantities, though it might in other cases, like Burgulence, be advisable to model structures as new fields explicitly, and observe their interactions in an effective field theory.
Applying the method discussed above to action (63), we arrive at the local action with Here, φ 1 , φ 2 , φ 3 denote the new scalar fields. The terms are already sorted by their order of derivatives. This is the result for the action S. Depending on how the determinant (35) is expressed, other non-local interactions may have to be rewritten in the same way.

Discussion
In this section we have shown how to transform non-local into local interactions: either by redefinition of unphysical fields, or by introduction of intermediate propagators.
A major drawback will be that we have to approximate this action to account for RG-transformations, and a common way is the derivative expansion. Now, due to our "localization" of the action, the original terms have been mixed concerning the order of derivatives. Moreover, the number of derivatives has increased for most interactions, which means that we would have to expand the RG-flow to a high order in the derivative expansion. Also, the number of fields involved makes it a cumbersome work to do, even to the lowest orders. Nevertheless, this expansion is feasible to any order in derivatives, as shown in [18] for the first two orders. The results are rather lengthy rate equations, which will not be elaborated on here.

Renormalization Group Considerations
The renormalization group has its origins in the work of Gell-Mann and Low in 1954 [15] 2 . The first works concerning the Wilsonian Exact Renormalization Group (ERG) have been published in the early 70s [32,33], based on Kadanoff's block-spin picture [19]. It is surprising that some very basic questions, e.g. concerning the renormalization step and the anomalous dimension, are still being discussed. We will consider this point in detail in paragraph 4.4, as it seems to be in order (at least to us) to address these matters, especially the anomalous dimension and the graphical representation of the flow.
In this section, we remind the reader of the origins of the ERG, and will try to give a basic understanding of the flow equations. We will not repeat the derivation of the equations, as this can be found in a number of articles, but quite a bit of the graphical representation of the different terms, as it will lay the foundations for our numerical investigations that closely follow the loop expansion.

Form of the Action
We are looking for a RG-flow of a given theory defined by its generating functional Z. To be definite, let us work with the theory of a vector field v i , and write Z in the following way: The action depends on two momentum scales Λ and Λ 0 . By Λ 0 we denote the scale on which we impose the initial renormalization condition -e.g. the value of the fourpoint-function is fixed to a certain value λ 4 if all external momenta equal Λ 0 . The term S 0 might look uncommon, but is necessary to pick up terms nonlinear in J that will be generated by the RG-flow. As initial condition, we set S 0 [Λ = Λ 0 ] = 0.
Starting from a renormalized action on scale Λ 0 , the flow is going to generate the renormalized action on all lower scaled Λ, which is the second momentum scale involved. From the RG-perspective, S[Λ = Λ 0 ] plays the role of the initial condition of the flow. The flow equations depend on the choice of the kinetic action, so we will define the kinetic term to be where we define C is the cutoff-function, which has the following properties: Though it is by no means necessary, one usually assumes that C −1 is monotonous, and that it is a smooth approximation of the step function, thus suppressing degrees of freedom on scales bigger than Λ, while not effectively altering those on scales below. The last equation (92) is ambiguous, but we define a value for C −1 (1), so that the role of Λ becomes definite. We will say that the degrees of freedom that are suppressed are "integrated out", as this part of the involved integrals can be interpreted as already being performed. Apart from the properties (90-92), we are free to define C in any way we like. It follows that not even (88) is enforced; other definitions of the propagator have been tried. In practice, some propagators will lead to simpler numerical calculations than others. A very special choice of C is the sharp cutoff C s : which leads to the Wegner-Houghton equation and will be treated separately.
In a similar matter, we define the following kinetic terms for anti-commuting Grassmann variables: where P −1 Ψµν is an antisymmetric matrix in the indices µ and ν. For completeness, we already mention here that we will expand the interaction part of the action, S int , in powers of the fields to illustrate some examples, where we will call a k-vertex. The derivation of the flow equations does not depend on this expansion; but it is useful in some definite calculations.

Integrating out degrees of freedom / Lowering the Cutoff
A simple way to derive the RGE is to calculate the effect of a change of the cutoff on an action, keeping in mind that both the generating functional and the correlation functions may not change. 3 In our case, we will lower the cutoff by lowering Λ, leaving Λ 0 as a unit of measurement unchanged. Applied to the vector theory, for example, we arrive at the following equation for the interaction term of the action: where we dropped a field-independent term. Taking the kinetic term into account, we find the simple equatioṅ The term 1 2 p δS δvjṖ vji δS δvi will from now on be called the link-term of the flow-equation, while we will call − 1 2 p δ δvjṖ vji δS δvi the loop-term. These names will be justified in the following subsection. In complete analogy, equations for theories involving Grassmann variables ψ * and ψ with propagator (95) can be derived: In case of anti-commuting fields it is important to keep track of all extra signs that arise. These equations describe the lowering of the cutoff, or integrating out of degrees of freedom. Before we proceed, let us remark on the graphical interpretation of the RGequations.

Graphical Representation
Let us begin with the interpretation of the link-term. For the time being, we assume that it is applied to a part of the interaction term of the form (97), a vertex with n 1 + 1 attached lines. Then the functional derivative of this gives us a vertex with n 1 lines; the missing line is linked by the part of the propagator that is integrated out,Ṗ , to a second, similar vertex with, say, n 2 + 1 lines. The graphical result is shown in Fig. 3. Observe that the functional derivatives automatically lead to the correct symmetry factor of the graph.  This graph gives a contribution to the (n 1 + n 2 )-vertex, proportional to (105) The two δ-functions will pose problems in the numerical process, as they will have to be approximated in order to be suitable for a derivative expansion. We can eliminate one δ-function by realizing that which gives us at least the overall conservation of momentum that we need. The loop-term is equally easy to understand: From a vertex with n + 2 attached lines, two are joined by a prop-agatorṖ (Fig. 4). Again, the symmetry factor is given correctly.  This graph gives a contribution to the (n)-vertex; The overall δ-function does not pose any problems in this case. So far, we explained the effect of the flow equation as only S int is concerned. Let us now investigate the contributions of the kinetic term. The loop-term generated from the kinetic term (Fig. 5) is trivial , as it is field-independent and can be dropped. Let us consider the terms arising when one field derivative in the link-term acts on the interaction, and the other one on the kinetic term; this one is compensated by another term in the RG-eq.: The remaining term to be considered is We can use1 and this is, as defined in (88), the change in the kinetic term. Let us summarize: We have seen that the RG-flow can be expressed graphically. Iteratively, we calculate the contributions from all graphs with one propagatorṖ , and all other propagators P , that are inner propagators which have been generated in RG-steps before. This is simply an application of the product rule: If we formally define G[f (p)] to be the sum of all possible Feynman-graphs of the theory with inner propagators f (p), we can even write down a formal solution to the Wilson equation: In the limit Λ 0 → ∞ this becomes This reveals the meaning of the changes to the action: The RG-flow sums up the part of the propagator that is cut off iteratively. Notice that we seem to subtract all graphs -this is because we are working with e −S rather than e S . Of course the derivation and graphical interpretation of the flow-equation does not apply directly to the case of a sharp cutoff as defined in (93, 94). As our numerical approach favours the Wegner-Houghton equation, we discuss it in the following paragraph.

The Wegner-Houghton Equation
For numerical purposes, it is easiest to work with the sharp cutoff function C s defined in (93, 94). This changes the form of the flow equation drastically. Again, we will not present a derivation here as it is found in the original literature [31], but present a short overview. We start from dividing degrees of freedom into a high-momentum part that is to be integrated out, and a low-momentum part that is kept. Expressing the integrated part of the functional as a change ∆S to the action, we find where the prime denotes integration over the momentum shell between Λ − dΛ and Λ. We now expand the action two second order in the fields and get where we completed the square. Integrating out the field v in the shell of momentum, which is assumed to be of the thickness ∆Λ Λ ≪ 1, we get and with the aid of (det A) α = exp{α ln det A} = exp{αTr ln A} (122) we can finally write down the Wegner-Houghton equatioṅ In the derivation it is used that it is sufficient to work in one-loop-order, as higher order contributions are also of higher order in dΛ. A proof of this is found in the original work [31].
Notice that (123) seems to differ from the Wilson equation (103) by an overall-sign; but this is explained asṖ is negative.
Eq. (123) is especially convenient for numerical applications, as the contributions to the integrals can be calculated explicitly. The logarithm looks problematic at first glance, but we shall demonstrate in the next paragraph that it has a very simple graphical interpretation, and is thus favourable for our graphically based program.

Graphical Representation
Once again we start with the interpretation of the linkterm, which for the Wegner-Houghton equation does not involve the bare propagator alone, but the quantity , which is more than just the inverse of the kinetic term. Using the geometric series, we can write This is the sum of all graphs with i vertices, linked into a line by propagators. The first and the last vertex in the line are also attached to propagators, which link them to terms δS δv . Again, let us assume first that these act on the interaction part of the action, thus the chain described above is linked to other vertices. We therefore find the graphical representation Fig. 6.  In a similar way, we derive the graphical representation of the loop-term. We re-write the logarithm as: The first term is field-independent and is dropped. The second logarithm is expanded as a Taylor-series, reading The corresponding graph is shown in Fig. 7. It is similar to the link-term before, but closed to a loop by the trace. The factor 1 n compensates the rotational symmetry of the graph.   As before, we still have to sort out the link-terms involving the kinetic action. Let us start with an example. P P . . . r 1 r 2 r n P −1 P −1 q 1 q 2 Fig. 8: Graph of the Wegner-Houghton flow, linking a vertex to two outer propagators in this case.
The graph in Fig. 8 is obviously one of those that arise from the link-term; the reader may focus his attention to one of the P P −1 legs. Integration is again over the momentum shell, so P P −1 = 1. The result looks like the vertex, but with the difference that we know that fields q 1 and q 2 are depending only on momenta less than Λ − dΛ.
We can thus interpret these terms as integrating out the momenta on remaining fields. Integrating out more outer fields at the same time would again be of higher order in dΛ, and is thus omitted. The change in the kinetic term itself is again simple, and not even a sign problem arises as in the Wilson-case. The corresponding graph is and as P −1 P P −1 = P −1 , this is exactly Again, this is precisely the change of the kinetic action, as expected from Eq. (88). As in the case of the Wilson equation, we are now able to give a formal solution to the WHE. The final result reads

Kinetic Action -Renormalization and Rescaling
The field renormalization step is completely artificial. Let us clearly discuss the renormalization step, as it seems to have lead to confusion. For definiteness, we will demonstrate the concept using φ 4 -theory in D-dimensions; for other theories the procedure works in exactly the same way. Let us emphasize that this step is not unique to the ERG, but also found in perturbative renormalization.

Field Strength Renormalization
We started our integration step with the kinetic term We thus defined S kin to be the only term quadratic in fields and momenta in the limit p → 0. After the integration step (lowering the cutoff from Λ 0 to Λ), new terms that, according to our definition, should belong to the kinetic term, can be found within the interaction part of the action. Such a term has to be counted as kinetic term by hand. In practice, the first contribution to the kinetic term arises in the second step of the RG-flow as it is of two-loop order. The simplest graph contributing to field strength renormalization is the so-called sunset graph, Fig. 10. In our case, a graph analogous to Fig. 10 is to be computed by our numerical approach in an iterative way later on, summing up contributions from every infinitesimal integration. The result depends on the used renormalization scheme; by means of a Taylor-expansion, we are always able to identify the contribution to the kinetic action; we will call this This term is transferred from the interaction term to the kinetic action by hand, which then becomes As our renormalization condition for the field strength, we will choose that the coefficient of the kinetic action in the limit p → 0 is equal to 1 2 . According to the definition of the cutoff-properties, we introduce the field-strength renormalization factor Z in a way that compensates for the new term in S ′ kin . If we write for the original action (133) we conclude that Z transforms as (The sign changes as in the integration step, we actually lower Λ.) The initial condition has to be so that (133) is fulfilled. We can integrate (137) easily to Z is now absorbed into the fields in the following way, which explains the name of field strength renormalization: We will use this in the next paragraph to determine the scaling of the field. This is in complete accordance with renormalization conditions met in perturbative renormalization, see for example [28,7]. The kinetic term now reads which, expressed in renormalized fields, is exactly of the same form as (133). The point stressed by Golner [16] and Bervillier [6] is that we have to do the rescaling step very carefully, to account for the renormalization step consistently.

Rescaling
The last step in the renormalization group process is the rescaling of the momenta, and the functions thereof. We define new momentap in the following way: so we replace We can directly see that this replacement also changes the cutoff function in the expected way: so we are ready to identify the new with the old cutoff. From the renormalization step, it is now easy to deduce the scaling of a field. Beginning with the original kinetic term (133), we conclude that, as S kin does not scale at all, the rescaled fieldφ has to be scaled as We call the canonical dimension D φ,can = D−2 2 . Eq. (140) tells us the anomalous exponent: We will now see that due to the effect of the rescaling of the cutoff-function, the renormalized kinetic action indeed does not scale. As the cutoff is a function of the ratio p 2 Λ 2 , a change in p has the inverse effect as the same change in Λ. By this, we reverse the effect of the integration step. As a part of this, terms are re-distributed back to the interaction of the theory. Now scaling back the kinetic term: The canonical scaling of the fields is compensated by the integration measure, and the anomalous scaling by the renormalization step. We thus trivially find again formally identic to (133) as a function of the rescaled quantities, as desired. We will drop the tildes and primes, formally getting back to (133).

The Interaction Terms
As we have explicitly discussed all steps when applied to the kinetic action, we now just have to apply them to the interaction terms. Let us, as an example, have a look at a four-field-interaction

Renormalization
Once we change φ to the renormalized field φ ′ , without changing the interaction term, we have to renormalize the coupling as follows: We need the change in λ for an infinitesimal integration step, which isλ or, for a general interaction S int with any number of vertices,Ṡ int, Ren = as the operator φ δ δφ counts the number of fields in a vertex.

Rescaling
This step is now straightforward. Let us sum up the contributions: -Integral: For each integration measure, we get a factor D, and there is one integration measure less than there are fields (because of the δ-function), so we get a contributioṅ (we are still lowering the cutoff). -Momentum: The vertex could and will depend explicitly on the momentum, so we introduce another operator φ(p)p ∂ ∂p ′ δ δφ(p) that counts the powers of momenta in each vertex. The prime at the derivative indicates that it is not acting upon the momentum conserving δ-function. We get the contributioṅ -Fields: As derived above, each field brings a contribution proportional to − D+2−η 2 , so in total we finḋ -Renormalized Coupling: Not to be forgotten, we renormalized the couplings according to (150), so it scales itself anomalously, exactly compensating the anomalous scaling of the fields: So, to sum up all contributions, we finally get the rescaling terṁ (155)

The RG-Equation
For the interaction term, we thus find the flow equatioṅ This equation clearly depends on the choice of propagator (89). To give an example: If we had defined the inverse propagator to bẽ as often found in literature, we would have ended up with the same equation as Bervillier [6] or Golner [16], namelẏ which in turn is completely equivalent to Wilson's equation.
We want to point out that both equations (156) and (158) are in a way correct, even though they seem to differ by a sign. The equivalence is obscured by a different choice of propagator functions, taking advantage of the reparametrization-invariance of the equation. Propagator (157) indeed has some advantages, as in principle other values for k are also possible, and simplify the implementation of K41, as we will see. On the other hand, the derivative of (157) is a bit uglier, and especially the derivative expansion gets mixed up.
We can now aim for the final assembly of the RG-equation, written down for vectorial and Grassmannian fields: This equation is general enough for all our applications, including different ways of considering the functional determinant Eq. (35).

Derivative expansion
Clearly, any action has to be approximated to be of numerical use, as it represents infinitely many degrees of freedom. A common way of approximation is the derivative expansion, see e.g. [17] and [26]. In principle, applied to the scalar theory, it amounts to expanding the action in powers of derivatives: in contrast to an expansion in powers of fields, which can be seen as expansion around a weak field. As a special case, in the Local Potential Approximation (LPA) the action is reduced to an interaction term depending only locally on the field φ(x) (and not on its derivatives) and a kinetic term whose coefficient Z is held constant throughout the flow: If we apply the loop-and link-terms to the action expanded in powers of fields, we are led to rate equations for the coefficients. Let us, as an example, apply the Local Potential Approximation (LPA) to Eqs. (105) and (107), expanded in powers of fields. Starting with Eq. (105) for the link-term, we see that the non-trivial part of graph 3 is proportional tȯ λ n+m ∝ Ṗ p 2 Λ 2 λ n+1 (p 1 , . . . , p n , p)λ m+1 (q 1 , . . . , q m , p) ×δ(p 1 + p 2 + . . . + p n + p)δ(q 1 + q 2 + . . . + q m − p)d D p.
As explained above, we transform one of the δ-functions into an overall momentum conservation, and use the other to evaluate the integral: In the LPA, we approximate the couplings to be momentum-independent, and develop the cutoff to zeroth order in the momenta, so: A difficulty is the momentum-dependence of the factoṙ P which we need to expand, according to the derivative expansion. The result obviously depends on the choice of the cutoff; if we apply it to the LPA, we can summarize the result into the constant If the cutoff is an approximation of the step function, we would expect the limit to converge, andP 1 = 0. This clearly is not an option, as it would suppress the nontrivial character of the RG-flow.
On the other hand, the loop-equation (107) leads in the LPA to Again, the integral depends on the choice of the cutoff; for the LPA we will write Rather than to specify a cutoff, in the LPA it is sufficient to define the constantsP 1 andP 0 . In higher orders of the derivative expansion, additional information concerning the cutoff will be required.
In the case of a vector theory in three dimensions, the situation is not that simple, as products of the type v i v i or any contraction with other three-component fields will be present. We need to keep track of this to calculate the contributions to a renormalization group flow, so we propose to expand the terms of the action in powers of fields and momenta in the following way: From now on, we will work with the coefficients In first order, the terms of the derivative expansion are even more complicated, as we also have to keep track of terms like p i v i (q).
As the overall number of momenta is fixed for each term, and the action itself is a scalar, we get the following possible values for the indices of V : and for Z equivalently. For completeness, we also have to keep track of the number of time derivatives Der in a term -so the quantities we will work with are the coefficients We worked with a predictor-corrector-, as well as a Runge-Kutta-integration algorithm, both with self-adjusting step width. We used two sets of algorithms -one of them involves explicitly programmed versions of the rate equations, while the others worked out the loop-and linkgraphs automatically, only needing the parameters of the physical system. Apart from the algorithm of the flow simulation itself, we developed a number of tools for the analysis of the resulting data. As the coupling space, in which we are working, is very abstract and high-dimensional, it is helpful to make first steps by the study of unphysical toy systems, i.e. simple and solvable physical systems, and reduced turbulent systems (Burgulence) first, to gain the confidence that the algorithm is working correctly, and to develop some intuition for the work with renormalized couplings.

Non-Turbulent Systems Analyzed
We started our investigations by analyzing unphysical (toy-) systems with arbitrary constants and dimensionality of space, to learn more about the detection and features of different sorts of fixed points. A main question was how structures in coupling space can be recognized, if the dimensionality of the coupling space is high, and whether terms of higher order in the field expansion contribute as a correction.
In a second step, we applied the algorithm to standard physical systems with known properties, as the scalar and the O(3)-symmetric field theory. The idea behind this was on the one hand to ensure that the algorithm works correctly and to see how closely we can reproduce analytic values for fixed point scalings; and on the other hand to approach turbulent hydrodynamics in a stepwise manner, interpreting it as a special case of the general 3-vectormodel.

Toy Systems
We worked with a number of unphysical systems for trials of the algorithm and analysis tools, thus merely looking for interesting structures. These systems were defined by an action consisting of a propagator, a two-field-and a four-field-interaction, where the field was a 3-vector-field. Parameters were deliberately adjusted to allow the presence of different fixed points. Investigations of the coupling space were mainly done using the shooting method, which is especially useful for finding fixed points. In practice, one initiates a number of RG-flows, starting from initial conditions sufficiently close to each other, and searches the topology of the flow for interesting structures. An example can be seen in Fig. 11, where the location of the fixed point is easy to guess.
To identify the location more precisely, one would repeat the method with initial conditions closer to the guessed fixed point couplings, leading to a picture similar to Fig. 12. In this way, one approaches the fixed point iteratively. Following this procedure, the simulated trajectories approach the ideal trajectories, i.e. the flows directly running into our out of the fixed point.
The shooting method is limited by the numerical accuracy of the computer program, and the stepsize adjustment of the flow integration, as the algorithm slows down drastically when a fixed point is approached.
In a simulation involving more than two couplings (as is usually the case) the projection of the flow onto a twodimensional subspace will in general not look so evident, but quite similar if the fixed point is approached closely enough.  Fig. 11 after adjusting the initial conditions. Position and ideal trajectories of the fixed point can clearly be identified.
Other interesting structures include focus fixed points; an example can be seen in Fig. 13. As explained by e.g. Sornette [30], these fixed points can be associated with complex eigenvalues, and can thus be identified with some discrete scale invariance. Observe that the simulated trajectory ends at the fixed point due to numerical errors, rather than to encircle it and approaching it asymptotically.

Simple Physical Systems
The renormalization group flows of the scalar theory and the O(3)-symmetric theory in D dimensions have been analyzed by P. Düben [9]. By reproducing known values of these theories like fixed point locations and scaling (also in the ǫ-expansion), we tried to get closer to the much more divert general three-vector-theory, and also to check the correctness of our algorithms. We found that we are able to accurately reproduce the values known from literature, to a given order of the ǫ-expansion. Our results will be published in a forthcoming article.

Local Potential Approximation of Hydrodynamics
We now return to the analysis of the action for hydrodynamics derived above in the LPA. The system is specified by the dimensionality of space and symmetries of the fields; action (82) merely serves as the initial condition of the flow. The simulations were performed using two distinct algorithms: The first one iterating the rate equations derived in the previous sections and doing the book-keeping of the terms involved explicitly, the second one finding the graphs to be computed automatically. The second formulation turned out to be not only more elegant, but a great deal faster than the cumbersome implementation of the book-keeping.
The advantage of our approach is the fast integration of a large number of couplings, and in that way circumventing the drawbacks of the expansions. Simulations were done with up to 100 couplings, though it has to be said that the identification of fixed points becomes nearly impossible in these high-dimensional spaces. Working with such a number of terms can only be done iteratively, meaning that one starts with a low number of couplings, identifies the fixed point and than changes to more and more terms, hoping that these act as corrections to the overall behaviour.  It is a simple calculation to show that for values η > 1.5 of the anomalous exponent, a non-trivial fixed point exists in the vicinity of the trivial one. We used the shooting method to determine the position of this non-trivial fixed point, depending on the anomalous exponent, as can be seen in Fig. 14 and Fig. 15. The distance to the origin of coupling space can be seen to grow linearly with η, but, alas, no physical property could be ascribed to it. Note that for η < 1.5, this fixed point does not exist. As we have already explained, in flow-simulations it is generally preferable to calculate η, rather than to search for it by means of the shooting method. In the LPA, on the other hand, one has η = 0 as no field renormalization is performed.

Scaling of the Trivial Fixed Point
It is straightforward to analyze the scaling of the trivial fixed point. The correlation functions of even orders are directly computed by the RG-flow; after Fouriertransformation to physical space we can read off the scaling, and find:

Order of the Correlation Function
Scaling Exponent 2 0, 666 ± 0, 017 4 1, 338 ± 0, 035 6 1, 999 ± 0, 052 Table 1. Scaling exponents at the trivial fixed point The correlation functions of odd orders are not explicit terms of the action and so have to be measured indirectly. The correlation function of order n can, for example, be derived from the term uv n , if the scaling of the field u is known. We suggest to measure the scaling of u from the two-point-function uu , and subtract it from uv n . In this way the following exponents can be measured:

Conclusion / Outlook
We have shown how to define a generating functional for hydrodynamic turbulence, including a strict treatment of the incompressibility condition. The non-local interactions have been transformed into local ones, by a re-definition of unphysical fields and constants, and by introduction of new fields. We have shown how to approximate the resulting action within a derivative expansion. On the side of the renormalization group, we discussed the procedure of renormalization and rescaling in detail in order to clarify some obscure point in the literature. Our work leads to a RG-equation for the turbulent action, and a set of rate equations after the application of the derivative expansion. These rate equations have been simulated numerically, and the resulting data have been analyzed. As will be published in another article, we tested the numerical algorithm by reproducing known values for nontrivial scalings of the scalar theory in 4 − ǫ dimensions, and the O(3)-symmetric field theory. The results are in agreement with values found in literature, so the numerical algorithm can be considered to be reliable. Our numerical work involves the computation of the RGflow for products of Grassmannian variables, as well as 3-vectors that have to be kept track of according to their combinations.
We were able to identify the trivial fixed-point with the scaling exponents predicted by the K41-theory.
On the other hand, we have to admit that we were not yet able to reproduce intermittent exponents for the structure functions of turbulence that would agree with the experimental values. The reason clearly is the complexity of the general 3-vector-model, including all theories that are based on hydrodynamics. Although the basic foundations of these theories are easily understood, all of them (including Navier-Stokes and Burgers turbulence) involve the same dimensionality of space and symmetry of the fields, while leading to different predictions for the intermittent exponents. Moreover, it is questionable whether the analysis of a fixed point will eventually lead to an understanding of intermittency. Data seem to hint to a cutoff-phenomenon; the probability distribution of the velocity increment looks, for small distances, like a Lévy-distribution; on large scales like normally distributed [12]. We are strongly convinced that this can be explained by a crossover between two fixed points, but cannot yet justify it by our simulations.