# Phase and amplitude evolution in the network of triadic interactions of the Hasegawa–Wakatani system

Ö. Gürcan, J. Anderson, S. Moradi, A. Biancalani, P. Morel

► To cite this version:

Ö. Gürcan, J. Anderson, S. Moradi, A. Biancalani, P. Morel. Phase and amplitude evolution in the network of triadic interactions of the Hasegawa–Wakatani system. Physics of Plasmas, 2022, 29 (5), pp.052306. 10.1063/5.0089073. hal-03865423

HAL Id: hal-03865423

https://hal.science/hal-03865423

Submitted on 22 Nov 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Phase and amplitude evolution in the network of triadic interactions of the Hasegawa-Wakatani system

Ö. D. Gürcan, J. Anderson, S. Moradi, A. Biancalani, P. Morel

CNRS, LPP, Ecole Polytechnique

Hasegawa-Wakatani system, commonly used as a toy model of dissipative drift waves in fusion devices is revisited with considerations of phase and amplitude dynamics of its triadic interactions. It is observed that a single resonant triad can saturate via three way phase locking where the phase differences between dominant modes converge to constant values as individual phases increase in time. This allows the system to have approximately constant amplitude solutions. Non-resonant triads show similar behavior only when one of its legs is a zonal wave number. However when an additional triad, which is a reflection of the original one with respect to the y axis is included, the behavior of the resulting triad pair is shown to be more complex. In particular, it is found that triads involving small radial wave numbers (large scale zonal flows) end up transferring their energy to the subdominant mode which keeps growing exponentially, while those involving larger radial wave numbers (small scale zonal flows) tend to find steady chaotic or limit cycle states (or decay to zero). In order to study the dynamics in a connected network of triads, a network formulation is considered including a pump mode, and a number of zonal and non-zonal subdominant modes as a dynamical system. It was observed that the zonal modes become clearly dominant only when a large number of triads are connected. When the zonal flow becomes dominant as a 'collective mean field', individual interactions between modes become less important, which is consistent with the inhomogeneous wave-kinetic picture. Finally, the results of direct numerical simulation is discussed for the same parameters and various forms of the order parameter are computed. It is observed that nonlinear phase dynamics results in a flattening of the large scale phase velocity as a function of scale in direct numerical simulations.

# I. INTRODUCTION

Two dimensional Hasegawa-Wakatani equations[1] with proper zonal response consist of an equation of plasma vorticity

$$ \frac{\partial}{\partial t} \nabla^2 \Phi + \hat{\mathbf{z}} \times \nabla \Phi \cdot \nabla \nabla^2 \Phi = C(\tilde{\Phi} - \tilde{n}) + D_{\Phi} (\nabla^2 \Phi) \quad (1) $$

and an equation of continuity

$$ \frac{\partial}{\partial t} n + \hat{\mathbf{z}} \times \nabla \Phi \cdot \nabla n + \kappa \partial_y \Phi = C(\tilde{\Phi} - \tilde{n}) + D_n(n), \quad (2) $$

with the $E \times B$ velocity defined as $v_E = \hat{\mathbf{z}} \times \nabla \Phi$ in normalized form and $\tilde{\Phi} = \Phi - \langle \Phi \rangle$ where $\langle \Phi \rangle$ denotes averaging in y (i.e. poloidal) direction. Here $n$ is the fluctuating particle density normalized to a background density $n_0$, $\Phi$ is the electrostatic potential normalized to $T/e$, $\kappa$ is the diamagnetic velocity normalized to speed of sound, $C$ is the so called adiabaticity parameter, which is a measure of the electron mobility and $D_\Phi$ and $D_n$ are dissipation functions for vorticity and particle density respectively. For fluctuations we have $D_\Phi(\nabla^2 \tilde{\Phi}) = \nu \nabla^4 \tilde{\Phi}$ from kinematic viscosity, whereas for the zonal flows $D_\Phi(\nabla^2 \tilde{\Phi}) = -\nu_{ZF} \nabla^2 \tilde{\Phi}$ from large scale friction. Unless the system represents a renormalized formulation, $D_n$ should actually be zero, however here we include it for completeness and numerical convenience and take it to have the same form as the vorticity dissipation with diffusion $D_n(\tilde{n}) = D \nabla^2 \tilde{n}$ and particle loss $D(\tilde{n}) = -D_{ZF} \tilde{n}$.

The Hasegawa Wakatani model, was initially devised as a simple, nonlinear model of dissipative drift wave tur-

bulence in tokamak plasmas. It has the same nonlinear structure as the passive scalar turbulence [2] -with vorticity evolving according to 2D Navier-Stokes equations- or more complex problems such as rotating convection [3, 4]. From a plasma physics perspective it can be considered as the minimum non-trivial model for plasma turbulence, since it has i) linear instability (e.g. Hasegawa-Mima model does not [5]), ii) finite frequency (so that resonant interactions are possible [6]), and iii) a proper treatment of zonal flows[7]. The model is well known to generate high levels of large scale zonal flows, especially for $C \gtrsim 1$ [8–10]. It has been studied in detail for many problems in fusion plasmas including dissipative drift waves in tokamak edge [11, 12], subcritical turbulence[13], trapped ion modes [14], intermittency [15, 16], closures [17–19], feedback control [20], information geometry [21] and machine learning [22]. Variations of the Hasegawa-Wakatani model are regularly used for describing turbulence in basic plasma devices[23–25].

Formation of large scale structures, in particular Zonal flows in drift wave turbulence is one of the key issues in the study of turbulence in fusion plasmas, which can be formulated in terms of modulational instability of either a gas of drift wave turbulence using the wave kinetic formulation [26] or a small number of drift modes[27], resulting in various forms of complex amplitude equations such as the celebrated nonlinear Schrödinger equation (NLS) [28]. It is also common to talk about zonal flows as resulting from a process of inverse cascade[29, 30] and their back reaction on turbulence[31, 32] results in predator-prey dynamics, possibly leading up to the low to high
confinement transition in tokamaks [33, 34]. While the role of the complex phases in nonlinear evolution of the amplitudes, especially in the context of structure formation, for example as in the case of soliton formation in NLS, was always well known, its particular importance for zonal flow formation in toroidal geometry has been underlined recently [35].

Here we revisit the Hasegawa-Wakatani system, with proper zonal response, as a minimum system that allows a description of zonal flow formation in drift wave turbulence, and study interactions between various number of modes from three wave interactions to the full spectrum of modes described by direct numerical simulations, focusing in particular on phase dynamics and the possibility of phase locking and synchronization. It turns out that while resonant three wave interactions involving unstable and damped modes favor phase locking (i.e. a state where the differences between individual phases remain roughly constant as they increase together), interactions involving zonal flows (i.e. four wave interactions including the triad reflected with respect to the y axis), seems to have a complicated set of possible outcomes depending on if the zonal flow wave number is larger or smaller than the pump wave-number. It therefore becomes critical to study a “network” of connected triads in order to see the collective effects of a number of triads on the evolution of zonal flows and of the relative phases between modes. Two different network configurations are considered: that with a single $k_y$ and many different q's, and that with a single q but many different $k_y$'s. Note that the algorithm that we use computes all possible interactions between the modes in a given collection of triads and then computes the interaction coefficients and evolves the system nonlinearly according to those.

Finally we consider the results from direct numerical simulations (DNS) using a pseudo-spectral 2D Hasegawa-Wakatani solver. The DNS and the network models correspond exactly, in the sense that if we consider an $N_x \times N_y$ grid and consider all the possible triads in such a grid and solve this problem using our network solver, we obtain exactly the same problem (including the boundary conditions that are periodic) as the DNS. The results of the DNS show qualitatively similar behavior to the two network models that we considered. However looking at the evolution of phases, we observe a nonlinear flattening of the phase velocity for large scales computed as a function of x, suggesting nonlinear structure formation in the classical sense of nonlinearity balancing dispersion resulting in a constant velocity propagation at least for large r scale structures. These vortex-like structures that move at a constant velocity are also clearly visible in the time evolution of density and vorticity fields.

The rest of the paper is organized as follows. In the remainder of the introduction, the Hasegawa-Wakatani system is reformulated in terms of its linear eigenmodes and the amplitude and phase equations for these eigen-

modes, writing out explicitly the nonlinear terms that appear in this formulation. In Section II, different types of interactions among dissipative drift waves are considered using these linear eigenmodes, starting with the basic three wave interaction. After showing that there is no qualitative difference between a near resonant and an exactly resonant (within numerical accuracy) triad, the details of the phase dynamics of such a single triad are discussed. In Section III, the interaction with zonal flows are considered. It is noted that when we consider a triad and its reflection with respect to its pump wave-number together as a pair, the behavior of the system is qualitatively different from the single triad case. After a discussion of order parameters for this system, a network formulation is considered and the results from such a network model is presented. Finally the results from direct numerical simulations of the Hasegawa-Wakatani system is discussed and compared with those earlier results based on reduced number of triads. Section IV is conclusion.

