Stochastic averaging, large deviations and random transitions for the dynamics of 2D and geostrophic turbulent vortices

For the barotropic equation (2), the regime corresponding to the emergence of large-scale jets is given by $\alpha \ll 1$ (dynamical time scale of the jets much larger than the time scale of the turbulent fluctuations) and $\nu \ll \alpha$ (turbulent regime). This is the regime we consider in the following.

The large scale zonal jets are charecterized by either a zonal velocity field ${{{\bf v}}_{{\rm z}}}({\bf r})=U(y){{{\bf e}}_{x}}$ or its corresponding zonal potential vorticity ${{q}_{{\rm z}}}(y)=-U^{\prime} (y)+h(y)$. For reasons that will become clear in the following discussion (we will explain that this is a natural hypothesis and prove that it is self-consistent in the limit $\alpha \ll 1$), the non-zonal perturbation to this zonal velocity field is of order $\sqrt{\alpha }$. We then have the decomposition

$$\begin{eqnarray}&&q({\bf r})={{q}_{{\rm z}}}(y)+\sqrt{\alpha }{{\omega }_{m}}({\bf r}),\quad {\bf v}({\bf r})=U(y){{{\bf e}}_{x}}+\sqrt{\alpha }{{{\bf v}}_{m}}({\bf r}),\end{eqnarray} \tag{ 3 }$$

where the zonal projection is defined by $\left\langle f \right\rangle (y)=\frac{1}{2\pi {{l}_{x}}}\int _{0}^{2\pi {{l}_{x}}}{\rm d}x\;f({\bf r})$.

We now project the barotropic equation (2) into zonal and non-zonal part, assuming for simplicity that the random forcing does not act directly on the zonal degrees of freedom¹ $(\left\langle C \right\rangle =0)$

$$\begin{eqnarray}&&\frac{\partial {{q}_{{\rm z}}}}{\partial t}=-\alpha \frac{\partial }{\partial y}\left\langle v_{m}^{(y)}{{\omega }_{m}} \right\rangle -\alpha {{\omega }_{{\rm z}}}+\nu \frac{{{\partial }^{2}}{{\omega }_{{\rm z}}}}{\partial {{y}^{2}}},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\frac{\partial {{\omega }_{m}}}{\partial t}+{{L}_{U}}\left[ {{\omega }_{m}} \right]=\sqrt{2}\eta -\sqrt{\alpha }{{{\bf v}}_{m}}.\nabla {{\omega }_{m}}+\sqrt{\alpha }\left\langle {{{\bf v}}_{m}}.\nabla {{\omega }_{m}} \right\rangle ,\end{eqnarray} \tag{ 5 }$$

where

$$\begin{eqnarray*}&&{{L}_{U}}[{{\omega }_{m}}]=U(y){{\partial }_{x}}{{\omega }_{m}}+{{\partial }_{y}}{{q}_{{\rm z}}}(y){{\partial }_{x}}{{\psi }_{m}}\end{eqnarray*}$$

is the linearized dynamics operator around the zonal base flow U. We see that the zonal potential vorticity is coupled to the non-zonal one through the zonal average of the advection term $\frac{\partial }{\partial y}\langle v_{m}^{(y)}{{\omega }_{m}}\rangle$. If our rescaling of the equations is correct, we clearly see that the natural time scale for the evolution of the zonal flow is $1/\alpha$. By contrast, the natural time scale for the evolution of the non-zonal perturbation is one. These remarks show that in the limit $\alpha \ll 1$, we have a time scale separation between the slow zonal evolution and a rapid non-zonal evolution. Our aim is to use this remark in order to describe precisely the stochastic behavior of the Reynold stress in this limit (by integrating out the non-zonal turbulence), and to prove that our rescaling of the equations and this time scale separation hypothesis is a self-consistent hypothesis.

We will use the remarks that we have a time scale separation between zonal and non-zonal degrees of freedom in order to average out the effect of the non-zonal turbulence. This amounts to treating the zonal degrees of freedom adiabatically. This kind of problems are described in the theoretical physics literature as adiabatic elimination of fast variables (Gardiner 1994) or may also be called stochastic averaging in the mathematics literature. Our aim is to perform the stochastic averaging of the barotropic flow equation and to find the equation that describes the slow evolution of zonal flows. In this stochastic problem, it is natural to work at the level of the probability density function (PDF) of the flow, $P[q]=P[{{q}_{{\rm z}}},{{\omega }_{m}}]$. Then, the dynamical equations (2) or (4) and (5) are equivalent to the so-called Fokker–Planck equation for P.

