%Last edit YL 12/15/2004

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

\def \BE {\begin{equation}}
\def \EE {\end{equation}}
\def \BEA {\begin{eqnarray}}
\def \EEA {\end{eqnarray}}
\def \CR {\nonumber \\}
\def \zer {^{(0)}}
\def \one {^{(1)}}
\def \two {^{(2)}}
\def \kp {{k^\prime}}
\def \e {{\epsilon}}

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

\begin{document}
\title{Anomalous probability of large amplitudes in wave turbulence\\
{\bf Physics Letters A,Volume 339, Pages 361-369, 2005}}

\author{Yeontaek Choi$^*$, Yuri V. Lvov$^\dagger$, Sergey Nazarenko$^*$
and Boris Pokorni$^\dagger$}
\address
{
$^*$
Mathematics Institute, The University of Warwick,  Coventry, CV4-7AL, UK \\
$^\dagger$
Department of Mathematical Sciences, Rensselaer Polytechnic Institute,
Troy, NY 12180 }

\maketitle



\begin{abstract}
{Time evolution equation for the Probability Distribution Function
(PDF) is derived for system of weakly interacting waves, dominated by
the four-wave process. It is shown that a steady state for such system
may correspond to strong intermittency.}
\end{abstract}


\begin{multicols}{2}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
{\bf Introduction ---} Wave Turbulence (WT) is a common name for the
fields of dispersive waves which are engaged in stochastic weakly
nonlinear interactions over a wide range of scales.  Numerous examples
of WT are found in oceans, atmospheres, plasmas and Bose-Einstein
condensates~\cite{ZLF,Ben,GS,Zakfil,hasselman,zaslav,OsbornePRL,JaLubimii,DNPZ}.
For a long time, describing and predicting the energy spectra was the
only concern in WT theory. More recently, some attention was given to
the study of turbulence intermittency.  WT intermittency, or
``burstiness'' of the turbulent signal, was observed experimentally
and numerically and was attributed, as in most turbulent systems, to
the presence of coherent structures. Examples include collapsing
filaments in Bose-Einstein condensates with attractive potentials
\cite{DNPZ,Bradley}, condensate quasi-solitons in systems with
repulsive potentials \cite{DNPZ,Pitaevsky,zgpd}, white caps of sea
waves at small scales \cite{rough}, freak ocean waves at larger scales
\cite{jansen}.  Often, such coherent structures are intense but quite
sparse so that in most of the space waves remain weakly nonlinear and
mostly unaffected by these structures.

Recent analysis of the higher order cumulants \cite{biven} showed that
WT becomes strongly non-Gaussian at the same length scale where it
fails to be weakly nonlinear. In scale invariant systems, the ratio of
nonlinear time to the linear wave period grows as a power-law either
in to small or toward large wavenumbers. When this growth coincides
with the cascade direction then one expects the WT breakdown if the
inertial range is large enough. Otherwise intermittency never occurs
provided that turbulence is weak at the forcing scale \cite{cnn}.
Further, even if a significant non-Gaussianity occurs, it does not in
itself imply intermittency because PDF may remain, in principle, of
the same order as Gaussian in all of its parts. This motivates us to
study PDFs in WT. Study of PDF in WT context can be traced back to as
early as the work of Pierls~\cite{peierls}, and, latter,
in~\cite{bp,prig}, who considered waves in anharmonic crystals, a
special case of 3-wave systems. Recently, equations for multimode and
one-mode PDF's where derived and analyzed for the general case of
3-wave systems~\cite{pla,physd}. PDF's of the three wave systems were
also studied to explain entropy production in three-wave turbulence
systems~\cite{Newell}. In the present paper, we will be concerned with
the 4-wave case.  We are also motivated by a puzzling numerical
evidence of a low-wavenumber intermittency in the system of
water-surface gravity waves \cite{yoko} whereas the analysis of
\cite{biven} predicts intermittency at high wave numbers only.
Explaining this fact could shed light on the phenomenon of freak waves
\cite{jansen}.