### A. Linear Eigenmodes

We can write the Hasegawa-Wakatani system in Fourier space for non-zonal modes (i.e. $k_y \neq 0$) as follows:

$$ \partial_t \Phi_k + (A_k - B_k) \Phi_k = \frac{C}{k^2} n_k + N_{\Phi k} \quad (3) $$

$$ \partial_t n_k + (A_k + B_k) n_k = (C - i\kappa k_y) \Phi_k + N_{nk} \quad (4) $$

where

$$ A_k = \frac{1}{2} \left[ (Dk^2 + C) + \left( \frac{C}{k^2} + \nu k^2 \right) \right] \quad (5) $$

and

$$ B_k = \frac{1}{2} \left[ (Dk^2 + C) - \left( \frac{C}{k^2} + \nu k^2 \right) \right] . \quad (6) $$

using the simpler notation $\Phi_k \rightarrow \Phi_k$ and defining:

$$ N_{nk} = \frac{1}{2} \sum_{\Delta} \hat{\mathbf{z}} \times \mathbf{p} \cdot \mathbf{q} (\Phi_p^* n_q^* - \Phi_q^* n_p^*) \quad (7) $$

and

$$ N_{\Phi k} = \frac{1}{2} \sum_{\Delta} \frac{\hat{\mathbf{z}} \times \mathbf{p} \cdot \mathbf{q} (q^2 - p^2) \Phi_p^* \Phi_q^*}{k^2} . \quad (8) $$

Diagonalizing the linear terms we can write:

$$ \partial_t \xi_k^\pm + i\omega_k^\pm \xi_k^\pm = N_{\xi k}^\pm \quad (9) $$
with the complex eigen-frequencies $\omega_k^{\pm} = \omega_{rk}^{\pm} + i\gamma_k^{\pm}$ that can be written as:

$$\omega_k^{\pm} = \Omega_k^{\pm} - iA_k$$

with

$$\Omega_k^{\pm} = \pm \left( \sigma_k \sqrt{\frac{H_k - G_k}{2}} + i \sqrt{\frac{H_k + G_k}{2}} \right) \quad (10)$$

where $\sigma_k = \text{sign}(\kappa k_y)$,

$$H_k = \sqrt{G_k^2 + C^2 \kappa^2 k_y^2 / k^4}, \quad (11)$$

and

$$G_k = \left( B_k^2 + \frac{C^2}{k^2} \right). \quad (12)$$

This allows us to write the two linear eigenmodes as:

$$\xi_k^{s_k} = n_k + \frac{k^2}{C} [B_k - i\Omega_k^{s_k}] \Phi_k . \quad (13)$$

where $s_k = \pm$. The nonlinear terms in (9) become:

$$N_{\xi_k}^{s_k} = N_{nk} + \frac{k^2}{C} (B_k - i\Omega_k^{s_k}) N_{\Phi_k}, \quad (14)$$

and the inverse transforms can be written as:

$$\Phi_k = \frac{i}{2} \frac{C}{k^2} \sum_{s_k} \frac{\xi_k^{s_k}}{\Omega_k^{s_k}} \quad (15)$$

$$n_k = -\frac{i}{2} \sum_{s_k} \frac{1}{\Omega_k^{s_k}} (B_k + i\Omega_k^{s_k}) \xi_k^{s_k} . \quad (16)$$

Considering the inviscid limit, $\{D, \nu\} \to 0$ and $k_y \to O(\epsilon)$, where we keep terms only up to $O(\epsilon)$ we obtain:

$$\xi_k^+ = n_k + \left( k^2 - i \frac{\kappa k_y}{2A_k^2} \right) \Phi_k$$

$$\xi_k^- = n_k - \left( 1 + i \frac{\kappa k_y}{2A_k^2} \right) \Phi_k,$$

which means that one could loosely refer to these two modes as the potential vorticity mode (i.e. $\xi_k^+ = n_k + k^2\Phi_k$) and the non-adiabatic electron density mode (i.e. $\xi_k^- = n_k - \Phi_k$), somewhat similar to the real space decomposition used in Ref. 36. Since the equations are already diagonal for $k_y = 0$ modes, we can use $\xi_k^+ = k^2\Phi_k$ and $\xi_k^- = n_k$ for those (or $\bar{\Phi}_k$ and $\bar{n}_k$ explicitly as we will do below).

Notice that the two eigenmodes in (13) are not orthogonal. They have the same frequencies (in opposite directions) but different growth rates with $\gamma_k^+ > \gamma_k^-$ (with

$\gamma_k^- < 0$, while $\gamma_k^+$ can be positive or negative depending on the wave-number). The full nonlinear initial value problem can be solved using linear eigenmodes by first computing $\xi_k^{s_k}(0)$ from (13), and then advancing those to $\xi_k^{s_k}(t)$ using (9), where the linear matrix is now diagonal (but the nonlinear coupling terms are rather complicated), and finally going back to compute $\Phi_k(t)$ and $n_k(t)$ using (15-16). Obviously, this approach does not involve any kind of approximation.

## B. Amplitude and Phase Equations

Substituting $\xi_k^{\pm} = \chi_k^{\pm} e^{i\phi_k^{\pm}}$ into (9), we get:

$$\partial_t (\chi_k^{\pm} e^{i\phi_k^{\pm}}) + i\omega_k^{\pm} \chi_k^{\pm} e^{i\phi_k^{\pm}} = \left|N_{\xi_k}^{\pm}\right| e^{i\phi_k^{\pm}} \quad (17)$$

taking the real part we obtain the amplitude equations:

$$(\partial_t - \gamma_k^{\pm}) \chi_k^{\pm} = \left|N_{\xi_k}^{\pm}\right| \cos \left( \phi_k^{N_{\xi_k}^{\pm}} - \phi_k^{\pm} \right) \quad (18)$$

and taking the imaginary part and dividing by $\chi_k^{\pm}$ we get the phase equations:

$$\partial_t \phi_k^{\pm} = -\omega_{kr}^{\pm} + \frac{\left|N_{\xi_k}^{\pm}\right|}{\chi_k^{\pm}} \sin \left( \phi_k^{N_{\xi_k}^{\pm}} - \phi_k^{\pm} \right) . \quad (19)$$

The form of the amplitude equation (18) means that the fixed point for the amplitude evolution is determined by the phase difference between $N_{\xi_k}^{s_k}$ and $\xi_k^{s_k}$ for each $s_k$. However such a fixed point keeps evolving since the phases themselves increase linearly with the linear frequency while being deformed by the nonlinear terms. Note that if the nonlinear phase is dominated by a slowly evolving mean phase (could be the case if the nonlinear interactions are dominated by the interactions with a zonal flow), the individual phases will be attracted to this nonlinear mean phase, since if the individual phase is behind the nonlinear phase, the $\sin(\phi_k^{N_{\xi_k}^{\pm}} - \phi_k^{\pm})$ will be positive, causing the individual phase to accelerate, whereas if the individual phase is ahead of the nonlinear phase it will be slowed down. However since we have linear frequencies it is impossible for individual phases to become phase locked directly with the slow nonlinear phase. Instead the nonlinear term plays a role akin to that of the ponderomotive force in parametric instability.

## C. Nonlinear Terms

In order to compute $N_{\xi_k}^{\pm}$ in terms of $\xi_k^{\pm}$, we need to go back to $\Phi_k$ and $n_k$ using (15-16), compute the nonlinear terms (7-8) using those and combine them as in (14). They can then be written in the form:
<figure><img src="image_2.png" /><figcaption>Figure 1. The resonance manifold $\Delta\omega = \omega_{kr}^+ + \omega_{pr}^+ + \omega_{qr}^+ = 0$ of the Hasegawa-Wakatani system for the case $C = 1.0$, $\kappa = 0.2$, $\nu = D = 10^{-3}$ is shown corresponding to the wave vector $\mathbf{k} = (5, 5)$ that is shown explicitly. Any $\mathbf{p}$ that falls onto the region inside the resonance manifold (shown here with a finite width of $\pm 0.04$ with $\Delta\omega > 0$ in red and $\Delta\omega < 0$ in blue if in color) gives $\Delta\omega \approx 0$ (with $\mathbf{q} = -\mathbf{k} - \mathbf{p}$). As discussed in the text, because of the fact that the (+) and (-) modes have the same frequency (but opposite direction of propagation in $y$ direction) all possible combinations of (+) and (-) modes resonate on the same manifold.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 2. The resonance manifold, shown on top of the growth rate where red corresponds to $\gamma_k^+ > 0$ and blue to $\gamma_k^+ < 0$ for a) the most unstable model on the grid $\mathbf{k} = (0, 1.125)$, b) a nearby mode with a small $k_x$ component $\mathbf{k} = (0.250, 1.125)$, c) a mode with $k_x = k_y$ that is $\mathbf{k} = (1.125, 1.125)$ and finally d) a mode that has $k_x \gg k_y$ with $\mathbf{k} = (1.125, 0.125)$.</figcaption></figure>

