A formal derivation of the local energy transfer (LET) theory of homogeneous turbulence

A statistical closure of the Navier–Stokes hierarchy which leads to equations for the two-point, two-time covariance of the velocity field for stationary, homogeneous isotropic turbulence is presented. It is a generalisation of the self-consistent field method due to Edwards (1964) for the stationary, single-time velocity covariance. The probability distribution functional P[u,t] is obtained, in the form of a series, from the Liouville equation by means of a perturbation expansion about a Gaussian distribution, which is chosen to give the exact two-point, two-time covariance. The triple moment is calculated in terms of an ensemble-averaged infinitesimal velocity-field propagator, and shown to yield the Edwards result as a special case. The use of a Gaussian zero-order distribution has been found to justify the introduction of a fluctuation-response relation, which is in accord with modern dynamical theories. In a sense this work completes the analogy drawn by Edwards between turbulence and Brownian motion. Originally Edwards had shown that the noise input was determined by the correlation of the velocity field with the externally applied stirring forces but was unable to determine the system response. Now we find that the system response is determined by the correlation of the velocity field with internal quasi-entropic forces. This analysis is valid to all orders of perturbation theory, and allows the recovery of the local energy transfer (LET) theory, which had previously been derived by more heuristical methods. The LET theory is known to be in good agreement with experimental results. It is also unique among two-point statistical closures in displaying an acceptable (i.e. non-Markovian) relationship between the transfer spectrum and the system response, in accordance with experimental results. As a result of the latter property, it is compatible with the Kolmogorov (K41) spectral phenomenology.


Introduction
It is sometimes remarked that the theory of turbulence is mired in controversy. In reality this is not so. Controversy implies activity. What we actually have is that perceptions of the subject appear to be quite static and dominated by a small number of issues which arose in the 1960s/70s. None of these issues is necessarily of any great significance and some have been effectively resolved. In this respect, the statistical theory is really just like the rest of fundamental research on turbulence. One thinks of topics such as the free decay of the total energy, the dissipation anomaly, and the Kolmogorov (K41) theory: all of these are characterised by disagreements and unresolved issues. For some remarks on this aspect of turbulence research, see the review by Sreenivasan [1].
The reason why there is such a static (and dismissive) perception of the subject is probably two-fold. First, there is the exotic nature of the theoretical physics approach, as it is seen by a turbulence community which is dominated by engineers and applied mathematicians. Indeed, it seems quite unfortunate that those who are best placed to understand the statistical field theoretic approach (e.g. particle theorists, condensed matter theorists) have little knowledge of the phenomenology of turbulence or indeed much interest in it. Conversely, those who are engaged in practical applications, generally do not have the appropriate background to appreciate renormalization methods or statistical field theory.
Secondly, there is the absence of consensus. If there were a clear consensus among theorists, then those concerned with practical applications might well find the statistical theory more accessible. But the small numbers working on this approach, allied to a rather sporadic rate of progress, mean that any sort of consensus is hard to find. Indeed, to those not in the field, going back to the 1960s there must appear to be disagreements between theorists over the reasons why their respective closures failed to yield the Kolmogorov energy spectrum (i.e. Kraichnan [2] and Edwards [3]). Moreover, the perturbative formalism of Wyld [4], which could offer an attractive pedagogic way into the subject, was criticised by Martin, Siggia and Rose [5], in the course of developing their functional formalism. It is only recently that this particular issue has been resolved by Berera, Salewski and McComb [6].
To add to these uncertainties, we have the fact that Kraichnan introduced a new way of representing the triangle condition in wavenumber space. This involved working with the triad of wavenumber amplitudes k, p, and q chosen such that their associated wavevectors formed the sides of a triangle. This leads to quite complicated forms for coefficients which require rather esoteric trigonometric identies in order to evaluate them. In contrast, Edwards followed the existing convention of using wavenumber amplitudes k and j, along with µ = cos θ where θ is the angle between the two wavevectors. The benefit of this formulation of the problem is its simplicity, with it being possible to deduce the conservation properties of the approximation to the inertial term by inspection. It is the Edwards formulation which we follow here.
In the hope of clearing up some of this confusion, we present here a new and fundamental derivation (in the sense of theoretical physics) of an existing theory which first saw the light of day in 1974 [7]. This is the local energy transfer (or LET) theory which was derived by modifying the Edwards single-time theory to conform to the (by that time) established experimental picture of turbulent energy transfer. The result was a theory which was compatible with the Kolmogorov (K41) theory 2 . In 1978, the LET theory was extended heuristically to the two-time case, and in the process led to a new, more unified view of the reasons why the Kraichnan and Edwards theories were not compatible with K41.
In the present work we follow in the footsteps of Edwards [3] and take the probability distribution functional (pdf) of fluctuating velocities to be the solution to the Liouville equation (expressing conservation of probability) and introduce a model pdf, as a non-interacting or quasi-particle representation of the actual pdf, which is constrained to give the same covariance as the exact distribution.
However, our present approach differs from that of Edwards in that we retain the full timedependence of the problem. Also, we do not introduce the Fokker-Planck equation, with associated self-consistent procedures for determining the system response as a dynamical friction. Instead, we argue that the system response can be identified as the ensemble average of the Liouville operator against the Gaussian ground-state pdf and can be calculated as a generalised fluctuation response relation (or FRR). It is perhaps also worth remarking that our approach lies somewhere between that of Edwards and the more general, abstract theory of Herring [8,9]. That is, unlike Herring, we take the model zero-order distribution to be Gaussian; but, unlike Edwards, we leave the base operator undetermined.
Significant new features of the work may be highlighted as follows: 1. In sections 3.2 and 3.3 we present a new analysis of the failure of Kraichnan's DIA, using recent advances in the theory of turbulent dissipation, and also show that this provides a unified explanation of the failure of both DIA and the Edwards SCF. Previously it had been assumed that the diagnosis of the Edwards theory must also apply to DIA as they are cognate theories. 2. In section 4.3 we introduce a more general base operator than that of Edwards to generate the zero-order Gaussian distribution. We show that the Edwards base operator is overdetermined by his assumption of the Fokker-Planck form. 3. In section 4.6, when evaluating the triple moment, we recover the Kraichnan-Wyld-Edwards (KWE) equations for the covariances. But this is no mere re-derivation. The crucial difference is that we obtain them in conjunction with a response function which arises naturally from the functional differentiation and which (unlike previous theories) is compatible with K41 phenomenology and does not have a Markovian relationship to the energy balance equation. 4. In section 4.7 we show that the renormalized response function arises from the covariance of the velocity field with the internal quasi-entropic forces, rather as the energy injection rate arises from the covariance of the velocity field with the external stirring forces. 5. In section 5.1, we provide a formal justification for using the linear fluctuation-response relationship, as this holds to all orders in perturbation theory. We note that in the past such relationships have been used to interpret theories and even on an ad hoc basis to provide a closure approximation. 6. In section 6 we make the new point that a mean-field approximation made by Kraichnan in the course of deriving the DIA is justified by a comparison with the Wyld-Lee formalism where this problem does not arise.
The paper is organised into the following sections: Section 2: we summarise the basic equations and formulate the closure problem in terms of the turbulence pdf, rather than employing the ensemble-based abstract notation which is usual in turbulence theory. We also introduce the velocity-field propagator which is the basis of our present approach. Section 3: We review statistical closures, with particular emphasis on the Eulerian formulations, and present a unified diagnosis of the failure of first-generation theories due respectively to Kraichnan, Wyld, Edwards and Herring. Section 4: A preliminary treatment which shows how the FRR may be applied to turbulence [10] is amplified and extended in order to provide a formal derivation of the two-time LET theory. Section 5: The relationship of this work to other work is discussed, along with a consideration of possible applications.