The idea of the present letter is based on the observation that even
if the ``hard'' breakdown (as in \cite{biven}) does not occur, there
will always be a part of the PDF tail for which the amplitudes are too
high for WT to work. Such a ``mild'' breakdown will modify the PDF
tail in a way that may correspond to intermittency.  In fact, this
case is easier to study analytically because WT still works for most
of the PDF and the wave breaking phenomenon can be modeled simply as a
phenomenological cutoff of the PDF tail reflecting the fact that no
waves exist above the breaking amplitude. The wave breaking causes
``leakage'' and, therefore, a flux in the {\em amplitude} space which
is the key phenomenon leading to deviations from the Gaussian
equilibrium and intermittency. Note an analogy with the well-known
{k-space} fluxes (cascades) corresponding to Kolmogorov turbulence
which is qualitatively different from the thermodynamic equilibrium
state.  In this paper we will derive an equation for the wave
amplitude PDF and we will find its steady state solutions
corresponding to the finite flux in the {\em amplitude} space.
Consequently, we will show that the resulting wave fields are
intermittent at each wavenumber with an anomalously large probability
of the large-amplitude waves.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
{\bf Definition of RPA fields --- }
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  Previously, the random phase approximation (RPA) has typically
assumed that the phases evolve much more rapidly than the amplitudes
and, therefore, there exist time intervals where the phases are random
but the amplitudes are deterministic \cite{ZLF}. However, numerical
simulations indicate that the phase and the amplitude vary at the same
time scale \cite{push}.  Thus, we need to generalize RPA to the case
where both the phases and the amplitudes are random quantities.  Such
generalization was done in \cite{physd,pla,ln} where 3-wave systems
were considered. In the present letter, we will be dealing with 4-wave
systems.

Let us consider a wavefield $a({\bf x}, t)$ in a periodic box of
volume ${\cal V}$ and let the Fourier transform of this field be $a_k$
where $ k \in {\cal Z}^d$ and $d$ is the space dimension.  Later we
take the large box limit in order to consider homogeneous wave
turbulence. Let us write complex $a_k$ as $a_k =A_k \psi_k $ where
$A_k \in {\cal R}^+ $ is the amplitude and $\psi_k \in {\cal S}^1$ is
a phase factor (${\cal S}^1$ being the unit circle in the complex
plane).  We say the wavefield $a_k$ is of the {\em RPA type} if all
variables in the set $\{ A_k, \psi_k ; k \in {\cal Z}^d \}$ are
statistically independent random variables and $\psi_k$'s are
uniformly distributed on ${\cal S}^1$.
%
Defined this way RPA refers not only to the phase but also the
amplitude statistics and therefore we suggest a slightly different
reading of this acronym: ``Random Phase and Amplitude''.

The above properties are sufficient for our WT analysis and yet such
fields may be strongly non-Gaussian.  Indeed, RPA allows any shape of
the PDF for amplitudes $A_{k}$ and, therefore, it will be a good tool
for describing intermittency.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
{\bf  Weakly nonlinear evolution ---}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Consider a weakly nonlinear wavefield dominated by the 4-wave
interactions, e.g.  the water-surface gravity
waves~\cite{ZLF,hasselman,OsbornePRL,rough}, Langmuir waves in
plasmas~\cite{ZLF,GS} and the waves described by the nonlinear
Schroedinger equation~\cite{DNPZ}.  In the finite box, we have the
following Hamiltonian equations for the Fourier modes of this field,
%
\begin{equation}
i \dot b_l = \epsilon \sum_{\alpha\mu\nu}
{ W^{l\alpha}_{\mu\nu}}\bar b_\alpha b_\mu b_\nu 
e^{i\omega^{l\alpha}_{\mu\nu}t}
\delta^{l\alpha}_{\mu\nu}
\label{FourWaveEquationOfMotionB}
\end{equation}
%
where $b_l$ is the wave action variable in the interaction
 representation, $l \in {\cal Z}^d$, $W^{l \alpha} _{\mu \nu} \sim 1$
 is an interaction coefficient,
%
$\omega^{l\alpha}_{\mu\nu} =
\omega_l+\omega_{\alpha}-\omega_{\mu}-\omega_{\nu}$,
%
$\omega_l$ is the frequency of mode $l$
%
and $\epsilon \ll 1$ is a nonlinearity parameter.  We are going expand
in $\epsilon$ and consider the long-time behavior of a wave field, but
in order to make such an analysis consistent we have to renormalize
the frequency of (\ref{FourWaveEquationOfMotionB}) as
%
%
\begin{equation}
i \dot a_l = \epsilon \sum_{\alpha\mu\nu}
{ W^{l\alpha}_{\mu\nu}}\bar a_\alpha a_\mu a_\nu 
e^{i\tilde\omega^{l\alpha}_{\mu\nu}t}
\delta^{l\alpha}_{\mu\nu} - \Omega_l a_l,
\label{FourWaveEquationOfMotionA}
\end{equation}
%
where $a_l = b_l e^{i \Omega_l t} $, $\Omega_l = 2 \epsilon \sum_\mu
 W^{l \mu}_{l \mu} |A_\mu(0)|^2 $ is the nonlinear frequency shift
 arising from self-interactions and 
%
$\tilde \omega^{l\alpha}_{\mu\nu} = \omega^{l\alpha}_{\mu\nu}
 +\Omega_l+\Omega_\alpha-\Omega_\mu-\Omega_\nu$.
%

