2. Fluid models

A two-field fluid model for toroidal ITG turbulence describing the nonlinear temporal evolution of electrostatic potential $\phi _k=\phi (k_x,k_y)$ and pressure $p_k=p(k_x,k_y)$ is used to analytically investigate ZF saturation mechanisms. The model may be expressed as a spatially Fourier transformed pair of coupled nonlinear PDEs of the form

(2.1) \begin{align} \frac {\partial p_k}{\partial t} + ik_y(1+\eta )\phi _k + \chi k_\perp ^4 p_k = - \sum _{\boldsymbol{k}^{\prime }} (\boldsymbol{k}^{\prime } \times \hat {z} \boldsymbol{\cdot }\boldsymbol{k})\phi _{\boldsymbol{k}^{\prime }}p_{\boldsymbol{k}^{\prime \prime }}, \end{align}

(2.2) \begin{align} \big[\delta (k_y)+k_\perp ^2\big]\frac {\partial \phi _k}{\partial t} & + ik_y \phi _k - i k_y \epsilon _n p_k + \nu k_\perp ^2 \phi _k = -\frac {1}{2} \sum _{k^{\prime }} (\boldsymbol{k}^{\prime } \times \hat {z} \boldsymbol{\cdot }\boldsymbol{k}) \big(\boldsymbol{k}_\perp ^{\prime \prime 2} - \boldsymbol{k}_\perp ^{{\prime } 2}\big)\phi _{k^{\prime }}\phi _{k^{\prime \prime }}\nonumber\\[4pt]& +\frac {1}{2} \sum _{k_x^{\prime }}k_y\big[k_x'\phi _{k^{\prime }}\phi _{k^{\prime \prime }}|_{k_y'=0}+(k_x-k_x')\phi _{k^{\prime }}\phi _{k^{\prime \prime }} |_{k_y'=k_y}\big] . \end{align}

Equations (2.1) and (2.2) describe two-dimensional turbulence that is strongly localised along a magnetic field line at the outboard midplane of a tokamak. The $x$ coordinate corresponds to the radial direction while $y$ is binormal. The model is spatially homogeneous because the density gradient is treated as constant and strong ballooning implies a localised mode structure to the outboard midplane (Horton, Choi & Tang Reference Horton, Choi and Tang1981). The parameter $\eta$ is the ratio of density to temperature gradient scale lengths, $\epsilon _n$ is the inverse scale length of magnetic field curvature normalised to the density-gradient scale length, $\delta (k_y)$ is the adiabatic electron response, and $\nu$ and $\chi$ are coefficients of collisional dissipation. Length scales are normalised to the ratio of sound speed $c_s = (T_e/m_i)^{1/2}$ to ion cyclotron frequency $\varOmega _{ci} = qB/m_i$ in the form $\rho _s=c_s/\varOmega _{ci}$ , such that $k_\perp \rho _s \rightarrow k_\perp$ . $T_e$ is the electron temperature, $m_i$ is the ion mass, $q$ is the charge of the ion and $B$ is the on-axis magnetic field strength. For modes with $k_y=0$ , there is no adiabatic electron response, hence,

(2.3) \begin{align} \delta (k_y)= \begin{cases} 0 & k_y =0 \\ 1 & k_y\neq 0 \end{cases}. \end{align}

This system is valid in the limit of low frequency relative to the ion cyclotron frequency, and strong ballooning of the fluctuation eigenfunction about the outboard midplane in a toroidal coordinate system transformed via the ballooning representation is assumed. For the electron dynamics, the system assumes adiabatic electrons. This model has been derived and refined in Horton et al. (Reference Horton, Choi and Tang1981), Horton, Hong & Tang (Reference Horton, Hong and Tang1988) and Holland et al. (Reference Holland, Diamond, Champeaux, Kim, Gurcan, Rosenbluth, Tynan, Crocker, Nevins and Candy2003), investigated mathematically and computationally in considerable depth and shown to make good agreement with features of gyrokinetic simulations (Holland et al. Reference Holland, Diamond, Champeaux, Kim, Gurcan, Rosenbluth, Tynan, Crocker, Nevins and Candy2003; Makwana et al. Reference Makwana, Terry and Kim2012, Reference Makwana, Terry, Pueschel and Hatch2014; Terry et al. Reference Terry, Faber, Hegna, Mirnov, Pueschel and Whelan2018, Reference Terry, Li, Pueschel and Whelan2021). (For more detailed discussions of this system, see Horton et al. (Reference Horton, Choi and Tang1981, Reference Horton, Hong and Tang1988), Holland et al. (Reference Holland, Diamond, Champeaux, Kim, Gurcan, Rosenbluth, Tynan, Crocker, Nevins and Candy2003) and Terry et al. (Reference Terry, Faber, Hegna, Mirnov, Pueschel and Whelan2018).) The second nonlinearity of (2.2) arises from the advection of the adiabatic electron response by the ZF. This nonlinearity was neglected in the aforementioned works; in an eigenmode decomposition its contribution to the nonlinearity is seen to be small for peaked density – of order $\epsilon _n$ relative to the dominant advection-of-ion-pressure nonlinearity. This nonlinearity is retained here for the sake of generality.

It can be shown that the linear eigenmodes of this system at $k_y\neq 0$ are stable and unstable modes with eigenfrequencies $\omega _j$ , whose amplitudes are $\beta _2=\beta _2(\boldsymbol{k})$ and $\beta _1=\beta _1(\boldsymbol{k})$ and whose eigenvector coefficients are $R_2$ and $R_1$ . The amplitudes $\beta _1$ and $\beta _2$ are related to $\phi _k$ and $p_k$ via the transformation

(2.4) \begin{align} \begin{pmatrix} p_k \\[2pt] \phi _k \end{pmatrix} = \begin{pmatrix} R_1(\boldsymbol{k}) & \quad R_2(\boldsymbol{k}) \\[2pt] 1 & \quad 1 \end{pmatrix} \begin{pmatrix} \beta _1(\boldsymbol{k}) \\[2pt] \beta _2(\boldsymbol{k}) \end{pmatrix}\!, \end{align}

in terms of eigenvector coefficients

(2.5) \begin{align} R_j = \frac {-\omega _j\big(1+k_\perp ^2\big)+k_y[1+\epsilon _n(1-\sqrt {8})]-\mathrm{i}\nu k_\perp ^2}{ 2k_y \epsilon _n}. \end{align}

It is fairly commonplace in studies of fluid mechanical systems to decompose governing equations into their respective eigenmode representations (Lamb Reference Lamb1911; Cowling Reference Cowling1941; Elsasser Reference Elsasser1950; Tassoul Reference Tassoul1980; Lien & Müller Reference Lien and Müller1992; Terry et al. Reference Terry, Baver and Gupta2006, Reference Terry, Faber, Hegna, Mirnov, Pueschel and Whelan2018) and interpret the nonlinearity as interactions among the constituent linear eigenmodes (Kumar & Goldreich Reference Kumar and Goldreich1989; Galtier et al. Reference Galtier, Nazarenko, Newell and Pouquet2001; Terry et al. Reference Terry, Baver and Gupta2006, Reference Terry, Faber, Hegna, Mirnov, Pueschel and Whelan2018; Azelis et al. Reference Azelis, Perez and Bourouaine2024; Remond-Tiedrez et al. Reference Remond-Tiedrez, Smith and Stechmann2024; Van Beeck et al. Reference Van Beeck, Van Hoolst, Aerts and Fuller2024).

Transforming a system to be expressed in terms of its linear eigenfunctions and eigenvalues often simplifies the mathematical structure of the governing equations. It is a standard technique (introduced in undergraduate physics texts) that is broad in applicability and can readily be applied to many different systems. Its mathematical implementation is not overly preoccupied with what particular physical processes are occurring. Despite this, it allocates considerable physical intuition with which to characterise nonlinear dynamics using the language of the more familiar and easily comprehensible linear physics of the system. For instance, the eigenmode decomposition proposed by Elsasser (Reference Elsasser1950) readily allows for the nonlinearity of the incompressible magnetohydrodynamic system to be understood as interactions between counter-propagating Alfvénic fluctuations (Galtier et al. Reference Galtier, Nazarenko, Newell and Pouquet2001). This is in contrast with potentially esoteric notions of considering all the ways in which velocity and magnetic fields may interact with themselves and each other.

The eigenmode decomposition used in this work casts the dynamics in terms of fluctuations that transport energy up and down the temperature gradient in the form of stable and unstable modes, respectively. In addition, casting the fluid model in terms of these modes allows for immediate inference and quantification of the rates at which energy is being added to and removed from the system of modes due to its interactions with a mean temperature gradient. In the context of fusion-relevant turbulence, where the central issues are stability and confinement degradation, this decomposition offers a particularly natural and physically transparent framework. Its utility has been tested in comparisons with the singular value decomposition, bases with and without dissipation, bases that evolve with relaxing profiles in decaying turbulence and bases that are intrinsically nonlinear. Its utility and transparency is evident in many situations.

It has been previously shown that the toroidal branch of the ITG instability saturates via ZF-catalysed energy transfer from the unstable mode to the stable mode (Hatch et al. Reference Hatch, Terry, Jenko, Merz and Nevins2011; Makwana et al. Reference Makwana, Terry and Kim2012, Reference Makwana, Terry, Pueschel and Hatch2014; Terry et al. Reference Terry, Faber, Hegna, Mirnov, Pueschel and Whelan2018; Li & Terry Reference Li and Terry2022). This particular mechanism also plays a prominent role in dictating other nonlinear dynamical aspects of these types of systems (Whelan et al. Reference Whelan, Pueschel and Terry2018, Reference Whelan, Pueschel, Terry, Citrin, McKinney, Guttenfelder and Doerk2019; Terry et al. Reference Terry, Li, Pueschel and Whelan2021; Li & Terry Reference Li and Terry2022). Saturation by ZF-catalysed nonlinear energy transfer to stable modes is also observed in turbulence driven by the Kelvin–Helmholtz (Fraser et al. Reference Fraser, Terry, Zweibel, Pueschel and Schroeder2021; Tripathi et al. Reference Tripathi, Fraser, Terry, Zweibel and Pueschel2022a , Reference Tripathi, Fraser, Terry, Zweibel and Pueschelb , Reference Tripathi, Fraser, Terry, Zweibel, Pueschel and Anders2023a , Reference Tripathi, Terry, Fraser, Zweibel and Pueschelb ), Goldreich–Schubert–Fricke (Tripathi et al. Reference Tripathi, Barker, Fraser, Terry and Zweibel2024) and semi-collisional trapped electron mode instabilities (Terry et al. Reference Terry, Gatto and Baver2002). More generally, in the absence of ZFs, saturation by nonlinear energy transfer to stable modes is a common feature of fluid models for instability-driven turbulence (Makwana et al. Reference Makwana, Terry, Kim and Hatch2011). This paradigm also arises in both trapped electron mode turbulence (Baver et al. Reference Baver, Terry, Gatto and Fernandez2002; Terry et al. Reference Terry, Baver and Gupta2006) and electron temperature gradient-driven turbulence (Kim & Terry Reference Kim and Terry2010).

In the present analysis, we are interested in the wavenumber range where $k_\perp ^2 \ll 1$ , and use this to further assume that $\chi k_\perp ^2 \ll \nu$ (Terry et al. Reference Terry, Li, Pueschel and Whelan2021). The limit of small $k_\perp$ is taken as a simplifying approximation that is valid because the instability is driven and saturated by stable modes in the regime $k_\perp ^2 \ll 1$ . However, mode dispersion, which enters from the ion polarisation drift, is necessary because otherwise the nonlinear triplet interaction time becomes infinite under the resonant vanishing of the frequency mismatch for non-dispersive waves. This means that the lowest-order contribution to the frequency mismatch for zero collisionality is order $k_\perp ^2$ . This limit is taken consistently, as shown by the ordering scheme given in Terry et al. (Reference Terry, Faber, Hegna, Mirnov, Pueschel and Whelan2018). This scheme yields both an internally consistent expansion, and is consistent with general numerical solutions of the many-wavenumber two-field fluid system as described by Makwana, Terry & Kim (Reference Makwana, Terry and Kim2012). Given the above, the eigenvalues $\omega _j$ can be written as