Complete Fokker–Planck equation

The evolution equation for the PDF reads

$$\begin{eqnarray}&&\frac{\partial P}{\partial t}={{\mathcal{L}}_{0}}P+\sqrt{\alpha }{{\mathcal{L}}_{n}}P+\alpha {{\mathcal{L}}_{{\rm z}}}P,\end{eqnarray} \tag{ 6 }$$

where

$$\begin{eqnarray}&&{{\mathcal{L}}_{0}}P\equiv \int {\rm d}{{{\bf r}}_{1}}\;\frac{\delta }{\delta {{\omega }_{m}}({{{\bf r}}_{1}})}\left[ {{L}_{U}}\left[ {{\omega }_{m}} \right]({{{\bf r}}_{1}})P+\int {\rm d}{{{\bf r}}_{2}}\;{{C}_{m}}({{{\bf r}}_{1}}-{{{\bf r}}_{2}})\frac{\delta P}{\delta {{\omega }_{m}}({{{\bf r}}_{2}})} \right]\end{eqnarray} \tag{ 7 }$$

is the Fokker–Planck operator that corresponds to the linearized dynamics close to the zonal flow U, forced by a Gaussian noise, white in time and with spatial correlations C. This Fokker–Planck operator acts on the non-zonal variables only and depends parametrically on U. This is in accordance with the fact that on time scales of order 1, the zonal flow does not evolve and only the non-zonal degrees of freedom evolve significantly. It should also be remarked that this term contains dissipation terms of order α and ν. These dissipation terms can be included in ${{\mathcal{L}}_{0}}$ because in the limit $\nu \ll \alpha \ll 1$, the non-zonal dynamics is dominated by the interaction with the mean flow, thanks to the so-called Orr mechanism (Orr 1907). This crucial point will be discussed in the following paragraph. At order $\sqrt{\alpha }$, the nonlinear part of the perturbation

$$\begin{eqnarray}&&{{\mathcal{L}}_{n}}P\equiv \int {\rm d}{{{\bf r}}_{1}}\;\frac{\delta }{\delta {{\omega }_{m}}({{{\bf r}}_{1}})}\left[ \left( {{{\bf v}}_{m}}.\nabla {{\omega }_{m}}({{{\bf r}}_{1}})-\langle {{{\bf v}}_{m}}.\nabla {{\omega }_{m}}({{{\bf r}}_{1}})\rangle \right)P \right]\end{eqnarray} \tag{ 8 }$$

describes the nonlinear interactions between non-zonal degrees of freedom. At order α, the zonal part of the perturbation

$$\begin{eqnarray}&&{{\mathcal{L}}_{{\rm z}}}P\equiv \int {\rm d}{{y}_{1}}\;\frac{\delta }{\delta {{q}_{{\rm z}}}({{y}_{1}})}\left[ \left( \frac{\partial }{\partial y}\left\langle v_{m}^{(y)}{{\omega }_{m}} \right\rangle ({{y}_{1}})+{{\omega }_{{\rm z}}}({{y}_{1}})-\frac{\nu }{\alpha }\Delta {{\omega }_{{\rm z}}}({{y}_{1}}) \right)P \right.\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&\left. +\int {\rm d}{{y}_{2}}\;{{C}_{{\rm z}}}({{y}_{1}}-{{y}_{2}})\frac{\delta P}{\delta {{q}_{{\rm z}}}({{y}_{2}})} \right]\end{eqnarray} \tag{ 10 }$$

contains the terms that describe the large-scale friction and the coupling between the zonal and non-zonal flow.

Stationary distribution of the fast variables

The goal of our approach is to get an equation that describes only the zonal, slowly evolving part of the PDF, but taking into account the fact that the non-zonal degrees of freedom have rapidly relaxed to their stationary distribution. The first step is then to determine this stationary distribution of the non-zonal, fastly evolving degrees of freedom. This stationary distribution is given by the stationary state of (6), retaining only the first order term: ${{\mathcal{L}}_{0}}P=0$. For the special case of a determined zonal flow $P\left[ {{q}_{{\rm z}}},{{\omega }_{m}} \right]=\delta \left( {{q}_{{\rm z}}}-{{q}_{0}} \right)Q({{\omega }_{m}})$, ${{\mathcal{L}}_{0}}$ is the Fokker–Planck operator that corresponds to the dynamics of the non-zonal degrees of freedoms, for quasi-geostrophic equations linearized around the base flow with potential vorticity q₀