For small nonlinearity, the linear time-scale $2 \pi / \omega$ is a
lot less than the nonlinear evolution time which (as will be evident
below, see e.g. (\ref{pdeofmoment})) is $2 \pi \epsilon^2/ \omega$.
Thus, to filter out fast oscillations at the wave period, let us seek
for the solution at an intermediate time $T$ such that $2 \pi / \omega
\ll T \ll 1/(\omega \epsilon^2)$.  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)}.
%\label{Expansion}
$
%
Substituting this in (\ref{FourWaveEquationOfMotionA}) we get in the
zeroth order a time independent result,
%
$ a_l^{(0)}(T)=a_l(0)\label{definitionofa} $.
%
For simplicity, we will write $a_l(0)= a_l$, understanding that a
quantity is taken at $T=0$ if its time argument is not mentioned
explicitly.  The first iteration of (\ref{FourWaveEquationOfMotionA})
gives
%
\begin{eqnarray}
a_l\one(T) = - i \sum_{\alpha\mu\nu} W^{l\alpha}_{\mu\nu} \bar
 a_\alpha a_\mu a_\nu \delta^{l\alpha}_{\mu\nu} \Delta^{l
 \alpha}_{\mu\nu} + i {\Omega_l \over \epsilon} a_l T.  \label{FirstIterate}
\end{eqnarray}
%
where 
%
$$\Delta^{l \alpha}_{\mu\nu} \equiv \Delta^{l \alpha}_{\mu\nu}(T) =
({e^{i\tilde\omega^{l\alpha}_{\mu\nu}T}-1})/{i \tilde\omega^{l\alpha}_{\mu\nu}}.
 \label{NewellsDelta}$$
%
Iterating one more time we get
\end{multicols}
\leftline{---------------------------------------------------------------------
----}
\begin{eqnarray}
a_l\two(T)=\sum_{\alpha\beta\mu\nu v u}\left( 
W^{\mu\nu}_{\alpha u} W^{lu}_{v\beta} \delta^{\mu\nu}_{\alpha u} \delta^{lu}_{v\beta}
 a_{\alpha} a_v a_\beta 
\bar a_{\mu} \bar a_{\nu} 
E(\tilde\omega^{l\mu\nu}_{\alpha v \beta}, \tilde\omega^{lu}_{v\beta}) 
- 
2 W^{\alpha v}_{\mu\nu} W^{l u}_{v \beta}
\delta^{\alpha v}_{\mu\nu} \delta^{l u}_{v \beta}
\bar a_{\alpha} \bar a_{u}
a_{\mu} a_{\nu} a_\beta 
E(\tilde \omega^{l\alpha u}_{\mu\nu\beta}, \tilde \omega^{lu}_{v\beta})\right) 
-\Omega_l^2 a_l \frac{T^2}{2\epsilon^2}
\nonumber  \\ 
+ \frac{1}{\epsilon}
\sum_{\alpha\mu\nu}\left(
\Omega_l 
W^{l\alpha}_{\mu\nu}
\delta^{l\alpha}_{\mu\nu}
\bar a_{\alpha} a_{\mu} a_{\nu} 
E(\tilde\omega^{l\alpha}_{\mu\nu},0)  
-W^{l \alpha}_{\mu\nu }
\delta^{l \alpha}_{\mu\nu }\bar{ a_{\alpha}} 
a_{\mu} a_{\nu}(\Omega_\alpha
-2\Omega_\nu) 
\int \limits_0^T
\tau e^{i\tilde\omega^{l \alpha}_{\mu\nu } \tau} d \tau
\right)
, \hbox{\ with\ }
E(x,y)=\int_0^T \Delta(x-y)e^{i y t} d t .
\label{SecondIterate}
\end{eqnarray}
%
\rightline{---------------------------------------------------------------------
----}
\begin{multicols}{2}

{\bf Evolution of statistics --- } We will now develop a statistical
description via averaging over the initial fields $a_k(0)$ which are
taken to be of the RPA type.  Of course, to have a non-trivial
description valid over the nonlinear evolution time, the fields must
remain of the RPA type over the nonlinear time in the leading order in
$\epsilon$. The proof of this statement in the 3-wave case was
presented in~\cite{pla,physd}. Also \cite{physd} contains an
announcement of the 4-wave equation for the multi-mode PDF; the
details of the relevant analysis will be published separately.  Let us
introduce a generating function $ Z(\lambda ,t)=\langle e^{\lambda
|a_k(t)|^2}\rangle, $
%
where $\lambda$ is a real parameter.  Then PDF of the wave intensities
$s= |a_k(t)|^2$ at each ${\bf k}$ can be written as an inverse Laplace
transform,
%
$
P(s,t) = 
\langle \delta (|a_k(t)|^2 -s) \rangle =
{1 \over 2 \pi i}  \int_{-i \infty}^{+i \infty} Z(\lambda, t)e^{-s
\lambda}d\lambda. \label{distribution}
$
%
For the one-point moments  we have
\begin{eqnarray} 
M_k^{(p)} \equiv \langle
|a_k|^{2p}\rangle 
=\langle
|a|^{2p}e^{\lambda |a|^2}\rangle|_{\lambda=0}=\nonumber\\Z_{\lambda \cdots \lambda}|_{\lambda=0}
= \int_{0}^{\infty} s^p P(s, t) \, ds \label{momPDF}, 
\end{eqnarray}
where $p \in {\cal N}$ and subscript $\lambda$ means differentiation
with respect to $\lambda$ $p$ times.