(2.6) \begin{align} \omega _{1,2} &= \frac {\epsilon _n k_y\big[2+(1+\sqrt {8})k_\perp ^2\big]+k_y - \mathrm{i}\nu k_\perp ^2}{2\big(1+k_\perp ^2\big)} \nonumber \\[4pt] &\pm \frac {\sqrt {-8 \epsilon _n k_y^2(1+\eta )\big(1+k_\perp ^2\big)+\big\{\epsilon _n k_y\big[2\sqrt {8}+(1+\sqrt {8})k_\perp ^2\big]-k_y+\mathrm{i}\nu k_\perp ^2\big\}^2}}{2\big(1+k_\perp ^2\big)}. \end{align}

Note that subscripts $1$ and $2$ refer to unstable/stable eigenmode labels throughout this work. The eigenvalues given by (2.6) follow from the linear dispersion relation (DR) for stable and unstable modes of this system. These modes have drift-wave oscillation, instability growth due to the temperature gradient $\eta$ and magnetic curvature $\epsilon _n$ , an ITG-stabilising effect induced by the density gradient $L_n^{-1}$ and collisional damping due to $\nu$ . It should be noted that $\nu$ acts on both the turbulence and ZFs. Furthermore, $\omega _j$ is often algebraically manipulated throughout this work using the generic complex-number form

(2.7) \begin{align} \omega _j = \varpi _j \pm \mathrm{i} \gamma _j \end{align}

for wave frequency $\varpi _j = \mathrm{Re}(\omega _j)$ and unstable/stable growth/damping rate $\gamma _j=\mathrm{Im}(\omega _j)$ .

The dominant mechanism for saturation of the ITG instability in the models presented herein (Horton et al. Reference Horton, Choi and Tang1981, Reference Horton, Hong and Tang1988; Holland et al. Reference Holland, Diamond, Champeaux, Kim, Gurcan, Rosenbluth, Tynan, Crocker, Nevins and Candy2003; Li & Terry Reference Li and Terry2022) is ZF-catalysed energy transfer from the unstable to stable mode, for which $\gamma _j\neq 0$ (Terry et al. Reference Terry, Faber, Hegna, Mirnov, Pueschel and Whelan2018; Li & Terry Reference Li and Terry2022). In the interest of preserving the key turbulent saturation physics while using the wavenumber-truncated model, we restrict our analysis to choices of eigenmode wavevectors $\boldsymbol{k}$ and $\boldsymbol{k}^{\prime \prime }$ for which $\gamma _j = \mathrm{Im}(\omega _j) \neq 0$ . For suitable choices of $\epsilon _n$ (approximately $.6 \lt \epsilon _n \lt 1.25$ ) and $\eta$ above the threshold of linear instability, these appropriate $\boldsymbol{k}$ s correspond to the radicand of (2.6) being negative. When $\nu =0$ , these are points in the $k_xk_y$ plane satisfying

(2.8) \begin{align} -8 \epsilon _n k_y^2(1+\eta )\big(1+k_\perp ^2\big)+\big\{\epsilon _n k_y\big[2\sqrt {8}+(1+\sqrt {8})k_\perp ^2\big]-k_y\big\}^2 \lt 0. \end{align}

It can be shown that this is the interior of a curve that is quartic in $\boldsymbol{k}_\perp$ . For $\boldsymbol{k}$ outside of this curve, the eigenmodes are at most (for $k_y \neq 0$ ) purely oscillatory dispersive waves and possess no linear behaviour otherwise ( $k_y=0$ ).

One particularly important feature of the DR is that $\omega _1$ and $\omega _2$ constitute a complex-conjugate pair when $\nu =0$ . Both the wave frequencies

(2.9) \begin{align} \varpi _{1,\nu =0} = \varpi _{2,\nu =0} = \varpi \end{align}

and growth/damping rates

(2.10) \begin{align} \gamma _{1,\nu =0} = \gamma _{2,\nu =0} = \gamma \end{align}

for pairs of eigenmodes become equal; their values become eigenmode label independent. For $\nu =0$ , the complex mode frequency can be expressed as

(2.11) \begin{align} \omega _{j,\nu =0} = \varpi \pm \mathrm{i} \gamma . \end{align}

It will prove useful to introduce the exchange matrix

(2.12) \begin{align} J_2 = \begin{bmatrix} 0 & \quad 1 \\[2pt] 1 & \quad 0 \end{bmatrix}\!, \end{align}

which acts on a column vector

(2.13) \begin{align} \boldsymbol{v} = \begin{bmatrix} a \\[2pt] b \end{bmatrix}\!, \end{align}

to exchange its elements

(2.14) \begin{align} J_2 \boldsymbol{v} = \begin{bmatrix} 0 & \quad 1 \\[2pt] 1 & \quad 0 \end{bmatrix} \begin{bmatrix} a \\[2pt] b \end{bmatrix} = \begin{bmatrix} b \\[2pt] a \end{bmatrix}\!. \end{align}

In addition, we define the transformation $\mathcal{T}$ , which acts on a square matrix

(2.15) \begin{align} \boldsymbol{A} = \begin{bmatrix} a & \quad b \\[2pt] c & \quad d \end{bmatrix} \end{align}

as

(2.16) \begin{align} \mathcal{T} (\boldsymbol{A}) = J_2 \boldsymbol{A} J_2 = \begin{bmatrix} 0 & \quad 1 \\[2pt] 1 & \quad 0 \end{bmatrix} \begin{bmatrix} a & \quad b \\[2pt] c & \quad d \end{bmatrix} \begin{bmatrix} 0 & \quad 1 \\[2pt] 1 & \quad 0 \end{bmatrix} = \begin{bmatrix} d & \quad c \\[2pt] b & \quad a \end{bmatrix}. \end{align}

Then from (2.10), for $\nu = 0$ , the linear growth/decay rate for unstable/stable modes possesses invariance under $J_2$ of the form

(2.17) \begin{align} J_2 \gamma _{j,\nu =0} = J_2\begin{bmatrix} \gamma \\[2pt]\gamma \end{bmatrix} =\begin{bmatrix} \gamma \\[2pt]\gamma \end{bmatrix} = \gamma _{j,\nu =0}, \end{align}

which is generally not true for finite collisionality, i.e.

(2.18) \begin{align} J_2\gamma _j =J_2 \begin{bmatrix} \gamma _{1} \\[2pt]\gamma _{2} \end{bmatrix} =\begin{bmatrix} \gamma _{2} \\[2pt] \gamma _{1} \end{bmatrix} \neq \gamma _j. \end{align}

The exchange operation may appear to be quite trivial and insignificant for the time being, but proves to be of fundamental importance in the subsequent analysis of this system for collisionless ZF saturation. The permutation symmetry (Sakurai Reference Sakurai1994; Artin Reference Artin1998) present in the linear dynamics results in an analogous property in the nonlinearity of the ZF, as it is ultimately constructed from arithmetic combinations of the linear eigenvalues of the interacting modes.

Applying the transformation (2.4) to $p_k$ and $\phi _k$ (2.1) and (2.2) result in the system of equations given by