$$\begin{eqnarray}&&\frac{\partial {{\omega }_{m}}}{\partial t}+{{L}_{U}}[{{\omega }_{m}}]=\sqrt{2}\eta .\end{eqnarray} \tag{ 11 }$$

It is a linear stochastic process (Orstein–Uhlenbeck process) with zero average value, so we know that its stationary distribution is a centered Gaussian, entirely determined by the variance of ω_m. The variance is the stationary value of the two-points correlation function of ω_m, $g({{{\bf r}}_{1}},{{{\bf r}}_{2}},t)={\bf E}\left[ {{\omega }_{m}}({{r}_{1}},t){{\omega }_{m}}({{r}_{2}},t) \right]$. The evolution of g is given by the so-called Lyapunov equation, which is obained by applying the Ito formula to (11)

$$\begin{eqnarray}&&\frac{\partial g}{\partial t}+\left( L_{U}^{(1)}+L_{U}^{(2)} \right)g=2C.\end{eqnarray} \tag{ 12 }$$

($L_{U}^{(i)}$ means that the operator is applied to the ith variable). We now understand that the asymptotic behavior of this equation is a crucial point for the whole theory. It can be proved (Bouchet and Morita 2010, Bouchet et al 2013) that g has a well-defined limit (in the distributional sense) for $t\to \infty$, even in the absence of any dissipation mechanism ($\alpha =\nu =0$). This may seem paradoxical as we deal with a linearized dynamics with a stochastic force and no dissipation mechanism. This is due to the Orr mechanism (Orr 1907, Bouchet and Morita 2010) (the effect of the shear through a non-normal linearized dynamics), that acts as an effective dissipation. The fact that (12) has a finite limit when $\alpha \to 0$ is the precise justification of the scaling (3), and it is thus the central point of the theory.

The average of an observable $A[{{q}_{{\rm z}}},{{\omega }_{m}}]$ over the stationary gaussian distribution is still a function of q_z, and it is an average over the non-zonal degrees of freedom, taking into account the fact that they have relaxed to their stationary distribution. In the following, we denote this average

$$\begin{eqnarray}&&{{{\bf E}}_{U}}\left[ A \right]=\int \mathcal{D}[{{\omega }_{m}}]\;A[{{q}_{{\rm z}}},{{\omega }_{m}}]G\left[ {{q}_{{\rm z}}},{{\omega }_{m}} \right],\end{eqnarray} \tag{ 13 }$$

the subscript U recalling that this quantity depends on the zonal flow. With this definition, the adiabatic reduction can be performed. The details of the computation, that follow Gardiner (1994), are reported in Bouchet et al (2013). Only the final result and its consequences are presented here.

Effective zonal Fokker–Planck equation

The final Fokker–Planck equation for the slowly evolving part of the zonal jets PDF $R[{{q}_{{\rm z}}}]$ reads

$$\begin{eqnarray}&&\frac{\partial R}{\partial \tau }=\int {\rm d}{{y}_{1}}\;\frac{\delta }{\delta {{q}_{{\rm z}}}({{y}_{1}})}\left\{ \left[ \frac{\partial F\left[ U \right]}{\partial {{y}_{1}}}+{{\omega }_{{\rm z}}}({{y}_{1}})-\frac{\nu }{\alpha }\frac{{{\partial }^{2}}{{\omega }_{{\rm z}}}}{\partial y_{1}^{2}} \right]R[{{q}_{{\rm z}}}] \right.\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&\left. +\int {\rm d}{{y}_{2}}\;\frac{\delta }{\delta {{q}_{{\rm z}}}({{y}_{2}})}\left[ {{C}_{{\rm R}}}({{y}_{1}},{{y}_{2}})\left[ {{q}_{{\rm z}}} \right]R\left[ {{q}_{{\rm z}}} \right] \right] \right\},\end{eqnarray} \tag{ 15 }$$

which evolves over the time scale $\tau =\alpha t$, with the drift term

$$\begin{eqnarray*}&&F\left[ U \right]={{{\bf E}}_{U}}\left[ \left\langle v_{m}^{(y)}{{\omega }_{m}} \right\rangle \right]({{y}_{1}})+\alpha {{F}_{1}}\left[ U \right],\end{eqnarray*}$$

