%Yuri Lvov 8/27/3

%PaperId: math-ph/0305028, PaperPassword: pw9md (access no longer password restricted)


\documentstyle[aps,prl,multicol,epsf]{revtex}

\topmargin 0.05in\textheight 8in

\def \BE {\begin{equation}}
\def \EE {\end{equation}}
\def \BEA {\begin{eqnarray}}
\def \EEA {\end{eqnarray}}
\def \CR {\nonumber \\}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%% START OF PAPER %%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{document}
\title{
 Noisy spectra, long correlations, and intermittency in wave turbulence\\
{\bf Phys. Rev. E 69, 066608 (2004) }}
\author{Yuri V. Lvov$^\dagger$ and Sergey Nazarenko$^*$ 
}
\address
{ 
$^\dagger$ 
Department of Mathematical Sciences, Rensselaer Polytechnic Institute,
Troy, NY 12180 \\ 
$^*$
Mathematics Institute, The University of Warwick,  Coventry, CV4-7AL, UK}

\maketitle

\pacs 04.30.Nk , 24.60.-k, 52.35.Ra, 92.10.Cg, 87.15.Ya 

 

%04.30.Nk Wave propagation and interactions

%24.60.-k Statistical theory and fluctuations

%52.35.Ra Plasma turbulence

%92.10.Cg Capillary waves 

%87.15.Ya Fluctuations 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%% ABSTRACT %%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{abstract}
We study the k-space fluctuations of the waveaction about its mean
spectrum in the turbulence of dispersive waves. We use a minimal model
based on the Random Phase Approximation (RPA) and derive evolution
equations for the arbitrary-order one-point moments of the wave
intensity in the wavenumber space. The first equation in this series
is the familiar Kinetic Equation for the mean waveaction spectrum,
whereas the second and higher equations describe the fluctuations
about this mean spectrum. The fluctuations exhibit a nontrivial
dynamics if some long coordinate-space correlations are present in the
system, as it is the case in typical numerical and laboratory
experiments. Without such long-range correlations, the fluctuations
are trivially fixed at their Gaussian values and cannot evolve even if
the wavefield itself is non-Gaussian in the coordinate space.  Unlike
the previous approaches based on smooth initial k-space cumulants, our
approach works even for extreme cases where the k-space fluctuations
are absent or very large and intermittent.  We show, however, that
whenever turbulence approaches a stationary state, all the moments
approach the Gaussian values.
\end{abstract}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%% MAIN DOCUMENT %%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{multicols}{2}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

The concept of Wave Turbulence (WT), which describes an ensemble of
weakly interacting dispersive waves, significantly enhanced our
understanding of the spectral energy transfer in complex systems like
the ocean, the atmosphere, or in
plasmas~\cite{ZLF,Ben,GS,Zakfil,hasselman}.  This theory also became a
subject of renewed interest recently, (see,
e.g. \cite{OsbornePRL,ZakharovPRL,SoomerePRL,LT}).  Traditionally, WT
theory deals with derivation and solutions of the Kinetic Equation
(KE) for the mean waveaction spectrum (see e.g. \cite{ZLF}).  However,
all experimentally or numerically obtained spectra are ``noisy'',
i.e. exhibit k-space fluctuations which contain a complimentary to the
mean spectra information.  Such fluctuations will be studied (for the
first time) in this manuscript.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{ Random phases vs Gaussian fields --- } The random phase
approximation (RPA) has been popular in WT because it allows a quick
derivation of KE \cite{ZLF,GS}.  We will use RPA in this paper because
it provides a minimal model for for description of the k-space
fluctuations of the waveaction about its mean spectrum, but we will
also discuss relation to the approach of  \cite{Ben} which does not
assume RPA.  By definition, RPA for an ensemble of complex fields $a_k
= A_k e^{i\phi_k} $ means that the phases $\phi_k$ are uniformly
distributed in $(0,2\pi]$ and are statistically independent of each
other and of the amplitude $A_k$, $\langle \phi_{k_1} \phi_{k_2}
\rangle = \pi \, \delta^1_2, \;\; \langle \phi_{k_1} A_{k_2} \rangle =
0$\footnote{ We start by considering fields in a periodic box which is
an essential intermediate step in the definition of RPA and the new
correlators $M^{(p)}_k$ introduced later in this work.  Therefore
$\delta^1_2$ here is the Kronecker symbol. Later, we take the large
box limit corresponding to homogeneous wave turbulence.}.  Thus, the
averaging over the phase and over the amplitude statistics can be
performed independently. In RPA, the fluctuations of the amplitudes
$A_k$ must also be decorrelated at different $k$'s\footnote{This
property is typically not mentioned explicitly (but used implicitly)
when RPA is employed. }, $\langle A_{k_1}^n A_{k_2}^m \rangle =\langle
A_{k_1}^n \rangle \langle A_{k_2}^m \rangle \;\; (m,n = 1,2,3,...)$.

