A linearly unstable, sinusoidal *E *×* B* shear flow is examined in the gyrokinetic framework in both the linear and nonlinear regimes. In the linear regime, it is shown that the eigenmode spectrum is nearly identical to hydrodynamic shear flows, with a conjugate stable mode found at every unstable wavenumber. In the nonlinear regime, turbulent saturation of the instability is examined with and without the inclusion of a driving term that prevents nonlinear flattening of the mean flow and a scale-independent radiative damping term that suppresses the excitation of conjugate stable modes. From a variety of analyses, the nonlinear state is found to have a significant component associated with stable modes. The role of these modes is investigated through a simple fluid model that tracks how momentum transport and partial flattening of the mean flow scale with the driving term. From this model, it is shown that, except at high radiative damping, stable modes play an important role in the turbulent state and yield significantly improved quantitative predictions when compared with corresponding models neglecting stable modes.

## I. INTRODUCTION

The prevalence of sheared flows in diverse systems has motivated their study for over a century. Their potential to drive instabilities and turbulence in fluids and plasmas is central to angular momentum transport in astrophysical disks,^{1,2} to transport in the Earth's magnetosphere,^{3} and possibly to the generation and saturation of confinement-modifying zonal flows in fusion devices.^{4} The linear stability of simple shear flow configurations has been thoroughly investigated from linear equations^{5,6} and provides a rough understanding of the nature of more complex flow profiles early in their development, before unstable perturbations grow too large.^{7,8} However, as these flows develop beyond the regime of validity of linearized models, and nonlinear interactions between different components of the fluctuation become important, capturing or understanding their behavior with any set of constructs based on linear analysis becomes problematic.^{9–11}

Instead, studies generally rely on direct numerical simulations to investigate relevant physical effects.^{3,8,12} In many cases of interest, these methods cannot produce solutions for physically relevant parameters, such as the high Reynolds numbers found in astrophysical systems. This motivates the development of scaling models that can inform how the system extrapolates to parameter regimes inaccessible to simulations. Valid scaling models require an understanding of the physics of all relevant phenomena, including turbulent responses that modify the unstable flow, like nonlinear fluctuation structures, cascades, and momentum transport.

Regarding nonlinear processes that become relevant as the linear growth phase ends, recent analytical work on shear-flow instability saturation has demonstrated the importance of fluctuation dissipation that arises at large scales due to excitation of stable modes.^{13} When an unstable shear flow is perturbed from equilibrium, these linear modes are generally a part of the initial perturbation, decaying from their small initial amplitude. Given this initial decay, stable modes are typically ignored in constructing reduced nonlinear models that draw from linear physics.^{9–11} However, nonlinear interactions with unstable modes can drive stable modes to large amplitude. Because they are linearly stable, they provide a route for energy to be removed from fluctuations at large scales, before it is able to cascade to small scales, thereby modifying the flow, its spectrum, and its transport.^{13} This represents a significant departure from the usual picture of instability-driven turbulence, where energy injection by unstable modes is assumed to be balanced by conservative nonlinear energy transfer to small, dissipative scales.

While it has been shown that nonlinear interactions with large-scale stable modes can be important in saturating shear-flow instabilities, their amplitude and contribution to the fluctuating flow and momentum transport in fully developed turbulence remains an open question, which we pursue in this paper. Additionally, we explore whether reduced models of shear-flow-driven turbulence that are based solely on the linear instability might be improved by including the effects of large-scale stable modes. This is a natural expectation given their importance in saturating the instability, their introduction of a large-scale linear energy sink, and their potential to modify momentum transport. This is also motivated by recent work in the context of instability-driven turbulence in fusion devices, where reduced turbulence models that include details of stable modes and instability saturation physics have been shown to be effective.^{14–16}

We address these questions by performing direct numerical simulations of an unstable shear flow that develops into turbulence and comparing the contribution of different linear modes to the turbulent flow and the Reynolds stress. Our simulations are performed using the gyrokinetic turbulence code Gene,^{17,18} which has previously been used to examine stable modes in other turbulent systems,^{19,20} and includes both initial value and eigenvalue solvers. This allows us to benchmark our calculated growth rates against previous gyrokinetic studies of the same system,^{21} as well as investigate differences between shear flow instabilities in hydrodynamics and gyrokinetics with regard to both the linear mode spectrum and instability saturation. In particular, while it is understood that all unstable, inviscid, incompressible, two-dimensional (2D) hydrodynamic flows include one stable mode for every unstable mode,^{6} and previous work has shown that these stable modes are nonlinearly driven in the fluid system,^{13} whether these results apply to the gyrokinetic case as well has not been explored. To allow for more direct comparisons with previous work, all simulations presented in this paper are effectively 2D, with no variations in the direction of the strong guide field (*k _{z}* = 0).

The flow we examine is a sinusoidally varying *E *×* B* parallel shear flow with periodic boundary conditions. The hydrodynamic counterpart to this flow is often referred to as Kolmogorov flow when it is maintained by a constant forcing term.^{22–24} This flow profile is particularly relevant to astrophysical disks, where its Kelvin-Helmholtz (KH) instability is studied as a saturation mechanism for the magnetorotational instability^{25–30} or its collisionless counterpart,^{31} and in fusion devices, where it is studied as a potential secondary and tertiary instability to streamers and zonal flows.^{4,32} In order to admit a quasi-stationary state of driven turbulence where energy dissipation is balanced by energy injection, we continually reinforce the mean flow using a Krook operator previously employed similarly to reinforce current gradients in tearing mode studies,^{33} and referred to as a linear relaxation term in studies of barotropic jets.^{34} With this forcing term, the system bears a strong resemblance to Kolmogorov flow,^{22–24} with the exception that it is not a constant forcing. From a numerical perspective, Kolmogorov flow presents a convenient choice of unstable shear flow to study due to its simple description in a Fourier basis and the lack of no-slip boundary conditions that could otherwise generate boundary layers. This also allows us to address whether the saturation physics active in the free shear layer^{13} is applicable to a driven periodic shear flow.

Our simulations also include damping terms in the form of hyperdissipation and scale-independent radiative damping. The form of the radiative damping term is such that it damps every mode equally. In systems with pairs of stable and unstable modes, this disproportionately affects the stable mode amplitude relative to the unstable one in the nonlinear state.^{35} Thus, varying the degree of radiative damping in our system allows us to assess whether different shear-driven turbulence regimes exist with significantly different stable mode effects, and how these regimes might differ.

The remainder of this paper is organized as follows: Section II starts with a brief review of hydrodynamic parallel shear flows for comparison with our gyrokinetic results, as well as some unique aspects of the particular flow profile studied here, followed by a discussion of the numerical implementation used in our work, including the specific forms of forcing and dissipation. In Sec. III, we show the full eigenmode spectrum for the gyrokinetic KH instability. A description of the nonlinear evolution of the flow is presented in Sec. IV, where we discuss saturation and decaying turbulence when forcing is absent, driven turbulence with external forcing, and turbulent momentum transport in this system. Section V examines the turbulence in terms of the role played by the linear eigenmodes and compares reduced descriptions and scaling models of the turbulence with and without stable modes. Conclusions are presented in Sec. VI.

