The Kelvin–Helmholtz (KH) instability of a shear layer with an initially uniform magnetic field in the direction of flow is studied in the framework of 2D incompressible magnetohydrodynamics with finite resistivity and viscosity using direct numerical simulations. The shear layer evolves freely, with no external forcing, and thus broadens in time as turbulent stresses transport momentum across it. As with hydrodynamic KH, the instability here features a conjugate stable mode for every unstable mode in the absence of dissipation. Stable modes are shown to transport momentum up its gradient, shrinking the layer width whenever they exceed unstable modes in amplitude. In simulations with weak magnetic fields, the linear instability is minimally affected by the field, but enhanced small-scale fluctuations relative to the hydrodynamic case are observed. These enhanced fluctuations coincide with increased energy dissipation and faster layer broadening, with these features more pronounced in simulations with stronger fields. These trends result from the magnetic field reducing the effects of stable modes relative to the transfer of energy to small scales. As field strength increases, stable modes become less excited, thus transporting less momentum against its gradient. Furthermore, the energy that would otherwise transfer back to the driving shear because of the stable modes is instead allowed to cascade to small scales, where it is lost to dissipation. Approximations of the turbulent state in terms of a reduced set of modes are explored. While the Reynolds stress is well-described using just two modes per wavenumber at large scales, the Maxwell stress is not.

## I. INTRODUCTION

Shear layers are ubiquitous in space and astrophysical systems, including Earth's magnetosphere,^{1} relativistic jets,^{2} and clouds passing by galactic and circumgalactic gas.^{3} These flows exhibit extremely large Reynolds numbers and, thus, are often susceptible to shear-flow instabilities that can give rise to turbulence. Among those instabilities is the Kelvin–Helmholtz (KH) instability, the canonical shear-flow instability,^{4–6} which is triggered in unstable flow profiles by sufficiently strong flow shear (i.e., it does not require other physical effects for it to be destabilized) and thus can exist in a wide range of systems. Despite the relative simplicity of the instability in idealized systems, the background magnetic fields present in many astrophysical systems can have significant and complex effects on its dynamics. A uniform field in the direction of flow can stabilize KH if the field strength exceeds some threshold.^{4} This stability threshold depends on fluid properties and the flow profile (see Ref. 7 and references therein) but, in the perfectly conducting case, is roughly characterized by the Alfvén velocity (in terms of the flow-aligned component of the field) exceeding the difference in flow velocity on either side of the layer. For weaker magnetic fields, the instability remains but with a reduction in growth rate that depends on the sonic Mach number (see Ref. 8 and references therein). However, despite a uniform magnetic field decreasing the growth rate relative to the hydrodynamic case, magnetohydrodynamic (MHD) simulations show even weak magnetic fields significantly enhance the generation of small-scale fluctuations and the rate at which momentum is transported across the shear layer, causing the layer to broaden at a faster rate.^{7,9,10} While the focus of the present work is on uniform equilibrium fields, it is worth noting that both the linear and nonlinear dynamics are significantly different in the case of a nonuniform magnetic field.^{8,11}

This instability-driven turbulence generally transports momentum, as well as heat and particles, far faster than viscosity and molecular diffusion would alone. Thus, this turbulence can have important effects on systems where it is found. Indeed, transport driven by shear-flow turbulence is often necessary to explain observations^{12} or similarly plays a key role in reduced dynamical models.^{13,14} This motivates the studies of shear-driven turbulence, particularly in pursuit of transport models that can be employed when considering systems too complex for direct numerical simulation (e.g., in stellar evolution codes^{13,14}) With the significant impacts that background magnetic fields can have on turbulent transport in shear-driven turbulence, it is important that such transport models account for magnetic effects. However, reduced models, where transport is assumed to scale with the growth rate of the driving instability and no details of the nonlinear saturation are included, are clearly inadequate in this case as increasing magnetic field strength increases transport while decreasing the instability's growth rate. A primary goal of this work is to explore what details of the dynamics are responsible for this scaling and, thus, might be necessary to include in reduced models to capture key trends accurately.

In the context of turbulence driven by gyroradius-scale instabilities in fusion plasmas, novel reduced transport models,^{15,16} as well as corrections to existing models,^{17,18} have been derived by accounting for the physical mechanisms that saturate the instability. Here, instability saturation refers to the arresting of the exponential growth of fluctuations that are seeded by initially small perturbations. For example, in systems where perturbations grow by drawing energy from an unstable momentum, density, or temperature gradient, exponential growth might cease once perturbations have drawn so much energy from the driving gradient that it relaxes and is no longer unstable; this is sometimes referred to as quasilinear flattening. Particularly in systems with fixed background gradients, saturation might instead occur when the injection of energy by the instability is balanced by the nonlinear transfer of energy to small, dissipative scales (see also Refs. 19 and 20, where magnetic field generation was noted as a saturation mechanism). A third, distinct saturation mechanism involves the transfer of energy to large-scale, linearly stable (damped) modes.^{21,22} These modes are eigenfunctions of the linearized governing equations, and they decay exponentially in the absence of nonlinear energy transfer from other modes. Signatures of their excitation due to nonlinear energy transfer have been measured in dipole-confined plasmas,^{23} and, in this paper, a previously observed feature of hydrodynamic shear layers in laboratory experiments^{24,25} and simulations,^{26} namely, counter-gradient momentum transport, will be identified as a consequence of stable-mode excitation.

In the fusion context, analytical calculations have shown that stable-mode excitation is almost universally a significant contributor to saturation.^{21,27} Subsequent direct numerical simulations have demonstrated that stable eigenmodes not only affect instability saturation, but also remain excited in the ensuing turbulence.^{28,29} Understanding the nonlinear interactions primarily responsible for energy transfer to stable modes^{30} has enabled the development of a variety of reduced models that incorporate stable mode effects.^{15–18} In the case of shear-flow instabilities, the same analytical saturation calculations of Refs. 21 and 27 were applied to a hydrodynamic, KH-unstable system in Ref. 22, showing that stable modes are important in saturating the KH instability and that when they are excited they can significantly affect turbulent momentum transport (see also Ref. 31, where the methods were extended to systems with eigenmodes that cannot be derived in closed form, like the system considered here, and thus the calculation must be done numerically). In Ref. 32, gyrokinetic simulations of an unstable shear flow revealed that stable modes are excited to significant amplitudes in shear-driven turbulence except when heavily damped by a radiative damping term (also called drag or friction in other contexts). There, an external forcing term partially maintained the unstable flow profile against quasilinear flattening, permitting a quasi-stationary state of driven turbulence. The simulations had only two spatial dimensions, and the dynamics were essentially hydrodynamic, with the forcing term^{33,34} reminiscent of the forcing considered in the well-studied Kolmogorov flow problem.^{35–37} In parameter regimes where stable modes were significantly excited, the relevant flow fluctuations could be approximated well by linear combinations of stable and unstable modes alone (neglecting the continuum of marginally stable modes^{38}). A scaling model for a quantity directly related to the Reynolds stress as a function of the forcing was derived in terms of the stable and unstable mode amplitudes and shown by comparison with simulations to be very accurate.

The present work explores the role of stable modes in shear-flow instability saturation and turbulent momentum transport for a system that differs from Ref. 32 in two key regards. First, no forcing terms are included. Thus, no quasi-stationary state is formed, quasilinear flattening is permitted, and the effects of layer broadening on saturation and the ensuing turbulence are investigated. Second, an initially uniform magnetic field in the direction of flow is included. The present study focuses on the weak-field regime, where the growth rate of the instability is only slightly reduced compared to the hydrodynamic case. The system is studied in the MHD framework via direct numerical simulations using the code Dedalus.^{39,40} The simulations are of a two-dimensional (2D), incompressible fluid with finite viscosity and resistivity included explicitly.