To illustrate the relation between the random phases and Gaussianity,
let us consider the fourth-order moment for which RPA gives
%
\begin{equation}
\langle a_{k_1} a_{k_2} {\overline a_{k_3}} {\overline a_{k_4}}
\rangle = n_{k_1} n_{k_2} (\delta^{1}_{3} \delta^{2}_{4} + 
\delta^{1}_{4}\delta^2_3) + Q_{k_1} \delta^1_{2}
\delta^{1}_{3} \delta^{1}_{4},
\label{wick}
\end{equation}
%
where $ n_k = \langle A_{k}^2 \rangle $ is the waveaction spectrum and
$ Q_k = \langle A_{k}^4 \rangle $ is a cumulants coefficient.  The
last term in this expression appears because the phases drop out for
$k_1=k_2=k_3=k_4 $ and their statistics poses no restriction on the
value of this correlator at this point. This cumulant part of the
correlator can be arbitrary for a general random-phased field whereas
for Gaussian fields $ Q_k$ must be zero. Such a difference between the
Gaussian and the random-phased fields occurs only at a vanishingly
small set of modes with $k_1=k_2=k_3=k_4 $ and it has been typically
ignored before because its contribution to KE is negligible.  
Therefore, if the mean waveaction spectrum was the only thing we were
interested in, we could safely ignore contributions from all
(one-point) moments $M^{(p)}_k = \langle |a_{k}|^{2p} \rangle \;\;
(p=1,2,3,..)$.

However, it is precisely moments $M^{(p)}_k$ that contain information
about fluctuations of the waveaction about its mean spectrum. For
example, the standard deviation of the waveaction from its mean is
%
$\sigma_k = (\langle |a_{k}|^4 \rangle - \langle |a_{k}|^2
\rangle^2)^{1/2} = (M^{(2)}_k - n_k^2)^{1/2} 
%= (n_k^2 + Q_k)^{1/2}
$.
%
This quantity can be arbitrary for a general random-phased field
whereas for a Gaussian wave field the fluctuation level $\sigma_k$ is
fixed, $\sigma_k = n_k$.  Note that different values of moments
$M^{(p)}_k$ can correspond to hugely different typical wave field
realizations.  In particular, if $M^{(p)} = n^p$ then there is no
fluctuations and $A_k$ is deterministic, $\sigma_k=0$.  For the
opposite extreme of large fluctuations we would have $M^{(p)} \gg n^p$
which means that the typical realization is sparse in the k-space and
is characterized by few intermittent peaks of $A_k$ and close to zero
values in between these peaks. Note that the information about the
spectral fluctuations of the waveaction contained in the one-point
moments $M^{(p)}$ is completely erased from the multiple-point moments
by the random phases and it is precisely why these new objects play a
crucial role for the description of the fluctuations.

Will the waveaction fluctuations appear if they were absent initially?
Will they saturate at the Gaussian level $\sigma_k = n_k$ or will they
keep growing leading to the k-space intermittency? To answer these
questions, we will use RPA to derive and analyze equations for the
moments $M^{(p)}_k$ for arbitrary orders $p$ and thereby describe the
statistical evolution of the spectral fluctuations.  Note that RPA,
without a stronger Gaussianity assumption, is totally sufficient for
the WT closure at any order.  This allows us to study wavefields with
moments $M^{(p)}_k$ very far from their Gaussian values, which may
happen, for example, because of the choice of initial conditions or a
non-Gaussianity of the energy source in the system.