Throughout this paper, we adopt the notation that $A\u0302(x,ky)$ denotes the Fourier transform in *y* of *A*(*x*, *y*), and $A\u0303(kx,ky)$ denotes the Fourier transform in *x* and *y*.

## II. SHEAR FLOW INSTABILITY

### A. Rayleigh's stability equation

The stability of parallel shear flows is generally investigated by examining infinitesimal perturbations about equilibrium solutions to the Navier-Stokes equation. When considering a 2D, inviscid, incompressible flow that is perturbed from an equilibrium, the vorticity equation becomes

where *V*(*x*) is the *y*-directed equilibrium shear flow, and $\varphi (x,y,t)$ is the streamfunction of the perturbation $v=\u2207\varphi \xd7z\u0302$. The linear dynamics can then be explored by dropping the nonlinearities and using the normal mode ansatz

yielding

Equation (3) is known as Rayleigh's stability equation, or as the Orr-Sommerfeld equation when the effect of viscosity on *ϕ* is included. It can be solved as an eigenvalue problem, yielding a set of eigenvalues *ω _{j}* and eigenmodes $\varphi \u0302j$, with

*j*enumerating the eigenmodes at a given

*k*. The eigenvalue

_{y}*ω*is complex, with real frequency $Re(\omega j)$ and growth rate $\gamma j=\u2212Im(\omega j)$. If any eigenmode has a positive growth rate, the flow is unstable. Furthermore, taking the complex conjugate of Eq. (3) shows that for each unstable solution there exists a stable solution with equal and opposite growth rate.

_{j}^{6}Previous work

^{13}demonstrated that nonlinear interactions with these stable modes play an important role in saturating the growth of unstable modes. In the present work, we perform nonlinear simulations of an unstable shear flow and examine the role played by stable modes beyond the onset of saturation.

### B. Kolmogorov flow

One unstable flow profile of relevance in fusion and astrophysical systems is a sinusoidal equilibrium flow with periodic boundary conditions.^{2,4,21,25–31} For a sinusoidal flow profile $V(x)=V0\u2009cos\u2009(kxeqx)$ in a periodic domain, Eq. (3) lends itself well to being solved using spectral methods. Defining $\varphi \u0303j(kx,ky)$ as the Fourier series expansion of $\varphi \u0302j(x,ky)$, the Fourier representation of Eq. (3) is

where $\varphi \u0303j\xb1\u2261\varphi \u0303j(kx\xb1kxeq,ky)$. Equation (4) immediately demonstrates that each eigenmode exhibits a discrete, comb-like structure when viewed through a Fourier transform: for a given eigenmode $\varphi \u0302j(x,ky)$, if its Fourier transform $\varphi \u0303j(kx,ky)$ is nonzero at some *k _{x}*, then it is also nonzero at $kx+nkxeq$ for every integer

*n*(though $\varphi \u0303j$ is still expected to fall off at large $|kx|$, so that calculations with a finite number of

*k*can be expected to capture the structure well). This property of the system will have important consequences in Secs. III A and IV C when we compare simulations with different box sizes, and in Sec. V B when we explore the possibility of approximating the turbulent state by truncating the summation over

_{x}*j*in Eq. (2) to a reduced number of modes.

### C. Numerical implementation and benchmarking

We perform simulations of a KH-unstable sinusoidal *E *×* B* flow using the gyrokinetic framework^{36} as implemented in the Gene code.^{17,18} The gyrokinetic framework applies to systems with a strong guide field, where the parallel length scale of fluctuations is much larger than the perpendicular length scale, and the relevant frequencies are much smaller than the ion cyclotron frequency. The use of gyrokinetics for this work is motivated by Gene's unique tools for performing eigenmode decompositions.^{19,20} We simulate a system with two spatial dimensions, with a *y*-directed flow that varies sinusoidally in *x*, a strong guide field in the *z* direction, and no variations in *z*. Our simulation domain is a periodic box of dimensions $Lx\xd7Ly$ with no curvature or magnetic shear. The flow arises from the *E *×* B* drift of the particles, allowing the electrostatic potential *ϕ* to serve as the streamfunction for the flow. We model the plasma with gyrokinetic ions and electrons with hydrogen mass ratio, ion and electron background temperatures $Ti=Te$, no collisions, and no electromagnetic fluctuations (plasma *β* = 0).

We drive instability with a potential and corresponding distribution function that vary sinusoidally in *x*. Gene uses a *δf* formalism, where the full distribution function is separated into equilibrium *F*_{0} and fluctuation *f*, with the code solving for the evolution of the fluctuation. We let $f(x,y,v\u2225,\mu ,s,t)$ and $f\u0303(kx,ky,v\u2225,\mu ,s,t)$ denote the (guiding-center) distribution function for species *s* in real and Fourier space. For the remainder of this paper, we will generally use notation that suppresses the species and velocity dependence of *f*, and instead focus on its dependence on the spatial coordinates and time.

For benchmarking against previous work,^{21} the instability is first examined by implementing the sinusoidal flow with low-amplitude white noise as an initial condition in the fluctuation, formally evolving the system nonlinearly, with a homogeneous Maxwellian equilibrium *F*_{0}. This corresponds to solving the equation

with a sinusoidal initial condition $f(t=0),\varphi (t=0)\u223c\u2009sin(kxeqx)$ and low-amplitude noise to seed instability. The only term on the right-hand side of Eq. (5)

is the *E *×* B* nonlinearity, whose Fourier transform becomes $\u2211kx\u2032,ky\u2032(kx\u2032ky\u2212kxky\u2032)\varphi \xaf\u0303(kx\u2032,ky\u2032)f\u0303(kx\u2212kx\u2032,ky\u2212ky\u2032)$. Here, $\varphi \xaf$ is the gyro-averaged *ϕ*, whose Fourier transform is given by $\varphi \xaf\u0303(kx,ky,\mu ,s)=J0(kx2+ky2\rho )\varphi \u0303(kx,ky)$, where *J*_{0} is a Bessel function, and *ρ* is the gyroradius of species *s* with magnetic moment *μ*. The code evolves *f* according to Eq. (5) and calculates *ϕ* using Gauss's law as described in Refs. 37 and 38. The normalizations used by Gene are described in Ref. 37; however, in this paper, we will follow the standard convention used in the fluids community and normalize quantities with respect to the equilibrium flow velocity *V*_{0} and its wavelength $kxeq$, which are normalized in the code by $Vphys=V\rho scs/Lref$ and $kxphys=kx/\rho s$, where *c _{s}* is the ion sound speed and

