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.

## I. INTRODUCTION

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 fluids^{32} 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 assumed^{15} 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 events^{12,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 interest^{42} 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-root^{46} 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 instability^{46} 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.

## II. MODEL AND SIMULATION SETUP

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

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. Background flow, magnetic field, and forcing

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\u0302$ and a flow-aligned magnetic field, initially uniform, as $B(x,z,t=0)=B0x\u0302$. The parameters $a,U0,and\u2009B0$ 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 *U*_{0} to the Alfvèn speed can be written as the Alfvènic Mach number $MA=U04\pi \rho /B0$. The viscosity and resistivity are quantified via fluid Reynolds number $Re=U0a/\nu $ and magnetic Reynolds number $Rm=U0a/\eta $, respectively.

In two dimensions, a more convenient formalism is available, using the streamfunction $\varphi $ and flux function *ψ*. Defining $u=y\u0302\xd7\u2207\varphi $ and $B=y\u0302\xd7\u2207\psi $, the vorticity and the current become $\u22072\varphi y\u0302$ and $\u22072\psi y\u0302$, respectively. Taking the curl of Eq. (1b) and rewriting Eq. (1d) in terms of the stream- and flux-functions yield^{49}

where the Poisson bracket is ${P,Q}=\u2202xP\xb7\u2202zQ\u2212\u2202zP\xb7\u2202xQ$, e.g., ${\varphi ,\psi}=\u2212u\xb7\u2207\psi $. Here, *k _{x}* 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 $5\u2009000$ 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=5\u2009000$ becomes $Rm=5\u2009000\xd7Lx\u22481.5\xd7105$, where

*L*represents the box-size along the mean flow direction. The external body force, $f=f(kx=0,z,t)x\u0302$, 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,

_{x}^{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 gravitation

^{50}tends to build shear layers. We assume here such a force, represented as a Krook operator,

^{51–53}as

where $DKrook$, sometimes also referred to as the profile relaxation rate,^{54} measures the forcing strength (in units of $U0/a$) and $\u27e8U(x,z,t)\u27e9x$ 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 *F*_{0} is implemented only to balance the viscous diffusion of the initial shear layer: $Re\u22121\u22072Uref(z)+F0=0$ ensures an initial equilibrium state to which small-amplitude perturbations are added before the system is evolved.

### B. Initial and boundary conditions

As in the unforced study,^{44} a simulation box of $Lx=10\pi $ is considered, but double the size along the *z*-axis ($Lz=20\pi $), 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\pi ,8\pi )$, 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=\xb1Lz/2$.^{44,45}

The initial equilibrium state is seeded with small-amplitude perturbations $(\varphi \u0303,\psi \u0303)$ at all Fourier wavenumbers, as^{44}

and

Here, $A\varphi $ and $A\psi $ 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\varphi (kx)$ and $r\psi (kx)$, forming a uniform distribution in $[0,2\pi )$, are issued for each different *k _{x}* 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\varphi =A\psi =10\u22123$. This set of parameters offers distinct linear and nonlinear phases of evolution.

## III. LINEAR EIGENMODES

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}

### A. Complete eigenspectrum

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

Fourier transforming along the *x*-axis and assuming time variation at each Fourier wavenumber take the form $\varphi \u0302(kx,z,\omega )ei\omega (kx)t$, Eqs. (6a) and (6b) become

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\gamma (kx)t$ and $e\u2212\gamma (kx)t$, respectively. This mode-pair is shown, for the first Fourier wavenumber $kx=2\pi /Lx=0.2$, in Fig. 1(a), along with all the purely real eigenvalues. The latter constitute the eigenmode continuum^{59} 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 $\omega /kx+Uref(z)\xb1vA,ref(z)=0$, where $vA,ref(z)$ is the Alfvèn speed along the reference magnetic field at the vertical coordinate *z*.

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 $\varphi 1(kx,z)$ into the stable mode $\varphi 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.

### B. Competing roles of unstable and stable modes

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.

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}

## IV. NONLINEAR EVOLUTION

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.

### A. Finite-amplitude Kelvin–Helmholtz instability

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.

### B. Momentum transport

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,

where $U=U(x,z,t)$ represents the instantaneous flow, $\u27e8\xb7\u27e9x$ signifies *x*-averaging operation, and *τ _{u}* and

*τ*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

_{b}*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

These stresses are evaluated from nonlinear simulations and shown in Fig. 4. Fluctuations of Reynolds stress are concentrated in the shear layer, near $z\u22480$. 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.

## V. DECOMPOSITION OF NONLINEAR SIMULATION ONTO LINEAR MODES

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 $\chi \u0303turb=(\varphi \u0303turb,\psi \u0303turb)$, which is expanded as

where the eigenmode basis $\chi j(kx,z)$ is employed along the *z*-axis at each wavenumber *k _{x}* to decompose the fluctuations. The complex mode-amplitude $\beta 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–45} *j *=* *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 *k _{x}*.