In \cite{Ben} non-Gaussian fields of a rather different kind were
considered.  Namely, statistically homogeneous wave fields were
considered in an infinite space which initially have decaying
correlations in the coordinate space and, therefore, smooth cumulants
in the k-space, e.g.
%
$$
\langle a_{k_1} a_{k_2} {\overline a_{k_3}} {\overline a_{k_4}}
\rangle = n_{k_1} n_{k_2} (\delta^{k_1}_{k_3} \delta^{k_2}_{k_4} + 
\delta^{k_1}_{k_4}\delta^2_3) + C_{123} \delta^{k_1+k_2}_{k_3 + k_4}, $$
%
where $C_{123}$ is a smooth function of $k_1, k_2, k_3$ and
$\delta$'s now mean Dirac deltas. On the other hand, by taking the
large box limit it is easy to see that our expression (\ref{wick})
corresponds to a {\em singular} cumulant $C_{123} = Q_{k_1}/{\cal V}
\, \delta^{k_1}_{k_2} \delta^{k_1}_{k_3}$. It tends to zero when the
box volume ${\cal V}$ tends to infinity and yet it gives a finite
contribution to the waveaction fluctuations in this limit.\footnote{
Thus, assuming a finite box is an important intermediate step when
introducing the relevant to the fluctuations objects like $Q_k$.}
This singular cumulant corresponds to a small component of the
wavefield which is long-correlated - the case not covered by the
approach of \cite{Ben}. On the other hand, it would be
straightforward to go beyond our RPA by adding a cumulant part of
the initial fields which tends to a smooth function of $k_1, k_2,
k_3$ in the infinite box limit (like in \cite{Ben}). However, such
cumulants would give a box-size dependent contribution to the
waveaction fluctuations which vanishes in the infinite box limit
(e.g.  it would change $\sigma_k^2$ by $C_{kkk}/{\cal V}$). Thus, in
large boxes the waveaction fluctuation for the fields with smooth
cumulants is fixed at the same value as the for Gaussian fields,
$\sigma_k = n_k$, and introduction of the singular cumulant is
essential to remove this restriction on the level of fluctuations.
On the other hand, the smooth part of the cumulants has no bearing on
the closure, as shown in~\cite{Ben} and on the large-box
fluctuation and, therefore, will be omitted here  for
brevity and clarity of the analysis.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{  Time-scale separation analysis --- }
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Consider weakly nonlinear dispersive waves in a periodic box.  
Here  we consider quadratic nonlinearity and the linear dispersion
relations $\omega_k$ which allow three-wave interactions. Example of
such systems include surface capillary waves~\cite{Zakfil} and
internal waves in the ocean~\cite{LT}.  In Fourier space, the general
form for the Hamiltonian systems with quadratic nonlinearity looks as
follows,\footnote{We will follow the RPA approach as presented by
Galeev and Sagdeev \cite{GS} but deal with a slightly more general
case where the wave field is not restricted by the condition
$\overline a(k) = a(-k)$. We will also use elements of the technique
and notations of \cite{Ben}. }
%
\BEA
{\cal H} &=& \sum_{n=1}^\infty \omega_n|c_n|^2 +\epsilon
\sum_{l,m,n=1}^\infty \left( 
V^l_{mn} \bar c_{l} c_m  c_n\delta^l_{m+n}+c.c.\right),\CR
i\dot c_l &=&\frac{\partial {\cal H}}{\partial \bar c_l}, \ \ 
c_l=a_l e^{-i \omega_l t}, \CR
i \, \dot a_l &=& \epsilon \sum_{m,n=1}^\infty \left( V^l_{mn} a_{m}
a_{n}e^{i\omega_{mn}^l t} \, \delta^l_{m+n} 
\right.\CR && \left. \hspace{3cm}
+ 2 \bar{V}^{m}_{ln} \bar a_{n}
a_{m} e^{-i\omega^m_{ln}t } \, \delta^m_{l+n}\right),
\label{Interaction} \EEA 
%
where $a_n=a(k_n)$ is the complex wave amplitude in the interaction representation, 
$k_n = 2 \pi n/L $, $L $ is the box side length,
$n=(n_1,n_2)$ for 2D, or $ n=(n_1,n_2, n_3)$ in 3D, (similar for $k_l
$ and $ k_m$), $
\omega^l_{mn}\equiv\omega_{k_l}-\omega_{k_m}-\omega_{k_m}$ and
$\omega_l=\omega_{k_l}$ is the wave linear dispersion relation.  Here,
$V^l_{mn} \sim 1$ is an interaction coefficient and $\epsilon$ is
introduced as a formal small nonlinearity parameter.

In order to filter out fast oscillations at the wave period, let us
seek for the solution at time $T$ such that $2 \pi / \omega \ll T \ll
1/\omega \epsilon^2$.  The second condition ensures that $T$ is a lot
less than the nonlinear evolution time.  Now let us use a perturbation
expansion in small $\epsilon$,
%
$$a_l(T)=a_l^{(0)}+\epsilon a_l^{(1)}+\epsilon^2 a_l^{(2)}.$$
%
Substituting this expansion in (\ref{Interaction}) we get in
the zeroth order
%
$ a_l^{(0)}(T)=a_l(0)\label{definitionofa} $,
%
i.e. the zeroth order term is time independent. This corresponds to
the fact that the interaction representation wave amplitudes are
constant in the linear approximation.  For simplicity, we will write
$a^{(0)}_l(0)= a_l$, understanding that a quantity is taken at $T=0$
if its time argument is not mentioned explicitly.  The first order is
given by
% 
\BEA
a^{(1)}_l (T) = -i \sum_{m,n=1}^\infty \left(   V^l_{mn}
a_m a_n \Delta^l_{mn} \delta^l_{m+n}\right.\CR\left.\hskip 4cm
+
2 \bar{V}^m_{ln}a_m\bar{a}_n \bar\Delta^m_{ln}\delta^m_{l+n}
\right),
\label{FirstIterate}
\EEA
%
where 
$ \Delta^l_{mn}=\int_0^T e^{i\omega^l_{mn}t}d t =
({e^{i\omega^l_{mn}T}-1})/{i \omega^l_{mn}}. \label{NewellsDelta}
$
%
Here we have taken into account that $a^{(0)}_l(T)= a_l$ and
$a^{(1)}_k (0)=0$.  