*ρ*is the ion sound Larmor radius.

_{s}Consistent with fluid theory, our system is unstable to perturbations of the same form as Eq. (2) for a range of perturbation wavenumbers *k _{y}*, with the growth rate scaling with the base flow amplitude

*V*

_{0}. Growth rates from this formally nonlinear setup are indicated by crosses in Fig. 1. For direct comparison with previous work,

^{21}the wavenumber of the equilibrium $kxeq$ was varied at fixed

*k*, where perturbations are unstable for $0<ky/kxeq<1$. For this reason, in the remainder of this paper we focus our discussion on modes that lie in this range.

_{y}As demonstrated in Fig. 1, nonlinear simulations with appropriate initial conditions can be used to investigate some of the linear dynamics of this system, such as the growth rate and mode structure of the most unstable mode at each *k _{y}*. However, to solve for other linear modes, which are known to exist in fluid models,

^{6}terms corresponding to interaction with the driving flow need to be implemented as a linear operator, so that Eq. (5) can be linearized similarly to Eq. (3), in the form

for a linear differential operator $LKH$. To that end, we have implemented the linear operator $LKH$ in the Gene code. This allows computations to be performed with $LKH[f]$ on the right-hand side of the equation for $\u2202tf$ for

where $\varphi 0$ is the electrostatic potential (streamfunction) for the sinusoidal base flow, and *f*_{0} is the self-consistent distribution function corresponding to $\varphi 0$. Specifically, we use

where $\delta kx,kx\u2032$ is the Kronecker delta, and $F0(s)$, *q _{s}*, and

*T*are the equilibrium Maxwellian, charge, and temperature of species

_{s}*s*. Here, $J0(k\u22a5\rho )$, a Bessel function, and $\Gamma 0(b)=e\u2212bI\u03020(b)$ (where $I\u03020$ is a modified Bessel function, $b(s)=vTs2k\u22a52\Omega s\u22122/2$, and

*v*and Ω

_{Ts}_{s}are the thermal velocity and Larmor frequency of species

*s*) relate to finite Larmor radius (FLR) effects as detailed in Refs. 37 and 38. This form of

*f*

_{0}is used for secondary instability tests in tokamak-relevant systems

^{39}and yields a sinusoidal $\varphi 0(x)$ corresponding to a sinusoidal equilibrium flow in the

*y*direction with amplitude

*V*

_{0}and wavenumber $kxeq$.

Note that $LKH$ has *x* dependence but no *y* dependence, so its eigenmodes will have Fourier dependence in *y* and more complex structure in *x*, similar to the hydrodynamic case discussed in Sec. II A. The dots in Fig. 1 are obtained by solving Eq. (7), and their agreement with the corresponding crosses demonstrates successful implementation of the linear drive. For both setups, convergence checks were performed in spatial and velocity coordinates. Well-converged growth rates generally require 33 grid points in *x*, although far fewer points were required for $ky/kxeq\u22720.5$. For the remainder of this paper, results are presented with $V0=10$ and $kxeq=0.1$ using the linearized $LKH$. A convenient consequence of these parameters is that times and frequencies have the same value when expressed in standard Gene normalizations as they do in typical normalizations used in calculations of unstable shear flow in the fluids community, where *t* is often measured in units of $(kxeqV0)\u22121$.

Cyan triangles in Fig. 1 are obtained from solving the Orr-Sommerfeld equation with the Dedalus code^{40,41} (where $kxeq$ is the only length scale in the system) with a Reynolds number *Re* = 400. Their agreement with the other curves supports the notion that kinetic effects do not play a significant role in determining the growth rate of this mode. Crosses and dots in Fig. 1 corresponding to lower values of *k _{y}* show especially good agreement with the fluid results. As each curve represents a fixed

*k*with varying $kxeq$, FLR effects become more important as $kxeq$ increases (i.e., as $ky/kxeq$ decreases), suggesting that the reduced growth rates in the $ky=0.2$ simulations relative to the fluid results are due to FLR effects. In non-periodic shear layers, such as $V\u223ctanh(x)$, it is observed that FLR effects can be stabilizing or destabilizing depending on the alignment of the equilibrium vorticity and magnetic field.

_{y}^{3,12}Due to the sinusoidal nature of the flow studied here, the simulation domain includes regions where vorticity and magnetic field are aligned and where they are anti-aligned, suggesting that the FLR stabilization observed in our system is qualitatively different from what is found in shear layers. We speculate that the FLR stabilization is due to a reduction in the gyro-averaged potential $\varphi \xaf$, as $\varphi \xaf/\varphi $ generally decreases with increasing

*k*.

### D. Forcing and dissipation terms

In this paper, nonlinear calculations often include additional terms corresponding to forcing and dissipation, which we introduce here. Hyperdissipation $\u2212D\u22a5(kx4+ky4)f\u0303$ (Ref. 42) is employed to provide small-scale dissipation in place of collisions, which are not expected to sufficiently dissipate small-scale fluctuations at achievable resolutions within valid limits of collision models. We note that our hyperdissipation term differs from what is more standard in the fluids community, where $kx4+ky4$ is replaced by $(kx2+ky2)2$. A second dissipative term $\u2212Dradf\u0303$ is spatially uniform and sometimes referred to as radiative damping or friction.^{43} It absorbs energy transferred to large scales,^{43,44} while also serving as a “symmetry-breaking” parameter that adjusts the relative growth rates of linear modes without modifying their structure.

Finally, we introduce a Krook operator $\u2212DKrook\delta kx,\xb1kxeq\delta ky,0f\u0303$, where $\delta i,j$ is the Kronecker delta, to represent forcing of the unstable equilibrium and prevent it from decaying due to turbulent fluctuations.^{33} Aside from being linear in *f* and therefore not constant in time, this is identical to the inhomogeneous body forcing used in studies of Kolmogorov flow.^{22–24} While the sign of the Krook operator seems to suggest that it removes energy from the system, that is merely a consequence of our separation between equilibrium and perturbation. As explained in Ref. 45, the kinetic energy of the full flow is $E=\u222b|V+v|2dxdy$, so that if the $(kx,ky)=(kxeq,0)$ component of **v** opposes that of **V** and is not larger in magnitude, as we will see to be the case in Sec. IV, terms that appear to dissipate the “perturbation energy” $\u222b|v|2dxdy$ at that wavenumber will actually increase the true energy *E*.

Having constructed a linear operator that yields consistent results for the most unstable eigenmode's growth rate, we now address the rest of the spectrum of eigenvalues.

## III. EIGENSPECTRUM

### A. Subdominant modes

At each *k _{y}* there exists a spectrum of eigenmodes $f\u0302j$ and eigenvalues

*ω*, with corresponding potential structures $\varphi \u0302j$. For $0<ky/kxeq<1$, we let

_{j}*j*=

