Visco-resistive magnetohydrodynamic turbulence, driven by a two-dimensional unstable shear layer that is maintained by an imposed body force, is examined by decomposing it into dissipationless linear eigenmodes of the initial profiles. The down-gradient momentum flux, as expected, originates from the large-scale instability. However, continual up-gradient momentum transport by large-scale linearly stable but nonlinearly excited eigenmodes is identified and found to nearly cancel the down-gradient transport by unstable modes. The stable modes effectuate this by depleting the large-scale turbulent fluctuations via energy transfer to the mean flow. This establishes a physical mechanism underlying the long-known observation that coherent vortices formed from nonlinear saturation of the instability reduce turbulent transport and fluctuations, as such vortices are composed of both the stable and unstable modes, which are nearly equal in their amplitudes. The impact of magnetic fields on the nonlinearly excited stable modes is then quantified. Even when imposing a strong magnetic field that almost completely suppresses the instability, the up-gradient transport by the stable modes is at least two-thirds of the down-gradient transport by the unstable modes, whereas for weaker fields, this fraction reaches up to 98%. These effects are persistent with variations in magnetic Prandtl number and forcing strength. Finally, continuum modes are shown to be energetically less important, but essential for capturing the magnetic fluctuations and Maxwell stress. A simple analytical scaling law is derived for their saturated turbulent amplitudes. It predicts the falloff rate as the inverse of the Fourier wavenumber, a property which is confirmed in numerical simulations.

Owing to their ubiquity in laboratory,1 geophysical,2,3 and astrophysical environments,4–9 shear layers have been extensively studied.10,11 Observations and analyses from experiments and direct numerical simulations have offered insights into the connection between large-scale vortical structures formed from the instability of a shear layer and turbulent transport across the layer.12–14 Properties like shape and scale of the nonlinearly saturated vortices, which dominate the transport, are generally attributed to the linearly unstable eigenmodes or closely related nonlinear fluctuations.15,16 The nonlinear saturation of the instability, however, can be more complex than just the finite-amplitude modifications of unstable modes, as emerging understanding in fusion plasma instability demonstrates.17–30 

Already in the late 1960s, using one of the early numerical simulations of shear instability,31 it was hinted that the nonlinear saturation of Kelvin–Helmholtz instability involves, contrary to finite-amplitude modifications of unstable modes, quasi-periodic oscillations in the fluctuations. Later, an intuitive understanding of how such a phenomenon occurs in sheared fluids32 has been reported by invoking vortex nutation:15 fluctuation-amplitude oscillations correlate with oscillations in the mean flow energy and lead to vortex nutation. Fluctuations, however, are usually not decomposed into the complete set of linear eigenmodes and are commonly assumed15 to be due to unstable mode structures. However, since unstable modes always drive a down-gradient momentum transport, they cannot explain the increase in kinetic energy in the mean flow.

Notably, occasional up-gradient momentum transport has been observed in several experimental and numerical studies where an unstable shear layer drives the turbulence.33–37 Analyses of these transient events12,38–41 do not address the underlying conditions producing this dynamics—whether the transient up-gradient transport is a part of an ongoing subdominant process with occasional breakthroughs or simply spontaneous fluctuations. The laboratory and prior numerical experiments alone are not sufficient to definitively answer this question. One way to expose the underlying process is to examine the turbulent fluctuations using a complete eigenmode decomposition and to assign roles and activities to each mode in the transport phenomena. Indeed, there can be modes other than the unstable modes that are important in the turbulent phase, as an insightful study hints: the dominant vortex in a turbulent background orients quasi-periodically against (or toward) the mean flow and drives the down-gradient (or up-gradient) momentum transport.12 To understand such behaviors in detail, it is instructive to also analyze how the instability saturates, a question that has long been of interest42 but for which understanding remains incomplete.

When turbulence is sustained via continuous energy injection from a large-scale instability, there exist two primary candidates for instability saturation. A common (but not necessarily justified) assumption is that energy injected by the instability is transferred conservatively to increasingly smaller scales in a forward, Kolmogorov-like cascade, where nonlinear interactions move energy between linearly unstable or marginal modes until a dissipation range is reached at small scales.5 An alternative process involves linearly stable eigenmodes at the large injection scales, which absorb and remove significant energy from scales that launch the inertial cascade. In several studies of microturbulence in fusion plasmas, linearly stable modes have been found to be excited to significant levels via nonlinear interactions and to drastically affect the saturated amplitudes and transport characteristics of the system.17–22,26–30

Stable modes in shear flow turbulence, however, have been studied only recently,43–45 and more remains to be understood, e.g., their role in mixing and magnetic field evolution and how they might affect reduced models of turbulence and transport. It was predicted in Ref. 43 that the Kelvin–Helmholtz instability in its nonlinear evolution excites a linearly stable conjugate-root46 of the inviscid instability, which affects the instability saturation even when viscosity is finite. This was later verified in numerical simulations of freely evolving shear layers.44 However, the rapid relaxation of the layer toward a stable profile on a time scale similar to that of stable-mode excitation prevented general conclusions from being reached, regarding how the turbulence and transport are affected by the stable modes. The issue is aggravated by the addition of a flow-aligned magnetic field, which causes the layer to relax even more rapidly. To circumvent this challenge, one may drive the mean flow toward the unstable profile and, thus, achieve a thorough statistical quantification of the stable modes. Note that driven profiles are quite common in astrophysical shear flows, with forces like gravity providing free energy for the drive. For these reasons, driven shear flow is studied here.

The principal result of this study is that significant up-gradient momentum transport is driven by nonlinearly excited (linearly-)stable modes, canceling a substantial portion of the down-gradient transport by unstable modes, and notably, this transport is present not just during turbulent momentum flux reversals, but is continuously at work at a slightly lower level than that of the unstable modes. This finding is robust even for variations of orders of magnitude in background magnetic field strength, magnetic Prandtl number (or resistivity), and forcing strength of the mean flow. Note that the stronger background magnetic field tends to suppress the instability46 and disrupt the large-scale vortices,47 while the larger magnetic Prandtl number (weaker resistivity for a fixed viscosity) extends the scale range of magnetic fluctuations, compared to the flow fluctuations.48 We also show, for astrophysical applications, that a turbulent viscosity can be defined, with the addition of stable modes, that can reliably capture the Reynolds stress: without stable modes, however, the stresses are greatly over-predicted by the unstable modes.

This article is organized in the following manner. Section II entails the magnetohydrodynamic (MHD) model of the shear flow and details the system setup. In Sec. III, the complete linear eigenspectrum is presented, along with a discussion on the roles of different eigenmodes. Section IV shows the full nonlinear evolution of MHD Kelvin–Helmholtz instability using direct numerical simulations. A decomposition of the turbulent fluctuations onto linear eigenmodes is performed in Sec. V, where a detailed study of imprints of stable modes in turbulence and transport is presented. Section VI summarizes the findings of this work.

An incompressible magneto-fluid is considered in this study, and standard MHD equations are adopted,

·u=0,
(1a)
tu+u·u=Pρ+(×B)×B4πρ+ν2u+f,
(1b)
·B=0,
(1c)
tB=×(u×B)+η2B,
(1d)

where u, B, P, ρ, ν, η, and f, respectively, denote the fluid velocity, magnetic field, pressure, fluid density, viscosity, Ohmic diffusivity, and externally supplied acceleration to the magneto-fluid.

A shear layer is examined on a two-dimensional (x, z) plane with the initial fluid velocity given by u(x,z,t=0)=U0tanh(z/a)x̂ and a flow-aligned magnetic field, initially uniform, as B(x,z,t=0)=B0x̂. The parameters a,U0,andB0 signify the half-width of the flow-shear, maximum initial fluid velocity, and initial magnetic field, respectively. These parameters are exploited to non-dimensionalize all the variables henceforth. Length, time, and energy (per unit mass) are hereafter measured in units of a, a/U0, and U02, respectively. Thus, the initial (or reference) flow and magnetic field are represented by Uref(z)=tanh(z) and Bref(z)=1 in the rest of this article. The ratio of the maximum fluid speed U0 to the Alfvèn speed can be written as the Alfvènic Mach number MA=U04πρ/B0. The viscosity and resistivity are quantified via fluid Reynolds number Re=U0a/ν and magnetic Reynolds number Rm=U0a/η, respectively.

In two dimensions, a more convenient formalism is available, using the streamfunction ϕ and flux function ψ. Defining u=ŷ×ϕ and B=ŷ×ψ, the vorticity and the current become 2ϕŷ and 2ψŷ, respectively. Taking the curl of Eq. (1b) and rewriting Eq. (1d) in terms of the stream- and flux-functions yield49 

t2ϕ+{2ϕ,ϕ}=MA2{2ψ,ψ}+Re14ϕ+zf(kx=0,z,t),
(2a)
tψ={ϕ,ψ}+Rm12ψ,
(2b)

where the Poisson bracket is {P,Q}=xP·zQzP·xQ, e.g., {ϕ,ψ}=u·ψ. Here, kx is the Fourier wavenumber along the x-axis. The parameters Re=Rm=500 are chosen for all simulations unless mentioned otherwise (where Rm is changed to 50 and 5000 in different simulations). It should be emphasized that these Reynolds numbers are defined using the initial scale a of the sharpest gradient in the flow as the characteristic length scale; however, as the system evolves nonlinearly via vortex merging, despite the forced mean flow, eddies of the size of the simulation box appear, which may be considered as the characteristic length scale of motion.11 When choosing this normalization, non-dimensional numbers should be scaled accordingly, e.g., Rm=5000 becomes Rm=5000×Lx1.5×105, where Lx represents the box-size along the mean flow direction. The external body force, f=f(kx=0,z,t)x̂, is applied to the mean flow only, which is highlighted in Eq. (2a) using the explicit mention of kx=0. As in a recent study,45 the forcing drives the instantaneous mean flow toward the initial unstable profile Uref(z). A similar forcing mechanism exists for geo- and astrophysical flows where gravitation50 tends to build shear layers. We assume here such a force, represented as a Krook operator,51–53 as

f=DKrook[Uref(z)U(x,z,t)x]+F0,
(3)

where DKrook, sometimes also referred to as the profile relaxation rate,54 measures the forcing strength (in units of U0/a) and U(x,z,t)x represents the instantaneous x-averaged flow. If DKrook=0, the shear layer evolves freely, and decaying turbulence is realized as a result of the Kelvin–Helmholtz instability and the turbulence it generates.

The time-independent force F0 is implemented only to balance the viscous diffusion of the initial shear layer: Re12Uref(z)+F0=0 ensures an initial equilibrium state to which small-amplitude perturbations are added before the system is evolved.

As in the unforced study,44 a simulation box of Lx=10π is considered, but double the size along the z-axis (Lz=20π), given that the quasi-stationary turbulence simulated herein is run for much longer time, which tends to create fully developed turbulent features that are larger in size. Thus, we adopt a larger domain to minimize their potential interactions with the boundaries in the z-axis. Note the forcing applied to the mean flow prevents profile relaxation, and the turbulence remains mostly in the vicinity of the shear layer. The numerical code Dedalus,55 a pseudospectral solver, is used in this study. Fourier modes along the x-axis and Chebyshev polynomials along the z-axis are employed with (Nx,Nz)=(2048,2048) spectral modes. We confirmed that the spectral energy density and dissipation are converged at this resolution. Note also that these high resolutions benefit the eigenmode projection of nonlinear data in the post-processing analysis. Only for the simulation with magnetic Prandtl number of 10, the box size was changed to (Lx,Lz)=(6π,8π), and the resolution was increased to (Nx,Nz)=(4096,4096); the same simulation was repeated with (Nx,Nz)=(4096,8192), but only for early times due to computational cost and found to reproduce, among others, the energy evolution. All simulations use 3/2 dealiasing rule, additionally. The boundary conditions used in all simulations are periodic along the x-axis, and perfectly conducting, no-slip, co-moving (with the initial flow) at the top and bottom boundaries, z=±Lz/2.44,45