where $\theta_{M_{kpq}}^{s_k s_p s_q}$ is the phase of the nonlinear interaction coefficient $M_{\xi_{kpq}}^{s_k s_p s_q}$.

$$N_{\xi k}^{s_k} = \frac{1}{2} \sum_{\Delta} \sum_{s_p, s_q} M_{kpq}^{s_k s_p s_q} \xi_p^{s_p*} \xi_q^{s_q*} \quad (20)$$

in terms of the linear eigenmodes, where the sum is over $s_p, s_q = \{(+, +), (+, -), (-, +), (-, -)\}$ for $s_k = (+, -)$. The nonlinear interaction coefficients in (20) can be written (i.e. between 3 non-zonal modes) as:

$$M_{\xi_{kpq}}^{s_k s_p s_q} = m_{kpq}^{s_k s_p s_q} \left[ q^2 (B_q - i\Omega_q^{s_q*}) - p^2 (B_p - i\Omega_p^{s_p*}) - (q^2 - p^2) (B_k - i\Omega_k^{s_k*}) \right] \quad (21)$$

where

$$m_{kpq}^{s_k s_p s_q} = \frac{C \hat{z} \times \mathbf{p} \cdot \mathbf{q}}{4 \Omega_p^{s_p*} \Omega_q^{s_q*} q^2 p^2}$$

and $\Omega_k^\pm$ is given in (10).

Note that these coefficients are complex, and have different phases in general. In other words the explicit forms of (19) can be written as:

$$\partial_t \phi_k^{s_k} = -\omega_{kr}^{s_k} + \sum_{\Delta} \sum_{s_p, s_q} \frac{|M_{\xi_{kpq}}^{s_k s_p s_q}| |\xi_p^{s_p}| |\xi_q^{s_q}|}{|\xi_k^{s_k}|} \times \sin \left( \theta_{M_{kpq}}^{s_k s_p s_q} - \phi_p^{s_p} - \phi_q^{s_q} - \phi_k^{s_k} \right) \quad (22)$$

## II. INTERACTIONS AMONG DRIFT WAVES

### A. Three wave interactions

Consider three separate modes $k, p$ and $q$ that satisfy the triadic interaction condition $\mathbf{k}+\mathbf{p}+\mathbf{q}=0$, possibly in the presence of other modes. The nonlinear term for the wave number $k$ can then be divided into the interaction with the pair $p$ and $q$, and the interaction with the rest of the modes (if they exist). If the three wave interaction that we consider is resonant, slightly off-resonance, or completely non-resonant, its evolution is likely to be different, which can be considered as different scenarios. It may also be possible to model the effects of rest of the modes as background forcing, modification of the linear terms ($\partial_t la$ eddy damping) or simply as stochastic noise. Thus separating the nonlinear term into the interaction with the pair $p$ and $q$ (i.e. $N_{\xi_{kpq}}^\pm$) and the interaction with the rest of the modes (i.e. $\delta N_{\xi_{kpq}}^\pm$), we can write:

$$\partial_t \xi_k^\pm + i\omega_k^\pm \xi_k^\pm = N_{\xi_{kpq}}^\pm + \delta N_{\xi_{kpq}}^\pm \quad (23)$$

where

$$N_{\xi_{kpq}}^{\pm} = M_{\xi_{kpq}}^{\pm+++} \xi_p^{+*} \xi_q^{+*} + M_{\xi_{kpq}}^{\pm+--} \xi_p^{+*} \xi_q^{-*} + M_{\xi_{kpq}}^{\pm-+-} \xi_p^{-*} \xi_q^{+*} + M_{\xi_{kpq}}^{\pm+-} \xi_p^{-*} \xi_q^{-*}$$
with $M_{\xi k p q}^{\pm \pm}$ being (complex) nonlinear interaction coefficients, and

$$ \delta N_{\xi k p q}^{\pm} = N_{\xi k}^{\pm} - N_{\xi k p q}^{\pm}. $$

The p and q modes evolve similarly:

$$ \partial_t \xi_p^{\pm} + i \omega_p^{\pm} \xi_p^{\pm} = N_{\xi pqk}^{\pm} + \delta N_{\xi pqk}^{\pm} \quad (24) $$

$$ \partial_t \xi_q^{\pm} + i \omega_q^{\pm} \xi_q^{\pm} = N_{\xi qkp}^{\pm} + \delta N_{\xi qkp}^{\pm} \quad (25) $$

Notice that, since there are two eigenmodes (23-25) represent 6 equations. One can therefore consider resonances between 3 growing modes, 2 growing modes and a damped mode, or a growing mode and 2 damped modes etc. However since the frequencies are the same with opposing signs, and due to the condition that the flow field is real, we have both $k_y$ and $-k_y$ components, whenever we have a resonance say of the form $\omega_k^+ + \omega_p^+ + \omega_q^+ = 0$, (with $k+p+q=0$), we also have $\omega_k^- + \omega_p^- + \omega_q^- = 0$, $\omega_k^+ - \omega_{-p}^- - \omega_{-q}^- = 0$ or $\omega_{-k}^- - \omega_p^+ - \omega_q^+ = 0$ etc. In other words, whenever we have a resonance for three + modes, we also have all the other combinations. The form of the resonance manifold can be seen in figures 1 and 2, for $C=1$, $\kappa=0.2$, and $\nu = D = 10^{-3}$, which we will refer to as the “C = 1 case”.

The three wave interaction system (23-25) can be implemented numerically without much difficulty by dropping the $\delta N_\xi$ terms above. One can also formulate the same three wave interaction problem in the original variables $\Phi_k$, $\Phi_p$, $\Phi_q$, $n_k$, $n_p$ and $n_q$ using the form (3-4) before the transformation, and then transform the result using (13). Obviously those two approaches are numerically equivalent and naturally they give exactly the same results. We used this to verify that the eigenmode computation was correct. While in general it is unclear if the eigenmode formulation provides any concrete advantage apart from diagonalizing the linear system, the advantage becomes self-evident if the resulting fluctuations have $\xi_k^+ \gg \xi_k^-$ and we can drop the $\xi_k^-$ mode for example.

1. Is there a difference between exact and near resonances?

We first pick a primary wave-number $k = (0, 1.125)$ which is the linearly most unstable mode on a grid with $dk_x = dk_y = 0.125$ for the $C = 1.0$ case and consider the resonance manifold as shown in figure 2a in order to pick a second wave-number $p = (-0.5, -1.0)$ as the point on the k-space grid that is closest to the resonance manifold. The third wave-number $q$ is computed from $k+p+q=0$. While a direct numerical simulation only has the wave-numbers on grid points, a three wave equation solver is not constrained to such a grid. We can instead compute $p$ to be exactly on the resonance manifold -at least

<figure><img src="image_4.png" /><figcaption>Figure 3. Comparison between exact or near resonances, with real parts of each eigenmode shown for each wave number as labeled on the left side of the figure. The solid line is the exact (i.e. $\Delta\omega \approx 2 \times 10^{-15}$) resonance of $k = (0, 1.125)$ with $p = (-0.5, -1.0632325265492)$ whereas the dashed line is the near resonance with $p = (-0.5, -1.0)$ and $\Delta\omega \approx 0.01$. While some details change, the overall behavior, and saturation levels are actually very similar.</figcaption></figure>

within some numerical precision- for example by choosing $p = (-0.5, -1.0632325265492)$. Solving the three wave equations numerically, using these slightly different sets of wave-numbers, we find that having exact resonance or near resonance (i.e. $\Delta\omega \approx 2 \times 10^{-15}$ vs. $\Delta\omega \approx 0.01$) does not make much difference in terms of time evolution (see figure ), while picking something like $p = (-0.5, -1.5)$, which gives $\Delta\omega \approx 0.07$ (with $\omega_k \approx 0.1$ for comparison) gives a completely different time evolution, where one of the modes keeps growing linearly without being able to couple to the other two. We verified this for a bunch of different sets of wave numbers, and while there are some differences in detail, generally both exactly resonant or near resonant triads lead to saturation but non-resonant triads can not saturate, possibly due to lack of efficient interactions. The boundary between what can be considered a near resonant vs. non-resonant interaction can actually be defined using this criterion. In particular, it seems that the triads with one of the frequencies much smaller than the other two (i.e. $\omega_q \ll \omega_p \sim \omega_k$) tend to support larger overall $\Delta\omega$, and nonetheless reach saturation. However it is not clear whether these observations from a single triad continue to hold when many triads are interacting with each other.

B. Phase Evolution