Basic equations and notation
We consider the solenoidal Navier-Stokes equation (NSE) for the velocity field u α (k, t) in wavenumber (k) space (see either [14] or [17]) as: where ν is the fluid kinematic viscosity, and f α (k, t) is an arbitrarily chosen stirring force which we have to specify. Note that we follow Edwards and use Greek letters to denote the usual Cartesian tensor indices, where these take the values 1, 2 or 3, as appropriate for a threedimensional space. There should be no confusion with the conventional use of Greek indices in Minkowski four-space. The inertial transfer operator M αβγ (k) is given by where i = √ −1, while the projector P αβ (k) is expressed in terms of the Kronecker delta as Note that the use of the projector ensures that the velocity field remains solenoidal. The covariance of the fluctuating velocity field may be introduced as and for isotropic, homogeneous turbulence we may write this as where we anticipate the effects of homogeneity in writing the left hand side. As is usual, the angle brackets . . . denote the operation of taking an average. We shall discuss this further in the next section, in connection with the introduction of the probability distribution functional. If we consider the case t = t , then we may introduce the spectral density function: (2.6) and the energy spectrum: At this stage, Kraichnan introduced his infinitesimal response tensor which related an infinitesimal fluctuation in the velocity field to a corresponding fluctuation in the stirring forces. Instead of doing this, we will introduce an infinitesimal velocity field propagator as determining the system response. In its most general form this is defined by the relationship between two infinitesimal fluctuations of the velocity field, as: for t t, (2.8) where the hat symbol indicates that the propagator R αβ is a random variable. Further, we rewrite equation (2.8) in terms of the functional derivative, in order to introduce the ensembleaveraged response function, thus: where we have in turn invoked homogeneity and then isotropy. The quantity R(k; t, t ) is called the response function.

The stirring forces
In order to define the ensemble (and, when required, to study stationary turbulence) Edwards [3] introduced stirring forces. These are denoted by f α (k, t) and can be added to the right hand side of the equation of motion, as we have done here in (2.1). These forces must be chosen to be isotropic, homogeneous and (in order to maintain incompressibility) solenoidal. We consider random forces with a multivariate normal probability distribution (that is, the pdf for each mode k is Gaussian) such that the associated functional integrals are analytically tractable. It is also usual to assume that the autocorrelation of the forces is instantaneous in time and we represent this by choosing the time autocorrelation to be a delta function. In signal analysis this case is known as 'white noise', but it was introduced into turbulence theory by Edwards.
A form of correlation which satisfies all these requirements is: Here F(k) is a spectral energy density which is related to the rate at which the force does work on the fluid. The rate at which the stirring forces do work on the fluid can be shown to take the form [3,18]: We may also introduce the energy injection spectrum of the stirring forces W(k) as: with associated rate of doing work ε W , where: For stationary flows ε W must be equal to the dissipation rate ε.

The probability distribution functional (pdf) and the moment hierarchy
We characterise the system by its probability distribution functional (pdf) which we denote by P[u, t], or P for short. It may be defined, in the language of statistical mechanics, as: P [u, t] = the probability that the 'phase' of the system lies between u(k, t) and u(k, t) + du(k, t). It is the pdf which is required to evaluate the covariance C αβ (k, k ; t, t ), and other associated statistical quantities. Hence, although the symbol . . . is normally used in turbulence theory to indicate an ensemble average, our use of angle brackets in this work will always imply an average against the system pdf.
We now consider evaluating moments in the context of the closure problem. In a simplified notation, with Du representing the functional integral, we have: and so on. Problems in many-body physics are generally (as in turbulence) mathematically intractable, and rely on an approximate technique, such as mean-field theory, or self-consistent field and perturbation methods. Often a combination of these approaches is required for success. It is universally accepted that the ideal starting point is to find a model system which is, in some sense, close to the actual system; but is, despite that, mathematically tractable.
In the case of Wyld's perturbative analysis of the NSE [4], only the second of these criteria is satisfied. The zero-order (or ground-state) model is characterised by the molecular viscosity and is therefore very remote from the turbulent system it is modelling. Of course it is an easy matter to generate the perturbation series in terms of a book-keeping parameter associated with the nonlinear term. The price one then pays is the enormous mathematical complexity of the resulting series and the need to find ways to sum classes of terms to all orders in order to achieve renormalization.
The technique pioneered by Edwards satisfies both criteria for a base model. But the way in which it does this is rather subtle, and relies on underlying symmetry considerations, which only emerge when one works through his analysis in some detail. As we are using a similar method, based on that of Edwards, it may be helpful to make the symmetry aspects manifest from the beginning, and this we now do.
Let us formally decompose the pdf into symmetric and antisymmetric parts, as denoted by subscripts, thus: (2.18) Then, from elementary symmetry considerations, it follows that equations (2.14)-(2.15) may be rewritten as: and so on. In other words, the symmetric part of the distribution P S determines the even-order moments, whereas the antisymmetric part P A determines the odd-order moments. This decomposition can be interpreted in terms of the perturbation expansion of the pdf, which we introduce later as (4.5), where P S consists of the even order terms P 0 , P 2 , . . . while P A comprises the odd-order terms.

Statement of the statistical closure problem
This material may be found in the literature but it will be convenient to repeat it here, partly to define some notation, but mainly to allow us to make some points about the assessment of closure theories which are often overlooked.
We form an equation for the covariance C(k; t, t ) in the usual way. Multiply each term in (2.1) by u α (−k, t ) and take the average, to obtain: where we have used equations (2.3) and (2.5), cancelled the factor δ(k + k ) across, and invoked isotropy, along with the property Tr P αβ (k) = 2. We can also write this in the compact form: where P(k; t, t ) is just the right hand side of (2.23). The problem of expressing this in terms of the covariance is the well-known statistical closure problem. Also using (2.1), we can obtain an equation describing the energy balance. On the timediagonal, the covariance equation takes the form for isotropic turbulence: where F(k) is the energy injection spectral density, as defined in terms of the stirring forces by equation (2.11). Again, we can write this in a more compact form as: (2.26) where Q(k; t, t) can be determined by comparison with the right hand side of (2.25). The energy spectrum is related to the spectral density by equation (2.7). Accordingly, if we multiply (2.25) across by 4πk 2 , and rearrange terms, the governing equation of the energy spectrum takes the form: (2.27) where the energy transfer spectrum is T(k, t) = 4πk 2 Q(k; t, t) ≡ 4πk 2 Q(k, t). This equation for the energy spectrum is well known, and is nowadays sometimes referred to as the Lin equation [15,17]. More detailed discussions of these equations can be found in chapter 3 of the book [17].

A concise review of statistical closures
The statistical closure problem was first formulated for turbulent shear flows by Reynolds in the 1890s. He showed that the equation for the mean velocity Ū contains the unknown covariance of two fluctuating velocities uu : that is, the Reynolds stress. It was later formulated for isotropic turbulence by Taylor in the 1930s. In this case the mean velocity is taken to be zero and the equation for the correlation of two velocities uu is found to contain the unknown correlation of three velocities uuu . The equation for uuu contains the unknown fourth-order moment uuuu , and so on. The problem is therefore seen as an open hierarchy of statistical equations (one-point for Reynolds, two-point for Taylor) which requires some means of closure. Theories used in engineering applications, such as the eddy-viscosity, mixing-length and n-equation models, are all effectively statistical closure approximations for the single-point problem. In contrast, the Heisenberg eddy viscosity and the quasi-normality theory are effectively closure approximations for the two-point problem. Further reading on this topic can be found in the following books [11,12,14,17]. The first formal treatment of the closure problem was the quasi-normality hypothesis (Proudman and Reid [19], Tatsumi [20]). The basic idea was: solve the next equation in the hierarchy for uuu by factorizing uuuu in terms of uu × uu . This then gives a closed set of equations for the covariance. However, quasi-normality failed when computed numerically in the 1960s, as it predicted negative spectra: it was not physically realizable [21,22].
The attempt at quasi-normality illustrates an important aspect of all statistical theories of turbulence and indeed of all statistical field theories. In all theories the averages are evaluated using the tractable properties of the Gaussian or normal distribution. In particular, even-order moments can be factorised into products of pair-correlations. But odd-order moments cannot be evaluated that way, because they vanish identically by symmetry. In quasi-normality the fourth-order moment is evaluated in terms of pair-correlations, whereas the triple moment is calculated from the resulting closed equation by inverting the linear operator of the Navier-Stokes equation. Some such manoeuvre is needed in any theory if it is to be successful, with the degree of success depending on the choice of manoeuvre.