The initial equilibrium state is seeded with small-amplitude perturbations (ϕ̃,ψ̃) at all Fourier wavenumbers, as44 

ϕ̃(x,z,t=0)=Aϕkx0kxaeirϕ(kx)ez2/σ2eikxx,
(4)

and

ψ̃(x,z,t=0)=Aψkx0kxaeirψ(kx)ez2/σ2eikxx.
(5)

Here, Aϕ and Aψ set the overall amplitudes of the perturbations that have a Gaussian width controlled by σ and the rate at which they falloff with the wavenumbers given by a. The random phases rϕ(kx) and rψ(kx), forming a uniform distribution in [0,2π), are issued for each different kx using a pseudo-random number generator. Different choices of these initial conditions were investigated in Ref. 44, motivating the choice here: a = −1, σ = 2, and Aϕ=Aψ=103. This set of parameters offers distinct linear and nonlinear phases of evolution.

Aiming to understand the nonlinear excitation of linear eigenmodes in the turbulent phase, first the nonlinear initial-value problem is solved to collect high-fidelity turbulent data. Afterward, a separate eigenvalue problem is solved to obtain a complete linear eigenspectrum and eigenmodes, which are used to expand the nonlinear data on this basis to track the amplitude of each eigenmode. Such a basis is obtained by linearizing the governing equations around the initial flow and magnetic field profiles, by dropping the dissipative terms. The eigenmodes, thus, obtained are of a dissipationless linear operator. Of course, the meaning and utility of this linear basis are a priori unknown. Nevertheless, when a basis forms a complete set, one can always expand an arbitrary fluctuation on that basis. As the non-dissipative equations of motion preserve Parity-Time (PT-) reversal symmetry, such a system is theoretically guaranteed to yield a complete basis as established recently in PT-symmetric quantum mechanics.56 Previous studies in gyrokinetic and MHD plasmas have also revealed the usefulness of dissipationless linear eigenmodes in interpreting dissipative nonlinear systems.27,44,57

With the intent of obtaining dissipationless linear eigenmodes, the variables (ϕ,ψ) in Eqs. (2a) and (2b) are decomposed into background and perturbations, (ϕ,ψ)=(ϕref,ψref)+(ϕ̃,ψ̃). The linearized dissipationless equations for the evolution of perturbations are

t2ϕ̃=[Urefx2(z2Uref)·x]ϕ̃+1MA2[Brefx2(z2Bref)·x]ψ̃,
(6a)
tψ̃=Urefxψ̃+Brefxϕ̃.
(6b)

Fourier transforming along the x-axis and assuming time variation at each Fourier wavenumber take the form ϕ̂(kx,z,ω)eiω(kx)t, Eqs. (6a) and (6b) become

ω(z2kx2)ϕ̂=kx[Uref(z2kx2)(z2Uref)]ϕ̂+1MA2kx[Bref(z2kx2)(z2Bref)]ψ̂,
(7a)
ωψ̂=Urefkxψ̂+Brefkxϕ̂.
(7b)

Solving Eqs. (7a) and (7b), the eigenvalues ω are found to be real except when 0<|kx|<1, where two of the real eigenvalues coalesce to produce imaginary eigenvalues,58 as complex conjugate to each other. These are the growth rates of the unstable eigenmode and its conjugate stable eigenmode, which evolve in time as eγ(kx)t and eγ(kx)t, respectively. This mode-pair is shown, for the first Fourier wavenumber kx=2π/Lx=0.2, in Fig. 1(a), along with all the purely real eigenvalues. The latter constitute the eigenmode continuum59 and are theoretically infinite in number, although numerical discretization always yields a finite but very large number of modes (>3000 for each wavenumber in this study). These eigenvalues are given by the relation ω/kx+Uref(z)±vA,ref(z)=0, where vA,ref(z) is the Alfvèn speed along the reference magnetic field at the vertical coordinate z.

FIG. 1.

(a) Linear eigenspectrum of the MHD shear-flow system (at kx=0.2 with MA=10). The unstable and stable modes are shown with thick (colored) crosses. Among the continuum modes that form a vertical line, a zero-frequency continuum mode is displayed with a green-colored star. (b)–(d) Eigenfunctions in z-space, with real (Re) and imaginary (Im) parts, for unstable (ϕ1), stable (ϕ2), and one continuum (ω = 0) mode. (e)–(g) Corresponding eigenmode structures in (x, z) space. Note that the eigenmodes ϕ1 and ϕ2 are complex conjugate to each other. Imaginary parts in their eigenfunctions induce relative tilt between them in (x, z) space, which will be consequential for momentum transport in Sec. III B. Each eigenmode is normalized to have unit total energy [following which the maximum values of ϕ in (b)–(g) are chosen].

FIG. 1.

(a) Linear eigenspectrum of the MHD shear-flow system (at kx=0.2 with MA=10). The unstable and stable modes are shown with thick (colored) crosses. Among the continuum modes that form a vertical line, a zero-frequency continuum mode is displayed with a green-colored star. (b)–(d) Eigenfunctions in z-space, with real (Re) and imaginary (Im) parts, for unstable (ϕ1), stable (ϕ2), and one continuum (ω = 0) mode. (e)–(g) Corresponding eigenmode structures in (x, z) space. Note that the eigenmodes ϕ1 and ϕ2 are complex conjugate to each other. Imaginary parts in their eigenfunctions induce relative tilt between them in (x, z) space, which will be consequential for momentum transport in Sec. III B. Each eigenmode is normalized to have unit total energy [following which the maximum values of ϕ in (b)–(g) are chosen].

Close modal

The eigenfunctions, normalized to have unit total energy, are also shown in Fig. 1: along the z-axis, see Figs. 1(b)–1(d), and in (x, z) space, see Figs. 1(e)–1(g). Note that complex conjugation transforms the unstable mode ϕ1(kx,z) into the stable mode ϕ2(kx,z) and vice versa. This is a direct consequence of spontaneous PT-symmetry breaking in the ideal shear-flow instability.58 (The spontaneous symmetry breaking does not imply that the equation of motion or the associated Hamiltonian breaks PT-symmetry; it is rather some of the eigenfunctions of such a PT-symmetry-preserving Hamiltonian that break PT-symmetry.)

A representative eigenfunction of a continuum mode, shown in Fig. 1(d), exhibits sharp and narrow structure. To what physics each type of eigenmode structure contributes will be explored in this article.

Shown in Fig. 2 is a schematic diagram, illustrating how the relative tilts in the eddies can transport momentum in opposite directions across the shear layer.14 It can be qualitatively observed from Figs. 1(e) and 1(f) that the unstable and stable modes drive down- and up-gradient momentum transport, respectively. Precise quantitative measures will be built and computed later in Sec. V D.

FIG. 2.

(a) The unstable mode of the flow transport momentum in the down-gradient direction: –x-directed momentum at A is carried to A′ and +x-directed momentum at B is carried to B′. Fluxes A A′ and B B′ act to relax the mean flow gradient (shown with the long horizontal arrows). (b) Oppositely tilted eddies, which correspond to a stable mode, transport momentum in the up-gradient direction: –x-directed momentum at C is carried to C′ and +x-directed momentum at D is carried to D′. Both of these fluxes replenish the mean flow. The direction of the streamlines (shown with gray arrows on the elliptic eddies) does not alter these properties, but the tilt does.

FIG. 2.

(a) The unstable mode of the flow transport momentum in the down-gradient direction: –x-directed momentum at A is carried to A′ and +x-directed momentum at B is carried to B′. Fluxes A A′ and B B′ act to relax the mean flow gradient (shown with the long horizontal arrows). (b) Oppositely tilted eddies, which correspond to a stable mode, transport momentum in the up-gradient direction: –x-directed momentum at C is carried to C′ and +x-directed momentum at D is carried to D′. Both of these fluxes replenish the mean flow. The direction of the streamlines (shown with gray arrows on the elliptic eddies) does not alter these properties, but the tilt does.

Close modal

Since the unstable and stable modes compete with each other to transport momentum in opposing directions, the excitation levels of these modes are crucial. In the linear phase of instability evolution, the transport by the unstable modes dominates over the transport by the stable modes. However, this need not be the case in the nonlinear phase, as nonlinear processes can excite the stable modes to appreciable levels. Whenever the stable modes surpass the unstable modes in amplitudes, net momentum is transported in the up-gradient direction.44,60 In extremely simplified models of transport, such as eddy viscosity models, this contributes to negative eddy viscosity. Computing the amplitude of each eigenmode in the nonlinear phase can, thus, be helpful to build improved reduced transport models. A recent investigation also demonstrated that this kind of competition between the two large-scale eigenmodes alters the magnetic cascade substantially.45 

Having provided a description of linear eigenmodes, we now turn to properties of the nonlinear system, before discussing how expressions of linear modes may be identified in turbulent fluctuations.

Small-amplitude perturbations in the flow and magnetic field evolve exponentially fast in the linear regime of the instability, giving rise to a chain of spiral vortices, as evident in Figs. 3(a) and 3(d). These structures then interact nonlinearly with nearby vortices to yield even larger turbulent structures as in Figs. 3(b) and 3(c). A contrast is to be made between forced and unforced simulations. In the latter, the gradient of the mean flow flattens out as the instability extracts energy. Decaying turbulence then ensues. Forcing the mean flow, however, leads to a quasi-stationary turbulence, as the energy in the gradient is replenished with the instability drawing on its energy. In the saturated stage, energy input through the unstable modes is balanced by energy removal via stable modes as well as dissipative channels.

FIG. 3.

Time evolution of vorticity in (a)–(c) in a simulation with a forced background flow, DKrook=2, and (d)–(f) in a simulation with a freely evolving shear layer, DKrook=0; both for MA=10. Panels (a)–(c) share the same colorbar, and (d)–(f) share another colorbar. The instantaneous mean flow profile in each of the subplots is shown with a black dotted curve, where the vertical axis represents the z-coordinate and the horizontal direction corresponds to the x-velocity U0(z,t), as exemplified in the inset of (a). Two arrows pointing in opposite directions show the direction of the flow in the regions z > 0 and z < 0. The initial flow profile U0(z,t=0)=tanh(z) is shown with a red dashed curve in (a) and (d). Rapid profile flattening is evident in (d). While the instability dies out in the unforced case, quasi-stationary turbulence is realized in the forced case in (c).

FIG. 3.

Time evolution of vorticity in (a)–(c) in a simulation with a forced background flow, DKrook=2, and (d)–(f) in a simulation with a freely evolving shear layer, DKrook=0; both for MA=10. Panels (a)–(c) share the same colorbar, and (d)–(f) share another colorbar. The instantaneous mean flow profile in each of the subplots is shown with a black dotted curve, where the vertical axis represents the z-coordinate and the horizontal direction corresponds to the x-velocity U0(z,t), as exemplified in the inset of (a). Two arrows pointing in opposite directions show the direction of the flow in the regions z > 0 and z < 0. The initial flow profile U0(z,t=0)=tanh(z) is shown with a red dashed curve in (a) and (d). Rapid profile flattening is evident in (d). While the instability dies out in the unforced case, quasi-stationary turbulence is realized in the forced case in (c).