### A. Nonlinear excitation of stable modes

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.

The energy in individual eigenmodes $|\beta j|2$, averaged over a turbulent state ($t=150\u22121000$), 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.

### B. Reduced representation of the turbulent flow

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 $\varphi \u0303approx$ can be constructed from a class of eigenmodes at each wavenumber, e.g., $\varphi \u0303approx(x,z,t)$ can be written as a sum of an unstable mode per wavenumber $\u2211kxeikxx\beta 1(kx,t)\varphi 1(kx,z)$ or as a sum of an unstable and a stable mode per wavenumber $\u2211kxeikxx[\beta 1(kx,t)\varphi 1(kx,z)+\beta 2(kx,t)\varphi 2(kx,z)]$. Respective short-hand notations $\beta 1\varphi 1$ and $\beta 1\varphi 1+\beta 2\varphi 2$ will be used hereafter, i.e.,

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,\u2009and\u20090.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 $\gamma \u22121\u22485$) 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 $t\u224830$). In this respect, the forced shear layer is markedly different from the freely evolving layer.

### C. Performance of reduced representations

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:

where $(\u2202x\varphi )2$ and $(\u2202z\varphi )2$ are the squared *z*- and *x*-components of velocities; $\varphi \u0303diff=\varphi \u0303exact\u2212\varphi \u0303approx$ with $\varphi \u0303exact$ and $\varphi \u0303approx$ 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., $t\u227230$). 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.

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.

### D. Competing up- and down-gradient momentum transport and their reduced models

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 *k _{x}*, the Reynolds stress from all the fluctuations $\varphi \u0302kx$ reads

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

and

where $\varphi \u0302kx$ is the Fourier transform of the streamfunction at wavenumber *k _{x}* and $\varphi j,kx$ represents the

*z*-dependent

*j*th 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., $\varphi 2,kx(z)=\varphi 1,kx\u2217(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),

where

Thus, we obtain

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.

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.

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,\u2009DKrook=2$, and $Re=500$.

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:

where $\u27e8A\u27e9t$ represents a time-averaging operation on *A*.

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

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, $|\beta 1|2\u2212|\beta 2|2$ can be written as $|\beta 1|2(1\u2212|\beta 2|2/|\beta 1|2)$. Since $|\beta 2|2$ has been found to on the same order of $|\beta 1|2$, e.g., see Fig. 12, where $|\beta 2|2/|\beta 1|2$ can range from $\u22480.8$ to $\u22480.95$, yielding $(1\u2212|\beta 2|2/|\beta 1|2)\u22480.05\u2009to\u20090.2$. Therefore, neglecting stable modes can overestimate the transport by a factor of 5 to 20.

### E. Reynolds vs Maxwell stresses

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

where $\psi \u0302kx$ is the Fourier transform of the flux function at a wavenumber *k _{x}* and $\psi j,kx$ represents the

*z*-dependent

*j*th 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 $kx\u227210$ 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}

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.

### F. Scaling law for continuum modes

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}

where $\beta j(kx)$ represents the complex amplitude of an eigenmode *j* at wavenumber *k _{x}* with

*ω*the associated mode-frequency; the nonlinear mode coupling coefficient

_{j}*C*measures the three-wave overlap, which dictates the strength of nonlinear beating between an eigenmode

_{jmn}*m*at wavenumber $kx\u2032$ and an eigenmode

*n*at wavenumber $kx\u2033$, driving an eigenmode

*j*at wavenumber

*k*(with the constraint $kx\u2032+kx\u2033=kx$).

_{x}For the continuum modes, as was mentioned in Sec. III A, their frequencies depend linearly on the wavenumber *k _{x}* as

^{59}$\omega /kx+Uref(z)\xb1vA,ref(z)=0$. This implies $\omega j=\omega \u221dkx$.

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,