with F₁ a functional of q_z, and the diffusion coefficient ${{C}_{{\rm R}}}({{y}_{1}},{{y}_{2}})\left[ {{q}_{{\rm z}}} \right]$, that also depends on the zonal flow q_z.

This Fokker–Planck equation is equivalent to a nonlinear stochastic partial differential equation for the potential vorticity q_z,

$$\begin{eqnarray}&&\frac{\partial {{q}_{{\rm z}}}}{\partial \tau }=-\frac{\partial F}{\partial y}[U]-{{\omega }_{{\rm z}}}({{y}_{1}})+\frac{\nu }{\alpha }\frac{{{\partial }^{2}}{{\omega }_{{\rm z}}}}{\partial {{y}^{2}}}+\zeta ,\end{eqnarray} \tag{ 16 }$$

where ζ is white in time Gaussian noise with spatial correlation C_R. As C_R depends itself on the velocity field U, this is a nonlinear noise. The main physical consequences and the numerical implementation of this equation are discussed in the following paragraphs.

At first order in α, we obtain a deterministic evolution equation for q_z

$$\begin{eqnarray}&&\frac{\partial {{q}_{{\rm z}}}}{\partial t}=-\alpha \frac{\partial }{\partial y}{{{\bf E}}_{U}}\left[ \left\langle v_{m}^{(y)}{{\omega }_{m}} \right\rangle \right]-\alpha {{\omega }_{{\rm z}}}+\nu \Delta {{\omega }_{{\rm z}}}.\end{eqnarray} \tag{ 17 }$$