At $t=T$ we have 
%
\begin{eqnarray}
Z(T) = \langle e^{\lambda|a_k^{(0)} +\epsilon a_k^{(1)} +\epsilon^2
a_k^{(2)}|^2}\rangle  = \hspace{1.5cm} \nonumber \\
\langle e^{\lambda
|a_k^{(0)}|^2} \langle 1+\lambda \epsilon (a_k^{(1)} \bar a_k^{(0)}+
{\rm cc})
\nonumber \\
 \lambda
\epsilon^{2}(|a_k^{(1)}|^{2} + (a_k^{(2)}\bar a_k^{(0)} +{\rm cc}) 
) + 
\nonumber \\
{\lambda^2 \e^2 \over 2} (a_k^{(1)}\bar a_k^{(0)} +{\rm cc}
% \bar a_k^{(1)}a_k^{(0)}
)^2
\rangle_{\psi} \rangle_{A}
\nonumber \\
=Z(0) + \e \lambda \langle e^{\lambda
|a_k^{(0)}|^2} \langle a_k^{(1)} \bar a_k^{(0)}+
{\rm cc} \rangle_{\psi} \rangle_{A} + \nonumber \\
\epsilon^{2} \langle \langle (\lambda + \lambda^2 A^2) |a_k^{(1)}|^{2} 
+ \lambda (a_k^{(2)}\bar a_k^{(0)} +{\rm cc} ) \nonumber \\ 
+ {\lambda^2  \over 2} (a_k^{(1)2} \bar a_k^{(0)2} +{\rm cc})
\rangle_{\psi} \rangle_{A} \hspace{1.5cm}
\label{zmoment} 
\end{eqnarray}
%
where ${\rm cc}$ stands for complex conjugate of the previous terms
and $\langle\dots\rangle_\psi$ and $\langle\dots\rangle_A$ denote
phase and amplitude averaging respectively. Note that in RPA fields
the phases and the amplitudes are statistically independent so that
these two averaging could be done independently. First let us
substitute $a_k^{(1)}$ and $a_k^{(2)}$ from (\ref{FirstIterate}) and
(\ref{SecondIterate}) respectively and perform the phase
averaging. For the terms proportional to $\epsilon$ we have
%
\begin{eqnarray}
\langle a_k\one\bar a_k^{(0)}\rangle_\psi= \nonumber \\
 - i \sum_{\alpha\mu\nu} W^{k\alpha}_{\mu\nu} 
\langle \bar a_k \bar
 a_\alpha a_\mu a_\nu\rangle_\psi
 \delta^{k\alpha}_{\mu\nu} \Delta^{k
 \alpha}_{\mu\nu} + i \Omega_l \langle |a_k|^2 \rangle  T\nonumber\\
 =- 2 i \sum_{\alpha} W^{k\alpha}_{k\alpha} 
 A_k^2 A_\alpha^2 \cdot T
+ i \Omega_l A_k^2 T.
\label{EpsilonTerms}
\end{eqnarray}
%
where we have used the fact that $\Delta(0)=T$ and we have used the
RPA's ``Wick's Theorem''
%
$$\langle \bar a_k \bar  a_\alpha a_\mu a_\nu\rangle_\psi =
A_k^2 A_\alpha^2 (
\delta^k_\mu \delta^\alpha_\nu +\delta^k_\nu\delta^\alpha_\mu).
$$ 
%
%
Note that the above expression should also contain the {\it singular
cumulant}, i.e. term $(A_\alpha^4 - 2A_\alpha^2) \delta^\alpha_\nu
\delta^\alpha_\mu \delta^\alpha_k$, see \cite{ln}. However we do not
write this term here, since its contribution is of the order of $1/N$
smaller because it has one extra delta function.

We see from (\ref{EpsilonTerms}) that the choice 
%
\begin{equation}
\Omega_k  =2 \sum_{\alpha} W^{k\alpha}_{k\alpha} A_\alpha^2
\label{FrequencyRenormalization}
\end{equation}
%
makes the contribution of $\langle a_k\one\bar a_k^{(0)}\rangle_\psi$ terms 
to be equal to zero. 