The DIA of Kraichnan and the SCF theories of Edwards and Herring
The beginning of the modern age of turbulence theory was signalled by the publication in 1959 of Kraichnan's direct interaction approximation (DIA) [23]. This was both a mean-field theory and a new kind of perturbation theory 3 . At the time its most important feature was seen as being its physical realizability. It did not have the catastrophic behaviour of quasinormality. Shortly afterwards, Wyld [4] showed that the DIA could be recovered from a conventional renormalized perturbation theory. This work was couched in the language (and used the techniques) of quantum field theory, including the introduction of diagrams analogous to Feynmann diagrams. Although of considerable pedagogical interest, it contained some procedural errors. A discussion of this aspect may be found in section 10.2 of the book [17], but these errors do not affect anything we say here.
In 1964 Edwards published his self-consistent field (SCF) theory [3], which he also described as a random phase approximation, and which was based on the probability distribution functional (pdf) of the fluctuating velocities. Employing a technique which has since become commonplace in statistical field theory, he began by introducing a model system with a Gaussian pdf and the same velocity covariance as the turbulence system. Then he introduced a perturbation expansion in terms of the difference between the model system and the real system. A particularly interesting feature of Edwards's work was the way in which he exploited the symmetries of the problem to compensate for the lack of a variational principle, such as the minimum free energy, as possessed by microscopic systems in thermal equilibrium. The result was a set of equations for the single-time covariance and a renormalized response function. It should be noted that this result was cognate to the DIA theory of Kraichnan and, like it, maintained the properties of the inertial term, such as conservation of energy.
Another attempt at a self-consistent field theory was made by Herring [8] in 1965. This was a more abstract approach than that of Edwards, although it also started from the Edwards form of the Liouville equation. In strategy it lay somewhere between the Edwards SCF and the Kraichnan DIA, in that it formally renormalised a 'bare' operator. Herring's basic approach involved the introduction of single-mode projectors and single-mode distributions. Then the perturbation expansion was in terms of the difference between single-mode and coupled forms. There are similarities and differences between the two forms of SCF, and they both lead to the same form of covariance equation. However, the response function in Herring's theory was the same as the time-independent DIA form and hence differed from that of Edwards: see section 4.3 of [14] for a fuller account of this work.
In the following year, Herring was able to extend his SCF theory to the two-time case [9], the result being a theory which was closely related to Kraichnan's DIA, rather then the Edwards theory. In 2002 Edwards (with Schwartz) [24] revisited the problem in the more general context of noise-driven nonlinear field equations, but did not succeed in generalising his theory to the two-time case in a way comparable to the DIA of Kraichnan and the SCF of Herring. In effect, such a generalisation is the subject of the present paper.

The pioneering closures and the Kolmogorov (K41) theory
Historically, the primary test of a spectral closure approximation lay in the answer to the question: does it possess the Kolmogorov '-5/3' spectrum as its asymptotic solution in the limit of infinite Reynolds numbers? This consideration dominated thinking in the subject during the 1960s. But later on in the 1980s, developments in the study of structure functions led to this view being challenged. In particular, the behaviour of power-law exponents with increasing order was believed by some to indicate the presence of intermittency effects, which in wavenumber space would change the '−5/3' spectrum to a slightly different value. A discussion of this topic will be found in section 3.2.1 of the book [14], and later work based on multifractal models is covered in the review by Boffetta, Mazzino and Vulpiani [25].
However, in recent years there as been a growing interest in finite-Reynolds-number corrections, to such an extent that the view that K41 is not a test of closures seems increasingly dated. Theoretical treatments have been given by Effinger and Grossmann [26], Barenblatt and Chorin [27], Qian [28], Gamard and George [29], and Lundgren [30]. The last-named is particularly interesting as he used matched asymptotic expansions to show that the second-and third-order exponents were of K41 form, with only finite-Reynolds-number corrections, and this was supported by a comparison with the experimental results of Mydlarski and Warhaft [31]. Experimental and numerical investigations of finite-Reynolds-number effects have been made by Gotoh, Fukayama and Nakano [32], Antonia and Burattini [33], Tchoufag, Sagaut and Cambon [34], McComb et al [35], and Antonia et al [36]. In particular, McComb et al [35], present results which suggest that systematic error with increasing order of structure function is more likely than anomalous exponents, while Antonia et al [36] conclude that a failure to take account of the effects of finite Reynolds numbers has resulted in misguided assessments of the Kolmogorov theory.
Here we separately discuss the nature of the difficulty for DIA, and then for the Edwards SCF, of showing compatibility with K41.

Incompatibility of DIA with K41.
Kraichnan had originally concluded [23] that the DIA predicted an inertial-range spectrum of the form where C is a constant, and U stands for the root-mean-square velocity. In 1964, he further concluded [2] that: Recent experimental evidence gives strong support to [the Kolmogorov '−5/3' form] and rules out [the '−3/2' form above] as a correct asymptotic law.
The experimental evidence that he cited was the iconic investigation of Grant, Stewart and Moillet [37]; and, of course, since that time there has been overwhelming evidence to support the K41 form of the energy spectrum: see the survey papers by Sreenivasan [38] and Yeung and Zhou [39].
In examining the inertial range behaviour of his theory, Kraichnan followed the example of Batchelor [40], and introduced the flux of energy through wavenumber κ, which he denoted by the symbol Π, as defined by Then the condition for an inertial range is given in terms of the maximum value of the flux by: This condition, when it holds for a range of wavenumbers, is known as scale invariance.
A version of Kraichnan's analysis leading to equation (3.1) can be found translated into the simpler Edwards formalism in section 6.1.6 of the book by McComb [14]. This analysis led to an asymptotic form of the response function G(k, τ ) as where τ = t − t and J 1 is the first-order Bessel function of the first kind. We should note the presence of the viscous time-scale (νk 2 ) −1 and the energy-containing range time-scale For large values of k the exponential factor may be treated as unity, and Kraichnan concluded that the two-time correlation and response functions scaled on (Uk) −1 ; also known as the convective or sweeping time-scale. This led Kraichnan to the conclusion that a satisfactory closure could not be obtained in the Eulerian coordinate system. This seems rather counter-intuitive, as this is the coordinate system which is universally employed throughout fluid dynamics. His way out of this 'impasse' was to apply the ideas of DIA in a Lagrangian-history coordinate system. Here we will consider some of the consequences of this conclusion.
In particular, it stimulated interest in the scaling of two-time correlations. This is an important subject in its own right. But, in certain quarters, the above analysis did tend to polarise matters into 'convective' versus 'Kolmogorov' scaling. Or indeed, an assumption that Eulerian correlations scale on (Uk) −1 in the inertial range, whereas Lagrangian correlations were supposed to scale on (ε 1/3 k 2/3 ) −1 . In reality, the situation is a great deal more complicated than that, as a glance at some recent review articles on the subject will show [41][42][43]. However, we will make some specific points here, as derived from our own earlier work.
In order to do this, it is helpful to put (3.1) into a form which is more immediately comparable with the Kolmogorov spectrum. We may eliminate the dependence on the rms velocity U by substituting from the asymptotic form of the Taylor surrogate for the dissipation rate (see McComb, Berera, Yoffe and Linkmann [44]), where L is the integral length-scale and C ε,∞ is a constant. We can then write (3.1) as where µ = 1/6. Kraichnan only computed the Eulerian DIA for free decay at low Reynolds numbers (we will discuss this later). However, in 1989 McComb, Shanmugasundaram and Hutchinson [45] reported calculations for free decay of both DIA and LET for Taylor-Reynolds numbers in the range 0.5 R λ (t f ) 1009, where t f is the final time of the computation. These results do not support the asymptotic form of the DIA energy spectrum, as given above by (3.1). It was found that (for example) at R λ (t f ) = 533, the two theories were virtually indistinguishable and both gave the Kolmogorov spectrum to within the accuracy of the numerical methods. It was shown that this result was not an artefact of the initial conditions, by taking k −3/2 as the initial spectrum, whereupon it was found that both theories evolved away from this form to once again give k −5/3 as the final spectrum.
In that investigation, two-time correlations for both theories, as well as propagator (LET) and response (DIA) functions, were tested with both types of scaling. It was found that both theories behaved in a very similar way. Preliminary conclusions were drawn as follows: 1. Both scaling methods were partially successful, with convective scaling tending to collapse data at lower wavenumbers, and Kolmogorov scaling being more effective at higher wavenumbers. 2. Overall, Kolmogorov scaling seemed to be the more effective of the two.
3. The difference between the two methods increased with increasing Taylor-Reynolds number, particularly at shorter lag times. 4. When a shorter reference time was taken (t f = 0.01 as opposed to t f = 0.06), the position was reversed, with convective scaling behaviour becoming more effective, even at large wavenumbers.
In considering Kraichan's approximate analysis of the inertial range behaviour of the DIA, we should also note the unphysical nature of the response function (3.4), which is oscillatory when a monotonic decline would be expected; and that the condition (3.2) is inconsistent with the −3/2 form, when it actually leads to −5/3 [46]. We make these points to emphasise the fact that Kraichnan's analysis of the failure of DIA to yield the Kolmogorov spectrum is neither definitive nor prescriptive for other theories. We may obtain a different view of this problem by considering the analogous analysis of the Edwards theory [47].

