% 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 $<t.$ It is clear how to construct such 
a functional approximately using a numerical integration routine on a spacetime grid. This solution is 
now formally introduced back into LE by substituting there
\be \btau_\ell=\bT_\ell(\vl_\ell+\Vs_\ell,\vl_\ell+\Vs_\ell). \lb{45h} \ee
As a result $\vs_\ell$ is formally eliminated from that equation, which is closed in terms of the large-scale
velocity $\vl_\ell,$ along with the random force $\fs_\ell$ acting at small scales. We therefore refer to this 
as the ``closed large-eddy equation''. It is clear that the $\vl_\ell$ obtained by filtering the solution 
of the full equation (with an appropriate random force) must also be a solution of the ``closed LE.'' 
The latter inherits important properties from the original Navier-Stokes dynamics such as energy conservation, 
which can be stated as the constancy of 
\be E(t)={{1}\over{2}}\int d\br\,\left(\vl_\ell(\br,t)+\Vs_\ell[\br,t;\vl]\right)
                          \bdot\left(\vl_\ell(\br,t)+\Vs_\ell[\br,t;\vl]\right)
\lb{71} \ee
for $\fs_\ell\rightarrow\bz$ and $\nu_0\rightarrow 0,$ that is, in the inviscid, unforced limit. 

The next step is to take the limit $\ven_0\rightarrow 0$ in the closed equation. We conjecture that in 
this limit the stress $\btau_\ell$ will not become deterministic but that, instead, it will remain a random 
functional of $\vl_\ell$. \footnote{The limit here must be considered in a ``weak'' sense, i.e. for distributions
on histories in path space.} This may naturally be called the {\em stochastic large-eddy hypothesis}. It is 
motivated by the known chaotic behavior of the Navier-Stokes dynamics, which should enhance the tiny random 
accelerations in the small scales to the degree that their effects do not vanish as their amplitude 
is taken to zero. If $\ell$ is chosen to lie within the inertial interval in a turbulent flow, then the 
statistics of $\btau_\ell$ ought to be determined by the ``universal equilibrium state'' postulated by Kolmogorov for 
the small-scales of a turbulent flow. We refer to the closed equation in that limit as the {\em stochastic large-eddy
equation}(SLE) and we write it out here explicitly, 
\be \partial_t \vl_\ell+\grad\cdot(\vl_\ell\vl_\ell+\btau_\ell[\vl])=-\grad \pl_\ell+\nu_0\bigtriangleup\vl_\ell, 
\lb{SLE} \ee
where, as the notation indicates, the stress $\btau_\ell$ is now a functional of $\vl_\ell,$ but a random one depending
upon the realization of the noise. The random force $\fs_\ell$ is simply used as a mechanism to select the correct 
small-scale statistics in the limit as $\ven_0\rightarrow 0.$ Such a weak-noise limit has been exploited in dynamical 
systems theory to select the ``physical'' stationary measure. The characterization of the ``physical measure'' as the 
stochastically stable one is known in many cases to agree with other reasonable characterizations based upon ``smoothness'' 
properties, when the dynamics has sufficiently strong ergodic properties: see the article of Ruelle \cite{Rue1} and 
Section IV.H  of the review of Eckmann and Ruelle \cite{ER}. The previous hypothesis extends this idea to evolving
states in which the small-scale modes are presumed, because of their faster dynamics, to achieve a ``local equilibrium''
relative to the instantaneous value of the large-scale modes. Because there is no sharp separation
of space and time scales in fluid turbulence, one can only hope for a universal law of the fluctuating turbulent
stress $\btau_\ell$ if the filtering length $\ell$ is taken well within the inertial interval. 

>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'<t.$ Also, $\bF'[x;\bv]$ is 
a zero-mean random force which may be ``history-dependent,'' so that its generating functional
\be \int {\cal D}\bF'[\bv]\,e^{-i\langle \hat{\bv},\bF'[\bv]\rangle}
    =\exp\left\{\sum_{p\geq 2}{{(-i)^p}\over{p!}}\int d^{d+1}x_1\cdots\int d^{d+1}x_p
                       D_{i_1\cdots i_p}^{(p)}[x_1,...,x_p;\bv]\hat{v}_{i_1}(x_1)\cdots\hat{v}_{i_p}(x_p)\right\}, \lb{20} \ee
contains cumulants ${\bf D}^{(p)}[\bv]$ which are causal functionals of $\bv.$ This force is state-independent only if the 
cumulants ${\bf D}^{(p)}$ are constants (independent of $\bv$) and Gaussian only if ${\bf D}^{(p)}=0$ for all $p\geq 3.$ 
Then, by the same argument as given above, this generalized Langevin dynamics has an equivalent field-theory formulation 
in terms of an action \footnote{It is known that this is {\em not} the correct continuum form of the action even for a 
general diffusion process with multiplicative noise: see Section 4 of \cite{Gr1}. We emphasize again that it is crucial 
to interpret all such formulae for the discretized dynamics.}
\begin{eqnarray}
S[\bv,\hat{\bv}]& = & \int d^{d+1}x \hat{\bv}(x)\left(\partial_t\bv(x)-{\bf K}[x;\bv]\right)
                                          +{{i}\over{2}}\int d^{d+1}x\,\,{\rm tr}\left({{\delta \bK(x;\bv)}\over{\delta\bv(x)}}\right) \cr
              \, &   & \,\,\,\,\,\,+\sum_{p\geq 2}{{(-i)^{p+1}}\over{p!}}\int d^{d+1}x_1\cdots\int d^{d+1}x_p
                       D_{i_1\cdots i_p}^{(p)}[x_1,...,x_p;\bv]\hat{v}_{i_1}(x_1)\cdots\hat{v}_{i_p}(x_p). \lb{21}
\end{eqnarray}
If this action is used in the formula Eq.(\ref{17}), then a generating functional is obtained for all of the correlation
and response functions of a general Langevin dynamics. In fact, the MSR field-theory formalism is not limited only to 
situations where the production mechanism or energy injection is by a stochastic force. Just as discussed above in the 
context of ``stochastic filtering,'' it is natural in the completely deterministic cases to add to the dynamics a weak 
stochastic force, which is then taken to zero after computation of averages. From the point of view of the previous 
manipulations, the addition of a spacetime white-noise force with covariance $\sim \ven_0$ contributes a term 
$-{{\ven_0}\over{2}}\int d^{d+1}x\,|\hat{\bv}(x)|^2$ to the exponent, which makes the functional integral over 
$\hat{\bv}$ absolutely convergent. 

Our analysis to this point has been formally exact and non-perturbative, even mathematically rigorous if we consider a 
discretization of the dynamics on a spacetime grid. However, for the dynamical problem of Eq.(\ref{11}) there is a 
naive perturbation expansion in the nonlinearity studied by Wyld \cite{Wyl} and Kraichnan \cite{Kr61}. See also
\cite{ZL} and Appendix I for a brief review.  Its basic elements are the ``bare'' response function $\bG^{(0)},$ 
correlation function $\bU^{(0)},$ and average field $\bar{\bv}^{(0)}.$ It is not obvious that this expansion will be 
quantitatively useful in the case of high Reynolds number turbulence, but it does help to build up some intuition about the 
physical meaning of the formal analysis and to give the point of departure of various approximation methods. Clearly, the
``bare expansion'' degenerates in the vanishing noise limit $\ven_0\rightarrow 0.$ In fact, $\bU^{(0)}\rightarrow\bz$ in 
that case. Since any closed loop of response lines vanishes by causality, all that remains are the ``tree graphs'' for 
$\bv$ terminating in $\bar{\bv}^{(0)}$'s, which represent the perturbative solution of the deterministic Navier-Stokes 
equation with a body-force. However, a nontrivial statistical limit may still be obtained by making a reversion of ``bare'' 
propagators and vertices for ``dressed'' or ``renormalized'' ones prior to taking $\ven_0\rightarrow 0.$ Such a reversion  
can be accomplished by non-perturbative functional Legendre transform techniques, as in the work of MSR \cite{MSR},
or by perturbative series reversion as in \cite{Kr77}. 

\noindent {\em (2.3) Statistical Filtering in the MSR Formalism}

We shall now explain the MSR approach to ``statistical filtering'' of the Navier-Stokes dynamics. It was earlier 
described without a detailed account in our paper \cite{I-RNG} which discussed the foundations of the renormalization 
group (RG) method in turbulence. \footnote{It appeared in that context because the definition of the RG operation always uses 
an elimination of highwavenumber degrees of freedom, or filtering, as one of its parts. However, RG supplements this 
with a rescaling of the remaining ``large-scale'' variables by a ``renormalization factor'' $Z(\ell)$ in such a way as
to get a limit as the filtering length $\ell\rightarrow +\infty.$ Essentially this strategy is a generalization of probabilistic 
``central limit theorem'' ideas to a functional context. Furthermore, rather than performing the filtering and 
subsequent rescaling in one step, RG proceeds iteratively, filtering out some fraction of scales (say, ${{1}\over{2}}$ of them) 
at each step and rescaling the remaining ones. This allows the evolution in scale to be visualized as a ``dynamical flow''
in a space of theories, in which complex behavior can build up in the ``long-time'' limit. None of these RG strategies will
be applied here, since we have not found them to be useful for justifying turbulence model equations.} The first step is
to add {\em independent} random forces $\fl$ and $\fs$ to the LE and SE, respectively:
\be \partial_t \vl+\grad\cdot(\vl\,\vl+\btau)=
                               -\grad \pl+\nu_0\bigtriangleup\vl+\fl, \lb{50} \ee