*1 denote the fastest-growing mode. The KH instability has been investigated in gyrokinetics before, but previous calculations did not address linear modes other than $f\u03021$ or their role in saturation. With the linear operator $LKH$ now implemented in Gene, its full spectrum of eigenmodes and eigenvalues can be obtained.*

^{19,20,46}Like the inviscid fluid analog, for each

*k*in the unstable range there exist one unstable mode, one stable (damped) mode with equal and opposite growth rate,

_{y}^{6}and a continuum of marginally stable modes,

^{47}shown in Fig. 2(a) for $ky/kxeq=2/3$. The additional degrees of freedom gained in gyrokinetics relative to a fluid calculation, by taking into account the velocity-space structure of multiple species, increases the rank of the discretized linear operator considerably and leads to many more marginally stable modes. A single point on the marginally stable continuum in Fig. 2 corresponds to hundreds of eigenmodes (depending on velocity-space resolution), each with similar electrostatic potentials but significantly different velocity-space structure.

Despite the added degrees of freedom in gyrokinetics, there are still only one stable and one unstable eigenmode per *k _{y}* for $0<ky/kxeq<1$ when the box size

*L*equals the wavelength of the equilibrium flow, denoted by $\lambda eq\u22612\pi /kxeq$. Consistent with magnetohydrodynamic (MHD) studies of a similar system,

_{x}^{25,27}we find that flows where multiple wavelengths of the equilibrium are present (i.e., setting $Lx=n\lambda eq$ where $n\u22652$ is an integer) exhibit pairs of subdominant unstable ($0<\gamma j<\gamma 1$) and stable ($\gamma 2<\gamma j<0$) modes, shown in Fig. 2(b). This means that simulations with larger boxes but with an equilibrium flow of the same wavelength are expected to have different dynamics than simulations with $Lx=\lambda eq$, as they include additional modes through which $LKH$ can inject or remove energy. In Sec. IV, we will demonstrate that including

*D*

_{Krook}and

*D*

_{rad}admits a system where, for sufficiently large

*L*, observables are converged with respect to a further increase in

_{x}*L*.

_{x}As stated above, for each *k _{y}* in $0<ky/kxeq<1$, we let

*j*=

*1 refer to the dominant unstable mode. We will further let*

*j*=

*2 refer to the corresponding stable mode, and*

*j*>

*2 to all other modes. Figure 3 shows the*

*x*-dependence of $\varphi \u03021$ at $ky/kxeq=1/2$ alongside the streamfunction for the equilibrium flow. Consistent with the fluid case,

^{6,13}we find $\gamma 2=\u2212\gamma 1$ (such that both $|\gamma 1,2|$ are reduced by FLR effects), and $\varphi \u03022(x,ky)=\varphi \u03021*(x,ky)$. Accordingly, we refer to

*f*

_{2}as a conjugate stable mode.

Consistent with the fluid case discussed in Sec. II B, the sinusoidal nature of the equilibrium gives eigenmodes a discrete, comb-like structure in *k _{x}*, where $f\u0303j$ is zero at every

*k*except for a countably infinite number of

_{x}*k*that are each separated by $kxeq$. All of the modes whose eigenvalues are plotted in Fig. 2(a), including $f\u03031$ and $f\u03032$, have nonzero amplitudes at integer multiples of $kxeq$. Many of the additional modes gained in Fig. 2(b) by extending

_{x}*L*to $2\lambda eq$, such as the modes with finite growth rate and real frequency, are nonzero at half-integer multiples of $kxeq$. This implies that arbitrary linear combinations of the modes in Fig. 2(a) are only nonzero at integer multiples of $kxeq$.

_{x}### B. Forcing and dissipation effects

The additional physics effects introduced in Sec. II D each modify the eigenmodes to varying degrees. The Krook operator enters the Vlasov equation only at *k _{y}* = 0, so it has no impact on the $ky>0$ eigenmode spectra. The radiative damping term reduces the growth rate of every eigenmode by $Drad$ without changing the mode structure. The hyperdissipation term has a more significant impact on the spectrum. It reduces the growth rate of the unstable mode with minor modifications to its structure and replaces both the stable mode and marginal continuum with a set of damped modes that does not include any mode resembling the conjugate stable mode.

We now turn our attention to the nonlinear saturation of this system.

## IV. INSTABILITY SATURATION

### A. Saturation and decaying turbulence

To investigate the saturation of this instability, we include in Eq. (7) the full *E *×* B* nonlinearity, yielding

Owing to the way in which the linear drive terms were derived and implemented, the evolution of Eq. (10) with some initial condition $finit$ is identical to the evolution of Eq. (5) with the initial condition $f0+finit$, presuming no dissipation or drive is included. When dissipation terms are added to Eq. (10), they do not act on *f*_{0}, unlike those in Eq. (5).

As the system evolves according to Eq. (10), the nonlinearity transfers energy across a range of scales, but with zero energy injection and nonzero dissipation, the initial energy eventually decays away. In terms of saturation of a linear instability, this can be understood as quasilinear flattening, where the fluctuations reduce mean gradients until the system is linearly stable. This is observed in simulations of Eq. (10) with added hyperdissipation, as shown in Figs. 4 and 5. Once unstable wavenumbers reach a sufficient amplitude, fluctuations at the wavenumbers of the equilibrium flow, i.e., $(kx,ky)=(\xb1kxeq,0)$, quickly grow to offset the unstable profile of the mean flow. From this point, the system exhibits features of decaying turbulence: the dynamics are highly intermittent, with long periods of coherent behavior punctuated by the merging of vortices. This is consistent with previous 2D KH simulations^{3} and can be expected given the lack of external forcing; the linear drive in Eq. (10) appears similar to an external forcing term, but as argued in the preceding paragraph, that is merely a consequence of the formal separation between the equilibrium and fluctuations.

### B. Driven turbulence

In many physical systems where shear-flow instability saturation and turbulence are of interest, the unstable shear flow is not some ideal initial condition but is brought about by a separate process. Examples include shear flows driven by boundary conditions,^{48} drift-wave instabilities^{4,21} in laboratory experiments, and jets, gravity, or another instability^{25} in astrophysical systems. We include a Krook operator, introduced in Sec. II, with the intent of capturing some of the effects of such continual forcing but without modeling the subtleties of any particular system where forcing produces a shear flow.