In other systems of instability-driven turbulence, stable modes are known to remove energy from fluctuations that would otherwise cascade to small scales.^{30,41} Here, the enhancement of turbulent momentum transport and small-scale fluctuations with increasing magnetic field strength is examined in the context of stable modes to identify whether the enhanced small-scale fluctuations are due to a reduction in stable-mode activity. At large scales, the linearized, dissipationless system features a pair of unstable and stable modes at each horizontal wavenumber.^{4} Unstable modes gain energy from the shear flow and their associated Reynolds and Maxwell stresses transport momentum down its gradient and broaden the layer.^{6} Stable modes return energy to the background, and their stresses transport momentum against its gradient and tend to shrink the layer. Transient instances of counter-gradient momentum transport, resembling those seen in experiments,^{24,25} are observed here and shown to occur whenever stable modes exceed unstable modes in amplitude. As field strength is increased between simulations, this counter-gradient momentum transport becomes weaker and eventually ceases partly because stable modes are less excited with stronger magnetic fields. With less stable mode activity, the energy that would otherwise be returned to the background flow instead cascades to small scales, producing more small-scale fluctuations and a significant increase in energy dissipation relative to the hydrodynamic case. The small-scale fluctuations also produce a down-gradient Maxwell stress that becomes significant for strong initial fields or sufficiently low resistivity. The enhanced layer broadening is, thus, a combination of reduced counter-gradient momentum transport by stable modes and enhanced down-gradient transport by small-scale magnetic fluctuations.

This paper is organized as follows: The system is described in Sec. II, including the equilibrium, governing equations, and the numerical implementation. The dissipationless linear modes are discussed in Sec. III. Nonlinear simulations are presented in Sec. IV, beginning with an overview of the nonlinear evolution of the system, followed by discussions of the effects of layer broadening in Sec. IV A, stable-mode excitation and momentum transport in Sec. IV B, and small-scale fluctuations and dissipation in Sec. IV C. Conclusions are presented in Sec. V.

## II. SYSTEM SETUP

### A. Equilibrium, governing equations

We study the evolution of a two-dimensional free shear layer in an incompressible fluid in MHD with finite viscosity and resistivity, with an initially uniform magnetic field in the direction of the flow. Specifically, we consider an initial flow in the horizontal direction $x\u0302$ that varies in the vertical direction $z\u0302$, i.e., $V\xaf0=U\xaf(z\xaf)x\u0302$, where $U\xaf(z\xaf)=U\xaf0tanh(z\xaf/d\xaf)$ is the initial flow profile, $U\xaf0$ is the flow speed away from the layer, and $d\xaf$ is the layer half-width, with an initial, uniform magnetic field $B\xaf0=B\xaf0x\u0302$. Here, an overbar denotes a dimensional quantity. Henceforth, we non-dimensionalize all speeds, distances, and fields according to $U=U\xaf/U\xaf0,\u2009(x,z)=(x\xaf/d\xaf,z\xaf/d\xaf)$, and $B=B\xaf/B\xaf0$, respectively, such that $V0=tanh(z)x\u0302$ and $B0=x\u0302$. All other physical quantities will be non-dimensionalized in terms of $U\xaf0,\u2009d\xaf,\u2009B\xaf0$, and combinations thereof.

We describe the flow velocity and magnetic field in terms of a streamfunction $\varphi $ and a flux function *ψ*, so that $v=y\u0302\xd7\u2207\varphi $ and $B=y\u0302\xd7\u2207\psi $. Under our chosen non-dimensionalization, we may write the governing equations as^{42}

and

Here, $MA$ is the Alfvén Mach number, or the ratio of the equilibrium flow speed to the Alfvén speed, and scales like $MA\u221dU\xaf0/B\xaf0$; the Reynolds number $Re$ and magnetic Reynolds number $Rm$ are defined as $Re=U\xaf0d\xaf/\nu \xaf$ and $Rm=U\xaf0d\xaf/\mu \xaf$, respectively, where $\nu \xaf$ is the kinematic viscosity and $\mu \xaf$ is resistivity; ${f,g}\u2261\u2202xf\u2202zg\u2212\u2202xg\u2202zf$. Equation (1) describes the evolution of the vorticity $\u2207\xd7v=\u22072\varphi y\u0302$. The second term on the left-hand side is the vorticity advection term, the first term on the right-hand side is the curl of the Lorentz force, and the second term on the right-hand side is standard viscous dissipation. The terms on the right-hand side of Eq. (2) correspond to flux advection and resistive diffusion. This system, with the above equilibrium, is known to be linearly unstable for $MA$ above a critical threshold that lies between 1 and 2.^{7}

### B. Perturbation equations

As will be described in Sec. II C, we solve Eqs. (1) and (2) numerically using the initial value problem capabilities in the Dedalus code.^{39} Additionally, we use Dedalus' eigenvalue problem capabilities to calculate the complex frequencies (eigenvalues) and eigenmodes of these equations linearized about an unstable equilibrium. Solving the linearized system as an eigenvalue problem allows the full set of eigenmodes at each *k _{x}*, including stable modes, to be calculated, whereas initial value calculations yield only the most unstable mode at each

*k*. This is necessary to track the amplitudes of these modes in the ensuing turbulence in solutions of Eqs. (1) and (2), which then informs how much energy they remove from fluctuations.

_{x}^{30}As with previous studies of stable modes in shear-flow turbulence,

^{22,32}we are specifically interested in the dissipationless modes of this system, so eigenmodes are calculated with viscosity and resistivity neglected. While eigenmodes could be calculated with dissipation included, such modes would mix together the physical effects of conservative energy transfer between the shear flow and fluctuations—a primary focus of this paper—and non-conservative dissipation at large scales. The stable modes of the dissipative system owe their stability to a combination of these two effects, and thus these modes do not lend themselves as conveniently to calculations involving the conservative effect alone. Furthermore, the dissipationless modes have previously been shown to still be relevant in the full, dissipative system.

^{32,43}

To derive linearized equations, we separate the system into a horizontal, uniform (in *x*) background flow *U*(*z*) and field $Bx(z)$, and perturbations $\varphi \u0303$ and $\psi \u0303$. This allows Eqs. (1) and (2) to be similarly separated into equations describing the background, and the following equations for the perturbations:

and

where primes denote derivatives with respect to *z*, and we have neglected viscosity and resistivity. Equations (3) and (4) are the MHD equivalent of Eq. (1) in Ref. 22 and describe how fluctuations interact linearly with the background flow and field, and nonlinearly with one another. When the fluctuations are small enough that the nonlinearities can be neglected, Fourier transforming Eqs. (3) and (4) in *x* and assuming solutions vary in time as $exp\u2009[i\omega (kx)t]$ yields

and

where $\varphi \u0302$ and $\psi \u0302$ are the Fourier transforms (in *x*) of $\varphi \u0303$ and $\psi \u0303$. Thus, at every *k _{x}*, we have a separate system of linear, ordinary differential equations in

*z*. For a given

*k*and $MA$,

_{x}*U*(

*z*) and $Bx(z)$, and appropriate choice of boundary conditions, this system forms a generalized eigenvalue problem that can be solved to obtain a spectrum of eigenvalues

*ω*and eigenmodes $f\u2192j\u2261(\varphi j(z),\psi j(z))$, where $j=1,2,\u2026$ enumerates the different solutions at each