To calculate the second iterate, write
%
\BEA
i\dot{ a}^{(2)}_l  = \sum_{m,n=1}^\infty \Big[  
V^l_{mn}\delta^l_{m+n} e^{i \omega^l_{mn} t}
\left(a_m^{(0)} a_n^{(1)}+ a_m^{(1)} a_n^{(0)}\right)
\CR \hspace{3cm}+
2\bar{V}^m_{ln}\delta^m_{l+n} e^{ -i \omega^m_{ln} t}
\left(a_m^{(1)} \bar{a}_n^{(0)}+ a_m^{(0)} \bar{a}_n^{(1)}\right)
\Big].\CR
\label{SecondIterateTimeDerivative}
\EEA 
%
We now have to substitute (\ref{FirstIterate}) into
(\ref{SecondIterateTimeDerivative}) and integrate over time to obtain
\end{multicols}
\leftline{-------------------------------------------------------------------------}
%
\BEA a_l^{(2)} (T)  &=& \sum_{m,n, \mu, \nu=1}^\infty \left[ 2 V^l_{mn} 
\left(
-V^m_{\mu \nu}a_n a_\mu a_\nu E[\omega^l_{n \mu \nu},\omega^l_{mn}]
\delta^m_{\mu + \nu} 
-2 
\bar V^\mu_{m \nu}a_n a_\mu \bar a_\nu \bar 
E[\omega^{l \nu}_{n \mu},\omega^l_{mn}]\delta^\mu_{m + \nu}\right)
\delta^l_{m+n} \right.\CR && \left.
+ 2 
\bar V^m_{ln} 
 \left(
-V^m_{\mu \nu}\bar a_n a_\mu a_\nu E[\omega^{ln}_{\mu \nu},-\omega^m_{ln}]
\delta^m_{\mu + \nu} 
-
2 \bar V^\mu_{m \nu}\bar a_n a_\mu \bar a_\nu 
E[-\omega^\mu_{n \nu l},-\omega^m_{l n}]  \delta^\mu_{m + \nu}
\right) \delta^m_{l+ n} \right. \CR && \left.
+ 2 
\bar V^m_{ln} 
 \left(
\bar V^n_{\mu \nu}a_m \bar a_\mu  \bar a_\nu \delta^n_{\mu + \nu}
E[-\omega^m_{l\nu\mu},-\omega^m_{ln}] 
+
2  V^\mu_{n \nu}a_m \bar a_\mu  a_\nu E[\omega^{\mu l}_{\nu m},
-\omega^m_{ln}]\delta^\mu_{n + \nu}\right)\delta^m_{l+n}
\right],\CR\label{SecondIterate}
\EEA
\rightline{---------------------------------------------------------------------
----}
%
\begin{multicols}{2}
\noindent
where we used $a^{(2)}_k (0)=0$ and introduced 
%
$E(x,y)=\int_0^T \Delta(x-y)e^{i y t} d t .$
%

\section{ Statistical description --- } 

Let us now develop a statistical
description applying RPA to the fields $a^{(0)}_k$. Since phases and
the amplitudes are statistically independent in RPA, we will first
average over the random phases (denoted as $\langle
... \rangle_{\phi}$) and then we average over amplitudes (denoted as
$\langle ... \rangle_{A}$) to calculate the moments,
%
$$ M^{(p)}_k(T)\equiv \langle |a_k(T)|^{2p}\rangle_{\phi,A}.\ \ \ 
p=1,\ 2, \, 3\ ...,  $$
%
First, let us calculate $|a_l(T)|^{2 p}$ as
%
\end{multicols}
\leftline{-------------------------------------------------------------------------}
%
\BEA|a_l(T)|^{2 p}
= \left(a_l^{(0)}+\epsilon a_l^{(1)}+\epsilon^2 a_l^{(2)}\right)^p
\left(\bar a_l^{(0)}+\epsilon \bar a_l^{(1)}+\epsilon^2
\bar a_l^{(2)}\right)^p= 
|a_l^{(0)}|^{2p} +
\epsilon p |a^{(0)}_l|^{2p-2}\left(a_l^{(0)}\bar a_l^{(1)}+  \bar a_l^{(0)}
a_l^{(1)} \right) +\CR
\epsilon^2|a_l|^{2p-4}\Big[ C^2_p(a_l^{(0)} \bar a_l^{(1)})^2 +
C^2_p(\bar a_l^{(0)} a_l^{(1)})^2 + p^2 |a_l^{(0)}|^2 |a_l^{(1)}|^2 +
p | a_l^{(0)} |^2 \left( a_l^{(0)} \bar a_l^{(2)}+
\bar a_l^{(0)} a_l^{(2)}\right)\Big] + ... ,\nonumber \EEA 
\rightline{---------------------------------------------------------------------
----}
%
\begin{multicols}{2}
\noindent
where $C^2_p$ is the binomial coefficient.


Up to the second power in $\epsilon$ terms, we have
\BEA \langle |a_l(T)|^{2p}\rangle_\phi= 
 |a_l|^{2p} + \CR