Close modal

It is now timely to discuss the turbulent transport of momentum in nonlinear simulations. To derive the turbulent stresses, the evolution equation of the mean flow can be written by x-averaging the momentum equation,

tUx=z(τu+τb)+DKrook[Uref(z)Ux]+1Re2z2[UxUref(z)],
(8)

where U=U(x,z,t) represents the instantaneous flow, ·x signifies x-averaging operation, and τu and τb are the Reynolds and Maxwell stresses, arising from the correlations of turbulent fluctuations of the flow and the magnetic fields, respectively. Note that in Fraser et al.,44 a negative sign was typographically missed in front of the first term on the right-hand side of Eq. (8). With the sign displayed in Eq. (8) above, the turbulent stresses are given by

τu=ũxũzx=zϕ̃·xϕ̃x,
(9a)
τb=1MA2b̃xb̃zx=1MA2zψ̃·xψ̃x.
(9b)

These stresses are evaluated from nonlinear simulations and shown in Fig. 4. Fluctuations of Reynolds stress are concentrated in the shear layer, near z0. Time histories of the Reynolds and Maxwell stresses, at z =0, where they are largest in magnitude, are compared in Figs. 4(b) and 4(d). Note the recurring dominant up-gradient transport via the Reynolds stress. The Maxwell stress, however, is almost always down-gradient. Figures 4(a) and 4(c) also convey that the Maxwell stress is generally broader along the z-axis than the Reynolds stress, which is more localized near the shear layer.

FIG. 4.

Time evolution of MHD stresses. (a) Reynolds stress τu(t,z). (b) Reynolds stress τu(t,z=0) at the middle of the layer, z =0. (c) Maxwell stress τb(z,t). (d) Maxwell stress τb(t,z=0) at z =0. All data shown are for a single simulation with MA=10,DKrook=2. The Reynolds stress in (b) reverses several times, in contrast to the Maxwell stress in (d), which is almost always down-gradient.

FIG. 4.

Time evolution of MHD stresses. (a) Reynolds stress τu(t,z). (b) Reynolds stress τu(t,z=0) at the middle of the layer, z =0. (c) Maxwell stress τb(z,t). (d) Maxwell stress τb(t,z=0) at z =0. All data shown are for a single simulation with MA=10,DKrook=2. The Reynolds stress in (b) reverses several times, in contrast to the Maxwell stress in (d), which is almost always down-gradient.

Close modal

To probe the nonlinear simulation data, the turbulent fluctuations are expanded on the linear eigenmode basis described in Sec. III. Consider an arbitrary turbulent fluctuation χ̃turb=(ϕ̃turb,ψ̃turb), which is expanded as

χ̃turb(x,z,t)=kx0eikxxjβj(kx,t)χj(kx,z),
(10)

where the eigenmode basis χj(kx,z) is employed along the z-axis at each wavenumber kx to decompose the fluctuations. The complex mode-amplitude βj(kx,t), defined for each eigenmode j, can then be computed using properties of the linear operator, described in  Appendix A, even when the eigenmodes of the operator are non-orthogonal, as is the case here.

Following earlier studies,27,29,43–45j =1 and 2 will be used to represent unstable and stable modes, respectively. The computations herein resolve as many as 3109 eigenmodes at a particular kx.

The amplitudes of the unstable and stable modes are tracked in the nonlinear simulations, and their time series are plotted in Fig. 5(a). As expected, the unstable mode grows, and the stable mode decays exponentially in the early phase. However, as the fluctuations increase due to the growth of the unstable modes, nonlinear interactions among them begin exciting the linearly stable mode,30 causing it to rise to almost the same level as the unstable mode at that wavenumber, see Fig. 5(a). Later, in the fully nonlinear stage, all eigenmodes can participate in the energy redistribution.

FIG. 5.

(a) Time traces of the eigenmode amplitudes are shown for kx=0.2 for a simulation with MA=10 and DKrook=2. Note, in the inset, the nonlinear excitation of linearly stable mode (|β2|) in instability saturation (|β1| is the unstable mode amplitude). (b) All 3109 eigenmodes at kx=0.2 are plotted with their squared excitation levels in the nonlinearly saturated phase, which represent the energy in each eigenmode. The diameter of each circle shown corresponds to the energy in each eigenmode, and modes with lower energy are plotted on top of more highly excited modes, such that all data points are (partially) visible. Note that the total fluctuation energy is composed of both the modal and non-modal energy because of the non-orthogonal modes. Evaluating total energy at a wavenumber, E=dz[|u|2+|B|2/MA2]/2=dz[(mβmum)·(nβnun)+(mβmBm)·(nβnBn)/MA2]/2=m,nEmn, where (um,Bm) represents the mth eigenmode. When m and n belong to discrete (d) modes, Edd is, upon time-averaging (t=1501000), around 72% of the total energy, whereas when m and n belong to continuum (c) modes, Ecc is 22%; Edc is 6%.

FIG. 5.

(a) Time traces of the eigenmode amplitudes are shown for kx=0.2 for a simulation with MA=10 and DKrook=2. Note, in the inset, the nonlinear excitation of linearly stable mode (|β2|) in instability saturation (|β1| is the unstable mode amplitude). (b) All 3109 eigenmodes at kx=0.2 are plotted with their squared excitation levels in the nonlinearly saturated phase, which represent the energy in each eigenmode. The diameter of each circle shown corresponds to the energy in each eigenmode, and modes with lower energy are plotted on top of more highly excited modes, such that all data points are (partially) visible. Note that the total fluctuation energy is composed of both the modal and non-modal energy because of the non-orthogonal modes. Evaluating total energy at a wavenumber, E=dz[|u|2+|B|2/MA2]/2=dz[(mβmum)·(nβnun)+(mβmBm)·(nβnBn)/MA2]/2=m,nEmn, where (um,Bm) represents the mth eigenmode. When m and n belong to discrete (d) modes, Edd is, upon time-averaging (t=1501000), around 72% of the total energy, whereas when m and n belong to continuum (c) modes, Ecc is 22%; Edc is 6%.

Close modal

The energy in individual eigenmodes |βj|2, averaged over a turbulent state (t=1501000), is displayed in Fig. 5(b). It is evident that the unstable and stable eigenmode pair contains a majority (>70%) of the energy in the system. The remaining eigenmodes share a wide spectrum of the remaining energy. This suggests that the turbulent system at hand may be amenable to a substantial dimensionality reduction.19,27 For the cases of the weaker magnetic fields, this finding is more prominent, as evidenced in  Appendix B. In addition, the weaker fields support more coherent amplitude-oscillations, unlike the large excursions in the amplitudes observed with the stronger fields, e.g., MA=10 in Fig. 5(a). In the latter case, the stronger Lorentz back-reaction acting on the large-scale turbulent flow causes strong oscillations in the eigenmode amplitudes.

To obtain a better understanding of turbulent dynamics, it is of interest to compare different components of eigenmodes in the turbulent flow. An approximate (reduced) representation of the turbulent flow ϕ̃approx can be constructed from a class of eigenmodes at each wavenumber, e.g., ϕ̃approx(x,z,t) can be written as a sum of an unstable mode per wavenumber kxeikxxβ1(kx,t)ϕ1(kx,z) or as a sum of an unstable and a stable mode per wavenumber kxeikxx[β1(kx,t)ϕ1(kx,z)+β2(kx,t)ϕ2(kx,z)]. Respective short-hand notations β1ϕ1 and β1ϕ1+β2ϕ2 will be used hereafter, i.e.,

β1ϕ10<|kx|<1eikxxβ1(kx,t)ϕ1(kx,z),
(11a)
β1ϕ1+β2ϕ20<|kx|<1eikxx[β1(kx,t)ϕ1(kx,z)+β2(kx,t)ϕ2(kx,z)].
(11b)

The nonlinear fluctuations of the flow are compared in Fig. 6, viewed at different levels of truncation in the eigenmode expansion. The leftmost panel, Fig. 6(a), shows the full turbulent fluctuations in the Kelvin–Helmholtz (KH-)unstable wavenumbers kx=0.2,0.4,0.6,and0.8, which appear similar to the full turbulent fluctuations that include all wavenumbers in the nonlinear simulation (not shown); Fig. 6(b) displays the sum of unstable eigenmodes at each of these KH-unstable wavenumbers; and Fig. 6(c) presents the sum of unstable and stable eigenmodes at the same wavenumbers, while omitting all continuum modes. Adding stable modes produces a substantial improvement in the reconstruction. Note that such a reconstruction was found to deteriorate quickly over time (i.e., a few instability e-folding times where one e-folding time for the fastest growing mode kx=0.4 is γ15) in the study of unforced shear layers,44 as the rapid relaxation of the layer toward a stable profile rendered the unstable and stable eigenmodes of the system to be less representative of the decaying turbulence. The turbulent fluctuation shown in Fig. 6 is at t =702, which lies well within the nonlinear phase (the linear phase ends around t30). In this respect, the forced shear layer is markedly different from the freely evolving layer.

FIG. 6.

(a) Full turbulent fluctuations in streamfunction in the Kelvin–Helmholtz-unstable wavenumbers 0<|kx|<1 (thus called ϕ̃filtered), observed in nonlinear simulations with MA=30,DKrook=2. The shown plot of fluctuations includes all types of eigenmodes—the unstable, stable, and continuum modes. (b) Reconstruction of the turbulent fluctuations by summing only the unstable modes at the same wavenumber range. (c) Reconstruction by adding stable and unstable modes, while omitting all continuum modes. The reconstruction in (c) is clearly much alike the turbulent fluctuations in (a), in contrast to the reconstruction in (b). Saturation theory of instability that considers the unstable modes only, at best, can produce (b), but with the inclusion of the stable modes, substantial improvement can be achieved.

FIG. 6.

(a) Full turbulent fluctuations in streamfunction in the Kelvin–Helmholtz-unstable wavenumbers 0<|kx|<1 (thus called ϕ̃filtered), observed in nonlinear simulations with MA=30,DKrook=2. The shown plot of fluctuations includes all types of eigenmodes—the unstable, stable, and continuum modes. (b) Reconstruction of the turbulent fluctuations by summing only the unstable modes at the same wavenumber range. (c) Reconstruction by adding stable and unstable modes, while omitting all continuum modes. The reconstruction in (c) is clearly much alike the turbulent fluctuations in (a), in contrast to the reconstruction in (b). Saturation theory of instability that considers the unstable modes only, at best, can produce (b), but with the inclusion of the stable modes, substantial improvement can be achieved.

Close modal

While the qualitative analysis of the turbulent-flow reconstruction in Sec. V B is instructive, a quantitative measurement is desirable. To this end, following Ref. 44, the reconstructive capability of reduced representations is quantified, at each time step in the simulation, using the standard energy norm that measures the fraction of kinetic energy lost when the eigenmode basis is truncated, compared to the kinetic energy in the full turbulent flow data—see the definition in Eq. (12). The energy norm is well-suited for studying large-scale structures. Small-scale phenomena, however, may not be amenable to such analysis, although one may be able to find ties between the small- and large-scale phenomena in some cases. This measure is also called a “truncation error.” Note that this error arises not in the nonlinear simulations but merely in the reduced representations of turbulent fluctuations, when truncating the eigenmode basis in post-processing analyses.