and 
\be \partial_t \vs+\grad\cdot(\vl\vs+\vs\vl+\vs\vs-\btau)=
                       -\grad \ps+\nu_0\bigtriangleup\vs+\fs. \lb{51} \ee
(Here and afterward we omit the subscript $\ell$ for simplicity.) If we used the sharp 
Fourier filter and any homogeneous random force in the full Navier-Stokes equation, 
then the associated LE and SE would automatically contain such independent random forces. However, with 
graded filters, the independence property would not be true. Hence, we introduce the independent random forces 
{\em after} filtering. The amplitudes of both forces will be taken to zero. In fact, the introduction of 
the force in the LE equation is really not necessary, but introduced just as a mathematical device to 
regularize the path integrals. Then the same method used previously may be employed based upon the formula
\begin{eqnarray}
\,& &\prod_{i}\delta\left((\vl_i-\vl_{i-1})/\tau-\overline{K}({{\vl_i+\vl_{i-1}}\over{2}},{{\vs_i+\vs_{i-1}}\over{2}})
-\fl_i\right)\cr
\,& &\,\,\,\,\,\,\,\,\,\,\,\,\times\delta\left((\vs_i-\vs_{i-1})/\tau-K'({{\vl_i+\vl_{i-1}}\over{2}},{{\vs_i+\vs_{i-1}}\over{2}})
-\fs_i\right) \cr
\,& &\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\approx \prod_{i}\delta\left(\vl_i-\overline{\bf V}_i(\vl_{i-1},\fl_i,\vs_{i-1},\fs_i)\right)
\delta\left(\vs_i-{\bf V}_i'(\vl_{i-1},\fl_i,\vs_{i-1},\fs_i)\right) \cr
\,& &\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\times\exp\left[{{1}\over{2}}\tau\cdot 
{{\delta\overline{\bK}}\over{\delta\vl_i}}\left({{\vl_i+\vl_{i-1}}\over{2}},{{\vs_i+\vs_{i-1}}\over{2}}\right)+
{{1}\over{2}}\tau\cdot {{\delta\bK^\prime}\over{\delta\vs_i}}\left({{\vl_i+\vl_{i-1}}\over{2}},{{\vs_i+\vs_{i-1}}\over{2}}
\right)\right], \lb{52}
\end{eqnarray}
where $\overline{{\bf K}}$ and ${\bf K}^\prime$ represent the dynamical terms in the Eqs.(\ref{50}),(\ref{51}). It is shown 
in Appendix II that, just as for the full NS equation, there is a Liouville theorem for the nonlinear terms of the coupled 
LE and SE. Thus, the the Jacobian factor may be neglected. The final result is a path-integral representation of the 
distribution over histories with the following action function
\begin{eqnarray}
S[\vl,\hat{\vl},\vs,\hat{\vs}] & = &  \langle \hat{\vl}, (\partial_t-\nu_0\bigtriangleup)\vl
                                       +\bP(\grad)\grad\cdot(\vl\,\vl+\btau)\rangle 
                                      +{{i}\over{2}}\langle\hat{\vl},\Fl\cdot\hat{\vl}\rangle+ \cr
                             \,&   & \langle \hat{\vs}, (\partial_t-\nu_0\bigtriangleup)\vs
                                       +\bP(\grad)\grad\cdot(\vl\vs+\vs\vl+\vs\vs-\btau)\rangle 
                                      +{{i}\over{2}}\langle\hat{\vs},\Fs\cdot\hat{\vs}\rangle, \lb{53}
\end{eqnarray}
and integration measure ${\cal D}\vl{\cal D}\hat{\vl}{\cal D}\vs{\cal D}\hat{\vs}=
\prod_x\,\,d\vl(x)d\hat{\vl}(x)d\vs(x)d\hat{\vs}(x).$ This measure only exists as a formal expression, unless we consider 
the discretized dynamics on a spacetime grid.

The ``stochastic filtering'' is then carried out in the MSR formalism by simply integrating out the SS fields 
$\vs,\hat{\vs}$ to give an ``effective action'' $S_{{\rm eff}}[\vl,\hat{\vl}]$ for the LS fields. In fact, we shall 
show that performing this integration gives an MSR action of the form which we previously derived starting from 
a ``generalized Langevin equation'':
\begin{eqnarray}
S_{{\rm eff}}[\vl,\hat{\vl}]& = & \langle\hat{\vl},(\partial_t\vl-\Kl_\ef[\vl])\rangle \cr
                             \, &   & \,\,\,\,\,\,+\sum_{p\geq 2}{{(-i)^{p+1}}\over{p!}}
                                         \bD^{(p)}_\ef[\vl](1\cdots p)\hat{\vl}(1)\cdots\hat{\vl}(p). \lb{54a}
\end{eqnarray}
\newpage

\noindent We shall show that the systematic part of the dynamics obeys a Liouville theorem, so that the functional 
trace vanishes which would appear from the evaluation of a Jacobian: 
\be {\rm Tr}\left({{\delta \Kl_\ef[\vl]}\over{\delta\vl}}\right)=0. \lb{55} \ee
Therefore, it will be established that the LS velocity field obeys an ``effective dynamics'' in the generalized 
Langevin form:
\be \partial_t\vl(x)=\Kl_\ef[x;\vl]+\fl_\ef[x;\bv]. \lb{56} \ee
Furthermore, the derivation will show that the terms in Eq.(\ref{56}) are given explicitly as
\be \Kl_\ef[x;\bv]=-\bP(\grad)\grad\bdot\left(\vl\,\vl+\btaul[x;\vl]\right), \lb{57} \ee
in which $\btaul[x;\vl]$ is a symmetric tensor and causal functional of $\vl$, and, likewise,
\be \fl_\ef[x;\bv]= -\bP(\grad)\grad\bdot\btaus[x;\vl], \lb{58} \ee
where $\btaus[x;\vl]$ is a random symmetric tensor field with a distribution dependent only on the past history
of $\vl.$ From these results the exact representation of $\btau[\vl]$ in Eq.(\ref{41}) is obtained. 

We begin by making the following decomposition of the full MSR action for the original problem Eq.(\ref{11})
into three separate terms:
\begin{eqnarray}
S[\vl,\hat{\vl},\vs,\hat{\vs}] & = &  \langle \hat{\vl}, (\partial_t-\nu_0\bigtriangleup)\vl
                                       +\bP(\grad)\grad\cdot(\vl\,\vl)\rangle 
                                       +{{i}\over{2}}\langle\hat{\vl},\Fl\cdot\hat{\vl}\rangle \cr
               \,& &  +\langle \hat{\vs}, (\partial_t-\nu_0\bigtriangleup)\vs
                                       +\bP(\grad)\grad\cdot(\vl\vs+\vs\vl+\vs\vs-\bT(\vl+\vs,\vl+\vs))\rangle 
                                       +{{i}\over{2}}\langle\hat{\vs},\Fs\cdot\hat{\vs}\rangle \cr
               \,& & +\langle \hat{\vl}, \bP(\grad)\grad\cdot\bT(\vl+\vs,\vl+\vs)\rangle. \lb{59} 
\end{eqnarray}
The principle of the decomposition is this: the first part contains all the terms involving {\em only}
LS fields; the second and third parts contain {\em all} of the SS fields including their couplings
to the LS fields; of the latter two, only the third contains $\hat{\vl},$ the LS response field. We could, if 
we wished, include the contribution from the Leonard stress $\bL(\vl,\vl)$ to the third part instead with 
the first part, but it is a little easier to keep all of the contributions to turbulent stress together.
Now, if the integration over SS fields is performed, then the terms in the first part come out of the integral. 
The second line is the MSR action of the SS field $\vs$ ``conditioned'' on a given history $\vl$ of the LS field, 
\be S[\vs,\hat{\vs}|\vl]= \langle \hat{\vs}, (\partial_t-\nu_0\bigtriangleup)\vs
                                       +\bP(\grad)\grad\cdot(\vl\vs+\vs\vl+\vs\vs-\bT(\vl+\vs,\vl+\vs))\rangle 
                                       +{{i}\over{2}}\langle\hat{\vs},\Fs\cdot\hat{\vs}\rangle. \lb{60} \ee
We can also introduce a random spacetime field  
\be \bA(x;\vl,\vs)= -\bP(\grad_\bx)\grad_\bx\cdot\bT(\vl+\vs,\vl+\vs), \lb{61} \ee
which, for a given solution of the conditional SS dynamics, represents the ``acceleration history'' on LS fields
from the interaction with and by SS's. Note that it is a {\em function} of $\bv(\cdot,t)=\vl(\cdot,t)+\vs(\cdot,t)$---
not a functional of the entire past history. Then, the LS ``effective action'' is given exactly as
\begin{eqnarray} 
S_{{\rm eff}}[\vl,\hat{\vl}] & = & \langle \hat{\vl}, (\partial_t-\nu_0\bigtriangleup)\vl
                                          +\bP(\grad)\grad\cdot(\vl\,\vl)\rangle 
                                          +{{i}\over{2}}\langle\hat{\vl},\Fl\cdot\hat{\vl}\rangle \cr
                               \,&   & \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,+\Delta_{{\rm eff}}[\vl,\hat{\vl}], \lb{62}
\end{eqnarray}
in which $\Delta_{{\rm eff}}$ is represented by the path-integral
\be \exp(i\Delta_{\rm eff}[\vl,\hat{\vl}])=\int {\cal D}\vs{\cal D}\hat{\vs}\exp\left(iS[\vs,\hat{\vs}|\vl]
                                                    -i\langle\hat{\vl},\bA(\bv)\rangle\right). \lb{63} \ee
Assuming that $\Delta_{{\rm eff}}[\vl,\hat{\vl}]$ is functionally analytic in $\hat{\vl},$ we may write out 
its Taylor series expansion  
\be \Delta_{{\rm eff}}[\vl,\hat{\vl}]=\sum_{p\geq 1}{{(-i)^{p+1}}\over{p!}}
                                         \bD^{(p)}_\ef[\vl](1\cdots p)\hat{\vl}(1)\cdots\hat{\vl}(p). \lb{64} \ee
In terms of the diagrammatic expansion discussed in the previous subsection the graphical interpretation of
the perturbation terms contributing to $\bD^{(p)}_\ef(1,...,p)$ is that it consists of all {\em connected vacuum graphs} 
of the field theory with vertices generated by the nonlinear terms in $S[\vs,\hat{\vs}|\vl]$ and only SS propagators, 
and with {\em exactly} the insertions $\hat{\vl}(1),...,\hat{\vl}(p)$ and arbitrary insertions of $\vl.$ If we then take 
\be \Kl_\ef[x;\vl]\equiv \nu_0\bigtriangleup_\bx\vl(x)-\bP(\grad_\bx)\grad_\bx\cdot[\vl(x)\vl(x)]+\bD^{(1)}_\ef[x;\vl], 
\lb{65} \ee
we obtain a total ``effective action'' of the LS fields exactly in the form of Eq.(\ref{21}). We must only show that 
the $\Kl_\ef$ has the Liouville property Eq.(\ref{55}) when $\nu_0=0.$ 

We show that $\Kl_\ef$ satisfies the Liouville theorem (for $\nu_0=0$) to all orders in the 
diagrammatic perturbation series. We should note that the existence of the 
Liouville theorem is a necessary result for consistency of two facts: (i) that the effective LS dynamics is given
by a generalized Langevin equation (GLE) and (ii) that no $\hat{\vl}$-independent term appears in the expansion of 
$\Delta_{{\rm eff}}$ (which would correspond to a $p=0$ term in its series expansion in $\hat{\bv}$). The perturbative 
demonstration of (ii) is based on causality. Indeed, consider any graph with no $\hat{\vl}$-insertions and only 
$\vl$-insertions. Starting at any vertex and following the unique path along response lines, one must eventually arrive 
back at a vertex previously visited, because there are only a finite number of vertices and there is no ``outlet'' to an 
external $\hat{\vl}.$ However, any graph with a closed loop of retarded response lines vanishes, QED. This argument and 
the previously established representation of LS dynamics by a GLE implies the Liouville theorem. It can also be verified 
directly, as we show in Appendix II. Thus we have now verified that the elimination procedure in the MSR path-integral 
yields a ``generalized Langevin equation'' for the LS modes.

To make connection with our previous derivation, let us consider how the elimination step is accomplished 
formally. Integrating over $\hat{\vs}$ recovers the representation  
\begin{eqnarray}
\exp(i\Delta_{\rm eff}[\vl,\hat{\vl}]) & = & \int {\cal D}\vs \int {\cal D}P[\fs] 
\delta\left[\partial_t\vs+\bP(\grad)\grad\cdot[\vl\vs+\vs\vl+\vs\vs-\bT(\bv,\bv)]-\nu_0\bigtriangleup\vs-\fs\right] \cr
                                   \,&   &\,\,\,\,\,\times\exp\left(-i\langle\hat{\vl},\bA(\bv)\rangle\right). \lb{66} 
\end{eqnarray}
Given the LS history $\vl,$ the SE equation may be integrated forward from the initial time to give the SS field
as a solution 
\be \vs(x)=\Vs[x;\vl,\fs], \lb{62a} \ee
where the dependence on both $\vl$ and $\bF'$ is causal, exactly as before. We may now perform the integration 
over $\vs,$ noting again that the Jacobian factor is constant and can be absorbed into the integration measure, giving:
\be \exp(i\Delta_{\rm eff}[\vl,\hat{\vl}])=\int {\cal D}P[\fs] 
                       \exp\left(-i\langle\hat{\vl},\bA(\vl+\Vs[\vl,\fs])\rangle\right). \lb{67} \ee
Notice that the functionals $\bD^{(p)}_{\rm eff}$ introduced above by the Taylor series expansion of $\Delta_{{\rm eff}}$ 
in $\hat{\bv}$ are exactly the $p$th cumulants of the random variable $\bA[x;\vl,\fs]$ distributed with 
respect to $P[\fs]$ {\em with} $\vl$ {\em considered fixed.} If we represent this average with $\langle\cdot |\vl\rangle$
as before, then we may separate the turbulent stress $\bT$ into the ``systematic part''
\be \btaul[x;\bv]=\langle \bT[x;\bv,\cdot]|\bv\rangle, \lb{66a} \ee
and ``random part''
\be \btaus[x;\bv]=\bT[x;\bv,\cdot]-\btaul[x;\bv], \lb{67a} \ee
with zero ``conditional mean.'' This gives the LS effective dynamics in exactly the form earlier claimed
as in Eqs.(\ref{56})-(\ref{58}). For example, 
\be \bD^{(1)}_{\rm eff}[x;\bv]=-\bP(\grad_\bx)\grad_\bx\cdot\btaul[x;\bv] \lb{68} \ee
while
\be \bF^\prime[x;\bv] =-\bP(\grad_\bx)\grad_\bx\cdot\btaus[x;\bv]. \lb{69} \ee
Our choice of division into ``random'' and ``systematic'' terms was made precisely because of 
its natural appearance in the MSR method. The fact that the functionals $\bD^{(p)}_{\rm eff}$ are the cumulants of an actual
random field $\bA$ has an importance which cannot be over-emphasized. For example, it implies that an infinite
number of {\em realizability inequalities} for the $\bD^{(p)}_{\rm eff}$'s are satisfied, of which the lowest-order is 
\be \bD^{(2)}_{\rm eff}\geq 0. \lb{70} \ee
The inequality here is to be interpreted as positive-definiteness in the operator-sense. 

Our analysis has been so far exact (and even, for a numerical discretization, rigorous.) However, to get some intuition 
about the physical terms which arise from our ``statistical filtering'' it is useful to consider the perturbation series 
evaluation of $\Delta_{{\rm eff}}.$ It is then convenient to set $\bv=\vl+\vs[\vl],$ as in Eq.(\ref{74}), and make the 
following change of variables in the path-integral formula, $\hat{\vl}\rightarrow\hat{\bv},\hat{\vs}\rightarrow \hat{\ws},
\vs\rightarrow \vs[\vl]+\ws,.$ This leads to the new MSR action
\begin{eqnarray}
S[\bv,\hat{\bv}] & = & \langle \hat{\bv}, (\partial_t-\nu_0\bigtriangleup)\bv
                                       +\bP(\grad)\grad\cdot(\bv\bv)\rangle 
                                       +{{i}\over{2}}\langle\hat{\bv},\Fl\cdot\hat{\bv}\rangle \cr
               \,& & +\langle \hat{\ws},(\partial_t-\nu_0\bigtriangleup)\ws
                                       +\bP(\grad)\grad\cdot\left(\vs\ws+\ws\vs+\ws\ws-\bT(\bv+\ws,\bv+\ws)\right)
                                       +{{i}\over{2}}\langle\hat{\ws},\Fs\cdot\hat{\ws}\rangle \cr
               \,& & +\langle \hat{\bv}, \bP(\grad)\grad\cdot\bT(\bv+\ws,\bv+\ws)\rangle, \lb{77} 
\end{eqnarray}
The Feynman rules are changed for the calculation of the functions $\bD^{(p)}_{\rm eff}$ by the simple modification that 
insertions of $\bv$ occur rather than of $\vl.$ The generalized Langevin equation is then also given in an altered form, as  
\be \partial_t\bv+\bP(\grad)\grad\cdot(\bv\bv+\bT[\bv])=0. \lb{78} \ee
As before, the separation of $\bT$ into systematic and random parts may be made. 

To second order in the formal coupling parameter, four terms arise in $\Delta_{{\rm eff}},$ which were already noted in 
\cite{I-RNG} (and independently in \cite{Car}). 
They are:

\vspace{1in}

\noindent which gives the linear damping on $\bv$ due to the eliminated high-wavenumber modes,

\vspace{1in}

\noindent and

\vspace{1.25in}

\noindent which represent the generated noise terms, and

\vspace{1.25in}

\noindent which is the lowest order cubic nonlinearity $\sim \bv^3$ in the effective LS dynamics. 
Except for the ``multiplicative'' noise term (the third), these contributions were already noted in 1977 by Rose 
in the problem of passive scalar diffusion. In that case, the fluctuating equation was the conservation of 
scalar concentration $\phi$, or, $\partial_t\phi+\grad\cdot\bj=0,$ and all of the terms were observed to give 
contributions (systematic or random) to the concentration flux $\bj.$ Rose referred to the first as ``eddy diffusivity'' 
(here, ``eddy viscosity'') and the second as ``eddy noise.'' He referred to the last term as ``eddy-mediated diffusion'' 
(here, diffusion of momentum rather than of scalar concentration), which is a nonlocal transport effect due to the ``disappearance''
of the conserved substance into the subgrid modes and its subsequent ``re-emergence'' into the resolved modes in another 
location at a later time. For inhomogeneous flow there is an additional fifth term generated at first-order in the coupling,
which is represented by the graph

\vspace{1.25in}

\noindent It is independent of the systematic velocity field and represents a ``turbulent body-stress,'' or pressure tensor
due to the small-scale turbulent fluctuations. 

It should be observed that higher-order contributions besides these five terms cannot be considered negligible. 
We see no good theoretical argument that Yakhot's conjecture about large-scale equivalence of KS and noisy Burgers
should be {\em exactly} correct. RG methods do not support ``irrelevance'' of additional contributions to the dynamics 
at asymptotically large scales. As we have discussed in detail in \cite{I-RNG}, additional terms such as the new cubic 
nonlinearity (4th diagram above) are {\em marginal by power-counting} for any RG that respects Galilei covariance. Hence, 
that term ought to remain at large length-scales in the KS dynamics. \footnote{Note we have shown also in \cite{I-RNG} that 
those terms may be legitimately neglected using an ``improved'' RG method, but that revised method only works in the 
weak-coupling regime for small $\en.$ It is not helpful in the case of 1D noisy Burgers, which involves a strong-coupling
fixed point.}  However, we agree with the essence of Yakhot's conjecture that the large-scale dynamics should be 
stochastic and {\em similar} to noisy Burgers.  

Another defect of perturbation theory has already been mentioned (end of Section (2.2)), namely, that the noise terms all 
{\em vanish} when $\ven_0\rightarrow 0$ if bare propagators are used. This can be seen to occur in the work of Yakhot, for 
his noise term, Eq.(12) in \cite{Yak}, goes to zero in the limit $D(k)\rightarrow 0$ as his initial data become 
deterministic (see his Eq.(8)). Line-reverted or high-order vertex-reverted expansions seem to be potentially useful here.
In fact, V. S. L'vov et al. \cite{LLPP} have attempted to use such a ``self-consistent expansion'' to establish 
scaling predictions of Yakhot's conjecture. However, their ``proof'' still uses a weak-coupling expansion without 
justification where the interaction is large. Furthermore, it does not seem that any of these perturbative methods are 
sophisticated enough to distinguish between ``chaotic'' and ``integrable'' dynamics. The perturbative argument of
\cite{LLPP} uses no more than the gross form of the nonlinearity and could be used also to ``prove'' equivalence
of deterministic and noisy Burgers! Clearly, KS and Burgers as deterministic equations should behave very differently 
under ``stochastic filtering,'' since only the former should remain noisy, with some universal law of fluctuations, in the 
limit $\ven_0\rightarrow 0,$ whereas the noise in the latter will depend upon the imposed random forces and presumably 
vanish entirely as they are taken to zero. 

A concern with any approximation scheme attempted is the realizability constraints on the noise cumulants $\bD^{(p)}.$ 
Notice from our {\em exact} expressions that all realizability inequalities for the $\bD^{(p)}$ will obtain in the limit
$\ven_0\rightarrow 0,$ although they may conceivably degenerate (i.e. all $\bD^{(p)}\rightarrow \bz.$) One way to
ensure such constraints are satisfied is, as Kraichnan has emphasized \cite{Kr61}, by constructing model realizations
for the approximations. In the last Section 3 we shall describe how a direct-interaction approximation may be 
considered in which only the terms for the 1-loop graphs above are legitimately retained in the SLE dynamics and which 
represents the {\em exact} SLE dynamics of a certain ``random-coupling model.'' In that case, the propagator lines are the 
full propagators determined from the solution of a self-consistent ``DIA-equation'' for the SS mode statistics, into which 
are coupled the mean-values of the LS modes. 

\section{Principle of Least Effective-Action}

\noindent {\em (3.1) Derivation and Comparison with the Onsager Principle}

The formulation by MSR \cite{MSR} of classical statistical dynamics as a formal quantum field theory 
allows one to carry over many techniques from that field. In particular, there is a minimum principle
for vacuum expectations of field operators in Euclidean quantum field theory, which is due to K. Symanzik (see
the Appendix of \cite{Sym}.) This principle involves the ``effective action'' of the field-theory, a concept which has 
its roots in the early work of Heisenberg and Euler \cite{HE} and J. Schwinger \cite{Sch} in QED. However, it turns out that 
the equivalent principle for nonequilibrium statistical dynamics was introduced even earlier in the famous 1931 work of 
L. Onsager on the reciprocal relations \cite{O} (see Section 5 of both his papers.) This idea was further 
developed in a later paper of Onsager with his student S. Machlup \cite{OM}. In the situation they
considered the random force corresponded to molecular noise obeying the fluctuation-dissipation relation
and the action function had the physical meaning of a ``dissipation function.'' Therefore, the action principle 
corresponded to a ``principle of least dissipation.'' Although originally developed only for linear transport
processes, the principle can be extended to the nonlinear regime \cite{Gr3,I-LD}. Furthermore, as pointed 
out by R. Graham (\cite{Gr3}, Section 7), the principle is then seen to be precisely the analogue of the 
``least effective action principle'' of Euclidean field-theory.  

We shall give here a self-contained discussion of the least-action principle, following closely 
the accounts in \cite{Sym,Gr3}. It will be seen that the principle has very general basis and that,
in fact, its origin is the same as that of the familiar equilibrium variational principles of 
maximum entropy, minimum free-energy, etc. The main requirement for its validity in a global sense
is {\em finite exponential moments} of the statistical distribution. That is, it is required that 
\be \int {\cal D}P(\bv)\,e^{\langle \bF,\bv\rangle}<\infty, \lb{V1} \ee
where $\bF$ is a real-vector valued test function. We have in mind here that the distribution is 
over spacetime histories of a turbulent velocity field, but $P$ could really be any probability 
measure whatsoever. We will see that a considerably weaker condition than finite exponential moments 
is sufficient for the principle to be valid in just a local sense. However, if Eq.(\ref{V1}) holds, 
then we may define
\be W[\bF]\equiv \log\left[\int {\cal D}P(\bv)\,e^{\langle \bF,\bv\rangle}\right], \lb{V2} \ee
which is a cumulant-generating functional of the distribution $P.$ It is a consequence of the 
positivity of the distribution and the H\"{o}lder inequality that 
\be \int {\cal D}P(\bv)\,e^{\langle \lambda\bF_1+(1-\lambda)\bF_2,\bv\rangle}\leq
     \left(\int {\cal D}P(\bv)\,e^{\langle \bF_1,\bv\rangle}\right)^\lambda
     \left(\int {\cal D}P(\bv)\,e^{\langle \bF_2,\bv\rangle}\right)^{1-\lambda}, \lb{V3} \ee
for $0<\lambda<1,$ or
\be W[\lambda\bF_1+(1-\lambda)\bF_2]\leq \lambda W[\bF_1]+(1-\lambda)W[\bF_2]. \lb{V4} \ee
In other words, $W[\bF]$ is a globally convex functional of its argument. Observe that this 
is a result just of a simple realizability inequality for the distribution. 

However, it is easy to see that  
\be {{\delta W[\bF]}\over{\delta f_i(\br t)}}= v_i[\br t;\bF], \lb{V4a} \ee
where the latter is the expectation of the velocity-field in the distribution weighted
by the exponential factor. Therefore, the corresponding conjugate convex functional is
\be \Gamma[\bv]=\sup_{\bF}(\langle \bF,\bv\rangle-W[\bF]). \lb{V5} \ee
This is the definition of the {\em effective action}. Observe that it satisfies
\be {{\delta \Gamma[\bv]}\over{\delta v_i(\br t)}}= f_i[\br,t;\bv], \lb{V6} \ee
where $\bF[\bv]$ is the inverse functional to $\bv[\bF].$ Since $\Gamma[\bv]$ is also
globally convex under the assumption Eq.(\ref{V1}), it follows that it has an absolute
minimum (possibly nonunique if $\Gamma$ is not strictly convex). However, since $\vl=\langle \bv\rangle,$
the ordinary expectation of the velocity field \footnote{Observe that our notation here differs from 
that of the other chapters, where $\vl$ represented a large-scale velocity obtained by filtering.},
is just equal to $\vl=\bv[\bz]$ corresponding to $\bF=\bz,$ it follows that
\be {{\delta \Gamma[\vl]}\over{\delta v_i(\br t)}}= \bz. \lb{V7} \ee
That is, $\vl$ is the stationary point of $\Gamma[\bv],$ which is the point at which $\Gamma$ attains
its absolute minimum. That the mean velocity field is characterized as the point at which 
$\Gamma$ achieves its minimum is just the precise statement of {\em the principle of least effective action}.
 
Even without the assumption Eq.(\ref{V1}) of finiteness of all exponential moments a {\em local} version
of the principle may be formulated, assuming that $W[\bF]$ exists in some small neighborhood of $\bF=\bz$
and {\em finiteness of second moments}.In that case,
\be {{\delta^2 W[\bF]}\over{\delta f_i(\br t)\delta f_i(\br' t')}}= U_{ij}[\br t,\br t;\bF] \lb{V8} \ee
where $\bU[\bF]$ is the 2nd-order velocity cumulant in the exponentially weighted distribution
\be U_{ij}[\br t,\br' t';\bF]=\langle v_i(\br t) v_j(\br' t')\rangle_\bF-\langle v_i(\br t)\rangle_\bF
                              \langle v_j(\br' t')\rangle_\bF. \lb{V9} \ee
Again, it is a simple realizability condition that $\bU\geq 0,$ as an operator defined by the kernel 
in Eq.(\ref{V9}). The conjugate $\Gamma$ may still be defined locally by values at the stationary point
of the quantity in parentheses in Eq.(\ref{V5}). Then Eq.(\ref{V6}) still holds and further 
\be {{\delta^2 \Gamma[\bv]}\over{\delta v_i(\br t)\delta v_i(\br' t')}}= (\bU^{-1})_{ij}[\br t,\br t;\bF[\bv]]. \lb{V10} \ee
By positivity of the latter kernel, it follows that the stationary point $\bv$ of $\Gamma[\bv]$ is at least
a local minimum. 

The comparison with the Onsager principle may be found in \cite{Gr3,I-LD}. As noted, that principle applies
to thermal or molecular noise, governed by a fluctuation-dissipation relation. Since this noise is 
very weak, the corresponding path-integral formula for $W[\bF]$ may be evaluated by steepest descent.
The resulting action is just the ``classical'' one, given by the Onsager-Machlup langrangian functional \cite{OM}.
In other words, there is no ``renormalization'' of the action in this case. On the contrary, in turbulence there
should be, as we have discussed, a strong renormalization of the noise and the effective action will differ
from the ``bare'' MSR action appearing in the path-integral. This makes it far more challenging to calculate
the effective action for turbulence noise than for molecular noise. However, turbulence is for the same
reason a far more promising area of application of the principle. Molecular noise seems too weak in 
general to have important influence in nonequilibrium pattern formation and the ``least dissipation principle''
has few important applications in that area. On the other hand, turbulent fluctuations are large and must play a leading 
role in determining mean turbulent velocity profiles in such inhomogeneous situations as channel flow, flow
past an obstacle, etc. Therefore, the least-action principle should be relevant in those cases. 

This principle may be extended to statistics higher-order than mean quantities, if some stronger 
moment conditions are satisfied. This was noted by Cornwall et al. \cite{CJT} in the context of relativistic 
field-theory. For example, if the following moment condition is satisfied 
\be \int {\cal D}P(\bv)\,e^{\langle \bF,\bv\rangle+{{1}\over{2}}\langle \BF,\bv\otimes\bv\rangle}<\infty, \lb{V11} \ee
with
\be \langle \BF,\bv\otimes\bv\rangle=\int d\br dt\int d\br' dt'\,\,F_{ij}(\br t,\br' t')v_i(\br t)v_j(\br' t'),\lb{V12} \ee
then a generalized generating function may be defined as
\be W[\bF,\BF]\equiv \log\left[\int {\cal D}P(\bv)\,e^{\langle \bF,\bv\rangle+{{1}\over{2}}\langle \BF,\bv\otimes\bv\rangle}
\right]. \lb{V13} \ee
The H\"{o}lder inequality still implies for $0<\lambda<1$ that
\be W[\lambda\bF_1+(1-\lambda)\bF_2,\lambda\BF_1+(1-\lambda)\BF_2]\leq \lambda W[\bF_1,\BF_1]
                                   +(1-\lambda)W[\bF_2,\BF_2]. \lb{V14} \ee
In other words, $W[\bF,\BF]$ is a globally convex functional of both arguments. Furthermore,
\be {{\delta W[\bF,\BF]}\over{\delta f_i(\br t)}}= v_i[\br t;\bF,\BF] \lb{V15} \ee
and 
\be {{\delta W[\bF,\BF]}\over{\delta F_{ij}(\br t,\br' t')}}= {{1}\over{2}}(v_i[\br t;\bF,\BF]v_j[\br' t';\bF,\BF]
                                                               +U_{ij}[\br t,\br' t';\bF,\BF]). \lb{V16} \ee 
Hence, a conjugate convex functional $\Gamma[\bv,\bU]$ may be introduced as 
\be \Gamma[\bv,\bU]=\sup_{\bF,\BF}\left(\langle \bF,\bv\rangle+{{1}\over{2}}\langle\BF,\bv\otimes\bv+\bU\rangle
                                                 -W[\bF,\BF]\right). \lb{V17} \ee
It satisfies
\be {{\delta \Gamma[\bv,\bU]}\over{\delta v_i(\br t)}}= f_i[\br,t;\bv,\bU]+\int d\br' dt'\,\,F_{ij}(\br t,\br' t')
                                                        \bv_j(\br' t'), \lb{V18} \ee
and
\be {{\delta \Gamma[\bv,\bU]}\over{\delta U_{ij}(\br t)}}= {{1}\over{2}}F_{ij}[\br t,\br' t';\bv,\bU]. \lb{V19} \ee
Thus, the average $\vl$ and cumulant $\overline{\bU}$ for the $\bF=\BF=\bz$ ensemble are just the minimum points
of the generalized action:
\be {{\delta \Gamma[\vl,\overline{\bU}]}\over{\delta v_i(\br t)}}
              ={{\delta \Gamma[\vl,\overline{\bU}]}\over{\delta U_{ij}(\br t)}}=0. \lb{V20} \ee
This principle characterizes by a minimum principle also second order statistics such as the 
Reynolds stress $\tau_{ij}(\br t)=\overline{U}_{ij}(\br t,\br t)$ and the turbulent intensity 
$K(\br t)={{1}\over{2}}\overline{U}_{ii}(\br t,\br t).$ It is easy to see how to formulate also restricted versions
of the previous construction just for these quantities, or, in the other direction, to extend the principle
to statistics of 3rd and higher order. 

We therefore see now that---in principle---all turbulence statistics will be characterized as the minimizers
of an appropriate action functional. Needless to say, this is practically useless without a means to calculate the
action. Fortunately, there exist in quantum field theory methods of doing so. One of the methods is the 
semiclassical loop expansion. However, this depends upon the existence of a small parameter (Planck's 
constant $h\!\!\!\!^-$ in the quantum mechanics context) and there is no such small parameter in turbulence. However,
there is another method in field theory which is nonperturbative that is based upon a theorem of Symanzik 
\cite{Sym}. The latter provides a variational characterization of the {\em static effective action}, more
commonly referred to as the {\em effective potential}. This is obtained from the full action by defining, for any 
time-independent c-number field $\phi(\br)$, the time-extended field $\phi_T(\br,t)$ by
\be \phi_T(\br,t)\equiv \left\{ \begin{array}{ll}
                                \phi(\br)& \mbox{if $|t|\leq T/2$} \cr
                                0 & \mbox{ otherwise}. 
                                \end{array} 
                        \right. \lb{V21} \ee 
Then the ``effective potential'' $V[\phi]$ is defined as the infinite-time limit
\be V[\phi]=\lim_{T\rightarrow +\infty}{{\Gamma[\phi_T]}\over{T}}. \lb{V22} \ee
The effective potential is appropriate to determine expected field configurations in the time-invariant ground
state of the theory. What Symanzik proved \cite{Sym} is that this quantity may be characterized as the 
expected value of the quantum Hamiltonian operator $\hat{H}$ over all Hilbert space vectors $\Psi$ 
with fixed field expectation $\phi(\br)$. That is,
\be V[\phi]=\inf\{ \langle\Psi,\hat{H}\Psi\rangle: \langle\Psi,\Psi\rangle=1,
                                     \langle\Psi,\hat{\Phi}(\br)\Psi\rangle=\phi(\br)\}. \lb{V23} \ee
Symanzik's theorem corresponds to a constrained version of the well-known variational principle
for the ground state energy in quantum mechanics, according to which
\be E_{{\rm ground}\,\,\,\,{\rm state}}=\inf\{ \langle\Psi,\hat{H}\Psi\rangle: \langle\Psi,\Psi\rangle=1\}. \lb{V24} \ee
The practical significance of this result is that it can be used to formulate Rayleigh-Ritz variational 
calculations of the effective potential. This strategy has been employed in \cite{CJT} to calculate $V[\phi]$ 
by choosing a family of trial states $\{\Psi_\lambda\}$ and then varying the parameters $\lambda$ to optimize the guess. 
This yields an approximate result both for $V$ and also for the ground state wavefunction $\Omega.$ The important 
feature of such variational methods is that they are systematically improvable, even convergent to the true answer,
by suitably enlarging the class of trial states.

We shall show in the following sections that Symanzik's theorem, as well as various generalizations of it to the full, 
time-dependent action and to action functionals for higher-order statistics, can all be extended to the MSR context. 
Furthermore, this will allow us to elaborate Rayleigh-Ritz methods applicable to turbulence. The standard statements
and proofs are not automatically valid in MSR field theory because of the peculiar feature---noted by MSR themselves---that 
their formal quantum theory is one with a ``non-Hermitian Hamiltonian.'' However, we shall see that a simple modification 
suffices to make the results all valid there. 

\newpage

\noindent {\em (3.2) Operator Formulation of MSR Field Theory}

As already noted, the original paper of MSR \cite{MSR} developed the field-theory method by simply 
postulating an appropriate set of commutation relations between the usual variables and a new 
``response operator.'' A more concrete study was then made by Phythian \cite{Phy1,Phy2}, who constructed
canonical representations of the MSR commutation relations. It was also observed by him that the operator
approach seems to be limited to the case of {\em Markov dynamics}. We shall describe here briefly the 
operator formulation, following a similar approach as Phythian (but with a few differences). 

To this end, let us consider again the case of a generalized Langevin equation 
\be \partial_t x_n=K_n[\bx]+f_n'[\bx]. \lb{V32} \ee
where $\bK[\bx]$ is an arbitrary causal functional of histories $\bx(t)$ and $\bF'[\bx]$ is 
a zero-mean random force with generating functional
\be \int {\cal D}\bF'[\bx]\,e^{-i\langle \bp,\bF'[\bx]\rangle}
    =\exp\left\{\sum_{m\geq 2}{{(-i)^m}\over{m!}}\int p_{i_1}\cdots p_{i_m}D_{i_1\cdots i_m}^{(m)}[\bx]\right\}, 
\lb{V32a} \ee
where the cumulants ${\bf D}^{(m)}[\bx]$ are again causal functionals of $\bx.$ Then, by the same argument as 
given previously, this generalized Langevin dynamics has an equivalent field-theory formulation 
in terms of an action 
\begin{eqnarray}
S[\bx,\bp]& = & \int dt\,\,p_n(t)\left(\dot{x}_n(t)- K_n[t;\bx]\right)
             +{{i}\over{2}}\int dt\,\,{{\delta K_n[t;\bx]}\over{\delta x_n(t)}} \cr
        \, &   & \,\,\,\,\,\,+\sum_{m\geq 2}{{(-i)^{m+1}}\over{m!}}\int dt\,\,p_{i_1}(t)\cdots p_{i_m}(t)
                D_{i_1\cdots i_m}^{(m)}[t;\bx]
\lb{V33}
\end{eqnarray}
Note that
\be S[\bx,\bp] = \int (p_n\cdot\dot{x}_n- H[\bx,\bp]), \lb{V34} \ee
where
\be H[\bx,\bp]=p_n K_n[\bx]-\sum_{m\geq 2}{{(-i)^{m+1}}\over{m!}}p_{i_1}\cdots p_{i_m}D_{i_1\cdots i_m}^{(m)}[\bx].
\lb{V35} \ee 
For simplicity, we consider now and hereafter the case where a Liouville theorem is satisfied (or else Ito rule is
used.) The Euler-Lagrange equations of motion of this action are
\be \dot{x}_n(t)- K_n[t;\bx]+\sum_{m\geq 2}{{(-i)^{m+1}}\over{(m-1)!}}p_{i_1}(t)\cdots p_{i_{m-1}}(t)
                               D_{n i_1\cdots i_{m-1}}^{(m)}[t;\bx]=0, \lb{V36} \ee
and 
\be -\dot{p}_n(t)- \int dt'\,\,p_m(t'){{\delta K_m}\over{\delta x_n(t)}}[t';\bx]
                   +\sum_{m\geq 2}{{(-i)^{m+1}}\over{m!}}\int dt'\,\,p_{i_1}(t')\cdots p_{i_m}(t')
                               {{\delta D_{i_1\cdots i_m}^{(m)}}\over{\delta x_n(t)}}[t';\bx]=0. \lb{V37} \ee

Now assume that 
\be \bK[t;\bx]=\bK(\bx(t)), \,\,\,\,\,\,\,\,\,\,\,\,\,\,\bD^{(m)}[t;\bx]=\bD^{(m)}(\bx(t)), \lb{V38} \ee
that is, these quantities depend only upon the instantaneous values of the state variables. This restricts one
to the case of Markov dynamics, as previously noted by Phythian. In that case it is 
possible to restate the problem as a {\em formal quantum theory} with Hamiltonian operator 	
\be \hH=\hP_n K_n[\hbX]-\sum_{m\geq 2}{{(-i)^{m+1}}\over{m!}}\hP_{i_1}\cdots \hP_{i_m}D_{i_1\cdots i_m}^{(m)}[\hbX], 
\lb{V39} \ee
where the usual canonical commutation relations hold between $\hX_n,\hP_m$:
\be [\hX_n,\hP_m]=i\delta_{n,m}. \lb{V40} \ee
The equations of motion then follow in operator form as {\em Heisenberg equations}
\be i\partial_t \hat{O}=[\hat{O},\hH], \lb{V41} \ee
or 
\be \partial_t \hX_n=K_n(\hbX)-\sum_{m\geq 2}{{(-i)^{m+1}}\over{(m-1)!}}\hP_{i_1}\cdots \hP_{i_{m-1}}
                               D_{n i_1\cdots i_{m-1}}^{(m)}(\hbX), \lb{V42} \ee
and 
\be \partial_t \hP_n= -\hP_m {{\partial K_m}\over{\partial x_n}}[\hbX]
                   +\sum_{m\geq 2}{{(-i)^{m+1}}\over{m!}}\hP_{i_1}\cdots \hP_{i_m}
                               {{\partial D_{i_1\cdots i_m}^{(m)}}\over{\partial x_n}}(\hbX). \lb{V43} \ee
Observe that $\hH$ is not formally self-adjoint: it is just the ``non-Hermitian Hamiltonian'' of MSR. 
Actually, MSR considered instead operator equations of motion $\partial_t \hat{O}=[\hat{O},\hL]$ with
the trivially related operator $\hL= -i\hH.$ For a general Markov process this is just
\be \hL= -{{\partial}\over{\partial x_n}}\left[ K_n(\bx)\cdot\right]+\sum_{m\geq 2}{{(-1)^m}\over{m!}}
                   {{\partial^m}\over{\partial x_{i_1}\cdots \partial x_{i_m}}}\left[D_{i_1\cdots i_m}^{(m)}(\bx)
          \cdot\right] 
\lb{V44} \ee
in the canonical representation. It coincides with the well-known {\em Kramers-Moyal expansion} of the Markov evolution
operator \cite{Ris}. In particular, for the specific case of a {\em diffusion process} 
\be \hL= -{{\partial}\over{\partial x_n}}\left[ K_n(\bx)\cdot\right]
 +{{1}\over{2}}{{\partial^2}\over{\partial x_l\partial x_m}}\left[D_{lm}(\bx)\cdot\right] 
\lb{V45} \ee
and coincides with the {\em Fokker-Planck operator} \cite{Ris}, in which $\bK$ is the drift vector and 
$\bD$ the diffusion tensor. 

In the operator formulation, the statistical correlation and response functions are represented as
vacuum expectation values of the time-ordered products of the operators $\hbX$ and $\hbP.$ The latter are
usually called ``Greens functions'' in field-theory. However, since $\hH$ is 
non-Hermitian, the ``right vacuum state''
\be \hH|\Omega^+\rangle=0, \lb{V46} \ee
and ``left vacuum state''
\be \hH^*|\Omega^-\rangle=0 \lb{V47} \ee
are distinct. In fact, 
\be \langle\bx |\Omega^+\rangle= \rho(\bx), \lb{V48} \ee
where $\rho$ is the density of the stationary measure of the process and 
\be \langle\bx |\Omega^-\rangle\equiv 1. \lb{V49} \ee
If time-dependent ``Heisenberg picture'' operators are introduced as
\be \hat{O}(t)\equiv e^{it\hH}\hat{O}e^{-it\hH}, \lb{V50} \ee
then it is easy to see that 
\be \langle x_{i_1}(t_1)\cdots x_{i_n}(t_n)\rangle=
       \langle\Omega^-|{\rm T}\left[\hX_{i_1}(t_1)\cdots \hX_{i_n}(t_n)\right]|\Omega^+\rangle, \lb{V51} \ee
in which the LHS is the n-point statistical correlation function of the process and the RHS is 
the vacuum expectation of the time-ordered (increasing right to left) coordinate operators $\hX_i(t).$
Furthermore, the time-ordered vacuum expectations of the momentum operators $\hP_i(t)$ along with the
coordinate operators give the linear and higher-order mean response functions. The linear response function is 
defined as usual by  
\be G_{nm}(t-t')=\langle {{\delta x_n(t)}\over{\delta f_m(t')}}\rangle, \lb{V52} \ee
where $f_n(t)$ is an external, deterministic force coupled additively into the Langevin dynamics. We shall just
show here that indeed 
\be G_{nm}(t)= -i\langle\Omega^-|{\rm T}\left[\hX_n(t)\hP_m(0)\right]|\Omega^+\rangle \lb{V53} \ee
and we leave the consideration of higher-order response functions to the reader. First we observe
that due to 
\be \hP_m|\Omega^-\rangle=0, \lb{V54} \ee
any Greens function whose latest time-variable is a momentum operator must vanish. Hence, the RHS of 
Eq.(\ref{V53}) vanishes for $t<0$ just as does the causal response function on the LHS. Furthermore, by
means of the canonical commutation relations
\begin{eqnarray}
 -i\hX_n(t)\hP_m(0) & = & i[\hP_m(0),\hX_n(t)] -i\hP_m(0)\hX_n(t) \cr
                      \,& = & {{\partial \hX_n(t)}\over{\partial x_m(0)}}-i\hP_m(0)\cdot \hX_n(t). \lb{V55} 