The result of including the Krook operator is readily seen in Figs. 4 and 6. When the Krook operator is added to Eq. (10), it suppresses the tendency for the $(kx,ky)=(\xb1kxeq,0)$ components of the fluctuation to cancel out *f*_{0}, thereby injecting energy into the system by reinforcing the unstable equilibrium. This in turn drives other Fourier modes via the KH instability, as is seen in the time trace of $\varphi \u0303(kxeq,0)$, shown in Fig. 4: the $(kx,ky)=(kxeq,0)$ component no longer reaches the amplitude necessary to cancel out the driving shear flow, and other Fourier modes no longer decay over time, leading to a quasi-stationary state of driven turbulence where the energy injected by the Krook drive is balanced by energy dissipation. As *D*_{Krook} increases, the saturated amplitude of $\varphi \u0303(\xb1kxeq,0)$ decreases, corresponding to an overall increase in $\varphi \u03030+\varphi \u0303(\xb1kxeq,0)$. The dominant balance that determines the amplitude of $\varphi \u0303(\xb1kxeq,0)$ in saturation is between the Krook drive and the Reynolds stress, which we explore further in Sec. IV C.

Also observed in Fig. 5 is the tendency for the system to form coherent vortices that gradually merge to the largest scale allowed by the simulation domain. This behavior is also observed in 2D shear layer simulations^{3} and is consistent with the inverse energy cascade to large scales in 2D hydrodynamics. The inverse cascade leads to a system with saturation properties that change as the box size is increased. The radiative damping term *D*_{rad} introduced in Sec. II damps low-*k* fluctuations, preventing energy from continuously building up at the largest scales, and thereby allowing fluctuation spectra to reach a stationary condition at low *k*. For this reason, and for the sake of presenting simulations where observables are converged with respect to the box size, the majority of our simulations were run with $Lx=12\lambda eq$ and *D*_{rad} = 0.05, a rate that is roughly 20% of the maximum linear growth rate in the dissipationless case. Figure 6 shows the results of a simulation with these parameters and *D*_{Krook} = 1. In contrast with Fig. 5, the system exhibits multiple excited scales in a quasi-stationary saturated state, providing the type of turbulence desired for studying momentum transport and eigenmode excitation.

### C. Momentum transport

We investigate the momentum transport driven by turbulent fluctuations in this system, examining the *xy* component of the Reynolds stress tensor, denoted as *τ*. We calculate τ from the average of the product of the *x* and *y* components of the fluctuating *E *×* B* flow in the homogeneous *y* direction

where $\u27e8A\u27e9q$ denotes an average of some quantity *A* over a domain in the variable *q*. Due to the sinusoidal variation in *x* of the equilibrium, *τ* changes sign along the *x* axis as the sign of the equilibrium flow changes, an expected feature of Kolmogorov flow.^{23}

In nonlinear gyrokinetic simulations, numerical convergence is typically tested by measuring changes of some scalar, time-averaged transport quantity with numerical parameters such as resolution and domain size. Due to the changes in sign of *τ*, the average of *τ* in the *x* direction and time, $\u27e8\tau \u27e9x,t$, is not appropriate for testing numerical convergence because it is typically 0. Instead, we calculate the root-mean-square (RMS) of *τ*, i.e., $\u27e8\tau 2\u27e9x$, and compare the time-average in the quasi-stationary state as resolution changes. For the simulation shown in Fig. 6, the time-averaged $\tau RMS$ in saturation changes by at most 2% when any spatial or velocity coordinate's domain size or resolution is doubled except *L _{x}* (expected due to the subdominant unstable modes and inverse cascade), where it changes by 9%. Therefore, despite creating additional unstable and stable eigenmodes as box size is increased, this simulation is numerically converged in

*L*with regard to $\tau RMS$.

_{x}Consistent with studies of Kolmogorov flow [where a constant force is typically used, while our forcing is proportional to $f\u0303(kxeq,0)$],^{23} we find that as the forcing increases, both the mean flow velocity and the Reynolds stress increase, such that at saturation the two are in balance. This is shown in Fig. 7, where the force on the mean flow applied by *D*_{Krook} is seen to balance the force due to *τ*. This can also be seen by considering the effect of a similar Krook operator on Eq. (1). When Eq. (1) is Fourier transformed in both *x* and *y*, our forcing term appears as $DKrook\delta kx,\xb1kxeq\delta ky,0(kx2+ky2)\varphi \u0303$. The $(kx,ky)=(kxeq,0)$ component of the equation then becomes

where the nonlinear term is the $kx=kxeq$ component of the Fourier-transformed Reynolds stress $\tau \u0303$. For a quasi-stationary, saturated state, the time-average of Eq. (12) yields a balance between the Reynolds stress and the *D*_{Krook} term. Figure 7 compares these terms for a range of simulations with different values of *D*_{Krook}, demonstrating good agreement with expectations. A small mismatch occurs because the effect of dissipation on the flow makes a small contribution to the force balance, but the other forces are clearly dominant. Because we only include dissipation on the fluctuation, not the equilibrium $\varphi 0$ which is independent of *D*_{Krook}, this contribution decreases as *D*_{Krook} increases.

## V. EIGENMODE ANALYSIS

### A. Eigenmode expansion

We investigate the role of stable modes in this system by expanding the turbulent state in a basis of the eigenmodes of the dissipationless operator $LKH$. We expand in eigenmodes of the dissipationless operator to focus on the role played by *f*_{2}, which, although it vanishes from the eigenmode spectrum in the dissipative system, may still characterize a component of the nonlinear state. This also allows for comparison with previous work,^{13} where the dissipationless modes were considered.

As discussed in Sec. III, the operator $LKH$ has a distinct set of $Nev$ eigenmodes ${f\u0302j}$ for each value of *k _{y}*. Therefore, an expansion of an arbitrary state $f(s,x,y,v\u2225,\mu )$ in a basis of eigenmodes $f\u0302j$ may be written as

As in Sec. III, the index *j* is a positive integer that enumerates the eigenmodes at a given *k _{y}*, and for $0<ky/kxeq<1,j=1$ and

*j*=

*2 label the most unstable mode and its stable conjugate, respectively. The number of eigenmodes*

*N*

_{ev}obtained by the eigenmode solver is equal to the number of degrees of freedom in the discrete numerical representation, i.e., the product of the number of grid points and the number of species, and the modes were verified to be linearly independent, so expansions of this form exist and are unique assuming the numerical resolutions of both sides of Eq. (13) are identical. Figure 8 shows time traces of $|\beta 1|$ in blue and $|\beta 2|$ in orange, as well as the time-averaged $|\beta j|$ in saturation for every

*j*at $ky/kxeq=0.25$ for the same simulation shown in Fig. 6.

The values *β _{j}*, which we refer to as the amplitudes of each eigenmode, can be understood as coordinates or components of

*f*in the basis ${fj}$. When such an expansion is performed at multiple time steps of a given simulation so that

*f*is a function of time, each

*β*becomes a function of time that indicates the relative contribution of eigenmode

_{j}*f*to the state of the system over time. In linear simulations, $\beta j(ky,t)=\beta j(ky,0)ei\omega jt$ for each

_{j}*j*and

*k*. Previous work showed how

_{y}*β*

_{1}and

*β*