_{j}*k*. While the equilibrium considered in this paper is specifically $U=tanh(z)$ and

_{x}*B*= 1, the more general equations are presented here because eigenmodes corresponding to other

_{x}*U*(

*z*) and $Bx(z)$ will be considered as well in this paper.

### C. Numerical implementation

Dedalus is a pseudo-spectral code with a variety of spectral bases available. We employ Fourier modes $exp\u2009[ikxx]$ in the *x* direction and Chebyshev polynomials $Tn(z)$ in *z*. Our simulation domain size is $Lx\xd7Lz=10\pi \xd710\pi $; thus, the minimum horizontal wavenumber is $kx=0.2$, with periodic boundaries at $x=\xb1Lx/2$ and perfectly conducting, no-slip, co-moving (with the equilibrium flow $V0$) walls at $z=\xb1Lz/2$. The simulations presented here use a resolution of $Nx\xd7Nz=512\xd72048$, with convergence tests performed at the highest values of $Rm$ by ensuring that changes in spectral energy density and dissipation with resolution are minimal. For many of the physical parameters studied here, this *z*-resolution is higher than necessary for well-resolved simulations.

Previous work has shown that the nonlinear development of KH-unstable flows depends sensitively on the choice of the initial perturbations that seed the instability.^{44} In studying free shear layers, a common choice of initial condition is a perturbation in one or more velocity fields that is sinusoidal in *x*, with a wavelength that matches the box size or the fastest-growing linear mode, and Gaussian in *z* centered about the shear layer,^{7,40} with lower-amplitude noise sometimes added to other horizontal wavenumbers.^{10} Here, we perturb both $\varphi $ and *ψ* at every nonzero *k _{x}* with Gaussians in

*z*that have randomly assigned,

*k*-dependent complex phases and amplitudes that decrease with

_{x}*k*as a power law. Thus, at

_{x}*t*=

*0, the streamfunction and flux function are*

and

where the *k _{x}* = 0 components are the unperturbed equilibrium profiles, $A\varphi $ and $A\psi $ set overall amplitudes for the perturbations,

*a*sets the steepness of the energy spectra of the perturbations,

*σ*sets the width of the Gaussian in

*z*, and at every nonzero wavenumber, $\Delta \varphi (kx)$ and $\Delta \psi (kx)$ are uniformly distributed pseudo-random numbers in $[0,2\pi )$. For the results presented here, we use

*σ*= 2,

*a*= −1, and $A\varphi =A\psi =5\xd710\u22124$, which allows for a clearly defined regime of linear growth before nonlinear interactions become important. For $MA=5$, setting $A\psi =0$ did not noticeably change how the instability saturated. This is likely a result of the flow-dominated nature of the instability at these values of $MA$ (as will be shown in Sec. III) and the well-defined linear growth regime permitted by our small value of $A\varphi $.

Previous work^{44} has shown that even when only two wavenumbers are perturbed, the details of the nonlinear stage after the instability saturates are sensitive to the complex phase differences and relative amplitudes between different *k _{x}*, the overall amplitude of the perturbation, and the structure in

*z*of the perturbations. In this work, we are interested in studying details of the saturated state as $MA$ and $Rm$ are varied. In an effort to ensure that our observed trends are not a unique feature of a particular choice of initial conditions, we perform multiple simulations at each $MA$ and $Rm$, with different ensembles of $\Delta \varphi $ and $\Delta \psi $. In practice, this is done by selecting different seeds for our pseudo-random number generator (we use numpy.random.RandomState

^{45,46}to ensure consistency across different computers) and using the same seeds for different $MA$ and $Rm$ so that $MA$ and $Rm$ can be varied independently with $\Delta \varphi $ and $\Delta \psi $ held fixed. For each value of $MA$ and $Rm$ presented here, at least five different sets of initial conditions were simulated. While the majority of this paper presents results from only one set of initial conditions, the trends we present were robust and broadly representative of the range of initial conditions we sampled.

## III. EIGENMODES AND EIGENVALUES FOR $U=tanh(z)$ AND *B*_{x} = 1

_{x}

For $U=tanh(z)$ and *B _{x}* = 1, unstable modes, solutions to Eqs. (5) and (6) with positive growth rates $\gamma j=\u2212Im[\omega j]$, are observed as expected for wavenumbers in the range $0<kx<1$ as long as $MA$ is above a critical threshold between 1 and 2 (the precise value depends on fluid properties and the flow profile, see Ref. 7 and references therein). The growth rate of the fastest-growing mode for this system is plotted against

*k*for a variety of $MA$ in Fig. 1. Magnetic tension provides a stabilizing influence that suppresses instability for $MA$ below the critical threshold and significantly reduces the growth rate for $MA$ slightly above the threshold, but only marginally affects the growth rate for $MA\u22738$.

_{x}Taking the complex conjugate of Eqs. (5) and (6) shows that, as in the hydrodynamic^{5} and gyrokinetic cases,^{32} for every eigenvalue *ω _{j}* and eigenmode $(\varphi j,\psi j)$ that is a solution of Eqs. (5) and (6), the complex conjugate $\omega j*$ and $(\varphi j*,\psi j*)$ is a solution as well (recent work has explored the connection between this conjugate symmetry and parity-time symmetry

^{47}). Following Refs. 22 and 32, when describing the eigenmodes of this system, we label the most unstable mode at each

*k*as

_{x}*j*=

*1 and its conjugate stable mode as*

*j*=

*2. Real-space contours corresponding to $\varphi j$ and*

*ψ*for

_{j}*j*=

*1 and*

*j*=

*2 at $MA=40,kx=0.4$ are shown in Fig. 2. While the flow component of the mode can be roughly described as a superposition of two waves of vorticity localized about the edges of the layer (a wealth of literature exists on this subject*

^{48–51}), the current density of the mode is more localized about the center of the layer.

The source of free energy that drives the exponential growth of the unstable mode is the equilibrium flow *U*(*z*). In terms of Eqs. (3) and (4), the growth of the mode is due to the (dissipationless) linear terms on the on right-hand side. These terms were derived from the energy-conserving nonlinearities in Eqs. (1) and (2) by separating interactions involving *U*(*z*) and $Bx(z)$ from nonlinear interactions between fluctuations. If *U*(*z*) and $Bx(z)$ are identified as the horizontally averaged flow and field and held fixed in time, so that only $kx\u22600$ perturbations are allowed to evolve, then this exponential growth does not conserve energy because the energy injected into $f\u21921$ is not self-consistently removed from the mean flow. This is the case in Ref. 22 as well as previous studies of stable modes in plasma turbulence driven by instabilities aside from KH due to the fixed background gradients that were considered (e.g., Refs. 21 and 27–30). However, in direct numerical simulations of Eqs. (1) and (2) (or in other systems where driving gradients are not held fixed), energy is conservatively transferred from the equilibrium to growing perturbations by the nonlinearities. These considerations are critical for the relationship between the horizontally averaged flow, eigenmodes, and Reynolds stress in counter-gradient transport events described in Sec. IV. Viewed in terms of a separation between the mean and $kx\u22600$ fluctuations, the removal of energy from *U*(*z*) occurs via the *xz* components of the Reynolds and/or Maxwell stress tensors, which we denote as

and

respectively, where $\u27e8\xb7\u27e9x$ indicates an average in *x*. These stresses transport horizontal momentum along the vertical axis and evolve the mean flow according to (neglecting viscosity)

with a transport of momentum down the gradient lowering the kinetic energy of the mean flow.