Using the energy norm, we define the relative truncation error, which may also be called a normalized residual, in the following manner:

Residual=||ϕ̃exactϕ̃approx||KE2||ϕ̃exact||KE2=||ϕ̃diff||KE2||ϕ̃exact||KE2=dxdz[(xϕ̃diff)2+(zϕ̃diff)2]dxdz[(xϕ̃exact)2+(zϕ̃exact)2],
(12)

where (xϕ)2 and (zϕ)2 are the squared z- and x-components of velocities; ϕ̃diff=ϕ̃exactϕ̃approx with ϕ̃exact and ϕ̃approx representing, respectively, the turbulent streamfunction from nonlinear simulation and its reduced representation—either a summation over the unstable modes alone or over the unstable and stable modes together—both spanning fluctuations over a range of wavenumbers. Here, this range, taken to be the same for both, is considered to be 0<|kx|<1, which corresponds to the wavenumber range of the instability. If the residual is less than unity, the truncation in the eigenmode expansion may be considered as a representative of the full system, thus a candidate for reduced-order model building. On the contrary, the residual being around unity or more signifies the failure of the reduced representation in effectively capturing the overall nonlinear fluctuations.

The time evolution of the residuals is compared in Fig. 7 for varying forcing strengths. As expected, the unstable modes entirely capture the fluctuations in the linear phase (i.e., t30). In the nonlinear phase, however, the unstable modes capture only a rather limited fraction of the turbulent fluctuations. This is greatly improved when the stable modes are added. This suggests that the success of quasilinear models in capturing key properties of the turbulence can crucially depend on whether stable modes are considered when constructing such models.

FIG. 7.

Time traces of residuals, i.e., the fraction of energy missed in truncated bases, normalized to the total energy in the turbulent flow at each time step. The reconstruction uses truncated bases with unstable modes alone, and unstable and stable modes together. The forcing strength is varied in three different simulations with MA=10: (a) DKrook=1, (b) DKrook=0.1, and (c) DKrook=0. The unforced shear layer in (c) rapidly flattens out, and, thus, instability no longer drives the turbulence. As long as the turbulence is driven by the instability, the unstable and stable modes together can reconstruct a large fraction of the turbulent flow features in (a) and (b).

FIG. 7.

Time traces of residuals, i.e., the fraction of energy missed in truncated bases, normalized to the total energy in the turbulent flow at each time step. The reconstruction uses truncated bases with unstable modes alone, and unstable and stable modes together. The forcing strength is varied in three different simulations with MA=10: (a) DKrook=1, (b) DKrook=0.1, and (c) DKrook=0. The unforced shear layer in (c) rapidly flattens out, and, thus, instability no longer drives the turbulence. As long as the turbulence is driven by the instability, the unstable and stable modes together can reconstruct a large fraction of the turbulent flow features in (a) and (b).

Close modal

It is also interesting to note that the turbulence in the unforced shear layer, see Fig. 7(c), is different from the forced cases. In the former, the shear layer quickly flattens out and nearly shuts off the instability, leading to a decaying turbulence. Regardless of whether the unstable and/or stable modes are considered, the corresponding reconstructions fail to model the turbulence with any degree of accuracy. By contrast, when the shear layer is forced, a reduced representation of the turbulent flow with two modes (per wavenumber) is found to perform well, recovering a substantial fraction of the full nonlinear system.

A similar reconstruction is shown for various strengths of magnetic fields MA in Fig. 8, where residuals are time-averaged over a quasi-stationary state of turbulence. With stronger magnetic fields (lower MA), the vortices begin disrupting due to stronger Lorentz force and consequently generate more fluctuations at scales beyond the Kelvin–Helmholtz-instability (KHI) range.47 This accounts for an increase in the residual for low MA, although it remains below 0.2 for MA=10. For MA=3, the improvement with the inclusion of the stable modes is modest. Momentum transport by large-scale structures, formed from the unstable and stable modes, within the KHI range, however, may still dominate over the transport contributed by much smaller scales or the remaining continuum modes; hence, a quantitative analysis of transport will be conducted next.

FIG. 8.

Shown are the time-averaged residuals for different simulations with MA. Note the residuals are the fractions of energy missed in the truncated bases, compared to the total energy in the instantaneous full turbulent flow. Time-averaging is performed over a quasi-stationary state of turbulence (t=6001000). The reconstruction uses truncated bases with unstable modes alone, and unstable and stable modes together, leaving all the continuum modes. All simulations use DKrook=2. Note the dramatic improvement with the inclusion of the stable modes. For MA=3, the improvement is modest, as the stronger Lorentz force back-reacts on the large-scale turbulent flow, producing more fluctuations in the continuum modes.

FIG. 8.

Shown are the time-averaged residuals for different simulations with MA. Note the residuals are the fractions of energy missed in the truncated bases, compared to the total energy in the instantaneous full turbulent flow. Time-averaging is performed over a quasi-stationary state of turbulence (t=6001000). The reconstruction uses truncated bases with unstable modes alone, and unstable and stable modes together, leaving all the continuum modes. All simulations use DKrook=2. Note the dramatic improvement with the inclusion of the stable modes. For MA=3, the improvement is modest, as the stronger Lorentz force back-reacts on the large-scale turbulent flow, producing more fluctuations in the continuum modes.

Close modal

The Reynolds stress can be expressed in terms of the contribution from each wavenumber, which can further be decomposed into the contribution from each eigenmode. At a wavenumber kx, the Reynolds stress from all the fluctuations ϕ̂kx reads

τu(allmodes)=2Im[kxϕ̂kx·zϕ̂kx],
(13)

whereas the contributions from an unstable mode alone and from an unstable and a stable mode alone, at that wavenumber are, respectively, given as

τu(unstable)=2Im[kx(β1ϕ1,kx)·z(β1ϕ1,kx)]=2|β1|2Im[kxϕ1,kx·zϕ1,kx],
(14)

and

τu(stable)=2|β2|2Im[kxϕ2,kx·zϕ2,kx]=2|β2|2Im[kxϕ1,kx·zϕ1,kx]=2|β2|2Im[kxϕ1,kx·zϕ1,kx],
(15)

where ϕ̂kx is the Fourier transform of the streamfunction at wavenumber kx and ϕj,kx represents the z-dependent jth complex eigenmode: j =1 and 2 for unstable and stable modes, respectively. The conjugate symmetry of unstable and stable modes, as shown in Figs. 1(b) and 1(c), is used in Eq. (15), i.e., ϕ2,kx(z)=ϕ1,kx(z). The negative sign of the last expression in Eq. (15) corresponds to the up-gradient nature of momentum transport by stable modes, which was physically analyzed in Sec. III B and visually demonstrated in Fig. 2.

The summed contributions of unstable and stable modes in transport, however, can have cross-terms—quadratic correlations between unstable and stable modes—that do not appear in Eqs. (14) and (15) where the contributions from individual modes is shown. However, the cross-terms vanish when the unstable and stable modes are exactly complex conjugates of each other, as is the case for the ideal Kelvin–Helmholtz instability (when this conjugate symmetry is broken, e.g., in resistive tearing instability or in ion-temperature-gradient instability,60 the cross-terms can have non-zero contribution),

τu(unstable+stable)=2Im[kx(β1ϕ1,kx+β2ϕ2,kx)·z(β1ϕ1,kx+β2ϕ2,kx)]=2(|β1|2|β2|2)Im[kxϕ1,kx·zϕ1,kx]+cross-terms,
(16)

where

cross-terms=2Im[kx(β1ϕ1,kx)·z(β2ϕ2,kx)]+2Im[kx(β2ϕ2,kx)·z(β1ϕ1,kx)]=2Im[β1β2kxϕ1,kx·zϕ2,kx]+2Im[β2β1kxϕ2,kx·zϕ1,kx]=2Im[β1β2kxϕ1,kx·zϕ2,kx]+2Im[β2β1kxϕ1,kx·zϕ2,kx]=2Im[β1β2kxϕ1,kx·zϕ2,kx]+2Im[(β1β2kxϕ1,kx·zϕ2,kx)]=0.
(17)

Thus, we obtain

τu(unstable+stable)=2(|β1|2|β2|2)Im[kxϕ1,kx·zϕ1,kx]=τu(unstable)+τu(stable).
(18)

These relations inform us about the z-profile of the Reynolds stress, contributed by each wavenumber and each eigenmode. As largest momentum transport happens in the region with the largest flow-gradient, it is instructive to compute, in the forced shear layers, the turbulent stresses at the middle of the layer at z =0 and compare the stress contributions from different eigenmodes at various wavenumbers.

The total Reynolds stress from all modes and all wavenumbers in the simulations is compared in Fig. 9 with the stress contributions from the wavenumber range 0<|kx|<1, which is decomposed further into eigenmodes to assess the contribution of the unstable modes, stable modes, and their sum. The subplots demonstrate that the stable modes are highly efficient in transporting momentum in the up-gradient direction, as compared to the down-gradient transport by the unstable modes. Even for the strongest magnetic field MA=3, close to the instability threshold, the stable modes contribute significantly to a continuous reduction of the turbulent momentum flux. In addition, the occasional breakthroughs in stable-mode activity cause reversals of the transport direction. This reversal can be observed when the total Reynolds stress in the system is computed, without decomposing the stress into contributions by each eigenmode. However, when the stable modes are not overtaking the unstable modes in transport, the resulting down-gradient transport observed in simulations or experiments is difficult to interpret, in regard to the contributions of stable modes in subdominantly reducing the transport. An eigenmode decomposition of turbulent fluctuations, however, uncovers a complete picture, as is shown here.

FIG. 9.

Time variations of Reynolds stress at the middle of the shear layer, z =0. The stress contributions from unstable modes (blue), stable modes (orange), their sum (red), and full nonlinear fluctuations, i.e., all modes and all wavenumbers in the simulations (black), are compared, for varying strengths of magnetic fields: (a) MA=60, (b) MA=10, and (c) MA=3. Thin green lines represent the zero level. All simulations use DKrook=2. Although with stronger magnetic fields, the up-gradient momentum transport by stable modes is reduced, the up- and down-gradient transport nearly cancel each other throughout all cases.

FIG. 9.

Time variations of Reynolds stress at the middle of the shear layer, z =0. The stress contributions from unstable modes (blue), stable modes (orange), their sum (red), and full nonlinear fluctuations, i.e., all modes and all wavenumbers in the simulations (black), are compared, for varying strengths of magnetic fields: (a) MA=60, (b) MA=10, and (c) MA=3. Thin green lines represent the zero level. All simulations use DKrook=2. Although with stronger magnetic fields, the up-gradient momentum transport by stable modes is reduced, the up- and down-gradient transport nearly cancel each other throughout all cases.

Close modal

Similar variations of momentum transport across the middle of the shear layer are compared in Fig. 10 for different forcing strengths. Note the unforced case differs from the forced cases, as the nearly flattened shear layer has less momentum to be transported across the layer. As reported in Ref. 44, despite the profile relaxation, the two eigenmodes per wavenumber describe well the temporal variation of the Reynolds stress across the shear layer, although the stress itself is very low (note its vertical scale). In all cases, the stress captured via the sum of unstable and stable modes almost completely follows the total stress from all modes.