\epsilon^2 |a_l|^{2p-2} \left( p^2 \langle |a_l^{(1)}|^2 \rangle_\phi + p \, \langle a_l^{(0)} \bar
a_l^{(2)}+ \bar a_l^{(0)} a_l^{(2)}\rangle_\phi \right)\CR
\label{pmoment}\EEA
Here, the terms
proportional to $\epsilon$ dropped out after the phase averaging.
Further, we assume that there is no coupling to the $k=0$ mode, i.e.
$V^{k=0}_{k_1 k_2} = V^{k1}_{k_1 k=0}=0$, so that there is no 
contribution of the
term like $ \langle (a_l^{(0)} \bar a_l^{(1)})^2 \rangle_\phi $.
We now use (\ref{FirstIterate}) and (\ref{SecondIterate})
and the  averaging over the phases to obtain
%
\BEA
 \langle |a^{(1)}_l|^2\rangle_\phi = 4 \sum_{m,n}^\infty \big[
|V^l_{mn}|^2 \delta^l_{m+n} |\Delta(\omega^l_{mn})|^2 
|a_m|^2 |a_n|^2  \nonumber \\ 
+2 |V^n_{lm}|^2 |\Delta^n_{l+m}|^2 \delta(\omega^n_{lm}) 
 |a_n|^2  |a_m|^2 
\big],  \CR 
 \langle a^{(0)}_l \bar a^{(2)}_l+\bar a^{(0)}_l a^{(2)}_l
\rangle_\phi =
-8{\bf |a_l|^2}  \sum_{m,n}^\infty \big[
|V^l_{mn}|^2 \delta^l_{m+n} E(0,\omega^l_{mn})  |a_l|^2
\CR
+|V^n_{lm}|^2 \delta^n_{l+m} E(0,\omega^n_{lm}) (|a_m|^2- |a_n|^2) 
\big].  \nonumber \EEA
%
Let us substitute these expressions into (\ref{pmoment}), perform
amplitude averaging, take the large box limit\footnote{The large box
limit implies that sums will be replaced with integrals, the Kronecker
deltas will be replaced with Dirac's deltas, $\delta^l_{m+n}\to\delta^l_{mn}/{\cal V}$, 
where we introduced short-hand notation,
$\delta^l_{mn}=\delta(k_l-k_m-k_n)$. Further we redefine
$M^{(p)}_k/{\cal V}^p \to M^{(p)}_k$.} and then large $T$ limit ($T \gg 1/ \omega$)\footnote{Note
that $\lim\limits_{T\to\infty}E(0,x)= (\pi
\delta(x)+iP(\frac{1}{x}))$, and
$\lim\limits_{T\to\infty}|\Delta(x)|^2=2\pi T\delta(x)$ (see e.g.
\cite{Ben}).}. We have 
%
\BE
 M^{(p)}_k(T) =  M^{(p)}_k(0) + T
\left(-p \gamma_k M^{(p)}_k + p^2\rho_k M^{(p-1)}_k\right),
\EE
%
with
%
\BEA \rho_k = 4\ \epsilon^2 \int d {\bf k_1} d {\bf k_2} 
( 
|V^k_{12}|^2 \delta^k_{12} \delta(\omega^k_{12}) n_{1} n_{2} 
 \nonumber\\  \left.
+2 |V^2_{k1}|^2 \delta^2_{k1} \delta(\omega^2_{k1}) n_{2} n_{1}
\right),  \label{RHO} \\
\gamma_k = 
8 \epsilon^2 \int  d {\bf k_1} d {\bf k_2}   
(
|V^k_{12}|^2 \delta^k_{12} \delta(\omega^k_{12}) n_{2} 
\CR
 +|V^2_{k1}|^2 \delta^2_{k1} \delta(\omega^2_{k1}) (n_{1}- n_{2}) 
).  \label{GAMMA}
\EEA
%
Now, assuming that $T$ is a lot less than the nonlinear time ($T \ll
1/\omega \epsilon^2$) we finally arrive at our main result, 
%
\BE 
\dot M^{(p)}_k = -p \gamma_k M^{(p)}_k + 
p^2 \rho_k M^{(p-1)}_k.\label{MainResultOne}\EE
%
In particular, for the waveaction spectrum $M^{(1)}_k=n_k $
(\ref{MainResultOne}) gives the familiar kinetic equation (KE)
%
\BE \dot n_{k} = -\gamma_k n_{k} +\rho_k=\epsilon^2 J(n_{k}),
\label{KE1}\EE
%
 where $ J(n_{k})$ is the ``collision'' term \cite{ZLF,GS}, 
%
$$ J(n_{k})=\int d k_2 d k_1 (R^k_{12}-R^1_{k2}-R^2_{1k}), $$
%
 with 
%
\BE R_{k12}=4\pi|V^k_{12}|^2 \delta^k_{12}\delta(\omega^k_{12})
\Big(n_{2}n_{1}-n_{k}(n_{2}+n_{1})\Big). \EE
%
The second equation in the series (\ref{MainResultOne}) allows to
obtain the r.m.s. $\sigma_k^2 = M^{(2)}_k - n^2_k$ of the fluctuations
of the waveaction $\langle |a_k|^2\rangle$.  We emphasize that
(\ref{MainResultOne}) is valid even for strongly intermittent fields
with big fluctuations.