Incompatibility of the Edwards SCF with K41.
In a later paper [47], Edwards examined the properties of his equations in different situations and with different inputs. A particularly interesting and important feature, was his formulation of the infinite Reynolds number limit. Batchelor [40] had in 1953 observed that taking a limit of the viscosity shrinking to zero while maintaining the dissipation rate constant, must necessarily concentrate the dissipation at infinite wavenmbers. Edwards formalised this idea by introducing the delta function δ(k − ∞) and balancing it by another at the origin, which acted as an input. Thus for stationary stationary turbulence at infinite Reynolds numbers, the transfer spectrum T(k) could be written as: where ε W is the rate of doing work by the stirring forces, and ε is the dissipation rate. For stationarity the two rates must be the equal; and Edwards just used the symbol for the dissipation, as is usual in the subject 4 . Edwards argued that under these circumstances the Kolmogorov spectrum E(k) ∼ ε 2/3 k −5/3 should apply at all wavenumbers, and demonstrated that this was the case for his covariance equation, with the transfer spectrum taking the form given in (3.7). The problem was that his equation for the renormalized response involved an integral which diverged as the wavenumber tended to zero. This is (by analogy) often referred to as an infra-red divergence. Thus, despite its other achievements, the Edwards SCF theory did not possess the correct limiting behaviour at infinite Reynolds numbers.
This was necessarily true also of the DIA, as they are cognate forms. It is a trivial matter to show that the Edwards criterion can be integrated to give (3.3). So, in that respect at least, the two approaches are mathematically equivalent. A fuller discussion can be found in section 4.3 of [17].
As we have seen, Kraichnan and Edwards analysed the failure of their respective theories in different ways. In Kraichnan's case this led to a much more elaborate approach than the simpler Edwards view, which merely required some physical basis for introducing a counterterm to the response equation in order to cancel the infra-red divergence. In contrast, Kraichnan examined the behaviour of the DIA in terms of its ability to distinguish between convection and energy transfer and on the basis of much approximate analysis concluded that theories should be invariant under random Galilean transformations. He also concluded that Eulerian two-time theories could not satisfy this requirement, although single-time theories could [12,14].
As mentioned earlier, Kraichnan's response to this difficulty was to develop Lagrangianhistory theories, and in this he was followed by others, especially Kaneda and Kida [48][49][50][51]. As such theories are random Galilean invariant, the implication is that they are compatible with the Kolmogorov spectrum. It should be emphasised that all turbulence theories are Galilean invariant in the ordinary sense. The concept of random Galilean invariance conflates the exact symmetry of Galilean invariance with the definition of the statistical ensemble, and affects the operation of the ensemble average. For further discussion, see chapter 5 of [17], and also the paper by Frederiksen and Davies [52].

The original derivation of the LET theory
In 1974 it was shown by McComb [7] that the relationship between the response function and the energy balance equation in the Edwards SCF theory was incompatible with the known phenomenology of turbulence. As we know, SCF and DIA are cognate theories, so this conclusion also applied to the DIA. We shall discuss these matters in more detail presently. But for the moment we note that the physically correct use of the energy balance equation to define the response function led to the local energy transfer or LET theory. This was later extended in a rather heuristic way to the two-time case [53].
In order to understand the significance of this result, and how it led on to the LET theory, it is helpful to recall the form taken by the spectral energy balance (or Lin equation) in the presence of an energy input, as given by equation (2.27). For the stationary case, we have: Our immediate purpose is to understand the Edwards theory in terms of this phenomenology, and that essentially amounts to writing down the Edwards expression for T(k). This is given in appendix, and from (A.3) applied to the stationary case this is: where D(k, j, |k − j|) is given by (A.5). Note that, due to the antisymmetry of the integrand under interchange of k and j, the transfer spectrum satisfies as required for conservation of energy. The Edwards closure approximation was originally derived as follows [3]. Drawing an analogy with the Fokker-Planck equation, Edwards introduced diffusion and dynamical friction terms, s(k) and r(k) respectively, as: and Our definition of the dynamical friction differs from that of Edwards, only insofar as it permits us to write the transfer spectrum in terms of the energy spectrum E(k), rather than the covariance, thus: The interpretation of this result is crucial to an understanding of the Edwards approach (and, by extension, those of Kraichnan and Herring). We see that the transfer spectrum has been divided into two parts. The first part s(k) transfers energy into mode k by coupling to all other modes. Whereas, the second part involves a coefficient r(k), which arises from a sum over modes, which multiplies E(k) and which acts like an effective viscosity. Edwards interpreted it as such, and wrote an expression for the inverse lifetime ω(k) of mode k as (3.14) A physical interpretation was given by Edwards [47]. For this we will find it helpful to restore the time derivative on the left hand side, as in equation (A.3). We can think of this as being a slow variation near the steady state, and write the spectral energy balance as: where we have substituted from (3.13) for the transfer spectrum. It may be interpreted as follows: The rate of change of energy in mode k = The direct injection of energy into mode k + the transfer into mode k from all other modes − the transfer out of mode k to other modes − the direct dissipation of energy in mode k due to molecular viscosity.
A key point is that the loss of energy from mode k is proportional to the amount of energy contained in that mode. It was pointed out by Herring [8], that there were many ways of choosing s(k) and r(k) such that the appropriate constraints were satisfied. Yet this particular decomposition must have seemed very natural at the time, because this is precisely the form that transport equations take in statistical physics. One has, for instance, the Fermi master equation, the Chapman-Kolmogorov equation and the Boltzmann equation. Indeed, Edwards drew a specific analogy between the above result and the Peierls-Boltzmann equation for phonon transport in solids, which also involves triadic mode interactions [47]. However, also at that time, the first measurements of the turbulence energy balance were just beginning to be published and it is only with hindsight that one can see that these interpretations of closures were incompatible with the phenomenology of turbulence. In a pioneering paper in 1963, Uberoi [54] obtained T(k), for grid turbulence, from the Lin equation, by measuring the local (in wavenumber) dissipation and decay rates. The resulting transfer spectrum was found to behave like a sink of energy at low wavenumbers and as a source of energy at high wavenumbers. This behaviour has since received abundant confirmation from direct numerical simulations.
The preceding interpretation of the single-time Edwards theory also applies to DIA in twotime form. Although it is easier to understand these matters in the simpler formulation due to Edwards, nevertheless the DIA was also explicitly interpreted in purely Markovian terms by Kraichnan [23]. We may quote the relevant passage as follows: The net flow is the resultant of these absorption and emission terms. It will be noticed that in contrast to the absorption term, the emission terms are proportional to E(k). This indicates that the energy exchange acts to maintain equilibrium. If the spectrum level were suddenly raised to much higher than the equilibrium value in a narrow neighbourhood of k, the emission terms would be greatly increased while the absorption term would be little affected, thus energy would be drained from the neighbourhood and equilibrium re-established. The structure of the emission and absorption terms is such that we may expect the energy flow to be from strongly to weakly excited modes, in accord with general statistical mechanical principles.
This argument is essentially a more elaborate version of that due to Edwards, and presents what is very much a Markovian picture of turbulence energy transfer. But in later years, numerical experiments based on high-resolution direct numerical simulations, did not bear that picture out. In particular, we note the investigation by Kuczaj et al [55].
We can sum up the situation regarding the failure of the pioneering closures as follows. Their form of T(k) is only valid for Markov processes (with some examples listed above), so it is incompatible with the nature of turbulence which is non-Markovian. It is also incompatible with the phenomenology of turbulence, where the entire T(k) acts as input (or output), depending on the value of k [7]. This is in fact the basic flaw in both the DIA and the SCF theories: the fault lies not in the covariance equations but in the relationship of the response function to them.
The LET theory was introduced with the hypothesis that ω(k) is determined by T(k) and can be defined by a local energy balance [7]. It was extended to the two-time case [53]; and, less heuristically, in subsequent papers by McComb and co-workers [56][57][58][59]. Essentially, the two-time LET theory comprises the DIA covariance equations plus the generalized fluctuation-response relation. It may be compared to Herring's SCF [9] which comprises the DIA response equation, single-time covariance equation and the generalized fluctuation-response equation. It may also be compared directly to DIA in terms of response equations. However, before doing that, we will go back to the simplest case, and show how LET arose in relation to the Edwards SCF.
It was anticipated by McComb [7], that a correct assignment of the system response in terms of T(k) (i.e. 'correct' in the sense of agreeing with the turbulence phenomenology of energy transfer) could lead to a response function which was compatible with K41. This was found to be the case and, citing the form given in [60], we may write for the turbulence viscosity ν T (k): where ω(k) = ν T (k)k 2 . The lower limit on the integral with respect to j arises when we consider the flux through mode k. It was used in [7] to justify wavenumber expansions leading to differential forms but is not needed here and can be dropped. The interesting point here is made by rewriting this in terms of the Edwards dynamical friction r(k). From equation (3.12), with the substitution of (A.5) for the memory function D(k, j, |k − j|), we may rewrite (3.16) as: It was shown [7] that the second term cancelled the divergence in r(k). The question now arises: how does one extend this idea to the two-time covariance? The first step in doing this was based on the introduction of a velocity field propagator [53]; but, although this led to the LET governing equation for the turbulence response function R(k; t, t ), the concept was in itself unsatisfactory, and was later replaced by a propagator acting directly on velocity-field covariance [56], thus: C(k; t, t ) = R(k; t, t )C(k; t, t). (3.18) This also takes the form of the fluctuation-response relationship, and this interpretation was further developed by McComb and Kiyani [58,59]. Just as we did for the stationary LET and Edwards SCF theories, it is of interest to compare the DIA and LET response equations. We may write the DIA response equation in its usual form [23]: where G is the DIA response function. This may be compared with the LET response equation, which takes the form [59]: (3.20) If we compare (3.20) with the DIA response equation (3.19), the additional terms on the right hand side of (3.20) give rise to those terms in the single-time theory which cancel the infra-red divergence and so ensure compatibility with the Kolmogorov K41 spectrum. However, a more significant aspect is that these additional terms mean that, for the two-time LET, the relationship between the transfer spectrum and the response equation is in agreement with the turbulence phenomenology as established by experiment and direct numerical simulation.
Is worth noting that the fluctuation-response relation can be used directly to calculate the LET theory instead of the above response equation. With this simplification LET is much easier to calculate than DIA. So, despite the extra complication, the two-time LET bears the same relationship to DIA as its single-time, stationary version does to the Edwards theory.