We therefore obtain
%
\begin{eqnarray}
Z(T) -Z(0) = \hspace{2.5cm} \nonumber\\ 
\lambda
\epsilon^{2}
\Big\langle
\big\langle
\lambda |a_k^{(1)}|^{2} (1+\lambda |a_k|^2)+
a_k^{(2)}\bar a_k^{(0)} + 
\bar a_k^{(2)} a_k^{(0)}+\nonumber\\ 
+\frac{\lambda}{2} ( (\bar a_k^{(0)} a_k\one)^2+{\rm cc}) 
\big\rangle_{\psi} \Big\rangle_{A}. \label{Zintermediate}
\end{eqnarray}
%
where as an
intermediate result we have
%
\begin{eqnarray}
\langle \bar a_k^{(0)} a_k\two\rangle_\psi=\hspace{2.5cm}
\nonumber\\
2\sum\limits_{\alpha \mu\nu}
\delta^{k\alpha}_{\mu\nu}
|W^{k\alpha}_{\mu\nu}|^2 A_k^2 
( A_\mu^2 A_\nu^2 - 2 A_\alpha^2 A_\mu^2)
 E(0,\tilde\omega^{k\alpha}_{\mu\nu})\nonumber\end{eqnarray}
%
and
%
\begin{eqnarray}
\langle |\bar a_k\one|^2\rangle_\psi=
\sum\limits_{\alpha \mu\nu}\delta^{k\alpha}_{\mu\nu}
|W ^{k\alpha}_{\mu\nu}|^2 A_\alpha^2 A_\mu^2 A_\nu^2 
|\Delta(\tilde\omega^{k\alpha}_{\mu\nu})|^2.
\nonumber
\end{eqnarray}
%
Here terms proportional to $T^2$ drop from $\langle \bar a_k^{(0)}
a_k\two\rangle_\psi$ and $\langle |\bar a_k\one|^2\rangle$ because of
the choice (\ref{FrequencyRenormalization}) of frequency
renormalization.  
Furthermore, 
%
\begin{eqnarray}
\langle
(a_k\one\bar a_k^{(0)})^2                                                              
\rangle_\psi=\nonumber\\
-2 W^{k k }_{k k } A_k^6 T^2 \Omega_k
-A_k^4 \Omega_k^2 T^2 - 4 \sum_\alpha(W^{k\alpha}_{k\alpha}
A_k A_\alpha T )^2\nonumber
\end{eqnarray}
%

To complete the derivation of the equation for the time evolution of
the generating function $Z(T)$ we have to take a large box limit,
which implies that sums will be replaced with integrals, the Kronecker
deltas will be replaced with Dirac's deltas, $\delta^{l \alpha}_{m
n}\to\delta^{\alpha l}_{mn}/{\cal V}$, where we introduced short-hand
notation, $\delta^{\alpha l}_{mn}=\delta(k_\alpha+k_l-k_m-k_n)$.  Then
(\ref{Zintermediate}) will still hold, but with
%
\begin{eqnarray}
\langle \bar a_k^{(0)} a_k\two\rangle_\psi= \hspace{2.5cm}
\nonumber\\
2\int d 123
\delta^{k1}_{23}
|W^{k1}_{23}|^2 A_k^2 
( A_2^2 A_3^2 - 2 A_1^2 A_2^2)
 E(0,\tilde\omega^{k1}_{23})\nonumber\end{eqnarray}
%
and
%%
\begin{eqnarray}
\langle |\bar a_k\one|^2\rangle_\psi=
\int d 123\delta^{k1}_{23}
|W ^{k1}_{23}|^2 A_1 A_2 A_3 
|\Delta(\tilde\omega^{k1}_{23})|^2.
\nonumber
\end{eqnarray}
%
We also have
%
\begin{eqnarray}
\langle
(a_k\one\bar a_k^{(0)})^2
\rangle_\psi = 0, 
\end{eqnarray}
%
because this terms will be $1/N$ times smaller than $\langle |\bar
a_k\one|^2\rangle_\psi$ and $\langle \bar a_k a_k\two\rangle_\psi$
terms because it has one less summation index. Therefore it vanishes
in the $N\to\infty$ limit.


Further we take a large $T$ limit, and take into account that
%
$$\lim\limits_{T\to\infty}E(0,x)= (\pi
\delta(x)+iP(\frac{1}{x})) T,$$
%
and
$$\lim\limits_{T\to\infty}|\Delta(x)|^2=2\pi T\delta(x),$$
%
(see e.g. \cite{Ben}).

Finally we perform amplitude averaging, noticing that 
%
$$Z(0)=\langle e^{\lambda |a_k|^2}\rangle_A,$$
%
and
%
$$\frac{\partial Z}{\partial \lambda} = \langle |a_k|^2 e^{\lambda |a_k|^2}
\rangle_A.$$
%
to obtain
%
\begin{equation}
Z(T) = Z(0) + \epsilon^2T\cdot( 
\lambda \eta Z +(\lambda^2 \eta - \lambda \gamma)
 Z_\lambda).\label{pdeofmoment2}