\end{eqnarray}
Thus, for $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|<K\}$
or the square box ${\cal K}=\{\bk: \max_{i=1}^d|k_i|<K\}.$ Note that $\bP(\grad)$ includes the projection
back onto the ${\cal K}$-modes. In our proof it is helpful to use the ``cut-off delta function'' 
\be \delta^d(\bx)=\sum_{\bk\in {\cal K}}e^{i\bk\bdot\bx}, \lb{72aa} \ee
which was discussed in Section 2 of \cite{ZM}. As far as integration with respect to other cutoff functions,
it has the same properties as the ordinary delta-function. Our proof will make use of symmetry properties 
which are assumed to hold for the cutoff domain ${\cal K}$, particularly,
\be {\cal K}=-{\cal K} \lb{73aa} \ee
and invariance under a discrete rotation group mapping any coordinate direction into another. 
Of course, these are satisfied by either of the explicit choices above, as well as other reasonable choices. 

It is then easy to calculate, using incompressibility, that 
\be \partial_{x_i}{{\delta\tau^E_{ij}(x)}\over{\delta v_k(y)}}= (\partial_{x_i}\delta^{d+1})(x-y)
     [\delta_{ik}v_j(x)+\delta_{jk}v_i(x)]+\delta^{d+1}(x-y)\partial_{x_k}v_j(x). \lb{73a}\ee