As with the exponential growth of the unstable mode $f\u21921$, the exponential decay of the conjugate stable mode $f\u21922$ does not conserve energy if the background flow and field are held fixed. Thus, in Ref. 22 and previous studies of stable modes in instability-driven turbulence, stable modes necessarily present a nonconservative energy sink. However, they do conserve energy when directly simulating Eqs. (1) and (2), with energy injection into the mean provided by the same stresses, in addition to a minimal amount of energy that is transferred into the mean field. This is illustrated in Fig. 3, where the Reynolds and Maxwell stresses are shown for $f\u21921$ and $f\u21922$ at $kx=0.4$ for both $MA=25$ and $MA=2$. For unstable modes, both *τ _{u}* and

*τ*transport momentum down the gradient, so that they both contribute to a transfer of energy from

_{b}*U*(

*z*) to $f\u21921$. Likewise, for stable modes, both stresses yield counter-gradient momentum transport, transferring energy from $f\u21922$ to

*U*(

*z*). The transport of the two modes is symmetric in the sense that $\tau u[\varphi 2]=\u2212\tau u[\varphi 1]$ and $\tau b[\psi 2]=\u2212\tau b[\psi 1]$. As $MA$ is decreased, corresponding to a stronger equilibrium field, the relative amplitudes of

*τ*and

_{u}*τ*change, with $|\tau b|$ exceeding $|\tau u|$ for only the strongest equilibrium fields, starting around $MA\u22482.5$.

_{b}## IV. NONLINEAR EVOLUTION

Consistent with previous work,^{7,10} the addition of even a weak magnetic field causes significant changes to the nonlinear evolution of this system despite only slight changes to the linear instability. This is readily seen by inspecting snapshots of the flow. The top row of Fig. 4 shows vorticity and streamlines for a hydrodynamic simulations at three different times, corresponding to precursor vortex formation, vortex merging, and deep in the nonlinear regime. For comparison, the middle and bottom rows of Fig. 4 show vorticity and streamlines, as well as current density and field lines, for an MHD simulation with the same initial conditions. This simulation has an initially weak magnetic field, with $MA=60$ and $Rm=250$. In the hydrodynamic case, the vorticity equation becomes an advection-diffusion equation, so no negative vorticity is produced aside from boundary layer effects. Vorticity conservation is broken by the Lorentz force in MHD, as can be seen by the regions of negative vorticity that form late in the simulation. From close inspection of the vorticity fields, one can infer that vortex reconnection is occurring,^{52} in addition to the magnetic reconnection responsible for changes in field line topology. A thorough comparison of these two varieties of reconnection and their influence on the dynamics is beyond the scope of this paper although the significant dissipation by resistivity that will be shown in Sec. IV C suggests magnetic reconnection is a dominant effect. This is consistent with our observation that magnetic energy has a greater tendency to reach small scales than does the kinetic energy.

Magnetic fields embedded in high-strain-rate flows, such as the narrow vortex sheets, known as braids, connecting the coherent vortices in this system, are known to be amplified by the strain provided $Rm\u226b1$.^{53,54} The flow strain pushes neighboring field lines together, forming a current sheet along the braid as seen in Fig. 4. This figure also demonstrates the familiar amplification of magnetic fields by coherent vortices,^{10,55} with current sheets forming at the edge of the post-merger vortex as it wraps up the magnetic field, with the field lines in the interior of the vortex reconnecting and drifting outwards via flux expulsion.^{56} The field amplification provided by both high-strain-rate flow and coherent vortices increases with $Rm$.^{10,53} Thus, when studying trends with magnetic field strength in this system, not only does field strength vary with $MA$, but also with $Rm$. Decreasing $MA$ corresponds to a stronger initial field, while increasing $Rm$ allows the field to become more amplified as time goes on.

Consistent with Ref. 44, the precursor vortex formation time and the vortex merging time correspond to the initial saturation of $kx=0.4$ and $kx=0.2$, respectively. This is demonstrated in Fig. 5, where various components of energy are plotted over time, with solid and dashed lines correspond to kinetic and magnetic energy for the MHD case, dotted lines correspond to the hydrodynamic case, and different colors correspond to different components of the 1D spectral energy density, defined so that $KE=\u2211kxKEkx$ and $ME=\u2211kxMEkx$. The precursor vortices form roughly when $KEkx=0.4$ reaches its first maximum, and they merge roughly when $KEkx=0.2$ reaches its first maximum. Consistent with Ref. 44, we found the saturation time for $kx=0.4$ to be independent of the complex phase $\Delta \varphi (kx)$ of the initial flow perturbation, while the saturation time for $kx=0.2$ does depend on $\Delta \varphi (kx)$. Varying $\Delta \psi $ had no significant effect on the simulations. The simulations shown here correspond to a choice of $\Delta \varphi (kx)$ with a relatively long merging time.

### A. Layer broadening

The shear layer broadens over time in this system as energy is transferred from the mean flow to fluctuations at $kx>0$. Previous work has shown that the layer broadens more quickly for stronger magnetic fields.^{7,10} This is consistent with the results shown in Fig. 6, where $KEkx=0$ is plotted vs time for simulations with the same initial conditions and $Rm$, but different $MA$. As field strength increases, $KEkx=0$ decays more rapidly overall. From Eq. (11), this implies an overall increase in momentum transport down the gradient in *U*, and hence a broadening of the layer. Close inspection shows brief intervals in time where momentum transport reverses and energy is transferred back to the mean flow, as indicated by transient increases in $KEkx=0$, e.g., near $t\u224855$ and $t\u224885$. While the overall down-gradient transport increases with field strength over long times, the counter-gradient transport in these phases, as well as the down-gradient transport immediately before them, decreases with increased field strength. This will be explored in greater detail in Sec. IV B.

The broadening of the layer has important consequences for the eigenmodes and their impact on transport. While the eigenmodes described in Sec. III transfer energy to and from the initial base flow and field, the modes governing energy transfer with the background change as the background flow and field change.

In unstable shear layers, the growth rate of the instability generally scales with the difference in flow velocity on either side of the layer divided by the layer width. The most-unstable wavenumber and the critical wavenumber above which modes are no longer unstable scales with one divided by the layer width. Thus, as the layer broadens in time, the critical wavenumber decreases from its initial value of *k _{x}* = 1, and the linear growth rate of fluctuations about the mean flow decreases (see also Ref. 25, where the inverse is observed in a system where layer thickness is made to decrease over time). Similarly, the Reynolds and Maxwell stresses corresponding to the linear modes also broaden with the layer.

These trends can be shown directly by solving Eqs. (5) and (6) with $\u27e8U\u27e9x$ and $\u27e8Bx\u27e9x$ in place of *U* and *B _{x}*, where the mean flow and field are taken from individual time steps in nonlinear simulations and assumed to be independent of time. Throughout this paper, eigenmodes, complex frequencies, and growth rates obtained in this manner are, respectively, denoted $f\u2192\u27e8j\u27e9\u2261(\varphi \u27e8j\u27e9(z),\psi \u27e8j\u27e9(z)),\u2009\omega \u27e8j\u27e9$, and $\gamma \u27e8j\u27e9$ to distinguish them from the eigenmodes of the equilibrium described in Sec. III, with

*j*=