\end{equation}
%
Approximating $(Z(T)-Z(0))/T$ by $\dot Z$, we have
%
\begin{equation}
\dot Z = \lambda \eta Z +(\lambda^2 \eta - \lambda n \gamma)
 Z_\lambda,\label{pdeofmoment}
\end{equation}
%
where
\begin{eqnarray}
\eta_k &=& 4 \pi \epsilon^2 \int
|W^{k1}_{23}|^2 \delta^{k1}_{23}
\delta(\omega^{k1}_{23}) n_1 n_2 n_3 \, d123,\nonumber
\\ 
\gamma_k &=&
8 \pi \epsilon^2 \int  |W^{k1}_{23}|^2 \delta^{k1}_{23}
\delta(\omega^{k3}_{12}) 
 \Big[ n_1 (n_2 + n_3) - n_2 n_3\Big] \, d123,\nonumber
\end{eqnarray}
here wavenumbers $k, k_1, k_2, k_3 \in {\cal R}^d$, $\delta$'s now
mean Dirac $\delta$-functions, $n_{1,2,3} \equiv n(k_{1,2,3}) $ and
$d123 = d {\bf k_1} d {\bf k_2} d {\bf k_3}$. Differentiating
(\ref{pdeofmoment}) with respect to $\lambda$ $p$ times we get the
evolution equation for the moments:
%
$$ \dot M^{(p)}_k = -p \gamma_k
M^{(p)}_k + p^2 \eta_k M^{(p-1)}_k,\label{MainResultOne} $$
%
which, for $p=1$ gives the standard kinetic equation, $ \dot n_k = -
\gamma_k n_k + \eta_k . $ First-order PDE (\ref{pdeofmoment}) can be
easily solved by the method of characteristics.  Its steady state
solution is $$ Z=(1 -\lambda n_k)^{-1} $$ which corresponds to the
Gaussian values of momenta $M^{(p)} = p! n_k^p$.  However, these
solutions are invalid at small $\lambda$ and high $p$'s because large
amplitudes $s=|a|^2$, for which nonlinearity is not weak, strongly
contribute in these cases. Due to the integral nature of definitions
of $M^{(p)}$ and $Z$ with respect to the $s=|a|^2$, the ranges of
amplitudes where WT is applicable are mixed with, and contaminated by,
the regions where WT fails. Thus, to clearly separate these regions it
is better to work with quantities which are local in $s=|a|^2$, in
particular the probability distribution $P(s)$.  Taking the inverse
Laplace transform of (\ref{pdeofmoment}) we have
\begin{equation}
\dot P+\partial_{s}F=0, \label{peqn}
\end{equation}
 where $ F=-s(\gamma P+\eta\partial_{s}P) $ is a probability flux in
the s-space.
%
Consider the steady state solutions, $\dot P =0$, 
%
\begin{equation}
-s(\gamma P+\eta\partial_{s}P) = F= \hbox{const}. \label{f} 
\end{equation}
%
Note that in the steady state $\gamma /\eta = n_k$ which follows from
kinetic equation.  The general solution to (\ref{f}) is
%
$$P=P_{hom} + P_{part}$$
%
 where 
%
$$P_{hom} = \hbox{const} \, \exp{(-s/n)}$$ 
%
is the general solution to the homogeneous equation (corresponding to
$F=0$) and $P_{part}$ is a particular solution, 
%
$$P_{part} = -({F}/{\eta}) Ei({s}/{n}) \exp{(-s/n)} $$
%
where $Ei(x)$ is the integral exponential function.

At the tail of the PDF, $s \gg n_k$, the solution can be represented
as series in $1/s$,
%
$
P_{part} = - F/( s\gamma)-\eta F/(\gamma s)^2+\cdots.
$
%
Thus, the leading order asymptotic of the finite-flux solution is
$1/s$ which describes strong intermittency.