Then, again using incompressibility, 
\be {\rm tr}\left({{\delta\bK^E(x)}\over{\delta\bv(y)}}\right)= 2(d-1)(\bv(x)\cdot\grad_\bx)\delta^{d+1}(x-y)
                     +\bP(\grad_\bx)\delta^{d+1}(x-y):(\grad_\bx\bv)(\bx) \lb{74a} \ee
is obtained. Note that the factor $(d-1)$ arises from the trace of the solenoidal projection. However, it is observed that
\be \grad_\bx\delta^d(\bz)=\sum_{\bk\in {\cal K}}\bk=\bz \lb{75a} \ee
by reflection symmetry of the cutoff domain ${\cal K}$ in Fourier space. Likewise, 
\be \bP(\grad_\bx)_{ij}\delta(\bz)= C\delta_{ij} \lb{80a} \ee
if the reflection and discrete rotation symmetries hold. This gives finally that 
\be {\rm Tr}\left({{\delta\bK^E}\over{\delta\bv}}\right)\equiv \int d^{d+1}x\,\,
                      {\rm tr}\left({{\delta\bK^E(x)}\over{\delta\bv(x)}}\right) =0. \,\,\,\,\,\,\,\,\Box
\lb{81a} \ee

It follows immediately from this that a corresponding result holds for the coupled LE and SE obtained by filtering
the Navier-Stokes equation
\be {\rm tr}\left({{\delta\Kl^E(x)}\over{\delta\vl(x)}}\right)
     ={\rm tr}\left({{\delta\bK^{E\prime}(x)}\over{\delta\vs(y)}}\right)=0. \lb{82a} \ee