Numerical computation of the Eulerian closures
The numerical study of renormalized closure approximations was begun in 1964 by Kraichnan [61], who used DIA to calculate the free decay of turbulence from a Gaussian initial state with an arbitrarily chosen energy spectrum. Three such initial spectra were considered, and in all cases the qualitative behaviour was impressive, with the nonlinear coupling transferring energy to higher wavenumbers than those initially excited, and the development of the transfer spectrum from zero to a form in agreement with experiment. The development of selfpreservation and self-similarity were both observed during the decay. Quantitative agreement with experiment was as good as the agreement between experiments, and overall the results were very satisfactory. As the values of Taylor-Reynolds number were less than 40, the question of K41 behaviour did not arise. Nevertheless, this was a truly historic investigation of a theory which was derived from general principles, and which did not possess any adjustable parameters.
A similar but more general investigation was carried out by Herring and Kraichnan [62], who considered not only DIA, but also Herring's SCF; a model derived from the Edwards SCF and denoted by EDW; along with Kraichnan's test-field model (or TFM), which was a single-time theory derived from the DIA and containing one adjustable parameter 5 . There are two features of this investigation which are particularly worthy of note.
First, there was the emphasis on the skewness as providing the most sensitive criterion for deciding between theories. Secondly, there was their comparison with the results from the first DNS of isotropic turbulence, which had just been performed by Orszag and Patterson [63]. As well as this, there were comparisons with results from laboratory experiments up to Taylor-Reynolds numbers of about 50. In general, all the theories performed well in producing qualitatively correct forms of the evolving spectra. For quantities like the dissipation spectrum, all theories lay within the scatter of the experiments. In the case of the skewness, there was a possibility of comparision with the numerical simulation for values of the Taylor-Reynolds number of around 18. All theories behaved much like the simulation but took a range of asymptotic values, from TFM (which was TFM with a different value of the adjustable parameter) and reached a value just above that of the simulation, to EDW which lay well below the simulation. The skewness factor can also be seen as a measure of the effectiveness of theories at transferring energy through wavenumber. Calculation of transfer spectra and transport power revealed the EDW transferred the least energy, DIA and SCF were intermediate, and TFM and TFM were the most efficient.
In 1984, some twenty years after Kraichnan's first numerical computation of the DIA [61], a similar calculation was carried out for the LET, in comparison with DIA, and employing identical methods to those of Kraichnan. At low to moderate Reynolds numbers, the DIA and the LET theories were found to be very similar. Comparison of the newer DIA results with the earlier computation by Kraichnan revealed very similar behaviour, but with slight systematic differences which were attributed to the use of the different formulations: (k, j, µ), due to Edwards, and (k, j, l), due to Kraichnan 6 . At very high Reynolds numbers, the LET results were found to agree well with experiment and also with the Lagrangian-history theories [64]. Indeed, surprisingly, the purely Eulerian LET theory was found to agree rather closely with the strain-based Lagrangian-history direct interaction approximation, and further comparisons showed that this agreement extended to low Reynolds numbers as well.
These results were for free decay, from a variety of arbitrarily chosen initial spectra, at evolved Taylor-Reynolds numbers up to R λ = 533. Subsequent investigations followed the same format but explored other aspects. In 1989 [45], the concentration was on the two-time correlation and response functions over the range 0.5 R λ 1009. It was found that convective scaling was less important than inertial-range scaling, and this led on to a critical analysis of the concept of random Galilean invariance. We have previously referred to these results in section 3.2.1.
In 1992, as well as introducing the fluctuation-response relationship as the defining step of the LET theory, and providing a new simplified derivation (along with making some minor corrections to the equations), the work was extended to the case of passive scalar convection [56]. Results were obtained for free decay for 5 R λ 1060 and for Prandtl numbers of 0.1, 0.5 and 1.0. In both velocity and scalar spectra, Kolmogorov-type power laws were found at the higher Reynolds numbers. In the same year, Shanmugasundaram [65] published the results of an investigation into the processes of energy transfer in LET. He found that energy transfer 5 They also calculated an abridged form of Kraichnan's Lagrangian theory at somewhat higher Reynolds numbers, but we will not discuss that here. 6 For a detailed discussion of these formulations, see section 3.5 of [17]. was overwhelmingly local in character. At higher Reynolds numbers there were some indications of reverse nonlocal transfers, but these did not make a significant contribution to the transfer spectrum T(k). A particularly interesting result was obtained from the computation of Kraichnan's energy transfer locality function [66]. This was found to agree well with results from DNS and from previous computations by Kraichnan of his test field model closure. The author concluded by saying: …our picture of the modal interactions and the consequent energy transfers is not very far from the classical picture …and is a clear demonstration of the ability of the LET theory not only to predict the spectra correctly but also to resolve the underlying physical mechanisms.
In a later investigation which was reported in 2003, there were two new features. First, there was a detailed comparison of free decay results with an in-house direct numerical simulation of the Navier-Stokes equations with, as far as possible, the same initial conditions. Secondly, there was a preliminary attempt to apply the DIA and LET closures to stationary turbulence. This raised problems regarding the computation of time-history integrals over a long transient period, and for this reason the results were seen as being rather tentative in nature. In short, for free decay both closures were found to agree with the simulation within experimental error, as found in earlier investigations. In the forced case, agreement was not so good, except for high-wavenumber spectra and quantities derived from this region such as the Taylor microscale. This alone would suggest that the form of the low-wavenumber forcing needs further attention.
In 2004 Kiyani and McComb [58] reported some technical improvements to the LET, mainly to do with the introduction of time-ordering to the fluctuation-response relation.
The following year, McComb and Kiyani [59] assumed exponential time-dependences and formally verified that there was no infra-red divergence. They also derived a single-time Markovianised form of the LET. This latter theory still has to be computed for a standard test problem, such as free decay from a given initial spectrum.