Note that if the weakly nonlinearity assumption was valid uniformly to
$s=\infty$ then we would have to put $F=0$ to ensure positivity of $P$
and the convergence of its normalization, $\int P \, ds =1$. In this
case $P=P_{hom} = n \, \exp{(-s/n)}$ which is a pure Rayleigh
distribution corresponding to the Gaussian wave field.  However, WT
approach fails for the amplitudes $s \ge s_{nl}$ for which the
nonlinear time is of the same order or less than the linear wave
period and, therefore, we can expect a cut-off of $P(s)$ at $s =
s_{nl}$. An estimate based on the dynamical equation
(\ref{FourWaveEquationOfMotionB}) gives\footnote{This estimate assumes
that if the wave amplitude at some $k$ happened to be of the critical
value $s_{nl}$ then it will also be of similar value for a range of
$k$'s of width $k$. In other words, strong nonlinearity widens the
$k$-space correlation from zero (RPA value) to $k$ (value for the
coherent structures involved in the wave breaking).}  $s_{nl}=
\omega/\epsilon W k^2 $.  This phenomenological cutoff can be viewed
as a wave breaking process which does not allow wave amplitudes to
exceed their critical value, $P(s) =0$ for $s > s_{nl}$.  Now the
normalization condition can be satisfied for the finite-flux
solutions.  However, having a constant negative flux $F<0$ corresponds
to a source at $s=s_{nl}$ which dictates the necessity of a sink for
some $s < s_{nl}$ to preserve the normalization of $P(s)$. Note
however that the probability sink does not have to correspond to any
physical ``removal'' of waves with certain amplitudes. The sink should
be present solely because the probability is diluted due to acceptance
of new members with $s=s_{nl}$ into the statistical ensemble. In this
case, the sink must be proportional to the probability and, taking
into account the normalization condition, we can write a modified
equation for the PDF in the presence of cutoff,
%
\begin{equation}
\dot P- \partial_{s} (s\gamma P+s\eta\partial_{s}P)= - F_*, \label{peqnM}
\end{equation}
%
with $F_* = - P(s_{nl}) \gamma /s_{nl} $.
The general solution solution to this equation is
%
$
 P = [C  -{{F_*} } Ei({s}/{n} - \log s)/\eta ] \exp{(-s/n)},
$
%
there constant $C$ is fixed by the normalization condition.  This
solution is close to the Rayleigh distribution in the PDF core, $s
\sim n$, and it has a $1/s$ tail at $n \ll s < s_{nl} $.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\section
{\bf Discussion ---}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
We found that the WT intermittency shows as an anomalously high ($\sim
1/s$) probability of the large-amplitude waves whereas at lower
amplitudes distribution appears to be close to Rayleigh ($\sim
e^{-s/n}$) which corresponds to Gaussian wave fields.  We showed that
wave breaking is essential for WT intermittency to be present in the
system, yet the details of wave breaking are not important. The role
of wave breaking is just to ensure that no wave can have amplitude
greater than critical value $s_{nl}$. This simple condition leads to
huge mathematical consequences as it generates the flux solutions in
the amplitude space and therefore creates the $ 1/s$ intermittency.
On the other hand, the amplitude of the $1/s$ tail is not prescribed
by WT and will depend on a particular wave breaking mechanisms in a
particular system.  However, some conclusions about the dependence of
the tail amplitude on the physical parameters can be reached using a
dimensional arguments.

Consider the classical example of the gravity waves on the surface of
deep water.  The linear dispersion relation is given by $\omega_k =
\sqrt{ g k}$, and the coefficient of nonlinear interaction
%
$W^{k{\bf \alpha}}_{\bf \mu\nu}$
%
is given in \cite{ZLF}. This system has two power-law steady state
solutions.  First one is the spectrum corresponding to the direct
cascade of energy toward high-wave numbers, $n_k \propto k^{-4}$
\cite{ZLF,Zakfil}. Second one is the spectrum corresponding to the
inverse cascade of wave action toward the small $k$ values,
$n_k\propto k^{-23/6}$.  In addition to the gravity constant $g$, the
only quantity which determines the state of the system in the direct
cascade range is the energy flux ${\cal P}$ whereas in the inverse
cascade range - the particle flux $\cal Q$. The PDF tail strength can
be characterized by its area which is a dimensionless number and,
therefore, has to depend on the relevant dimensionless combinations in
the direct and the inverse cascade ranges, ${\cal P}^{2/3} k^{1/3} /g$
and ${\cal Q} k/g$ respectively. Thus, the PDF tail thickness grows
with $k$ but its length decreases until it completely disappears at $k
\sim k_{nl}$ (equal to $g^3/{\cal P}^2$ and $g/{\cal Q}$
respectively).

This effect is illustrated in Figure 1 which shows PDF's obtained by
numerical simulations of the dynamical equation for surface gravity
waves on deep water forced at low $k$'s and dissipated at high
$k$'s. Pseudospectral numerical method similar to that of
\cite{yoko},\cite{OsbornePRL} was used on a 256x256 grid. Unlike our
study, however, the works in \cite{yoko},\cite{OsbornePRL} considered
unforced turbulence which was freely decaying and unsteady. Note that
reaching the long-time steadiness was important for us, in order to
accumulate the statistical data necessary for measuring the PDF.

At moderate wavenumber ($k=15 k_{min}$) one can see a PDF tail in the
range $4 n_k < s < 10 n_k$ characterized by an order of magnitude
enhanced probabilities with respect to the Rayleigh distribution.
Unfortunately the range of $s$ where PDF converged to a stable value
in this experiment was not large enough to reach $s \gg n$ values and,
therefore, for an asymptotic scaling to develop.  To increase this
range a much longer computing to gain good statistics of very rare
events at the PDF tail is necessary, which we can not perform with our
resources.
 