where $P=P0+P\u0303$ and $Q=Q0+Q\u0303$ represent two fields (e.g., $\u22072\varphi $ or *ψ* for the present problem), with *P*_{0} representing the *x*-averaged mean component of *P* and $P\u0303$ standing for perturbations. The linear term, e.g., $ikxP\u0302\xb7\u2202zQ0$, which is in spectral space, contains only one perturbed field, whereas the nonlinear term, e.g., $ikx\u2032P\u0302\u2032\xb7\u2202zQ\u0302\u2033$, has two perturbed fields, with $P\u0302\u2032$ and $Q\u0302\u2033$ representing the Fourier-transformed quantities at wavenumbers $kx\u2032$ and $kx\u2033$, respectively. It may be supposed that the derivative $\u2202z$ on the perturbed quantities is roughly on the scale of $|\u2202z|\u223ckx$. (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 $\u2202z$ acting on *Q*_{0} clearly does not produce a factor of *k _{x}*, which is zero for the mean component

*Q*

_{0}. 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\u0302\xb7\u2202zQ0$ and $ikx\u2032P\u0302\u2032\xb7ikx\u2033Q\u0302\u2033$, respectively.

In Eq. (24), expanding the perturbations in the eigenmode basis, e.g., $P\u0302=\u2211l\beta lP\u0302l$ with $P\u0302l$ representing the *l*th 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., $ikx\u2032P\u0302\u2032\xb7ikx\u2033Q\u0302\u2033$ becomes $\u2212kx2P\u0302Q\u0302$; note the linear term has the form $ikxP\u0302\xb7\u2202zQ0$.

In assuming local interaction between the continuum modes spectral space in *k _{x}*, 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\u2032+kx\u2033$ does not allow two (un-)stable modes to beat together to drive a continuum mode at large

*k*, e.g., $kx>2$.

_{x}Continuing with the above assumption, the nonlinear term $\u2212kx2P\u0302Q\u0302$ has one extra *k _{x}* compared to the linear term $ikxP\u0302\xb7\u2202zQ0$. This implies that, for the continuum modes, the nonlinear mode coupling coefficients

*C*are expected to scale as

because the linear term for the continuum modes in Eq. (23) has the eigenfrequency that depends linearly on *k _{x}*, i.e.,

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 $\omega \beta \u223cC\beta 2$. Inserting their asymptotic dependences on *k _{x}*, the amplitudes of continuum modes are found to follow:

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 *k _{x}*, are plotted in Fig. 14(a) as functions of

*k*and eigenmode index

_{x}*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.

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 $\omega /kx=U00$ ($\omega /kx=\u2212U00$), where $U00=1$. These thick bands represent all eigenmodes that have phase speeds $\omega /kx=U00+cvA,0$ ($\omega /kx=\u2212U00+cvA,0$), where $\u22121<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 *k _{x}* and construct a

*j*-averaged spectral index. To this end, we note the amplitude $\beta j(kx=0.2)$ for each mode

*j*at $kx=0.2$ (the first wavenumber in the simulation) and compute a scaled mode-amplitude $\beta \u0302j(kx)$ as

which is expected to falloff with *k _{x}* as $\u223ckx\alpha $. In principle, the spectral index

*α*can depend on the eigenmode index

*j*, but a

*j*-averaged spectral index is sought now, following the procedure:

This *j*-averaged spectral distribution of the mode-amplitudes informs how, on average, each eigenmode amplitude depends on *k _{x}*.

A plot of $\u27e8\u2009log\u2009\u2009\beta \u0302j(kx)\u27e9j$ vs $log\u2009\u2009kx$ 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 *k _{x}* are shown in Fig. 14.

A finding in Fig. 14 is the identification of self-similar cascade of energy to smaller scales (larger *k _{x}*) 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).

## VI. CONCLUSIONS

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 $|\beta 1|$ and $|\beta 2|$ was recently made for ion-temperature-gradient-driven fusion plasma turbulence^{17} 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 approximation^{61} 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-driven^{6} 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.

## ACKNOWLEDGMENTS

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.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**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).

## DATA AVAILABILITY

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

### APPENDIX A: ORTHOGONALITY OF RIGHT EIGENMODE AND MODIFIED LEFT EIGENMODE

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\xi j=\lambda j\xi j$, where *L* is a linear operator whose *j*-th eigenmode is *ξ _{j}* with eigenvalue

*λ*. The generalized eigenvalue problem that we encounter here is

_{j}where *L* is a linear operator, *M* is another linear operator containing the Laplacian operation for our problem, and *X _{j}* is the

*j*th (right) eigenmode with corresponding (right) eigenvalue

*ω*. (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,

_{j}^{66,67}

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,

In the eigenvalue solver in Dedalus, the matrices $L\u2020$ and $M\u2020$ for each wavenumber are passed, and their eigenmodes $Yj\u2217$ and eigenvalues $\sigma j\u2217$ are found. It can be shown that the eigenvalues *σ _{j}* and

*ω*are the same (i.e., $\sigma j=\omega j$), by analyzing Eqs. (A1) and (A2). A modified orthogonality relation between the left and right eigenmodes can now be derived,

_{j}^{66,67}

which means the left and right eigenmodes are orthogonal to each other with a weight factor *M*, as long as their eigenvalues differ ($\sigma j\u2260\omega i$). For numerical computation, it is convenient to define $YjTM=(MTYj)T=ZjT$, where *Z _{j}* is the modified left eigenmode, which is—by construct—orthogonal to the right eigenmode without any weight factor: $ZjTXi\u221d\delta i,j$. Using this relation, the eigenmode coefficients $\beta j(kx,t)$ in the eigenmode expansion of turbulent fluctuations are computed at each wavenumber and at each time.