FIG. 10.

Time variations of Reynolds stress at the middle of the shear layer, z =0. The stress contributions from unstable modes (blue), stable modes (orange), their sum (red), and full nonlinear fluctuations, i.e., all modes and all wavenumbers in the simulations (black), are compared, for varying forcing strengths: (a) DKrook=25, (b) DKrook=1, and (c) DKrook=0. All simulations use MA=10. Thin green lines represent the zero level. Qualitative differences can be observed in unforced (DKrook=0) and forced cases (DKrook0): as instability extracts energy from the mean flow, the profile relaxation in the unforced layer leads to a decaying turbulence, and the transport rates become very small [note the vertical axis labels in (c)]. However, in all cases, the summed stable modes producing up-gradient transport nearly cancel the down-gradient transport by unstable modes. The addition of these two contributions produces a stress that is almost identical to the stress from all modes.

FIG. 10.

Time variations of Reynolds stress at the middle of the shear layer, z =0. The stress contributions from unstable modes (blue), stable modes (orange), their sum (red), and full nonlinear fluctuations, i.e., all modes and all wavenumbers in the simulations (black), are compared, for varying forcing strengths: (a) DKrook=25, (b) DKrook=1, and (c) DKrook=0. All simulations use MA=10. Thin green lines represent the zero level. Qualitative differences can be observed in unforced (DKrook=0) and forced cases (DKrook0): as instability extracts energy from the mean flow, the profile relaxation in the unforced layer leads to a decaying turbulence, and the transport rates become very small [note the vertical axis labels in (c)]. However, in all cases, the summed stable modes producing up-gradient transport nearly cancel the down-gradient transport by unstable modes. The addition of these two contributions produces a stress that is almost identical to the stress from all modes.

Close modal

In Fig. 11, the momentum transport by the unstable and stable modes is presented as a function of magnetic Prandtl number Pm=Rm/Re. All simulations until this point used Rm=500, which is now changed to Rm=50 and Rm=5000. In both cases of Pm=0.1 and Pm=10, the stable modes still substantially offset the turbulent momentum transport of the unstable modes. The shorter time trace for Pm=10 is due to the higher simulation cost. It should be noted that the quasi-stationary state in this simulation is still undergoing changes, unlike in the case of Pm=0.1 in Fig. 11(a) or Pm=1 in Fig. 9(b), all with the same MA=10,DKrook=2, and Re=500.

FIG. 11.

Time variations of Reynolds stress at the middle of the shear layer, z =0. The stress contributions from unstable modes (blue), stable modes (orange), their sum (red), and full nonlinear fluctuations, i.e., all modes and all wavenumbers in the simulations (black), are compared, for varying magnetic Prandtl numbers (resistivities): (a) Pm=0.1 and (b) Pm=10. All simulations use MA=10,DKrook=2, and Re=500. Thin green lines represent the zero level. It can be observed that the stable modes begin driving up-gradient momentum transport at around t30 when the nonlinear phase of evolution begins. By varying Pm by two orders of magnitude, around unity, the stable modes are found to substantially reduce the down-gradient transport; note the case of Pm=1 is shown in Fig. 9(b).

FIG. 11.

Time variations of Reynolds stress at the middle of the shear layer, z =0. The stress contributions from unstable modes (blue), stable modes (orange), their sum (red), and full nonlinear fluctuations, i.e., all modes and all wavenumbers in the simulations (black), are compared, for varying magnetic Prandtl numbers (resistivities): (a) Pm=0.1 and (b) Pm=10. All simulations use MA=10,DKrook=2, and Re=500. Thin green lines represent the zero level. It can be observed that the stable modes begin driving up-gradient momentum transport at around t30 when the nonlinear phase of evolution begins. By varying Pm by two orders of magnitude, around unity, the stable modes are found to substantially reduce the down-gradient transport; note the case of Pm=1 is shown in Fig. 9(b).

Close modal

The efficiency of time-averaged up-gradient momentum transport due to stable modes is compared in Fig. 12 with the time-averaged down-gradient transport due to unstable modes, via a measure, defined as follows:

Transportreductionefficiency=Up-gradienttransportbystablemodestDown-gradienttransportbyunstablemodest,
(19)

where At represents a time-averaging operation on A.

FIG. 12.

Parameter dependence of transport reduction efficiency, which is defined as the ratio of time-averaged up-gradient Reynolds stress due to stable modes and time-averaged down-gradient Reynolds stress due to unstable modes. The stress is measured at the middle of the shear layer, z =0, where the momentum transport is at its maximum. (a) Variations in MA=3,10,30,60,and120 with DKrook=2,Pm=1, and linear x-scale. (b) Variations in DKrook=0.1,1,2,and6 with MA=10,Pm=1, and logarithmic x-scale. (c) Variations in Pm=0.1,1,and10 (or, equivalent changes in resistivities) with MA=10,DKrook=2, and logarithmic x-scale. All plots have the same y-axis. The time-average for (a) and (b) is taken over a long quasi-stationary state of turbulence t=350900, while for (c), it is t=137237 where the quasi-stationary state is still undergoing changes. In all cases, substantial reduction of transport by stable modes is evident, which cancel, via their up-gradient transport, more than half of the down-gradient transport by unstable modes, and this fraction reaches up to 98%, see (a), for MA=60 and MA=120.

FIG. 12.

Parameter dependence of transport reduction efficiency, which is defined as the ratio of time-averaged up-gradient Reynolds stress due to stable modes and time-averaged down-gradient Reynolds stress due to unstable modes. The stress is measured at the middle of the shear layer, z =0, where the momentum transport is at its maximum. (a) Variations in MA=3,10,30,60,and120 with DKrook=2,Pm=1, and linear x-scale. (b) Variations in DKrook=0.1,1,2,and6 with MA=10,Pm=1, and logarithmic x-scale. (c) Variations in Pm=0.1,1,and10 (or, equivalent changes in resistivities) with MA=10,DKrook=2, and logarithmic x-scale. All plots have the same y-axis. The time-average for (a) and (b) is taken over a long quasi-stationary state of turbulence t=350900, while for (c), it is t=137237 where the quasi-stationary state is still undergoing changes. In all cases, substantial reduction of transport by stable modes is evident, which cancel, via their up-gradient transport, more than half of the down-gradient transport by unstable modes, and this fraction reaches up to 98%, see (a), for MA=60 and MA=120.

Close modal

Variations in magnetic field strength, forcing strength, and magnetic Prandtl number all demonstrate that the stable modes cancel an appreciable amount of the turbulent momentum flux associated with the unstable modes. On average, around 80% of the down-gradient flux is offset in this manner.

A remark should be made now regarding the use of unstable and stable modes for building a reliable reduced mode of transport for geo- and astro-physical problems. One approach would be to relate the activity of these two modes with a coefficient of diffusive flux (although the unstable and stable modes offer spatial profiles of transport as well, with both diffusive and non-diffusive fluxes, because they do not rely on an ad hoc eddy-viscosity model, which is an explicit diffusive-flux-based model). In the middle of the shear layer, the diffusive flux, however, dominates because of the maximum in the flow-gradient. The ad hoc turbulent viscosity can, thus, be defined,27 more importantly without a “free-parameter,” using Eq. (18) as

νturb=τu(unstable+stable)(dU0/dz)|z=0=0<kx<12(|β1|2|β2|2)Im[kxϕ1,kx·zϕ1,kx]|z=0.
(20)

Note that the denominator is unity for the shear-flow that has a linear profile in the vicinity of z =0. To assess the importance of stable modes in this construct, |β1|2|β2|2 can be written as |β1|2(1|β2|2/|β1|2). Since |β2|2 has been found to on the same order of |β1|2, e.g., see Fig. 12, where |β2|2/|β1|2 can range from 0.8 to 0.95, yielding (1|β2|2/|β1|2)0.05to0.2. Therefore, neglecting stable modes can overestimate the transport by a factor of 5 to 20.

With the above successful low-order representation of Reynolds stress, we now examine the fluctuations in the magnetic field that give rise of Maxwell stress. The stress can be quantified as

τb(allmodes)=2MA2Im[kxψ̂kx·zψ̂kx],
(21)
τb(unstable+stable)=2MA2(|β1|2β2|2)Im[kxψ1,kx·zψ1,kx],
(22)

where ψ̂kx is the Fourier transform of the flux function at a wavenumber kx and ψj,kx represents the z-dependent jth complex eigenmode (j =1 and 2 for unstable and stable modes, respectively). Again, cross-terms arising from the correlation between the unstable and stable modes can be shown to vanish, exactly as it was shown for the Reynolds stress in Eq. (17).

As can be seen in Fig. 13, the Reynolds stress is dominated by large scales, while the Maxwell stress involves a large number of different scales. Figure 13(a), using axes with linear scales, shows the dominance of Reynolds stress in the entire system, which the two-eigenmodes-per-wavenumber decomposition (unstable and stable modes) captures not only qualitatively but also quantitatively with great accuracy. In Fig. 13(b), a logarithmic scale is used to expose the range of small scales that contribute significantly to the magnetic fluctuations. Wavenumbers kx10 have major contributions, as opposed to kx<1 for the fluctuations of the flow. The fact that a large amount of flow energy resides at large scales suggests that the shear-flow turbulence may be amenable to some form of quasilinear modeling. Homogeneous isotropic turbulence, on the other hand, would not be reliably captured with such models, as no scale separation exists therein. Recent studies have highlighted that improved quasilinear models such as the generalized quasilinear approximation are realizable in systems with length or timescale separation.61 

FIG. 13.

Time-averaged turbulent stresses split into their contributions from different wavenumbers. (a) Stresses on a linear scale. (b) Log –log representation of the absolute value of the stresses. Note that only the wavenumbers kx<1 are Kelvin–Helmholtz-unstable. The simulation parameters are MA=60 and DKrook=2; the time-average is taken over a quasi-stationary state of turbulence, t=3501000. The total turbulent stress is dominated by the range |kx|<1, which is captured by the unstable and stable modes at those wavenumbers to a high precision. The small amount of stresses that are contributed by smaller scales of fluctuations spans a broad range of wavenumbers, due to the smaller scales in magnetic fields generated via straining by the flow.

FIG. 13.

Time-averaged turbulent stresses split into their contributions from different wavenumbers. (a) Stresses on a linear scale. (b) Log –log representation of the absolute value of the stresses. Note that only the wavenumbers kx<1 are Kelvin–Helmholtz-unstable. The simulation parameters are MA=60 and DKrook=2; the time-average is taken over a quasi-stationary state of turbulence, t=3501000. The total turbulent stress is dominated by the range |kx|<1, which is captured by the unstable and stable modes at those wavenumbers to a high precision. The small amount of stresses that are contributed by smaller scales of fluctuations spans a broad range of wavenumbers, due to the smaller scales in magnetic fields generated via straining by the flow.

Close modal

The magnetic fluctuations, on the other hand, span a broad range of scales. This can be physically interpreted as a result of the straining of the magnetic fields by the turbulent flow, which generates small scales in the magnetic fields.48,62–64 The straining process by the large-scale turbulent eddies converts the large-scale kinetic energy into the intermediate-scale magnetic energy.65 Magnetic fluctuations at such scales can then, via Lorentz force, feedback on the flow, although mostly at smaller scales. A comprehensive analysis of energy transfer for the present system will be reported in a forthcoming publication where nonlinear mode-coupling and energy transfer between fluctuations of discrete and continuum modes of velocity and magnetic fields are also analyzed.