*1 and 2 continuing to denote the most-unstable and conjugate stable modes, respectively, at each wavenumber. Note that Eqs. (5) and (6) still admit a conjugate stable mode for every unstable mode even when using the mean flow and field; thus, $\gamma \u27e82\u27e9=\u2212\gamma \u27e81\u27e9$ still holds. Figure 7 shows how growth rates $\gamma \u27e8j\u27e9$ evolve over time for the most unstable mode at the four initially unstable wavenumbers for $MA=40$. Even when the mean flow has hardly evolved in the first few time steps, the growth rates begin to decline noticeably, particularly at the higher wavenumbers. The highest wavenumbers stabilize first as the critical wavenumber decreases from*

*k*= 1. The linear growth regime for $KEkx=0.4$ is seen in Fig. 5 to end at about the same time that $\gamma \u27e81\u27e9(kx=0.4)$ approaches zero, suggesting that quasilinear flattening is a dominant saturation mechanism for this mode. The same is true for $kx=0.6$ and 0.8 (not shown in Fig. 5). The structure in

_{x}*z*of these modes and their corresponding Reynolds and Maxwell stresses are largely the same as $f\u21921$ and $f\u21922$, except that they broaden with the shear layer.

In freely evolving shear layers, the mean flow has been shown to depart from a simple, broadened $tanh$ profile in both MHD simulations^{7,10} and hydrodynamic experiments [see Eq. (5.2) and Fig. 2 in Ref. 57], and Ref. 58 (see their Table I) has shown that even minor departures from a broadened $tanh$ profile cause significant changes to eigenmode structures and the critical wavenumber. In solving for $f\u2192\u27e8j\u27e9$ and $\omega \u27e8j\u27e9$, minor features in $\u27e8U\u27e9x$ and $\u27e8Bx\u27e9x$ can have a significant impact in the full set of linear modes, including introducing new unstable and conjugate stable modes with finite real frequency. These modes are localized to regions other than *z *=* *0, including to new inflection points in $\u27e8U\u27e9x$. At times, changes in $\u27e8Bx\u27e9x$ separately introduce new modes as well, consistent with the understanding that nonuniform fields can destabilize shear flows in the absence of inflection points.^{59} These new modes can form in addition to the existing unstable KH mode, can replace the original conjugate pair of modes with two conjugate pairs of finite-frequency modes, or can emerge after the KH modes have already stabilized at that wavenumber, such as the $kx\u22650.4$ unstable modes that emerge around $t\u224860$ in Fig. 7. Each mode's complex frequency is well within the bounds of the modified semicircle theorems derived in Ref. 60, but these modes often bear little resemblance to the modes described in Sec. III. A detailed investigation of these modes, such as their scaling with different system features, their effects on transport, or their use in reduced models, is beyond the scope of this work. We note the existence of these modes, however, to point out that while the eigenmodes $f\u2192\u27e8j\u27e9$ of the broadened system correspond more directly to energy transfer to/from the mean than the eigenmodes $f\u2192j$ of the equilibrium, analyses in terms of these modes are often unwieldy because of their complexity. Furthermore, as Subsection IV B will show, the modes of the equilibrium lend themselves to analyses later into the simulation than the modes of the broadened flow.

### B. Large-scale eigenmode excitation and momentum transport

At every time *t* and wavenumber *k _{x}*, the Fourier-transformed system state $f\u0302(kx,z,t)\u2261(\varphi \u0302(kx,z,t),\psi \u0302(kx,z,t))$ can be expressed as

provided the eigenmodes of the equilibrium, ${f\u2192j(kx,z)}$, form a complete basis. The complex-valued $\beta j(kx,t)$ is the amplitude of mode $f\u2192j(kx,z)$ at time *t* and can be understood as the coefficient of the state vector $f\u0302$ expressed in this basis. The eigenmodes of the horizontally averaged system, ${f\u2192\u27e8j\u27e9(kx,z)}$, can also be used as a basis, and the amplitudes of these modes will be denoted as $\beta \u27e8j\u27e9(kx,t)$ here. Note that both bases are complete (specifically, over the finite-dimensional space of fluctuations in this numerically discretized system); thus, the corresponding mode amplitudes are uniquely defined and independent of choice of inner product.^{31} The procedure for calculating the mode amplitudes *β _{j}* is essentially identical to that employed in Refs. 28, 29, and 32, except that the Laplacian on the left-hand side of Eq. (3) must be taken into account. The same methods are also employed in Ref. 61.

Figure 8 shows mode amplitudes over time for the most unstable and conjugate stable mode for the four initially unstable wavenumbers for the same simulation as in Fig. 7, calculated using both sets of modes. At each time, $\beta \u27e8j\u27e9$ is obtained by expanding $f\u0302$ in terms of the eigenmodes ${f\u2192\u27e8j\u27e9}$ of the instantaneous mean flow and field. Leading into saturation, both sets of amplitudes, ${\beta j}$ and ${\beta \u27e8j\u27e9}$, evolve as expected for systems where unstable modes nonlinearly drive stable modes:^{21,22,27,32} the stable modes decay linearly before being nonlinearly driven while unstable modes are still growing linearly. Comparing the $|\beta j|$ and $|\beta \u27e8j\u27e9|$ curves demonstrates some of the differences between the two eigenmode bases. Two notable features indicate that $f\u2192j$ fails to capture the dynamics as precisely as the more relevant $f\u2192\u27e8j\u27e9$.

First, there are periods where $|\beta 2|$ grows at the same rate as $|\beta 1|$, when it is expected to and eventually does grow faster. These can be understood as follows: suppose the true stable mode amplitude in some system evolves as $\beta 2(t)=\beta 2(0)e\u2212|\gamma |t+\beta 12(t)$, and the unstable mode amplitude as $\beta 1(t)=\beta 1(0)e|\gamma |t$ [see, for example, Eqs. (22)–(25) in Ref. 21]. If a mode that differs slightly from $f\u21922$ is used to calculate $\beta 2(t)$ from simulation data, then an error of the form $\u03f5\beta 1$ will almost always be introduced. This will cause the apparent stable mode amplitude to briefly evolve as $\beta 2(t)\u223c\u03f5\beta 1(t)$ before the $\beta 2(t)\u223c\beta 12$ term becomes dominant, provided that $|\u03f5|\u2273|\beta 1(0)|$. For even small differences in $f\u21922$ (measured by some inner product or its corresponding norm), *ϵ* can still be quite large provided the eigenmodes are nonorthogonal (under this choice of inner product). Hence, while the differences between the two sets of modes are initially extremely small, they can cause the parametric driving of stable modes by unstable modes to be overlooked in these analyses.

A second subtle effect that is overlooked by the equilibrium modes can be seen when $|\beta \u27e82\u27e9|$ decreases briefly and dramatically before quickly increasing again. This can be attributed to a difference in the complex phases of different interactions (either with the background flow, or with other modes) that determine $\u2202\beta \u27e82\u27e9/\u2202t$:^{21,22,27} when one interaction overtakes another and becomes dominant, if the two have approximately opposite complex phases, then $\beta 2(t)$ briefly passes through or near 0 in the complex plane before growing in amplitude.

Note that these differences between the two sets of amplitudes begin quite early in the simulation (before *t *=* *10 for $kx=0.4$, see Fig. 8). This is despite almost unnoticeable departures of $\u27e8U\u27e9$ and $\u27e8Bx\u27e9$ from $U=tanh(z)$ and *B _{x}* = 1, respectively, at these early times. Thus, the mode amplitudes are sensitive to small perturbations in the linear system defined by Eqs. (5) and (6), as are the growth rates (noted in the discussion of Fig. 7). This sensitivity is a characteristic feature of non-normal linear operators (see Ref. 62 for a comprehensive review on this subject). Indeed, in the viscous, hydrodynamic case, the linear system is known to have significant properties that can be traced back to its non-normality,