Considering the (unwrapped) phase evolution of each of the modes of the near resonant triad with $k = (0, 1.125)$ and $p = (-0.5, -1.0)$, we observe that while
<figure><img src="image_5.png" /><figcaption>Figure 4. Time evolution of the amplitudes of the eigenmodes for $C = 1$ case with $\mathbf{k} = (0, 1.125)$ and $\mathbf{p} = (-0.5, -1.0)$. We have a “saturated” state with oscillating amplitudes. It seems that as $k$ and $p$ (the two unstable modes and the two larger legs of the triads) exchange energy, $q$ plays the role of the mediator.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5. Time evolution of the phases $\varphi_k^{\pm}$ and their sums $\psi_{kpq}^{s_p s_q s_k}$ for $C = 1$ case with $\mathbf{k} = (0, 1.125)$ and $\mathbf{p} = (-0.5, -1.0)$. Saturation of the amplitudes as seen in figure 4 is accompanied by a nonlinear frequency shift as shown in the top plot and the saturation of the $\psi_{kpq}^{s_p s_q s_k}$'s as shown in the bottom plot. Note that $\psi_{kpq}^{s_p s_q s_k} = \text{const.}$ would correspond to phase locking.</figcaption></figure>

some amplitude evolution continues, the phases converge towards straight lines, implying more or less constant frequencies in the final stage. These nonlinear frequencies are substantially shifted with respect to the initial linear frequencies due to the effect of nonlinear terms. However it appears that the system remains in resonance as the sum of the final nonlinear frequencies remain very close to zero. In fact, it appears that the “near resonant” system evolves towards resonance as a result of these nonlinear corrections, since $\Delta\omega$ decreases from the beginning towards the end.

Using (18) and (19) with the assumption that $\partial_t \chi_k^{\pm} \approx 0$

<figure><img src="image_7.png" /><figcaption>Figure 6. Time derivatives of the phases $\varphi_k^{\pm}$ for $C = 1$ case with $\mathbf{k} = (0, 1.125)$ and $\mathbf{p} = (-0.5, -1.0)$, corresponding to nonlinear frequencies. Notice that while $d\varphi_q^- / dt$ appears to oscillate wildly, since its amplitude $\chi_q^-$ is vanishingly small, as can be seen in figure 4, these oscillations are not important for the rest of the dynamics.</figcaption></figure>

and $\partial_t \varphi_k^{\pm} = -\omega_{k,nl}^{\pm}$ is a constant, we obtain the nonlinear frequency shift, i.e. $\delta\omega_{kr}^{\pm} = \omega_{k,nl}^{\pm} - \omega_{kr}^{\pm}$ as:

$$ \delta\omega_{kr}^{\pm} = \text{sign}(\omega_{kr}) \sqrt{\frac{|N_{\xi k}^{\pm}|^2}{|\xi_k^{\pm}|^2} - \gamma_k^{\pm 2}}, \quad (26) $$

which can be computed given the final amplitudes and the nonlinear interaction coefficients (21). For example for the case above the smoothed saturated amplitudes are shown in the table:

In order to elucidate dynamics of the phases in a triad, we define the sums of phases as a separate variable following Ref. 37:

$$ \psi_{kpq}^{s_k s_p s_q} \equiv \varphi_k^{s_k} + \varphi_p^{s_p} + \varphi_q^{s_q}. \quad (27) $$

We observe that while the phases keep increasing in time, for a steady state, the phase differences should remain bounded. We can write the equations for the amplitudes as

<table><caption>Table I. Saturated amplitudes, linear frequencies, linear growth rates, the final nonlinear frequencies and the $\delta\omega$'s that are computed from (26), rounded to two significant figures for the $C = 1$ case with $\mathbf{k} = (0, 1.125)$ and $\mathbf{p} = (-0.5, -1.0)$. Note that the basic assumption of (26) works only for linearly unstable modes, and for those $\delta\omega$ is not far from $\omega_{nl} - \omega_r$.</caption><thead><tr><th></th><th>k, +</th><th>p, +</th><th>q, +</th><th>k, −</th><th>p, −</th><th>q, −</th></tr></thead><tbody><tr><th>|ξ|</th><td>0.89</td><td>0.93</td><td>0.52</td><td>0.041</td><td>0.040</td><td>0.0017</td></tr><tr><th>ω<sub>r</sub></th><td>0.099</td><td>−0.088</td><td>−0.020</td><td>−0.099</td><td>0.088</td><td>0.020</td></tr><tr><th>γ</th><td>4.2 × 10<sup>−3</sup></td><td>3.1 × 10<sup>−3</sup></td><td>−1.8 × 10<sup>−4</sup></td><td>−1.8</td><td>−1.8</td><td>−4.8</td></tr><tr><th>ω<sub>nl</sub></th><td>0.20</td><td>−0.19</td><td>−0.016</td><td>0.20</td><td>−0.19</td><td>−0.016</td></tr><tr><th>δω</th><td>0.11</td><td>−0.11</td><td>−0.075</td><td>−1.06</td><td>1.12</td><td>23.2</td></tr></tbody></table>
$$ 
\begin{aligned}
\partial_t \chi_k^{s_k} & - \gamma_k^{s_k} \chi_k^{s_k} \\
& = \sum_{\sigma_p, \sigma_q} m_{kpq}^{s_k \sigma_p \sigma_q} \cos \left( \delta_{kpq}^{s_k \sigma_p \sigma_q} - \psi_{kpq}^{s_k \sigma_p \sigma_q} \right) \chi_p^{\sigma_p} \chi_q^{\sigma_q}
\end{aligned} 
\quad (28) $$

which contain the phases only through their sums (i.e. $\psi$ variables). We can also write an equation for the $\psi_{kpq}^{s_k s_p s_q}$ explicitly as:

$$ 
\begin{aligned}
\partial_t \psi_{kpq}^{s_k s_p s_q} &+ (\omega_k^{s_k} + \omega_p^{s_p} + \omega_q^{s_q}) \\
&= \sum_{\sigma_p, \sigma_q} m_{kpq}^{s_k \sigma_p \sigma_q} \sin \left( \delta_{kpq}^{s_k \sigma_p \sigma_q} - \psi_{kpq}^{s_k \sigma_p \sigma_q} \right) \frac{\chi_p^{\sigma_p} \chi_q^{\sigma_q}}{\chi_k^{s_k}} \\
&+ \sum_{\sigma_q, \sigma_k} m_{pqk}^{s_p \sigma_q \sigma_k} \sin \left( \delta_{pqk}^{s_p \sigma_q \sigma_k} - \psi_{pqk}^{s_p \sigma_q \sigma_k} \right) \frac{\chi_q^{\sigma_q} \chi_k^{\sigma_k}}{\chi_p^{s_p}} \\
&+ \sum_{\sigma_k, \sigma_p} m_{qkp}^{s_q \sigma_k \sigma_p} \sin \left( \delta_{qkp}^{s_q \sigma_k \sigma_p} - \psi_{qkp}^{s_q \sigma_k \sigma_p} \right) \frac{\chi_k^{\sigma_k} \chi_p^{\sigma_p}}{\chi_q^{s_q}}
\end{aligned} 
\quad (29) $$

While the form of (29) looks terribly complicated (e.g. when we expand the sums we have 8 equations, each of whom having 12 terms on their right hand side) it is useful for insight into phase locking. For example by setting $\partial_t \psi_{kpq}^{s_k s_p s_q} = 0$ in (29), and $\partial_t \chi_k^{s_k} = 0$ in (28), we can obtain constant amplitude, phase locked solutions, if such solutions exist. Unfortunately, even the computation of this “fixed point” requires numerical analysis. We can also integrate (28-29) numerically, which gives exactly the same result as the system in terms of $\xi_k^{\pm}$.

### III. INTERACTIONS WITH ZONAL FLOWS

When two non-zonal modes interact with a zonal one the evolution equations and the nonlinear interaction coefficients are different from non-zonal three wave interactions discussed in the previous section. Using the original variables $\Phi_k$ and $n_k$ as in (3-4), zonal and non-zonal modes interact with the same nonlinear interaction coefficients but different linear propagators. However, when we diagonalize the linear propagator, the nonlinear interaction coefficients for zonal and non-zonal modes differentiate.

In particular we have

$$ M_{kpq}^{\phi s_p s_q} = - \frac{\hat{\mathbf{z}} \times \mathbf{p} \cdot \mathbf{q} (q^2 - p^2) C^2}{4 \Omega_p^{s_p *} \Omega_q^{s_q *} k^2 p^2 q^2} \quad (30) $$

$$ M_{kpq}^{ns_p s_q} = \frac{\hat{\mathbf{z}} \times \mathbf{p} \cdot \mathbf{q} C}{4 \Omega_p^{s_p *} \Omega_q^{s_q *} p^2 q^2} \left[ (B_q - i \Omega_q^{s_q *}) q^2 - (B_p - i \Omega_p^{s_p *}) p^2 \right] \quad (31) $$