Let us now consider the stationary solution of (\ref{MainResultOne}),
$ \dot M^{(p)}_k =0$ for all $p$. Then for $p=1$ from (\ref{KE1}) we
have $\rho_k=\gamma n_{k}$. Substituting this into
(\ref{MainResultOne}) we have 
%
$$ M^{(p)}_k = p M^{(p-1)}_k n_{k}\label{ResultTwo}, $$
%
with the solution
%
$ M_k^{(p)}=p!\ n_{k}^p$.
%
Such a set of moments correspond to a Gaussian wavefield $a_k$.  To
see how such a Gaussian steady state forms in time, let us rewrite
(\ref{MainResultOne}) in terms of the deviations of $M^{(p)}_k$ from
their Gaussian values,
%
$Q^{(p)}_k = M^{(p)}_k - p! \ n^p_k,\ p=1,2,... .$
%
Then use (\ref{MainResultOne}) to obtain \BE \dot Q^{(p)}_k = -
\gamma_k p Q^{(p)}_k + p^2 \rho_k Q^{(p-1)}_k,
\label{MainResultTwo}
\EE 
for $p=2,3,4,...$  This results has a particularly simple form for $p=2$,
because  $Q^{(1)}_k\equiv 0$ and we get a decoupled equation 
%
$$Q^{(2)}_k = - 2\gamma_k Q^{(2)}_k.$$  According to this equation the
deviations from Gaussianity can grow or decay depending on the sign of
$\gamma$. However, according to KE near the steady state $\gamma =
\rho/n >0$ (because $\rho$ is always positive), the deviations from
Gaussianity decay. A similar picture arises also for the higher
moments. The easiest way to see this is to choose initial conditions
where $n$ is already at its steady value (but not the higher
moments). Then (\ref{MainResultTwo}) becomes a linear system which can
be immediately solved as an eigenvalue problem.  For large time, the
largest of these eigenvalues, $\lambda = 2 \gamma_k$, will dominate
and the solutions tend to $Q^{(p)}_k = C^{(p)}_k \exp(-2 \gamma_k t)$
where $C^{(p)}_k$ satisfy a recursive relation $C^{(p)}_k = p^2 n_k
C^{(p-1)}_k/(p-2)$ and $C^{(2)}_k$ is arbitrary (determined by the
initial conditions). Thus we conclude that the steady state
corresponding to a Gaussian wavefield is stable.

Predictions of equation (\ref{MainResultOne}) about the behavior of
fluctuations of the waveaction spectra can be tested by modern
experimental techniques which allow to produce surface water waves
with random phases and a prescribed shape of the amplitude $|a_k|$
\cite{lev}.  It is even easier to test (\ref{MainResultOne})
numerically.  Consider for example capillary waves on deep water.  If
a Gaussian forcing at low $k$ values is present, the steady state
solution of the kinetic equation corresponds to the Zakharov-Filonenko
(ZF) spectrum of Kolmogorov type~\cite{ZLF,Zakfil}. It is given by 
%
\BE n_k = A k^{-17/4} \label{ZF},\EE
%
with $A=\sqrt{P} \rho^{3/2} C/\sigma^{1/4}$, where $P$
is the value of flux of energy toward high wavenumbers, $\rho$ and
$\sigma$ are the density and surface tension of water, and $C\simeq
13.98$. The simplest experiment would be to start with a
zero-fluctuation (deterministic) spectrum and to compare the
fluctuation growth with the predictions of (\ref{MainResultOne}).
Note that such no-fluctuations initial conditions were used in
\cite{OsbornePRL,ZakharovPRL}.  

Let us calculate the rate at which a fluctuations grow for such an
initial conditions. To do that let us assume that the spectrum $n_k$
is isotropic, that is it depends only on the modulus of the vector,
not on its directions. We then can make an angular averaging of
(\ref{GAMMA}) obtaining: 
%
\BEA
\gamma_k = 
8 \epsilon^2 \int  d k_1 d k_2 \Delta_{k k_1 k_2}^{-1}  
(
|V^k_{12}|^2 \delta^k_{12} \delta(\omega^k_{12}) n_{2} 
\CR
 +|V^2_{k1}|^2 \delta^2_{k1} \delta(\omega^2_{k1}) (n_{1}- n_{2}) 
).,
\nonumber\\ 
\Delta_{k 1 2} = \left< \delta({\bf k}-{\bf
k_1}-{\bf k_2})\right>\equiv \int \delta({\bf
k}-{\bf k_1}-{\bf k_2}) \, d \theta_1 d \theta_2 \, , \nonumber\\
\Delta _{k 1 2} = \frac{1}{2}\sqrt{
2 \left( (k k_1)^2 +(k k_2)^2 +(k_1 k_2)^2 
\right)-k^4-k_1^4 -k_2^4} \, . \nonumber
 \label{GAMMAdeltaAVERAGE}
\EEA