(2.19) \begin{align} \partial _t \beta _{l}^{} + \mathrm{i} \omega _l \beta _{l}^{} &= \sum _{k^{\prime },k^{\prime }_y \neq 0,k_y} C_{lmn}^{({k},{k^{\prime }},{k^{\prime \prime }})} \beta _{m}^{{\prime }}\beta _{n}^{\prime\prime} \nonumber \\ & \quad + \sum _{k_x^{\prime }}\big[\big(C_{lzn}^{({k},{k^{\prime }},{k^{\prime \prime }})}v_z^{{\prime }} + C_{lpn}^{({k},{k^{\prime }},{k^{\prime \prime }})}p_z^{{\prime }}\big)\beta _{n}^{\prime\prime}\delta _{k_y^{\prime }} \nonumber \\ &\quad +\big(C_{lmz}^{({k},{k^{\prime }},{k^{\prime \prime }})}v_z^{\prime\prime} + C_{lmp}^{({k},{k^{\prime }},{k^{\prime \prime }})}p_z^{\prime\prime}\big)\beta _{m}^{{\prime } } \delta _{k_y',k_y}\big] ,\\[-10pt]\nonumber\end{align}

(2.20) \begin{align} \partial _t v_z^{} + \nu v_z^{} &= \sum _{k^{\prime }} C_{zmn}^{({k},{k^{\prime }},{k^{\prime \prime }})} \beta _{m}^{{\prime } }\beta _{n}^{\prime\prime}\delta _{k_y} ,\\[-10pt]\nonumber\end{align}

(2.21) \begin{align} \partial _t p_z + \chi k_x^2 p_z &= \sum _{k^{\prime }} C_{pmn}^{({k},{k^{\prime }},{k^{\prime \prime }})} \beta _{m}^{{\prime } }\beta _{n}^{\prime\prime}\delta _{k_y}, \end{align}

where $\beta _m'=\beta _m(k^{\prime })$ , $\beta _n^{\prime \prime }=\beta _n(k^{\prime \prime })$ , $v_z(k_x) = \mathrm{i} k_x\phi _{k_y=0}=\mathrm{i}k_x\phi _z$ is the $k_y=0$ component of the $E \times B$ flow and $p_z$ is the $k_y=0$ component of the pressure. The nonlinear coupling coefficients used in this work are defined as

(2.22) \begin{align} C_{lzn}^{({k},{k^{\prime }},{k^{\prime \prime }})} = \frac {({-}1)^{l-1}}{R_1(\boldsymbol{k})-R_2(\boldsymbol{k})}\frac {({-}\mathrm{i}k_y)}{2}\biggl [1+R_n(\boldsymbol{k}^{\prime \prime })-\frac {R_{3-l}(\boldsymbol{k})\big(\boldsymbol{k}^{{\prime\prime}2}_\perp -\boldsymbol{k}^{{\prime } 2}_\perp \big)}{1+k_\perp ^2}\biggl ],\\[-10pt]\nonumber \end{align}

(2.23) \begin{align} C_{lmz}^{({k},{k^{\prime }},{k^{\prime \prime }})} = \frac {({-}1)^{l-1}}{R_1(\boldsymbol{k})-R_2(\boldsymbol{k})}\frac {({-}\mathrm{i}k_y')}{2}\biggl [1+R_m(\boldsymbol{k}^{\prime })-\frac {R_{3-l}(\boldsymbol{k})\big(\boldsymbol{k}^{{\prime\prime}2}_\perp -\boldsymbol{k}^{{\prime } 2}_\perp \big)}{1+k_\perp ^2}\biggl ] \end{align}

and

(2.24) \begin{align} C_{zmn}^{({k},{k^{\prime }},{k^{\prime \prime }})} = \frac {-\mathrm{i} k_y'\big(k^{{\prime\prime}2}_\perp -k^{{\prime } 2}_\perp \big)}{2}. \end{align}

Expressions for the remaining coefficients can be found in Terry et al. (Reference Terry, Faber, Hegna, Mirnov, Pueschel and Whelan2018). The system of equations (2.19)–(2.21) reformulates (2.1) and (2.2) in terms of the nonlinearly coupled time-evolution equations for the stable/unstable modes ( $k_y \neq 0$ ), ZF $v_z$ ( $k_y=0$ ) and zonal pressure $p_z$ (also $k_y = 0$ ), which interact with one another. In particular, (2.19) separates the nonlinearity of the evolution equation for stable/unstable modes into interactions that involve a ZF, and those that do not. Equations (2.19)–(2.21) have been extensively studied in the presence of finite collisionality, notably making agreement with the well-known property of gyrokinetic solutions in which turbulence levels are directly proportional to the rate of ZF damping (Lin et al. Reference Lin, Hahm, Lee, Tang and Diamond1999; Terry et al. Reference Terry, Faber, Hegna, Mirnov, Pueschel and Whelan2018). They have also been assessed for nonlinear saturation scalings and used to explain heat flux reduction above the linear critical gradient of toroidal ITG-driven turbulence (Terry et al. Reference Terry, Faber, Hegna, Mirnov, Pueschel and Whelan2018, Reference Terry, Li, Pueschel and Whelan2021) in the phenomenon known as the Dimits shift (Dimits et al. Reference Dimits, Cohen, Mattor, Nevins, Shumaker, Parker and Kim2000a , Reference Dimitsb ).

More recently, the system of equations (2.19)–(2.21) has been studied for a truncated wavenumber space limited to three interacting wavenumbers and five eigenmodes $\beta _l=\beta _l(\boldsymbol{k})$ , $\beta _l^{\prime \prime }=\beta _l(\boldsymbol{k}^{\prime \prime })$ and $v_z'=v_z(k^{\prime })$ . The three-wavenumber system takes the form

(2.25) \begin{align} \dot {\beta _l} +\mathrm{i}\omega _l \beta _l = C_{lzn}^{({k},{k^{\prime }},{k^{\prime \prime }})} v_z^{{\prime } }\beta _{n}^{\prime\prime} ,\\[-5pt]\nonumber\end{align}

(2.26) \begin{align} \dot {\beta _{l}^{\prime\prime}} +\mathrm{i}\omega _{l}^{\prime\prime} \beta _{l}^{\prime\prime} = C_{lzn}^{({k^{\prime \prime }},{-k^{\prime }},{k})} v_z^{{\prime }*}\beta _{n}^{} ,\\[-5pt]\nonumber\end{align}

(2.27) \begin{align} \dot {v_z^{{\prime }}} + \nu v_z^{{\prime }} = C_{zmn}^{({k^{\prime }},{k},{-k^{\prime \prime }})}\beta _{m}^{}\beta _{n}^{{\prime\prime}*}. \end{align}

Assumptions were made in deriving (2.25)–(2.27) based on an ordering scheme inferred from gyrokinetic simulations (Terry et al. Reference Terry, Faber, Hegna, Mirnov, Pueschel and Whelan2018; Li & Terry Reference Li and Terry2022), which allows for discarding the nonlinearities in (2.19) that do not involve a ZF. Equations (2.25)–(2.27) have been shown to reproduce fundamental scalings of toroidal ITG-driven turbulence (Li & Terry Reference Li and Terry2022) and have been used to investigate limit cycle oscillations in these systems (Li & Terry Reference Li and Terry2024). While recovering the scalings of the turbulence with model parameters, the three-wavenumber truncation does not reproduce realistic saturation amplitudes (Li & Terry Reference Li and Terry2022). Given its relative simplicity and demonstrated utility, this system is used herein as a reduced model for considering collisionless ZF saturation. It is not difficult to show that the results generalise in applicability to the case of many interacting wavenumbers, as described by (2.19)–(2.21).

The collisionless truncated system presents a simple model for directly demonstrating that turbulence saturation via stable-mode damping is entirely distinct phenomenologically from the usual standardised picture presented by Kolmogorov (Reference Kolmogorov1941). The latter formalism assumes the existence of separate sets of length scales over which energy is injected, nonlinearly transferred between fluctuations in the form of a cascade and dissipated. A notable feature of these types of systems is that the energy injection rate is, in general, independent of the level of dissipation.

Turbulence described by systems of interactions among stable and unstable eigenmodes differs dramatically. For instance, there is not a clear distinction between the characteristic length scales over which injection, nonlinear interaction and dissipation occur; in general, they are overlapping and occur simultaneously at each scale. In addition, consider the complex mode frequency equation (2.6). Two systems with the same values of $\eta$ and $\epsilon _n$ but different $\nu$ values will in general experience different overall energy injection rates, as would be described approximately by $\sum _{\boldsymbol{k}} \gamma _1(\boldsymbol{k})$ .

Furthermore, in the wavenumber-truncated system, a cascade does not occur; the energy is only exchanged among the five modes. This trait then reveals a significant way in which the turbulent dynamics of the model differ from and fail to be characterised by a cascade. Given that, for $\nu =0$ , there is no collisional damping to thermalise energy, there must be some other process present to balance the overall injection of energy and bring about saturation. The wavenumber-truncated collisionless dynamics illustrate an alternative process whose essence is entirely distinct physically from a Kolmogorov cascade. This process is the ZF-catalysed energy transfer from the unstable to stable mode, a nonlinear process, which in the limit of $\nu \rightarrow 0$ puts energy back into the temperature profile at exactly the rate in which it was extracted by the instability.

Before proceeding, it is briefly acknowledged that the nonlinear coupling coefficients in (2.25)–(2.27) constitute tensors, e.g.

(2.28) \begin{align} C_{lzn}^{({k},{k^{\prime }},{k^{\prime \prime }})} & = \begin{bmatrix} C_{1z1}^{({k},{k^{\prime }},{k^{\prime \prime }})} & \quad C_{1z2}^{({k},{k^{\prime }},{k^{\prime \prime }})} \\[8pt] C_{2z1}^{({k},{k^{\prime }},{k^{\prime \prime }})} & \quad C_{2z2}^{({k},{k^{\prime }},{k^{\prime \prime }})} \end{bmatrix} \nonumber \\[5pt]& = \frac {-\mathrm{i}k_y}{2(R_1(\boldsymbol{k})-R_2(\boldsymbol{k}))} \nonumber \\[5pt]& \quad \times \begin{bmatrix} \left[1+R_1(\boldsymbol{k}^{\prime \prime })-\dfrac {R_{2}(\boldsymbol{k})(\boldsymbol{k}^{{\prime\prime}2}_\perp -\boldsymbol{k}^{{\prime }2}_\perp )}{1+k_\perp ^2}\right] & \quad\!\!\! \left[1+R_2(\boldsymbol{k}^{\prime \prime })-\dfrac {R_{2}(\boldsymbol{k})(\boldsymbol{k}^{{\prime\prime}2}_\perp -\boldsymbol{k}^{{\prime }2}_\perp )}{1+k_\perp ^2}\right] \\[20pt] -\left[1+R_1(\boldsymbol{k}^{\prime \prime })-\dfrac {R_{1}(\boldsymbol{k})(\boldsymbol{k}^{{\prime\prime}2}_\perp -\boldsymbol{k}^{{\prime }2}_\perp )}{1+k_\perp ^2}\right] & \quad\!\!\! -\left[1+R_2(\boldsymbol{k}^{\prime \prime })-\dfrac {R_{1}(\boldsymbol{k})(\boldsymbol{k}^{{\prime\prime}2}_\perp -\boldsymbol{k}^{{\prime }2}_\perp )}{1+k_\perp ^2}\right] \end{bmatrix}\!,\nonumber\\[5pt] \end{align}

and accordingly, terms with repeated indices imply sums

(2.29) \begin{align} \dot {\beta _l} +\mathrm{i}\omega _l \beta _l = \big(C_{lz1}^{({k},{k^{\prime }},{k^{\prime \prime }})} \beta _{1}^{\prime\prime} + C_{lz2}^{({k},{k^{\prime }},{k^{\prime \prime }})} \beta _{2}^{\prime\prime}\big)v_z^{{\prime }}, \end{align}

indicating that the modes $\beta _{l}^{}$ evolve due to ZF-catalysed nonlinear interactions with both $\beta _{1}^{\prime\prime}$ and $\beta _{2}^{\prime\prime}$ .

3. Weak-turbulence closure for ZF WKE

Analytical prediction of collisionless ZF saturation necessitates an expression for the time evolution of the energy $\langle v_z^{{\prime }2}\rangle$ of the ZF. Note that in this formalism, quantities within angle brackets $\langle . \rangle$ describe ensemble averages over many realisations of the system (Malfliet Reference Malfliet1972). As such, the starting point for analysis in this work is application of an instability-driven weak-turbulence closure technique (Malfliet Reference Malfliet1972) to the $3-\boldsymbol{k}$ system of (2.25)–(2.27). The closure calculates a WKE characterising ZF energy time evolution $\partial _t \langle v_z^{{\prime }2} \rangle$ . The details of this calculation are somewhat lengthy, but otherwise largely trivial to execute, and are given in Appendix A.

The ZF WKE is found to be

(3.1) \begin{align} \partial _t \langle v_z^{{\prime } 2} \rangle = - 2 \nu \langle v_z^{{\prime } 2} \rangle + 2 \mathrm{Re} \big\{ C_{zij}^{({k^{\prime }},{k},{-k^{\prime \prime }})} \big[ \langle v_z^{{\prime }2} \rangle \big( &C_{izm}^{({k},{k^{\prime }},{k^{\prime \prime }})} \tau _{izm}^{kk^{\prime }k^{\prime \prime }}\langle \beta _{j}^{{\prime\prime}*} \beta _{m}^{\prime\prime} \rangle \nonumber \\& + C_{jzm}^{({k^{\prime \prime }},{-k^{\prime }},{k})*} \tau _{jzm}^{k^{\prime \prime }k^{\prime }k}\langle \beta _{m}^{*}\beta _{i}^{} \rangle \big) \big]\big\}, \end{align}

where $\tau$ are the triplet correlation times (TCTs)

(3.2) \begin{align} \tau _{izm}^{kk^{\prime }k^{\prime \prime }} = \frac {1}{\mathrm{i}\big(\omega _{i}^{} - \omega _z - \omega _{m}^{\prime\prime}\big)} \end{align}

and

(3.3) \begin{align} \tau _{jzm}^{k^{\prime \prime }k^{\prime }k*} = \frac {1}{-\mathrm{i}(\omega _{j}^{\prime\prime}-\omega _z-\omega _{m}^{})^*}, \end{align}

and the ZF damping rate $\nu$ has been written as a wavenumber-independent complex frequency $\omega _z = -\mathrm{i}\nu$ .

Equation (3.1) describes how the ZF energy $\langle v_z^{{\prime }2} \rangle$ fluctuates with time according to linear collisional damping $\nu$ and all possible interactions with the eigenmodes $\beta _{l}^{}$ and $\beta _{j}^{\prime\prime}$ . Interactions occur with strength in proportion to the associated nonlinear coupling coefficients, TCTs (characteristic time scales over which the nonlinear three-wave processes transpire) and second-order correlation functions specific to each process. Equation (3.1) is valid for systems satisfying the weak-turbulence assumption made in its derivation. This corresponds to statistics obeying (A19)–(A21). Essentially, the statistical moment hierarchy is ordered by increasing polynomial degree in a small dimensionless parameter $\epsilon \ll 1$ ; the nth-order moment in statistics is $O(\epsilon ^{n-1})$ in perturbation. Through the lens of $n$ -point spatial correlations, this is in part a statement regarding the ‘randomness’ of the turbulence in real space; three points in the flow will be on average much less ‘similar’ to each other than two points, but far more than four. In addition, there is an implicit requirement that the various underlying asymptotic expansions (e.g. (A22), (A24)–(A26) and (A29)) used in the derivation of (3.1) are well ordered.

It is important for the ensuing analysis to illuminate the tensorial nature of the second-order correlations, e.g.

(3.4) \begin{align} \big\langle \beta _{j}^{{\prime\prime}*} \beta _{m}^{\prime\prime} \big\rangle = \begin{bmatrix} \langle \beta _{1}^{{\prime\prime} 2} \rangle & \quad \big\langle \beta _{1}^{{\prime\prime}*} \beta _{2}^{\prime\prime} \big\rangle \\[5pt] \big\langle \beta _{2}^{{\prime\prime}*} \beta _{1}^{\prime\prime} \rangle & \quad \langle \beta _{2}^{{\prime\prime}2} \big\rangle \end{bmatrix}, \end{align}

as well as $\tau$ s:

(3.5) \begin{align} \tau _{izm}^{kk^{\prime }k^{\prime \prime }} =\frac {1}{\mathrm{i}} \begin{bmatrix} (\omega _{1}^{} - \omega _z - \omega _{1}^{\prime\prime})^{-1} & \quad (\omega _{1}^{} - \omega _z - \omega _{2}^{\prime\prime})^{-1} \\[3pt] (\omega _{2}^{} - \omega _z - \omega _{1}^{\prime\prime})^{-1} & \quad (\omega _{2}^{} - \omega _z - \omega _{2}^{\prime\prime})^{-1} \end{bmatrix}\!. \end{align}

The time-evolution operator on the right-hand side of (3.1) is a mapping from the system’s space of second-order correlations to the instantaneous rate of growth or decay in ZF energy. The overall operator is a product of individual tensors expressing both the system’s present location in phase space (e.g. (3.4)) and the model’s constants of proportionality corresponding to each coordinate (e.g. (2.28) and (3.5)). Rather than being literal positions in configuration space, the coordinates describe the instantaneous amounts of energy moving both up and down gradient ( $\langle \beta _1^2\rangle$ and $\langle \beta _2^2\rangle$ , respectively) and the correlation between the opposing fluctuations (e.g. $\langle \beta _1^*\beta _2\rangle$ ). The presence of terms with repeated indices in (3.1) then implies sums over the different nonlinear processes present in the dynamics.

Equation (3.1) can be written as

(3.6) \begin{align} \partial _t \langle v_z^{{\prime } 2} \rangle = - 2( \nu - \kappa ) \langle v_z^{{\prime } 2} \rangle , \end{align}

where $\kappa =\kappa (k,k^{\prime })$ is a turbulent rate of energy transfer to and/or from the ZF. From (3.1) the rate $\kappa$ is given by

(3.7) \begin{align} \kappa = \mathrm{Re} \big[ C_{zij}^{({k^{\prime }},{k},{-k^{\prime \prime }})} \big ( C_{izm}^{({k},{k^{\prime }},{k^{\prime \prime }})} \tau _{izm}^{kk^{\prime }k^{\prime \prime }}\langle \beta _{j}^{{\prime\prime}*} \beta _{m}^{\prime\prime} \rangle + C_{jzm}^{({k^{\prime \prime }},{-k^{\prime }},{k})*} \tau _{jzm}^{k^{\prime \prime }k^{\prime }k*}\langle \beta _{m}^{*}\beta _{i}^{} \rangle \big ) \big ]. \end{align}

Note the weak-turbulence closure’s slow time variation of second-order correlations described by (A37). To lowest non-trivial order in perturbation, solutions of the WKE initialised from some state at time $t=0$ are of the form

(3.8) \begin{align} \langle v_z^{{\prime }2} \rangle = \langle v_z^{{\prime }2} \rangle _{t=0}\textrm {e}^{- 2( \nu - \kappa )t}. \end{align}

Conditions for ZF saturation correspond to stationary solutions of the WKE (3.6) and arguments in the exponential of (3.8) for which the time dependence vanishes. In general, this condition is

(3.9) \begin{align} \nu = \kappa , \end{align}

which in the absence of collisionality becomes

(3.10) \begin{align} \kappa = 0. \end{align}

The collisionless saturation condition (3.10) suggests that the ZF of the three-wavenumber system achieves saturation in the limit $\nu =0$ via balances between sources of drive and damping that are purely nonlinear in nature. This raises the question whether terms can be found in $\kappa$ that cause the ZF to grow and others that instigate decay. Noting that $v_z'$ interacts with unstable modes $\beta _{1}^{}$ and $\beta _{1}^{\prime\prime}$ and stable modes $\beta _{2}^{}$ and $\beta _{2}^{\prime\prime}$ , a related questions is whether a corresponding symmetric balance can be found between the sources of growth and decay to bring about a saturation. What follows is an analysis of the mathematical properties of (3.6), in particular, the structure of the nonlinear rate $\kappa$ , to investigate possible means of collisionless saturation and to compare with the collisional case. This investigation answers the questions just raised affirmatively.

4. Nonlinear time-evolution operator

The rate $\kappa$ quantifies nonlinear energy transfer to the ZF in terms of a set of resonant interactions between triads of fluctuations. A rigorous understanding of these details is necessary for predicting explicit stationary solutions of the ZF WKE and, in the limit $\nu \rightarrow 0$ , collisionless saturation. In addition, analysis is performed to clarify the workings of the nonlinearity, by expressing $\kappa$ entirely in terms of the more familiar and workable linear physics of the individual modes participating in a given interaction. An important point is that the eigenvalue or mode frequency, (2.6), exudes exchange symmetry of the form represented by (2.17) in the limit of zero collisionality. As a result, the nonlinearity (which is shown to be arithmetic combinations of the linear complex eigenfrequencies of individual modes) possesses centrosymmetry.

Investigation of the saturation physics contained in the ZF WKE begins with analysis of the properties of $\kappa$ as described by (3.7). Recalling the form of $\kappa$ , it is first observed that the ZF coupling coefficient,

(4.1) \begin{align} C_{zij}^{({k^{\prime }},{k},{-k^{\prime \prime }})} = \frac {\mathrm{i}k_y(k^2-k^{\prime \prime 2})}{2}, \end{align}

is purely imaginary. Therefore, only the imaginary portions of the terms of the form $C\langle \beta \beta \rangle \tau$ within the $\{ \}$ braces of (3.7) need be considered. Anticipating the general structure of the individual terms in braces, a complex-number representation of the form

(4.2) \begin{align} C\langle \beta \beta \rangle \tau = (C_R + \mathrm{i} C_I) (\langle \beta \beta \rangle _R+\mathrm{i}\langle \beta \beta \rangle _I)(\tau _R+\mathrm{i}\tau _I) \end{align}

is suggested. The imaginary portion of either of these terms $C\langle \beta \beta \rangle \tau$ will then resemble

(4.3) \begin{align} \mathrm{Im}(C\langle \beta \beta \rangle \tau ) & = (C_R \tau _I +C_I \tau _R)\langle \beta \beta \rangle _R+ (C_R\tau _R-C_I\tau _I)\langle \beta \beta \rangle _I \nonumber \\& = T_I\langle \beta \beta \rangle _R+T_R\langle \beta \beta \rangle _I, \end{align}

where

(4.4) \begin{align} T_I = C_R \tau _I +C_I \tau _R,\\[-7pt]\nonumber \end{align}

(4.5) \begin{align} T_R = C_R\tau _R-C_I\tau _I, \end{align}

and the shorthand notation

(4.6) \begin{align} \boldsymbol{T} \langle \beta \beta \rangle \equiv T_I\langle \beta \beta \rangle _R+T_R\langle \beta \beta \rangle _I \end{align}

is introduced. The condition for collisionless ZF saturation, $\kappa =0$ , is that terms of the form $\boldsymbol{T} \langle \beta \beta \rangle$ in $\kappa$ individually equate to zero as

(4.7) \begin{align} \boldsymbol{T}_{izm}^{k, k^{\prime }, k^{\prime \prime }} \langle \beta _{j}^{{\prime\prime}*} \beta _{m}^{\prime\prime}\rangle = 0 ,\\[-7pt]\nonumber \end{align}

(4.8) \begin{align} \boldsymbol{T}_{jzm}^{k^{\prime \prime }, -k^{\prime }, k (*)} \langle \beta _{m}^{*} \beta _{i}^{}\rangle = 0 , \end{align}

where the tensors $\boldsymbol{T}_{izm}^{k, k^{\prime }, k^{\prime \prime }}$ and $\boldsymbol{T}_{jzm}^{k^{\prime \prime }, -k^{\prime }, k (*)}$ are defined as

(4.9) \begin{align} \boldsymbol{T}_{izm}^{k, k^{\prime }, k^{\prime \prime }} &= C_{izm}^{({k},{k^{\prime }},{k^{\prime \prime }})} \tau _{izm}^{kk^{\prime }k^{\prime \prime }} ,\\[-7pt]\nonumber \end{align}

(4.10) \begin{align} \boldsymbol{T}_{jzm}^{k^{\prime \prime }, -k^{\prime }, k (*)} &= C_{jzm}^{({k^{\prime \prime }},{-k^{\prime }},{k})*} \tau _{jzm}^{k^{\prime \prime }k^{\prime }k*}. \end{align}

The notation $\boldsymbol{T}_{jzm}^{k^{\prime \prime }, -k^{\prime }, k (*)}$ represents complex conjugation of the corresponding nonlinear coupling coefficient and not the entire tensor. Alternatively to (4.7) and (4.8), the sum of the two terms may combine to zero, suggesting a different, more general saturation of the form

(4.11) \begin{align} \boldsymbol{T}_{izm}^{k, k^{\prime }, k^{\prime \prime }} \langle \beta _{j}^{{\prime\prime}*} \beta _{m}^{\prime\prime}\rangle + \boldsymbol{T}_{jzm}^{k^{\prime \prime }, -k^{\prime }, k (*)} \langle \beta _{m}^{*} \beta _{i}^{}\rangle = 0. \end{align}

Clearly, either criterion necessitates understanding of the structures of the tensors $\boldsymbol{T}_{izm}^{k, k^{\prime }, k^{\prime \prime }}$ and $\boldsymbol{T}_{jzm}^{k^{\prime \prime }, -k^{\prime }, k (*)}$ .

We begin by defining their real and imaginary parts with the expressions

(4.12) \begin{align} T_{R,izm}^{k, k^{\prime }, k^{\prime \prime }} &= \mathrm{Re}\big(C_{izm}^{({k},{k^{\prime }},{k^{\prime \prime }})} \tau _{izm}^{kk^{\prime }k^{\prime \prime }}\big),\\[-7pt]\nonumber \end{align}

(4.13) \begin{align} T_{I,izm}^{k, k^{\prime }, k^{\prime \prime }} &= \mathrm{Im}\big( C_{izm}^{({k},{k^{\prime }},{k^{\prime \prime }})} \tau _{izm}^{kk^{\prime }k^{\prime \prime }} \big),\\[-10pt]\nonumber \end{align}

(4.14) \begin{align} T_{R,jzm}^{k^{\prime \prime }, -k^{\prime }, k (*)} &= \mathrm{Re}\big(C_{jzm}^{({k^{\prime \prime }},{-k^{\prime }},{k})*} \tau _{jzm}^{k^{\prime \prime }k^{\prime }k*}\big), \\[-10pt]\nonumber\end{align}

(4.15) \begin{align} T_{I,jzm}^{k^{\prime \prime }, -k^{\prime }, k (*)} &= \mathrm{Im}\big(C_{jzm}^{({k^{\prime \prime }},{-k^{\prime }},{k})*} \tau _{jzm}^{k^{\prime \prime }k^{\prime }k*}\big). \end{align}

The above tensors are calculated in Appendix B. Their collisional forms are given by (B23), (B24), (B32) and (B33). They express the products of coupling coefficients and TCTs found in the nonlinear portion of the ZF spectral time-evolution operator in terms of the linear eigenvalues and eigenvector components of the interacting modes; the nature of the turbulence and associated resonant nonlinear energy transfer follows from the particular way the system combines individual wave frequencies and growth/damping rates. A notable commonality among all the matrix elements is clearly visible in their denominators. The factors of the form

(4.16) \begin{align} T \propto \frac {1}{(\varpi _a^{\prime \prime }-\varpi _b)^2} \end{align}

suggest that the ZF is more strongly influenced by interactions with oscillating modes $\beta$ and $\beta ^{\prime \prime }$ that are closer together in wave frequency, lending to a denominator that is closer to vanishing. Meanwhile, the sums and differences among $\gamma$ s and $\nu$ s are of the form

(4.17) \begin{align} T^{k, k^{\prime }, k^{\prime \prime }}_{azb} \propto \frac {1}{(\mp \gamma _a - \nu \pm \gamma _b^{\prime\prime})^2} \end{align}

and

(4.18) \begin{align} T^{k, k^{\prime }, k^{\prime \prime }}_{azb} \propto \frac {1}{(\mp \gamma _a^{\prime \prime } - \nu \pm \gamma _b)^2}. \end{align}

These functions are eigenmode label dependent. While $\nu$ is always subtractive, a stable or unstable mode may be additive or subtractive from the denominator. Hence, stable-mode damping and unstable-mode growth may contribute to or diminish the strength of an interaction, depending on how they relate to $\nu$ via the sign that precedes them.

The other significant feature is the consistent presence of $\nu$ in the numerators of the matrix elements in the eigenmode-label-independent forms $-\nu (\hat {\omega }^{\prime \prime }_\perp -k_y^2 \omega _y)$ and $-\nu (\hat {\omega }_\perp -\hat {\omega }_{\perp }^{\prime\prime} \omega _{\Delta }^{k^{\prime \prime },-k^{\prime }})$ . These $\nu$ -related terms always exhibit the same arithmetic behaviour of subtraction, regardless of matrix element, and therefore, help to break the permutation symmetries among pairs of diagonal or off-diagonal terms in a given $T_{izm}$ or $T_{jzm}$ . As such, the resulting energetic dynamics exhibited by solutions of the reduced two-field fluid models of (2.19)–(2.21) and (2.25)–(2.27) are not symmetric. To whit, in the quasi-stationary regime, the unstable-mode energy spectrum saturates at a larger value than the corresponding stable spectrum does at each $\boldsymbol{k}$ . The discrepancy between these averages and associated broken permutation symmetries occur in direct proportion to $\nu$ (Terry et al. Reference Terry, Faber, Hegna, Mirnov, Pueschel and Whelan2018, Reference Terry, Li, Pueschel and Whelan2021; Li & Terry Reference Li and Terry2022).

For systems with non-zero collisionality, $\nu \neq 0$ $\implies$ $ \gamma _1 \neq \gamma _2$ , $\gamma _1^{\prime \prime } \neq \gamma _2$ and the complex linear eigenmode frequency (2.6) exudes broken symmetry via $\nu$ . Consequentially, the tensors of (B23), (B24), (B32) and (B33) (made of linear eigenvalues) are non-trivial convoluted combinations of symmetric and antisymmetric parts (as is true for any matrix). Such a decomposition into these respective parts would likely result in production of expressions of increased complexity with reduced analytical tractability. Contrastingly, for $\nu =0$ , the unstable/stable growth/damping rates are symmetric under eigenmode-label permutation such that $\varpi _1=\varpi _2$ , $\gamma _1=\gamma _2$ , $\varpi _1^{\prime \prime } = \varpi _2^{\prime \prime }$ and $\gamma _1^{\prime\prime}=\gamma _2^{\prime\prime}$ . The tensors in this limit are described by (B26), (B27), (B35) and (B36). In matrix form, their structures are

(4.19) \begin{align} T_{R,izm}^{k, k^{\prime }, k^{\prime \prime }} &= \dfrac {k_y }{4\hat {\gamma }} \nonumber \\[5pt] &\quad \times \begin{bmatrix} \dfrac {\hat {\delta }^{k,k^{\prime \prime }}( \gamma - \gamma ^{\prime\prime}) + \big( \hat {\gamma }^{\prime \prime } + \hat {\gamma } \omega _{\Delta }^{k,k^{\prime }} \big)(\varpi ^{\prime \prime } - \varpi )}{( \gamma - \gamma ^{\prime\prime})^2 + (\varpi ^{\prime \prime } - \varpi )^2} & \quad \dfrac {\hat {\delta }^{k,k^{\prime \prime }}( \gamma + \gamma ^{\prime\prime}) - \big( \hat {\gamma }^{\prime \prime } - \hat {\gamma } \omega _{\Delta }^{k,k^{\prime }} \big)(\varpi ^{\prime \prime } - \varpi )}{( \gamma + \gamma ^{\prime\prime})^2 + (\varpi ^{\prime \prime } - \varpi )^2} \\[13pt] \dfrac {\hat {\delta }^{k,k^{\prime \prime }}( \gamma + \gamma ^{\prime\prime}) - \big( \hat {\gamma }^{\prime \prime } - \hat {\gamma } \omega _{\Delta }^{k,k^{\prime }} \big)(\varpi ^{\prime \prime } - \varpi )}{( \gamma + \gamma ^{\prime\prime})^2 + (\varpi ^{\prime \prime } - \varpi )^2} & \quad \dfrac {\hat {\delta }^{k,k^{\prime \prime }}( \gamma - \gamma ^{\prime\prime}) + \big( \hat {\gamma }^{\prime \prime } + \hat {\gamma } \omega _{\Delta }^{k,k^{\prime }} \big)(\varpi ^{\prime \prime } - \varpi )}{( \gamma - \gamma ^{\prime\prime})^2 + (\varpi ^{\prime \prime } - \varpi )^2} \end{bmatrix},\nonumber\\[5pt] \end{align}

(4.20) \begin{align} T_{I,izm}^{k, k^{\prime }, k^{\prime \prime }}& = \dfrac {k_y }{4\hat {\gamma }} \nonumber \\ &\quad\times \begin{bmatrix} - \dfrac {\left(\hat {\delta }^{k,k^{\prime \prime }}(\varpi ^{\prime \prime } - \varpi ) - \left( \hat {\gamma }^{\prime \prime } + \hat {\gamma } \omega _{\Delta }^{k,k^{\prime }} \right)( \gamma - \gamma ^{\prime\prime})\right) }{( \gamma - \gamma ^{\prime\prime})^2 + (\varpi ^{\prime \prime } - \varpi )^2} & - \dfrac {\big(\hat {\delta }^{k,k^{\prime \prime }}(\varpi ^{\prime \prime } - \varpi ) - \big({-}\hat {\gamma }^{\prime \prime } + \hat {\gamma } \omega _{\Delta }^{k,k^{\prime }} \big)( \gamma + \gamma ^{\prime\prime}) \big)}{(\gamma + \gamma ^{\prime\prime})^2 + (\varpi ^{\prime \prime } - \varpi )^2} \\[10pt] \dfrac {\hat {\delta }^{k,k^{\prime \prime }}(\varpi ^{\prime \prime } - \varpi ) - \big( -\hat {\gamma }^{\prime \prime } + \hat {\gamma } \omega _{\Delta }^{k,k^{\prime }} \big)( \gamma +\gamma ^{\prime\prime}) }{( \gamma + \gamma ^{\prime\prime})^2 + (\varpi ^{\prime \prime } - \varpi )^2} & \dfrac {\hat {\delta }^{k,k^{\prime \prime }}(\varpi ^{\prime \prime } - \varpi ) - \big( \hat {\gamma }^{\prime \prime } + \hat {\gamma } \omega _{\Delta }^{k,k^{\prime }} \big)( \gamma -\gamma ^{\prime\prime}) }{( \gamma - \gamma ^{\prime\prime})^2 + (\varpi ^{\prime \prime } - \varpi )^2} \end{bmatrix},\nonumber\\[5pt] \end{align}

(4.21) \begin{align} T_{R,jzm}^{k^{\prime \prime }, -k^{\prime }, k (*)} &= \dfrac {k^{\prime\prime}_y}{4\hat {\gamma }^{\prime\prime}} \nonumber \\ &\quad\times \begin{bmatrix} \dfrac {\hat {\delta }^{k^{\prime \prime },-k^{\prime }} ( \gamma ^{\prime\prime} - \gamma ) - \big( \hat {\gamma } + \hat {\gamma }^{\prime \prime } \omega _{\Delta }^{k^{\prime \prime },-k^{\prime }} \big)(\varpi ^{\prime \prime }-\varpi ) }{( \gamma ^{\prime\prime} - \gamma )^2 + (\varpi ^{\prime \prime }-\varpi )^2 } & \dfrac {\hat {\delta }^{k^{\prime \prime },-k^{\prime }} ( \gamma ^{\prime\prime} + \gamma ) - \big( -\hat {\gamma } + \hat {\gamma }^{\prime \prime } \omega _{\Delta }^{k^{\prime \prime },-k^{\prime }} \big)(\varpi ^{\prime \prime }-\varpi ) }{( \gamma ^{\prime\prime} + \gamma )^2 + (\varpi ^{\prime \prime }-\varpi )^2 } \\[10pt] \dfrac {\hat {\delta }^{k^{\prime \prime },-k^{\prime }} ( \gamma ^{\prime\prime} + \gamma ) - \big( -\hat {\gamma } + \hat {\gamma }^{\prime \prime } \omega _{\Delta }^{k^{\prime \prime },-k^{\prime }} \big) (\varpi ^{\prime \prime }-\varpi )}{( \gamma ^{\prime\prime} + \gamma )^2 + (\varpi ^{\prime \prime }-\varpi )^2 } & \dfrac {\hat {\delta }^{k^{\prime \prime },-k^{\prime }} ( \gamma ^{\prime\prime} - \gamma ) - \big(\hat {\gamma } + \hat {\gamma }^{\prime \prime } \omega _{\Delta }^{k^{\prime \prime },-k^{\prime }} \big)(\varpi ^{\prime \prime }-\varpi ) }{( \gamma ^{\prime\prime} - \gamma )^2 + (\varpi ^{\prime \prime }-\varpi )^2 } \end{bmatrix}\nonumber\\[5pt]\end{align}

and

(4.22) \begin{align} T_{I,jzm}^{k^{\prime \prime }, -k^{\prime }, k (*)} &= \frac {k^{\prime\prime}_y}{4\hat {\gamma }^{\prime\prime}} \nonumber \\ &\quad \times \begin{bmatrix} - \dfrac {\big(\hat {\delta }^{k^{\prime \prime },-k^{\prime }} (\varpi ^{\prime \prime }-\varpi ) + \big( \hat {\gamma } +\hat {\gamma }^{\prime \prime } \omega _{\Delta }^{k^{\prime \prime },-k^{\prime }} \big)( \gamma ^{\prime\prime} - \gamma )\big)}{( -\gamma ^{\prime\prime} +\gamma )^2 + (\varpi ^{\prime \prime }-\varpi )^2} & -\dfrac {\big(\hat {\delta }^{k^{\prime \prime },-k^{\prime }} (\varpi ^{\prime \prime }-\varpi ) + \big( -\hat {\gamma } + \hat {\gamma }^{\prime \prime } \omega _{\Delta }^{k^{\prime \prime },-k^{\prime }} \big)( \gamma ^{\prime\prime} + \gamma )\big)}{( \gamma ^{\prime\prime} + \gamma )^2 + (\varpi ^{\prime \prime }-\varpi )^2} \\[10pt] \dfrac {\hat {\delta }^{k^{\prime \prime },-k^{\prime }} (\varpi ^{\prime \prime }-\varpi ) + \big( -\hat {\gamma } + \hat {\gamma }^{\prime \prime } \omega _{\Delta }^{k^{\prime \prime },-k^{\prime }} \big)( \gamma ^{\prime\prime} + \gamma )}{( \gamma ^{\prime\prime} + \gamma )^2 + (\varpi ^{\prime \prime }-\varpi )^2} & \dfrac {\hat {\delta }^{k^{\prime \prime },-k^{\prime }} (\varpi ^{\prime \prime }-\varpi ) + \big( \hat {\gamma } + \hat {\gamma }^{\prime \prime } \omega _{\Delta }^{k^{\prime \prime },-k^{\prime }} \big)( \gamma ^{\prime\prime} - \gamma )}{( \gamma ^{\prime\prime} - \gamma )^2 + (\varpi ^{\prime \prime }-\varpi )^2} \end{bmatrix}.\nonumber\\[5pt] \end{align}

It is observed that the real tensors $T_R$ satisfy

(4.23) \begin{align} T_{R,1z1} = T_{R,2z2}, \\[-10pt]\nonumber\end{align}

(4.24) \begin{align} T_{R,1z2} = T_{R,2z1} \end{align}

and are centrosymmetric:

(4.25) \begin{align} \mathcal{T}(\textsf{T}_{\textsf{R}}) = \textsf{T}_{\textsf{R}}. \end{align}

In contrast, the imaginary $T_I$ obey

(4.26) \begin{align} T_{I,1z1} &= -T_{I,2z2},\\[-10pt]\nonumber \end{align}

(4.27) \begin{align} T_{I,1z2} &= -T_{I,2z1} \end{align}

and are skew centrosymmetric:

(4.28) \begin{align} \mathcal{T} ({\textsf{T}}_{\textsf{I}}) = -{\textsf{T}}_{\textsf{I}}. \end{align}

These observations imply that the constants of proportionality dictating nonlinear drive and damping of the ZF energy $\langle v_z^{{\prime } 2} \rangle$ become highly symmetric/antisymmetric in the collisionless limit. Correlations that are related through index permutation, e.g. $\langle \beta _{1}^{2} \rangle$ and $\langle \beta _{2}^{2} \rangle$ or $\langle \beta _{1}^{{\prime\prime}*} \beta _{2}^{\prime\prime}\rangle$ and $\langle \beta _{2}^{{\prime\prime}*} \beta _{1}^{\prime\prime}\rangle$ , will make contributions of equal magnitude to ZF spectral fluctuations, ‘per unit correlation’. Of particular importance is the skew-centrosymmetric property of the imaginary tensors, as they multiply the energy spectra (see (4.6)). This means that spectra for the stable/unstable mode at a given $\boldsymbol{k}$ are balanced on average, e.g. $\langle \beta _{1}^{2} \rangle =\langle \beta _{2}^{2} \rangle$ or $\langle \beta _{1}^{{\prime\prime}2} \rangle =\langle \beta _{2}^{{\prime\prime}2} \rangle$ , and therefore, add and remove energy from the ZF at equal and opposite rates, thus, nonlinear growth is saturated by an equal and opposite nonlinear decay. (This is precisely what is found in the subsequent section for one particular set of stationary solutions to (3.6)). In what follows, these properties of (4.25) and (4.28) are utilised to demonstrate several possible means by which the ZF $v_z'$ may collisionlessly saturate in the system under consideration. Before proceeding to the following section, however, the collisional and collisionless systems are graphically depicted and contrasted as a supplement to the analytic analysis and commentary.

Figures 1 and 2 geometrically compare the fundamental dynamical differences of the limits $\nu =0$ and $\nu \neq 0$ , as applied to the system of three wavevectors $\boldsymbol{k}=(0,k_y)$ , $\boldsymbol{k}'=({-}k_x',0)$ and $\boldsymbol{k}^{\prime \prime }=(k_x',k_y)$ associated with the five modes $\beta _{1}^{}$ , $\beta _{2}^{}$ , $\beta _{1}^{\prime\prime}$ , $\beta _{2}^{\prime\prime}$ and $v_z^{{\prime }}$ described by (2.25)–(2.27). The physical existence of these fluctuations in either limit ultimately arises from the constant temperature gradient $\eta$ , as depicted in the ellipses on the left-hand side of either figure. One function of the gradient is to act as a free-energy source for the encircled unstable modes $\beta _{1}^{}$ and $\beta _{1}^{\prime\prime}$ , as depicted by the solid red arrows. The circles and ellipses are positioned in space so that the reciprocal geometric length of a given red or blue solid arrow alludes to how strongly influenced one of these modes is by the gradient; unstable modes closer to the $\eta$ oval will grow faster under the linear instability than those farther away. This follows because $\boldsymbol{k}$ corresponds to a streamer mode ( $k_x=0$ ) while $\boldsymbol{k}^{\prime \prime }$ is a sideband with $k_x^{\prime \prime } = k_x' \neq 0$ . Hence, for the case $k_y=k_y^{\prime \prime }$ , $\gamma _1 \gt \gamma _1^{\prime\prime}$ (see (2.6)). Therefore, in either figure 1 or 2, the red arrow connecting $\eta$ to $\beta _{1}^{}$ is appreciably shorter than the one linking the former to $\beta _{1}^{\prime\prime}$ . With the $3-\boldsymbol{k}$ system continuously receiving free energy via transference from $\eta$ to the unstable modes, quasi-stationarity is enforced in the nonlinear saturated regime (predominantly) via the linear behaviour of stable modes $\beta _{2}^{}$ and $\beta _{2}^{\prime\prime}$ . These modes continuously deposit their energy back into the gradient at rates $\gamma _2$ and $\gamma _2^{\prime\prime}$ , as depicted by the solid blue arrows of figures 1 and 2. The same geometric length-based convention is used, implying that $\gamma _2 \gt \gamma _2^{\prime\prime}$ , with the linear behaviour of $\beta _{2}^{}$ more strongly influenced by the presence of $\eta$ than that of $\beta _{2}^{\prime\prime}$ . What remains is the question of how energy migrates from locations in the eigenmode space where it is injected to those where it is removed.

Figure 1.

Schematic depiction of the collisionless $3-\boldsymbol{k}$ system. The temperature gradient gives (red arrows) energy to unstable modes $\beta _{1}^{}$ and $\beta _{1}^{\prime\prime}$ while receiving it (blue arrows) from $\beta _{2}^{}$ and $\beta _{2}^{\prime\prime}$ , the ZF $v_z^{{\prime } }$ (yellow lines) is generated by nonlinear interactions among $\beta$ s. The absence of ZF damping results in the interactions being eigenmode-label-permutation symmetric.

Figure 1. Long description

A schematic diagram of the collisionless system. The diagram includes an oval labeled with the Greek letter eta. There are four circles labeled with mathematical expressions: |β₂|², |β₂|², |β₁|², and |β₂|². These circles are connected by arrows of different colors: red, blue, and yellow. The red arrows indicate the flow of energy from the temperature gradient to unstable modes. The blue arrows show the flow of energy from the unstable modes to the zonal flows (ZFs). The yellow arrows represent the nonlinear interactions among the modes. The absence of ZF damping results in eigenmode-label-permutation symmetric interactions.

Figure 2.

Schematic depiction of the collisional $3-\boldsymbol{k}$ system. The temperature gradient gives (red arrows) energy to unstable modes $\beta _{1}^{}$ and $\beta _{1}^{\prime\prime}$ while receiving it (blue arrows) from $\beta _{2}^{}$ and $\beta _{2}^{\prime\prime}$ , the ZF $v_z^{{\prime }}$ (yellow lines) is generated by nonlinear interactions among $\beta$ s. All modes linearly damp energy irreversibly in proportion to $\nu$ , as indicated by the waste receptacle. Non-zero $\nu$ breaks eigenmode-label-permutation symmetry in the strength of the nonlinear interactions accessible to the system.

The primary saturation mechanism of the toroidal ITG instability for the system at hand is described by ZF-catalysed resonant energy transfer from unstable modes to stable modes (Makwana et al. Reference Makwana, Terry and Kim2012, Reference Makwana, Terry, Pueschel and Hatch2014; Terry et al. Reference Terry, Faber, Hegna, Mirnov, Pueschel and Whelan2018; Li & Terry Reference Li and Terry2022). The ZFs are generated and sustained by the turbulence in an entirely nonlinear manner. To lowest non-trivial order in perturbation, the mode $v_z'$ is responsible for facilitating all three-wave interactions between pairs of $\beta _{l}^{}$ and $\beta _{j}^{\prime\prime}$ in the truncated system (Li & Terry 2022, Reference Li and Terry2024) described by (2.25)–(2.27). The aforementioned ideas are characterised by the solid yellow lines in figures 1 and 2; energy may flow through the ZF between the other eigenmodes in either ‘direction’. Note that all four yellow lines correspond to the same single mode $v_z'$ . The lengths of these lines are not intended to directly indicate the magnitude of nonlinear interaction strength, but rather to contrast the amount of symmetry present in the nonlinear ZF dynamics suggested by structures of the tensors $\boldsymbol{T}_{izm}^{k k^{\prime } k^{\prime \prime }}$ and $\boldsymbol{T}_{jzm}^{k^{\prime \prime } -k^{\prime } k (*)}$ in the collisional and collisionless limits.

In the collisionless system of figure 1, the growth and decay rates of unstable and stable modes are permutation symmetric, i.e. $\gamma _1=\gamma _2$ and $\gamma _1^{\prime\prime}=\gamma _2^{\prime\prime}$ . It was shown by the tensorial identities of (4.23)–(4.27) that the aforementioned symmetry results in the constants of proportionality governing nonlinear ZF time evolution possessing their own corresponding symmetries and antisymmetries. This idea manifests in figure 1 as follows. Given that $\gamma _1=\gamma _2$ , $\beta _{1}^{}$ and $\beta _{2}^{}$ are equidistant from the $\eta$ oval. A similar argument follows for the $k^{\prime \prime }$ fluctuations. Then, the nonlinearity, which has been shown to merely be arithmetic combinations of the individual modal linearities, appears as two pairs of permutation-symmetric yellow lines. In the collisionless diagram, the yellow line that connects $\beta _{1}^{}$ to $\beta _{2}^{\prime\prime}$ may be exchanged with the line linking $\beta _{2}^{}$ with $\beta _{1}^{\prime\prime}$ , and the diagram is invariant. The same argument holds true for the line-pair combinations $\beta _{1}^{},\beta _{1}^{\prime\prime}$ and $\beta _{2}^{},\beta _{2}^{\prime\prime}$ . This is a graphical analogue of the centrosymmetries symmetries of the nonlinear time-evolution operator as suggested by (4.25) and (4.28). As is already known from a mathematical perspective, the situation is fundamentally different for the collisional system.

For the case $\nu \neq 0$ , the linear mode frequency given by (2.6) suggests that all five eigenmodes of the system will exhibit linear damping due to $\nu$ . This idea is depicted by the dotted green arrows of figure 2; heat can be directly and irreversibly removed from any mode of the system to be externally dissipated and disposed of (as indicated by the waste receptacle), as opposed to always eventually being returned back to the gradient with the potential for being recycled by the system. The presence of $\nu$ breaks symmetry in the mode frequency such that stable-mode damping is stronger than unstable-mode growth for both pairs of $k_y\neq 0$ modes, i.e. $\gamma _1\lt \gamma _2$ and $\gamma _1^{\prime\prime}\lt \gamma _2^{\prime\prime}$ . This is depicted in the figure by $\beta _{2}^{}$ and $\beta _{2}^{\prime\prime}$ being closer in proximity to $\eta$ than their corresponding conjugate unstable modes. Broken symmetry in the linear dynamics due to $\nu$ lends to corresponding features in the nonlinear dynamics as suggested by the collisional time-evolution tensor of (B23)–(B33). In general, for $\nu \neq 0$ , the elements of a given tensor are not equal in magnitude to their index-cyclically permuted counterparts, i.e. the identities of (4.23)–(4.27) do not hold true. Hence, the diagram is not invariant with respect to the exchange of the yellow lines connecting $\beta _{1}^{}$ to $\beta _{2}^{\prime\prime}$ and $\beta _{2}^{}$ to $\beta _{1}^{\prime\prime}$ , or those linking $\beta _{1}^{}$ to $\beta _{1}^{\prime\prime}$ and $\beta _{2}^{}$ to $\beta _{2}^{\prime\prime}$ . (Geometrically, the yellow arrowheads would no longer terminate on modal circles under such an operation.) The constants of proportionality in the nonlinear ZF time-evolution operator are all different in magnitude from one another in the limit $\nu \neq 0$ .

5. Saturation theory and phenomenology

In this section, stationary solutions of the ZF WKE (3.6) are found from the preceding mathematical analysis, and a phenomenology is developed in connection to the thermal energy transport associated with these systems. We proceed by reviewing previous results for non-zero collisionality in the interest of making contact with pre-existing theory (Dimits et al. Reference Dimits2000b ).

5.1. Collisional systems

In general, for non-zero collisionality, the saturation condition (3.9) is given by

(5.1) \begin{align} \nu &= \kappa \nonumber \\ &= \mathrm{Re} \big[ C_{zij}^{({k^{\prime }},{k},{-k^{\prime \prime }})} \big( C_{izm}^{({k},{k^{\prime }},{k^{\prime \prime }})} \tau _{izm}^{kk^{\prime }k^{\prime \prime }}\langle \beta _{j}^{{\prime\prime}*} \beta _{m}^{\prime\prime} \rangle + C_{jzm}^{({k^{\prime \prime }},{-k^{\prime }},{k})*} \tau _{jzm}^{k^{\prime \prime }k^{\prime }k}\langle \beta _{m}^{*}\beta _{i}^{} \rangle \big ) \big]\nonumber \\[5pt] &= \mathrm{Re} \big[ C_{zij}^{({k^{\prime }},{k},{-k^{\prime \prime }})} \big( \boldsymbol{T}_{izm}^{k, k^{\prime }, k^{\prime \prime }}\langle \beta _{j}^{{\prime\prime}*} \beta _{m}^{\prime\prime} \rangle + \boldsymbol{T}_{jzm}^{k^{\prime \prime }, -k^{\prime }, k (*)}\langle \beta _{m}^{*}\beta _{i}^{} \rangle \big ) \big]\nonumber \\[5pt] &= C_{zij}^{({k^{\prime }},{k},{-k^{\prime \prime }})} \big( T_{I,izm}^{k, k^{\prime }, k^{\prime \prime }}\langle \beta _{j}^{{\prime\prime}*} \beta _{m}^{\prime\prime} \rangle _R+T_{R,izm}^{k, k^{\prime }, k^{\prime \prime }}\langle \beta _{j}^{{\prime\prime}*} \beta _{m}^{\prime\prime} \rangle _I + T_{I,jzm}^{k^{\prime \prime }, -k^{\prime }, k (*)}\langle \beta _{m}^{*}\beta _{i}^{} \rangle _R \nonumber \\[5pt] &\quad + T_{R,jzm}^{k^{\prime \prime }, -k^{\prime }, k (*)}\langle \beta _{m}^{*}\beta _{i}^{} \rangle _I \big) , \end{align}

and the ZF damping is balanced by the overall nonlinear drive of the ZF. The above may be further expanded, investigated and reduced, but there are already known results (Terry et al. Reference Terry, Faber, Hegna, Mirnov, Pueschel and Whelan2018; Li & Terry Reference Li and Terry2022) for this system suggesting that the ZF amplitude is set by the growth rate of the linear instability. From Terry et al. (Reference Terry, Faber, Hegna, Mirnov, Pueschel and Whelan2018),

(5.2) \begin{align} |v_z'|^2 = \frac {4.08 \gamma _1}{\mathrm{Re}(\tau _{1z2})|C_{izj}|^2}, \end{align}

while $\nu$ sets the level of the turbulent energy-like quantity $E_T$ , i.e.

(5.3) \begin{align} E_T = \frac {\nu }{|C_{z12}||C_{1z2}|\mathrm{Re}(\tau _{1z2})}, \end{align}

where $\tau _{1z2}$ is the TCT

(5.4) \begin{align} \tau _{1z2} = \frac {1}{\mathrm{i}(\omega _2^{\prime \prime } + \omega _z - \omega _1^*)}. \end{align}

In addition, it has been shown that, for the $3-\boldsymbol{k}$ system, the ratio of stable to unstable energy spectra in the saturated state is (Li & Terry Reference Li and Terry2022)

(5.5) \begin{align} \frac {\langle \beta _{2}^{2} \rangle }{\langle \beta _{1}^{2} \rangle } = \frac {\gamma _1}{\gamma _2} + O(\nu ). \end{align}

Further evaluation of (5.1) in terms of its indices is not particularly insightful or relevant to the goal of the current work; the algebra is cumbersome because the tensor identities of (4.23)–(4.27) do not have the simplifications possible in the limit $\nu =0$ . The primary point regarding (5.1) is that in the presence of collisionality, a combination of $\nu$ damping and transfer to stable modes balances the ZF, which is nonlinearly driven by transfer from the unstable modes. In the following subsection, it is shown that, for the collisionless case, the saturation theory is a limit of the aforementioned argument; a clear balance of nonlinear drive from unstable modes and nonlinear damping via energy transfer to stable modes brings about saturation. Importantly, the collisional saturation balance is in general more mathematically complicated because of the presence of $\nu$ -related terms in the tensor (B23)–(B33) that break the symmetry embodied in (4.23)–(4.27). The broken symmetry plays a significant role in the resulting thermal flux $Q$ .

The thermal flux is defined as

(5.6) \begin{align} Q = - \sum _k k_y \mathrm{Im} \langle \phi _{-k}p_k\rangle , \end{align}

or given in the eigenmode representation by

(5.7) \begin{align} Q = -\sum _k k_y \left [\mathrm{Im} R_1 \langle \beta _{1}^{2} \rangle + \mathrm{Im}R_2 \langle \beta _{2}^{2} \rangle + \mathrm{Im}(R_1+R_2)\mathrm{Re}\left \langle \beta _{1}^{}\beta _{2}^{*}\right \rangle \right ]. \end{align}

For collisional systems, it is known from gyrokinetic simulations (Lin et al. Reference Lin, Hahm, Lee, Tang and Diamond1999) and prior analysis (Terry et al. Reference Terry, Faber, Hegna, Mirnov, Pueschel and Whelan2018; Li & Terry Reference Li and Terry2022) of the reduced two-field fluid model of (2.19)–(2.21) that $Q$ scales in proportion to $\nu$ . In particular, it has been shown (Terry et al. Reference Terry, Faber, Hegna, Mirnov, Pueschel and Whelan2018) that

(5.8) \begin{align} Q \propto \sum _k \frac {\gamma _1}{k_\perp ^2} \frac {\nu }{\mathrm{Re}(\tau _{1z2})k_y^2\epsilon _n^{1/2}}, \end{align}

by which $\nu$ appears explicitly in the numerator but also resides with the TCT $\tau _{1z2}$ . In general, collisional stationary states imply that $Q\gt 0$ . This can be attributed to unstable modes on average saturating at amplitudes above conjugate stable modes in proportion to $\nu$ , as described be (5.5).

5.2. Collisionless systems

We consider application of $\nu =0$ to the results presented thus far. This case is explicitly evaluated to reveal the collisionless system’s saturation mechanism and accompanying phenomenology that arise due to the underlying dynamical symmetries. The exact means through which these properties arise are otherwise unclear if one simply applies $\nu =0$ to the prior results of (5.5) and (5.8).

We begin analysis of the collisionless system by solving the pair of saturation conditions described by (4.7) and (4.8). The equalities $\boldsymbol{T}_{izm}^{k k^{\prime } k^{\prime \prime }} \langle \beta _{j}^{{\prime\prime}*} \beta _{m}^{\prime\prime}\rangle = 0$ and $\boldsymbol{T}_{jzm}^{k^{\prime \prime } -k^{\prime } k (*)} \langle \beta _{m}^{*}\beta _{i}^{} \rangle = 0$ , which occur simultaneously as formally independent constraints, only depend on energy and cross-correlation at a single wavenumber $k$ or $k^{\prime \prime }$ . As such, (4.7) and (4.8) are referred to as monochromatic saturation (MS) conditions. Recalling the tensor-correlation expansion of (4.6), the $k^{\prime \prime }$ saturation condition of (4.7) is expanded as

(5.9) \begin{align} \boldsymbol{T}_{izm}^{k, k^{\prime }, k^{\prime \prime }} \langle \beta _{j}^{{\prime\prime}*} \beta _{m}^{\prime\prime}\rangle = T_{I,{i}{z}{m}}^{{k}{k^{\prime }}{k^{\prime \prime }}}\langle \beta _{j}^{{\prime\prime}*} \beta _{m}^{\prime\prime} \rangle _R + T_{R,{i}{z}{m}}^{{k}{k^{\prime }}{k^{\prime \prime }}}\langle \beta _{j}^{{\prime\prime}*} \beta _{m}^{\prime\prime} \rangle _I =0. \end{align}

It is not hard to show (see Appendix C), that this condition reduces to

(5.10) \begin{align} \boldsymbol{T}_{izm}^{k, k^{\prime }, k^{\prime \prime }} \langle \beta _{j}^{{\prime\prime}*} \beta _{m}^{\prime\prime}\rangle = \big(T_{I,{1}{z}{1}}^{{k}{k^{\prime }}{k^{\prime \prime }}}+T_{I,{2}{z}{1}}^{{k}{k^{\prime }}{k^{\prime \prime }}}\big)\big(\langle \beta _{1}^{{\prime\prime}2} \rangle - \langle \beta _{2}^{{\prime\prime}2} \rangle \big) =0, \end{align}

which implies that $\langle \beta _{1}^{{\prime\prime}2} \rangle = \langle \beta _{2}^{{\prime\prime}2} \rangle$ . Similarly, the MS condition of (4.8) for fluctuations with wavevector $\boldsymbol{k}$ may be expanded as

(5.11) \begin{align} \boldsymbol{T}_{jzm}^{k^{\prime \prime }, -k^{\prime }, k (*)} \langle \beta _{m}^{*}\beta _{i}^{} \rangle = T_{I,{j}{z}{m}}^{{k^{\prime \prime }}{-k^{\prime }}{k(*)}}\langle \beta _{m}^{ *} \beta _{i}^{} \rangle _R + T_{R,{j}{z}{m}}^{{k^{\prime \prime }}{-k^{\prime }}{k(*)}}\langle \beta _{m}^{ *} \beta _{i}^{} \rangle _I =0, \end{align}

and by the same reasoning that simplified the $k^{\prime \prime }$ MS condition reduces to (see Appendix C)

(5.12) \begin{align} \boldsymbol{T}_{jzm}^{k^{\prime \prime }, -k^{\prime }, k (*)} \langle \beta _{m}^{*}\beta _{i}^{} \rangle = \big(T_{I,{1}{z}{1}}^{{k^{\prime \prime }}{-k^{\prime }}{k(*)}}+T_{I,{2}{z}{1}}^{{k^{\prime \prime }}{-k^{\prime }}{k(*)}}\big)\big(\langle \beta _{1}^{2} \rangle -\langle \beta _{2}^{2} \rangle \big) =0 . \end{align}

This implies that $\langle \beta _{1}^{2} \rangle = \langle \beta _{2}^{2} \rangle$ . As such, the ZF of the $3-\boldsymbol{k}$ system of (2.25)–(2.27) with WKE (3.6) may saturate monochromatically when the energies of the unstable and conjugate stable mode at each wavenumber are equal. Mathematically, this corresponds to

(5.13) \begin{align} \langle \beta _{1}^{2} \rangle &= \langle \beta _{2}^{2} \rangle ,\\[-10pt]\nonumber\end{align}

(5.14) \begin{align} \langle \beta _{1}^{{\prime\prime}2} \rangle &= \langle \beta _{2}^{{\prime\prime}2} \rangle . \end{align}

Equations (5.13) and (5.14) coincide with the limit $\nu \rightarrow 0$ applied to the stable/unstable spectral ratio of (5.5). These conditions are not possible in the collisional case because the matrix elements multiplying the correlations are in general sufficiently different from one another, i.e. they do not satisfy the centrosymmetry identities of (4.23)–(4.27), and the simplifications in the calculations of Appendix C can not transpire.

At present, the collisionless saturation condition of (5.13) and (5.14) appear to suggest an ambiguity in the total amount of turbulent energy present in the system when $\langle v_z^{{\prime } 2} \rangle$ is fixed. Consider time-evolution equations for stable and unstable modes with wavenumber $k$ and generic nonlinearities of the forms

(5.15) \begin{align} \partial _t \langle \beta _{1}^{2} \rangle - 2 \gamma _1 \langle \beta _{1}^{2} \rangle = NL_1 \end{align}

and

(5.16) \begin{align} \partial _t \langle \beta _{2}^{2} \rangle + 2 \gamma _2 \langle \beta _{2}^{2} \rangle = NL_2. \end{align}

It is not hard to show that the steady-state collisionless behaviour of these modes when the ZF $v_z^{{\prime }}$ is saturated corresponds to the condition

(5.17) \begin{align} NL_1 = - NL_2, \end{align}

for which the nonlinearities are antisymmetric under permutation. Equation (5.17) only constrains the nonlinear energy-transfer rates and not the absolute turbulent energy level. We will address these ambiguities in turbulent saturation level for the collisionless $3-\boldsymbol{k}$ system in a subsequent work through a combination of closure theory and numerical simulation. Application of the limit $\nu = 0$ to the turbulent energy-like quantity given by (5.3) allows for considerable simplification using the permutation symmetries of the nonlinear coupling coefficients. The expression reduces to combinations of the eigenmode energy in the form

(5.18) \begin{align} E_{T,\nu =0} = 2\big(\langle \beta _{1}^{{\prime }2}\rangle - \langle \beta _{2}^{{\prime }2}\rangle + \langle \beta _{2}^{{\prime\prime}2} \rangle - \langle \beta _{1}^{{\prime\prime}2} \rangle \big), \end{align}

from which it is observed that $E_{T,\nu =0} =0$ for a state of MS. This is in agreement with (5.3).

Now we briefly consider the more general collisionless saturation condition described by (4.11). Using (5.10) and (5.12), (4.11) may be expanded to yield

(5.19) \begin{align} &\boldsymbol{T}_{izm}^{k, k^{\prime }, k^{\prime \prime }} \langle \beta _{j}^{{\prime\prime}*} \beta _{m}^{\prime\prime}\rangle + \boldsymbol{T}_{jzm}^{k^{\prime \prime }, -k^{\prime }, k (*)} \langle \beta _{m}^{*} \beta _{i}^{}\rangle \nonumber \\[3pt] &\quad =\big(T_{I,{1}{z}{1}}^{{k}{k^{\prime }}{k^{\prime \prime }}}+T_{I,{2}{z}{1}}^{{k}{k^{\prime }}{k^{\prime \prime }}}\big)\big(\langle \beta _{1}^{{\prime\prime}2} \rangle - \langle \beta _{2}^{{\prime\prime}2} \rangle \big)\\[-10pt]\nonumber \end{align}

(5.20) \begin{align} & \quad +\big(T_{I,{1}{z}{1}}^{{k^{\prime \prime }}{-k^{\prime }}{k(*)}}+T_{I,{2}{z}{1}}^{{k^{\prime \prime }}{-k^{\prime }}{k(*)}}\big)\big(\langle \beta _{1}^{2} \rangle - \langle \beta _{2}^{2} \rangle \big) =0, \end{align}

implying stationary solutions of the ZF WKE (3.6) described by the set of states satisfying the constraint

(5.21) \begin{align} \frac {\big(\langle \beta _{1}^{2} \rangle - \langle \beta _{2}^{2} \rangle \big)}{\big(\langle \beta _{1}^{{\prime\prime}2} \rangle - \langle \beta _{2}^{{\prime\prime}2} \rangle \big)} = -\frac {\big(T_{I,{1}{z}{1}}^{{k}{k^{\prime }}{k^{\prime \prime }}}+T_{I,{2}{z}{1}}^{{k}{k^{\prime }}{k^{\prime \prime }}}\big)}{\big(T_{I,{1}{z}{1}}^{{k^{\prime \prime }}{-k^{\prime }}{k(*)}}+T_{I,{2}{z}{1}}^{{k^{\prime \prime }}{-k^{\prime }}{k(*)}}\big)}. \end{align}

This saturation balance is characterised by a ratio of spectral symmetry broken in proportion to the ratio of the nonlinear time-evolution tensors. The tensor sums $T_{I,{1}{z}{1}}^{{k}{k^{\prime }}{k^{\prime \prime }}}+T_{I,{2}{z}{1}}^{{k}{k^{\prime }}{k^{\prime \prime }}}$ and $T_{I,{1}{z}{1}}^{{k^{\prime \prime }}{-k^{\prime }}{k(*)}}+T_{I,{2}{z}{1}}^{{k^{\prime \prime }}{-k^{\prime }}{k(*)}}$ are in general non-trivial and do not appear to reduce in a simple transparent manner algebraically. It can be shown in the collisionless limit that all second-order correlations of eigenmodes ( $\beta 's$ ) in this system will saturate under the MS conditions. In contrast, the asymmetric solution suggested by (5.21) only corresponds to $\partial _t \langle v_z^2\rangle =0$ , but does not yield a steady state for the eigenmodes. We present the relevant saturation calculations for the eigenmodes in a subsequent work (Azelis et al. Reference Azelis, Terry, Sprott and Li2026) and also confirm observations of MS solutions in a numerical simulation of the system. As such, the remainder of this analysis is focused on the symmetric solutions suggested by the MS condition of (5.13) and (5.14), which are a subset of the general case. We proceed by considering the collisionless thermal flux.

In the collisionless limit, the unstable/stable eigenvector expression of (2.5) becomes

(5.22) \begin{align} R_j = \frac {\omega _j(1+k_\perp ^2)-k_y}{ k_y \epsilon _n}, \end{align}

with imaginary part

(5.23) \begin{align} \mathrm{Im} R_j = \frac {\pm \gamma (1+k_\perp ^2)}{k_y \epsilon _n}. \end{align}

Therefore,

(5.24) \begin{align} \mathrm{Im} R_2 = - \mathrm{Im} R_1 \end{align}

and

(5.25) \begin{align} \mathrm{Im}(R_1 +R_2) = 0. \end{align}

The collisionless heat flux is then given by

(5.26) \begin{align} Q = -\sum _k k_y \mathrm{Im} R_1 \big( \langle \beta _{1}^{2} \rangle -\langle \beta _{2}^{2} \rangle \big) = \sum _k \frac {\gamma (1+k_\perp ^2)}{ \epsilon _n} \big( \langle \beta _{1}^{2} \rangle -\langle \beta _{2}^{2} \rangle \big), \end{align}

which after expanding the sum for the $3-\boldsymbol{k}$ system becomes

(5.27) \begin{align} Q = \frac {1}{\epsilon _n}\big[\gamma (1+k_y^2) \big( \langle \beta _{1}^{2} \rangle -\langle \beta _{2}^{2} \rangle \big) + \gamma ^{\prime \prime }(1+k_x^{{\prime }2} +k_y^2) \big( \langle \beta _{1}^{{\prime\prime}2} \rangle -\langle \beta _{2}^{{\prime\prime}2} \rangle \big)\big]. \end{align}

The terms in (5.27) clearly indicate that energy in an unstable mode corresponds to positive $Q$ and down-gradient transport, whereas the converse is true for energy in a stable mode. The MS condition of (5.13) and (5.14) then implies turbulence characterised by perfectly balanced up/down-gradient transport over each individual fluctuation length scale ( $k^{-1}$ ), i.e. states of zero net transport despite a finite fluctuation level. This matches the limiting behaviour of the collisional thermal flux (5.8) for $\nu \rightarrow 0$ . Remarkably, a saturated ZF with zero collisionality appears to suggest ideal confinement in the presence of the ITG instability, given this system. Interestingly, this is in contradiction with quasilinear theories of transport, which having been developed from $Q \propto \gamma /k_\perp ^2$ with no explicit accounting for stable modes, predict a finite flux whenever an instability is present (Kotschenreuther et al. Reference Kotschenreuther, Dorland, Beer and Hammett1995; Jenko, Dannert & Angioni Reference Jenko, Dannert and Angioni2005). This result raises the question of whether quasilinear models, such as Pueschel et al. (Reference Pueschel, Faber, Citrin, Hegna, Terry and Hatch2016), can be adapted to account for stable mode damping.

The result $Q=0$ for $\nu =0$ can be predicted without explicit consideration of ZF dynamics or the eigenmode decomposition utilised herein. One way to demonstrate this is by multiplying (2.2) by $\phi _k^*$ and performing a sum over all $\boldsymbol{k}$ to produce the expression

(5.28) \begin{align} \sum _{\boldsymbol{k}} [\delta (k_y)+k_\perp ^2]\frac {\partial |\phi _k|^2}{\partial t} - \epsilon _n Q + \sum _{\boldsymbol{k}} \nu k_\perp ^2 |\phi _k|^2 =0. \end{align}

We note that the nonlinearity vanishes via integration of an odd function over symmetric bounds. For $\nu =0$ , steady-state behaviour of $|\phi _k|^2$ corresponds to $Q=0$ , and this conclusion is reached without the lengthy calculations presented herein. However, when presented in this context, the property $Q=0$ merely takes the form of a constraint imposed upon the motions of the fluid, with little indication of the statistical properties of said motions or any significant phenomenology in terms of the fields being modelled. In the absence of the eigenmode decomposition, the constraint simply suggests dynamics that on average satisfy

(5.29) \begin{align} \sum _k k_y \mathrm{Im} \langle \phi _{-k}p_k\rangle = 0. \end{align}

While this does indeed suggest on average that the overall thermal energy being transported across all fluctuation length scales in the $+\hat {\boldsymbol{x}}$ is cancelled by fluxes in the $-\hat {\boldsymbol{x}}$ , it is not clear at all how the system achieves such a state and provides minimal description of the governing physics. Furthermore, there is no indication of whether the system has any particular preference for how the zero sum is achieved. It is an ambiguous linear combination of cross-correlations for which there are infinitely many possible $\phi _k$ and $p_k$ products that may satisfy it. Our result suggests that the system has a particular preference for equipartition of energy between eigenmodes at each fluctuation length scale $k^{-1}$ and demonstrates this within the context of the underlying turbulent saturation theory physics.

6. Generalisation to the many $-\boldsymbol{k}$ problem

In can be shown that the ZF of the many $-\boldsymbol{k}$ system described by (2.19)–(2.21) has a WKE at lowest non-trivial order in perturbation described by

(6.1) \begin{align} \partial _t \langle v_z^{2}\rangle = -2\left [\nu - K(k) \right ]\langle v_z^{2}\rangle , \end{align}

where

(6.2) \begin{align} K(k) = \mathrm{Re} \biggl [\sum _{k^{\prime }} C_{zmn}^{({k},{k^{\prime }},{k^{\prime \prime }})} \biggl (&\frac {C_{mzp}^{({k^{\prime }},{k},{-k^{\prime \prime }})}\langle \beta _{p}^{{\prime\prime}*}\beta _{n}^{\prime\prime}\rangle }{\mathrm{i}(\omega _{m}^{{\prime } }-\omega _{z}^{}+\omega _{p}^{{\prime\prime}*})} +\frac {C_{mqz}^{({k^{\prime }},{-k^{\prime \prime }},{k})}\langle \beta _{q}^{{\prime\prime}*}\beta _{n}^{\prime\prime}\rangle }{\mathrm{i}(\omega _{m}^{{\prime } }+\omega _{q}^{{\prime\prime}*}-\omega _{z}^{})} \nonumber \\[4pt] &+ \frac {C_{nzp}^{({k^{\prime \prime }},{k},{-k^{\prime }})}\langle \beta _{p}^{{\prime }*}\beta _{m}^{{\prime } }\rangle }{\mathrm{i}(\omega _{n}^{\prime\prime}-\omega _{z}^{}+\omega _{p}^{{\prime }*})} +\frac {C_{nqz}^{({k^{\prime \prime }},{-k^{\prime }},{k})}\langle \beta _{q}^{{\prime }*}\beta _{m}^{{\prime } }\rangle }{\mathrm{i}(\omega _{n}^{\prime\prime}+\omega _{q}^{{\prime }*}-\omega _{z}^{})}\biggl ) \biggl ]. \end{align}

The solution is given by

(6.3) \begin{align} \langle v_z^{2}\rangle (k_x,t) = \langle v_z^{2}\rangle _{t=0} \textrm {e}^{-2[\nu - K(k)]t}. \end{align}

The calculation is similar to Appendix A but with several subtle differences. One in particular is associated with the fourth-order moments not automatically being closed, but this is handled by the perturbation theory (Malfliet Reference Malfliet1972) with ease. This derivation is outside the scope of the current paper and would require many more pages, hence, it is omitted. However, it is not hard to show the result by following the methods of Malfliet (Reference Malfliet1972) in combination with the gyrokinetics-informed perturbative ordering scheme of Terry et al. (Reference Terry, Faber, Hegna, Mirnov, Pueschel and Whelan2018).

Striking similarities are observed between $K$ and $\kappa$ , with the key differences being the sum $\text{d}\boldsymbol{k}'$ over fluctuations and the fact that there are four terms instead of two. Nearly identical analysis can show nonetheless that the WKE of (6.1) also exhibits MS solutions in the collisionless limit with comparable up/down-gradient transport cancellations leading to the same saturation phenomenology. This is not surprising as the systems of equations (2.19)–(2.21) and (2.25)–(2.27) are akin, and it has previously been demonstrated that the $3-\boldsymbol{k}$ system reproduces similar saturation physics as the many $-\boldsymbol{k}$ system (Li & Terry Reference Li and Terry2022). The primary differences between the systems arise from the vastly increased number of dynamical degrees of freedom in the many $-\boldsymbol{k}$ system. In particular, this system does not impose any explicit constraints on the kurtosis. Moreover, the absence of wavenumber truncation enables the development of a full turbulence cascade, though its character differs from the standard forward-Kolmogorov picture in the collisionless limit.

For $\chi = \nu = 0$ , the imaginary portion $\gamma _j$ of the complex mode frequency (2.6) is finite and positive within a contour in the $k_xk_y$ plane satisfying (2.8). Outside this contour, $\gamma _j=0$ and all modes are marginally stable, exhibiting only oscillatory behaviour. (In contrast, for finite collisionality, $\gamma _j$ remains finite but asymptotically decays to zero with increasing wavenumber, i.e. $\lim _{k\rightarrow \infty } \gamma _j = 0$ .) Hence, in the collisionless limit, only a fraction of the modes are capable of either injecting or removing energy from the system, and this behaviour only occurs in a bounded and connected region of the Fourier domain. Saturation of the full wavenumber system therefore requires a mechanism by which energy, nonlinearly transferred via forward cascade to wavenumbers beyond the contour in (2.8), is eventually returned to the stable modes within that contour, where it can then be removed.

It has been observed that the plasma wave echo phenomenon (Gould, O’Neil & Malmberg Reference Gould, O’Neil and Malmberg1967; Malmberg et al. Reference Malmberg, Wharton, Gould and O’Neil1968) can inversely cascade energy from high to low wavenumbers in long-wavelength electrostatic turbulence in magnetised, weakly collisional plasmas, including (but not limited to) ITG drift-wave turbulence (Schekochihin et al. Reference Schekochihin, Parker, Highcock, Dellar, Dorland and Hammett2016). It has also been demonstrated that in trapped electron mode turbulence, where saturation occurs via energy transfer to stable modes (Terry et al. Reference Terry, Gatto and Baver2002), wavenumber anisotropy gives rise to inverse energy transfer that condenses onto zonal modes (Terry Reference Terry2004). However, a full investigation of energy-transfer mechanisms in the collisionless many $-\boldsymbol{k}$ system lies beyond the scope of this study and is left for future work.