$$ M_{kpq}^{s_k \phi s_q} = i \frac{\hat{\mathbf{z}} \times \mathbf{p} \cdot \mathbf{q}}{2 \Omega_q^{s_q *} q^2} \left[ (B_q - i \Omega_q^{s_q *}) q^2 - (B_k - i \Omega_k^{s_k}) (q^2 - p^2) \right] \quad (32) $$

$$ M_{kpq}^{s_k n s_q} = i \frac{\hat{\mathbf{z}} \times \mathbf{p} \cdot \mathbf{q} C}{2 \Omega_q^{s_q *} q^2} \quad (33) $$

so that for three waves $k, p$ and $q$ with $q_y = 0$, we can write:

$$ \partial_t \bar{\Phi}_q + \nu_{ZF} \bar{\Phi}_q = \sum_{s_k, s_p} M_{qkp}^{\phi s_k s_p} \xi_k^{s_k *} \xi_p^{s_p *} \quad (34) $$

$$ \partial_t \bar{n}_q + D_{ZF} \bar{n}_q = \sum_{s_k, s_p} M_{qkp}^{n s_k s_p} \xi_k^{s_k *} \xi_p^{s_p *} \quad (35) $$

$$ \partial_t \xi_k^{s_k} + i \omega_k^{s_k} \xi_k^{s_k} = \sum_{s_p} M_{kpq}^{s_k s_p \phi} \xi_p^{s_p *} \Phi_q^* + M_{kpq}^{s_k s_p n} \xi_p^{s_p *} n_q^* \quad (36) $$

$$ \partial_t \xi_p^{s_p} + i \omega_p^{s_p} \xi_p^{s_p} = \sum_{s_p} M_{pqk}^{s_p \phi s_k} \Phi_q^* \xi_k^{s_k *} + M_{pqk}^{s_p n s_k} n_q^* \xi_k^{s_k *} \quad (37) $$

We can write these in the form (23-25) by letting $\xi_q^+ = \Phi_q$ and $\xi_q^- = n_q$ and paying attention to the form of the interaction coefficient $M_{\xi_{kpq}}^{s_k s_p s_q}$ when one of the legs is zonal.

In order to study the interactions between two modes with a zonal flow in the Hasegawa-Wakatani system numerically, we pick a primary wave-number $\mathbf{k} = (0, 1.125)$ which is the linearly most unstable mode on a grid with $dk_x = dk_y = 0.125$ for the $C = 1.0$ case. We choose $\mathbf{p} = (-0.5, -1.125)$ so that $\mathbf{q} = (0.5, 0)$ is a zonal wave number. The 6 field variables are now $\xi_k^{\pm}$, $\xi_p^{\pm}$, $\Phi_q$ and $n_q$ whose evolutions are shown in figure 7 for the case $C = 1$, $\nu_Z = D_Z = 0$ and $\gamma_k \gtrsim \gamma_p > 0$. In the final state, the system finds a fixed point characterized by constant nonlinear frequency shifts, constant amplitudes and constant $\psi_{kpq}$'s. However this kind of steady state solution seems to be exclusive to the single triad case.

#### A. Triad pairs

Because of the symmetry of the system, if we consider two wave-numbers $\mathbf{p}_1 = -\mathbf{k} - \mathbf{q}$ and $\mathbf{p}_2 = -\mathbf{k} + \mathbf{q}$ with $\mathbf{k}$ in $\hat{\mathbf{y}}$ and $\mathbf{q}$ in $\hat{\mathbf{x}}$ directions, we get two triads that are reflections of one another with respect to the axis defined by $\mathbf{k}$. Such a system involves four different wave-numbers connected with two different triads. Including the $p \leftrightarrow q$ transformation we have four triads as shown in figure 8. However as long as we use symmetric forms for the interaction coefficients, we can drop the two triads we obtain
<figure><img src="image_8.png" /><figcaption>Figure 7. Time evolution of the three wave equations involving a zonal mode $q$, for the case $C = 1$, $\nu_Z = D_Z = 0$ and $\gamma_k \gtrsim \gamma_p > 0$ with $k_y = 1.125$ and $q = 0.5$ [i.e. $\mathbf{k} = (0, k_y)$, $\mathbf{p} = (-q, -k_y)$ and $\mathbf{q} = (q, 0)$]. The system reaches a steady state by introducing nonlinear frequencies in order to arrive at a state where the sums of phases $\psi_{kpq}$'s are constant. Note that it is $p$, which becomes the dominant mode in the final state and the existence of zonal flows does not lead to a complete suppression of turbulence. Instead the zonal flow acquires a constant nonlinear frequency.</figcaption></figure>

from the $p \rightleftharpoons q$ transformation and count only two triads. Since the two triads of such a pair are reflections of one another, the nonlinear interaction coefficients differ only in sign while the complex frequencies are the same, and as there are two eigenmodes for each wave-number, we have 8 equations. The equations for zonal modes can be written from (34-35) as:

$$ \partial_t \Phi_q + \nu_Z \Phi_q = \sum_{s_k, s_p} M_{qkp_1}^{\phi s_k s_p} \left( \xi_k^{s_k*} \xi_{p_1}^{s_{p_1}*} - \xi_k^{s_k} \xi_{p_2}^{s_{p_2}} \right) \quad (38) $$

$$ \partial_t n_q + D_Z n_q = \sum_{s_k, s_p} M_{qkp}^{ns_k s_p} \left( \xi_k^{s_k*} \xi_{p_1}^{s_{p_1}*} - \xi_k^{s_k} \xi_{p_2}^{s_{p_2}} \right), \quad (39) $$

which is possible since $M_{\xi_k p_2 q}^{s_k s_p \{n, \phi\}} = -M_{\xi_k p_1 q}^{s_k s_p \{n, \phi\}}$ because $p_2^2 = p_1^2$ and $p_{2y} = p_{1y}$ while $p_{2x} = -p_{2x}$. The equation for the primary mode, can be written as:

$$ \partial_t \xi_k^{s_k} + i \omega_k^{s_k} \xi_k^{s_k} = \sum_{s_p} \left[ M_{\xi_k p q}^{s_k s_p \phi} (\Phi_q^* \xi_{p_1}^{s_{p_1}*} + \Phi_q \xi_{p_2}^{s_{p_2}*}) + M_{\xi_k p q}^{s_k s_p n} (n_q^* \xi_{p_1}^{s_{p_1}*} + n_q \xi_{p_2}^{s_{p_2}*}) \right], \quad (40) $$

<figure><img src="image_9.png" /><figcaption>Figure 8. All the four triads involved in the interaction between the most unstable mode with $\mathbf{k} = k_y \hat{\mathbf{y}}$ with $k_y = 1.125$ and a given zonal mode with $q = 1.0$, obtained by reflection with respect to $\mathbf{k}$ and the exchange of $p$ and $q$ of the primary triad, which is shaded. The existence of the reflected triad is indeed important as it changes the qualitative behavior with respect to the single triad case.</figcaption></figure>

and the remaining two equations are the same as (37) but with different signs and conjugations:

$$ \partial_t \xi_{p_1}^{s_p} + i \omega_{p_1}^{s_p} \xi_{p_1}^{s_p} = \sum_{s_k} \left( M_{\xi_{p_1 k q}}^{s_p s_k \phi} \Phi_q^* + M_{\xi_{p_1 k q}}^{s_p s_k n} n_q^* \right) \xi_k^{s_k*} \quad (41) $$

$$ \partial_t \xi_{p_2}^{s_p} + i \omega_{p_2}^{s_p} \xi_{p_2}^{s_p} = - \sum_{s_k} \left( M_{\xi_{p_1 k q}}^{s_p s_k \phi} \Phi_q + M_{\xi_{p_1 k q}}^{s_p s_k n} n_q \right) \xi_k^{s_k*} \quad (42) $$

where $\omega_{p_2}^{s_p} = \omega_{p_1}^{s_p}$. Notice that this is also equivalent to one of the radial Fourier modes of a quasi-linear (e.g. zonostrophic) interaction, where for each field one would consider a single $p_y$ but the full spatial dependence in $x$.