Let us substitute ZF spectrum (\ref{ZF}) into
(\ref{GAMMAdeltaAVERAGE}), take the values of $\omega_k$ and
$V^k_{12}$ appropriate for the capillary waves on deep
water(\cite{ZLF}, eqs (5.2.1-2)). By changing the variables of
integrations via $k_1= k \xi_1, \ \ k_2 = k \xi_2$ we can factor out
the $k$ dependence of $\gamma_k$. Performing one of $\xi$ integrals
analytically with the use of the delta function in $\omega$'s, we
perform the remaining single integral numerically to obtain (all the
integrals converge):
%
$${\gamma= \frac{4.30 A\sqrt{\sigma}}{16 \pi \rho^{3/2}} 
k^{3/4}},
$$
%
where the dimensionless constant $4.30$ was obtained by numerical
integration. Consequently, our prediction for the fluctuations growth
is 
%
\BEA
Q^{(2)}_k = Q^{(2)}_{k0} e^{-2 \gamma_k t}, \nonumber \\ Q^{(3)}_k = 9
Q^{(2)}_{k0} n_k e^{-2 \gamma_k t},
\EEA
%
etc.  Note that fluctuations stabilize at Gaussian values faster for
high $k$ values.  It is also interesting to test equation
(\ref{MainResultOne}) when the forcing (and therefore the turbulence)
is non-Gaussian, as in most practical situations. 

\section{ Discussion --- } 

In this manuscript, we derived a hierarchy of
equations~(\ref{MainResultOne}) for the one-point moments $M^{(p)}_k$
of the waveaction $|a_k|^2$. This system of equations has a
``triangular'' structure: the time derivative of the $p$-th moment
depends only on the moments of order $p, p-1$ and 1 (spectrum). Their
evolution is not ``slaved'' to the spectrum or any other low moments
and it depends on the initial conditions. RPA allows the initial
conditions to be far from Gaussian and deviation of n'th moment from
its Gaussian level may even increase in a transient time dependent
state. Among two allowed extreme limits are the wavefield with a
deterministic amplitude $|a_k|$ (for which $M^{(p)}_k = n^p_k$) and
the intermittent wavefields characterized by sparse k-space
distributions of $|a_k|$ (for which $M^{(p)}_k \gg n^p_k$).  However,
the level of intermittency (and non-Gaussianity in general)
exponentially decreases as WT approaches a statistically steady state,
as given by (\ref{MainResultTwo}).  Importantly, in some situations WT
might never reach a steady state (for example because of the
nonstationary pumping) or it might spend a long time in a transient
non-stationary state. In this case WT can be highly intermittent and
yet equations (\ref{MainResultOne}) are still valid for description of
such turbulence. In the other words, the type of intermittency
discussed in the present manuscript appears within the weakly nonlinear
closure and not as a result of its breakdown as in
\cite{NazarenkoNewell}). 

{\bf Acknowledgments} Authors thank
anonymous referees for constructive comments and Alan Newell for
enlightening discussions. YL is supported by NSF CAREER grant DMS
0134955 and by ONR YIP grant N000140210528. SN thanks ONR for the
support of his visit to RPI.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%% THE BIBLIOGRAPHY %%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{thebibliography}{190}
\bibitem{ZLF}V.E.  Zakharov, V.S.  L'vov 
and G.Falkovich, "Kolmogorov Spectra of Turbulence", Springer-Verlag, 1992.
\bibitem{Ben} D.J. Benney and P.Saffman, Proc Royal. Soc, A(1966),
289, 301-320; B.  J . Benney and A.C.  Newell, Studies in Appl.
Math.  {\bf 48} (1) 29 (1969).
\bibitem{GS} A.A. Galeev and R.Z. Sagdeev, in Reviews of Plasma Physics 
Vol. 6 (Ed. M A Leontovich) (New York: Consultants Bureau, 1973) 
\bibitem{Zakfil} V.E.Zakharov and Filonenko, J. Appl. Mech. Tech.
Phys. {\bf 4} 506-515 (1967).
\bibitem{hasselman} K. Hasselmann, J. Fluid Mech {\bf 12} 481 (1962).
% Freely decaying weak turbulence for sea surface gravity waves 
\bibitem{OsbornePRL}
M. Onorato {\it et al},   
PRL {\bf 89}, 144501, (2002).
\bibitem{ZakharovPRL}
%Turbulence of capillary waves 
A.N. Pushkarev, V.E. Zakharov, PRL {\bf 76}, 3320-3, (1996). Physica
D, {\bf 135}, 98, (2000).
\bibitem{SoomerePRL}
%Generation of zonal flow and meridional anisotropy in two-layer weak geostrophic turbulence 
T. Soomere, PRL {\bf 75}, 2440-3, (1995).
\bibitem{LT}
Yu. V. Lvov and E.G.  Tabak, 
PRL {\bf 87}, 168501, (2001) 
\bibitem{NazarenkoNewell} L. Biven, S.V. Nazarenko and A.C. Newell, 
%''Breakdown of wave turbulence and the onset of intermittency'', 
Phys Lett A, {\bf 280},  28-32, (2001)
A.C. Newell, S.V. Nazarenko and L. Biven, Physica
D, {\bf 152-153}, 520-550, (2001).
\bibitem{lev} L. Shemer, H. Jiao, E. Kit and Y. Agnon 
%Evolution of a nonlinear wave field along a tank: experiments and numerical simulations based on the spatial Zakharov equation. 
J. Fluid Mech. {\bf 427}, 107-129, (2001)
\end{thebibliography}
\end{multicols}
\newpage

Letter to the editor:


Dear Editor, 