At a higher wavenumber ($k=35 k_{min}$) one can see that the large
amplitude waves are less probable than the ones predicted by the
Rayleigh distribution.  This is because the wave breaking happens now
closer to the PDF core causing the PDF cut-off seen at the figure.

In this letter we considered WT which is weak on average so that the
wave breaking occurs only in the PDF tail, i.e. $s_{nl} \gg n$.  It
does not apply to the cases when, at some large $k$, the wave breaking
may become so strong that it occurs for most of the waves in the PDF
core.  These cases where predicted and discussed in \cite{biven}, but
their statistics would be hard to describe analytically because of the
strong nonlinearity.

{\bf Acknowledgments ---} Yeontaek Choi's work is supported by KOSEF
M07-2003-000-10003-0.  YL is supported by NSF CAREER grant DMS 0134955
and by ONR YIP grant N000140210528

\begin{figure}
\epsfxsize=8cm\epsffile{Pokorb.eps}
\caption[]{\footnotesize\baselineskip11pt PDF of $|a_k|^2$ for $k=15
 k_{min}$ (thick curve) and $k=35 k_{min}$ (thin curve) and their
 comparison with Rayleigh distribution (dotted line). Amplitude s is
 normalized so that the two curves have the same slope at $s=0$.}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%% 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 (New York: Consultants Bureau, 1973)
\bibitem{Zakfil} V.E.Zakharov and N.Filonenko, J. Appl. Mech. Tech.
Phys. {\bf 4} 506-515 (1967).
\bibitem{hasselman} K. Hasselmann, J. Fluid Mech {\bf 12} 481 (1962).
\bibitem{zaslav} G.M. Zaslavski and R.Z. Sagdeev, Sov. Phys. JETP {\bf 25}
718 (1967).
\bibitem{OsbornePRL}
M. Onorato {\it et al},
PRL {\bf 89}, 144501, (2002).
\bibitem{JaLubimii} Y.V. Lvov and E.G. Tabak, PRL {\bf  87}, 168501,(2001). 
Y.V. Lvov, K.L. Polzin and E.G.Tabak, PRL, {\bf 92}, 128501 (2004).
\bibitem{DNPZ} S. Dyachenko, A.C. Newell, A. Pushkarev and V.E. Zakharov,
Physica D {\bf 57} 96-160 (1992).
\bibitem{Bradley}
 C. C. Bradley, C. A. Sackett, and R. G. Hulet, PRL  {\bf 78}, 985 (1997).
\bibitem{Pitaevsky} Dalfovo F., Giorgini S., Pitaevsky L., Stringari S., Rev Modern Phys, {\bf 71}, 463, (1999).
\bibitem{zgpd} 
V.E. Zakharov et al,  Physica
D, {\bf 152-153}, 573-619, (2001).
\bibitem{rough} 
A.C. Newell and V.E. Zakharov, 
%``Rough sea foam'', 
PRL {\bf 69}, 1149 (1992).
%\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{jansen} 
P.A.E.M. Janssen, 
%``Nonlinear four-wave interactions and freak waves''
J. Phys Oceanography, {\bf 33} 864 (2003).
\bibitem{biven} 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{cnn}
C. Connaughton S. Nazarenko and A.C. Newell, Physica D, {\bf 184} 86-97 (2003) 
\bibitem{peierls} R. Peierls, Annalen Physik {\bf 3} (1029) 1055.
\bibitem{bp} R. Brout and I. Prigogine, Physica {\bf 22} (1956) 621-636.
\bibitem{prig} I. Prigogine, ``Non-equilibrium Statistical Mechanics'' chap. 2,
Wiley, 1962.
\bibitem{physd} Y. Choi, Y. Lvov, S. Nazarenko, to appear in Physica D, 2005.
\bibitem{pla} Y. Choi, Y. Lvov, S. Nazarenko, Physics Letters A, 332,
230, (2004).
\bibitem{Newell} Per Jakobsen and Alan C Newell J. Stat. Mech.:
Theor. Exp. (2004) L10002
\bibitem{yoko}
N. Yokoyama, J. Fluid Mech. 501 (2004), 169-178.
\bibitem{ln} Y. Lvov and S. Nazarenko, 
%``Statistics of spectral fluctuations in wave turbulence'', 
Phys. Rev. E {\bf 69}, 066608 (2004)      
\bibitem{push}                                                                                                         
V.E. Zakharov, P. Guyenne, A.N. Pushkarev, F. Dias, 
%Wave turbulence in one-dimensional models, 
Physica D {\bf 152-153} 573-619 (2001)
\end{thebibliography}



\end{multicols}

\end{document}