A formal derivation of the LET as a self-consistent field theory
Essentially we now develop our approach in parallel with an account of the Edwards method for an invariant exact pdf, which we wish to extend to the non-invariant case.

The Liouville equation for the exact pdf
As is well known, the Liouville equation expresses conservation of probability in the form where P is the pdf of the system, and the evaluation of the total time-derivative depends on the nature of the system. It is most familiar in the context of statistical mechanics where it is usually applied to the N-body Hamiltonian system, the Liouvillian operator L being evaluated with the aid of Hamilton's equations. It was shown by Edwards [3] that the analogous equation describing the evolution of the probability density functional P [u, t] of the velocity field for stirred fluid motion could be written as where L = V + L, for this system, and all the operators are defined by their action on P. Thus: and (4.4) It is important to note that the Liouville equation is linear in the pdf.
In addition to the original treatment by Edwards [3], derivations can be found in the books by Leslie [12] and by McComb [14] where, with some changes of notation, it appears as equation (6.71). It can also be derived by consideration of the characteristic (or generating) functional: see the book by Beran [11]. For sake of completeness, we note that Beran [67] has also demonstrated that the method of Edwards can be used to recover the familiar result for the N-body Hamiltonian system.
At this stage, as we recognize that P [u, t] is non-Gaussian, with non-zero skewness and a flatness factor that differs from the Gaussian value of 3, we have to calculate the effects of the nonlinear coupling which should lead to corrections to the zero-order, or base, Gaussian model distribution. We do this iteratively in the next section.

Perturbation expansion of the Liouville equation
Historically, the concept of perturbation theory has been associated with a small perturbation. However, with the growth of statistical field theory from the 1970s onwards, the term has been more generally associated with an iterative calculation based on an expansion parameter, which is used to group terms of the same order, and which is subsequently put equal to one. In particular, this has been the standard usage in turbulence theory for many years.
In order to follow the general strategy of Edwards, we expand the exact pdf about a Gaussian model distribution P 0 [u] in a perturbation series, thus: where ε is a book-keeping parameter which is later set equal to one 7 . As P 0 is an invariant distribution, the time-dependence of the pdf enters through the higher-order terms. Since the model distribution is both normalised, and chosen to recover the exact covariance, the higherorder coefficients are constrained, as follows, and 7 The use of this type of ordering parameter, which is treated as if small during the calculation and set to one at the end, is discussed in the books by Leslie [12] and McComb [14].
for integer i = 1, 2, 3, . . . ∞. We recall that the odd orders cannot contribute to these integrals since they are antisymmetric in u.
We introduce the zero-order operator associated with the model which is the basis of our perturbation theory through the relationship (4.8) We shall refer to the zero-order operator and distribution as the base operator and the base distribution respectively. Further details of their forms are given in the next section. The operator L 0 is now introduced into the Liouville equation, as given by equation (4.2), by simultaneously adding and subtracting L 0 P. With some rearrangement it becomes: For later convenience, we will rewrite this equation as: It should be noted that we can also write This follows from the invariance of the zero-order distribution and the definition of the L 0 operator, as given by (4.8).
Since the operator V stems from the non-linear term and is antisymmetric in u, we assign order ε to it. The operator L − L 0 generates a correction to the flatness factor of the pdf and is taken to be of order 2 . These two terms represent the coupling of the model (i.e. Gaussian) system to the exact Navier-Stokes system.
The perturbation expansion of the pdf is then substituted for P [u, t], such that Collecting terms with the same order of ε, we find 0 : L t,0 P 0 = 0, (4.14) 1 : 16) and so on, to all orders. However, we shall only work to order 2 . Recall that equation (4.14) at O 0 is satisfied trivially, because both the partial time derivative and the operator L 0 vanish independently when acting on P 0 . Formally we may obtain the first-and second-order coefficients from (4.15) and (4.16) as: (4.18) where in the last line we invoke L 0 P 0 = 0. From this we may write the exact pdf as an expansion, thus: (4.19) where the coefficients are given by: (4.21) and so on, to all orders in powers of ε.

The model system
Following Edwards [3], we now define our model system in terms of a Gaussian distribution P 0 [u], which is chosen such that it is normalised to unity and recovers the exact covariance. That is: and respectively. This is the introduction of the two-time covariance. However, in order to explain the Edwards model, we revert temporarily to the single-time case. For the sake of simplicitly, we employ the reduced notation of Herring [8], as used extensively by Leslie [12] and others. In this notation we represent the velocity field by X i , where the index is a combined wave-vector and cartesian tensor index (i.e. our k and α). Accordingly, we introduce the Edwards-Fokker-Planck operator as the sum of single-mode operators, in the form: where ω i is a renormalized eddy decay rate and φ i is the covariance of the velocity field, such that and P is the exact distribution. Then it is readily verified that the model equation: has the Gaussian solution However, it is important to note, and is also readily verified, that a more general form of the operator L 0 , which is given by where H(X i ) is an arbitrarily chosen well behaved function, also yields the same Gaussian solution for the zero-order equation: Hence at this stage the operator L 0 is not fully determined. Edwards [3] was guided by an analogy with the theory of Brownian motion and in effect made the choice in order to generate a base operator which could be inverted in terms of an eigenfunction expansion of Hermite polynomials. In this process, the {ω i } appeared as eigenvalues.

Evaluation of the coefficients P 0 , P 1 , P 2 ,
The practical problem now is how to actually perform the operations needed to evaluate the coefficients as given by equations (4.17) and (4.18); or, alternatively from equations (4.20) and (4.21). In short, we need to be able to invert the operator L t,0 . We should remind ourselves at this point that Edwards worked with the invariant exact distribution P[u], so that he was able to set ∂/∂t = 0 in the Liouville equation (4.2). Hence, for the case studied by Edwards, we have (4.31) Thus, in the Edwards theory, the first-order coefficient given by (4.17) becomes As mentioned earlier, this inversion was performed using an expansion in Hermite polynomials: for details see [3,14] or [12]. At lowest nontrivial order the resulting equation for the covariance was This result may be found as equation (7.60) in the book by Leslie [12]. The operator M ijm is just a compact form of the usual inertial transfer operator as defined by equation (2.2). Note that Leslie followed Herring [8] in including a negative damping term ω i to sustain the turbulence, but we do not need that here. The calculation can be carried on to higher orders, with increasing numbers of eigenvalues {ω i } appearing in the denominator. Further details can be found in the references cited above. We shall close this brief discussion with a quotation from Leslie (see page 121 of the book [12]). Arguing that Kraichnan's DIA was more general than the Edwards theory, because it allowed other time dependences than just the exponential form, he stated the following: It should be possible to derive alternatives to the Fokker-Planck operator which also permit other time-dependences while retaining the Gaussian form for the basis function P F . So far as we know, the present work is the first attempt to do this. Now in order to explain our present approach, we return to equation (4.15) for P 1 . Substituting for L t,0 , and re-arranging gives If we choose to neglect the second term on the right hand side, then our starting point for the iterative calculation is the assumption That is, inversion of L t,0 is taken to be just a time integral. Evidently the second line of (4. 35) suggests that this step can be analysed in terms of the commutator of the operators L 0 and Θ 1 , with the freedom to choose the functional H [u] being exploited to ensure that this commutator is no lower in order than 2 . Also, for the iteration to be carried on to higher orders, we would need the additional condition for all values of n. However, we will not pursue that here, but instead will explore the consequences of our basic assumption for the calculation at lowest nontrivial order.
With the preceding assumptions, equation (4.15) allows us to calculate the leading order contribution to odd moments of the velocity field. It can be rewritten as which is then integrated to give a form for the first-order coefficient P 1 , thus: (4.39) Similarly we may calculate the second-order coefficient P 2 from equation (4.16) as where we have used the condition L 0 P 0 = 0. Therefore our second-order approximation to the full pdf P[u, t] is non-Gaussian, with the first-order antisymmetric coefficient P 1 [u, t] expressed in terms of an operator acting on the Gaussian zero-order distribution. In turn, the second-order symmetric coefficent P 2 [u, t] is expressed in terms of operators acting on P 1 and P 0 ; and, as we have seen earlier, can be further reduced to operators acting on the zero-order distribution only.