_{2}interact nonlinearly in a fluid system, derived equations for $\u2202\beta j/\u2202t$ by inserting expansions of the form Eq. (13) into the governing equations of the system, and compared the relative sizes of different terms leading into instability saturation.

^{13,49,50}Here, we directly calculate the evolution of each

*β*over time in nonlinear simulations, extending analysis beyond the onset of saturation. Our procedure for calculating each

_{j}*β*relies on the left eigenmodes of $LKH$ and is described in Refs. 19 and 20.

_{j}Similar analyses have been performed for gyroradius-scale instabilities in reduced fluid models,^{50} and gyrokinetic models.^{19,51} These eigenmode expansions are related to the “projections” calculated in related work,^{16,51,52} defined as

(where $fNL$ is the nonlinearly evolved distribution function, and the summations are over each species), but the two are generally quite different. Identifying $\u27e8g,h\u27e9\u2261\u222bdxdvg*h$ as an inner product on the space of distribution functions *f*, projections *p _{j}* are inner products normalized by the lengths of

*f*and

_{j}*f*so that

*p*= 0 if they are orthogonal (under this inner product) and

_{j}*p*= 1 if they are parallel. The eigenvectors of an arbitrary linear operator are not guaranteed to be mutually orthogonal under a given inner product (we have verified that the eigenmodes of our system are not mutually orthogonal under the above inner product), which leads to the possibility that the projection onto one eigenvector depends on the amplitudes of every eigenvector. For example, one could find that the projection

_{j}*p*onto a stable mode counterintuitively grows over time in a linear simulation due to nonorthogonality, even though the amplitude

_{j}*β*of the stable mode decays. Likewise, if the projection onto a stable mode is large in the saturated state, it is not immediately clear whether this is due to a large stable mode amplitude, significant nonorthogonality with the dominant unstable mode, or even due to nonorthogonality with an entirely different mode that has a large amplitude. This situation is avoided if the linear operator has mutually orthogonal eigenvectors (e.g., if it is Hermitian), if the set of modes

_{j}*f*are replaced by an orthogonal set, such as from a proper orthogonal decomposition,

_{j}^{19}or by applying an orthogonalization procedure like Gram-Schmidt.

^{51}However, the relationship between the eigenmode amplitudes and the orthogonalized mode amplitudes is not immediately clear. We focus our attention on the eigenmode amplitudes

*β*rather than projections

_{j}*p*because linear energy transfer due to $LKH$ is directly related to

*β*, not

_{j}*p,*

^{49}and to facilitate comparison with Ref. 13.

For the simulation shown in Fig. 6, we use the parameters $Lx=12\lambda eq,DKrook=1,Drad=0.05$, and $D\u22a5=1.6$, with 512 grid points in the *x* direction. Calculating every eigenmode of the system at that resolution is prohibitively expensive. Instead, to generate Fig. 8 we perform eigenvalue computations with $Lx=\lambda eq$. Due to the discrete, comb-like eigenmode structure in *k _{x}* discussed in Secs. II B and III B, this reduced set of modes does not describe the full state of Eq. (13) because it lacks modes obtained when $Lx>\lambda eq$ [see Fig. 2]. But this does allow for a full expansion of the components of $f\u0303(kx,ky)$ given by $kx=nkxeq$ for integer

*n*, and this does not affect the obtained values of

*β*

_{1}and

*β*

_{2}.

Consistent with analytical calculations and reduced models,^{13,49,50} Fig. 8 shows that $|\beta 2|$ decays before being nonlinearly driven at a rate faster than the unstable mode's concurrent exponential growth. We stress that the evolution of $|\beta 2|$ is remarkably consistent with the inviscid fluid problem^{13} despite the influence of nonzero $D\u22a5$ in the nonlinear simulation, which modifies the structure of *f*_{1} and eliminates the conjugate stable mode *f*_{2} from the eigenmode spectrum of the dissipative operator. A similar observation was made in studies of ion temperature gradient (ITG) pseudospectra, where a similar conjugate stable mode vanished in the dissipative case, but was nonetheless a part of the pseudospectrum and was significantly excited in saturation when dissipation was included.^{53} Figure 8 only shows amplitudes for the $ky/kxeq=0.25$ eigenmodes, but every other *k _{y}* in $0<ky/kxeq<1$ exhibits similar results. The amplitude of

*f*

_{2}in saturation nearly matches that of

*f*

_{1}both at saturation onset and for the rest of the simulation. Since the two modes are nearly conjugate symmetric, this suggests that the linear energy dissipation due to

*f*

_{2}is a significant fraction of the linear energy injection due to

*f*

_{1}at the onset of saturation and throughout the quasi-stationary state. This suggests that the predictive capabilities of the threshold parameter

*P*analysis studied in Refs. 49 and 50 carry over to systems more general than plasma microturbulence, and that a significant amount of the energy transferred to

_{t}*k*> 0 fluctuations via $LKH$ makes its way back into the mean flow rather than smaller scales.

_{y}### B. Truncated eigenmode expansions

In turbulence models, it is common practice to separate the flow into mean and fluctuating parts. If there is further separation between large and small scale structures, the former are often approximated by the most unstable eigenmode.^{7,9,10} Here we demonstrate the potential for improving such models by including the stable mode in the approximation for the large scales.

Figure 9 compares part of the flow structure at $t\u2248501(V0kxeq)\u22121$ to three different expansions. The top-left contours show the electrostatic potential *ϕ* for the simulation described in Fig. 8. To focus on the components of *ϕ* where the eigenmodes discussed in Figs. 1 and 2(a) can be used to approximate the flow, a filtering procedure has been applied in Fig. 9 to remove all but a subset of Fourier components (*k _{x}*,

*k*). Only

_{y}*k*in $0<ky/kxeq<1$ are included, and only $kx=nkxeq$ for integer

_{y}*n*are included. This allows the eigenmodes in Fig. 8, and the equivalent eigenmodes at other unstable

*k*, to be used as a basis in the sense of Eq. (13). The top-right contours show the

_{y}*ϕ*structure obtained from summing over these eigenmodes at each unstable

*k*, verifying that they indeed serve as a basis. The excellent agreement helps demonstrate that the wavenumber filtering only affects the amplitudes

_{y}*β*of eigenmodes that arise from having $Lx>\lambda eq$, and fully captures the structure and amplitudes of the $Lx=\lambda eq$ eigenmodes. Extremely minor differences between the top-left and top-right contours arise due to the higher

_{j}*x*resolution in the nonlinear simulation than in the linear eigenmode calculations. To investigate the differences between these large-scale flows and the results of approximating them using just the unstable mode, the bottom-left contours show the result of excluding every eigenmode in Eq. (13) except the most unstable at each

*k*, as is often done in reduced models. The bottom-right contours are obtained similarly, but both the most unstable mode

_{y}*ϕ*

_{1}and the conjugate stable mode

*ϕ*

_{2}at each

*k*are included. Including

_{y}*ϕ*