To model any aspect of magnetic fluctuations, one must, thus, rely on tools such as statistical theories to obtain scaling laws that can offer insights into these fluctuations. One such approach is detailed next.

Until this point, the discrete modes—unstable and stable modes—which describe the turbulent flow well, have been our focus. The magnetic fluctuations, on the other hand, result from the straining of field lines by the flow, exciting the remaining continuum modes. Hence, these modes are necessary for a successful reconstruction of the magnetic fluctuations. Thus, we seek a simple scaling law for the saturated turbulent amplitudes of the continuum modes.

1. Analytical prediction for continuum mode amplitudes

It is instructive to write the nonlinear MHD equations in the eigenmode basis, arriving at what is also referred to as the mode-amplitude evolution equation,30,43,45,66,67

tβj(kx)=iωj(kx)βj(kx)+kx,kx,m,nkx+kx=kxCjmn(kx,kx)βm(kx)βn(kx),
(23)

where βj(kx) represents the complex amplitude of an eigenmode j at wavenumber kx with ωj the associated mode-frequency; the nonlinear mode coupling coefficient Cjmn measures the three-wave overlap, which dictates the strength of nonlinear beating between an eigenmode m at wavenumber kx and an eigenmode n at wavenumber kx, driving an eigenmode j at wavenumber kx (with the constraint kx+kx=kx).

For the continuum modes, as was mentioned in Sec. III A, their frequencies depend linearly on the wavenumber kx as59ω/kx+Uref(z)±vA,ref(z)=0. This implies ωj=ωkx.

Heuristically, the scaling of the nonlinear mode coupling coefficient with wavenumber can be obtained in the following manner: in Eqs. (2a) and (2b), the separation of linear and nonlinear terms arises in Poisson brackets. Consider a prototype equation,

tP̃={P,Q}+={P̃,Q̃}+{P̃,Q0}+{P0,Q̃}+,
(24)

where P=P0+P̃ and Q=Q0+Q̃ represent two fields (e.g., 2ϕ or ψ for the present problem), with P0 representing the x-averaged mean component of P and P̃ standing for perturbations. The linear term, e.g., ikxP̂·zQ0, which is in spectral space, contains only one perturbed field, whereas the nonlinear term, e.g., ikxP̂·zQ̂, has two perturbed fields, with P̂ and Q̂ representing the Fourier-transformed quantities at wavenumbers kx and kx, respectively. It may be supposed that the derivative z on the perturbed quantities is roughly on the scale of |z|kx. (This can be shown analytically for all the eigenmodes, where the background flow is approximately uniform, see Ref. 43.) Notice, however, that this argument applies only to the perturbations: the operator z acting on Q0 clearly does not produce a factor of kx, which is zero for the mean component Q0. We now use this distinction to make a prediction for the amplitudes of perturbations, in particular the continuum mode-amplitudes. The linear and nonlinear terms, thus, assume the forms ikxP̂·zQ0 and ikxP̂·ikxQ̂, respectively.

In Eq. (24), expanding the perturbations in the eigenmode basis, e.g., P̂=lβlP̂l with P̂l representing the lth eigenmode, and diagonalizing the linear terms (operator), one finds an equation of the form given in Eq. (23). We can now attempt to understand the behavior of the nonlinear mode coupling coefficients that drive the continuum modes. Assuming nonlinear interactions between the continuum modes are local in spectral space—interaction of three wavenumbers of similar scales—the nonlinear term in Eq. (24) simplifies, e.g., ikxP̂·ikxQ̂ becomes kx2P̂Q̂; note the linear term has the form ikxP̂·zQ0.

In assuming local interaction between the continuum modes spectral space in kx, the involvement of unstable and stable modes in nonlinear interactions is ignored, which, otherwise, could bring in non-local effects. This may be a valid assumption for continuum modes at scales much above the Kelvin–Helmholtz-unstable wavenumber range, i.e., kx>1, as the wavenumber convolution constraint of kx=kx+kx does not allow two (un-)stable modes to beat together to drive a continuum mode at large kx, e.g., kx>2.

Continuing with the above assumption, the nonlinear term kx2P̂Q̂ has one extra kx compared to the linear term ikxP̂·zQ0. This implies that, for the continuum modes, the nonlinear mode coupling coefficients C are expected to scale as

Ckx2,
(25)

because the linear term for the continuum modes in Eq. (23) has the eigenfrequency that depends linearly on kx, i.e.,

ωkx.
(26)

Such a property of nonlinear coupling coefficient is common in other turbulence calculations, as well.26 

In order to obtain a phenomenological scaling law, we now make no distinction between different continuum modes and, thus, balance the linear and nonlinear terms of Eq. (23) in the quasi-stationary state as ωβCβ2. Inserting their asymptotic dependences on kx, the amplitudes of continuum modes are found to follow:

βkx1.
(27)

Note that the assumptions made in arriving at this simple scaling law are crude. The next step will be to determine from nonlinear simulations whether this scaling can be recovered or whether a number of assumptions made above render the result inapplicable.

2. Numerical verification

Time-averaged eigenmode amplitudes from nonlinear simulations, after multiplying with kx, are plotted in Fig. 14(a) as functions of kx and eigenmode index j, arranged in order of increasing real frequency of the eigenmodes. The appearance of vertical near-equicontour lines signifies that eigenmodes are excited in a similar pattern across a large range of scales.

FIG. 14.

(a) Dependence of time-averaged eigenmode amplitudes on kx. The indices js of the eigenmodes are arranged in increasing order of their real frequencies. Vertical lines signify the self-similar cascade of energy to small scales in kx. Mode amplitudes are averaged over a quasi-stationary state of turbulence t=400800 for a simulation with MA=120 and DKrook=2. (b) kx spectra of mode amplitudes in a nonlinear simulation (shown with empty circles) in a log10log10 plot. Shown with a solid line is the analytical prediction, made for the wavenumbers that lie beyond the Kelvin–Helmholtz-unstable range, i.e., for kx>1. The spectral index, predicted based upon a number of simple assumptions, can be seen to fit the data reasonably well.

FIG. 14.

(a) Dependence of time-averaged eigenmode amplitudes on kx. The indices js of the eigenmodes are arranged in increasing order of their real frequencies. Vertical lines signify the self-similar cascade of energy to small scales in kx. Mode amplitudes are averaged over a quasi-stationary state of turbulence t=400800 for a simulation with MA=120 and DKrook=2. (b) kx spectra of mode amplitudes in a nonlinear simulation (shown with empty circles) in a log10log10 plot. Shown with a solid line is the analytical prediction, made for the wavenumbers that lie beyond the Kelvin–Helmholtz-unstable range, i.e., for kx>1. The spectral index, predicted based upon a number of simple assumptions, can be seen to fit the data reasonably well.

Close modal

The eigenmodes that lie within the yellow bands are localized in space (z-axis), but the band spans a range of heights, outside the shear layer |z|<1. Empirically, we note that the center of the rightmost (leftmost) band corresponds to ω/kx=U00 (ω/kx=U00), where U00=1. These thick bands represent all eigenmodes that have phase speeds ω/kx=U00+cvA,0 (ω/kx=U00+cvA,0), where 1<c<1 and vA,0=1/MA; note that |c|=1 is not included in these bands. All of these eigenmodes have peaks and oscillations in their eigenfunctions outside of the shear layer. In the layer, the unstable and stable modes maintain their dominance, and, thus, these two discrete modes alone almost completely regulate the momentum transport across the layer, as was noted in Sec. V D.

It is of interest now to compute from numerical simulation data how the amplitude of each eigenmode j falls off with kx and construct a j-averaged spectral index. To this end, we note the amplitude βj(kx=0.2) for each mode j at kx=0.2 (the first wavenumber in the simulation) and compute a scaled mode-amplitude β̂j(kx) as

β̂j(kx)=βj(kx)βj(kx=0.2),
(28)

which is expected to falloff with kx as kxα. In principle, the spectral index α can depend on the eigenmode index j, but a j-averaged spectral index is sought now, following the procedure:

β̂j(kx)(kx/0.2)α,logβ̂j(kx)αlogkx,logβ̂j(kx)jαjlogkx.
(29)

This j-averaged spectral distribution of the mode-amplitudes informs how, on average, each eigenmode amplitude depends on kx.

A plot of logβ̂j(kx)j vs logkx is shown in Fig. 14(b), along with the analytical prediction of inverse-in-wavenumber falloff of the mode amplitudes, at scales above the Kelvin–Helmholtz-unstable wavenumber range. It should be highlighted that the computation of all the eigenmode amplitudes at each wavenumber at each simulation time is computationally demanding, as the process requires the computation of modified left eigenmodes for each right eigenmode at each wavenumber, apart from the mode projection calculation at each time step. Therefore, only the first 24 Fourier modes in kx are shown in Fig. 14.

A finding in Fig. 14 is the identification of self-similar cascade of energy to smaller scales (larger kx) in eigenmode space. This result also hints that the interaction involving the continuum modes may be reasonably simplified and potentially be valuable in estimating the amplitudes of unstable and stable modes, using closure theories (see Ref. 26 for a recent example).

We have investigated MHD turbulence in two-dimensions, driven by a forced unstable shear flow, using a complete eigenmode decomposition of fluctuations in nonlinear simulations, which exposes the nonlinearly saturated excitation level of each eigenmode and its role. Intrinsic to linear instability, the unstable modes derive fluctuation energy from the mean flow gradient. The linearly decaying stable modes, however, contain almost the same amount of energy as the unstable mode, which they receive via nonlinear excitation. This truncated basis of two eigenmodes per wavenumber is found to reconstruct essential large-scale features of turbulent flow and the associated momentum transport via Reynolds stress. Quantifying transport due to unstable modes alone shows an overestimation up to an order of magnitude higher relative to the actual flux. The reduction in the flux is identified to be due to the continuous up-gradient transport by the stable modes, which causes a near-cancellation of down-gradient transport driven by unstable modes.

The continuum modes, on the other hand, describe small-scale fluctuations of the flow and magnetic field, where the above large-scale unstable and stable modes manifest themselves as a quasi-coherent vortex. To predict the mode amplitudes of the continuum, a simple scaling law is derived from the governing nonlinear equations, and the predicted inverse-in-wavenumber falloff rate is found reasonably agree with the simulation data.

Although both the momentum transport and fluctuation energy are largely described by the discrete modes, the former is more efficiently captured (Figs. 9–13) as almost all the momentum transport occurs near z =0, which is the region where the discrete modes dominate [Figs. 1(b) and 1(c)]. The fluctuation energy, on the other hand, is related to fluctuations that are scattered in and around the large-scale eddies; a portion of this energy is claimed by the continuum modes, although a large fraction still belongs to the discrete modes [Figs. 5(b) and 6–8].

Transport reduction by stable modes can also be used to improve phenomenological constructs like eddy viscosity, which are generally agnostic as to the nonlinear excitation of stable modes. By predicting the turbulent amplitudes of the unstable and stable modes, for astrophysically relevant parameters, e.g., very large Rm, Re, and Pm compared to unity, the simple relation between turbulent viscosity and eigenmode amplitudes [in Eq. (20)] can be exploited to reliably model transport processes in astrophysical objects, which, otherwise, cannot be solved using current state-of-the-art direct numerical simulations. It should be noted that such a prediction for the mode amplitudes |β1| and |β2| was recently made for ion-temperature-gradient-driven fusion plasma turbulence17 using statistical closure theory.68 Undertaking such a task for the present system is interesting, but beyond the purview of this work and will, thus, be left for future investigations.