Evaluation of the triple moment
We now confront the full closure problem in the form of the two-time covariance equation, given in this work as (2.23). Replacing the angle brackets by expectation values obtained by averaging against the exact pdf P[u, t], we have: (4.41) Then, substituting the perturbation expansion for P [u, t] into both sides leads to an approximate equation at lowest non-trivial order, thus: (4.42) which, from the self-consistency of the model and the symmetry of the odd-order coefficients, holds to order 2 . It is at this point that the Edwards method uses symmetry to compensate for the lack of a variational principle, such as the Bogoliubov principle in the theory of critical phenomena (e.g. see section 7.5 of the book [68]).
Where convenient, we can restore the angle-bracket notation in modified form as: (4.43) Here · · · 0 indicates that the average has been taken against P 0 ; and similarly · · · 1 is the average against P 1 . The left hand side of (4.43) yields the covariance tensor. Taking the trace of both sides, and invoking isotropy, leads to the form ∂ ∂t + νk 2 C(k; t, t ) = P(k; t, t ), which we had previously as (2.24), but where now 8 (4.44) We recall that repeated tensor indices are summed. Substituting the expression for P 1 [u, t] as given by equation (4.39) yields 8 There should be no confusion between the symbol P(k; t, t ) and the use of P to represent the pdf.
where, in going to the last line, we performed an integration by parts with respect to the velocity field, and dropped the resulting boundary terms as they are zero, such that the derivative now acts on velocity components involved in the triple moment. It may be helpful at this point to recall that Fourier coefficients of the velocity belonging to different wavevectors are independent variables. Accordingly, functional derivatives evaluated at the same time give rise to delta functions, in the usual way. This is not the case for nonsimultaneous functional derivatives, as the Fourier coefficients are correlated in time. In that case we can replace the non-simultaneous functional derivatives in terms of the infinitesimal velocity-field propagator, as defined by equation (2.8). Then, using the model distribution to perform the average, and acting on each of the three velocity components in turn, we find Following Kraichnan, we assume that the instantaneous infinitesimal response function can be averaged independently of the velocity field. This was described as a weak dependence principle: see page 507 of the article by Kraichnan [23] or page 371 of the book by Beran [11]. In effect it is also a mean-field approximation (see the books by Leslie [12] or McComb [14]), and allows us to factor the averages as (4.47) Since the model distribution is Gaussian, we can then factor the fourth-order moment into products of covariances in the usual way. Recalling that the correlation tensor takes the isotropic form as given by equation (2.5), we note that pairings of the velocities such as q and p − q violate the triangle condition, since they give rise to δ(p) which forces the vertex operator M(p) = M(0) = 0. We first evaluate I αβγ : (4.48) where in the last line the factor 2 comes from the symmetry M ρµν (j) = M ρνµ (j). The terms J αβγ and K αβγ are evaluated in a similar way to give (4.50) These results, along with equation (4.46), complete our calculation of P(k; t, t ).

The equations for the covariance
We substitute from equation (2.9) for the ensemble-averaged infinitesimal propagator into (4.46), along with equations (4.48)-(4.50), and invoke isotropy to write are well known (e.g. see [14] or [17]) and satisfy the relationship Using this relationship to replace A(k, j, k − j), along with the change of dummy variables j → k − j for the two terms which now contain L(k, k − j), we arrive at the two-time covariance equation, We may also derive an equivalent equation for C(k; t, t), thus: Equations (4.55) for C(k; t, t ) and (4.56) for C(k; t, t) are identical in form to the Kraichnan-Wyld-Edwards covariance equations [17]; but with one crucial difference. Here the response of the system is not determined by the variation of the velocity field due to an infinitesimal change in the stirring force. Instead the response function arises quite naturally when we carry out the functional differentiation with respect to changes in the velocity field. Note that we have not added an input term to (4.56), but this can be done at any stage, just as in equation (2.25). We should note in passing that our formulation of these covariance equations is different from that used by Kraichnan [23] and was originally due to Edwards [3] who introduced the L(k, j) coefficients. Also, it is a simple matter to recover the Edwards theory from (4.56) and this may be found in appendix.

Equation for the response function
From equation (2.9) we have the definition of the response function. We adopt a simplified, temporary notation in which only the time arguments appear explicitly and re-write this as: which may be further written as the last step following by partial integration and the assumption that the boundary terms vanish as P 0 → 0 as u → ∞. Restoring full notation, our expression for the response function is therefore: , (4.59) and evaluation of this only requires our explicit form for P 0 . It is worth noting that there are some points of resemblance between this response function and the analogous form g(k; t, t ) in Herring's time-dependent self-consistent field theory. This latter is defined by where U(k; t, t ) satisfies the integral equation for t > t and zero otherwise, and L(t, t ) is a particular form of Liouvillian which arises in Herring's theory [9].
Let us assume a general (exponential) form for the base distribution, thus: (4.62) where N is the normalization and S has to be specified. If this were a problem in microscopic physics at equilibrium, then S would be the thermodynamic entropy (divided by the Boltzmann constant). However it has been argued by Edwards and McComb [69], that the entropy is available as a concept in the statistical theory of turbulence, provided it is defined in terms of the information. Accordingly, we may then introduce a 'force like' object f , which it is natural to refer to as a quasi-entropic force, thus: (4.63) Hence from (4.59) we may express the response function in terms of f as It may be noted that the tilde distinguishes f from the stirring force f. Edwards showed that uf was the rate of doing work by the stirring forces on the velocity field, whereas the new quantity uf 0 determines the two-time response. Carnevale et al [70] showed that an expression of this form leads to a linear relationship if the initial distribution is Gaussian. In our case, P 0 [u] is the zero-order distribution which is Gaussian for all times t t. Accordingly, putting we havef and hence (4.67) Then, multiplying across, and with wavenumber labels restored, our final result is: This is the general form of the fluctuation-response relation and, taken with the Kraichnan-Wyld-Edwards equations for C(k; t, t ) and C(k; t, t), as derived earlier, constitutes the LET theory, as previously obtained by more ad hoc methods: see [59] and references therein.

Discussion
Three topics seem to merit further discussion. First there is the question of the relevance of a fluctuation-response relation to turbulence, a subject which has received a certain amount of attention over the years. Then there is the question of how the local energy transfer (LET) theory relates to the earlier Eulerian theories. This might be subsumed under the heading: how relevant is the infra-red divergence? And lastly there is the rather open question of practical applications. We will discuss these three topics in turn.