The results of the system (38-42) are shown in figure 9 for the $C = 1$ case with $k_y = 1.125$ [i.e. $\mathbf{k} = (0, k_y)$, $\mathbf{p}_1 = (-q, -k_y)$, $\mathbf{p}_2 = (q, -k_y)$ and $\mathbf{q} = (q, 0)$] for $q = (1.0, 1.5, 2.0, 4.0)$ from top to bottom respectively. For $q \le k_y$ we have instability and $p$ keeps growing exponentially whereas for $q > k_y$ we get some sort of steady or limit cycle state. Performing a scan of $k_y$ and $q$ for this two triad system (keeping in mind that for $k_y > 2$ we have no instability and therefore the pump mode decays) we observe that we can define a four wave interaction condition of the form $\omega_{kr}^{s_k} + \omega_{p1r}^{s_{p1}} + \omega_{p2r}^{s_{p2}} + \omega_{qr}^{s_q} = 0$, which turns into $\Omega_k^{s_k} + 2\Omega_p^{s_p} = 0$ since $\omega_{qr} = 0$ and $\omega_{p1r} = \omega_{p2r} = \omega_{pr}$. There seems to be 3 distinct regions in figure X: for $q < 1$, the $\xi_p^+$ modes grow exponentially as in the top plot of figure 9, for the central region where $q \approx 1$, we have saturation and then somewhat chaotic evolution, and finally for $q \gg 1$, we observe limit cycle oscillations between $\xi_k^+$
<figure><img src="image_10.png" /><figcaption>Figure 9. Evolution of a triad pair with the same parameters as figure 7, no zonal flow damping $\nu_Z = D_Z = 0$ and $k_y = 1.125$ [i.e. $\mathbf{k} = (0, k_y)$, $\mathbf{p}_1 = (-q, -k_y)$, $\mathbf{p}_2 = (q, -k_y)$ and $\mathbf{q} = (q, 0)$] for four different values of $q = (1.0, 1.5, 2.0, 4.0)$ from top to bottom for which the growth rates of the subdominant modes are $\gamma_p = (0.00099, -0.0016, -0.0042, -0.017)$ respectively. Note that apart from the second plot, which displays some chaotic behavior, the curves for $\xi_{p1}^+$ and $\xi_{p2}^+$ overlap almost exactly.</figcaption></figure>

and $\xi_p^+$ modes, mediated by zonal flows.

One is tempted to argue that since the $p$ with $p_x < p_y$ wins the competition to attract more energy, the cascade will proceed in this direction, and in the next step we can consider the interaction of this $\xi_p^+$ as the pump mode for the next triad etc. However, since each mode interacts with many triads simultaneously, the fact that $\xi_p^+$ wins the competition in the single triad (or one triad and its reflection) configuration does not really mean the energy will indeed go this way.

## B. Triad Networks

In order to study the fate of the cascade, we need to consider multiple triads that are connected to one another. However as we add more zonal and non-zonal modes, it becomes quite complicated to keep track of all the interaction coefficients, conjugations etc. In order to simplify this task, we can divide the problem into two steps i) construction of a network of three body interactions and ii) computation of the evolution of the field variables on this network. For example for the above problem we need to consider a network of $N_k = 4$ wave

<figure><img src="image_11.png" /><figcaption>Figure 10. The structure of the network with a single $k_y$ with $k_y = 1.125$ shown as a filled (red if in color) node. A reduced version with $q$ values that only go up to 0.5 is shown for clarity. Notice that in this network while all of the 26 triads involve one of the zonal modes, only 8 of them involve the $q = 0$ mode.</figcaption></figure>

number nodes, coupled to $N_t = 2$ triads, with $N_f = 2$ fields in each node, with an interaction coefficient of the size $N_f \times N_f \times N_f$ for each connection. Since a network in Fourier space is made up of three body interactions, for each node, we can compute a list of interacting pairs and the interaction coefficients, so that we can write

$$ \partial_t \xi_{\ell}^{i} + L_{\ell}^{ij} \xi_{\ell}^{j} = \frac{1}{2N} \sum_{\ell', \ell''=i_{\ell}} M_{\ell\ell'\ell''}^{ijk} (\xi_{\ell'}^{j})^{c_{\ell'}} (\xi_{\ell''}^{k})^{c_{\ell''}} \quad (43) $$

where $i_\ell$ is the list of precomputed interaction pairs for the node $\ell$. The indices $i, j$ and $k$ correspond to different fields (eigenmodes or $\Phi_k$ and $n_k$), the matrix $L_\ell^{ij}$ is the linear matrix in $k$ space (i.e. diagonal with the elements $i\omega_\ell^\pm$ for the eigenmodes), the $M_{\ell\ell'\ell''}^{ijk}$ is the interaction coefficient for each interaction and $N$ is the number of independent wave number nodes so that when we reach the full grid, we have exactly the same interaction coefficients as the system formulated using discrete fast Fourier transforms (i.e. divided by $N_x \times N_y$). Finally if we write the triad interaction condition in the form $\mathbf{k}_\ell + \sigma_{\ell'} \mathbf{k}_{\ell'} + \sigma_{\ell''} \mathbf{k}_{\ell''} = 0$, where $\sigma$ are $\pm 1$, the $(\xi_{\ell'}^j)^{c_{\ell'}}$ are defined as:

$$ (\xi_{\ell'}^j)^{c_{\ell'}} = \begin{cases} \xi_{\ell'}^j & \sigma_{\ell'} = -1 \\ \xi_{\ell'}^{j*} & \sigma_{\ell'} = +1 \end{cases} $$

This is necessary unless we have the negative of each wave number vector as a separate node in the network.

Notice that when computing the nonlinear interaction coefficients for the eigenmodes, we would use (21) if all the nodes have nonzero $k_y$. In contrast we would use (30)
and (31) if the receiving node (i.e. node $\ell$) is zonal or (32) and (33) if one of the interacting pairs (i.e. $\ell'$ or $\ell''$) are zonal. Two or more zonal mode do not interact because of the geometric factor $\hat{\mathbf{z}} \times \mathbf{p} \cdot \mathbf{q}$, which appear in front of all the interaction coefficients.

Finally, if it makes sense to zero out some of the fields at a given wave-number (e.g. in eigenmode formulation, we may decide to throw away some damped modes), one may switch to a formulation where each node corresponds to a wave-number/field variable combination via $\{k_x, k_y, s_k\} \to \ell$. In this case, assuming that the linear matrix $L_{\ell}^{ij}$ in (43) diagonal takes the form:

$$ \partial_t \xi_{\ell} + i\omega_{\ell}\xi_{\ell} = \frac{1}{N} \sum_{\ell', \ell''=i_{\ell}} M_{\ell\ell'\ell''} \xi_{\ell'}^{c_{\ell'}} \xi_{\ell''}^{c_{\ell''}} \quad (44) $$

### C. Order Parameters

The phases of wave-number nodes in Hasegawa-Wakatani turbulence evolve according to (19) or written explicitly as (22). This suggests that one can possibly define some kind of order parameter for this system. The usual definition of the Kuramoto order parameter can be written for the network formulation of (44) as:

$$ z = re^{i\psi} = \frac{1}{N} \sum_{\ell} e^{i\varphi_{\ell}} \quad (45) $$

without explicitly distinguishing + or – modes. However this order parameter based on an unweighted sum is probably relevant only if all the oscillators were identical with all-to-all, unweighted couplings of the Kuramoto type. Instead we can use an amplitude filtered Kuramoto order parameter (i.e. the sum is computed only over the oscillators with an amplitude larger than a threshold), or define a weighted version of (45) as:

$$ z = re^{i\psi} = \frac{\sum_{\ell} \chi_{\ell} e^{i\varphi_{\ell}}}{\sum_{\ell} \chi_{\ell}} \quad (46) $$

whose absolute value would tend towards 1 if the relevant phases (i.e. those that have large amplitude) are the same. However note that the weighted order parameter tends towards 1 also when one of the modes dominate over the others, while $\psi$ as defined in (46), can still be used as a mean phase.

It would also make sense to look at the net effect on the nonlinear term on the phases instead. As discussed in Section IB, since we can write:

$$ \partial_t \varphi_{\ell} = -\omega_{\ell} + \frac{1}{N\chi_{\ell}} \mathrm{Im} \left[ \sum_{\ell', \ell''=i_{\ell}} M_{\ell\ell'\ell''} \xi_{\ell'}^{c_{\ell'}} \xi_{\ell''}^{c_{\ell''}} e^{-i\varphi_{\ell}} \right] \quad (47) $$

<figure><img src="image_12.png" /><figcaption>Figure 11. Time evolution for a number of triad pairs (as defined in section III A) with different values of $q$ in the network of interacting triads for $C = 1$ case with $\nu_{ZF} = D_{ZF} = 10^{-3}$. A steady state turbulence level is observed, with elevated levels of zonal flows at large scales.</figcaption></figure>

for the evolution of the phase, we can define:

$$ Z_{\ell} = R_{\ell} e^{i\psi_{\ell}} = \frac{1}{N\chi_{\ell}} \left( \sum_{\ell', \ell''=i_{\ell}} M_{\ell\ell'\ell''} \xi_{\ell'}^{c_{\ell'}} \xi_{\ell''}^{c_{\ell''}} \right) \quad (48) $$

with $d_{\ell}$ being the number of interactions for the node $\ell$ (i.e. length of $i_{\ell}$), as some kind of local order parameter for the node $\ell$, allowing us to write the phase equation as:

$$ \partial_t \varphi_{\ell} = -\omega_{\ell} + R_{\ell} \sin(\psi_{\ell} - \varphi_{\ell}), \quad (49) $$

which attracts the system towards $\varphi_{\ell} = \psi_{\ell} + 2n\pi$.