Note that the LE and SE equations are regarded as a system for {\em autonomous} variables $\vl$ and $\vs.$
Consider the first case. Since
\be \overline{K}_i^E[x;\vl,\vs]=\int d^d\by\,\,G(\bx-\by)K_i^E[\by;\vl+\vs], \lb{83a} \ee
it follows that in the equation corresponding to Eq.(\ref{74a}), $\bv$ is replaced everywhere by $\vl+\vs$ 
and $\delta^{d+1}(x)$ is replaced by $G(\bx)\delta(t)$:
\begin{eqnarray}
{\rm tr}\left({{\delta\Kl^E(x)}\over{\delta\vl(y)}}\right)& = & 2(d-1)[(\vl(x)+\vs(x))\cdot\grad_\bx]G(\bx-\by)\delta(t-s)\cr
      \,& &\,\,\,\,\,\,\,\,\,\,\,\,+\bP(\grad_\bx)G(\bx-\by)\delta(t-s):[(\grad_\bx\vl)(\bx)+(\grad_\bx\vs)(\bx)]. \lb{84a} 
\end{eqnarray}
If one assumes that $\hat{G}$ has the necessary symmetries---invariance under reflections and at least discrete rotation
invariance---then $\grad G(\bz)=\bz$ and $P_{ij}(\grad)G(\bz)=C\cdot\delta_{ij}.$ These facts imply that the first part
of Eq.(\ref{82a}) holds. The demonstration of the second part is the same, simply replacing $G$ by $H.$