^{63}and insight has similarly been gained into the magnetorotational instability by considering consequences of non-normality.

^{64,65}While Figs. 7 and 8 show some consequences of non-normality, we stress that the excitation of stable modes is not a linear effect: it is a result of nonlinear energy transfer from unstable modes and is a consequence of significant coupling between these modes through the nonlinearity,

^{21,22,27}rather than a consequence of non-normality. Furthermore, the non-orthogonality of a set of eigenmodes, and common metrics for the consequences of non-normality (non-trivial

*ε*-pseudospectra and nonmodal growth), depend on choice of inner product and norm, while the mode amplitudes we present here do not. We encourage readers interested in further discussion of stable modes and non-normal operators to refer to Hatch

*et al.*,

^{43}where connections between these topics are more thoroughly explored.

While the mean-flow amplitudes $\beta \u27e8j\u27e9$ capture subtler details leading into saturation than the equilibrium amplitudes *β _{j}*, and the corresponding eigenmodes $f\u2192\u27e8j\u27e9$ capture the layer broadening, they become unreliable once $\gamma \u27e81\u27e9\u21920$ for a given wavenumber. Furthermore, when additional modes emerge as discussed in Sec. IV A, identifying broadened modes corresponding to the original stable and unstable modes, and tracking their amplitudes, can become difficult or impossible. Despite the physical effects overlooked by the equilibrium modes, they can still be used to assess stable mode activity and connect it to counter-gradient momentum transport. Their utility in assessing stable mode activity is demonstrated by following similar methods to Ref. 32: construct two approximations for $f\u0302$ by truncating the summation in Eq. (12) to either include only

*j*=

*1, or*

*j*=

*1 and*

*j*=

*2. Comparing the accuracy of these two approximations provides a measure of how significant stable modes are in the turbulent state in a way that is normalized to the overall amplitude of the turbulence (whereas $|\beta 2|$ alone is not a normalized measure).*

Figures 9 and 10 show, as functions of time, the resulting truncation errors (conceptually similar to, but not to be confused with, the familiar truncation errors in numerical simulations with spectral methods^{66}) when fluctuations at $kx=0.2$ and 0.4, respectively, are approximated in this manner for a simulation with *M _{A}* = 40 and $Rm=500$. Black and gray curves compare approximations using

*β*and $f\u2192j$, while blue and orange curves compare approximations using $\beta \u27e8j\u27e9$ and $f\u2192\u27e8j\u27e9$, where $f\u2192\u27e81\u27e9$ and $f\u2192\u27e82\u27e9$ are defined as the mode with the largest growth rate and its conjugate. The truncation errors are calculated separately for $\varphi \u0302$ and $\psi \u0302$ according to

_{j}and

where $||\varphi \u0302||KE2$ and $||\psi \u0302||ME2$ are the kinetic and magnetic energy of a given $\varphi \u0302$ and $\psi \u0302$. Unsurprisingly, the $f\u2192\u27e8j\u27e9$ modes describe flow fluctuations better than the $f\u2192j$ ones until quasilinear flattening stabilizes modes at a given wavenumber, at which point such approximations become ill-defined. However, even when quasilinear flattening has set in, flow approximations using eigenmodes of the initial equilibrium retain some fidelity. More importantly, the observation that $\u03f5T[\varphi \u0302]$ is reduced when stable modes are included holds true when either set of modes is used. Thus, the equilibrium mode amplitudes *β _{j}* serve as sufficient indicators of stable mode activity despite the physical effects they overlook. Figure 11 is identical to Fig. 9 but for a simulation with $MA=7.5$. Here, including stable modes does not reduce $\u03f5T[\varphi \u0302]$ as dramatically as in the $MA=40$ case (using either set of modes), suggesting stable modes play less of a role as field strength is increased. Note that in all three figures, each approximation fails to describe the fluctuating field well, particularly after the linear regime.

Figures 9–11 assess the ability of the equilibrium modes $f\u2192j$ to describe the state $f\u0302$ over the entire domain in *z*, demonstrating that they do assess stable mode activity despite the evolving mean flow and field. Focusing on the center of the shear layer, *z *=* *0, demonstrates that they describe fluctuations there particularly well and ties stable mode activity to counter-gradient momentum transport. This is shown in Fig. 12, where dashed curves show the Reynolds stress at the center of the layer due to fluctuations at different wavenumbers, $\tau u(kx,z=0)$, over time for three different $MA$. The quantity $|\beta 2|2\u2212|\beta 1|2$ [cf. Eq. (23) in Ref. 22 and Eq. (15) in Ref. 32] is plotted for comparison. In each case, the two curves essentially overlap and are nearly identical up to a constant of proportionality that depends on *k _{x}* but not $MA$. The equilibrium modes $f\u2192j$ capture the Reynolds stress at the center of the layer well, and doing so requires only the stable and unstable mode at each wavenumber. Furthermore, whenever stable mode amplitudes exceed unstable mode amplitudes at a given wavenumber, the mid-layer Reynolds stress at that wavenumber changes sign and the stable mode drives counter-gradient momentum transport. Conversely, whenever the mid-layer Reynolds stress changes sign, stable modes exceed unstable modes in amplitude.

From Fig. 12 alone, it is unclear whether the general decrease in $||\beta 2|2\u2212|\beta 1|2|$ with decreasing $MA$ is due to a decrease in both mode amplitudes or whether they are instead trending toward equipartition. The first case might reflect a decrease in transport due to an overall decrease in fluctuation amplitudes, while the second reflects a decrease in transport that is independent of fluctuation amplitudes. Figure 13 compares $|\beta 1|$ and $|\beta 2|$ at $kx=0.2$ for $MA=40$ and $MA=15$, showing that while stable mode amplitudes broadly decrease with increasing field strength, unstable mode amplitudes do as well. As field strength increases, the two modes tend toward equipartition. Thus, the reduction in large-scale Reynolds stresses with decreasing $MA$ in Fig. 12 is not purely a result of decreased fluctuation amplitudes. With stronger magnetic fields, the large-scale eddies arrange themselves in a manner that reduces their associated Reynolds stress.

We have established that stable modes are directly responsible for reversals in the mid-layer Reynolds stress of low-*k _{x}* fluctuations. To demonstrate that

*τ*is the dominant contributor to the total momentum transport at high $MA$ and solely responsible for the local minima in $KEkx=0$ identified in Sec. IV A, we compare

_{u}*τ*and

_{u}*τ*. This is done at four different times for a simulation with $MA=60$ and $Rm=250$ in Fig. 14. The four times shown are shortly before and after the first two local minima in $KEkx=0$. At each time, the transport is separated into the Reynolds stress

_{b}*τ*and Maxwell stress

_{u}*τ*, which are themselves separated into contributions from different

_{b}*k*. At all times, $|\tau u|$ exceeds $|\tau b|$ by over an order of magnitude. The Reynolds stress is dominated first by $kx=0.4$ and then by $kx=0.2$, the same scales that contain the majority of the kinetic energy in $kx>0$ fluctuations. Thus, not only are the local minima in $KEkx=0$ due to reversals in

_{x}*τ*and not

_{u}*τ*, they are specifically due to

_{b}*τ*reversals at small wavenumbers, where we have identified stable-mode excitation as the effect that causes these reversals. This definitively ties counter-gradient momentum transport in this system to stable-mode excitation.

_{u}In contrast to the Reynolds stress, the Maxwell stress *τ _{b}* is dominated by scales beyond the initially unstable range, aside from the very first panel in Fig. 14. This is similarly consistent with the distribution in