Fluctuation-response relations (FRR)
If we apply the fluctuation-response relation (4.68) to a microscopic system in thermal equilibrium (e.g. Brownian motion), then it is a simple matter to show that it can be reduced to the well known fluctuation-dissipation theorem. We refer the interested reader to the classic review by Kubo [71] for further information on its relevance to microscopic problems. Here we concentrate on the relevance of this kind of approach to fluid turbulence. It is for this reason that we follow the recent trend and refer to it in this way as the fluctuation-response relation (FRR), in order to keep in mind that we are referring to the more general form.
It should be mentioned that we originally introduced the relationship (4.68) as a propagator relationship, with some heuristic justification coming from a consideration of Wyld's analysis [4]. However, one can interpret theories in terms of such a relationship. For instance, as pointed out by Frederiksen and Davies [72], both Herrings's SCF theory and the LET theory can by related to Kraichnan's DIA, by means of the fluctuation-response relation. In the case of Herring's SCF, the single-time covariance and response function are calculated from the DIA, and the FRR is used to calculate the two-time covariance. Whereas, for LET, the DIA covariance equations 9 are used to calculate C(k; t, t ) and C(k; t , t ), and the FRR is then used to fix R(k; t, t ).
This is not the same as simply assuming the validity of the fluctuation-response relation. Although that too has a history, as in the past it has proved popular in atmospheric turbulence and climate research. The seminal work of Leith [73] and Bell [74], was inspired by Kraichnan's demonstration that the FRR could be applied to classical nonlinear systems provided that they were in thermal equilibrium [75]. A recent discussion of the use of the FRR in climate modelling can be found in the article by Zidikheri and Frederiksen [76], along with many references.
However, from our present point of view, we note that the subject of fluctuation-response relations has also received much attention from the point of view of dynamical systems theory over the last two decades. A significant development has been a generalisation of Kraichnan's result [75] to macroscopic systems that are both chaotic and mixing [70]. Under these circumstances the system response can be expressed in terms of its invariant probability measure at the initial time. Further, Carnevale et al [70] have shown that if this invariant measure is Gaussian, then the linear result equivalent to equation (4.68) is recovered. The interested reader is also referred to the paper by Boffetta et al [77] and references therein.
Our present work can be compared with that of Carnevale et al [70]. For instance, our equation (4.59) is essentially of the same form as their equation (5). The main difference is that they average over trajectories against an initial distribution, whereas we average over our model (or zero-order) pdf which is valid for all times. Bearing in mind that the coefficients in the expansion of the exact pdf in our case are evaluated perturbatively against the zero-order distribution, it is clear that in this context the FRR applies at all times and to all orders of perturbation theory.
As ours is a continuum (statistical field theoretic) problem, which is why our analysis is in terms of functionals, we do not have any concerns about the differentiability of distributions. In contrast, when averaging over trajectories, Carnevale et al [70] invoke the role of noise in smoothing the probability distribution of trajectories (see also Boffetta et al [77] for a further discussion of this point). Interestingly, the Edwards theory [3] had a similar step, in that the pdf of the velocities was obtained by averaging against the distribution of the stirring forces. But arguably this is only necessary to facilitate the many-body theoretical derivation of the Liouville equation [3], rather than to impose any additional smoothness. However, irrespective 9 Now referred to by us as the Kraichnan-Wyld-Edwards or KWE covariance equations. of such considerations, when Carnevale et al [70] consider their initial invariant distribution to be Gaussian, the result is the linear relationship which is their equation (21), while our use of the Gaussian zero-order distribution ensures that our equation (4.68) takes the same form.

How relevant is the infra-red divergence?
The analysis by Edwards of the failure of his SCF theory to yield the K41 spectrum, as discussed in section 3.2.2, is both elegant and rigorous. It is also of its time, and is mainly of historical significance. Yet the associated idea of the infra-red divergence as an index of problems with statistical closures has taken root, and is still often quoted. This can lead to misunderstandings and we will address that aspect here.
First we should clear one possible source of confusion out of the way. In later work on the application of renormalization group methods to stirred fluid motion, references to infra-red and ultra-violet divergences started to occur. It should be understood that these divergencs have nothing to do with the infinite Reynolds number limit of Edwards, but are entirely to do with the arbitrarily chosen stirring forces. For instance, where Edwards chose an input term involving the delta function δ(k) in order to give the Kolmogorov spectrum, later workers chose the power law k −3 for the same purpose; and obviously this latter form could lead to both types of divergence. A general discussion of this topic has been given elsewhere [78].
The analysis by Edwards showed that his covariance equation was well behaved at the possible points of divergence, due to cancellations of pairs of individually divergent terms. It was the lack of such cancellations that was the problem with his response equation. This analysis was so transparent and easy to follow, that there was a tendency to use it for pedagogic reasons when analysing other theories. This was particularly true of the book by Leslie [12], which is still influential in this subject.
However, it needs to be clearly understood that in applying the Edwards analysis to a fully two-time theory such as DIA or LET, approximations have to be made. In particular, the assumption of exponential time-dependences; and the further assumption of the same decay time for both correlation and response functions, are both incorrect. This too was discussed by Leslie. It is therefore the case that, although the concept of the infra-red divergence in the limit of infinite Reynolds numbers can offer some guidance and possible insight, it is not in itself a prescriptive test for the more general theories.
Having said that, we should nevertheless note that the assumption of specific time-depend ences can still tell us something about the structure of theories like DIA and LET, even if this may be seen as heuristic in a mathematical sense. The reason for this is that a oneto-one correspondence is maintained between terms in the two formulations. This fact was exploited in the original derivation of the two-time form of the LET theory [53].
Also, while on this topic, we may mention that some aspects of the work by Edwards and McComb [69] foreshadowed the later local energy transfer theory [7]. In particular they noted that the avoidance of the infra-red divergence in the covariance equation relied on cancellations of pairs of terms making up T(k). Hence at j = 0 and |k − j| = 0 these terms can not be treated as independent modes.

Potential for practical applications
In view of the great importance of turbulent flow in engineering applications, it is suprising that the two-point statistical closures have been so little exploited for practical applications.
To our knowledge, there have only been two exceptions to this. First there is the activity on climate and atmospheric modelling by Frederiksen and co-workers (see [79] for a recent account), who have used DIA, Herring's SCF and LET theories over several decades for this purpose. This work also includes sub-grid modelling.
Secondly, there is the eddy-damped quasi-normal Markovian (EDQNM) theory. This arose from Orszag's analysis of the failure of the original quasi-normality hypothesis and reflects a later and better understanding of turbulence phenomenology. The EDQNM theory is really a model (compare Kraichnan's test field model [80]), in that it makes a particular assumption which requires the fitting of an adjustable constant. In stationary form it reduces to the Edwards SCF theory. Because it is a single-time theory, it is relatively easy to compute and has been widely applied to practical problems (e.g. [81][82][83][84]). For a recent discussion see the book by Sagaut and Cambon [15].
In contrast, the LET theory, consisting as it does of the KWE equations (4.55) and (4.56), plus the fluctuation-relaxation relation (4.68), is a fully two-time theory and has no adjustable parameters. Although rather more complicated than the EDQNM theory, and more demanding to compute, it nevertheless could be applied to non-isotropic and even inhomogeneous problems. Even if presenting too formidable a computational challenge, it could provide a basis for single-point modelling, which might offer a useful extension to existing approaches.

Conclusion
The LET theory consists of the two covariance equations: (4.55), for t = t ; and (4.56), for t = t ; along with an equation for the response function. This can take the form of either the velocity propagator relation as given by equation (3.18); or the relation derived from it as equation (3.20), which may be compared with the DIA form as given by (3.19). The two forms of LET are equivalent. On the one hand, the relationship (3.20) has pedagogic value in showing clearly the relationship between LET and DIA, and in particular showing that the LET possesses the non-Markovian structure which is implied by experimental results. While, on the other hand, the propagator relationship (3.18) is much easier to compute, and in fact its use makes the LET easier to compute than the DIA.
Originally, the use of the propagator relation (3.18) was hypothesised from a consideration of the Wyld-Lee formalism (see Berera et al [6]). It has now been derived and identified formally as a fluctuation-response relation. In deriving this, we have made three approximations, which we list as follows: The LET theory has been computed for the test problems of isotropic turbulence (including for diffusion of a passive scalar) with good results, both qualitative and quantitative (see section 3.4 of the present paper). It is arguable that it now has a formal derivation which puts it on the same footing as the other two-time theories, particularly Kraichnan's DIA and Herring's SCF. In view of the fact that it has the unique characteristic of a relationship between the transfer spectrum and the system response which is qualitatively in accord with experiment, one might hope that this will stimulate interest in this class of approach to the turbulence problem.

Acknowledgments
SRY acknowledges the support of the UK EPSRC (grant number EP/N028694/1) and the EC's H2020 LASERLAB-EUROPE (grant number 654148).

Appendix. Recovering the Edwards SCF theory
Edwards derived his closure approximation by inverting the Fokker-Planck form of L 0 in terms of an eigenfunction expansion, with the inverse modal lifetimes appearing as eigenvalues. However, instead of this, it is a simple matter to recover the Edwards SCF theory directly from equation (4.56).
Our starting point for this, is the recognition by Kraichnan in 1964 [85] that the DIA and Edwards covariance equations were equivalent if one assumed the same exponential time dependence for the covariance and the response function. That is, where ω(k) is the inverse lifetime for correlation and response of mode k. Its prescription is the second part of the Edwards closure approximation. Later, Kraichnan showed that the Edwards theory could be extended to statistically non-stationary states [86], while Herring and Kraichnan [62] demonstrated that the generalized Edwards covariance equation could be obtained from (4.56), with C(k, t) ≡ C(k; t, t), as: where they introduced a memory function D, given by For the stationary case, this could be written [85]: D(k, j, |k − j|) = 1 ω(k) + ω( j) + ω(|k − j|) . (A.5)