We now show that the Liouville theorem for the SS contribution to $\Kl_{\rm eff},$
\be {\rm Tr}\left({{\delta\bD^{(1)}_{\rm eff}}\over{\delta\vl}}\right)=0, \lb{76a} \ee
follows from a causality argument. In fact, note that in the term $\delta\bD^{(1)}_{\rm eff}[x;\vl]/\delta\vl(y)$ 
(or in the similar functional derivative of $\btaul[x;\bv]$) there must be a line of response functions leading 
from the $y$-vertex to the $x$-vertex. Otherwise, starting with the $y$-vertex one must encounter a closed loop of 
response lines according to our argument in Section 2.2, and the graph will vanish. However, this implies that
\be \left.{{\delta\taul_{lk}[x;\bv]}\over{\delta\bv_j(y)}}\right|_{t=s}=\delta^d(\bx-\by)B_{jkl}[y;\bv], \lb{77a} \ee
since the response functions at equal time arguments are spatial delta-functions. In that case, it follows from 
Eq.(\ref{68}) in Section 2.2 that
\begin{eqnarray} 
\left.{{\bD^{(1)}_{{\rm eff},i}[x;\bv]}\over{\delta v_j(y)}}\right|_{t=s} & = & -P_{il}(\grad_\bx)(\grad_\bx)_k
                                                                  \left.{{\delta\taul_{lk}[x;\bv]}\over{\delta\bv_j(y)}}\right|_{t=s} \cr
                                                               \,& = & -\left. P_{il}(\grad_\bx)(\grad_\bx)_k\delta^d(\bx-\by)
                                                                          \cdot B_{jkl}[y;\bv]\right|_{t=s}. \lb{78a} 
