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 \vert a_k(T)\vert^{2p}\rangle_{\phi,A}. \ p=1, 2, 3 ...,$$

First, let us calculate $\vert a_l(T)\vert^{2 p}$ as

-------------------------------------

$$\vert a_l(T)\vert^{2 p} = \left(a_l^{(0)}+\epsilon a_l^{(1)}+\eps... ..._l\vert^{2p-2}\left(a_l^{(0)}\bar a_l^{(1)}+ \bar a_l^{(0)} a_l^{(1)} \right) +$$

$$\epsilon^2\vert a_l\vert^{2p-4}\Big[ C^2_p(a_l^{(0)} \bar a_l^{(1... ...2 \left( a_l^{(0)} \bar a_l^{(2)}+ \bar a_l^{(0)} a_l^{(2)}\right)\Big] + ... ,$$

----------------------------------- --

where $C^2_p$ is the binomial coefficient.

Up to the second power in $\epsilon$ terms, we have

$$\langle \vert a_l(T)\vert^{2p}\rangle_\phi= \vert a_l\vert^{2p} +$$

$$\epsilon^2 \vert a_l\vert^{2p-2} \left( p^2 \langle \vert a_l^{(1... ... \langle a_l^{(0)} \bar a_l^{(2)}+ \bar a_l^{(0)} a_l^{(2)}\rangle_\phi \right)$$

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 (3) and (5) and the averaging over the phases to obtain

$$\langle \vert a^{(1)}_l\vert^2\rangle_\phi = 4 \sum_{m,n}^\infty ... ...lta^l_{m+n} \vert\Delta(\omega^l_{mn})\vert^2 \vert a_m\vert^2 \vert a_n\vert^2$$

$$+2 \vert V^n_{lm}\vert^2 \vert\Delta^n_{l+m}\vert^2 \delta(\omega^n_{lm}) \vert a_n\vert^2 \vert a_m\vert^2 \big],$$

$$\langle a^{(0)}_l \bar a^{(2)}_l+\bar a^{(0)}_l a^{(2)}_l \rangle... ... \big[ \vert V^l_{mn}\vert^2 \delta^l_{m+n} E(0,\omega^l_{mn}) \vert a_l\vert^2$$

$$+\vert V^n_{lm}\vert^2 \delta^n_{l+m} E(0,\omega^n_{lm}) (\vert a_m\vert^2- \vert a_n\vert^2) \big].$$

Let us substitute these expressions into (6), perform amplitude averaging, take the large box limit and then large $T$ limit ( $T \gg 1/ \omega$). We have

$$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),$$ (4)

with

$$\rho_k = 4 \epsilon^2 \int d {\bf k_1} d {\bf k_2} ( \vert V^k_{12}\vert^2 \delta^k_{12} \delta(\omega^k_{12}) n_{1} n_{2}$$

$$\left. +2 \vert V^2_{k1}\vert^2 \delta^2_{k1} \delta(\omega^2_{k1}) n_{2} n_{1} \right),$$ (5)

$$\gamma_k = 8 \epsilon^2 \int d {\bf k_1} d {\bf k_2} ( \vert V^k_{12}\vert^2 \delta^k_{12} \delta(\omega^k_{12}) n_{2}$$

$$+\vert V^2_{k1}\vert^2 \delta^2_{k1} \delta(\omega^2_{k1}) (n_{1}- n_{2}) ).$$ (6)

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,

$$\dot M^{(p)}_k = -p \gamma_k M^{(p)}_k + p^2 \rho_k M^{(p-1)}_k.$$ (7)

In particular, for the waveaction spectrum $M^{(1)}_k=n_k$ (10) gives the familiar kinetic equation (KE)

$$\dot n_{k} = -\gamma_k n_{k} +\rho_k=\epsilon^2 J(n_{k}),$$ (8)

where $J(n_{k})$ is the ``collision'' term [1,3],

$$J(n_{k})=\int d k_2 d k_1 (R^k_{12}-R^1_{k2}-R^2_{1k}),$$

with

$$R_{k12}=4\pi\vert V^k_{12}\vert^2 \delta^k_{12}\delta(\omega^k_{12}) \Big(n_{2}n_{1}-n_{k}(n_{2}+n_{1})\Big).$$ (9)

The second equation in the series (10) allows to obtain the r.m.s. $\sigma_k^2 = M^{(2)}_k - n^2_k$ of the fluctuations of the waveaction $\langle \vert a_k\vert^2\rangle$. We emphasize that (10) is valid even for strongly intermittent fields with big fluctuations.

Let us now consider the stationary solution of (10), $\dot M^{(p)}_k =0$ for all $p$. Then for $p=1$ from (11) we have $\rho_k=\gamma n_{k}$. Substituting this into (10) we have

$$M^{(p)}_k = p M^{(p-1)}_k n_{k} ,$$

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 (10) 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 (10) to obtain

$$\dot Q^{(p)}_k = - \gamma_k p Q^{(p)}_k + p^2 \rho_k Q^{(p-1)}_k,$$ (10)

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 (13) 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 (10) 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 $\vert a_k\vert$ [11]. It is even easier to test (10) 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 [1,4]. It is given by

$$n_k = A k^{-17/4} ,$$ (11)

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 (10). Note that such no-fluctuations initial conditions were used in [6,7].

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 (9) obtaining:

$$\gamma_k = 8 \epsilon^2 \int d k_1 d k_2 \Delta_{k k_1 k_2}^{-1} ( \vert V^k_{12}\vert^2 \delta^k_{12} \delta(\omega^k_{12}) n_{2}$$

$$+\vert V^2_{k1}\vert^2 \delta^2_{k1} \delta(\omega^2_{k1}) (n_{1}- n_{2}) ).,$$

$$\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 ,$$

$$\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} .$$

Let us substitute ZF spectrum (14) into (), take the values of $\omega_k$ and $V^k_{12}$ appropriate for the capillary waves on deep water([1], 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

$$Q^{(2)}_k = Q^{(2)}_{k0} e^{-2 \gamma_k t},$$

$$Q^{(3)}_k = 9 Q^{(2)}_{k0} n_k e^{-2 \gamma_k t},$$ (12)

etc. Note that fluctuations stabilize at Gaussian values faster for high $k$ values. It is also interesting to test equation (10) when the forcing (and therefore the turbulence) is non-Gaussian, as in most practical situations.