we would like to resubmit our paper for PRE. As explained below, the
main objectives of the referee's are coming from poor understanding of
the issues involved in our paper. We ahve sligtly expanded the
manuscript, and we hope that it could be published in PRE with out
further signinficant delays. (kak skazat; eto pomagche?).

The main point of our paper - we present NEW results (not technical
improvement of the old results), as we give a QUANTITATIVE and
VERIFIABLE predictions for the quantity, which has been traditionally
ignored in the theory of weakly interacting dispersive waves. 

We answer referee's questions below:


We disagree with the referee A that there are mistakes
in our work, - as you will see below this is a result 
of not too careful reading. We also disagree that
our work is just a technical improvement and not an
essentially new theory with respect to the classical
weak turbulence. Indeed, the classical weak turbulence
deals with short correlated wave fields close to gaussian
whereas we consider the long correlated fields the deviation
from gaussianity for which is of order one. These fields
are widespread in numerical simulations, including the
operational wave weather forecasts which always start
with fields with random phases and deterministic
amplitudes.

\small {`` ... Indeed, the signature in Fourier space of long-range
correlations in physical turbulence is generally a power law
singularity for small values of wavenumber, rather than a
delta-function singularity. If there is some reasonable model of
turbulence which does have delta-function singularities in the
cumulant spectrum, it should be cited in the manuscript; I am
unaware of any such example.''}

Examples of such wave fields are widespread in numerical
simulations. Indeed, the typical start is a field with
random phases and deterministic amplitudes, which is
long correlated in the x-space, strongly non-gaussian
and is a special case of the fileds we are dealing with.
Note that the referee A wrongly calls such field deterministic
and this is the main problem with his objections (see
also below).


{\small{``
... So if one stays instead within the context of random
models of turbulence, what kind of experiments and numerical
simulations are contemplated? Is weak turbulence non-Gaussian? I
know strong turbulence has non-Gaussian features, but there should
be some target experiment or numerical simulation cited which
motivates the theory.''}

The referee has not read our paper carefully. We are
not talking about the deterministic initial data.
Instead, we talk about deterministic absolute values
and  random phases. This is a typical statistical
initial condition for numerics and it is obviously 
a special case of the fields we are dealing with
(RPA requirements are described on the first page).
This mix-up seems to be the main basis for the
referee A objections.

{\small{``The discussion about the Fourier spectrum being pinned to
the Gaussian value in the large-box limit is misleading.  Indeed, the
Fourier-space statistics are nearly Gaussian, but the small deviation
(proportional to some inverse power of the volume of the box) plays a
physically important role reflecting the deviations from Gaussianity
in physical space, and the third order cumulant plays an important
role in giving the correct kinetic equation even for the
nearly-Gaussian standard weak turbulence theory!}}


Yes, in standard weak turbulence there are small 
but important deviations from gaussianity. However,
our analysis includes fields which are order-one 
different from being gaussian and this is what we 
wanted to emphasize.

{\small{``
The claim (on p. 4) that the
one-point moments contain information wiped out in the multiple-point
moments is misleading. Indeed, the one-point moments are just the
multiple-point moments evaluated at coalescing values of wavenumber.
There is no ``fusion'' problem in computing the multiple-point moments
and then evaluating them with coalescing arguments in wavenumber space;
the ``fusion'' problem has to do with coalescing points in phsyical
space within the context of a high Reynolds number theory (which is not
weak turbulence!).  The source of confusion is that the cumulants of
non-coincident values of the wavenumber are simply assumed to be zero by
an ad hoc invocation of RPA; so the loss of information is due to the
RPA assumption and not any physical mechanism.''
}}

This is indeed a classical "fusion" effect, although
in the Fourier rather than the coordinate space.
It also has classical roots of non-cummuting limits,
in our case, these are the limits of large time and
of zero phase correlation length (in the k-space).
And of course there are deep physical reasons for this
phase decorrelation to appear in the dispersive wave
systems.

{\small{``
On p. 8, the time scale $ T $ of interest is assumed to be much
less than the nonlinear time $ 1/\omega \epsilon^2 $. But the
differential equations for the cumulants in (8) have the right
hand sides proportional to $ \varepsilon^2 $, so with this
restriction, not much can be said because the time scale is
restricted so that the $ Q_k^{(p)} $ change negligibly from their
initial values. Now, one might object that when $ \omega $ is
small, there are nontrivial predictions, but often the coupling
coefficients for nonlinear waves vanish along with $ \omega $ so
if $ \omega $ is small, so will $ \gamma_k $ and $ \rho_k $ be
small and it still seems like the restriction of the prediction in
equation (8) to time scales $ T \ll 1/\omega \varepsilon^2$ means
simply that not much happens to the cumulants on these time
scales, which is not surprising since one would expect the
cumulants to evolve on the actual nonlinear time scale.''}}


Typical timescale in our paper is $1/epsilon^2$ and not $T$.
$T$ is just a technical intermediate time needed for
the time-scale separation technique (same as in the
derivation of the usual Kinetic equation). 
Our final equations show that all qumulants evolve at
the same $1/epsilon^2$ timescale as the spectrum itself
and they can experience order 1 changes.

In summary, we 

\end{document}



       