\end{eqnarray}
Therefore, 
\begin{eqnarray} 
{{\bD^{(1)}_{{\rm eff},i}[x;\bv]}\over{\delta v_j(x)}} & = & -P_{il}(\grad)(\grad)_k\delta^d(\bz) \cdot B_{jkl}[x;\bv] \cr
                                         \, & = & 0, \lb{79a} 
\end{eqnarray}
by our previous formula Eq.(\ref{75a}) from reflection symmetry. This gives at once the Liouville theorem (Eq.(\ref{76}).) 
                         
This result may seem paradoxical in view of the well-known term in the effective dynamics associated to the 
``eddy-viscosity.'' Such a term has the same form as the usual viscous dynamics in Navier-Stokes (from molecular viscosity)
and apparently violates the Liouville theorem! However, it should be clear that this violation results from
various approximate evaluations of the exact expressions for the eddy-damping. In fact, the 1-loop perturbation diagram
for the eddy-damping term---from which eddy-viscosity contributions are ordinarily derived---is easily seen to 
satisfy the Liouville theorem according to the argument given here. Note that a distinction must be made between
``dissipative dynamics'', which in the current usage of dynamical systems means those without a Liouville theorem,
and ``non-conservative dynamics,'' designating those for which conservation laws, such as energy conservation, are violated.
There is no question that the eddy-damping term violates energy conservation of the LS modes and that its principal effect 
at low wavenumbers is to drain energy from those modes.

\vspace{.4in}

\noindent {\bf Acknowledgements.} This work has grown out of discussions with a large number of scientists over a 
period of years, including D. Carati, U. Frisch, N. Goldenfeld, R. H. Kraichnan, J. Krug, J. L. Lebowitz, P. C. Hohenberg,
Y. Oono, R. Maier, U. Piomelli, H. Rose, Z.-S. She, H. Spohn, D. Stein, V. Yakhot, and X. Zeng. I am grateful to
all for their insights and information. 

\newpage
                                          
\begin{thebibliography}{99}
\bibitem[1]{Kr76} R. H. Kraichnan, ``Eddy viscosity in two and three dimensions,'' J. Atmos. Sci.
               {\bf 33} 1521 (1976). 
\bibitem[2]{Kr87} R. H. Kraichnan, ``Eddy viscosity and diffusivity: exact formulas and approximations ,'' 
               Complex Systems {\bf 1} 805 (1987).
\bibitem[3]{LL} L. D. Landau and E. M. Lifshitz, {\em Fluid Mechanics}. (Pergamon Press, New York, 1959),
               Chapter 17.
\bibitem[4]{Rose} H. A. Rose, ``Eddy diffusivity, eddy noise, and sub-grid scale modelling,'' J. Fluid. Mech.
               {\bf 81} 719 (1977). 
\bibitem[5]{Lei} C. E. Leith, ``Stochastic backscatter in a subgrid-scale model: plane shear mixing layer,''
               Phys. Fluids. A {\bf 2} 297 (1990).
\bibitem[6]{Chas} J. R. Chasnov, ``Simulation of the Kolmogorov inertial subrange using an improved subgrid 
               model,'' Phys. Fluids A {\bf 3} (1991). 
\bibitem[7]{MT}P. J. Mason and D. J. Thomson, ``Stochastic backscatter in large-eddy simulations of boundary
               layers,'' J. Fluid Mech. {\bf 242} 51 (1992).
\bibitem[8]{Car} D. Carati, ``Iterative filtering of the forced Navier-Stokes equation,'' preprint (1994). 
\bibitem[9]{YO} V. Yakhot and S. A. Orszag, ``Renormalization group analysis of turbulence, I. basic theory,'' 
               J. Sci. Comp. {\bf 1} 3 (1986).
\bibitem[10]{Kr85} R. H. Kraichnan, ``Decimated amplitude equations in turbulence dynamics,'' in:
                 {\em Theoretical Approaches to Turbulence}, ed. D. L. Dwoyer, M. Y. Hussaini,
                 and R. G. Voigt (Springer, New York, 1985). 
\bibitem[11]{MSR} P. C. Martin, E. D. Siggia, and H. A. Rose, ``Statistical dynamics of classical systems,'' 
               Phys. Rev. A {\bf 8} 423 (1973).
\bibitem[12]{O} L. Onsager, ``Reciprocal relations in irreversible processes. I,II,'' 
                 Phys. Rev. {\bf 37} 405 (1931); {\bf 38} 2265 (1931).