The reduced representation of turbulent flow and transport presented here is also useful for building sub-grid-scale models, which can allow performing nonlinear simulations at extreme parameters with less-intensive computational demands. Progress can, thus, be made in seeking models that reduce the number of degrees of freedom while capturing essential features of the turbulent system. Techniques like proper orthogonal decomposition, dynamic mode decomposition, etc. also exist for such purposes,69 but they operate on output from nonlinear simulations, and it can be difficult to assign intuitive physical meaning to the characteristic mode structures. Here, the truncated eigenmode basis, composed of the unstable and stable modes, has been demonstrated to reconstruct nonlinear fluctuations to an appreciable degree, thus suggesting that they can be leveraged as a physically motivated basis for extreme parameter studies, without having to first perform a direct numerical simulation. These modes may also be useful for generating, via their nonlinear interactions with continuum modes, fluctuations associated with the continuum modes. Such a test could be performed to analyze magnetic fluctuations. The reduced basis, composed of unstable and stable mode alone, can also serve in direct statistical simulations,54 which have shown promises toward simulating the slowly evolving turbulent statistics, e.g., two-point two-time correlations, three-point correlations between fluctuating fields, etc., rather than the fast-evolving field variables themselves, e.g., flow velocities. Other improved forms of quasilinear models like the generalized quasilinear approximation61 may also benefit from using this truncated basis. This possibility will be explored in a separate publication.

In the future, procedures similar to that employed here can be used to examine the properties of other forms of instability-driven turbulence, such as magneto-rotational-instability-driven6 and stratified-shear-flow-driven turbulence.70 Building nonlinear energy transfer diagnostics in shear-flow turbulence to study the physical processes and scales that impact the difference in unstable- and stable-mode amplitudes are another possible avenue. Such investigations constitute steps toward deployment in service to one-dimensional stellar transport models.5 Central to improved predictiveness are the stable modes, whose properties will similarly require additional studies.

The authors appreciate K. Burns, E. H. Anders, and the Dedalus developers for assistance in several aspects of leveraging the numerical code. Thanks also to J. Fuller and other participants of the program “Transport in Stellar Interiors, 2021” at the Kavli Institute of Theoretical Physics for useful discussions. This material is based upon work funded by the Department of Energy (No. DE-SC0022257) through the NSF/DOE Partnership in Basic Plasma Science and Engineering. We also gratefully acknowledge support from NSF Grant Nos. AST-1814327 and AST-1908338. The simulations reported herein were performed using the XSEDE supercomputing resources via Allocation No. TG-PHY130027.

The authors have no conflicts to disclose.

Bindesh Tripathi: Writing – review and editing (equal). Adrian Everett Fraser: Writing – review and editing (equal). Paul Terry: Writing – review and editing (equal). Ellen G. Zweibel: Writing – review and editing (equal). MJ Pueschel: Writing – review and editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Due to the non-normality of the linear operator of the shear-flow instability, the eigenmodes are not orthogonal. This presents a significant challenge in the computation of mode amplitudes. An additional challenge is the generalized eigenvalue nature of the problem at hand, when written in vorticity formalism, as in Eq. (6a). This differs from the standard eigenvalue problem, Lξj=λjξj, where L is a linear operator whose j-th eigenmode is ξj with eigenvalue λj. The generalized eigenvalue problem that we encounter here is

LXj=ωjMXj,
(A1)

where L is a linear operator, M is another linear operator containing the Laplacian operation for our problem, and Xj is the jth (right) eigenmode with corresponding (right) eigenvalue ωj. (Often times, the distinction between left and right eigenmodes of a linear operator is not made as they happen to be the same; however, this is not the case here for the non-normal operator.) The right eigenmodes, although non-orthogonal to each other, can be made orthogonal with an appropriate weight factor to the left eigenmodes, which are solutions to another eigenvalue problem,66,67

YjTL=σjYjTM.
(A2)

Here, YjT is the transpose of the left eigenmode with its left eigenvalue σj. A slight reformulation is possible to this equation by taking Hermitian-transpose,

LYj=σjMYj.
(A3)

In the eigenvalue solver in Dedalus, the matrices L and M for each wavenumber are passed, and their eigenmodes Yj and eigenvalues σj are found. It can be shown that the eigenvalues σj and ωj are the same (i.e., σj=ωj), by analyzing Eqs. (A1) and (A2). A modified orthogonality relation between the left and right eigenmodes can now be derived,66,67

(YjTL)Xi=YjT(LXi)=YjT(ωiMXi)(σjYjTM)Xi=YjT(ωiMXi)YjTMXi(σjωi)=0YjTMXiδi,j,
(A4)

which means the left and right eigenmodes are orthogonal to each other with a weight factor M, as long as their eigenvalues differ (σjωi). For numerical computation, it is convenient to define YjTM=(MTYj)T=ZjT, where Zj is the modified left eigenmode, which is—by construct—orthogonal to the right eigenmode without any weight factor: ZjTXiδi,j. Using this relation, the eigenmode coefficients βj(kx,t) in the eigenmode expansion of turbulent fluctuations are computed at each wavenumber and at each time.

Same as in Fig. 5, but for MA=120 (Fig. 15).

FIG. 15.

Same as in Fig. 5, but for MA=120. The amplitude of the nonlinearly excited stable mode is almost exactly the same as that of the unstable mode. The nature of their oscillations is also similar, although the oscillations in the stable-mode-amplitude lags behind that of the unstable mode. The lag is likely an outcome of a time-delay in energy transfer from the unstable to the stable mode at the same wavenumber, which, thus, requires a series of nonlinear interactions with fluctuations at other wavenumbers.

FIG. 15.

Same as in Fig. 5, but for MA=120. The amplitude of the nonlinearly excited stable mode is almost exactly the same as that of the unstable mode. The nature of their oscillations is also similar, although the oscillations in the stable-mode-amplitude lags behind that of the unstable mode. The lag is likely an outcome of a time-delay in energy transfer from the unstable to the stable mode at the same wavenumber, which, thus, requires a series of nonlinear interactions with fluctuations at other wavenumbers.