### D. Specific network configurations

#### 1. Network with a single $k_y$:

We consider a network of triad pairs as discussed in section III A, with a single value of $k_y$ and $q$ values that go from 0.125 to 4.0 in steps of 0.125. Notice that such a network has many different types of interactions as shown in figure 10, but all of those involve one of the zonal modes, which means that if we compute the inverse Fourier transform in the $x$ direction, the network can be seen to be equivalent to the single $k_y$, full-$x$, quasi-linear model [38, 39], since in both cases we have full spatial evolution but only nonlinear coupling is with the zonal flow.

For the case $C = 1$, without zonal flow damping (not shown) we observe that the zonal flows dominate and all the other modes decay to zero. This may well be what happens also in direct numerical simulations (DNS) eventually: what we observe in numerical simulations without zonal flow damping is a continual increase of zonal flows even for very long simulations.

In contrast, when we introduce zonal flow damping by letting $\nu_{ZF} = D_{ZF} = 10^{-3}$, we get dynamics and k-
<figure><img src="image_13.png" /></figure>

Figure 12. The top plot shows the order parameter $r$ defined in (45) or (46) as a function of time for a network with single $q$ and multiple $k_y$. The two definitions are in reasonable agreement apart from the peak around $t = 2500$ for the weighted order parameter, which corresponds to the linear growth phase, where only a few modes around the most unstable mode dominate. This can be seen at the bottom plot where the amplitudes of a triad pair with $q = 0.5$ and $k_y = 1.125$ are shown. Around $t = 2500$ the blue curve clearly dominates.

<figure><img src="image_14.png" /></figure>

Figure 14. The top plot shows the order parameter $r$ defined in (45) or (46) as a function of time for a network with single $q$ and multiple $k_y$. The two definitions are in reasonable agreement apart from the peak around $t = 2500$ for the weighted order parameter, which corresponds to the linear growth phase, where only a few modes around the most unstable mode dominate. This can be seen at the bottom plot where the amplitudes of a triad pair with $q = 0.5$ and $k_y = 1.125$ are shown. Around $t = 2500$ the blue curve clearly dominates.

<figure><img src="image_15.png" /></figure>

Figure 13. The structure of the network with a single $q = 0.5$ zonal mode, shown as a filled (red if in color) node. A reduced version with $k_y$ values that only go up to 0.5 is shown for clarity. Only 8 of the full 26 triads involve the zonal flow.

spectra which look more like fully developed Hasegawa-Wakatani turbulence, as shown in figure 11, with high levels of zonal flows at large scales.

## 2. Network with a single $q$:

Here, we consider a network of triad pairs with a single $q$, and a grid of values of $k_y$ going from 0.125 to 4.0 in steps of 0.125. A reduced version of such a network is shown in figure 13. Physically this network corresponds to the opposite case where we consider a single $q$ with the whole $y$ dynamics if we compute the inverse Fourier transform in $y$. Since it involves bunch of oscillators with different frequencies (as $\omega$ is mostly a function of $k_y$) that are coupled to each other and to a zonal mode that may play the role of a dominant mean field, it has the basic ingredients that may lead to synchronization.

Nonetheless numerical observations suggest that there is no obvious route to global synchronization in the three body network of interacting triads consisting of a zonal mode and drift waves of different $k_y$ either. The weighted order parameter shows a brief increase during the nonlinear saturation phase as the energy is transferred to the zonal flow, but otherwise remain close to zero, while the Kuramoto order parameter simply remains close to zero the whole time as can be seen in figure 14. Since we observed no qualitative difference between the runs with or without zonal flow damping for this case, we only show those with $\nu_{ZF} = D_{ZF} = 10^{-3}$.
<figure><img src="image_16.png" /></figure>

Figure 15. The top plot shows the order parameter $r$ defined in (45) or (46) as a function of time for a DNS. The bottom plot shows the amplitudes of a triad pair with $q = 0.5$ and $k_y = 1.125$ in order to compare with the earlier plots. The saturation levels for the amplitudes are different because of the normalization factor $N_x^{-1}N_y^{-1}$ in front of the nonlinear term implied in discrete Fourier transforms.

E. Direct numerical simulations

One can think of direct numerical simulation (DNS) on a regular rectangular grid as a “network” in Fourier space, in the sense that it consists of a collection of wave number nodes connected to each other through triadic interactions. In contrast to the networks that we considered that contain a single zonal mode, or a single $q = 0$ mode, a regular rectangular grid has all the possible wave-numbers in a particular range, and it allows using more efficient methods for computing the convolution sums. In practice, the high resolution direct numerical simulations that we discuss here were performed with a standard pseudo-spectral solver (i.e. with periodic boundary conditions in both directions) using 2/3 rule for dealiasing and adaptive time stepping.

As with all the previous examples of single or multiple triads, or networks with a particular selection of nodes and triads, we use $C = 1$, $\kappa = 0.2$. Since we have a larger range of wave-numbers, we choose $\nu = D = 10^{-4}$, with a box size of $L_x = L_y = 16\pi$ and a padded resolution of $1024 \times 1024$. The results show (see figures 15 and 16):

i. Initial linear growth followed by nonlinear saturation.

ii. Formation and finally suppression of nonlinear of convective cells that transfer vorticity radially.

iii. Consequent stratification of vorticity leading to a state dominated by zonal flows (as in figure 16).

iv. Coherent nonlinear structures (e.g. vortices) that are advected by the zonal flows in regions of weak zonal shear, get sheared apart if they fall into a region of strong zonal shear.

<figure><img src="image_17.png" /></figure>

Figure 16. Snapshots of vorticity and density at $t=5000$ from DNS. The blue curve in both plots shows the zonal velocity whose values are given on the right hand axes. An example coherent vortex, that was moving upwards is encircled.

Since the wave-like dynamics seems to be primarily in $y$ direction and reasonably localized in $x$, we can compute the Fourier transform in $y$ and plot phase of $\xi_{k_y}^{\pm} = \chi_{k_y}^{\pm} e^{i\phi_{k_y}^{\pm}}$ at each $x$, compute $\partial_t \phi_{k_y}^{\pm} (x,t)$ in order to compute the phase speeds (see figure 17). We can also compute an order parameter as a function of $x$ and $t$ from this data.

While it is clear from 15 that there is no global synchronization in direct numerical simulations, the plateau form of the phase velocity as a function of $k_y$ at the radii where it is positive for large scales, suggest that a process of phase locking similar to soliton formation in nonlinear Schrödinger equation, where nonlinearity would balance dispersion is at play for a range of $k_y$ values around the linearly unstable mode. While $\omega/k_y$ being the same across a range of $x$ and $k_y$ values is obviously very different from $\omega$ being the same. However if we note that the nonlinear dispersion relation takes the form $\omega(x, k_y) = \bar{v}_\phi(x) k_y$, at the lowest order we can see that the frequency in the frame moving with the zonal flow velocity becomes zero. This is roughly consistent with what we see in time evolution, where coherent structures like rotating vortices are advected by zonal flows. In order for such a detailed structure

IV. CONCLUSION

A detailed analysis of triadic interactions formulated in terms natural frequencies reveals the complex nature of the dynamics of the phases and amplitudes in the Hasegawa Wakatani system. In particular, it is observed that a single resonant (or near resonant) triad, including a pump mode and two other modes, can saturate by adjusting the sums of phases of its legs ($\psi_{kpq}^{s_k s_p s_q} = \phi_k^{s_k} + \phi_p^{s_p} + \phi_q^{s_q}$) to be asymptotically constant, resulting in a set of nonlinearly shifted frequencies and constant amplitudes. When the interactions with zonal flows are considered, a similar saturation is possible for a single
triad even without the condition of resonance. However this solution breaks down when we add the triad, which is the reflection of the original one with respect to the y axis (or the wave-vector **k** vector). Instead we observe three different behavior for these triad pairs as a function of the radial wave number.

i. For smaller radial wave numbers, we find that the subdominant mode becomes the dominant one and grows exponentially. We call those unstable triads. They are associated with unstable subdominant modes.

ii. For medium radial wave numbers, after an initial growth phase, the system saturates with a more or less chaotic evolution, where the energy goes back and forth between the modes. We call these saturated triads. They are associated with weakly unstable, or weakly damped subdominant modes.