To summarize, we found that at leading order in α, the zonal flow is forced by the average of the advection term due to the non-zonal fluctuations (Reynolds' stress), and that this quantity is computed from the linearized dynamics for the fluctuations. In other words, we could have applied the same stochastic reduction technique to the quasi-linear dynamics ((4) and (5) without nonlinear terms), and we would have obtained at leading order the same deterministic equation (17). The system (17) and (12) is a quasi-Gaussian (or second-order) closure of the dynamics. Working directly at the level of the PDF, and using the tools of the stochastic reduction, we have been able to justify the closure of this problem. This quasi-Gaussian closure has been already studied in numerical works (SSST in Farrell and Ioannou (2003) and CE2 in Marston et al (2008)) and analytical works (Srinivasan and Young 2012), and is known to give very good results.

Using again the results about the Orr mechanism (Orr 1907, Bouchet and Morita 2010, Bouchet et al 2013) some important facts about equation (17) can be proved. First, we can make sure that the Reynolds' stress is well-behaved, even in the inertial limit $\alpha ,\nu \to 0$, so that the zonal flow equation (17) is always well-defined. We can also show that the energy in the non-zonal degrees of freedom is of order α. As a consequence, a vanishing amount of energy is dissipated in the fluctuations and almost all the energy injected by the stochastic forcing goes to the zonal degrees flow. Moreover, the dynamics defined by (17) and (12) are much simpler to solve numerically than the full nonlinear dynamics (2). A very efficient algorithm used to compute the forcing term $-\alpha \frac{\partial }{\partial y}{{{\bf E}}_{U}}[\langle v_{m}^{(y)}{{\omega }_{m}}\rangle ]$ is presented in the next section 3.2.

In this section we present an efficient algorithm used to compute the forcing term ${{{\bf E}}_{U}}[\langle v_{m}^{(y)}{{\omega }_{m}}\rangle ]$ appearing in the first order equation for the slow evolution of zonal jets (17). We recall that this quantity is the limit for infinite time of the statistical and zonal average of $v_{m}^{(y)}{{\omega }_{m}}$, where ω_m evolves according to the linearized stochastic dynamics (11), with the base flow U held fixed. Equivalently, this quantity can be computed as a linear transform of the infinite-time limit of the solution of the Lyapunov equation (12). In the limit $\nu \ll \alpha \ll 1$, the classical numerical resolution of (12) obtained discretizing the operators $L_{U}^{(1)}$ and $L_{U}^{(2)}$ is a very difficult task. Indeed, the discretization step should be taken smaller and smaller as ν goes to 0. On the contrary, with the algorithm presented below, we can directly take the limit $\nu =0$.

The forcing correlation function can be expanded as ${{C}_{m}}({\bf r})={{\sum }_{k\gt 0,l}}{{c}_{k,l}}{\rm cos} (kx+ly)$, with ${{c}_{k,l}}\geqslant 0$. Then we can decompose the vorticity field as ${{\omega }_{m}}({\bf r},t)=\sum _{k,l=-\infty }^{+\infty }\sqrt{{{c}_{k,l}}}\;{{\omega }_{k,l}}(y,t){{{\rm e}}^{{\rm i}kx}}$ where the coefficients satisfy

$$\begin{eqnarray}&&{{\partial }_{t}}{{\omega }_{k,l}}+{{L}_{U,k}}[{{\omega }_{k,l}}]={{{\rm e}}^{{\rm i}ly}}{{\eta }_{k,l}}(t),\end{eqnarray} \tag{ 18 }$$

with

$$\begin{eqnarray}&&{{L}_{U,k}}={\rm i}kU-{\rm i}k(U^{\prime\prime} -\beta )\Delta _{k}^{-1}+\alpha ,\end{eqnarray} \tag{ 19 }$$

the kth component of the linear operator L_U. The white noises satisfy $\eta _{k,l}^{*}={{\eta }_{-k,-l}}$ and ${\bf E}[{{\eta }_{{{k}_{1}},{{l}_{1}}}}({{t}_{1}})\eta _{{{k}_{2}},{{l}_{2}}}^{*}({{t}_{2}})]={{\delta }_{{{k}_{1}},{{k}_{2}}}}{{\delta }_{{{l}_{1}},{{l}_{2}}}}\delta ({{t}_{1}}-{{t}_{2}})$, and ${{c}_{k,l}}$ is defined for $k\lt 0$ by ${{c}_{k,l}}=c_{-k,-l}^{*}={{c}_{-k,-l}}$. With this decomposition, the Reynolds' stress divergence reads

$$\begin{eqnarray}&&{{{\bf E}}_{U}}\left[ \left\langle v_{m}^{(y)}{{\omega }_{m}} \right\rangle \right]=\mathop{\sum }\limits_{k,l=-\infty }^{+\infty }-{\rm i}k{{c}_{k,l}}\;{{h}_{k,l}}(y)=\mathop{\sum }\limits_{k\gt 0,l}2k{{c}_{k,l}}\;{\rm Im}[{{h}_{k,l}}(y)],\end{eqnarray} \tag{ 20 }$$

with the vorticity-stream function correlation function ${{h}_{k,l}}={{{\bf E}}_{U}}[{{\omega }_{k,l}}\psi _{k,l}^{*}]$. In (Bouchet et al 2013), we presented a numerical scheme for the computation of ${{h}_{k,l}}$ based on an integral equation. This algorithm turns out to be efficient only in very specific cases. Here we present a much more efficient way to compute ${{h}_{k,l}}$, based on the computation of the resolvent of the operator ${{L}_{U,k}}$. From the formal solution of the Ornstein–Uhlenbeck process (18), we have

$$\begin{eqnarray}&&{{h}_{k,l}}(y)=\int _{-\infty }^{\infty }{\rm d}t\;{{\tilde{\omega }}_{k,l}}(y,t)\tilde{\psi }_{k,l}^{*}(y,t),\end{eqnarray} \tag{ 21 }$$

where ${{\tilde{\omega }}_{k,l}}(y,t)$ is the vorticity field obeying the deterministic initial value problem

$$\begin{eqnarray*}&&\forall t\gt 0,\quad {{\partial }_{t}}{{\tilde{\omega }}_{k,l}}(y,t)+{{L}_{U,k}}\left[ {{\tilde{\omega }}_{k,l}} \right](y,t)=0,\end{eqnarray*}$$

$$\begin{eqnarray*}&&{{\tilde{\omega }}_{k,l}}(y,0)={{{\rm e}}^{{\rm i}ly}},\end{eqnarray*}$$

$$\begin{eqnarray*}&&\forall t\lt 0,\quad {{\tilde{\omega }}_{k,l}}(y,t)=0,\end{eqnarray*}$$

and ${{\tilde{\psi }}_{k,l}}(y,t)$ is the associated stream function.

We move now to the frequency domain. The Laplace transform of the vorticity is denoted with

$$\begin{eqnarray}&&{{\hat{\omega }}_{k,l}}(y,c)=\int _{0}^{\infty }{\rm d}t\;{{\tilde{\omega }}_{k,l}}(y,t){{{\rm e}}^{{\rm i}kct}}=\int _{-\infty }^{\infty }{\rm d}t\;{{\tilde{\omega }}_{k,l}}(y,t){{{\rm e}}^{{\rm i}kct}}.\end{eqnarray} \tag{ 22 }$$

For all $\alpha \gt 0$, the deterministic vorticity field ${{\tilde{\omega }}_{k,l}}$ decays to 0 for $t\to \infty$. Then, the Laplace transform defined by (22) for real values of c coincides with the Fourier transform with respect to t. As a consequence, it can be simply inverted as

$$\begin{eqnarray}&&{{\tilde{\omega }}_{k,l}}(y,t)=\frac{|k|}{2\pi }\int _{-\infty }^{\infty }{\rm d}c\;{{\hat{\omega }}_{k,l}}(y,c){{{\rm e}}^{-{\rm i}kct}}.\end{eqnarray} \tag{ 23 }$$

We stress that this property is valid only for decaying fields ${{\tilde{\omega }}_{k,l}}$, and thus only when $\alpha \ne 0$. It is also useful to define the Laplace transform of the stream function ${{\hat{\psi }}_{k,l}}$; this quantity is usually referred in literature as the resolvent of the operator ${{L}_{U,k}}$. It is related to ${{\hat{\omega }}_{k,l}}$ through ${{\hat{\omega }}_{k,l}}(y,c)=\left( \frac{{{{\rm d}}^{2}}}{{\rm d}{{y}^{2}}}-{{k}^{2}} \right){{\hat{\psi }}_{k,l}}(y,c)$ and is the solution of the linear ordinary differential equation

$$\begin{eqnarray}&&\left( \frac{{{{\rm d}}^{2}}}{{\rm d}{{y}^{2}}}-{{k}^{2}} \right){{\hat{\psi }}_{k,l}}(y,c)-\frac{U^{\prime\prime} (y)-\beta }{U(y)-c-{\rm i}\frac{\alpha }{k}}{{\hat{\psi }}_{k,l}}(y,c)=\frac{{{{\rm e}}^{{\rm i}ly}}}{{\rm ik}\left( U(y)-c-{\rm i}\frac{\alpha }{k} \right)}.\end{eqnarray} \tag{ 24 }$$

Using (23) and (24) in (21), and performing the integration over t, we get

$$\begin{eqnarray}&&{{h}_{k,l}}(y)=\frac{|k|}{2\pi }\int _{-\infty }^{\infty }{\rm d}c\;\frac{(U^{\prime\prime} (y)-\beta ){{{\hat{\psi }}}_{k,l}}(y,c)+{{{\rm e}}^{{\rm i}ly}}/{\rm i}k}{U(y)-c-{\rm i}\frac{\alpha }{k}}\;\hat{\psi }_{k,l}^{*}(y,c).\end{eqnarray} \tag{ 25 }$$

The numerical computation of the resolvant ${{\hat{\psi }}_{k,l}}$ is a very easy task, an algorithm is detailed in Bouchet and Morita (2010). For a given value of c, computing ${{\hat{\psi }}_{k,l}}(y,c)$ takes a few seconds on a laptop. Then, the computation of the integral (25) is a matter of a few minutes, with a very good accuracy (of the order of 10⁻³). The fact that we are able to compute the infinite-time limit of a statistical average, in the limit of no viscosity, so fast and with such accuracy, is an important result in itself. An example of application of this method is given in figure 2. Using (25) for smaller and smaller values of α allows to check direclty the theoretical statements of the previous paragraph 3.1: the Reynolds' stress divergence ${{{\bf E}}_{U}}[\langle v_{m}^{(y)}{{\omega }_{m}}\rangle ]$ converges to a well-defined function in the limit $\alpha \to 0$, and all the energy is transferred from the forcing to the zonal flow in this limit.

Zoom In Zoom Out Reset image size

From the full Fokker–Planck equation (6), we expect the non-linear operator ${{\mathcal{L}}_{n}}$ to produce terms of order ${{\alpha }^{1/2}}$ and ${{\alpha }^{3/2}}$ in the zonal Fokker–Planck equation (14). The detailed computation shows that these terms exactly vanish (Bouchet et al 2013). As a consequence, we have proved that the quasi-Gaussian closure (17,12) is correct in the limit $\alpha \ll 1$, with correction only at order α².

We then have a correction F₁ to the drift $F[U]$ due to the nonlinear interactions. At this order, the quasilinear dynamics and nonlinear dynamics differ. We also see the appearance of the noise term, which has a qualitatively different effect than the drift term. For instance if one is interested in large deviations from the most probable states, correction of order α to F₀ will still be vanishingly small, whereas the effect of the noise will be essential. This issue is important for the description of the bistability of zonal jets and phase transitions.