% The paper is formatted in LaTex Version 2.09. % Figures are available upon request from the author. % However, the first four figures are identical with % the first four figures in ref.[39], Phys. Fluids, % vol.6, 3063 (1994) and can be obtained there. % The macro "seceq.sty" is included in the following % file, which numbers equations sequentially by section. % It should be included in the directory when formatting % the LaTex file for the printer. If numbering by section % is not desired, then the "seceq" may simply be omitted % from the "documentstyle" command below. % %%%%%%%%%%%%%%%%%%%%%%%%% %% seceq.sty macro %% %%%%%%%%%%%%%%%%%%%%%%%%% % % ripped off by jonathan from euthesis.sty % % numbers equations consecutively within sections % **************************************** % * MISCELLANEOUS * % **************************************** % % EQUATION and EQNARRAY -- put here because it must follow \chapter definition % % \newcounter{equation} % \@addtoreset{equation}{section} % Makes \chapter reset 'equation' counter. \def\theequation{\thesection.\arabic{equation}} % \jot = 3pt % Extra space added between lines of an eqnarray environment % The macro \@eqnnum defines how equation numbers are to appear in equations. % \def\@eqnnum{(\theequation)} %\def\@eqnnum{{\rm (\thesection .\theequation)}} % BODY \documentstyle[11pt,seceq]{article} \newcommand{\be}{\begin{equation}} \newcommand{\ee}{\end{equation}} \newcommand{\lb}{\label} \newcommand{\en}{\epsilon} \newcommand{\ven}{\varepsilon} \newtheorem{Th}{Theorem} \newtheorem{Lem}{Lemma} \newtheorem{Conj}{Conjecture} \newtheorem{Prop}{Proposition} \newtheorem{Ex}{Example} \newcommand{\bv}{{\bf v}} \newcommand{\br}{{\bf r}} \newcommand{\bk}{{\bf k}} \newcommand{\bl}{{\bf l}} \newcommand{\bh}{{\bf h}} \newcommand{\bj}{{\bf j}} \newcommand{\bm}{{\bf m}} \newcommand{\bn}{{\bf n}} \newcommand{\bpd}{{\bf \nabla}} \newcommand{\bE}{\hat{{\bf e}}} \newcommand{\bz}{{\bf 0}} \newcommand{\bR}{{\bf R}} \newcommand{\bP}{{\bf P}} \newcommand{\bp}{{\bf p}} \newcommand{\bD}{{\bf D}} \newcommand{\bZ}{{\bf Z}} \newcommand{\bV}{{\bf V}} \newcommand{\bK}{{\bf K}} \newcommand{\bT}{{\bf T}} \newcommand{\bA}{{\bf A}} \newcommand{\bU}{{\bf U}} \newcommand{\bG}{{\bf G}} \newcommand{\ef}{{\rm eff}} \newcommand{\bC}{{\bf C}} \newcommand{\bB}{{\bf B}} \newcommand{\bL}{{\bf L}} \newcommand{\bI}{{\bf I}} \newcommand{\med}{{\rm med}} \newcommand{\bu}{{\bf u}} \newcommand{\bF}{{\bf f}} \newcommand{\fs}{{\bf f}^\prime} \newcommand{\fl}{\overline{{\bf f}}} \newcommand{\fS}{f^\prime} \newcommand{\fL}{\overline{f}} \newcommand{\Fs}{{\bf F}^\prime} \newcommand{\Fl}{\overline{{\bf F}}} \newcommand{\BF}{{\bf F}} \newcommand{\ba}{{\bf a}} \newcommand{\bs}{{\bf s}} \newcommand{\hX}{{\hat{X}}} \newcommand{\hP}{{\hat{P}}} \newcommand{\hH}{{\hat{H}}} \newcommand{\hL}{{\hat{L}}} \newcommand{\hbX}{{\hat{{\bf X}}}} \newcommand{\hbP}{{\hat{{\bf P}}}} \newcommand{\bb}{{\bf b}} \newcommand{\bc}{{\bf c}} \newcommand{\bx}{{\bf x}} \newcommand{\bX}{{\bf X}} \newcommand{\by}{{\bf y}} \newcommand{\bw}{{\bf w}} \newcommand{\grad}{{\mbox{\boldmath $\nabla$}}} \newcommand{\bun}{{\mbox{\boldmath $1$}}} \newcommand{\st}{$\,_{*}$} \newcommand{\vl}{\overline{{\bf v}}} \newcommand{\cvl}{\overline{{\bf V}}} \newcommand{\taus}{{\tau^r}} \newcommand{\taul}{{{\tau^s}}} \newcommand{\vs}{{\bf v}^{\prime}} \newcommand{\ws}{{\bf w}^{\prime}} \newcommand{\Vs}{{\bf V}^{\prime}} \newcommand{\Ws}{{\bf W}^{\prime}} \newcommand{\ul}{\overline{{\bf u}}} \newcommand{\Kl}{\overline{{\bf K}}} \newcommand{\us}{{\bf u}^{\prime}} \newcommand{\Ks}{{\bf K}^{\prime}} \newcommand{\pl}{\overline{p}} \newcommand{\ps}{p^\prime} \newcommand{\zl}{\overline{z}} \newcommand{\zs}{z^\prime} \newcommand{\Zl}{\overline{Z}} \newcommand{\Zs}{Z^\prime} \newcommand{\Gl}{\overline{G}} \newcommand{\Gs}{G^\prime} \newcommand{\Bl}{\overline{b}} \newcommand{\Bs}{b^\prime} \newcommand{\BL}{\overline{B}} \newcommand{\BS}{B^\prime} \newcommand{\FL}{\overline{F}} \newcommand{\FS}{F^\prime} \newcommand{\Sigmal}{\overline{\Sigma}} \newcommand{\Sigmas}{\Sigma^\prime} \newcommand{\phil}{\overline{\phi}} \newcommand{\Phil}{\overline{\Phi}} \newcommand{\Phis}{\Phi^\prime} \newcommand{\el}{\overline{e}} \newcommand{\btau}{{\mbox{\boldmath $\tau$}}} \newcommand{\bdot}{{\mbox{\boldmath $\cdot$}}} \newcommand{\btimes}{{\mbox{\boldmath $\times$}}} \newcommand{\bleta}{{\mbox{\boldmath $\eta$}}} \newcommand{\btaul}{{\mbox{\boldmath $\tau$}}^s} \newcommand{\btaus}{{\mbox{\boldmath $\tau$}}^r} \newcommand{\bphi}{{\mbox{\boldmath $\phi$}}} \newcommand{\ben}{{\mbox{\boldmath $\epsilon$}}} \newcommand{\bsigma}{\overline{\mbox{\boldmath $\sigma$}}} \textwidth6.25in \textheight8.5in \oddsidemargin.25in \topmargin0in \renewcommand{\baselinestretch}{1.7} \begin{document} \title{Turbulence Noise} \author{Gregory L. Eyink\\{\em Department of Mathematics, University of Arizona}\\ {\em Building No. 89, Tuscon, AZ, 85721}} \date{ } \maketitle \begin{abstract} We show that the large-eddy motions in turbulent fluid flow obey a modified hydrodynamic equation with a stochastic turbulent stress whose distribution is a causal functional of the large-scale velocity field itself. We do so by means of an exact procedure of ``statistical filtering'' of the Navier-Stokes equations, which formally solves the closure problem, and we discuss relations of our analysis with the ``decimation theory'' of Kraichnan. We show that the statistical filtering procedure can be formulated using field-theoretic path-integral methods within the Martin-Siggia-Rose formalism for classical statistical dynamics. We also establish within the MSR formalism a ``least-effective-action principle'' for mean turbulent velocity profiles, which generalizes Onsager's principle of least dissipation. This minimum principle is a consequence of a simple realizability inequality and therefore holds also in any realizable closure. Symanzik's theorem in field-theory---which characterizes the static effective action as the minimum expected value of the quantum Hamiltonian over all state vectors with prescribed expectations of fields---is extended to MSR theory with non-Hermitian Hamiltonian. This allows stationary mean velocity profiles and other turbulence statistics to be calculated variationally by a Rayleigh-Ritz procedure. Finally, we develop approximations of the exact Langevin equations for large eddies, e.g. a random-coupling DIA model, which yield new stochastic LES models. These are compared with stochastic subgrid modelling schemes proposed by Rose, Chasnov, Leith, and others, and various applications are discussed. \newpage {\em Key words:} Navier-Stokes, turbulence noise, LES, generalized Langevin equations, variational principles {\em PACS numbers:} 47.27.Sd, 03.40.Gc, 47.25.Cg, 64.60.Ak \end{abstract} \section{Introduction} The analogy of turbulent motion with the hydrodynamics of molecular fluids has played a central role in attempts to understand and to model the dynamics of large, turbulent eddies. In ordinary Newtonian hydrodynamics, at least in the low Mach number, incompressible regime, the local fluid velocity obeys an equation of the form \be \partial_t \bv+\grad\cdot(\bv\bv+\btaul)=-\grad p, \lb{1} \ee in which $p$ is the (kinematic) pressure, required to enforce the incompressibility constraint $\grad\cdot\bv=0,$ and $\btaul$ is the {\em viscous stress tensor} \be \btaul= -\nu_0[(\grad\bv)+(\grad\bv)^\top]. \lb{2} \ee The so-called {\em molecular viscosity} $\nu_0$ represents the residual, dissipative effects of the ``graininess'' of matter at the fluid level. Likewise, in the case of turbulent motion, the ``large-scale velocity'' $\bv_\ell,$ which may be conveniently defined by a smooth filtering \be \vl_\ell(\br)\equiv \int d^d\br'\,\,G_\ell(\br-\br')\bv(\br'), \lb{3} \ee obeys an equation of the form \be \partial_t \vl_\ell+\grad\cdot(\vl_\ell\vl_\ell+\btau_\ell)=-\grad \pl_\ell, \lb{4} \ee in which $\pl_\ell$ is the filtered pressure field and $\btau_\ell$ is a tensor representing the {\em turbulent stress} of the eliminated small-scale eddies $<\ell.$ If the Reynolds number of the fluid is high and $\ell$ is chosen in the long ``inertial range'' of scales, then the molecular viscosity is believed to play a negligible role in the evolution of the large eddies $>\ell.$ Instead the primary effect is believed to be due to the smaller scale turbulence, which is modelled, by analogy with the molecular fluids, as \be \btau_\ell= -\nu_\ell[(\grad\vl_\ell)+(\grad\vl_\ell)^\top], \lb{5} \ee in which $\nu_\ell$ is now a so-called {\em eddy viscosity}. This representation of the effect of small-scales as a simple damping is, however, not nearly so accurate as in the molecular case, where there is generally a large separation of scales between the atomic and fluid degrees of freedom. In contrast, in the case of turbulence, there is no such scale-separation and the eddy-viscosity representation is flawed, leading to a variety of viscoelastic, nonlinear effects \cite{Kr76,Kr87}. Here we shall be concerned with another over-simplification of the eddy-viscosity representation. In fact, even in the case of molecular fluids there is an additional influence of the microscopic degrees of freedom, which is a stochastic effect due to the chaotic molecular motions. Because of such effects, Landau and Lifshitz \cite{LL} proposed in 1959 that there should be also in the Navier-Stokes fluid equation Eq.(1) a {\em random stress} $\btaus,$ as \be \partial_t \bv+\grad\cdot(\bv\bv+\btaul+\btaus)=-\grad p. \lb{6} \ee Since the molecular degrees of the fluid are, by mixing properties of the dynamics, locally in a state of thermal equilibrium, the random stress field $\btaus(\br,t)$ should have statistics which are those of the universal Gibbs states. Motivated by such considerations, Landau and Lifshitz hypothesized that $\btaus$ for the fluid at temperature $T$ should be a Gaussian random field with zero mean and covariance \be \langle \tau^r_{ij}(\bx,t)\tau^r_{kl}(\by,s)\rangle=(2k_BT\nu_0/\rho)(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk}) \delta^d(\bx-\by)\delta(t-s), \lb{7} \ee where the delta-functions arise from the fast decay of correlations in the Gibbs states. This ``fluctuation-dissipation relation'' (FDR) relating the strength of the noisy stresses and the molecular viscosity makes clear that the two contributions to the momentum flux have the same origin in the microscopic degrees-of-freedom. Generally speaking, the molecular noise gives a negligible effect on the gross fluid motions, but it does have observable consequences for fluctuation behavior seen, for example, in light-scattering experiments. By the same analogy with molecular fluids which motivated the eddy-viscosity model Eq.(\ref{5}), one therefore should expect that the stress tensor in the equation of motion of large turbulent eddies will consist of both a systematic (or, deterministic) and a fluctuating part: \be \partial_t \vl_\ell+\grad\cdot(\vl_\ell\vl_\ell+\btaul_\ell+\btaus_\ell)=-\grad \pl_\ell. \lb{8} \ee In fact, an exact formula for the turbulent stress can be deduced from the filtering procedure, as \be \btau_\ell=\overline{(\bv\bv)}_\ell-\vl_\ell\vl_\ell. \lb{9} \ee This representation of the stress is as a function of the {\em full} velocity field $\bv$, that is, of the small-scale (or, sub-grid) component $\vs_\ell,$ defined as \be \vs_\ell(\br)\equiv \bv(\br)-\vl_\ell(\br). \lb{10} \ee as well as of the large-scale velocity $\vl_\ell$ itself. Therefore, one should expect that, when only the large-scale field $\vl_\ell$ is specified, the stress $\btau_\ell$ will not be fully determined but that, in fact, a distribution of values will occur with some given frequency depending upon the precise values of the small-scale velocity $\vs_\ell.$ Unlike the case of laminar fluids where molecular noise represents only a small perturbation to deterministic motion, in turbulent flow the fluctuations are large and the importance of the random component of the stress, or {\em eddy noise}, is likely to be fully as significant as the systematic part. On the other hand, if the statistics of the small-scale component $\vs_\ell$ are {\em universal}, as usually believed, then the distribution law of the fluctuating stress ought to be fixed by the Kolmogorov ``universal equilibrium state'' when $\ell$ is chosen in the inertial-interval. Nothing so simple as the FDR in Eq.(\ref{7}) fixing the molecular noise can be expected, again due to lack of scale-separation, but it ought to be possible to construct a generally applicable model of the statistical distribution. The recognition of the probable importance of ``eddy noise'' in fully-developed turbulence does not originate with this work but was pointed out already in the first cited paper of Kraichnan \cite{Kr76} and in the contemporary work of Rose \cite{Rose}. More recently, a number of stochastic large-eddy models have been proposed to account for such effects \cite{Lei,Chas,MT,Car}. Our purpose here is to set up a general formalism for {\em statistical filtering} of the Navier-Stokes turbulence dynamics, leading to exact ``fluctuating hydrodynamic equations'' for the large eddies of the form proposed. We should point out that in the present method no finite ``stirring forces'' are added, as in the so-called ``RNG'' method of Yakhot and Orszag \cite{YO}. Instead, the ``eddy noise'' we study must be dynamically generated by the chaotic fluid motion of the small-scale eddies. In fact, we believe that a careful analysis of the existence and properties of this turbulent noise must be made {\em before} any modelling is attempted, for otherwise serious confusions in the physics will result. This work develops the framework for such an analysis. Our procedure is most similar to that used by Rose \cite{Rose}, but no iteration schemes motivated by renormalization group ideas will be employed (since we have not found them to present any significant advantages). Although our stochastic equations are exact results for the Navier-Stokes dynamics, they are rather formal and not directly useful for practical modelling. However, subsequent approximations based upon standard ideas (random-coupling model, multiple-scale expansion and/or Markovianization, functional reversion techniques, etc.) yield simplified forms which are suitable for numerical solution by computer. The precise contents of this work are as follows: in the following Section 2 we shall set up the ``statistical filtering'' scheme leading to stochastic equations for the large turbulent eddies and discuss its statistical-dynamical foundations. We compare there as well the relations of this method with Kraichnan's ``decimation theory'' \cite{Kr85}. His theory leads also to generalized Langevin equations by a procedure of elimination or ``weeding'' of modes that are then modelled by Langevin forces subject to realizability constraints. Afterward we develop from first principles the Martin-Siggia-Rose (MSR) field theory formulation of classical statistical dynamics \cite{MSR}, emphasizing the non-perturbative and exact basis of the approach. The MSR formalism is shown to yield an elegant and efficient formulation of the statistical filtering scheme. In the next Section 3 we discuss a rather different---but closely related---subject of variational principles for mean velocity profiles based upon the MSR effective action. It is shown that a principle of ``least effective-action'' characterizes the mean velocity as a consequence of a simple realizability inequality and thereby generalizes the famous least-dissipation principle of Onsager \cite{O} for most probable states in nonequilibrium thermodynamics. The least-action principle is made the basis to calculate mean turbulent profiles and other turbulence statistics by variational strategies rather than by the solution of effective equations of motion. Finally, in the last Section 4 we briefy describe approximations of the exact large-eddy equations, more practical for LES, and we compare these models with those which have been introduced previously \cite{Lei,Chas,MT,Car}. We also briefly discuss applications of such equations to the problem of atmospheric predictability, where eddy noise can play a significant role in LES, and to simulation of inhomogeneous turbulence in a complex geometry. \section{The Statistical Filtering Method and MSR Field Theory} \noindent{\em (2.1) Formulation of Statistical Filtering} Let us first recall the usual (deterministic) filtering method, of which an interesting recent account can be found in the paper of Germano \cite{Ger1} (see also our paper \cite{I-LES}). It was already briefly reviewed in the Introduction. The basic idea of the method is to convolute the turbulent velocity field with a smooth ``filter function'' $G_\ell(\br)=\ell^{-d}G(\br/\ell),$ as in Eq.(\ref{3}) there, in order to define a ``large-scale'' velocity $\vl_\ell.$ It is then possible to derive an infinite hierarchy of equations for ``generalized central moments,'' of which the turbulent stress $\btau_\ell$ in Eq.(\ref{9}) is the 2nd order example. The lowest equation in the hierarchy is precisely the Eq.(\ref{4}) in the Introduction. If the hierarchy is truncated at this order, then one encounters the ``closure problem'' that $\btau_\ell$ is not a functional of $\vl_\ell$ alone but also of $\vs_\ell$ and a similar problem occurs at every higher-order truncation of the hierarchy as well. In the current large-eddy simulation (LES) literature, second-order modelling is not attempted and only the lowest equation in the hierarchy is retained. Hence a model representation of $\btau_\ell$ must be sought as a (causal) functional of $\vl_\ell$ in order to obtain an autonomous equation for this field. We shall show that a ``stochastic filtering'' method solves the closure problem in principle (but not in practice) by yielding an {\em exact} representation of $\btau_\ell$ as \be \btau_\ell[\vl_\ell]=\btaul_\ell[\vl_\ell]+\btaus_\ell[\vl_\ell], \lb{41} \ee in which $\btaul_\ell[\vl_\ell]$ is a ``systematic stress,'' some causal functional of $\vl_\ell,$ and $\btaus_\ell[\bv_\ell]$ is a Langevin force, or ``random stress,'' whose distribution given the past history of $\vl_\ell$ is specified. We shall refer to this exact Langevin dynamical equation as the ``stochastic large-eddy equation'' (SLE). Although the functional $\btaul_\ell[\vl_\ell]$ and the distribution of $\btaus_\ell[\vl_\ell]$ are calculable in principle from formulae given below, in practice the exact expressions cannot be evaluated and approximation schemes must be employed. This is the subject of the last Section 4. Just as for the deterministic filtering procedure, a rather wide class of choices for filter function $G$ may be made. The main requirement is that the Fourier transform $\widehat{G}(\bk)$ must decay sufficiently rapidly as $k\rightarrow+\infty,$ so that the convolution reallys represents an ``elimination'' of high-wavenumber modes. The most common choices are the {\em Gaussian filter}, usually defined in the LES literature as \be G(\bx)=\left({{6}\over{\pi}}\right)^{d/2}\exp\left[-{{6|\bx|^2}\over{\pi}}\right], \lb{42} \ee the {\em tophat filter} \be G(\bx)= \left\{ \begin{array}{ll} 1 & \mbox{ if $\max_i|x_i|<{{1}\over{2}}$} \cr 0 & \mbox{ otherwise,} \end{array} \right. \lb{43} \ee and the {\em sharp Fourier cutoff filter}, which is most easily defined by its Fourier transform, as: \be \widehat{G}(\bk)=\left\{ \begin{array}{ll} 1 & \mbox{ if $\max_i|k_i|<\pi$} \cr 0 & \mbox{ otherwise.} \end{array} \right. \lb{44} \ee In physical space the latter gives \be G(\bx)= 2^d\prod_{i=1}^d {{\sin\pi x_i}\over{x_i}}. \lb{45} \ee However, as discussed at length in our work \cite{I-LES}, the sharp Fourier filter gives rise to pathological features in the deterministic method and these can be seen to appear as well in the stochastic formulation. Therefore, we will always consider a ``graded'' filter obeying the modest condition, Eq.(31) in \cite{I-LES}, which allows us to avoid the bad features of the sharp Fourier cutoff. The first steps in the ``statistical filtering'' method are the same as those in the usual deterministic scheme. The ``large-scale (LS) velocity'' \be \vl_\ell(\br)=(G_\ell*\bv)(\br), \lb{45a} \ee and ``small-scale (SS) velocity'' \be \vs_\ell(\br)=(H_\ell*\bv)(\br), \lb{45b} \ee are introduced by convolution with low-pass filter $G_\ell$ and high-pass filter $H_\ell$ satisfying \be \widehat{G}(\bk)+\widehat{H}(\bk)=1. \lb{45c} \ee The latter guarantees that $\bv=\vl_\ell+\vs_\ell.$ If the filters are applied to the dynamical equations themselves, then there result coupled equations: \be \partial_t \vl_\ell+\grad\cdot(\vl_\ell\vl_\ell+\btau_\ell)=-\grad \pl_\ell+\nu_0\bigtriangleup\vl_\ell, \lb{40} \ee which we call the ``large-eddy equation'' (LE), and \be \partial_t \vs_\ell+\grad\cdot(\vl_\ell\vs_\ell+\vs_\ell\vl_\ell+\vs_\ell\vs_\ell-\btau_\ell)= -\grad \ps_\ell+\nu_0\bigtriangleup\vs_\ell, \lb{40a} \ee which we call the ``small-eddy equation'' (SE). Here $\btau_\ell$ is the ``turbulent stress-tensor,'' which is given by an explicit quadratic function of the total velocity, $\btau_\ell=\bT_\ell(\bv,\bv),$ \be \bT_\ell(\bv,\bv)=\overline{(\bv\bv)}_\ell-\vl_\ell\vl_\ell. \lb{45d} \ee A useful decomposition of this function was made by Germano \cite{Ger2}, according to the following straightforward procedure. One substitutes $\bv=\vl_\ell+\vs_\ell$ into $\bT_\ell$ to obtain: \be \bT_\ell(\bv,\bv)=\bL_\ell(\bv,\bv)+\bC_\ell(\bv,\bv)+\bR_\ell(\bv,\bv), \lb{112} \ee where \be \bL_\ell(\bv,\bv)\equiv \bT_\ell(\vl_\ell,\vl_\ell), \lb{113} \ee is the so-called {\em Leonard stress}, \be \bC_\ell(\bv,\bv)\equiv \bT_\ell(\vl_\ell,\vs_\ell)+\bT_\ell(\vs_\ell,\vl_\ell), \lb{114} \ee is the {\em cross stress}, and \be \bR_\ell(\bv,\bv)\equiv\bT_\ell(\vs_\ell,\vs_\ell), \lb{115} \ee is denoted as {\em Reynolds stress}. Note this is a modification of the traditional decomposition of Leonard \cite{Leo}. It is clear that the LE and SE are coupled equations for the fields $\vl_\ell$ and $\vs_\ell.$ The fundamental new step in ``statistical filtering'' is to add to the SE a weak random force $\fs_\ell$: \be \partial_t \vs_\ell+\grad\cdot(\vl_\ell\vs_\ell+\vs_\ell\vl_\ell+\vs_\ell\vs_\ell-\btau_\ell)= -\grad \ps_\ell+\nu_0\bigtriangleup\vs_\ell+\fs_\ell, \lb{45e} \ee where $\fs_\ell$ has zero mean and spectral support at wavenumbers greater than $\sim 2\pi/\ell.$ For convenience, we choose $\fs_\ell$ Gaussian. A possible choice of the covariance is \be \langle \fs_\ell(\br t)\fs_\ell(\br' t')\rangle= 2\ven_0 H_\ell(\br-\br')\delta(t-t'). \lb{45f} \ee The parameter $\ven_0$, when multiplied by the density $\rho$ of the fluid, has units of energy per time. The role of this force is simply to add some random perturbations to the otherwise deterministic SE and thus $\rho\ven_0$ should be very small compared to the total energy transfer per time to the small-scale motions. In fact, $\ven_0$ will momentarily be taken to zero. First, we formally solve the SE part of the dynamics {\em with the large-scale velocity field considered fixed.} \footnote{ This is similar to Kraichnan's ``clamping method'' in Section 8 of \cite{Kr85}.} The solution may be written as \be \vs_\ell(\br t)=\Vs_\ell[\br t;\vl_\ell,\fs_\ell], \lb{45g} \ee in which $\Vs[\br,t; \vl_\ell,\fs_\ell]$ is a causal functional of time-histories of $\vl_\ell$ and $\fs_\ell,$ i.e a functional depending only upon their values for times $From a more physical perspective, one can consider the random force to represent the molecular noise Eq.(\ref{7}). Just as the molecular viscosity is believed to play no direct role in the dynamics of the large-scale turbulence, so that the limit $\nu_0\rightarrow 0$ is well-defined, likewise the molecular noise is believed to play no essential role in that limit except to select the correct measure. This idea was already proposed some time ago by Ruelle \cite{Rue2}. In the language of field-theory, the bare, molecular noise will be replaced by a ``renormalized'' eddy noise. This is due to the same physics by means of which molecular viscosity is replaced by an effective or ``renormalized'' eddy viscosity. The field-theoretic point of view will be exploited later on to come up with approximation schemes and calculational methods for ``statistical filtering.'' Stochastic descriptions similar to the SLE have been postulated by Hohenberg and Shraiman \cite{HoSh} in the context of general spatiotemporally chaotic systems and, specifically, for the Kuramoto-Sivashinsky equation. In the latter example there is an important conjecture of Yakhot \cite{Yak} that the chaotic behavior generated by small-scale instabilities will produce an ``effective dynamics'' of the large-scales which is just the noisy Burgers equation. See also \cite{Zal,SKJJB}. Our method generalizes and systematizes this earlier work, which was based upon weak-coupling perturbation expansions (without a small parameter!). As we shall see shortly, the ``stochastic filtering'' we have formulated non-perturbatively may also be carried out in weak-coupling expansions and then contains the earlier results. In fact, even for weak-coupling our systematic technique obtains terms missed in the earlier work. The ideas proposed here have also some similarities with the ``decimation theory'' of Kraichnan \cite{Kr85}, but they are really essentially distinct. Kraichnan proposed a general strategy for economical computation of many-mode systems by a procedure in which modes that are redundant due to underlying statistical symmetries are eliminated and replaced by suitable Langevin forces. The latter are constructed to enforce certain statistical constraints imposed by the exact dynamics. In this way efficient numerical approximations might be obtained and, by adding more and more constraints from the true dynamics, convergence of the approximations to the true values might even be found. The two methods have in common the elimination of modes with ``generalized Langevin equations'' as the output. However, the ``statistical filtering'' we have formulated is exact and without any approximation whatsoever. Of course, the output SLE are, at this stage, completely formal and useless as a tool for numerical computations. Approximations must be developed for any practical application. Nevertheless, the idea is not at its basis an attempt to develop {\em statistical approximations}, as is decimation, but is rather an attempt to characterize a physical {\em statistical law}. This is seen most clearly by comparing the two procedures for ``chaotic'' vs. ``integrable'' systems. Decimation theory might be successfully applied to both cases, yielding generalized Langevin models as statistical approximations of the true dynamics. In fact, Kraichnan's first example of his decimation method \cite{Kr85} was for the three-mode dynamical model \be \partial_t x= A_x yz,\,\,\,\,\,\partial_t y= A_y xz\,\,\,\,\,\partial_t z= A_z xy, \lb{45i} \ee with $A_x+A_y+A_z=0.$ This dynamics has two quadratic integrals and is thus integrable by quadratures in terms of elliptic functions, as discussed long ago by Lorenz \cite{Lor}. On the contrary, if the ``stochastic filtering'' were applied to an integrable system, such as the above three-mode model or Burgers' equation, then the SLE for the explicit modes would degenerate, or become deterministic, in the vanishing noise limit $\ven_0\rightarrow 0.$ The chaotic properties of the dynamics are necessary to create a true statistical law as described by the SLE equations. After these somewhat philosophical remarks, we return to the formal development. It is useful to separate the turbulent stress $\btau_\ell$ into a ``systematic'' part $\btaul_\ell$ and a ``random'' part $\btaus_\ell,$ \be \btau_\ell[x;\vl_\ell]=\btaul_\ell[x;\vl_\ell]+\btaus_\ell[x;\vl_\ell]. \lb{45j} \ee (We represent here the space-time point as $x=(\br,t).$ ) Any such division is essentially arbitrary and can be made according to various schemes. The one we have found most convenient is as follows: consider the functional $\Vs_\ell[x;\vl_\ell,\fs_\ell]$ obtained by solving the SE for fixed past histories of $\vl_\ell$ and $\fs_\ell.$ Substituted into the expression for $\bT_\ell$, this can then be averaged over the ensemble of random forces $\fs_\ell$ {\em treating } $\vl_\ell$ {\em as a deterministic quantity} (or, equivalently, as an independent random variable). Denoting this average as $\langle\cdot |\vl_\ell\rangle,$ \footnote{Note that for any functional $F$ of the large-scale velocity alone \be \langle F[\vl_\ell]G[\vl_\ell,\fs_\ell]|\vl_\ell\rangle= F[\vl_\ell]\langle G[\vl_\ell,\fs_\ell]|\vl_\ell\rangle. \lb{45m} \ee Hence, the average $\langle\cdot|\vl_\ell\rangle$ has the appearance of a ``conditional average'' with large-scale velocity fixed, as its notation also suggests. However, this is not true and rather misleading. The {\em solution} $\vl_\ell$ {\em of the SLE} for a prescribed past history of random forces $\fs_\ell$ will develop a functional dependence upon them, so that \be \langle F[\vl_\ell]G[\vl_\ell,\fs_\ell]\rangle\neq F[\vl_\ell]\langle G[\vl_\ell,\fs_\ell] \rangle, \lb{45n} \ee when $\langle\cdot\rangle$ denotes the ordinary average over the ensemble of forces $\fs_\ell.$ It would be possible to consider a true conditional average over $\fs_\ell$ with the {\em solutions} $\vl_\ell$ fixed, but such a conditioning would change the statistics of the forces $\fs_\ell$ from their {\em a priori} Gaussian statistics prescribed by Eq.(\ref{45f}). Therefore, the property Eq.(\ref{45m}) is really a consequence just of the {\em definition} of the linear operation $\langle\cdot |\vl_\ell\rangle$ applied to arbitrary functionals $F[\vl_\ell,\fs_\ell]$ of the variables $\vl_\ell$ and $\fs_\ell.$ Its usefulness will appear in the context of the MSR formalism developed later.} we then set \be \btaul_\ell[x;\vl_\ell]=\langle \bT_\ell[x;\vl_\ell,\cdot]|\vl_\ell\rangle, \lb{45k} \ee which no longer depends upon $\fs_\ell,$ and \be \btaus_\ell[x;\bv_\ell,\fs_\ell]=\bT_\ell[x;\vl_\ell,\fs_\ell]-\btaul[x;\vl_\ell]. \lb{45l} \ee Note also that the ``conditional mean'' of the SS field, \be \vs[x;\vl]\equiv \langle \Vs[x;\vl]|\vl\rangle, \lb{72} \ee will not vanish in general, since the past history of the LS field will set up a SS flow. Therefore, it is natural to separate this mean contribution from the SS field, by making the definition \be \Ws[x;\bv]\equiv \Vs[x;\bv]-\vs[x;\vl]. \lb{73} \ee We can now define a net ``systematic field,'' \be \bv[x;\vl]\equiv \vl(x)+\vs[x;\vl] \lb{74} \ee as a deterministic functional of the LS field, whereas $\Ws[\vl]$ represents the random part. As a final remark, let us note that it is possible in principle to ``defilter'' the solutions of the SLE if \be \fs_\ell=H_\ell*\bF, \lb{45o} \ee where $\bF$ is a random force whose spectral support is disjoint from the support of $\widehat{G}_\ell(\bk).$ In that case, there is a solution of the full Navier-Stokes equation with the force $\bF$ which, when filtered by $G_\ell,$ is a solution of the SLE given above. When the filter functions are graded this is not true for the explicit choice of $\fs_\ell$ suggested in Eq.(\ref{45f}) above. The latter corresponds to filtering with $H_\ell$ the spacetime white-noise force $\bF$ whose covariance is \be \langle \bF(\br t)\bF(\br' t')\rangle= 2\ven_0 \delta^3(\br-\br')\delta(t-t'). \lb{45p} \ee Since $\widehat{G}_\ell$ and $\widehat{H}_\ell$ must have overlapping supports for graded filters and the spatially white-noise force $\bF$ has a flat wavenumber spectrum, a corresponding term $\fl_\ell=G_\ell*\bF$ must then appear in the LE. However, since one does not believe that the results in the limit $\ven_0\rightarrow 0$ will depend too essentially upon the precise form of the force, we expect the same results by using the ``non-defilterable'' force Eq.(\ref{45f}) as a ``defilterable'' one. The possibility of using different types of forcing in the limit $\ven_0\rightarrow 0$ limit, with hopefully equivalent results, will be exploited further later on. \noindent{\em (2.2) Martin-Siggia-Rose Formalism for Statistical Dynamics} We review here the field-theory method of MSR \cite{MSR}. Its fundamental motivation, without technicalities, is the observation that it is difficult to characterize the probability distribution describing a stationary ensemble of turbulent flows. In that case, it may be more profitable to work instead with the ensemble of {\em histories} of the Navier-Stokes dynamics. This can be specified more concretely and, in many respects, provides an effective substitute for the stationary measure. To explain the formal principles, we consider first the example of a randomly forced Navier-Stokes fluid. In this case, the statistical problem to be solved is given by the dynamical equation \be \partial_t\bv +\bP(\grad)(\bv\cdot\grad)\bv=\nu_0\bigtriangleup \bv+\bF, \lb{11} \ee where $\bP(\grad)$ is a solenoidal projection required to maintain the incompressibility condition (replacing a pressure term) and $\bF$ is a random solenoidal force with a Gaussian distribution and covariance \be \langle f_i(\bx,t)f_j(\bx',t')\rangle = P_{ij}(\grad_\bx)F(\bx-\bx')\delta(t-t'). \lb{12} \ee As we shall discuss later, the stochastic nature of the force $\bF$ is not really necessary to our discussion, or any restriction of the method. The statistical problem posed by Eq.(\ref{11}) may also be considerably generalized by employing other classes of random forces besides those Gaussian and white-noise in time. The models with random forcing subsume the problems with random initial data, since those may be represented by an impulsive random force \be \bF(\br t)=\bv_0(\br)\delta(t-t_0) \lb{13} \ee imposed at the initial instant $t_0$ on a quiescent fluid, with the initial datum $\bv_0$ chosen from some selected random ensemble. Although we do not need to assume that $\langle \bF(\br t)\rangle=0,$ it is more convenient to assume this and to add, if desired, an explicit mean force $\bar{\bF}(\br t)$ to the righthand side of Eq.(\ref{11}). In such cases \be \bar{\bv}(\br t)\equiv \langle \bv(\br t)\rangle \lb{13a} \ee will not be zero. If the forcing spectrum $\widehat{F}(\bk)$ is supported only in a small interval near $k_0=1/L,$ then Eqs.(\ref{11}),(\ref{12}) are a relatively realistic model of homogeneous and isotropic turbulence, where the force is just a convenient way of injecting energy at large scales, i.e. stirring the fluid at length-scale $L.$ Rather than following the somewhat mysterious operator formulation originally developed by MSR \cite{MSR} (see also Phythian \cite{Phy1,Phy2}), we shall briefly describe the path-integral formulation of Janssen \cite{Jan} and DeDominicis \cite{DeD}. The simple idea is to write a path-integral representation for the generating functional $Z[\bleta,\hat{\bleta}]$ of correlation and response functions by incorporating the dynamics through a delta functional in its representation by an integral over an exponential. The objects playing the role of ``momentum'' $p$ in the functional integral analogue of $\delta(x)={{1}\over{2\pi}}\int dp\,e^{ipx},$ are the fields $\hat{\bv}(\br t),$ whose joint correlations with $\bv$'s turn out to be the response functions. The expression for the generating functional is just \begin{eqnarray} Z[\bleta,\hat{\bleta}]&=&\int{\cal D}\bF \exp\left[-{{1}\over{2}}\langle \bF, F^{-1}\bF\rangle \right] \cr \,& &\times \int{\cal D}\bv{\cal D}\hat{\bv}\exp\left[ i\int d^d\br\int dt\hat{\bv}(\br t) \left( (\partial_t-\nu_0\bigtriangleup_\br)\bv(\br t) \right. \right.\cr \,& &\,\,\,\,\,\,\,\,\,\,\,\left.+\bP(\grad_\br)(\bv(\br t)\cdot \grad_\br)\bv(\br t)-\bar{\bF}(\br t)\right) \cr \,& &\,\,\,\,\,\,\,\,\,\,\,\left.-i\langle \hat{\bv},\bF\rangle -i\int d^d\br dt \left(\bleta(\br t)\cdot\bv(\br t) +\hat{\bleta}(\br t)\cdot\hat{\bv}(\br t)\right) \right]J[\bv], \lb{14} \end{eqnarray} where $J[\bv]$ is a Jacobian factor which appears in the transformation \be \delta[\bv-{\bf V}(\bF)]=\delta[(\partial_t-\nu_0\bigtriangleup)\bv-\bar{\bF}+\bP(\grad)(\bv\cdot\grad)\bv-\bF]J[\bv], \lb{15} \ee and ${\bf V}(\bF)$ is the exact solution of the dynamical equation Eq.(\ref{11}) with specified force $\bF.$ We note that a discretization of the dynamics in time will always be implicitly assumed in this work to give meaning to such formulae as those above, even where we use continuum notation as a convenient shorthand. Along with this, a space-regularization by means of a spatial grid or a Fourier-Galerkin truncation will be employed. Among various choices of time-discretization, we could use the explicit Euler scheme, \be (\bv_i-\bv_{i-1})/\tau=\bK(\bv_{i-1})+\bF_i, \lb{15a} \ee or an implicit numerical scheme, such as: \be (\bv_i-\bv_{i-1})/\tau=\bK({{\bv_i+\bv_{i-1}}\over{2}})+\bF_i. \lb{15aa} \ee It is easy to check that the form of the Jacobian $J[\bv]$ and the resulting path-integral expression depend upon the discretization adopted. In fact, there is a substantial literature on this point \cite{Gr1,Gr2,Wei,Wis}. If the force in our equation Eq.(\ref{11}) were multiplied by some state-dependent coefficient $\bB(\bv),$ then it is known that the symmetrical time-splitting as in the implicit rule Eq.(\ref{15aa}) must be used in the noise coefficient in order to yield the continuum equation in Stratonovich form \cite{Arn}. For a stochastic Langevin equation the Stratonovich interpretation is usually implicitly assumed since it guarantees validity of the familiar rules of calculus. On the other hand, the explicit or ``non-anticipating'' scheme leads to the Ito form. We calculate here the Jacobian only for these two schemes. Beginning with the symmetrical time-splitting, we note if ${\bf V}_i(\bv_{i-1},\bF_i)$ is the solution of the discretized nonlinear equation Eq.(\ref{15aa}) for $\bv_i$ in terms of $\bv_{i-1}$ and $\bF_i$, it follows from \begin{eqnarray} \,& &\prod_{i}\delta\left((\bv_i-\bv_{i-1})/\tau-\bK({{\bv_i+\bv_{i-1}}\over{2}})-\bF_i\right) \cr \,& &\,\,\,\,\,\,\,\,\,\,\,\,=\prod_{i}\delta\left(\bv_i-{\bf V}_i(\bv_{i-1},\bF_i)\right) \left| {{1}\over{\tau}}-{{1}\over{2}}{{\delta\bK}\over{\delta\bv_i}}\left({{\bv_i+\bv_{i-1}}\over{2}}\right) \right|^{-1} \cr \,& &\,\,\,\,\,\,\,\,\,\,\,\,\approx \prod_{i}\delta\left(\bv_i-{\bf V}_i(\bv_{i-1},\bF_i)\right) {{1}\over{\tau}}\exp\left[{{1}\over{2}}\tau\cdot {{\delta\bK}\over{\delta\bv_i}}\left({{\bv_i+\bv_{i-1}}\over{2}}\right)\right], \lb{16} \end{eqnarray} that the Jacobian factor is \be J[\bv]=\exp\left[-{{1}\over{2}}\int d^d\br dt\,\, {\rm tr}{{\delta \bK(\br t;\bv)}\over{\delta\bv(\br t)}}\right]. \lb{16a} \ee On the other hand, it is clear by the same method that $J[\bv]\equiv 1$ for the Ito discretization. We should remark that for situations where the nonlinear part of the dynamics satisfies a Liouville theorem, the Jacobian in fact is only a field-independent factor and may always be neglected. It was shown long ago by T. D. Lee \cite{Lee} that the Fourier-Galerkin truncation of the Navier-Stokes dynamics falls in this category. We shall rederive this result in Appendix II, where it is shown more generally that the ``effective dynamics'' generated by elimination of the small-scale modes satisfies such a Liouville property. It thus turns out that in all the cases we consider the nonlinear dynamics satisfies the Liouville theorem and we are justified in ignoring the Jacobian factor. Performing finally the Gaussian integration over force $\bF$ by completion of squares, we obtain \be Z[\bleta,\hat{\bleta}]=\int{\cal D}\bv{\cal D}\hat{\bv}\exp\left(iS[\bv,\hat{\bv}]-i\langle\bv,\bleta\rangle -i\langle\hat{\bv},\hat{\bleta}\rangle\right), \lb{17} \ee a path-integral over histories with the ``action'' \begin{eqnarray} S[\bv,\hat{\bv}] &= & \int d^d\br\int dt \hat{\bv}(\br t)\left( (\partial_t-\nu_0\bigtriangleup_\br)\bv(\br t) \right. \cr \,& &\,\,\,\,\,\,\,\,\,\,\,\,\left.+\bP(\grad_\br)(\bv(\br t)\cdot \grad_\br)\bv(\br t)-\bar{\bF}(\br t)\right) \cr \,& &\,\,\,\,\,\,\,\,\,\,\,\,+{{i}\over{2}}\int d^d\br dt\int d^d\br'dt' \hat{v}_i(\br t)F_{ij}(\br t,\br't')\hat{v}_j(\br't'). \lb{18} \end{eqnarray} This expression may be used to calculate arbitrary multitime correlations of the velocity field by functional differentiation with respect to $\bleta(\br t).$ The source $\hat{\bleta}(\br t),$ on the other hand, appears in the action exactly as an external force in the Eq.(\ref{11}), so that the derivatives ${{\delta^{p+1}Z[\bleta,\hat{\bleta}]}\over {\delta\eta_i(\br t)\delta\hat{\eta}_{j_1}(\br_1t_1)\cdots\delta\hat{\eta}_{j_p}(\br_pt_p)}}$ are the higher-order average response functions. It is not difficult to extend the formalism to non-Gaussian or ``colored'' random forces and to forces multiplicative in the random velocity $\bv.$ Indeed, one can consider any ``generalized Langevin dynamics'' \be \partial_t\bv(x)={\bf K}[x;\bv]+\bF'[x;\bv]. \lb{19} \ee Here $x=(\br,t)$ is a space-time point and ${\bf K}[x;\bv]$ is an arbitrary functional of $\bv$ which is``non-anticipating'' or ``causal,'' i.e. that depends only upon the ``past'' values of $\bv(x')$ with $t'0$ \be -i\langle\Omega^-|{\rm T}\left[\hX_n(t)\hP_m(0)\right]|\Omega^+\rangle =\langle {{\partial x_n(t)}\over{\partial x_m(0)}}\rangle. \lb{V56} \ee Now making a small perturbation $\bx(0)\rightarrow \bx(0)+\ben$ in the initial data is the same as making a small perturbation $\bF(t)\rightarrow \ben\cdot\delta(t)$ in the external force. By the definition of functional differentiation it follows that \be {{\partial x_n(t)}\over{\partial x_m(0)}}= {{\delta x_n(t)}\over{\delta f_m(0)}}, \lb{V57} \ee the usual instantaneous response function $\hat{G}_{nm}(t,0)$. The stated identity Eq.(\ref{V53}) immediately follows. \noindent {\em (3.3) Symanzik-Type Theorems for MSR} We can now state and prove the version of Symanzik's theorem which holds in MSR field theory. It requires a very simple modification associated to the non-self-adjoint character of the formal Hamiltonian. More precisely, we have \begin{Th} The effective potential \be V[\bx]=\lim_{T\rightarrow +\infty}{{1}\over{T}}\Gamma[\bx_T], \lb{V58} \ee for a stationary Markov process is the value at the extremum point of the functional \be V[\Psi^+,\Psi^-]=-\langle\Psi^-|\hL|\Psi^+\rangle \lb{59a} \ee varying over all pairs of state vectors $\Psi^+,\Psi^-$ subject to the constraints \be \langle\Psi^-|\Psi^+\rangle=1 \lb{V60} \ee and \be \langle\Psi^-|\hbX|\Psi^+\rangle=\bx. \lb{V61} \ee \end{Th} Whereas the original version of the theorem required just one trial state, now there must be {\em two independent trial states.} Nevertheless, the proof is similar to the original one of Symanzik \cite{Sym}. Observe first that the generating functional $W[\bh]$ introduced earlier may be represented in the operator formulation by \be W[\bh]=\log\langle\Omega^-| {\rm T}\exp\left(\int_0^T dt\,\hL_\bh(t)\right)|\Omega^+\rangle, \lb{V62} \ee where \be \hL_\bh(t)=\hL+\bh(t)\bdot\hbX . \lb{V63} \ee No time-dependence is required for the coordinate operators because the exponential factors automatically introduce the correct Heisenberg picture operators after differentiating and setting $\bh$ to zero. We note then that for a {\em static} field $\bh$ in the limit $T\rightarrow +\infty,$ \begin{eqnarray} \exp(W[\bh_T]) & = & \langle\Omega^-|\exp\left(T\cdot \hL_\bh \right)|\Omega^+\rangle \cr \,& \approx & \langle\Omega^-|\Omega^+[\bh]\rangle\langle\Omega^-[\bh]|\Omega^+\rangle \times \exp(T\cdot\lambda[\bh]), \lb{V64} \end{eqnarray} where $\lambda[\bh]$ is the eigenvalue of the ``perturbed operator'' \be \hL_\bh=\hL+\bh\bdot\hbX \lb{V65} \ee with the {\em largest real part} and $\Omega^+[\bh],\Omega^-[\bh]$ are the associated right and left ``ground state'' eigenvectors: \be \hL_\bh|\Omega^+[\bh]\rangle=\lambda[\bh]|\Omega^+[\bh]\rangle, \lb{V65a} \ee and \be \hL_\bh^*|\Omega^-[\bh]\rangle=\lambda^*[\bh]|\Omega^-[\bh]\rangle. \lb{V66} \ee Furthermore, we can see that \be {{\partial W[\bh_T]}\over{\partial h_n}}= T\cdot x_n[\bh]+o(T), \lb{V66a} \ee with \be x_n[\bh]= \langle\Omega^-[\bh]|\hX_n|\Omega^+[\bh]\rangle. \lb{V67} \ee This can be obtained from the formula \begin{eqnarray} \exp(W[\bh_T]){{\partial W[\bh_T]}\over{\partial h_n}} & = & \langle\Omega^-|{{\partial}\over{\partial h_n}}\exp\left(T\cdot \hL_\bh \right)|\Omega^+\rangle \cr \,& = & \langle\Omega^-|\Omega^+[\bh]\rangle\langle\Omega^-[\bh]|\Omega^+\rangle \langle\Omega^-[\bh]|{{\partial}\over{\partial h_n}}\exp\left(T\cdot \hL_\bh \right)|\Omega^+[\bh]\rangle \cr \,& & \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, +O\left(e^{-T\cdot\Delta\lambda}\right), \lb{V67a} \end{eqnarray} where $\Delta\lambda$ is the spectral gap between the real parts of the ``ground state'' eigenvalue and the next highest eigenvalue. We also use the well-known fact that, for any one-parameter family of operators $\hL(h)$ depending smoothly on a parameter $h,$ \be {{\partial}\over{\partial h}}\exp(\hL(h))=\exp(\hL(h))\varphi(-{\rm Ad}\hL(h)) \left[{{\partial \hL(h)}\over{\partial h}}\right], \lb{V68} \ee where ${\rm Ad}\hL$ denotes the ``adjoint super-operator'' defined by commutator, \be ({\rm Ad}\hL)[\hat{O}]=[\hL,\hat{O}], \lb{V69} \ee and $\varphi(z)$ is the entire function \be \varphi(z)=(e^z-1)/z=1+{{1}\over{2!}}z+{{1}\over{3!}}z^2\cdots. \lb{V70} \ee See \cite{HS}. Since \be \langle\Omega^-[\bh]|[\hL_\bh,\hat{O}]|\Omega^+[\bh]\rangle=0, \lb{V71} \ee for any operator $\hat{O},$ only the first term survives in the expansion of $\varphi$ when substituted into the first term of formula Eq.(\ref{V67a}). This yields Eq.(\ref{V66a}). Now let us consider the variational problem. If we incorporate the constraints by suitable Lagrange multipliers, then the variational equation is just \be \delta\left[ -\langle\Psi^-|\hL|\Psi^+\rangle-\bh\bdot\langle\Psi^-|\hbX|\Psi^+\rangle +\lambda\langle\Psi^-|\Psi^+\rangle\right]=0, \lb{V72} \ee or \be \langle\delta\Psi^-|\hL_\bh-\lambda|\Psi^+\rangle +\langle\Psi^-|\hL_\bh-\lambda|\delta\Psi^+\rangle=0. \lb{V72c} \ee In other words, there are infinitely many stationary points of the functional $V[\Psi^+,\Psi^-]$ subject to the constraints. They consist precisely of pairs $(\Psi_\alpha^+[\bh],\Psi_\alpha^-[\bh])$ of eigenvectors of $\hL_\bh$, \be \hL_\bh|\Psi_\alpha^+[\bh]\rangle=\lambda_\alpha[\bh]|\Psi^+_\alpha[\bh]\rangle, \lb{V72a} \ee and \be \hL_\bh^*|\Psi_\alpha^-[\bh]\rangle=\lambda^*_\alpha[\bh]|\Psi_\alpha^-[\bh]\rangle, \lb{V72b} \ee corresponding to different branches of eigenvalues $\lambda_\alpha[\bh],\alpha=0,1,2,...$ To be precise, we should consider the stationary point corresponding to the branch with largest real part for each $\bh,$ that is, the pair of ``ground state'' eigenvectors $(\Omega^+[\bh],\Omega^-[\bh])$ introduced above. Applying the left eigenvector to the eigen-equation of the right vector and using the constraints gives \be \langle\Omega^-[\bh]|\hL|\Omega^+[\bh]\rangle+ \bh\bdot\bx[\bh]=\lambda[\bh], \lb{V73} \ee and thus \begin{eqnarray} -\langle\Omega^-[\bh]|\hL|\Omega^+[\bh]\rangle & = & \bh\bdot\bx[\bh]-\lambda[\bh] \cr \, & = & {{1}\over{T}}\left[ \langle\bh_T,{{\delta W}\over{\delta\bh}}[\bh_T]\rangle-W[\bh_T]\right]+o(1), \cr \, & = & {{1}\over{T}}\Gamma[\bx_T]+o(1). \lb{V74} \end{eqnarray} The first quantity is independent of $T,$ so that we see taking the limit $T\rightarrow +\infty$ that \be -\langle\Omega^-[\bh]|\hL|\Omega^+[\bh]\rangle = V[\bx], \lb{V75} \ee as was claimed. We have given only a formal proof of the theorem without a careful statement of the conditions, which would certainly involve spectral properties of the ``Liouville operator'' $\hL,$ etc. For deterministic dynamics the existence of a spectral gap in (properly speaking) the Perron-Frobenius operator has been established only for a few special cases, such as the work of Pollicot and Ruelle on Axiom A systems. See \cite{Rue3}. The assumption of a spectral gap is probably stronger than required and a fast polynomial decay, rather than exponential, should suffice. We just make one remark here on the mathematical aspects, which is that trial states $\Psi^+,\Psi^-$ clearly must be taken from different spaces. In fact, $\Psi^+$ varies over the space $L^1,$ of integrable functions of coordinates, while $\Psi^-$ varies over the space $L^\infty$ of bounded functions. Since $L^\infty$ is the Banach space dual to $L^1,$ the Dirac ``bra-ket'' notation, $\langle\cdot|\cdot\rangle$ above, must be interpreted as the canonical dual space action of $L^\infty$ vectors on $L^1$ vectors. Later on we will formulate another more symmetrical version of the result in which both $\Psi^+,\Psi^-$ vary over $L^2.$ The theorem can also be generalized in two important ways. First, there is an analogous variational characterization of the full, time-dependent effective action $\Gamma[\bx].$ Let us state the result formally as \begin{Th} The effective action $\Gamma[\bx]$ for a stationary Markov process is the value at the extremum point of the functional \be \Gamma[\Psi^+,\Psi^-]= \int^{+\infty}_{-\infty}dt\,\,\langle\Psi^-(t),(\partial_t-\hL)\Psi^+(t)\rangle, \lb{V25} \ee when that is independently varied over all pairs of time-dependent state vectors subject to the constraints for each time $t$: \be \langle\Psi^-(t),\Psi^+(t)\rangle=1 \lb{V26} \ee and \be \langle\Psi^-(t),\hbX\Psi^+(t)\rangle=\bx(t), \lb{V27} \ee and also to the boundary conditions \be \lim_{t\rightarrow\mp\infty}|\Psi^\pm(t)\rangle=|\Omega^\pm\rangle. \lb{76} \ee \end{Th} We shall not give here the proof of this theorem because it is essentially the same as the proof of a corresponding result in quantum field theory due to Jackiw and Kerman \cite{JK}. They have shown similarly that the effective action $\Gamma[\phi]$ of field-theory is the stationary point of the functional \be \Gamma[\Psi^+,\Psi^-]= \int^{+\infty}_{-\infty}dt\,\,\langle\Psi^-(t),(ih\!\!\!\!^-\partial_t-\hat{H})\Psi^+(t)\rangle, \lb{V77} \ee varied over pairs $(\Psi^-(t),\Psi^+(t))$ with constraints \be \langle\Psi^-(t),\Psi^+(t)\rangle=1 \lb{V78} \ee and \be \langle\Psi^-(t),\hat{\Phi}(\br)\Psi^+(t)\rangle=\phi(\br,t). \lb{V79} \ee Just as the Symanzik theorem is a constrained version of the familiar quantum variational principle for energy eigenvalues and eigenvectors, the Jackiw-Kerman theorem can be seen as a constrained version of Dirac's \cite{Dir} variational formulation of the Schr\"{o}dinger equation (a quantum analogue of Hamilton's principle). According to that principle, the solutions of the Schr\"{o}dinger equation \be ih\!\!\!\!^-\partial_t\Psi^{\pm}(t)=\hat{H}\Psi^{\pm}(t), \lb{V80} \ee are obtained as the stationary points of $\Gamma[\Psi^+,\Psi^-]$ subject to independent variations of the two wavevectors. In addition to providing a basis for time-dependent Rayleigh-Ritz calculations, the Jackiw-Kerman-type theorem establishes the existence of a Lagrangian functional for the effective action. A second generalization of Symanzik's theorem has been made by Cornwall, Jackiw, and Tomboulis (CJT)\cite{CJT} to the static effective action for higher-order statistics. For example, if $\Gamma[\phi,G]$ is the quantum analogue of $\Gamma[\bv,\bU]$ defined previously, then CJT defined a static version appropriate to discuss those 2nd-order statistics in the time-invariant ground state. Their definition corresponds to substituting into $\Gamma[\phi,G]$ time-invariant versions $\phi(\br)$ and $G(\br,t-t';\br',0),$ with this $G$ then eliminated in terms of the equal-time correlation $G(\br,\br')=G(\br,0;\br',0)$ in a way discussed in detail in \cite{CJT}. After this, the result was divided by length of the total time interval $T$ and the limit $T\rightarrow +\infty$ taken. Then, CJT showed that the resulting quantity $V[\phi,G]$ is also given by the following variational prescription: \be V[\phi,G]=\langle\Psi,\hat{H}\Psi\rangle \lb{V29} \ee at the minimum point subject to constraints \be \langle\Psi,\Psi\rangle=1,\langle\Psi,\hat{\Phi}(\br)\Psi\rangle=\phi(\br), \lb{V30} \ee as before and also \be \langle\Psi,\hat{\Phi}(\br)\hat{\Phi}(\br')\Psi\rangle=\phi(\br)\phi(\br')+G(\br,\br'). \lb{V31} \ee This is the straightforward generalization of Symanzik's original result. It can also be generalized to MSR field theory in a way which should be now obvious. We will refrain here from giving the precise statement and, instead, turn to the practical implementation of the results. \newpage \noindent {\em (3.4) Formulation of Rayleigh-Ritz Variational Methods} We will just outline here a simple variational method of Rayleigh-Ritz type to approximate the effective action and, thereby, the ensemble means. To initiate the method a {\em trial weight} must be chosen, \be w(\bx)\geq 0,\,\,\,\,\,\,\,\,\,\,\int d\bx\,\,w(\bx)=1, \lb{V81} \ee as a plausible {\em a priori} guess for the density $\rho(\bx)$ of stationary measure. This weight will contain a number of parameters which can be optimized by means of the variational principle in Theorem 1. The most straightforward implementation uses expansions in orthogonal polynomials with respect to the trial weight, \be \int P_n(\bx)P_m(\bx) \,w(\bx)\,d\bx=\delta_{n,m}, \lb{V82} \ee which are assumed to form a complete set in the weighted $L^2$-space associated to $w.$ General convergence properties of these polynomials are discussed by Kraichnan in \cite{Kr70} and Section 2 of \cite{Kr85}. In the present application, the pair of state vectors to be varied are expanded to some finite order $N$ as \be \Psi^+= w\cdot \sum_{n=0}^{N-1} c^+_n P_n, \lb{V83} \ee and \be \Psi^-= \sum_{n=0}^{N-1} c^-_n P_n. \lb{V84} \ee Then, incorporating constraints by Lagrange multipliers, the variational equation becomes \be \delta\left[\sum_{n,m=0}^{N-1}c_n^- c_m^+(L_{nm}+\bh\bdot\bX_{nm})-\lambda\sum_{n=0}^{N-1}c_n^-c_n^+\right]=0, \lb{V85} \ee where for all $n,m=0,1,2,....$ \be L_{nm}=\langle P_n,\hL (wP_m)\rangle,\,\,\,\,\,\,\,\,\,\,\bX_{nm}=\langle P_n,\hbX (wP_m)\rangle. \lb{V86} \ee The variation may be instituted in two stages. Varying first with respect to the expansion coefficients $\bc^+,\bc^-,$ the $N$-rank eigenvalue equations are obtained \be \sum_{m=0}^{N-1}(L_{nm}+\bh\bdot\bX_{nm})c_m^+=\lambda c_n^+, \lb{V87} \ee for $n=0,1,...,N-1,$ and \be \sum_{n=0}^{N-1}(L_{nm}+\bh\bdot\bX_{nm})c_n^- =\lambda c_m^- \lb{V88} \ee for $m=0,1,...,N-1.$ Then calculate, for each $\bh,$ the eigenvalue of largest real part and the associated eigenvectors. With $\bx$ fixed, the value of $\bh$ is determined, yielding values $\lambda_N(\bx),\bc^\pm_N(\bx).$ The intermediate approximation of rank $N$ to the effective potential $V_N(\bx)$ is \be V_N(\bx)= -\sum_{n,m=0}^{N-1}\,c_{N,n}^-(\bx)c_{N,m}^+(\bx)L_{nm}. \lb{V89} \ee However, this potential may next be optimized with respect to the remaining parameters in the weight function itself, giving the final approximation. A natural way to implement this procedure in calculating mean turbulent velocity profiles would be to make a Gaussian ansatz for the velocity fluctuations. The mean velocities in the Gaussian trial weight would be taken as the variational parameters whereas the velocity covariance could be fixed by hypothesis. The covariance could be determined, for example, from a model energy spectrum which is $k^{-5/3}$ in the inertial range. In this way, the physics of the K41 theory could be incorporated into the trial state. The orthogonal polynomials with such a Gaussian trial weight are just appropriate multidimensional Hermite polynomials, or polynomials ``Wick-ordered'' with respect to the model covariance. The procedure outlined above would produce a set of variational equations for the mean velocity field at every stage of expansion of state vectors into Hermite polynomials up to the $N$th order. This scheme is essentially a first-order closure method for mean quantities by postulating 2nd-order statistics but with the advantage that it may be systematically improved by increasing $N.$ Alternatively, 2nd-order closures could be implemented variationally as well in which the entire velocity covariance (or some partial 2nd-order statistics such as a local value of the Kolmogorov constant, or the local turbulence intensity, etc.) would be allowed to range over some trial set and then varied to optimize the guess. Clearly, similar methods could also be used in time-dependent problems, such as turbulence growth under a mean constant shear, or freely-decaying turbulence, etc. based upon the variational characterization of the full time-dependent effective action. We do not present any such realistic applications at this point. We shall only examine here a couple of very simple, exactly soluble stochastic models, in order to allow easy comparison of the variational method with reality. First, let us consider the standard Ornstein-Uhlenbeck process (cf. \cite{Ris}) associated to the linear Langevin equation \be \partial_t x= -\gamma (x-x_0)+ \sqrt{D}\eta, \lb{V90} \ee with $\langle\eta(t)\eta(t')\rangle=2\delta(t-t'),$ and Fokker-Planck operator \be \hL= \gamma{{\partial}\over{\partial x}}\left( (x-x_0)\cdot\right)+ D{{\partial^2}\over{\partial x^2}}. \lb{V91} \ee As is well-known, the stationary density is just a Gaussian \be \rho(x)= {{1}\over{\sqrt{2\pi\sigma^2}}}\exp[-(x-x_0)^2/2\sigma^2], \lb{V92} \ee with $\sigma^2=D/\gamma.$ In particular, the mean position is $\bar{x}=x_0.$ However, suppose we ignore this fact and instead use a reasonable guess of the 2nd-order statistics to try to obtain the mean value variationally. A natural choice of trial weight in this case is also Gaussian \be w(x)= {{1}\over{\sqrt{2\pi\sigma^2}}}\exp[-(x-a)^2/2\sigma^2], \lb{V93} \ee with arbitrary center $a.$ Furthermore, let us make the simple {\em ansatz} that \be \Psi^+(x)= w(x)\cdot c_0^+, \lb{V94} \ee and \be \Psi^-(x)= 1+ c_1^-(x-a). \lb{V95} \ee Thus, we use a ``mixed-order'' Hermite expansion in which the series for $\Psi^+$ is taken only to 0th-order but that for $\Psi^-$ is carried to 1st-order. It is then easy to calculate that \be \langle\Psi^-,\hL\Psi^+\rangle= -\gamma(a-x_0)\cdot c_1^-c_0^+, \lb{V96} \ee that \be \langle\Psi^-,\Psi^+\rangle= c_0^+, \lb{V97} \ee and that \be \langle\Psi^-,\hX\Psi^+\rangle= a\cdot c_0^++\sigma^2 \cdot c_1^-c_0^+ . \lb{V98} \ee Imposing the constraints $\langle\Psi^-,\Psi^+\rangle=1$ and $\langle\Psi^-,\hX\Psi^+\rangle=x,$ it is easy to determine that \be c_0^+=1 \,\,\,\,\,\,\,\,\,\, c_1^-={{x-a}\over{\sigma^2}}, \lb{V99} \ee so that \be \langle\Psi^-,\hL\Psi^+\rangle= -{{\gamma^2}\over{D}}(a-x_0)(x-a). \lb{V100} \ee Taking $V(x)= -\langle\Psi^-,\hL\Psi^+\rangle$ and varying with respect to $a,$ the value $a={{x+x_0}\over{2}}$ is obtained. Substituting this stationary point value, the final estimate is obtained: \be V_{{\rm approx}}(x)= {{\gamma^2}\over{4D}}(x-x_0)^2. \lb{V101} \ee Clearly, the minimum is just $x_0,$ so that the method has produced the correct mean value. Of course, this example is very simple and getting the exact mean value at first try is not expected in general. Notice that the same approximation may be obtained from the ``classical'' action for the process, \be \Gamma_{OM}[x]={{1}\over{4D}}\int_{-\infty}^{+\infty}dt\,[\dot{x}+\gamma(x-x_0)]^2, \lb{V102} \ee i.e. the Onsager-Machlup action. For a general diffusion process this is only the exact effective action in the zero-noise limit $D\rightarrow 0,$ and otherwise it will be just a ``leading'' term in a diagrammatic loop-expansion. It is therefore similar to the ``Hartree-Fock approximation'' used in \cite{CJT}. However, since \be {{\Gamma_{OM}[x_T]}\over{T}}={{\gamma^2}\over{4D}}(x-x_0)^2 \lb{V102a} \ee identically at each finite $T,$ the previous result for $V$ is recovered. It may be worthwhile giving just one more simple example, to see whether the method can succeed when the exact statistics are {\em non-Gaussian}. To that end, we take an exactly soluble 1-dimensional gradient dynamics \be \partial_t x= -\varphi'(x)+ \sqrt{D}\eta, \lb{V103} \ee with potential \be \varphi(x)={{\lambda}\over{4!}}(x-x_0)^4. \lb{V104} \ee As is well-known, the Fokker-Planck operator is now \be \hL= {{\lambda}\over{6}}{{\partial}\over{\partial x}}\left((x-x_0)^3\cdot\right) + D{{\partial^2}\over{\partial x^2}}. \lb{V105} \ee and the stationary density is just given by \be \rho(x)\propto \exp\left[-{{\lambda\cdot(x-x_0)^4}\over{4!\cdot D}}\right]. \lb{V106} \ee Although the fluctuations are non-Gaussian, the mean is still given simply as $\bar{x}=x_0.$ Let us use the same Gaussian trial guess as above for $w$ and the same {\em ansatz} for $\Psi^+,\Psi^-$ as in Eqs.(\ref{V94}),(\ref{V95}). Then Eqs.(\ref{V97}) and (\ref{V98}) are unchanged, while \begin{eqnarray} -\langle\Psi^-,\hL\Psi^+\rangle & = & \lambda\left((a-x_0)\cdot{{1}\over{2}}\sigma^2+{{1}\over{6}}(a-x_0)^3\right) c_1^-c_0^+, \cr \,& = &{{\lambda}\over{\sigma^2}} \left((a-x_0)\cdot{{1}\over{2}}\sigma^2+{{1}\over{6}}(a-x_0)^3\right)(x-a), \lb{V107} \end{eqnarray} imposing the constraints. The approximate effective potential $V$ can be obtained by substituting the value $a$ at the stationary point, which is determined by the variational equation $\partial V_a(x)/\partial a=0,$ or \be \left[ {{1}\over{2}}\sigma^2+{{1}\over{2}}(a-x_0)^2\right](x-a)= \left[(a-x_0)\cdot{{1}\over{2}}\sigma^2+{{1}\over{6}}(a-x_0)^3\right]. \lb{V108} \ee It is possible to solve this cubic polynomial equation for $a,$ but it is easier to determine the approximate mean value without doing so explicitly. In fact, using Eq.(\ref{V108}), it is straightforward to see that \be V_{{\rm approx}}(x)={{\lambda}\over{\sigma^2}}\left({{1}\over{2}}\sigma^2+{{1}\over{2}}(a-x_0)^2\right)(x-a)^2 \geq 0. \lb{V109} \ee Thus, the minimum $V_{{\rm approx}}=0$ is achieved if and only if $x=a.$ Substituting back into the cubic Eq.(\ref{V108}), it is found that \be {{1}\over{6}}(a-x_0)^3+{{\sigma^2}\over{2}}(a-x_0)=0 \lb{V110} \ee at the minimum. This equation has only one real-valued solution $a=x_0.$ Hence, we conclude that \be \bar{x}=x_0, \lb{V111} \ee which is again the exact value. Notice that the variance $\sigma^2$ of the Gaussian trial state never needed to be specified in this argument. In fact, the Gaussian {\em ansatz} gives a rather poor representation of the fluctuations for any choice of $\sigma^2,$ but all choices lead here to the same mean value. Let us make a few remarks on the convergence properties of the systematic expansion procedure described previously in the limit as $N\rightarrow +\infty.$ Although we give here no rigorous proofs, it is clear that it should give a convergent scheme if some reasonable properties are satisfied by the dynamics and the trial weight. A natural condition is that the operator $\hL_w$ defined by the following formal similarity transformation \be \hL_w\equiv {{1}\over{\sqrt{w}}}\hL\sqrt{w}, \lb{V112} \ee be a (generally unbounded) operator on $L^2$ such that \be \{\sqrt{w}P_n: n=0,1,2,...\}\subset {\rm dom}(\hL_w)\bigcap{\rm dom}(\hL_w^*), \lb{V113} \ee and that the perturbed operators $\hL_w+\bh\bdot\hbX$ have a complete set of biorthogonal left and right eigenvectors $\Psi_{w,\alpha}^\pm[\bh],\,\,\,\alpha=0,1,2,...$ in $L^2.$ Notice that if we take $\Omega_w^\pm[\bh]$ to be the ``ground state'' vectors ($\alpha=0$) with largest real part, then $\Omega_w^+[\bz]=\rho/\sqrt{w}$ and $\Omega_w^-[\bz]=\sqrt{w}$. \footnote{It follows from this that, for $w=\rho,$ the two ground states at $\bh=\bz$ of the transformed operator coincide: $\Omega^\pm_\rho=\sqrt{\rho}.$ However, this is probably not helpful, since one must expect the left and right ground states to split under the perturbation $\bh\bdot\hbX.$} In particular, the first implies the usual condition for orthogonal expansion methods that \be \int d\bx\,{{\rho^2(\bx)}\over{w(\bx)}}<\infty. \lb{V114} \ee The importance of the similarity transformation is that it removes the asymmetry of the problem, replacing the two different spaces of trial vectors, $L^1$ and $L^\infty,$ by the single Hilbert space $L^2.$ The variational method may then be reformulated equivalently in terms of the functional \be V_w[\Psi^+,\Psi^-]= -\langle\Psi^-,\hL_w\Psi^+\rangle, \lb{V115} \ee varied subject to the same constraints as before but now with both $\Psi^\pm\in L^2.$ It is clear that the Rayleigh-Ritz method outlined above is equivalent to one for the present principle if the trial states are taken as \be \Psi^\pm(\bx)=\sqrt{w}\sum_{n=0}^{N-1} c_n^\pm P_n(\bx). \lb{V116} \ee This reformulation is necessary to justify the procedure, since the trial states in Eqs.(\ref{V83}),(\ref{V84}) do not belong to the proper spaces $L^1,L^\infty.$ It also allows some proofs of convergence of the $N$th-order approximation $V_N$ pointwise to the true effective potential $V,$ under suitable hypotheses. Note that we consider here convergence with a fixed choice of trial weight $w,$ although one would expect convergence to be improved by optimizing the weight at each stage. Let us just close this section by comparing the variational method we have proposed with some previous ones. Kraichnan (see Section 4.3 of \cite{Kr58}, also \cite{Kr79}) and Qian \cite{Qi} have devised variational schemes, which involve satisfying equations of motion in mean-square sense. Instead, we use only the linear ``Liouville operator'' of the dynamics. More essentially, the scheme advocated by Kraichnan was intended as a scheme to approximate the entire statistical distribution, whereas the principle discussed here is constructed to obtain just mean-values or other low-order statistics. The calculations of Qian \cite{Qi} do show that the $5/3$-spectrum and a reasonable value of the Kolmogorov constant can be obtained by a single-time variational calculation. There may be a connection of our ideas with those of Castaing \cite{Cas}. The action-principle we have established here differs from the ``optimum theory'' of Busse \cite{Bus}, in that it characterizes variationally the true ensemble-averages. Instead, Busse derived variational bounds for turbulent transport properties and he has proposed that the extremalizing vector fields achieving these bounds will be ``similar'' to the ensemble-averaged turbulent fields. However, in general, they will be distinct. Detailed comparison of the virtues and failings of these different principles must await future work. \newpage \section{Stochastic LES Models and Their Applications} \noindent {\em (4.1) A Subgrid Random Coupling Model and Comparison With Other Models} A number of stochastic LES models have already been proposed and implemented for turbulence simulation \cite{Lei,Chas,MT,Car}. In principle, they might be regarded as simplified approximations of the {\em exact} SLE derived here. However, there are certain differences. The model of Chasnov \cite{Chas} was derived by applying the sharp Fourier cutoff filter to the momentum equation of the ``generalized Langevin equation'' model \footnote{ Properly speaking, these equations differ in type from what we have called ``generalized Langevin equations'' because they have the property that noise and damping terms are determined from the past statistics of an infinite ensemble of realizations of the dynamics. A better term to describe this type of stochastic equation would be ``self-consistent Langevin equation''.} for the EDQNM closure equations. As we have discussed elsewhere \cite{I-LES}, the use of the sharp Fourier filter is unwise and introduces non-universal features in the ``subgrid stresses'' which are due to large-scale convection. A more basic difficulty with constructing LES models in this way was noted a long time ago by Kraichnan and Herring \cite{HKr} (p.162). The self-consistent Langevin models were introduced to establish realizability for the closures and, while they are expected to approximate well the low-order {\em statistics} of the turbulence dynamics (e.g. mean energy spectra and transfer), they cannot be expected to provide faithful approximations for individual realizations. In particular, Herring and Kraichnan noted that the models ``scramble the dynamics in a way that makes it implausible that the model could reproduce the build-up of complicated correlations among large numbers of wavevector modes which occurs in solutions of Navier-Stokes equations.'' This is a particularly serious concern if the LES model is to be used to study the coherent structures which evolve spontaneously in the turbulent flow, since production of these structures will be inhibited and their lifetimes reduced. Another type of model was introduced by Leith \cite{Lei} and applied in a somewhat modified form to simulation of a turbulent boundary layer by Mason and Thompson \cite{MT}. Those authors modelled the random accelerations given in our Eq.(\ref{61}) by an application of the same energy-balance ideas used in deriving the Smagorinsky model and by Kolmogorov-style dimensional reasoning. More precisely, a heuristic dimensional argument was given that the rate of energy backscatter from turbulent eddies of size $\ell_e$ to the resolved scale at the filter length $\ell_f$ is of the form $C_B\left({{\ell_e}\over{\ell_f}}\right)^5\ven.$ Here $\ven$ is the local value of energy dissipation, which was estimated by assuming local energy balance of subgrid flux (from a Smagorinsky stress model) with dissipation and backscatter: \be \nu \overline{S}^2=\ven+C_B\left({{\ell_e}\over{\ell_f}}\right)^5\ven. \lb{X1} \ee The backscatter was then implemented in the scheme by generating a field of independent random vectors $\bphi=(\phi_x,\phi_y,\phi_z)$ of zero mean at each point of the spacetime grid. This random vector field was used to generate a ``vector potential'' whose curl gave the random acceleration. The random vector $\bphi$ was filtered on the scale $\ell_f$ to obtain the components of the potential and rescaled appropriately to give the estimated contribution to backscatter rate in root mean square sense. Compared with the exact Eq.(\ref{61}), this model acceleration field has one clearly unnatural feature, that it is given as the curl of a vector $\bA=\grad\btimes\bphi$ rather than the divergence of a symmetric stress tensor $\bA=-\grad\bdot\btaus$ plus the associated contribution to pressure. Mason and Thompson argued that only the divergence-free part of the random acceleration needed to be considered. However, the random backscatter contribution to pressure could make an important contribution, e.g. to isotropization of the small scales, and their procedure makes the pressure term unrealistically a deterministic quantity (i.e. a function just of resolved fields). Another difference compared to the exact expression in Eq.(\ref{61}) is that the model acceleration in this scheme is completely uncorrelated in space-time, whereas the exact expression allows long-range correlations in space and memory of past history in time. Note that a heuristic ``multifractal model'' suggests the presence of such correlation effects in energy dissipation \cite{MC}. A model which contains some of these effects can be constructed by applying Kraichnan's ``random coupling'' method \cite{Kr61} in conjunction with the exact stochastic filtering scheme. \footnote{However, these models will have the same deficiency as \cite{Chas} in omitting coherent subgrid structures.} We shall not give full details here but just outline some of the ideas in the simple context of coupled random oscillators. The dynamical variables are now complex numbers $\zl,\zs$ governed by \be \dot{\zl}(t)+i\Bl\zl(t)+ic\zs(t)=\fL(t), \lb{X2} \ee and \be \dot{\zs}(t)+ic^*\zl(t)+i\Bs\zs(t)=\fS(t), \lb{X3} \ee Here $\Bs,\Bl$ are real-valued random frequencies with independent distributions of zero means and variances $\BL=\langle(\Bl)^2\rangle,\BS= \langle(\Bs)^2\rangle,$ also $c$ is a random complex coupling coefficient with mean $\langle c\rangle$ and variance $C=\langle|c|^2\rangle-|\langle c\rangle|^2,$ and $\fL,\fS$ are external random forces with means $\langle\fL\rangle,\langle\fS\rangle$ and covariances \be \langle(\fL(t)-\langle\fL(t)\rangle)(\fL(t')-\langle\fL(t')\rangle)^*\rangle=\FL(t,t')\,\,\,\,\,\,\,\,\,\, \langle(\fS(t)-\langle\fS(t)\rangle)(\fS(t')-\langle\fS(t')\rangle)^*\rangle=\FS(t,t'). \lb{X4} \ee This system can be considered a caricature of a passive scalar, in which the $z$ variables mimic the scalar concentration, $b,c$ mimic the convecting velocity with prescribed statistics, and $f$'s represent scalar sources. A pair of oscillators are considered to represent the division into ``large-scale'' modes $\zl$ and ``small-scale'' modes $\zs.$ Because the dynamics is linear, it is trivial to implement the stochastic filtering with \be \Zs[t;\zl,\fS]=\int_{t_0}^t ds\,\,e^{-i\Bs(t-s)}[-ic^*\zl(s)+\fS(s)]. \lb{X5} \ee Substituted back into the ``LE'', Eq.(\ref{X2}), this yields the exact stochastic equation for $\zl.$ Note that in this example, which is clearly an ``integrable'' case, the exact SLE dynamics becomes deterministic as $\BS, C,$ and $\FS$ all vanish. In that case, the filtering just produces the new term \be \overline{K}[t,\zl]= -|\langle c\rangle|^2\int_{t_0}^t ds\,\,\zl(s)\cdot\zl(t) \lb{X6} \ee in the effective dynamics. To make a closure for the general case, we may apply Kraichnan's random coupling scheme. We introduce as in \cite{Kr61} $N$ independent copies $\zl_{[n]},\zs_{[n]}$ of the above system,$n=0,1,...,N-1,$ and rewrite them coupled together in ``collective variables'' introduced by a discrete ${\bf Z}(N)$-Fourier transform: \be \zl_\alpha={{1}\over{\sqrt{N}}}\sum_{n=0}^{N-1}e^{2\pi i\alpha n/N}\zl_{[n]}\,\,\,\,\,\,\,\,\,\, \zs_\alpha={{1}\over{\sqrt{N}}}\sum_{n=0}^{N-1}e^{2\pi i\alpha n/N}\zs_{[n]}. \lb{X7} \ee Into the coupled equations are then introduced random interaction phases $\theta_{\alpha,\beta,\gamma}, \theta'_{\alpha,\beta,\gamma}$ {\em but only into the terms involving the SS modes (alone or coupled to the LS modes)}. The dynamics of the LS modes are left unaltered. For obvious reasons, we refer to this new system as the ``subgrid random coupling model.'' The result upon transferring back to the ``individual variables'' $\zl_{[n]},\zs_{[n]}$ and taking $N\rightarrow +\infty$ is that the individuals become again uncorrelated in the limit and that each has LS modes now governed by a stochastic effective dynamics of the form \be \partial_t\zl(t)+i\Bl\zl(t)+\int_{t_0}^t ds\,\,\Sigmal(t,s)\zl(s)=\fL(t)+\phil(t), \lb{X8} \ee in which $\Sigmal$ is a kernel representing additional damping and $\phil$ is a new random force with mean $\langle\phil(t)\rangle$ and covariance \be \langle(\phil(t)-\langle \phil(t)\rangle)(\phil(t')-\langle\phil(t')\rangle)^*\rangle=\Phil(t,t'). \lb{X9} \ee These additional terms represent the effects of the eliminated SS modes and can be characterized by the statistics of those modes. Precisely, \be \Sigmal(t,s)= C\cdot \Gs(t,s) \lb{X10} \ee and \be \Phil(t,s)= C\cdot\left(\Zs(t,s)+\langle\zs(t)\rangle\langle\zs(s)\rangle^*\right), \lb{X11} \ee and \be \langle\phil(t)\rangle= -i\langle c\rangle\langle\zs(t)\rangle. \lb{X11a} \ee In addition to the mean $\langle\zs(t)\rangle$ these involve \be \Gs(t,s)=\langle {{\delta\zs(t)}\over{\delta \fS(s)}}\rangle, \lb{X12} \ee the average response function of the SS modes, and \be \Zs(t,s)=\langle (\zs(t)-\langle\zs(t)\rangle)(\zs(s)-\langle\zs(s)\rangle)^*\rangle, \lb{X12a} \ee the SS correlation function. Note that the new terms in the SLE dynamics Eq.(\ref{X8}) are precisely those appearing at the 1-loop level in the line-reverted expansion (corresponding to the first and second diagrams in Section 2.3). All other contributions from the SS modes are eliminated by averaging over the random phases in the coupling interactions in the large-N limit. The SS statistics are then determined by a set of self-consistent equations, which we call ``subgrid DIA equations.'' They involve the SS statistics as well as the corresponding statistical quantities $\langle\zl(t)\rangle, \Gl(t,s), \Zl(t,s)$ of the LS modes. Precisely, the equations are \be \partial_t\langle\zs(t)\rangle+i\langle c\rangle^*\langle\zl(t)\rangle+ \int_{t_0}^t ds\,\,\Sigmas(t,s)\langle\zs(s)\rangle=\langle \fS(t)\rangle, \lb{X13} \ee and \be \partial_t \Gs(t,t')+\int_{t'}^t ds\,\,\Sigmas(t,s)\Gs(s,t')=\delta(t-t'), \lb{X14} \ee and \be \partial_t \Zs(t,t')+\int_{t_0}^t ds\,\,\Sigmas(t,s)\Zs(s,t')= \int_{t_0}^{t'}ds\,\,\left[\FS(t,s)+\Phis(t,s)\right]\Gs(t',s)^*. \lb{X15} \ee In these equations the new quantities are \be \Sigmas(t,s)= \BS\cdot\Gs(t,s)+C\cdot\Gl(t,s), \lb{X16} \ee and \be \Phis(t,s)= \BS\cdot\left[\Zs(t,s)+\langle\zs(t)\rangle\langle\zs(s)\rangle^*\right]+ C\cdot\left[\Zl(t,s)+\langle\zl(t)\rangle\langle\zl(s)\rangle^*\right]. \lb{X17} \ee These subgrid DIA equations may be obtained from a set of exact Schwinger-Dyson integral equations for the SS dynamics ``conditioned'' on a fixed past history of the LS modes. Because of averaging over random phases and the limit $N\rightarrow +\infty,$ they depend only upon the low-order statistics of the LS modes and only the 1-loop terms in the ``self-energies'' $\Sigmas$ and $\Phis$ survive. It may be verified that the combined system of SLE dynamics and subgrid DIA equations satisfy necessary consistency properties, such as conservation in the mean \be \partial_t\langle|\zl(t)|^2\rangle+\partial_t\left[\Zs(t,t)+|\langle\zs(t)\rangle|^2\right]=0 \lb{X18} \ee for $\langle f\rangle,F=0,$ which follow from the existence of a model representation. This approximation to the exact SLE dynamics clearly does retain some memory effects of the SS modes depending upon the past history of the mean statistics of the LS modes. A similar ``subgrid random coupling model'' can be constructed for the nonlinear Navier-Stokes equation, but it is far more complicated and will not be written here in detail. In this case there is a serious defect of the DIA approximation to the exact SLE equations, which is the failure of {\em random Galilei invariance} (RGI). This problem is exactly the same as that in ordinary DIA, discussed at length by Kraichnan in \cite{Kr87,Kr85,Kr77}. The failure of RGI gives rise to qualitatively bad approximations to the damping, noise, etc. calculated from the DIA equations, or, more generally, from any truncation of the perturbation series in the ``line-reverted'' form. These defects can be cured by using Lagrangian reformulations of the line-reverted expansions \cite{Kr77} but only at the price of losing the model representation. Probably the exact ``subgrid DIA'' equations for Navier-Stokes are too complicated anyway to be of practical use in simulations and further approximations must be made along the lines leading to EDQNM-type closures. There is one feature of the ``subgrid DIA'' equations for the coupled random oscillator system which should be mentioned. In the limit as $\BS,C,\FS\rightarrow 0$ these equations degenerate, so that $\Gs(t,t')\rightarrow\theta(t-t'), \Zs(t,t')\rightarrow 0,$ and the mean obeys the simple equation \be \partial_t\langle\zs(t)\rangle+i\langle c\rangle^*\langle\zl(t)\rangle=\langle \fS(t)\rangle. \lb{X19} \ee These results correctly account for the fact that the exact SS dynamics becomes deterministic in that limit. Likewise, when $C\rightarrow 0,$ the covariance $\Phil\rightarrow 0$ and the force $\phil(t)$ generated by the SS modes becomes deterministic. Therefore, no spurious stochasticity is generated in the LS dynamics by the DIA approximation. However, this is possibly an artefact of the linearity of the problem. In the case of nonlinear dynamics, there is no obvious distinction between the subgrid DIA equations for an ``integrable'' case such as Burgers and a ``chaotic'' case such as Kuramoto-Shivashinsky in the limit $\ven_0\rightarrow 0,$ where external noise is removed. However, it may be that a qualitative distinction emerges from their detailed solution. The issue may be formulated as follows: with small randomness $\sim \ven_0$ in the {\em initial data} of the SS modes, the KS equation should generate an $O(1)$ effective noise for the LS modes in a time proportional to a local ``mixing time'' of the dynamics, with a coefficient which is only weakly dependent on $\ven_0$ (e.g. logarithmically). However, the Burgers equation will presumably generate an $O(1)$ effective noise for the LS modes from randomness $\sim\ven_0$ in the initial data of the SS modes in a time which goes rapidly to infinity in the limit as $\ven_0\rightarrow 0.$ It would be interesting to know whether the DIA approximation to the exact SLE equations is sophisticated enough to make this distinction between ``chaotic'' and ``integrable'' dynamics. \noindent {\em (4.2) Atmospheric Predictability and Complex Flows} Let us just conclude this work by a brief mention of some situations of physical interest where ``turbulent eddy noise'' may play a significant role. It was already argued by Mason and Thompson \cite{MT} that stochastic backscatter is necessary to correct errors in the {\em mean} velocity profiles of the turbulent boundary layer produced by a deterministic LES model of Smagorinsky type. Similarly, the importance of turbulent fluctuations generally in determining the means lead us to believe that the least-action principle may have some use in calculating those mean statistics. Mason and Thompson also found in their study significant qualitative differences between individual realizations in LES with and without stochastic backscatter. For example, with backscatter the velocity fields were much ``rougher'' or ``irregular.'' Thus, eddy noise may also have important effects on the character of large-scale structures that evolve in the flow. Another obvious area of interest is the role of eddy noise in the predictability problem for weather forecasting and climate change. An early study of this issue within the TFM closure was made by Leith and Kraichnan \cite{LKr}. It is clear that deterministic LES models of atmospheric flow make a spurious prediction that LS velocities are completely predictable---in principle---given exact information on just the initial LS fields. A stochastic LES model with eddy noise from the unresolved SS fields corrects this obvious defect. It is of interest how such effects of turbulence noise interact with the ``deterministic chaos'' view on the predictability problem \cite{Lor2}. Since even deterministic LES models show generally ``sensitive dependence'' to initial perturbations, adding eddy noise may not produce any qualitative differences. However, a quantitative difference, even by a factor of two or so, would have significant impact on forecasting ability. These issues are currently investigated in the meteorological literature: see \cite{PNNWZ} for a review. Improved physical understanding of the characteristics of turbulence noise would be valuable for this inquiry. \vspace{.3in} \noindent {\bf Appendix I: The Kraichnan-Wyld Perturbation Theory} The perturbation expansion for the Navier-Stokes equation Eq.(\ref{11}) can be developed by writing an exact integral solution, as: \be \bv(\br t)=\bv^{(0)}(\br t)+\int_{-\infty}^tdt'\int d^d\br'\,G^{(0)}(\br t,\br't')\bP(\grad_{\br'}) (\bv(\br't')\cdot\grad_{\br'})\bv(\br't'), \lb{22} \ee where $G^{(0)}$ is the bare response function \be G^{(0)}(\br t,\br't')=(\partial_t-\nu\bigtriangleup)^{-1}(\br t,\br't'), \lb{23} \ee and $\bv^{(0)}$ is the solution of the linearized problem \be \bv^{(0)}(\br t)=\int_{-\infty}^tdt'\int d^d\br'\,G^{(0)}(\br t,\br't')\bF(\br't'). \lb{24} \ee An exact integral solution can also be obtained for the full response tensor, \be \widehat{G}_{ij}(\br t,\br't')\equiv {{\delta v_i(\br t)}\over{\delta f_j(\br't')}}, \lb{25} \ee by taking the functional derivative of both sides of Eq.(\ref{22}) with respect to $\bF(r't')$, yielding \begin{eqnarray} \,& & \widehat{G}_{ij}(\br t,\br't')=G_{ij}^{(0)}(\br t,\br't') \cr \,& &\,\,\,\,\,\,\, +\int_{-\infty}^tdt''\int d^d\br''\,G^{(0)}(\br t,\br''t'')P_{il}(\grad_{\br''}) \cr \,& &\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\times\left[\widehat{G}_{kj}(\br''t'',\br't')(\grad_{\br''})_k v_l(\br''t'') +(\bv(\br't')\cdot\grad_{\br''})\widehat{G}_{lj}(\br''t'',\br't')\right]. \lb{26} \end{eqnarray} Now, iterating the two expressions, Eq.(\ref{22}) and Eq.(\ref{26}), one obtains expansions of the exact $\bv$ and $\widehat{{\bf G}}$ as power series in $\bv^{(0)}$ and ${\bf G}^{(0)}.$ One may substitute the series into the ensemble averages for the covariance \be U_{ij}(\br t,\br't')=\langle v_i(\br t)v_j(\br't')\rangle-\langle v_i(\br t)\rangle \langle v_j(\br't')\rangle, \lb{27} \ee average response tensor \be G_{ij}(\br t,\br't')=\langle \widehat{G}_{ij}(\br t,\br't')\rangle, \lb{28} \ee and mean field $\bar{\bv}(\br t)$ (Eq.(\ref{13a}) above). In the case of Gaussian force $\bF$, the averages over products of the $\bv^{(0)}$'s may be resolved into products of $\bar{\bv}^{(0)}$'s and bare correlation functions \be U^{(0)}_{ij}(\br t,\br't')=\langle v_i^{(0)}(\br t)v_j^{(0)}(\br't')\rangle-\langle v_i^{(0)}(\br t)\rangle \langle v_j^{(0)}(\br't')\rangle, \lb{29} \ee by the use of Wick's theorem. In this way, series expansions for $\bar{\bv},{\bf U},$ and ${\bf G}$ may be developed in terms of the corresponding bare quantities. A graphical interpretation of this iteration procedure was introduced in \cite{Wyl,Kr61}, in which the basic elements are the ``propagators,'' which are the bare response and correlation function, along with the mean field and the bare interaction vertex \begin{eqnarray} \gamma_{3;ijk}(\br t,\br't',\br''t'') & = & -{{1}\over{2}}\left(P_{ij}(\grad_\br)(\grad_\br)_k+P_{ik}(\grad_\br)(\grad_\br)_j\right) \cr \,& &\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\times \delta^d(\br-\br')\delta(t-t')\delta^d(\br-\br'')\delta(t-t''). \lb{30} \end{eqnarray} Diagrammatically one represents the bare field as: \be v^{(0)}_{i}\equiv\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \lb{31}\ee the average bare response function as \be G^{(0)}_{ij}\equiv\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \lb{32} \ee and the bare vertex as \be \gamma_{3;ijk}\equiv\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \lb{33}.\ee In terms of these graphical expressions, one generates the series for $\bv$ by first writing down all tree graphs whose ``limbs'' are given by ${\bf G}^{(0)}$'s and whose ``top branches'' are given by $\bv^{(0)}$'s. To form the series for $\bar{\bv}$ one eliminates all the $\bv^{(0)}$'s in the single trees for $\bv$, either by replacing them singly by $\bar{\bv}^{(0)}$'s or joining them in pairs to yield ${\bf U}^{(0)}$'s, which are represented as \be \bar{v}^{(0)}_{i}\equiv\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\lb{34}\ee and \be U^{(0)}_{ij}\equiv\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \lb{35}\ee The averaging elimination procedure performed on all pairs of trees (joined by at least one ${\bf U}^{(0)}$) gives the series for ${\bf U}.$ Finally, the series for ${\bf G}$ is obtained if one first replaces one ``top branch'' in single trees by a ${\bf G}^{(0)}$ and then performs the averaging elimination of the $\bv^{(0)}$'s. A convenient compendium of some of the topological properties of graphs in the resulting diagrammatic perturbation series is contained in \cite{ZL}. The same perturbation expansion can be obtained as a Feynman diagram expansion for the MSR field-theory. It is convenient to introduce a doublet field $\Phi_\sigma,\,\sigma=\pm$, following \cite{MSR}, as: \be \Phi(\br t)=\left( \begin{array}{c} \Phi_+(\br t) \\ \Phi_-(\br t) \end{array} \right) = \left( \begin{array}{c} \bv(\br t) \\ i\hat{\bv}(\br t) \end{array} \right) \lb{36}. \ee Then the path-integral expression for the generating functional $Z[\eta]=Z[-i\eta_+,\eta_-],$ $\eta=(\eta_+,\eta_-)^\top,$ is just \begin{eqnarray} Z[\eta]&=&\int{\cal D}\Phi \exp\left[-{{1}\over{2}}\Gamma^{(0)}_2(12) \Phi(1)\Phi(2) \right. \cr \,& &\,\,\,\,\,\,\,\,\,\,\,\,\left.+\gamma_1(1)\Phi(1)+{{1}\over{3!}}\gamma_3(123)\Phi(1)\Phi(2)\Phi(3)+\eta(1)\Phi(1)\right], \lb{37} \end{eqnarray} where we have employed a compact notation $(\sigma_1i_1\br_1t_1)\equiv 1$ along with a summation convention on repeated indices. Here, $\gamma_1=(0,\bar{\bF})^\top.$ $\Gamma^{(0)}_2$ is the ``wave operator'' given by the $2\times 2$ matrix in doublet space \be \Gamma^{(0)}_{2;ij}(\br t,\br't')=\left(\begin{array}{ll} 0 & -(\partial_t+\nu_0\bigtriangleup_\br)\delta_{ij}\delta^d(\br-\br')\delta(t-t') \cr (\partial_t-\nu_0\bigtriangleup_\br)\delta_{ij}\delta^d(\br-\br')\delta(t-t') & -F_{ij}(\br t,\br't') \end{array} \right), \lb{38} \ee and $\gamma_3(123)$ is a completely symmetric vertex which is equal to the $\gamma_3$ already defined in Eq.(\ref{30}) when the doublet indices are $(-++),(+-+),$ or $(++-)$ and zero otherwise. Of course, it must be the case---and it is easy to check directly---that $\Gamma^{(0)}_2=(G^{(0)}_2)^{-1},$ where $G^{(0)}_2$ is the ``bare'' matrix correlation function $G^{(0)}_{2;\sigma\sigma'}=\langle \Phi_\sigma^{(0)}\Phi_{\sigma'}^{(0)}\rangle-\langle \Phi_\sigma^{(0)}\rangle\langle \Phi_{\sigma'}^{(0)}\rangle$ which has three non-vanishing entries \be G^{(0)}_2(i\br t,j\br't')=\left(\begin{array}{ll} U^{(0)}_{ij}(\br t,\br't') & G^{(0)}_{ij}(\br t,\br't') \cr \bar{G}^{(0)}_{ij}(\br t,\br't') & 0 \end{array} \right), \lb{39} \ee with $\bar{G}^{(0)}_{ij}(\br t,\br't')=G^{(0)}_{ji}(\br't',\br t)$ a ``(bare) anti-response function'' and with $\bG^{(0)}$ and $\bU^{(0)}$ already defined in Eqs.(\ref{23}) and (\ref{29}). \noindent {\bf Appendix II: Proof of the Liouville Property} It is helpful to begin by giving a proof of Lee's result \cite{Lee}, that the Liouville theorem is valid for the cutoff Euler dynamics associated to the ``dynamical field'' \be \bK^E(x)=-\bP(\grad_\bx)\grad_\bx\cdot\btau^E(x), \lb{71a} \ee with \be \btau^E(x)=\bv(x)\bv(x). \lb{72a} \ee All of the velocity fields are now defined with a highwavenumber cutoff \be \bv(\bx,t)=\sum_{\bk\in {\cal K}}\,\widehat{\bv}(\bk,t)e^{i\bk\bdot\bx}, \lb{71aa} \ee associated to some bounded domain of wavevectors ${\cal K},$ e.g. the sphere ${\cal K}=\{\bk: |\bk|