_{2}produces a flow structure that is remarkably similar to the top-left and top-right flow structures, unlike what one obtains when only

*ϕ*

_{1}is included. Unsurprisingly, the more accurate flow structure leads to a more accurate Reynolds stress (not shown).

To compare the efficacy of these three eigenmode expansions over time, rather than the one timestep shown in Fig. 9, we calculate the error, given by $||\varphi \u2212\u2211j\beta j\varphi j||/||\varphi ||$, of each relative to the filtered nonlinear data (the top-left plot in Fig. 9), shown in Fig. 10. Here *ϕ* refers to the filtered nonlinear *ϕ*, and $||.||$ is the standard *L*_{2} norm. Due to differences in *x* resolution, the full expansion (in green) has minor errors that decay away as the simulation progresses. Errors in both the unstable-only expansion (blue) and the combined unstable-stable expansion (orange) start large due to choice of initial condition, gradually decay as the most unstable mode grows in the linear phase, and peak at the onset of saturation before fluctuating about an average value in the quasi-stationary state, with the inclusion of the stable mode reducing the average error in the saturated state by a factor of three.

### C. Influence of forcing and dissipation

Figure 8 shows significant excitation of the stable mode in the saturated state for a simulation with $DKrook=1,\u2009Drad=0.05$, and $D\u22a5=1.6$, with Figs. 9 and 10 demonstrating its importance in describing the large-scale fluctuations in *ϕ*. To investigate the role of these parameters in determining the influence of stable modes in saturation, we vary them between different simulations. In particular, we pay close attention to the relative amplitudes of *β*_{1} and *β*_{2} as *D*_{rad} and *D*_{Krook} are varied. Because *D*_{rad} is a symmetry-breaking term in the sense that it decreases the growth rate of *f*_{1} and increases the damping of *f*_{2} without changing their mode structures, it reduces the parametric driving of *f*_{2} by *f*_{1}. (The parametric driving of *f*_{2} by *f*_{1} depends on their mode structures, the form of the nonlinearity, and *γ*_{2}/*γ*_{1}.^{13,35,49,50} Of those, only *γ*_{2}/*γ*_{1} is affected by *D*_{rad}, making its influence on $|\beta 1/\beta 2|$ more transparent.) Figure 11 shows that this leads to significantly smaller $|\beta 2|$ relative to $|\beta 1|$ in the saturated state. This also suggests that for unstable shear flow in systems without radiative damping $|\beta 2|\u2248|\beta 1|$ is expected, consistent with Ref. 13. The *k _{y}* dependence of the ratio $|\beta 1/\beta 2|$ roughly follows that of

*γ*

_{1}, except that it approaches 1, rather than 0, at

*k*= 0 and $ky=kxeq$.

_{y}Figure 12 shows how the time-averaged values of $|\beta 1/\beta 2|$ in saturation vary with *D*_{Krook}. The shape of the curves remains fairly consistent as *D*_{Krook} changes. Two regimes are apparent: below *D*_{Krook} = 4, increasing *D*_{Krook} drives the ratio $|\beta 1/\beta 2|$ closer to unity, while above *D*_{Krook} = 4 the ratio is significantly less affected. This behavior is consistent with the notion that reinforcement of the unstable profile by larger *D*_{Krook} allows *β*_{2} to be nonlinearly pumped to its maximal level, whereas for smaller *D*_{Krook} the quasilinear depletion of the profile cuts off the pumping of *β*_{2} before it reaches its maximal level. Note that *f*_{2} tends to reduce the Reynolds stress *τ*, suggesting that the increase in *τ* with *D*_{Krook} must be due to an increase in overall fluctuation level, rather than a change in $|\beta 1/\beta 2|$.

### D. Influence of stable modes in analytical models

To better understand how the unstable and stable modes affect the mean flow in saturation, we develop a reduced model that expresses the mean flow amplitude in terms of *β*_{1} and *β*_{2}. The model considers a 2D inviscid, incompressible fluid, assumes $\varphi \u2248\beta 1\varphi 1+\beta 2\varphi 2$ for $0<ky/kxeq<1$, and assumes that the force applied by the Krook operator balances the turbulent Reynolds stress in saturation. These assumptions are consistent with the findings presented in Figs. 7, 9, and 10.

For perturbations about a sinusoidal equilibrium flow, the linearized system becomes Eq. (4), which was derived in Sec. II but is repeated here

Equation (4) can be expressed as a matrix equation $\omega j\varphi \u2192j+M\varphi \u2192j=0$ where the components of $\varphi \u2192j$ are $\varphi \u0303j$ at different *k _{x}* and the dimension of $\varphi \u2192j$ and

**M**is infinite. Reasonable approximations of the eigenmodes and eigenvalues can be obtained by solving the matrix equation with $\varphi \u0303j\u22600$ for some finite number of

*k*values, and $\varphi \u0303j=0$ for all other

_{x}*k*. This has previously been found useful in similar KH and tearing mode calculations.

_{x}^{28}For example, solving the system with $kx=0,\xb1kxeq$ yields

where the vectors are written in the form $(\varphi \u0303(\u2212kxeq),\varphi \u0303(0),\varphi \u0303(kxeq))T$ and $\kappa (ky)\u22612((kxeq)2+ky2)/((kxeq)2\u2212ky2)$.

To arrive at an expression of force balance between the Reynolds stress and Krook drive, we return to Eq. (12), which we repeat here

Considering a steady state where $\u2202\varphi \u0303(kxeq,0)/\u2202t=0$, and assuming $\varphi \u0303=\beta 1\varphi \u03031+\beta 2\varphi \u03032$ with just the $kx=0,\xb1kxeq$ Fourier modes considered, Eq. (12) can be manipulated to yield

From here, values for $|\beta 1|$ and $|\beta 2|$ can be inserted to arrive at values for the mean flow amplitude. In other systems, $|\beta j|$ have been calculated using statistical closures.^{15,54,55} Extending the above approach with such a calculation would yield a complete model, but is outside the scope of this paper. Instead, we insert values of $|\beta 1|$ and $|\beta 2|$ from nonlinear simulations into Eq. (15). Our interest is in the scaling of $\varphi \u0303(kxeq,0)$ with *D*_{Krook}, and what role *β*_{2} plays in that scaling. Thus, we perform simulations with a range of *D*_{Krook} and compare three quantities: the time-averaged value of $\varphi \u0303(kxeq,0)$ in saturation, the result of Eq. (15) using both *β*_{1} and *β*_{2}, and the result obtained when *β*_{2} = 0 is assumed. In Fig. 13, this comparison is made for three values of *D*_{rad}, corresponding to three systems with varying degrees of stable mode excitation (recall stable modes are more excited at lower values of *D*_{rad}, see Fig. 11). For *D*_{rad} = 0.05, where stable modes were shown in Fig. 8 to be significant, the *β*_{2} = 0 model captures general trends, but is significantly improved when *β*_{2} is properly included. For *D*_{rad} = 0.025, where stable modes are even more important, the *β*_{2} = 0 model fails to even capture the decrease in $|\varphi \u0303(kxeq,0)|$ with *D*_{Krook}, while the model that includes *β*_{2} does capture the correct qualitative and often even quantitative behavior. At *D*_{rad} = 0.1, where stable modes are significantly weakened, their inclusion does not have a significant impact on the model. For each value of *D*_{rad}, the two models are scaled by a constant so that they agree with the simulation results at *D*_{Krook} = 4, which is where the change in scaling with respect to *D*_{Krook} was noted in Fig. 12. (It is the scaling properties of the models that we are assessing, not the absolute values.) Note that this model neglects all eigenmodes except *ϕ*_{1} and *ϕ*_{2}, including the modes with nonzero $\varphi \u0303$ at noninteger multiples of $kxeq$.

Comparing the two models at different values of *D*_{rad} demonstrates that when stable modes are excited in this system as in Fig. 8, they not only modify the shape of the flow, as shown in Fig. 9, but have an important impact on how the system responds to forcing. By nonlinearly transferring energy into large-scale stable modes, the fluctuating flow adjusts in a way that changes the feedback onto the large-scale mean flow, thus affecting how the system is forced. We also note that Eq. (15) was derived assuming an inviscid fluid, while the inserted values for $|\beta 1|$ and $|\beta 2|$ were obtained from gyrokinetic simulations with finite $D\u22a5$, suggesting gyrokinetic effects may not play a significant role in the eigenmode excitations and force balance in this system.

## VI. CONCLUSIONS

We have studied an unstable gyrokinetic shear flow, finding that a conjugate stable eigenmode is nonlinearly driven to a large amplitude leading into saturation, and continues to make important contributions to the Reynolds stress in the quasi-stationary turbulent state, except at high values of radiative damping. This demonstrates that previous findings on the role of stable eigenmodes in shear-flow instability saturation of a fluid shear layer are consistent with the quasi-stationary turbulent state of a gyrokinetic periodic shear flow. Furthermore, our results point to the potential for significantly improving reduced models of shear-driven turbulence by including stable-mode physics.

We have investigated the saturation of a linearly unstable *E *×* B* shear flow in gyrokinetics as it relates to the full eigenmode spectrum. We find that the gyrokinetic system compares well with its hydrodynamic counterpart with regard to the unstable mode, as well as the rest of the spectrum. Specifically, the dissipationless linear operator includes a single conjugate stable eigenmode for every unstable eigenmode, along with a continuum of marginally stable modes. Nonlinear simulations characterize the behavior of the flow in saturation, and we examine cases both with and without an external driving term. The drive is implemented in the form of a Krook operator and reinforces the unstable mean flow in a manner similar to Kolmogorov flow.

In simulations without the drive term, the system lacks any energy injection to maintain the unstable equilibrium. This causes fluctuations to quickly relax the unstable flow shear once nonlinear interactions become significant, and the turbulence subsequently decays. In simulations with forcing, we include a scale-independent radiative damping term that prevents accumulation of energy at the largest scales, and allows a quasi-stationary state of driven turbulence. In driven simulations, a partial relaxation of the mean flow is still observed, with the final state mostly determined by a force balance between the Krook drive and the turbulent Reynolds stress.

With a well-resolved system of quasi-stationary, driven turbulence, we investigate the role of different linear eigenmodes by performing an eigenvalue decomposition, where the turbulent state is expressed as a linear combination of the eigenmodes of the dissipationless linear operator. The evolution of the dominant pair of stable and unstable modes leading into saturation compares well with previous analytic calculations of an inviscid fluid shear layer,^{13} and the ensuing excitation of the stable mode in the turbulent state is consistent with previous findings in plasma microturbulence.^{50} By demonstrating that the role of stable modes in shear-flow instability saturation is consistent with their role in the fully developed turbulent system, we have extended the set of systems in which instability saturation analyses have proven to be predictive of the turbulent state to include fully global fluid instabilities, further motivating these sorts of analyses in other global instabilities where stable modes exist, such as the magnetorotational instability.^{56}

The significant excitation of linearly stable modes in the saturated state indicates that an important aspect of shear-driven turbulence is this previously neglected tendency for large-scale fluctuations to lose energy back into the mean flow via the linear operator. This idea is in contrast with the standard picture of instability-driven turbulence, where it is assumed that the largest scales are dominated by a balance between linear energy injection and nonlinear energy transfer to smaller scales. While many other modes are also excited in the saturated state, we have shown that the stable/unstable pair of modes is sufficient to capture many aspects of the flow. This also presents a significant modification to the existing understanding of shear-driven turbulence, where reduced models generally assume that large-scale fluctuations are dominated by unstable modes alone.^{7,9–11}

Consistent with previous work where the conjugate symmetry between unstable/stable pairs of modes was broken with dissipative terms,^{35} we find that the added radiative damping term, which increases the damping rate of the stable mode and reduces the growth rate of the unstable mode, suppresses the importance of the stable mode relative to the unstable one. This is observed by comparing the amplitudes of the two modes for a range of radiative damping values. Making use of the observations that the gyrokinetic and fluid systems behave similarly, that the mean flow amplitude at saturation is determined by force balance between driving and Reynolds stress, and that the stable and unstable modes alone describe large-scale fluctuations well despite the conjugate stable mode vanishing from the eigenmode spectrum when dissipation is included, we construct a reduced model that allows us to examine the role of stable modes in determining the mean flow in saturation. The model results in an equation where the contributions from stable modes can be isolated from unstable modes. We find that lower values of radiative damping, where stable modes exhibit higher amplitudes, require the inclusion of stable modes in the model in order for it to be even qualitatively correct. At higher radiative damping, where stable modes are suppressed, their inclusion in the model has no significant impact on its performance. Thus, in shear-flow systems where stable modes play an important role in instability saturation, they may also be expected to play an important role in understanding how fluctuations affect the mean flow, and thus how the system responds to external forcing. We further conclude that when effects observed to change turbulence characteristics also break the conjugate symmetry of an unstable/stable eigenmode pair,^{8} the change in turbulence may be related to differences in stable mode excitation.

## ACKNOWLEDGMENTS

The authors would like to thank B. N. Rogers, W. Dorland, M. E. Pessah, and F. Waleffe for valuable discussions and insights. Partial support for this work was provided by the National Science Foundation under Award No. PHY-1707236, the Wisconsin Alumni Research Foundation, the Vilas Trust, and the U.S. Department of Energy, Office of Science, Fusion Energy Sciences, under Award Nos. DE-FG02-89ER53291 and DE-FG02-04ER-54742. Computing resources were provided by the National Science Foundation through XSEDE computing resources, Allocation No. TG-PHY130027.