*k*of magnetic energy. While the low-

_{x}*k*Reynolds stresses resemble the Reynolds stress of the eigenmodes $\tau u[\varphi j]$—albeit broadened, consistent with the broadening of the flow profile and eigenmodes discussed in Sec. IV A—the low-

_{x}*k*Maxwell stresses bear less of a resemblance to $\tau b[\psi j]$, or to $\tau b[\psi \u27e8j\u27e9]$, particularly as time goes on. Furthermore,

_{x}*τ*was never observed to change sign for the simulations considered here. When counter-gradient momentum transport occurs, it is always due to

_{b}*τ*changing sign. This is consistent with the truncation error for the field, $\u03f5T[\psi \u0302]$, generally exceeding the truncation error for the flow, $\u03f5T[\varphi \u0302]$: if stable and unstable modes alone captured the Maxwell stress as well as they capture the Reynolds stress, then

_{u}*τ*at low

_{b}*k*would reverse sign whenever

_{x}*τ*does. Instead,

_{u}*τ*consistently transports momentum down the gradient. This

_{b}*τ*-driven transport is not captured by a simple Fick's law although more sophisticated models, such as a magnetic eddy viscosity model,

_{b}^{67}which applies to the same parameter regime studied here, may perform better.

While counter-gradient momentum transport has been firmly tied to stable mode activity in this paper, the weakening and eventual suppression of counter-gradient transport events with increased field strength is only partly due to magnetic fields reducing stable mode activity. This consistent, down-gradient *τ _{b}* becomes stronger for stronger magnetic fields and can reduce or even cancel the effect of any counter-gradient

*τ*. This is seen in Fig. 15, where

_{u}*τ*and

_{u}*τ*are shown for different $MA$, focusing on the first two counter-gradient transport phases. As field strength increases, $|\tau b|$ increases. This trend is seen for both a stronger initial field from a decreased $MA$ or with more field amplification from an increased $Rm$ (not shown, however trends are qualitatively similar, but the high-

_{b}*k*contributions become progressively more dominant at higher $Rm$). This reduces the net counter-gradient transport during these phases and increases down-gradient transport at other times. In the first phase of counter-gradient transport,

_{x}*τ*remains mostly unchanged with $MA$, as the relative amplitudes between stable and unstable modes does not significantly change with $MA$ at this time (see Fig. 12). So while stable mode activity does become less prominent at this time as field strength increases, it is not due to a clear suppression of stable modes relative to unstable ones (cf. Ref. 32) and a corresponding decrease in counter-gradient

_{u}*τ*. Instead, the influence of stable modes has become less prominent due to an increase in energy transfer to small scales leading to increased $|\tau b(kx)|$ for large

_{u}*k*—that is, the stable modes play a less prominent role only because small-scale magnetic fluctuations are playing more of a role. For $Rm=500$, this holds true for all simulations with $MA\u22735$. For lower $MA$, the counter-gradient $\tau u(kx=0.4)$ at this time does become partially reduced and then fully removed as $MA$ decreases, and

_{x}*τ*remains down-gradient but becomes dominated by larger scales. In the second phase of counter-gradient transport, the dominant contributor to the Reynolds stress, $kx=0.2$, becomes noticeably weaker with stronger fields. This trend in

_{b}*τ*amplitudes is similar to what happens during the down-gradient transport phases at

_{u}*t*=

*40 and 85: $|\tau u|$ becomes weaker as field strength increases at later times, but does not change with $MA$ at earlier times.*

From these results, the overall increase in down-gradient transport and decrease in $KEkx=0$ with field strength can be interpreted as a combination of two factors: first, reduced counter-gradient transport from large-scale flow fluctuations, and second, enhanced down-gradient transport from small-scale field fluctuations. The former is a direct result of decreased stable mode amplitudes. The latter can be understood as an indirect result of stable modes playing a less prominent role in this system. When stable modes activity is significant, energy that would otherwise cascade to small scales is instead returned to the mean flow. As previous work has noted,^{30} a reduction in stable mode activity allows more energy to reach small scales. The enhanced small-scale fluctuations are further explored in Subsection IV C.

### C. Small-scale fluctuations, dissipation

The enhanced generation of small-scale fluctuations as magnetic field strength increases has been noted in terms of the distribution in *k _{x}* of the Maxwell stress (Fig. 15) and of energy (Fig. 5), particularly magnetic energy. When stable mode activity becomes less prominent, the transfer of energy to small, dissipative scales is enhanced. This increases the energy dissipation rate.

Figure 5 shows that the total energy in the system decreases more rapidly in the MHD simulations than in the hydrodynamic one, and that it decreases more rapidly as field strength increases. The enhanced transfer of energy to small scales coincides with an overall increase in energy dissipation. While the magnetic field does increase the kinetic energy at small scales, the increased dissipation is almost entirely due to the addition of resistive dissipation. This is seen in Fig. 16, where viscous and resistive dissipation are plotted for a variety of simulations, each with the same initial condition but different $MA$ and $Rm$. For simulations where the amplified field is not sufficiently strong, the viscous dissipation remains essentially identical to that of the hydrodynamic case. This threshold depends on both $MA$ and $Rm$. For $MA=60$, viscous dissipation does not increase until $Rm$ reaches 1000, while for $MA=40$ the $Rm=500$ simulation does exhibit enhanced viscous dissipation while the $Rm=250$ case (not shown) does not. The precise threshold for this change in behavior, and a detailed investigation of its cause, is beyond the scope of this paper. Resistive dissipation, on the other hand, matches or even exceeds viscous dissipation, consistent with magnetic energy exceeding kinetic energy at small scales. As expected from Ref. 53 [see Eq. (6) there], the local maxima in resistive dissipation are greater and occur later with increasing $Rm$ as the current sheets formed by the high-strain-rate flow take longer to reach small enough scales that resistivity overtakes flux advection. This trend is mostly robust to changes in initial conditions. However, for initial conditions where the merging of the two initial vortices occurs rapidly, the first peak in resistive dissipation may be determined by the time when the merging process disrupts the first current sheet along the braid between the two vortices (e.g., the current sheet seen at *t *=* *40 in Fig. 4), rather than being determined by resistive effects disrupting the sheet.^{53} For these reasons, care must be taken when studying dissipative processes in KH-unstable systems in MHD to consider appropriate initial conditions, box sizes that permit mergers if relevant, and explicit viscosity^{40} *and* resistivity^{7} wherever possible. While $MA$ has no clear effect on the time of these local maxima in high-$MA$ simulations, it does affect the overall level of resistive dissipation, with stronger fields (lower $MA$), yielding greater resistive dissipation. This can be understood as a lower $MA$ increasing the amplitude of the magnetic fluctuations that reach small scales, thus increasing the amplitude of the small-scale, dissipative currents.

## V. CONCLUSIONS

In turbulence driven by unstable, freely evolving shear layers in MHD, the addition of a magnetic field has been observed to increase the layer broadening rate due to enhanced turbulent momentum transport,^{7,10} and to increase energy transfer to small scales, despite stabilizing the driving instability.^{4} This work has investigated the role of large-scale, dissipationless stable modes in this system to identify whether these trends are caused by a reduction in stable mode activity. These modes transfer energy from large-scale fluctuations back to the driving momentum gradient, shrinking the layer width, removing energy from fluctuations, and impeding the cascade to small scales.^{30,41} The results presented here show that the enhanced transfer to small scales and increased broadening rate with stronger magnetic fields do coincide with reduced stable mode activity. Furthermore, in allowing the shear layer to evolve freely without additional forcing terms, this presents the first investigation into the role of stable modes in a system where the source of instability is not fixed^{21,22,28} or quasi-stationary.^{32}