Close modal
1.
E. C.
Harding
,
J. F.
Hansen
,
O. A.
Hurricane
,
R. P.
Drake
,
H. F.
Robey
,
C. C.
Kuranz
,
B. A.
Remington
,
M. J.
Bono
,
M. J.
Grosskopf
, and
R. S.
Gillespie
, “
Observation of a Kelvin-Helmholtz instability in a high-energy-density plasma on the omega laser
,”
Phys. Rev. Lett.
103
,
045005
(
2009
).
2.
H.
Hasegawa
,
M.
Fujimoto
,
T.-D.
Phan
,
H.
Rème
,
A.
Balogh
,
M. W.
Dunlop
,
C.
Hashimoto
, and
R.
TanDokoro
, “
Transport of solar wind into Earth's magnetosphere through rolled-up Kelvin-Helmholtz vortices
,”
Nature
430
,
755
(
2004
).
3.
D. W.
Waugh
,
A. H.
Sobel
, and
L. M.
Polvani
, “
What is the polar vortex and how does it influence weather?
,”
Bull. Am. Meteorol. Soc.
98
,
37
(
2017
).
4.
P. L.
Read
,
R. M. B.
Young
, and
D.
Kennedy
, “
The turbulent dynamics of Jupiter's and Saturn's weather layers: Order out of chaos?
,”
Geosci. Lett.
7
,
10
(
2020
).
5.
J.
Fuller
,
A. L.
Piro
, and
A. S.
Jermyn
, “
Slowing the spins of stellar cores
,”
Mon. Not. R. Astron. Soc.
485
,
3661
(
2019
).
6.
M. E.
Pessah
,
C.-K.
Chan
, and
D.
Psaltis
, “
The signature of the magnetorotational instability in the Reynolds and Maxwell stress tensors in accretion discs
,”
Mon. Not. R. Astron. Soc.
372
,
183
(
2006
).
7.
J.
Goodman
and
G.
Xu
, “
Parasitic instabilities in magnetized, differentially rotating disks
,”
Astrophys. J.
432
,
213
(
1994
).
8.
J.
Alves
,
C.
Zucker
,
A. A.
Goodman
,
J. S.
Speagle
,
S.
Meingast
,
T.
Robitaille
,
D. P.
Finkbeiner
,
E. F.
Schlafly
, and
G. M.
Green
, “
A Galactic-scale gas wave in the solar neighbourhood
,”
Nature
578
,
237
(
2020
).
9.
R.
Fleck
, “
The ‘Radcliffe Wave’ as a Kelvin–Helmholtz instability
,”
Nature
583
,
E24
(
2020
).
10.
A.
Miura
, “
Self-organization in the two-dimensional Kelvin-Helmholtz instability
,”
Phys. Rev. Lett.
83
,
1586
(
1999
).
11.
D.
Lecoanet
,
M.
McCourt
,
E.
Quataert
,
K. J.
Burns
,
G. M.
Vasil
,
J. S.
Oishi
,
B. P.
Brown
,
J. M.
Stone
, and
R. M.
O'Leary
, “
A validated non-linear Kelvin-Helmholtz benchmark for numerical hydrodynamics
,”
Mon. Not. R. Astron. Soc.
455
,
4274
(
2016
).
12.
C.
Ho
and
P.
Huerre
, “
Perturbed free shear layers
,”
Annu. Rev. Fluid Mech.
16
,
365
(
1984
).
13.
F. K.
Browand
and
C. M.
Ho
, “
The mixing layer, an example of quasi two-dimensional turbulence
,”
J. Mec. Theor. Appl.
2
,
99
(
1983
); available at https://ui.adsabs.harvard.edu/abs/1983JMTAS.......99B/abstract.
14.
V. P.
Starr
and
N. E.
Gaut
, “
Negative viscosity
,”
Sci. Am.
223
,
72
(
1970
).
15.
A.
Miura
and
T.
Sato
, “
Theory of vortex nutation and amplitude oscillation in an inviscid shear instability
,”
J. Fluid Mech.
86
,
33
(
1978
).
16.
W.
Horton
,
T.
Tajima
, and
T.
Kamimura
, “
Kelvin–Helmholtz instability and vortices in magnetized plasma
,”
Phys. Fluids
30
,
3485
(
1987
).
17.
P. W.
Terry
,
P.-Y.
Li
,
M. J.
Pueschel
, and
G. G.
Whelan
, “
Threshold heat-flux reduction by near-resonant energy transfer
,”
Phys. Rev. Lett.
126
,
025004
(
2021
).
18.
G. G.
Whelan
,
M. J.
Pueschel
, and
P. W.
Terry
, “
Nonlinear electromagnetic stabilization of plasma microturbulence
,”
Phys. Rev. Lett.
120
,
175002
(
2018
).
19.
M. J.
Pueschel
,
B. J.
Faber
,
J.
Citrin
,
C. C.
Hegna
,
P. W.
Terry
, and
D. R.
Hatch
, “
Stellarator turbulence: Subdominant eigenmodes and quasilinear modeling
,”
Phys. Rev. Lett.
116
,
085001
(
2016
).
20.
K. D.
Makwana
,
P. W.
Terry
,
M. J.
Pueschel
, and
D. R.
Hatch
, “
Subdominant modes in zonal-flow-regulated turbulence
,”
Phys. Rev. Lett.
112
,
095002
(
2014
).
21.
D. R.
Hatch
,
F.
Jenko
,
A. B.
Navarro
, and
V.
Bratanov
, “
Transition between saturation regimes of gyrokinetic turbulence
,”
Phys. Rev. Lett.
111
,
175001
(
2013
).
22.
D. R.
Hatch
,
P. W.
Terry
,
F.
Jenko
,
F.
Merz
, and
W. M.
Nevins
, “
Saturation of gyrokinetic turbulence through damped eigenmodes
,”
Phys. Rev. Lett.
106
,
115003
(
2011
).
23.
M. J.
Pueschel
,
P.-Y.
Li
, and
P. W.
Terry
, “
Predicting the critical gradient of ITG turbulence in fusion plasmas
,”
Nucl. Fusion
61
,
054003
(
2021
).
24.
P.-Y.
Li
and
P. W.
Terry
, “
Assessing physics of ion temperature gradient turbulence via hierarchical reduced-model representations
,”
Phys. Plasmas
29
,
042301
(
2022
).
25.
G. G.
Whelan
,
M. J.
Pueschel
,
P. W.
Terry
,
J.
Citrin
,
I. J.
McKinney
,
W.
Guttenfelder
, and
H.
Doerk
, “
Saturation and nonlinear electromagnetic stabilization of ITG turbulence
,”
Phys. Plasmas
26
,
082302
(
2019
).
26.
P. W.
Terry
,
B. J.
Faber
,
C. C.
Hegna
,
V. V.
Mirnov
,
M. J.
Pueschel
, and
G. G.
Whelan
, “
Saturation scalings of toroidal ion temperature gradient turbulence
,”
Phys. Plasmas
25
,
012308
(
2018
).
27.
A. E.
Fraser
,
M. J.
Pueschel
,
P. W.
Terry
, and
E. G.
Zweibel
, “
Role of stable modes in driven shear-flow turbulence
,”
Phys. Plasmas
25
,
122303
(
2018
).
28.
K. D.
Makwana
,
P. W.
Terry
,
J.-H.
Kim
, and
D. R.
Hatch
, “
Damped eigenmode saturation in plasma fluid turbulence
,”
Phys. Plasmas
18
,
012302
(
2011
).
29.
K. D.
Makwana
,
P. W.
Terry
, and
J.-H.
Kim
, “
Role of stable modes in zonal flow regulated turbulence
,”
Phys. Plasmas
19
,
062310
(
2012
).
30.
P. W.
Terry
,
D. A.
Baver
, and
S.
Gupta
, “
Role of stable eigenmodes in saturated local plasma turbulence
,”
Phys. Plasmas
13
,
022307
(
2006
).
31.
R. H.
Levy
and
R. W.
Hockney
, “
Computer experiments on low-density crossed-field electron beams
,”
Phys. Fluids
11
,
766
(
1968
).
32.
N. J.
Zabusky
and
G. S.
Deem
, “
Dynamical evolution of two-dimensional unstable shear flows
,”
J. Fluid Mech.
47
,
353
(
1971
).
33.
L.-S.
Huang
and
C.-M.
Ho
, “
Small-scale transition in a plane mixing layer
,”
J. Fluid Mech.
210
,
475
(
1990
).
34.
R. D.
Moser
and
M. M.
Rogers
, “
The three-dimensional evolution of a plane mixing layer: Pairing and transition to turbulence
,”
J. Fluid Mech.
247
,
275
(
1993
).
35.
J. J.
Riley
and
R. W.
Metcalfe
, “
Direct numerical simulation of a perturbed turbulent mixing layer
,” AIAA Paper No. 1980-0274,
1980
.
36.
D.
Oster
and
I.
Wygnanski
, “
The forced mixing layer between parallel streams
,”
J. Fluid Mech.
123
,
91
(
1982
).
37.
Y.
Ito
,
K.
Nagata
,
Y.
Sakai
, and
O.
Terashima
, “
Momentum and mass transfer in developing liquid shear mixing layers
,”
Exp. Therm. Fluid Sci.
51
,
28
(
2013
).
38.
A.
López Zazueta
and
L.
Zavala Sansón
, “
Self-oscillations of a two-dimensional shear flow with forcing and dissipation
,”
Phys. Fluids
30
,
044101
(
2018
).
39.
A.
VanDine
,
H. T.
Pham
, and
S.
Sarkar
, “
Turbulent shear layers in a uniformly stratified background: DNS at high Reynolds number
,”
J. Fluid Mech.
916
,
A42
(
2021
).
40.
A. K. M. F.
Hussain
and
K. B. M. Q.
Zaman
, “
An experimental study of organized motions in the turbulent plane mixing layer
,”
J. Fluid Mech.
159
,
85
(
1985
).
41.
A. K. M. F.
Hussain
, “
Coherent structures and turbulence
,”
J. Fluid Mech.
173
,
303
(
1986
).
42.
L.
Landau
, “
On the problem of turbulence
,”
C. R. Acad. Sci. U. R. S. S.
44
,
311
(
1944
).
43.
A. E.
Fraser
,
P. W.
Terry
,
E. G.
Zweibel
, and
M. J.
Pueschel
, “
Coupling of damped and growing modes in unstable shear flow
,”
Phys. Plasmas
24
,
062304
(
2017
).
44.
A. E.
Fraser
,
P. W.
Terry
,
E. G.
Zweibel
,
M. J.
Pueschel
, and
J. M.
Schroeder
, “
The impact of magnetic fields on momentum transport and saturation of shear-flow instability by stable modes
,”
Phys. Plasmas
28
,
022309
(
2021
).
45.
B.
Tripathi
,
A. E.
Fraser
,
P. W.
Terry
,
E. G.
Zweibel
, and
M. J.
Pueschel
, “
Mechanism for sequestering magnetic energy at large scales in shear-flow turbulence
,” arXiv:2205.01298 (
2022
).
46.
S.
Chandrasekhar
,
Hydrodynamic and Hydromagnetic Stability
(
Clarendon Press
,
Oxford
,
1961
).
47.
J.
Mak
,
S. D.
Griffiths
, and
D. W.
Hughes
, “
Vortex disruption by magnetohydrodynamic feedback
,”
Phys. Rev. Fluids
2
,
113701
(
2017
).
48.
A. A.
Schekochihin
,
J. L.
Maron
,
S. C.
Cowley
, and
J. C.
McWilliams
, “
The small-scale structure of magnetohydrodynamic turbulence with large magnetic Prandtl numbers
,”
Astrophys. J.
576
,
806
(
2002
).
49.
D.
Biskamp
,
Magnetohydrodynamic Turbulence
(
Cambridge University Press
,
Cambridge
,
2003
).
50.
F.
Ebrahimi
,
S. C.
Prager
, and
D. D.
Schnack
, “
Saturation of magnetorotational instability through magnetic field generation
,”
Astrophys. J.
698
,
233
(
2009
).
51.
M. J.
Pueschel
,
D.
Told
,
P. W.
Terry
,
F.
Jenko
,
E. G.
Zweibel
,
V.
Zhdankin
, and
H.
Lesch
, “
Magnetic reconnection turbulence in strong guide fields: Basic properties and application to coronal heating
,”
Astrophys. J., Suppl. Ser.
213
,
30
(
2014
).
52.
J. B.
Marston
,
E.
Conover
, and
T.
Schneider
, “
Statistics of an unstable barotropic jet from a cumulant expansion
,”
J. Atmos. Sci.
65
,
1955
(
2008
).
53.
K. M.
Smith
,
C. P.
Caulfield
, and
J. R.
Taylor
, “
Turbulence in forced stratified shear flows
,”
J. Fluid Mech.
910
,
A42
(
2021
).
54.
A.
Allawala
,
S. M.
Tobias
, and
J. B.
Marston
, “
Dimensional reduction of direct statistical simulation
,”
J. Fluid Mech.
898
,
A21
(
2020
).
55.
K. J.
Burns
,
G. M.
Vasil
,
J. S.
Oishi
,
D.
Lecoanet
, and
B. P.
Brown
, “
Dedalus: A flexible framework for numerical simulations with spectral methods
,”
Phys. Rev. Res.
2
,
023068
(
2020
).
56.
C. M.
Bender
,
PT Symmetry: In Quantum and Classical Physics
(
World Scientific Publishing
,
2019
).
57.
D. R.
Hatch
,
F.
Jenko
,
A. B.
Navarro
,
V.
Bratanov
,
P. W.
Terry
, and
M. J.
Pueschel
, “
Linear signatures in nonlinear gyrokinetics: Interpreting turbulence with pseudospectra
,”
New J. Phys.
18
,
075018
(
2016
).
58.
Y.
Fu
and
H.
Qin
, “
The physics of spontaneous parity-time symmetry breaking in the Kelvin-Helmholtz instability
,”
New J. Phys.
22
,
083040
(
2020
).
59.
K. M.
Case
, “
Stability of inviscid plane Couette flow
,”
Phys. Fluids
3
,
143
(
1960
).
60.
P. W.
Terry
,
D. A.
Baver
, and
D. R.
Hatch
, “
Reduction of inward momentum flux by damped eigenmodes
,”
Phys. Plasmas
16
,
122305
(
2009
).
61.
J. B.
Marston
,
G. P.
Chini
, and
S. M.
Tobias
, “
Generalized quasilinear approximation: Application to zonal jets
,”
Phys. Rev. Lett.
116
,
214501
(
2016
).
62.
G. K.
Batchelor
, “
On the spontaneous magnetic field in a conducting liquid in turbulent motion
,”
Proc. Roy. Soc. London, Ser. A
201
,
405
(
1950
).
63.
G. K.
Batchelor
and
I.
Proudman
, “
The effects of rapid distortion of a fluid in turbulent motion
,”
Q. J. Mech. Appl. Math
7
,
83
(
1954
).
64.
A. A.
Townsend
,
The Structure of Turbulent Shear Flow
, 2nd ed. (
Cambridge University Press
,
Cambridge
,
1976
).
65.
A.
Alexakis
,
P. D.
Mininni
, and
A.
Pouquet
,
Phys. Rev. E
72
,
046301
(
2005
).
66.
K. J.
Burns
, “
Flexible spectral algorithms for simulating astrophysical and geophysical flows
,” Ph.D. thesis (
Massachusetts Institute of Technology
,
2018
).
67.
A. E.
Fraser
, “
Role of stable eigenmodes in shear-flow instability saturation and turbulence
,” Ph.D. thesis (
University of Wisconsin-Madison
,
2020
).
68.
S. A.
Orszag
, “
Analytical theories of turbulence
,”
J. Fluid Mech.
41
,
363
(
1970
).
69.
K.
Taira
,
S. L.
Brunton
,
S. T. M.
Dawson
,
C. W.
Rowley
,
T.
Colonius
,
B. J.
McKeon
,
O. T.
Schmidt
,
S.
Gordeyev
,
V.
Theofilis
, and
L. S.
Ukeiley
, “
Modal analysis of fluid flows: An overview
,”
AIAA J.
55
,
4013
(
2017
).
70.
P.
Garaud
, “
Double-diffusive convection at low Prandtl number
,”
Annu. Rev. Fluid Mech.
50
,
275
(
2018
).