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.

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 equations5,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 (kz = 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 instability25–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 layer13 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 Â(x,ky) denotes the Fourier transform in y of A(x, y), and Ã(kx,ky) denotes the Fourier transform in x and y.

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

t2ϕ+Vy2ϕd2Vdx2ϕy+ϕxy2ϕϕyx2ϕ=0,
(1)

where V(x) is the y-directed equilibrium shear flow, and ϕ(x,y,t) is the streamfunction of the perturbation v=ϕ×ẑ. The linear dynamics can then be explored by dropping the nonlinearities and using the normal mode ansatz

ϕ(x,y,t)=kyjϕ̂j(x,ky)ei(kyy+ωjt),
(2)

yielding

(ωj+kyV)(2x2ky2)ϕ̂jkyϕ̂jd2Vdx2=0.
(3)

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 ϕ̂j, with j enumerating the eigenmodes at a given ky. The eigenvalue ωj is complex, with real frequency Re(ωj) and growth rate γj=Im(ω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.6 Previous work13 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.

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)=V0cos(kxeqx) in a periodic domain, Eq. (3) lends itself well to being solved using spectral methods. Defining ϕ̃j(kx,ky) as the Fourier series expansion of ϕ̂j(x,ky), the Fourier representation of Eq. (3) is

ωj(kx2+ky2)ϕ̃j+kyV02[(kx22kxkxeq+ky2)ϕ̃j+(kx2+2kxkxeq+ky2)ϕ̃j+]=0,
(4)

where ϕ̃j±ϕ̃j(kx±kxeq,ky). Equation (4) immediately demonstrates that each eigenmode exhibits a discrete, comb-like structure when viewed through a Fourier transform: for a given eigenmode ϕ̂j(x,ky), if its Fourier transform ϕ̃j(kx,ky) is nonzero at some kx, then it is also nonzero at kx+nkxeq for every integer n (though ϕ̃j is still expected to fall off at large |kx|, so that calculations with a finite number of kx 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 j in Eq. (2) to a reduced number of modes.

We perform simulations of a KH-unstable sinusoidal E × B flow using the gyrokinetic framework36 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×Ly 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 F0 and fluctuation f, with the code solving for the evolution of the fluctuation. We let f(x,y,v,μ,s,t) and f̃(kx,ky,v,μ,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 F0. This corresponds to solving the equation

ft={f,ϕ¯}
(5)

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

{f,ϕ¯}fxϕ¯yϕ¯xfy,
(6)

is the E × B nonlinearity, whose Fourier transform becomes kx,ky(kxkykxky)ϕ¯̃(kx,ky)f̃(kxkx,kyky). Here, ϕ¯ is the gyro-averaged ϕ, whose Fourier transform is given by ϕ¯̃(kx,ky,μ,s)=J0(kx2+ky2ρ)ϕ̃(kx,ky), where J0 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 V0 and its wavelength kxeq, which are normalized in the code by Vphys=Vρscs/Lref and kxphys=kx/ρs, where cs is the ion sound speed and ρs is the ion sound Larmor radius.

Consistent with fluid theory, our system is unstable to perturbations of the same form as Eq. (2) for a range of perturbation wavenumbers ky, with the growth rate scaling with the base flow amplitude V0. 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 ky, 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.

FIG. 1.

Dispersion relation for the KH instability of a sinusoidal flow V=V0cos(kxeqx)ŷ. Growth rate γ is plotted against the wavenumber ky of the perturbation, with γ normalized to the equilibrium shear kxeqV0 and ky to the equilibrium wavenumber kxeq. Crosses are obtained from nonlinearly evolving a perturbed sinusoidal flow in Gene according to Eq. (5), while dots are from solving the linear Eq. (7). Results compare well with both previous gyrokinetic simulations (red curve, see Ref. 21) and hydrodynamic simulations of an equivalent system (magenta triangles). The stabilization of the ky=0.2 points at low ky/kxeq (i.e., high kxeq) can be attributed to finite Larmor radius effects. All modes have zero real frequency.

FIG. 1.

Dispersion relation for the KH instability of a sinusoidal flow V=V0cos(kxeqx)ŷ. Growth rate γ is plotted against the wavenumber ky of the perturbation, with γ normalized to the equilibrium shear kxeqV0 and ky to the equilibrium wavenumber kxeq. Crosses are obtained from nonlinearly evolving a perturbed sinusoidal flow in Gene according to Eq. (5), while dots are from solving the linear Eq. (7). Results compare well with both previous gyrokinetic simulations (red curve, see Ref. 21) and hydrodynamic simulations of an equivalent system (magenta triangles). The stabilization of the ky=0.2 points at low ky/kxeq (i.e., high kxeq) can be attributed to finite Larmor radius effects. All modes have zero real frequency.

Close modal

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 ky. 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

ft=LKH[f]
(7)

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 tf for

LKH[f]{f0,ϕ¯}+{f,ϕ¯0},
(8)

where ϕ0 is the electrostatic potential (streamfunction) for the sinusoidal base flow, and f0 is the self-consistent distribution function corresponding to ϕ0. Specifically, we use

f0(s)=V0kxeqδkx,kxeqδkx,kxeq2iF0(s)qsTsJ01Γ0Γ0,
(9)

where δkx,kx is the Kronecker delta, and F0(s), qs, and Ts are the equilibrium Maxwellian, charge, and temperature of species s. Here, J0(kρ), a Bessel function, and Γ0(b)=ebÎ0(b) (where Î0 is a modified Bessel function, b(s)=vTs2k2Ωs2/2, and vTs and Ω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 f0 is used for secondary instability tests in tokamak-relevant systems39 and yields a sinusoidal ϕ0(x) corresponding to a sinusoidal equilibrium flow in the y direction with amplitude V0 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/kxeq0.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)1.

Cyan triangles in Fig. 1 are obtained from solving the Orr-Sommerfeld equation with the Dedalus code40,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 ky show especially good agreement with the fluid results. As each curve represents a fixed ky 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 Vtanh(x), it is observed that FLR effects can be stabilizing or destabilizing depending on the alignment of the equilibrium vorticity and magnetic field.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 ϕ¯, as ϕ¯/ϕ generally decreases with increasing k.

In this paper, nonlinear calculations often include additional terms corresponding to forcing and dissipation, which we introduce here. Hyperdissipation D(kx4+ky4)f̃ (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 Dradf̃ 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 DKrookδkx,±kxeqδky,0f̃, where δ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=|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” |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.

At each ky there exists a spectrum of eigenmodes f̂j and eigenvalues ωj, with corresponding potential structures ϕ̂j. For 0<ky/kxeq<1, we let 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̂1 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 ky in the unstable range there exist one unstable mode, one stable (damped) mode with equal and opposite growth rate,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.

FIG. 2.

Eigenvalue spectra for ky/kxeq=2/3 with Lx=λeq (a) and Lx=2λeq (b). With Lx=λeq, at each unstable ky, the spectrum includes one unstable and one stable mode with equal and opposite growth rate γ, as well as a continuous spectrum of marginal modes corresponding to resonances between the phase velocity ω/ky and equilibrium flow.47 As the box size is increased to fit multiple wavelengths of the equilibrium, additional stable and unstable modes are introduced,25,27 and additional marginal eigenvalues appear due to an increase in number of values of V0cos(kxeqx) sampled by the extended grid (thus additional resonances with ω/ky).

FIG. 2.

Eigenvalue spectra for ky/kxeq=2/3 with Lx=λeq (a) and Lx=2λeq (b). With Lx=λeq, at each unstable ky, the spectrum includes one unstable and one stable mode with equal and opposite growth rate γ, as well as a continuous spectrum of marginal modes corresponding to resonances between the phase velocity ω/ky and equilibrium flow.47 As the box size is increased to fit multiple wavelengths of the equilibrium, additional stable and unstable modes are introduced,25,27 and additional marginal eigenvalues appear due to an increase in number of values of V0cos(kxeqx) sampled by the extended grid (thus additional resonances with ω/ky).

Close modal

Despite the added degrees of freedom in gyrokinetics, there are still only one stable and one unstable eigenmode per ky for 0<ky/kxeq<1 when the box size Lx equals the wavelength of the equilibrium flow, denoted by λeq2π/kxeq. Consistent with magnetohydrodynamic (MHD) studies of a similar system,25,27 we find that flows where multiple wavelengths of the equilibrium are present (i.e., setting Lx=nλeq where n2 is an integer) exhibit pairs of subdominant unstable (0<γj<γ1) and stable (γ2<γ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=λeq, as they include additional modes through which LKH can inject or remove energy. In Sec. IV, we will demonstrate that including DKrook and Drad admits a system where, for sufficiently large Lx, observables are converged with respect to a further increase in Lx.

As stated above, for each ky 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 ϕ̂1 at ky/kxeq=1/2 alongside the streamfunction for the equilibrium flow. Consistent with the fluid case,6,13 we find γ2=γ1 (such that both |γ1,2| are reduced by FLR effects), and ϕ̂2(x,ky)=ϕ̂1*(x,ky). Accordingly, we refer to f2 as a conjugate stable mode.

FIG. 3.

The equilibrium potential ϕ0 considered throughout this paper, which generates an E × B flow along the y-axis that varies sinusoidally in the x direction with wavenumber kxeq, alongside the real and imaginary parts of the potential corresponding to the unstable eigenmode ϕ̂1(x,ky) plotted with respect to x at ky/kxeq=1/2. The stable eigenmode's potential is the complex conjugate of the unstable eigenmode's potential, ϕ̂2=ϕ̂1*.

FIG. 3.

The equilibrium potential ϕ0 considered throughout this paper, which generates an E × B flow along the y-axis that varies sinusoidally in the x direction with wavenumber kxeq, alongside the real and imaginary parts of the potential corresponding to the unstable eigenmode ϕ̂1(x,ky) plotted with respect to x at ky/kxeq=1/2. The stable eigenmode's potential is the complex conjugate of the unstable eigenmode's potential, ϕ̂2=ϕ̂1*.

Close modal

Consistent with the fluid case discussed in Sec. II B, the sinusoidal nature of the equilibrium gives eigenmodes a discrete, comb-like structure in kx, where f̃j is zero at every kx except for a countably infinite number of kx that are each separated by kxeq. All of the modes whose eigenvalues are plotted in Fig. 2(a), including f̃1 and f̃2, have nonzero amplitudes at integer multiples of kxeq. Many of the additional modes gained in Fig. 2(b) by extending Lx to 2λ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.

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 ky = 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.

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

ft=LKH[f]+{f,ϕ¯}.
(10)

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 f0, 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)=(±kxeq,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 simulations3 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.

FIG. 4.

Left: Nonlinear simulations with dissipation quasilinearly flatten, then decay. Quasilinear flattening is measured by investigating ϕ̃ at (kx,ky)=(±kxeq,0). The perturbation cancels the drive once ϕ̃(kxeq,0) (blue) reaches a magnitude of 0.5V0/kxeq (black dashed line). Linearly unstable Fourier modes then turbulently decay over time. Right: Introducing a Krook operator (DKrook/(kxeqV0)=1 here) partially suppresses the Fourier mode responsible for quasilinear flattening, driving the system and leading to a quasi-stationary state of driven turbulence.

FIG. 4.

Left: Nonlinear simulations with dissipation quasilinearly flatten, then decay. Quasilinear flattening is measured by investigating ϕ̃ at (kx,ky)=(±kxeq,0). The perturbation cancels the drive once ϕ̃(kxeq,0) (blue) reaches a magnitude of 0.5V0/kxeq (black dashed line). Linearly unstable Fourier modes then turbulently decay over time. Right: Introducing a Krook operator (DKrook/(kxeqV0)=1 here) partially suppresses the Fourier mode responsible for quasilinear flattening, driving the system and leading to a quasi-stationary state of driven turbulence.

Close modal
FIG. 5.

Contours of the full (equilibrium and fluctuation) electrostatic potential for a nonlinear simulation with DKrook=0. From left to right, plots correspond to t(kxeqV0)46,101, and 502. Center and right plots show the tendency for small-scale fluctuations to dissipate, leaving coherent vortices that merge to progressively larger scales.

FIG. 5.

Contours of the full (equilibrium and fluctuation) electrostatic potential for a nonlinear simulation with DKrook=0. From left to right, plots correspond to t(kxeqV0)46,101, and 502. Center and right plots show the tendency for small-scale fluctuations to dissipate, leaving coherent vortices that merge to progressively larger scales.

Close modal

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 instabilities4,21 in laboratory experiments, and jets, gravity, or another instability25 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)=(±kxeq,0) components of the fluctuation to cancel out f0, 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 ϕ̃(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 DKrook increases, the saturated amplitude of ϕ̃(±kxeq,0) decreases, corresponding to an overall increase in ϕ̃0+ϕ̃(±kxeq,0). The dominant balance that determines the amplitude of ϕ̃(±kxeq,0) in saturation is between the Krook drive and the Reynolds stress, which we explore further in Sec. IV C.

FIG. 6.

Contours of the full (equilibrium and fluctuation) electrostatic potential for a nonlinear simulation with DKrook/(kxeqV0)=1 and Drad/(kxeqV0)=0.05. From left to right, plots correspond to t(kxeqV0)46,102, and 501. Comparing with Fig. 5 shows that the system no longer tends towards large-scale coherent vortices with gradual decay of energy. Instead, multiple scales are excited and form a quasi-stationary state.

FIG. 6.

Contours of the full (equilibrium and fluctuation) electrostatic potential for a nonlinear simulation with DKrook/(kxeqV0)=1 and Drad/(kxeqV0)=0.05. From left to right, plots correspond to t(kxeqV0)46,102, and 501. Comparing with Fig. 5 shows that the system no longer tends towards large-scale coherent vortices with gradual decay of energy. Instead, multiple scales are excited and form a quasi-stationary state.

Close modal

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 simulations3 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 Drad 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λeq and Drad = 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 DKrook = 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.

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

τϕxϕyy,
(11)

where Aq 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, τx,t, is not appropriate for testing numerical convergence because it is typically 0. Instead, we calculate the root-mean-square (RMS) of τ, i.e., τ2x, and compare the time-average in the quasi-stationary state as resolution changes. For the simulation shown in Fig. 6, the time-averaged τRMS in saturation changes by at most 2% when any spatial or velocity coordinate's domain size or resolution is doubled except Lx (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 Lx with regard to τRMS.

Consistent with studies of Kolmogorov flow [where a constant force is typically used, while our forcing is proportional to f̃(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 DKrook 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δkx,±kxeqδky,0(kx2+ky2)ϕ̃. The (kx,ky)=(kxeq,0) component of the equation then becomes

tϕ̃(kxeq,0)+kkykxeq[(kxeqkx)2+ky2]ϕ̃(kx,ky)ϕ̃(kxeqkx,ky)=DKrookϕ̃(kxeq,0),
(12)

where the nonlinear term is the kx=kxeq component of the Fourier-transformed Reynolds stress τ̃. For a quasi-stationary, saturated state, the time-average of Eq. (12) yields a balance between the Reynolds stress and the DKrook term. Figure 7 compares these terms for a range of simulations with different values of DKrook, 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 ϕ0 which is independent of DKrook, this contribution decreases as DKrook increases.

FIG. 7.

Comparison between average Krook drive amplitude |ϕ̃(±kxeq,0)t|DKrook and the average amplitude of the corresponding Fourier component of the Reynolds stress τ in saturation across a range of driving frequencies DKrook. Other simulation parameters are the same as in Fig. 6. In the saturated state, the mean flow is governed by a competition between the external forcing and the turbulent Reynolds stress. A small contribution is made by the influence of dissipation on ϕ̃(±kxeq,0), evidenced by the minor mismatch between the two curves at the lowest values of DKrook.

FIG. 7.

Comparison between average Krook drive amplitude |ϕ̃(±kxeq,0)t|DKrook and the average amplitude of the corresponding Fourier component of the Reynolds stress τ in saturation across a range of driving frequencies DKrook. Other simulation parameters are the same as in Fig. 6. In the saturated state, the mean flow is governed by a competition between the external forcing and the turbulent Reynolds stress. A small contribution is made by the influence of dissipation on ϕ̃(±kxeq,0), evidenced by the minor mismatch between the two curves at the lowest values of DKrook.

Close modal

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 f2, 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̂j} for each value of ky. Therefore, an expansion of an arbitrary state f(s,x,y,v,μ) in a basis of eigenmodes f̂j may be written as

f(s,x,y,v,μ)=kyj=1Nevβj(ky)f̂j(s,x,ky,v,μ)eikyy.
(13)

As in Sec. III, the index j is a positive integer that enumerates the eigenmodes at a given ky, 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 Nev 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 |β1| in blue and |β2| in orange, as well as the time-averaged |βj| in saturation for every j at ky/kxeq=0.25 for the same simulation shown in Fig. 6.

FIG. 8.

Left: Amplitudes of the unstable and stable eigenmodes, |β1| (blue) and |β2| (orange), respectively, as functions of time on a logarithmic scale (top, horizontal axis reduced to highlight the parametric growth of |β2|) and a linear scale (bottom) at horizontal wavenumber ky/kxeq=0.25 for the simulation with DKrook=1, Drad = 0.05, and D=1.6. Right: Full spectrum of eigenvalues ωj (real part ωr on the y-axis and growth rate γ on the x-axis), with color indicating time-averaged (starting from t =300) amplitudes |βj|, and dot size scaled proportionally to allow multiple |βj| with the same ωj to be shown. The decay of β2 is followed by nonlinear growth much faster than β1, while β1 continues its linear growth, consistent with Ref. 13. The remarkable similarity of the values of |β1| and |β2| in the saturated state was predicted in Ref. 13. Results are qualitatively similar for other unstable ky.

FIG. 8.

Left: Amplitudes of the unstable and stable eigenmodes, |β1| (blue) and |β2| (orange), respectively, as functions of time on a logarithmic scale (top, horizontal axis reduced to highlight the parametric growth of |β2|) and a linear scale (bottom) at horizontal wavenumber ky/kxeq=0.25 for the simulation with DKrook=1, Drad = 0.05, and D=1.6. Right: Full spectrum of eigenvalues ωj (real part ωr on the y-axis and growth rate γ on the x-axis), with color indicating time-averaged (starting from t =300) amplitudes |βj|, and dot size scaled proportionally to allow multiple |βj| with the same ωj to be shown. The decay of β2 is followed by nonlinear growth much faster than β1, while β1 continues its linear growth, consistent with Ref. 13. The remarkable similarity of the values of |β1| and |β2| in the saturated state was predicted in Ref. 13. Results are qualitatively similar for other unstable ky.

Close modal

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 βj becomes a function of time that indicates the relative contribution of eigenmode fj to the state of the system over time. In linear simulations, βj(ky,t)=βj(ky,0)eiωjt for each j and ky. Previous work showed how β1 and β2 interact nonlinearly in a fluid system, derived equations for βj/t 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 βj 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.

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

pj=|dxdvsfj*fNL|(dxdvs|fj|2dxdvs|fNL|2)1/2
(14)

(where fNL is the nonlinearly evolved distribution function, and the summations are over each species), but the two are generally quite different. Identifying g,hdxdvg*h as an inner product on the space of distribution functions f, projections pj are inner products normalized by the lengths of fj and f so that pj = 0 if they are orthogonal (under this inner product) and pj = 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 pj 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 fj are replaced by an orthogonal set, such as from a proper orthogonal decomposition,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 βj rather than projections p because linear energy transfer due to LKH is directly related to βj, not p,49 and to facilitate comparison with Ref. 13.

For the simulation shown in Fig. 6, we use the parameters Lx=12λeq,DKrook=1,Drad=0.05, and D=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=λeq. Due to the discrete, comb-like eigenmode structure in kx 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>λeq [see Fig. 2]. But this does allow for a full expansion of the components of f̃(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 |β2| decays before being nonlinearly driven at a rate faster than the unstable mode's concurrent exponential growth. We stress that the evolution of |β2| is remarkably consistent with the inviscid fluid problem13 despite the influence of nonzero D in the nonlinear simulation, which modifies the structure of f1 and eliminates the conjugate stable mode f2 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 ky in 0<ky/kxeq<1 exhibits similar results. The amplitude of f2 in saturation nearly matches that of f1 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 f2 is a significant fraction of the linear energy injection due to f1 at the onset of saturation and throughout the quasi-stationary state. This suggests that the predictive capabilities of the threshold parameter Pt 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 ky > 0 fluctuations via LKH makes its way back into the mean flow rather than smaller scales.

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 t501(V0kxeq)1 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 (kx, ky). Only ky in 0<ky/kxeq<1 are included, and only kx=nkxeq for integer n are included. This allows the eigenmodes in Fig. 8, and the equivalent eigenmodes at other unstable ky, to be used as a basis in the sense of Eq. (13). The top-right contours show the ϕ structure obtained from summing over these eigenmodes at each unstable ky, verifying that they indeed serve as a basis. The excellent agreement helps demonstrate that the wavenumber filtering only affects the amplitudes βj of eigenmodes that arise from having Lx>λeq, and fully captures the structure and amplitudes of the Lx=λeq eigenmodes. Extremely minor differences between the top-left and top-right contours arise due to the higher 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 ky, as is often done in reduced models. The bottom-right contours are obtained similarly, but both the most unstable mode ϕ1 and the conjugate stable mode ϕ2 at each ky are included. Including ϕ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).

FIG. 9.

Comparison between the components of ϕ that are spanned by the eigenmodes in Fig. 8 (top left), a summation over all of the eigenmodes in Fig. 8 at every unstable ky (top right), summation over just the most unstable mode at every unstable ky (bottom left), and summation over the most unstable and conjugate stable mode at every unstable ky (bottom right) at t501(V0kxeq)1 for the same simulation as Fig. 6. Due to only integer multiples of kxeq contributing to these eigenmodes, they are unable to effectively reproduce the full flow profile, plotted in Fig. 6. However, those components of ϕ that can be expressed as a linear combination of the eigenmodes in Fig. 8 are very well-described even by just the unstable ϕ1 and stable ϕ2.

FIG. 9.

Comparison between the components of ϕ that are spanned by the eigenmodes in Fig. 8 (top left), a summation over all of the eigenmodes in Fig. 8 at every unstable ky (top right), summation over just the most unstable mode at every unstable ky (bottom left), and summation over the most unstable and conjugate stable mode at every unstable ky (bottom right) at t501(V0kxeq)1 for the same simulation as Fig. 6. Due to only integer multiples of kxeq contributing to these eigenmodes, they are unable to effectively reproduce the full flow profile, plotted in Fig. 6. However, those components of ϕ that can be expressed as a linear combination of the eigenmodes in Fig. 8 are very well-described even by just the unstable ϕ1 and stable ϕ2.

Close modal

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 ||ϕjβjϕj||/||ϕ||, 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 L2 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.

Figure 8 shows significant excitation of the stable mode in the saturated state for a simulation with DKrook=1,Drad=0.05, and D=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 Drad and DKrook are varied. Because Drad is a symmetry-breaking term in the sense that it decreases the growth rate of f1 and increases the damping of f2 without changing their mode structures, it reduces the parametric driving of f2 by f1. (The parametric driving of f2 by f1 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 Drad, making its influence on |β1/β2| more transparent.) Figure 11 shows that this leads to significantly smaller |β2| relative to |β1| in the saturated state. This also suggests that for unstable shear flow in systems without radiative damping |β2||β1| is expected, consistent with Ref. 13. The ky dependence of the ratio |β1/β2| roughly follows that of γ1, except that it approaches 1, rather than 0, at ky = 0 and ky=kxeq.

FIG. 10.

Error in ϕ of each of the eigenmode expansions of Fig. 9 relative to the filtered nonlinear data. Including ϕ2 significantly improves fluctuation estimates in the quasi-stationary state.

FIG. 10.

Error in ϕ of each of the eigenmode expansions of Fig. 9 relative to the filtered nonlinear data. Including ϕ2 significantly improves fluctuation estimates in the quasi-stationary state.

Close modal
FIG. 11.

Time-averaged ratio of unstable mode amplitude to stable mode amplitude in saturation at each ky for a range of Drad with DKrook = 1. The growth rate dependence of the Pt analysis35,49,50 suggests that higher Drad causes β2 to be driven less leading into saturation. Here we see that this is reflected in the eigenmode amplitudes in the saturated state. The ky dependence of the ratio |β1/β2| roughly follows that of γ1.

FIG. 11.

Time-averaged ratio of unstable mode amplitude to stable mode amplitude in saturation at each ky for a range of Drad with DKrook = 1. The growth rate dependence of the Pt analysis35,49,50 suggests that higher Drad causes β2 to be driven less leading into saturation. Here we see that this is reflected in the eigenmode amplitudes in the saturated state. The ky dependence of the ratio |β1/β2| roughly follows that of γ1.

Close modal

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

FIG. 12.

Time-averaged ratio of unstable to stable mode amplitude in saturation at each ky for a range of DKrook. Between DKrook = 0.5 and 4, increasing DKrook generally pushes the ratio of amplitudes closer to unity. Above DKrook=4, increasing DKrook has a less pronounced impact on the ratio. We stress that, in both regimes, the mean flow amplitude and Reynolds stress increase with DKrook.

FIG. 12.

Time-averaged ratio of unstable to stable mode amplitude in saturation at each ky for a range of DKrook. Between DKrook = 0.5 and 4, increasing DKrook generally pushes the ratio of amplitudes closer to unity. Above DKrook=4, increasing DKrook has a less pronounced impact on the ratio. We stress that, in both regimes, the mean flow amplitude and Reynolds stress increase with DKrook.

Close modal

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 ϕβ1ϕ1+β2ϕ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

ωj(kx2+ky2)ϕ̃j+kyV02[(kx22kxkxeq+ky2)ϕ̃j+(kx2+2kxkxeq+ky2)ϕ̃j+]=0.

Equation (4) can be expressed as a matrix equation ωjϕj+Mϕj=0 where the components of ϕj are ϕ̃j at different kx and the dimension of ϕj and M is infinite. Reasonable approximations of the eigenmodes and eigenvalues can be obtained by solving the matrix equation with ϕ̃j0 for some finite number of kx values, and ϕ̃j=0 for all other kx. This has previously been found useful in similar KH and tearing mode calculations.28 For example, solving the system with kx=0,±kxeq yields

ω1=ikyV0κϕ1=(1,iκ,1)T,
ω2=ikyV0κϕ2=(1,iκ,1)T,
ω3=0ϕ3=(1,0,1)T,

where the vectors are written in the form (ϕ̃(kxeq),ϕ̃(0),ϕ̃(kxeq))T and κ(ky)2((kxeq)2+ky2)/((kxeq)2ky2).

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

tϕ̃(kxeq,0)+kkykxeq((kxeqkx)2+ky2)ϕ̃(kx,ky)ϕ̃(kxeqkx,ky)=DKrookϕ̃(kxeq,0).

Considering a steady state where ϕ̃(kxeq,0)/t=0, and assuming ϕ̃=β1ϕ̃1+β2ϕ̃2 with just the kx=0,±kxeq Fourier modes considered, Eq. (12) can be manipulated to yield

ϕ̃(kxeq,0)=2ikxeqDKrookky>0kyκ(ky)(|β1|2|β2|2).
(15)

From here, values for |β1| and |β2| can be inserted to arrive at values for the mean flow amplitude. In other systems, |β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 |β1| and |β2| from nonlinear simulations into Eq. (15). Our interest is in the scaling of ϕ̃(kxeq,0) with DKrook, and what role β2 plays in that scaling. Thus, we perform simulations with a range of DKrook and compare three quantities: the time-averaged value of ϕ̃(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 Drad, corresponding to three systems with varying degrees of stable mode excitation (recall stable modes are more excited at lower values of Drad, see Fig. 11). For Drad = 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 Drad = 0.025, where stable modes are even more important, the β2 = 0 model fails to even capture the decrease in |ϕ̃(kxeq,0)| with DKrook, while the model that includes β2 does capture the correct qualitative and often even quantitative behavior. At Drad = 0.1, where stable modes are significantly weakened, their inclusion does not have a significant impact on the model. For each value of Drad, the two models are scaled by a constant so that they agree with the simulation results at DKrook = 4, which is where the change in scaling with respect to DKrook 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 ϕ̃ at noninteger multiples of kxeq.

FIG. 13.

A comparison between the true value of |ϕ(kxeq,0)| (blue dots) and values predicted by Eq. (15) with (orange diamonds) and without (green crosses) the stable mode contribution β2 included over multiple values of DKrook and Drad. In each frame, the models are scaled by a constant coefficient to match the true value at DKrook=4 so that the scaling with DKrook can be investigated, rather than the absolute agreement. For the base case Drad = 0.05 (center frame), the scaling of |ϕ(kxeq,0)| with DKrook is qualitatively captured by the model with β2 neglected; however, the scaling is significantly improved when stable modes are included. For the larger value of Drad (right frame), where stable modes are largely suppressed in saturation (cf. Fig. 11), including β2 in the model produces little change, and decent quantitative agreement is observed by the model both with and without β2 included. For the smaller value of Drad (left frame), where stable modes are more important, the model fails to even qualitatively agree with simulations unless β2 is included, in which case the overall trend is captured.

FIG. 13.

A comparison between the true value of |ϕ(kxeq,0)| (blue dots) and values predicted by Eq. (15) with (orange diamonds) and without (green crosses) the stable mode contribution β2 included over multiple values of DKrook and Drad. In each frame, the models are scaled by a constant coefficient to match the true value at DKrook=4 so that the scaling with DKrook can be investigated, rather than the absolute agreement. For the base case Drad = 0.05 (center frame), the scaling of |ϕ(kxeq,0)| with DKrook is qualitatively captured by the model with β2 neglected; however, the scaling is significantly improved when stable modes are included. For the larger value of Drad (right frame), where stable modes are largely suppressed in saturation (cf. Fig. 11), including β2 in the model produces little change, and decent quantitative agreement is observed by the model both with and without β2 included. For the smaller value of Drad (left frame), where stable modes are more important, the model fails to even qualitatively agree with simulations unless β2 is included, in which case the overall trend is captured.

Close modal

Comparing the two models at different values of Drad 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 |β1| and |β2| were obtained from gyrokinetic simulations with finite D, suggesting gyrokinetic effects may not play a significant role in the eigenmode excitations and force balance in this system.

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.

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.

1.
S. A.
Balbus
and
J. F.
Hawley
,
Rev. Mod. Phys.
70
,
1
(
1998
).
2.
R.
Vanon
and
G. I.
Ogilvie
,
Mon. Not. R. Astron. Soc.
463
,
3725
(
2016
).
3.
M.
Faganello
and
F.
Califano
,
J. Plasma Phys.
83
,
535830601
(
2017
).
4.
B. N.
Rogers
,
W.
Dorland
, and
M.
Kotschenreuther
,
Phys. Rev. Lett.
85
,
5336
(
2000
).
5.
S.
Chandrasekhar
,
Hydrodynamic and Hydromagnetic Stability
(
Oxford University Press
,
1961
).
6.
P. G.
Drazin
and
W. H.
Reid
,
Hydrodynamic Stability
(
Cambridge University Press
,
1981
).
7.
M.
Gaster
,
E.
Kit
, and
I.
Wygnanski
,
J. Fluid Mech.
150
,
23
(
1985
).
8.
M. L.
Palotti
,
F.
Heitsch
,
E. G.
Zweibel
, and
Y.-M.
Huang
,
Astrophys. J.
678
,
234
(
2008
).
9.
W. W.
Liou
and
P. J.
Morris
,
Phys. Fluids A
4
,
2798
(
1992
).
10.
D. E.
Nikitopoulos
and
J. T. C.
Liu
,
Phys. Fluids
13
,
966
(
2001
).
11.
W.
Horton
,
T.
Tajima
, and
T.
Kamimura
,
Phys. Fluids
30
,
3485
(
1987
).
12.
P.
Henri
,
S. S.
Cerri
,
F.
Califano
,
F.
Pegoraro
,
C.
Rossi
,
M.
Faganello
,
O.
Šebek
,
P. M.
Trávníček
,
P.
Hellinger
,
J. T.
Frederiksen
,
A.
Nordlund
,
S.
Markidis
,
R.
Keppens
, and
G.
Lapenta
,
Phys. Plasmas
20
,
102118
(
2013
).
13.
A. E.
Fraser
,
P. W.
Terry
,
E. G.
Zweibel
, and
M. J.
Pueschel
,
Phys. Plasmas
24
,
062304
(
2017
).
14.
P. W.
Terry
,
B. J.
Faber
,
C. C.
Hegna
,
V. V.
Mirnov
,
M. J.
Pueschel
, and
G. G.
Whelan
,
Phys. Plasmas
25
,
012308
(
2018
).
15.
C. C.
Hegna
,
P. W.
Terry
, and
B. J.
Faber
,
Phys. Plasmas
25
,
022511
(
2018
).
16.
G. G.
Whelan
,
M. J.
Pueschel
, and
P. W.
Terry
,
Phys. Rev. Lett.
120
,
175002
(
2018
).
17.
F.
Jenko
,
W.
Dorland
,
M.
Kotschenreuther
, and
B. N.
Rogers
,
Phys. Plasmas
7
,
1904
(
2000
).
18.
See http://www.genecode.org for code details and access.
19.
D. R.
Hatch
,
P. W.
Terry
,
F.
Jenko
,
F.
Merz
,
M. J.
Pueschel
,
W. M.
Nevins
, and
E.
Wang
,
Phys. Plasmas
18
,
055706
(
2011
).
20.
P. W.
Terry
,
K. D.
Makwana
,
M. J.
Pueschel
,
D. R.
Hatch
,
F.
Jenko
, and
F.
Merz
,
Phys. Plasmas
21
,
122303
(
2014
).
21.
B. N.
Rogers
and
W.
Dorland
,
Phys. Plasmas
12
,
062511
(
2005
).
22.
N.
Platt
,
L.
Sirovich
, and
N.
Fitzmaurice
,
Phys. Fluids A
3
,
681
(
1991
).
23.
S.
Musacchio
and
G.
Boffetta
,
Phys. Rev. E
89
,
023004
(
2014
).
24.
D.
Lucas
and
R.
Kerswell
,
J. Fluid Mech.
750
,
518
(
2014
).
25.
J.
Goodman
and
G.
Xu
,
Astrophys. J.
432
,
213
(
1994
).
26.
H. N.
Latter
,
P.
Lesaffre
, and
S. A.
Balbus
,
Mon. Not. R. Astron. Soc.
394
,
715
(
2009
).
27.
M. E.
Pessah
and
J.
Goodman
,
Astrophys. J.
698
,
L72
(
2009
).
28.
M. E.
Pessah
,
Astrophys. J.
716
,
1012
(
2010
).
29.
H. N.
Latter
,
S.
Fromang
, and
O.
Gressel
,
Mon. Not. R. Astron. Soc.
406
,
848
(
2010
).
30.
P.-Y.
Longaretti
and
G.
Lesur
,
Astron. Astrophys.
516
,
51
(
2010
).
31.
J.
Squire
,
E.
Quataert
, and
M. W.
Kunz
,
J. Plasma Phys.
83
,
905830613
(
2017
).
32.
J.-H.
Kim
and
P. W.
Terry
,
Phys. Plasmas
17
,
112306
(
2010
).
33.
M. J.
Pueschel
,
D.
Told
,
P. W.
Terry
,
F.
Jenko
,
E. G.
Zweibel
,
V.
Zhdankin
, and
H.
Lesch
,
Astrophys. J. Suppl. Ser.
213
,
30
(
2014
).
34.
J. B.
Marston
,
E.
Conover
, and
T.
Schneider
,
J. Atmos. Sci.
65
,
1955
(
2008
).
35.
P. W.
Terry
,
D. A.
Baver
, and
D. R.
Hatch
,
Phys. Plasmas
16
,
122305
(
2009
).
36.
A. J.
Brizard
and
T. S.
Hahm
,
Rev. Mod. Phys.
79
,
421
(
2007
).
37.
F.
Merz
, Ph.D. thesis,
University of Münster
,
2009
.
38.
M. J.
Pueschel
,
F.
Jenko
,
D.
Told
, and
J.
Büchner
,
Phys. Plasmas
18
,
112102
(
2011
).
39.
M. J.
Pueschel
,
D. R.
Hatch
,
T.
Görler
,
W. M.
Nevins
,
F.
Jenko
,
P. W.
Terry
, and
D.
Told
,
Phys. Plasmas
20
,
102301
(
2013
).
40.
K. J.
Burns
,
G. M.
Vasil
,
J. S.
Oishi
,
D.
Lecoanet
, and
B.
Brown
,
Astrophys. Source Code Library
,
2011
, http://ascl.net/1603.015.
41.
42.
M. J.
Pueschel
,
T.
Dannert
, and
F.
Jenko
,
Comput. Phys. Commun.
181
(
8
),
1428
1437
(
2010
).
43.
S. M.
Tobias
,
K.
Dagon
, and
J. B.
Marston
,
Astrophys. J.
727
,
12
(
2011
).
44.
J. M.
Reynolds-Barredo
,
D. E.
Newman
,
P. W.
Terry
, and
R.
Sanchez
,
Europhys. Lett.
155
,
34002
(
2016
).
45.
F.
Waleffe
,
Phys. Fluids
7
,
3060
(
1995
).
46.
M.
Kammerer
,
F.
Merz
, and
F.
Jenko
,
Phys. Plasmas
15
,
052102
(
2008
).
47.
K. M.
Case
,
Phys. Fluids
3
,
143
(
1960
).
48.
A.
Brandstäter
,
J.
Swift
,
H. L.
Swinney
,
A.
Wolf
,
J. D.
Farmer
,
E.
Jen
, and
P. J.
Crutchfield
,
Phys. Rev. Lett.
51
,
1814
(
1983
).
49.
P. W.
Terry
,
D. A.
Baver
, and
S.
Gupta
,
Phys. Plasmas
13
,
022307
(
2006
).
50.
K. D.
Makwana
,
P. W.
Terry
,
J.-H.
Kim
, and
D. R.
Hatch
,
Phys. Plasmas
18
,
012302
(
2011
).
51.
M. J.
Pueschel
,
B. J.
Faber
,
J.
Citrin
,
C. C.
Hegna
,
P. W.
Terry
, and
D. R.
Hatch
,
Phys. Rev. Lett.
116
,
085001
(
2016
).
52.
B. J.
Faber
,
M. J.
Pueschel
,
P. W.
Terry
,
C. C.
Hegna
, and
J. E.
Roman
,
J. Plasma Phys.
84
,
905840503
(
2018
).
53.
D. R.
Hatch
,
F.
Jenko
,
A.
Bañón Navarro
,
V.
Bratanov
,
P. W.
Terry
, and
M. J.
Pueschel
,
New J. Phys.
18
,
075018
(
2016
).
54.
D. A.
Baver
,
P. W.
Terry
,
R.
Gatto
, and
E.
Fernandez
,
Phys. Plasmas
9
,
3318
(
2002
).
55.
P. W.
Terry
and
R.
Gatto
,
Phys. Plasmas
13
,
062309
(
2006
).
56.
S. E.
Clark
, Ph.D. thesis,
Columbia University
,
2017
.