Allowing the layer to broaden introduces quasilinear flattening as a saturation mechanism and yields two distinct sets of eigenmodes: the modes obtained by linearizing about the equilibrium flow $U=tanh(z)$ and field *B _{x}* = 1, denoted by $f\u2192j$, and those obtained by linearizing about the instantaneous mean flow $\u27e8U\u27e9x$ and field $\u27e8Bx\u27e9x$ at each time step, denoted by $f\u2192\u27e8j\u27e9$. The $f\u2192\u27e8j\u27e9$ modes correspond more directly to energy transfer to and from the mean flow, and they reflect the broadening of fluctuations as the layer broadens. However, they quickly vanish or become ill-defined as the layer evolves and develops small-scale features in

*z*. On the other hand, while the $f\u2192j$ modes do not correspond to the evolved system as directly as the $f\u2192\u27e8j\u27e9$ modes, they remain useful for assessing stable mode activity, and they capture the mid-layer Reynolds stress almost exactly when using only two modes at each

*k*. Thus, the same eigenmode decompositions used in systems with a fixed

_{x}^{28,29}or otherwise quasi-stationary

^{32}unstable profile driving the turbulence are also effective in this system, despite the evolution of the mean flow that drives instability.

These eigenmode decompositions were used to track the amplitudes of stable and unstable modes through the evolving turbulence. In simulations with stronger magnetic fields, stable mode activity becomes less prominent, and the layer broadening rate increases for two reasons. First, stable mode amplitudes are reduced relative to unstable ones in simulations with stronger magnetic fields. This reduces the counter-gradient momentum transport these stable modes drive via turbulent Reynolds stresses. Second, the down-gradient Maxwell stress due to small-scale magnetic fluctuations becomes stronger with stronger fields. This can be understood as a result of stable-mode excitation playing less of a role in instability saturation relative to energy transfer to small scales. When stable modes are less impactful, the energy they would remove from fluctuations instead makes its way to small scales.^{30,41} While the Maxwell stress produced by these small-scale fluctuations does not fit a simple Fick's Law in terms of the mean flow gradient, it is possible that more sophisticated models such as a magnetic eddy viscosity,^{67} might capture this effect.

These enhanced small-scale fluctuations, which are generally magnetic fluctuations, significantly boost the energy dissipation rate even for very weak magnetic fields (more than doubling the dissipation rate relative to the hydrodynamic case for *M _{A}* = 100), with dissipation increasing as $Rm$ increases or

*M*decreases. This increased dissipation is almost entirely due to resistivity although viscous dissipation increases with field strength as well, provided the field is stronger than some threshold that depends on $Rm$ and whose precise value for different $Rm$ is beyond the scope of this paper to determine. Thus, even in the presence of a magnetic field much too weak to significantly affect the linear instability, shear layers in non-ideal MHD can dissipate energy much faster than the hydrodynamic counterpart, with the reduction of stable mode activity providing an underlying explanation.

_{A}The observed variation in dynamics as stable-mode activity changes invites the question of how astrophysical shear flow instabilities saturate, where stable modes play a role, and how their presence or absence affects the system. For example, in galaxy halos, the entrainment of cold gas into galactic outflows is a topic of active research^{68} and is thought to occur in part via shear-flow instabilities. How these instabilities saturate in different parameter regimes likely affects the rate at which passive scalars are mixed, in addition to the momentum mixing studied here. Recent work^{69} found that background magnetic fields dramatically affect the drag between a cold cloud and hot medium, suggesting that this is a scenario where stable-mode activity or lack thereof may play a role. However, that work also noted the importance of 3D effects, which have been neglected in the present study. Furthermore, the dependence on magnetic Reynolds number we observed precludes reliable speculation of the background field strength necessary to suppress stable modes at the large $Rm$ (and large magnetic Prandtl number) found in these systems.

Fingering (or “thermohaline”) convection in stellar interiors^{70} presents a system that does not suffer from these concerns, but where the effects of stable modes are likely more subtle. Simulations show that the fingering instability saturates by driving so-called “parasitic” KH instabilities in both the hydrodynamic^{71,72} and MHD^{73} cases (similar to channel modes of the magnetorotational instability^{74–77}). The length scales of these KH modes are small enough that the Reynolds and magnetic Reynolds numbers are often similar to the values studied here, and turbulent transport models based on these KH parasites are remarkably accurate even when neglecting density gradients and 3D effects.^{71,73} This suggests that the results of the present paper are applicable to fingering convection and that the saturation mechanisms of these KH modes may have significant bearing on turbulent transport. In particular, recent work in the hydrodynamic case^{78} showed that adding a horizontal shear body forcing to fingering convection yields counter-gradient momentum transport in some parameter regimes. The results of the present paper suggest this counter-gradient transport may be suppressed by the addition of a magnetic field similar in strength to the weakest fields considered in Ref. 73: as stable modes become subdominant, smaller-scale magnetic fluctuations may produce significant down-gradient momentum transport. Accurately capturing these effects may require finer numerical resolutions than in the hydrodynamic case.

This work motivates two related possible future directions of inquiry. First, while the Reynolds stresses in this system can be reasonably modeled using stable and unstable modes alone, the Maxwell stresses cannot. Not only do these modes fail to describe large-scale field fluctuations well, but the Maxwell stress is primarily driven by small-scale fluctuations, while the modes considered in this paper exist at large scales. Recent work^{67} involving a quasi-stationary system of driven, shear-flow turbulence has suggested these stresses can be described by a magnetic eddy viscosity model. It is plausible that for the same system considered here, with the addition of a shear forcing term to permit a quasi-stationary state, the momentum transport could be well-described by a combination of stable and unstable modes for the Reynolds stress and a magnetic eddy viscosity model for the Maxwell stress. Furthermore, a quasi-stationary system enables a thorough investigation of the dominant nonlinear interactions in the saturated state, potentially informing saturation theories for predicting eigenmode amplitudes without first performing direct numerical simulations (similar to the saturation theories^{15} that were informed by studies of dominant nonlinear interactions^{30} in the context of fusion plasmas).

Provided a reduced model that is informed by stable mode activity, a second direction for this work is the application of such a model to physical systems. For example, if the effect of density stratification is included, such a model might be used to improve predictions of shear-driven transport in stellar evolution codes.^{13}

## ACKNOWLEDGMENTS

The authors thank K. Burns, D. Lecoanet, C. Sovinec, and J. Parker for insightful discussions, and N. Hurst for insightful discussions and a thorough reading of the manuscript. The authors also thank the Dedalus developers and the Dedalus user group for assistance with many aspects of installing and running Dedalus. Partial support for this work was provided by the National Science Foundation under Award No. PHY-1707236 and by the U.S. Department of Energy, Office of Science, Fusion Energy Sciences under Award No. DE-FG02-04ER54742. Computing resources were provided by the National Science Foundation through XSEDE computing resources, Allocation Nos. TG-PHY130027 and TG-PHY180047. This research was performed using the computer resources and assistance of the UW-Madison Center for High Throughput Computing (CHTC) in the Department of Computer Sciences. The CHTC is supported by UW-Madison, the Advanced Computing Initiative, the Wisconsin Alumni Research Foundation, the Wisconsin Institutes for Discovery, and the National Science Foundation, and is an active member of the Open Science Grid, which is supported by the National Science Foundation and the U.S. Department of Energy's Office of Science.

## DATA AVAILABILITY

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