iii. For large radial wave numbers the system decays to a steady state solution after a number of limit cycle oscillations. In some cases, these limit cycle oscillations can continue until the end of the simulation time. We call these decaying triads (even though they don't decay to zero but to a constant). They are associated with strongly damped subdominant modes.

In order to study the dynamics when those triads are connected to one another, we considered a network formulation where the wave numbers (or wave number eigen-mode combinations) are considered as nodes, and each triad represents a three body interaction. It is shown while the zonal flow is almost never dominant in a single triad, when the whole triad network with a large number of triads is considered, the zonal modes become dominant almost in each triad. Thus, the system can reach a steady state where the zonal flow dominates as the other modes decay.

In terms of triadic interactions, as the zonal flow becomes dominant, it plays the role of a collective mean field, in the sense that for each mode individual interactions with non-zonal modes start to become less important compared to the interaction with the zonal flow. This happens only when the number of triads is large enough so that the collective wins over the individual. It is interesting to note that this picture is qualitatively consistent with that of inhomogeneous wave-kinetic formulation, where the zonal flow is treated as a collective mean field, and the direct interaction between the modes are either dropped or modeled with a diffusion operator. This suggests that the wave-kinetic formulation may hold beyond its range of validity.

Playing with the range of radial wave-numbers of the network model, we observe that when the range includes only unstable triads [i.e. (i) above], or unstable and sat-

<figure><img src="image_18.png" /><figcaption>Figure 17. Profiles of phase velocity as a function of $k_y$, at three different values of $x$ (i.e. 14.74, 29.48 and 40.54) averaged over $t = [4500, 5000]$ shown at the top plot. The three plots that follow show the detailed time evolution (on the left y axes) of phase velocity as a function of $x$ for three different values of $k_y$ (i.e. 1, 3 and 5), together with the mean velocity profile shown for reference (on the right y axes). The phase velocity is computed using $v_\phi = -\partial_t \phi_{k_y}^+(x,t)/k_y$. The $k_y$'s for which the time evolution is given and the $x$'s for which the phase velocities are shown are marked with horizontal lines in the corresponding figures.</figcaption></figure>

urated triads [i.e. (i) and (ii) above] the network system remains unstable. It saturates only when we include a sufficient range of decaying triads, with subdominant modes with $\gamma_p^+ < 0$. This means that 'local coupling to damped modes' (i.e. $\gamma_p^-$ modes even though $\gamma_p^+ > 0$) is not a real mechanism for turbulent saturation. However since the fact that $\gamma_p^+ < 0$ for those modes do not come directly from dissipation but rather the detailed form of the linear growth/damping whose form is determined by various parameters including dissipation, it is correct to argue that in contrast to the Kolmogorov picture where there is an injection scale, a dissipation scale and the inertial range in between, plasma turbulence can generate and dissipate energy in much closer scales, even though one may observe clear power law scalings.

One of the goals of the current paper was to study the effect of nonlinear synchronization of drift waves[40] on the turbulent cascade using a framework similar to the Kuramoto model[41], which has already been attempted using simple models in fusion plasmas[42, 43]. We hoped by considering a network of connected triads interacting with zonal flows we could setup a system that would tend
toward synchronization through slight nonlinear modifications of the frequencies through their interactions with the zonal flow, playing the role of the control parameter. However due to particular form of the systematic dependency of the frequencies to the wave-numbers through the dispersion relation, such a system does not seem to tend towards synchronization. It should be checked whether or not the discretization resulting from boundary conditions, for example in cylindrical geometry change this picture drastically by impeding resonant interactions[44, 45] especially among large scale modes.

[1] A. Hasegawa and M. Wakatani, Pys. Rev. Lett 50, 682 (1983).

[2] M. Lesieur and J. Herring, Journal of Fluid Mechanics 161, 77 (1985).

[3] F. H. Busse and K. E. Heikes, Science 208, 173 (1980).

[4] L. K. Currie and S. M. Tobias, Physics of Fluids 28, 017101 (2016).

[5] A. Hasegawa and K. Mima, Phys. Fluids 21, 87 (1978).

[6] C. Connaughton, S. Nazarenko, and B. Quinn, Physics Reports 604, 1 (2015).

[7] P. H. Diamond, S.-I. Itoh, K. Itoh, and T. S. Hahm, PPCF 47, R35 (2005).

[8] C. Holland, G. R. Tynan, J. H. Y. A. James, D. Nishijima, M. Shimada, and N. Taheri, Plasma Physics and Controlled Fusion 49, A109 (2007).

[9] A. V. Pushkarev, W. J. T. Bos, and S. V. Nazarenko, Physics of Plasmas 20, 042304 (2013).

[10] Y. Zhang and S. I. Krasheninnikov, Physics of Plasmas 27, 122303 (2020).

[11] B. D. Scott, Journal of Computational Physics 78, 114 (1988).

[12] A. E. Koniges, J. A. Crotinger, and P. H. Diamond, Physics of Fluids B: Plasma Physics 4, 2785 (1992).

[13] B. Friedman and T. A. Carter, Physics of Plasmas 22, 012307 (2015).

[14] D. D. Sarto and A. Ghizzo, Fluids 2 (2017), 10.3390/fluids2040065.

[15] W. Bos, B. Kadoch, S. Neffaa, and K. Schneider, Physica D: Nonlinear Phenomena 239, 1269 (2010), at the boundaries of nonlinear physics, fluid mechanics and turbulence: where do we stand? Special issue in celebration of the 60th birthday of K.R. Sreenivasan.

[16] J. Anderson and B. Hnat, Physics of Plasmas 24, 062301 (2017).

[17] F. Y. Gang, P. H. Diamond, J. A. Crotinger, and A. E. Koniges, Physics of Fluids B: Plasma Physics (1989-1993) 3, 955 (1991).

[18] G. Hu, J. A. Krommes, and J. C. Bowman, Physics of Plasmas (1994-present) 4, 2116 (1997).

[19] R. Singh and P. H. Diamond, Plasma Physics and Controlled Fusion 63, 035015 (2021).

[20] I. R. Goumiri, C. W. Rowley, Z. Ma, D. A. Gates, J. A. Krommes, and J. B. Parker, Physics of Plasmas 20, 042501 (2013).

[21] J. Anderson, E.-j. Kim, B. Hnat, and T. Rafiq, Physics of Plasmas 27, 022307 (2020).

[22] R. A. Heinonen and P. H. Diamond, Phys. Rev. E 101, 061201 (2020).

[23] N. Kasuya, M. Yagi, M. Azumi, K. Itoh, and S.-I. Itoh, Journal of the Physical Society of Japan 76, 044501 (2007).

[24] P. Vaezi, C. Holland, S. C. Thakur, and G. R. Tynan, Physics of Plasmas 24, 092310 (2017).

[25] P. Donnel, P. Morel, C. Honoré, Ö. Gürcan, V. Pisarev, C. Metzger, and P. Hennequin, Physics of Plasmas 25, 062127 (2018).

[26] A. I. Smolyakov and P. H. Diamond, Phys. Plasmas 6, 4410 (1999).

[27] L. Chen, Z. Lin, and R. White, Physics of Plasmas 7, 3129 (2000).

[28] S. Champeaux and P. Diamond, Physics Letters A 288, 214 (2001).

[29] P. Manz, M. Ramisch, and U. Stroth, Phys. Rev. Lett. 103, 165004 (2009).

[30] U. Stroth, P. Manz, and M. Ramisch, Plasma Physics and Controlled Fusion 53, 024006 (2011).

[31] H. Biglari, P. H. Diamond, and P. W. Terry, Phys. Fluids B 2, 1 (1990).

[32] P. W. Terry, Rev. Mod. Phys. 72, 109 (2000).

[33] M. A. Malkov, P. H. Diamond, and M. N. Rosenbluth, Phys. Plasmas 8, 5073 (2001).

[34] E.-j. Kim and P. H. Diamond, Phys. Rev. Lett. 90, 185006 (2003).

[35] Z. B. Guo and P. H. Diamond, Phys. Rev. Lett. 117, 125002 (2016).

[36] T. Stoltzfus-Dueck, B. D. Scott, and J. A. Krommes, Physics of Plasmas 20, 082314 (2013).

[37] M. D. Bustamante and E. Kartashova, Europhysics Letters 85, 34002 (2009).

[38] N. Bian, S. Benkadda, O. E. Garcia, J.-V. Paulsen, and X. Garbet, Physics of Plasmas 10, 1382 (2003).

[39] Y. Sarazin, G. Dif-Pradalier, X. Garbet, P. Ghendrih, A. Berger, C. Gillot, V. Grandgirard, K. Obrejan, R. Varennes, L. Vermare, and T. Cartier-Michaud, Plasma Physics and Controlled Fusion 63, 064007 (2021).

[40] D. Block, A. Piel, C. Schröder, and T. Klinger, Phys. Rev. E 63, 056401 (2001).

[41] Y. Kuramoto, *Chemical Oscillations, Waves, and Turbulence* (Springer-Verlag, New York, 1984).

[42] S. Moradi, J. Anderson, and Ö. D. Gürcan, Physical Review E 92, 062930 (2015).

[43] S. Moradi, B. Teaca, and J. Anderson, AIP Advances 7, 115213 (2017).

[44] E. A. Kartashova, Phys. Rev. Lett. 72, 2013 (1994).

[45] E. Kartashova, Nonlinear Resonance Analysis, by Elena Kartashova, Cambridge, UK: Cambridge University Press, 2010 1 (2010).