\bibitem[13]{I-LES} G. L. Eyink, ``Large-eddy simulation and the `multifractal model' of turbulence:
                 {\em a priori} estimates on subgrid flux and locality of energy transfer,'' 
                 submitted to Phys. Fluids (1994). 
\bibitem[14]{Ger1}M. Germano, ``Turbulence: the filtering approach,'' J. Fluid. Mech. {\bf 238} 325 (1992).
\bibitem[15]{Ger2} M. Germano, ``A proposal for a redefinition of the turbulent stresses in the filtered 
                   Navier-Stokes equations,'' Phys. Fluids {\bf 29} 2323 (1986).
\bibitem[16]{Leo}A. Leonard, ``On the energy cascade in large-eddy simulations of turbulent flows,''
                 Adv. Geophys. {\bf 18A} 237 (1974).
\bibitem[17]{Rue1} D. Ruelle, ``Measures describing a turbulent flow,'' Ann. N. Y. Acad. Sci. {\bf 357} 1 (1980).
\bibitem[18]{ER}J.-P. Eckmann and D. Ruelle, ``Ergodic theory of chaos and strange attractors,'' 
                Rev. Mod. Phys. {\bf 57} 617 (1985).
\bibitem[19]{Rue2} D. Ruelle, ``Microscopic fluctuations and turbulence,'' Phys. Lett. A {\bf 72} 81 (1979).
\bibitem[20]{HoSh}P. C. Hohenberg and B. Shraiman, ``Chaotic behavior of an extended system,'' 
                Physica D {\bf 37} 109 (1989). 
\bibitem[21]{Yak}V. Yakhot, ``Large-scale properties of unstable systems governed by the Kuramoto-Sivashinsky 
                 equation,'' Phys. Rev. A {\bf 24} 642 (1983).
\bibitem[22]{Zal}S. Zaleski, ``A stochastic model for the large scale dynamics of some fluctuating interfaces,''
                 Physica D {\bf 34} 427 (1989).
\bibitem[23]{SKJJB}K. Sneppen et al., ``Dynamic scaling and crossover analysis for the Kuramoto-Sivashinsky 
                   equation,'' Phys. Rev. A {\bf 46} R7351 (1992).
\bibitem[24]{Lor}E. Lorenz, ``Maximum simplification of the dynamic equations,'' Tellus {\bf 12} 243 (1960). 
\bibitem[25]{Phy1}R. Phythian, ``The operator formalism of classical statistical dynamics,''
                 J. Phys. A {\bf 8} 1423 (1975).
\bibitem[26]{Phy2}R. Phythian, ``Further application of the Martin, Siggia, Rose formalism,''
                 J. Phys. A {\bf 9} 269 (1976).
\bibitem[27]{Jan}H. K. Janssen, ``On a Lagrangean for classical field dynamics and renormalization group calculations of
                dynamical critical properties,'' Z. Phys. B {\bf 23} 377 (1976).
\bibitem[28]{DeD}C. DeDominicis, ``Techniques de renormalisation de la th\'{e}orie des champs et dynamique des
                ph\'{e}nomen\`{e}s critiques,''J. Phys. (Paris) C {\bf 1} 247 (1976).
\bibitem[29]{Gr1}R. Graham, ``Path integral formulation of general diffusion processes,'' Z. Phys. B {\bf 26} 281 (1977).
\bibitem[30]{Gr2}R. Graham, ``Short-time propagators in Riemannian geometries,'' in: {\em Statphys 13}, Proceedings of the 
                 13th IUPAP Conference on Statistical Physics, Haifa, 1977.
\bibitem[31]{Wei}U. Weiss, ``Operator ordering schemes and covariant path integrals of quantum and stochastic 
                 processes in curved space,'' Z. Phys. B {\bf 30} 429 (1978).
\bibitem[32]{Wis}C. Wissel, ``Manifolds of equivalent path integral solutions of the Fokker-Planck equations,''
                 Z. Physik B {\bf 35} 185 (1979).
\bibitem[33]{Arn}L. Arnold, {\em Stochastic Differential Equations: Theory and Applications.}
                 (Wiley, New York, 1974).
\bibitem[34]{Lee}T. D. Lee, ``On some statistical properties of the hydrodynamical and magneto-hydrodynamical fields,''
                 Q. Appl. Math. {\bf 10} 69 (1952).
\bibitem[35]{Wyl}H. W. Wyld, ``Formulation of the theory of turbulence in an incompressible fluid,'' Ann. Phys. (N. Y.) 
                {\bf 14} 143 (1961).
\bibitem[36]{Kr61}R. H. Kraichnan, ``Dynamics of nonlinear stochastic systems,'' J. Math. Phys. {\bf 2} 124 (1961).
\bibitem[37]{ZL}V. E. Zakharov and V. S. L'vov, ``Statistical description of nonlinear wave fields,'' Izv. Vyss. Uch. Zav., 
                Radiofiz. {\bf 18} 1470 (1975).
\bibitem[38]{Kr77}R. H. Kraichnan, ``Eulerian and Lagrangian renormalization in turbulence theory,'' J. Fluid Mech. 
                {\bf 83} 349 (1977).
\bibitem[39]{I-RNG}G. L. Eyink, ``The renormalization group method in statistical hydrodynamics,'' Phys. Fluids {\bf 6}
                   3063 (1994).
\bibitem[40 ]{LLPP}V. S. L'vov et al., ``Proof of scale invariant solutions in the Kardar-Parisi-Zhang and 
                 Kuramoto-Sivashinsky equations in $1+1$ dimensions: analytical and numerical results,''
                 Nonlinearity {\bf 6} 25 (1993). 
\bibitem[41]{Sym}K. Symanzik, ``Renormalizable models with simple symmetry breaking,'' Commun. Math. Phys. {\bf 16} 48
                (1970).
\bibitem[42]{HE} W. Heisenberg and H. Euler, ``Folgerungen aus der Diracshen Theorie des Positrons,'' Z. Physik {\bf 98} 
                 714 (1936).
\bibitem[43]{Sch}J. Schwinger, ``On gauge invariance and vacuum polarization,'' Phys. Rev. {\bf 82} 664 (1951). 
\bibitem[44]{OM}L. Onsager and S. Machlup, ``Fluctuations and irreversible processes,'' Phys. Rev. {\bf 91}
                1505 (1953).
\bibitem[45]{Gr3}R. Graham, ``Path-integral methods in nonequilibrium thermodynamics and statistics,'' in: 
                {\em Stochastic Processes in Nonequilibrium Systems.} Sitges School, June 1978, eds. L. Garrido,
                P. Seglar, and P. J. Shepherd. Lecture Notes in Physics, vol. 84. (Springer-Verlag, Berlin, 1978). 
\bibitem[46]{I-LD}G. L. Eyink, ``Dissipation and large thermodynamic fluctuations,'' J. Stat. Phys. {\bf 61} 533 (1990).
\bibitem[47]{CJT}J. M. Cornwall, R. Jackiw, and E. Tomboulis, ``Effective action for composite operators,''
                 Phys. Rev. D {\bf 10} 2428 (1974).
\bibitem[48]{Ris}H. Risken, {\em The Fokker-Planck Equation.} (Springer-Verlag, Berlin, 1984).
\bibitem[49]{HS}M. Hausner and J. T. Schwartz, {\em Lie Groups and Their Lie Algebras}. (New York, Gordon and 
                Breach, 1968).
\bibitem[50]{Rue3}D. Ruelle, ``Locating resonances for Axiom A dynamical systems,'' J. Stat. Phys. {\bf 44} 281 (1986).
\bibitem[51]{JK}R. Jackiw and A. Kerman, ``Time-dependent variational principle and effective action,''
                Phys. Lett. A {\bf 71} 158 (1979).
\bibitem[52]{Dir}P. A. M. Dirac, ``Note on exchange phenomena in the Thomas atom,'' Proc. Camb. Philo. Soc. {\bf 26} 
                 (1930) 376.
\bibitem[53]{Kr70}R. H. Kraichnan, ``Convergents to turbulence functions,'' J. Fluid Mech. {\bf 41} 189 (1970).
\bibitem[54]{Kr58}R. H. Kraichnan, ``Irreversible statistical mechanics of incompressible hydromagnetic 
                  turbulence,'' Phys. Rev. {\bf 109} 1407 (1958).
\bibitem[55]{Kr79}R. H. Kraichnan, ``Variational method in turbulence theory,'' Phys. Rev. Lett. {\bf 42} 1263 (1979).
\bibitem[56]{Qi}J. Qian, ``Variational approach to the closure problem of turbulence theory,'' Phys. Fluids
                {\bf 26} 2098 (1983). 
\bibitem[57]{Cas}B. Castaing, ``Cons\'{e}quences d'un principe d'extremum en turbulence,'' J. de Phys. (Paris)
                 {\bf 50} 147 (1989).
\bibitem[58]{Bus}F. H. Busse, ``The optimum theory of turbulence,'' Adv. Appl. Math. {\bf 18} 77 (1978).
\bibitem[59]{HKr} J. R. Herring and R. H. Kraichnan, ``Comparison of some approximations for isotropic 
                 turbulence,'' in: {\em Statistical Models and Turbulence}, ed. M. Rosenblatt and 
                 C. Van Atta (Springer, New York, 1972). 
\bibitem[60]{MC}C. Meneveau and A. Chhabra, ``Two-point statistics of multifractal measures,'' 
                Physica A {\bf 164} 564 (1990).
\bibitem[61]{LKr}C. E. Leith and R. H. Kraichnan, ``Predictability of turbulent flows,'' J. Atmos. Sci.
                 {\bf 29} 1041 (1972).
\bibitem[62]{Lor2}E. N. Lorenz, ``Deterministic non-periodic flow,'' J. Atmos. Sci. {\bf 20} 130 (1963). 
\bibitem[63]{PNNWZ}R. A. Pielke et al. ``Several unresolved issues in numerical modelling of geophysical 
                   flows,'' submitted to Atmosphere-Oceans. 
\bibitem[64]{ZM}D. N. Zubarev and V. G. Morozov, ``Statistical mechanics of nonlinear hydrodynamic fluctuations,''
                Physica {\bf 120A} 411 (1983). 

\end{thebibliography}

\end{document}

\end
ENDBODY
