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.

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 observations12 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 codes13,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 experiments24,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 modes30 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 term33,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 modes38). 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.

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̂ that varies in the vertical direction ẑ, i.e., V¯0=U¯(z¯)x̂, where U¯(z¯)=U¯0tanh(z¯/d¯) is the initial flow profile, U¯0 is the flow speed away from the layer, and d¯ is the layer half-width, with an initial, uniform magnetic field B¯0=B¯0x̂. Here, an overbar denotes a dimensional quantity. Henceforth, we non-dimensionalize all speeds, distances, and fields according to U=U¯/U¯0,(x,z)=(x¯/d¯,z¯/d¯), and B=B¯/B¯0, respectively, such that V0=tanh(z)x̂ and B0=x̂. All other physical quantities will be non-dimensionalized in terms of U¯0,d¯,B¯0, and combinations thereof.

We describe the flow velocity and magnetic field in terms of a streamfunction ϕ and a flux function ψ, so that v=ŷ×ϕ and B=ŷ×ψ. Under our chosen non-dimensionalization, we may write the governing equations as42 

(1)

and

(2)

Here, MA is the Alfvén Mach number, or the ratio of the equilibrium flow speed to the Alfvén speed, and scales like MAU¯0/B¯0; the Reynolds number Re and magnetic Reynolds number Rm are defined as Re=U¯0d¯/ν¯ and Rm=U¯0d¯/μ¯, respectively, where ν¯ is the kinematic viscosity and μ¯ is resistivity; {f,g}xfzgxgzf. Equation (1) describes the evolution of the vorticity ×v=2ϕŷ. 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 

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 kx, including stable modes, to be calculated, whereas initial value calculations yield only the most unstable mode at each kx. 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.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 ϕ̃ and ψ̃. This allows Eqs. (1) and (2) to be similarly separated into equations describing the background, and the following equations for the perturbations:

(3)

and

(4)

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[iω(kx)t] yields

(5)

and

(6)

where ϕ̂ and ψ̂ are the Fourier transforms (in x) of ϕ̃ and ψ̃. Thus, at every kx, we have a separate system of linear, ordinary differential equations in z. For a given kx and MA, 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 ωj and eigenmodes fj(ϕj(z),ψj(z)), where j=1,2, enumerates the different solutions at each kx. While the equilibrium considered in this paper is specifically U=tanh(z) and Bx = 1, the more general equations are presented here because eigenmodes corresponding to other U(z) and Bx(z) will be considered as well in this paper.

Dedalus is a pseudo-spectral code with a variety of spectral bases available. We employ Fourier modes exp[ikxx] in the x direction and Chebyshev polynomials Tn(z) in z. Our simulation domain size is Lx×Lz=10π×10π; thus, the minimum horizontal wavenumber is kx=0.2, with periodic boundaries at x=±Lx/2 and perfectly conducting, no-slip, co-moving (with the equilibrium flow V0) walls at z=±Lz/2. The simulations presented here use a resolution of Nx×Nz=512×2048, 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 ϕ and ψ at every nonzero kx with Gaussians in z that have randomly assigned, kx-dependent complex phases and amplitudes that decrease with kx as a power law. Thus, at t =0, the streamfunction and flux function are

(7)

and

(8)

where the kx = 0 components are the unperturbed equilibrium profiles, Aϕ and Aψ 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, Δϕ(kx) and Δψ(kx) are uniformly distributed pseudo-random numbers in [0,2π). For the results presented here, we use σ = 2, a = −1, and Aϕ=Aψ=5×104, which allows for a clearly defined regime of linear growth before nonlinear interactions become important. For MA=5, setting Aψ=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ϕ.

Previous work44 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 kx, 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 Δϕ and Δψ. In practice, this is done by selecting different seeds for our pseudo-random number generator (we use numpy.random.RandomState45,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 Δϕ and Δψ 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.

For U=tanh(z) and Bx = 1, unstable modes, solutions to Eqs. (5) and (6) with positive growth rates γj=Im[ω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 kx 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 MA8.

FIG. 1.

Growth rate γ for the fastest-growing mode at each kx. Each curve corresponds to a different Alfvén Mach number MA. While stronger magnetic fields (lower MA) provide a stabilizing influence, γ varies little except when MA4. Horizontal dashed lines indicate the kx present in our nonlinear simulations.

FIG. 1.

Growth rate γ for the fastest-growing mode at each kx. Each curve corresponds to a different Alfvén Mach number MA. While stronger magnetic fields (lower MA) provide a stabilizing influence, γ varies little except when MA4. Horizontal dashed lines indicate the kx present in our nonlinear simulations.

Close modal

Taking the complex conjugate of Eqs. (5) and (6) shows that, as in the hydrodynamic5 and gyrokinetic cases,32 for every eigenvalue ωj and eigenmode (ϕj,ψj) that is a solution of Eqs. (5) and (6), the complex conjugate ωj* and (ϕj*,ψj*) is a solution as well (recent work has explored the connection between this conjugate symmetry and parity-time symmetry47). Following Refs. 22 and 32, when describing the eigenmodes of this system, we label the most unstable mode at each kx as j =1 and its conjugate stable mode as j =2. Real-space contours corresponding to ϕj and ψj for 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 subject48–51), the current density of the mode is more localized about the center of the layer.

FIG. 2.

The initial U=tanh(z) unstable equilibrium flow (top left) and uniform field (bottom left) are shown alongside contours (with arbitrary units) of the streamfunction (center column) and flux function (right column) for the unstable (top) and stable (bottom) modes at kx=0.4 for MA=40.

FIG. 2.

The initial U=tanh(z) unstable equilibrium flow (top left) and uniform field (bottom left) are shown alongside contours (with arbitrary units) of the streamfunction (center column) and flux function (right column) for the unstable (top) and stable (bottom) modes at kx=0.4 for MA=40.

Close modal

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 kx0 perturbations are allowed to evolve, then this exponential growth does not conserve energy because the energy injected into f1 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 kx0 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

(9)

and

(10)

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

(11)

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 f1, the exponential decay of the conjugate stable mode f2 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 f1 and f2 at kx=0.4 for both MA=25 and MA=2. For unstable modes, both τu and τb transport momentum down the gradient, so that they both contribute to a transfer of energy from U(z) to f1. Likewise, for stable modes, both stresses yield counter-gradient momentum transport, transferring energy from f2 to U(z). The transport of the two modes is symmetric in the sense that τu[ϕ2]=τu[ϕ1] and τb[ψ2]=τb[ψ1]. As MA is decreased, corresponding to a stronger equilibrium field, the relative amplitudes of τu and τb change, with |τb| exceeding |τu| for only the strongest equilibrium fields, starting around MA2.5.

FIG. 3.

xz components of the Reynolds stress τu (blue) and Maxwell stress τb (orange) for the unstable mode f1 (solid lines) and stable mode f2 (dashed lines) at kx=0.4 for MA=25 (left) and MA=2 (right). For MA=25, τb is rescaled by a factor of 5 to improve visibility. Modes are normalized to have unit total energy.

FIG. 3.

xz components of the Reynolds stress τu (blue) and Maxwell stress τb (orange) for the unstable mode f1 (solid lines) and stable mode f2 (dashed lines) at kx=0.4 for MA=25 (left) and MA=2 (right). For MA=25, τb is rescaled by a factor of 5 to improve visibility. Modes are normalized to have unit total energy.

Close modal

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.

FIG. 4.

Snapshots at three different times for a hydrodynamic simulation (top row) and an MHD simulation with MA = 60 and Rm=250 (middle and bottom rows). In the top and middle rows, the color shows vorticity 2ϕ and black lines show contours of the streamfunction ϕ (streamlines of the flow). In the bottom row, color shows current density 2ψ and black lines show contours of the flux function ψ (field lines). Current sheets form along the braids between vortices at early times and at the edges of the vortex at later times. For vorticity plots, the colorbar has been rescaled with 2ϕ=1 as the maximum and −1 as the minimum with white as 0 to demonstrate that, in the hydrodynamic case, all of the vorticity is into the page, consistent with the conservative, advection-diffusion nature of the vorticity equation, and that this is broken in MHD, with negative vorticity appearing at late times even at this high MA and low Rm.

FIG. 4.

Snapshots at three different times for a hydrodynamic simulation (top row) and an MHD simulation with MA = 60 and Rm=250 (middle and bottom rows). In the top and middle rows, the color shows vorticity 2ϕ and black lines show contours of the streamfunction ϕ (streamlines of the flow). In the bottom row, color shows current density 2ψ and black lines show contours of the flux function ψ (field lines). Current sheets form along the braids between vortices at early times and at the edges of the vortex at later times. For vorticity plots, the colorbar has been rescaled with 2ϕ=1 as the maximum and −1 as the minimum with white as 0 to demonstrate that, in the hydrodynamic case, all of the vorticity is into the page, consistent with the conservative, advection-diffusion nature of the vorticity equation, and that this is broken in MHD, with negative vorticity appearing at late times even at this high MA and low Rm.

Close modal

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 Rm1.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=kxKEkx and ME=kxMEkx. 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 Δϕ(kx) of the initial flow perturbation, while the saturation time for kx=0.2 does depend on Δϕ(kx). Varying Δψ had no significant effect on the simulations. The simulations shown here correspond to a choice of Δϕ(kx) with a relatively long merging time.

FIG. 5.

Energy components vs time for the simulations shown in Fig. 4 with MA=60 and Rm=250 (top row), and a similar simulation with MA=40 and Rm=500 (bottom row), both compared to the same hydrodynamic simulation. The left column plots energy on a log scale and the right on a linear scale. Solid, dashed, and dotted-dashed curves correspond to kinetic, magnetic, and total energy, respectively. Dotted curves correspond to the hydrodynamic case. Colors indicate different kx contributions, with total energy summed over all kx shown in brown, the contribution from all nonzero kx in red, and contribution from small scales, given by summing over kx1, shown in black. The initial saturation and merger stages show no obvious changes from the hydrodynamic case, but local minima in KEkx=0 (blue) become less pronounced with stronger fields. Small-scale fluctuations become enhanced with increased field strength.

FIG. 5.

Energy components vs time for the simulations shown in Fig. 4 with MA=60 and Rm=250 (top row), and a similar simulation with MA=40 and Rm=500 (bottom row), both compared to the same hydrodynamic simulation. The left column plots energy on a log scale and the right on a linear scale. Solid, dashed, and dotted-dashed curves correspond to kinetic, magnetic, and total energy, respectively. Dotted curves correspond to the hydrodynamic case. Colors indicate different kx contributions, with total energy summed over all kx shown in brown, the contribution from all nonzero kx in red, and contribution from small scales, given by summing over kx1, shown in black. The initial saturation and merger stages show no obvious changes from the hydrodynamic case, but local minima in KEkx=0 (blue) become less pronounced with stronger fields. Small-scale fluctuations become enhanced with increased field strength.

Close modal

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 t55 and t85. 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.

FIG. 6.

Kinetic energy over time in the mean flow, KEkx=0, for a variety of MA. Each simulation has Rm=500 and the same initial conditions. As field strength increases (MA decreases), KEkx=0 decreases more rapidly, equivalent to a faster layer broadening rate.

FIG. 6.

Kinetic energy over time in the mean flow, KEkx=0, for a variety of MA. Each simulation has Rm=500 and the same initial conditions. As field strength increases (MA decreases), KEkx=0 decreases more rapidly, equivalent to a faster layer broadening rate.

Close modal

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 kx = 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 Ux and Bxx in place of U and Bx, 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 fj(ϕj(z),ψj(z)),ωj, and γj 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, γ2=γ1 still holds. Figure 7 shows how growth rates γj 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 kx = 1. The linear growth regime for KEkx=0.4 is seen in Fig. 5 to end at about the same time that γ1(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 z of these modes and their corresponding Reynolds and Maxwell stresses are largely the same as f1 and f2, except that they broaden with the shear layer.

FIG. 7.

Growth rates for the eigenmodes of Ux and Bxx taken from a simulation with MA=40 and Rm=500. Different colors correspond to different wavenumbers. For each wavenumber and at each time, only the most unstable growth rate is plotted. Even when perturbations are small and the energy removed from the mean flow appears negligible at early times, the growth rates at higher kx decline rapidly.

FIG. 7.

Growth rates for the eigenmodes of Ux and Bxx taken from a simulation with MA=40 and Rm=500. Different colors correspond to different wavenumbers. For each wavenumber and at each time, only the most unstable growth rate is plotted. Even when perturbations are small and the energy removed from the mean flow appears negligible at early times, the growth rates at higher kx decline rapidly.

Close modal

In freely evolving shear layers, the mean flow has been shown to depart from a simple, broadened tanh profile in both MHD simulations7,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 fj and ωj, minor features in Ux and Bxx 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 Ux. At times, changes in Bxx 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 kx0.4 unstable modes that emerge around t60 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 fj of the broadened system correspond more directly to energy transfer to/from the mean than the eigenmodes fj 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.

At every time t and wavenumber kx, the Fourier-transformed system state f̂(kx,z,t)(ϕ̂(kx,z,t),ψ̂(kx,z,t)) can be expressed as

(12)

provided the eigenmodes of the equilibrium, {fj(kx,z)}, form a complete basis. The complex-valued βj(kx,t) is the amplitude of mode fj(kx,z) at time t and can be understood as the coefficient of the state vector f̂ expressed in this basis. The eigenmodes of the horizontally averaged system, {fj(kx,z)}, can also be used as a basis, and the amplitudes of these modes will be denoted as βj(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, βj is obtained by expanding f̂ in terms of the eigenmodes {fj} of the instantaneous mean flow and field. Leading into saturation, both sets of amplitudes, {βj} and {βj}, 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 |βj| and |βj| curves demonstrates some of the differences between the two eigenmode bases. Two notable features indicate that fj fails to capture the dynamics as precisely as the more relevant fj.

FIG. 8.

Mode amplitudes of the most unstable (blue) and conjugate stable (orange) modes for a simulation with MA = 40 and Rm=500. Dashed lines are mode amplitudes using a basis of eigenmodes of perturbations about the mean flow and field, solid lines are using eigenmodes of perturbations about the initial equilibrium. For kx=0.4, dashed lines terminate when that wavenumber first stabilizes.

FIG. 8.

Mode amplitudes of the most unstable (blue) and conjugate stable (orange) modes for a simulation with MA = 40 and Rm=500. Dashed lines are mode amplitudes using a basis of eigenmodes of perturbations about the mean flow and field, solid lines are using eigenmodes of perturbations about the initial equilibrium. For kx=0.4, dashed lines terminate when that wavenumber first stabilizes.

Close modal

First, there are periods where |β2| grows at the same rate as |β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 β2(t)=β2(0)e|γ|t+β12(t), and the unstable mode amplitude as β1(t)=β1(0)e|γ|t [see, for example, Eqs. (22)–(25) in Ref. 21]. If a mode that differs slightly from f2 is used to calculate β2(t) from simulation data, then an error of the form ϵβ1 will almost always be introduced. This will cause the apparent stable mode amplitude to briefly evolve as β2(t)ϵβ1(t) before the β2(t)β12 term becomes dominant, provided that |ϵ||β1(0)|. For even small differences in f2 (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 |β2| 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 β2/t:21,22,27 when one interaction overtakes another and becomes dominant, if the two have approximately opposite complex phases, then β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 U and Bx from U=tanh(z) and Bx = 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 βj capture subtler details leading into saturation than the equilibrium amplitudes βj, and the corresponding eigenmodes fj capture the layer broadening, they become unreliable once γ10 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̂ 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 |β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 methods66) when fluctuations at kx=0.2 and 0.4, respectively, are approximated in this manner for a simulation with MA = 40 and Rm=500. Black and gray curves compare approximations using βj and fj, while blue and orange curves compare approximations using βj and fj, where f1 and f2 are defined as the mode with the largest growth rate and its conjugate. The truncation errors are calculated separately for ϕ̂ and ψ̂ according to

(13)

and

(14)

where ||ϕ̂||KE2 and ||ψ̂||ME2 are the kinetic and magnetic energy of a given ϕ̂ and ψ̂. Unsurprisingly, the fj modes describe flow fluctuations better than the fj 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 ϵT[ϕ̂] 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 ϵT[ϕ̂] 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.

FIG. 9.

Error over time for approximations of the flow (left) and field (right) at kx=0.2 using a truncated eigenmode expansion for a simulation with MA=40 and Rm=500. Black and gray curves use βj and fj, while other curves use βj and fj. The horizontal dashed line indicates an error of 1. For a given approximation, ϵT1 implies that the difference between the exact state and the approximation is at least as large, energetically speaking, as the state itself, and so the approximation is unreliable. The importance of stable modes can be seen by noting that flow approximations are improved when stable modes are included in both sets of models, with models using βj and fj performing better than those that use βj and fj until quasilinear flattening stabilizes the flow, at which point approximations using fj become unreliable. Each approximation describes magnetic fluctuations poorly, particularly after the linear regime.

FIG. 9.

Error over time for approximations of the flow (left) and field (right) at kx=0.2 using a truncated eigenmode expansion for a simulation with MA=40 and Rm=500. Black and gray curves use βj and fj, while other curves use βj and fj. The horizontal dashed line indicates an error of 1. For a given approximation, ϵT1 implies that the difference between the exact state and the approximation is at least as large, energetically speaking, as the state itself, and so the approximation is unreliable. The importance of stable modes can be seen by noting that flow approximations are improved when stable modes are included in both sets of models, with models using βj and fj performing better than those that use βj and fj until quasilinear flattening stabilizes the flow, at which point approximations using fj become unreliable. Each approximation describes magnetic fluctuations poorly, particularly after the linear regime.

Close modal
FIG. 10.

Identical to Fig. 9, but the kx=0.4 fluctuations are approximated. Unlike the kx=0.2 fluctuations, here the fj modes describe ϕ̂ well for significantly longer than the fj modes. This is to be expected as γj approaches 0 much sooner at this kx.

FIG. 10.

Identical to Fig. 9, but the kx=0.4 fluctuations are approximated. Unlike the kx=0.2 fluctuations, here the fj modes describe ϕ̂ well for significantly longer than the fj modes. This is to be expected as γj approaches 0 much sooner at this kx.

Close modal
FIG. 11.

Identical to Fig. 9, but for MA=7.5. Here, including stable modes does not improve approximations as noticeably, consistent with their diminished role in the dynamics at lower MA.

FIG. 11.

Identical to Fig. 9, but for MA=7.5. Here, including stable modes does not improve approximations as noticeably, consistent with their diminished role in the dynamics at lower MA.

Close modal

Figures 9–11 assess the ability of the equilibrium modes fj to describe the state f̂ 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, τu(kx,z=0), over time for three different MA. The quantity |β2|2|β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 kx but not MA. The equilibrium modes fj 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.

FIG. 12.

Dashed curves show the contribution to Reynolds stress τu made by fluctuations at kx=0.2 (left) and 0.4 (right) evaluated at z =0 vs time for three of the simulations shown in Fig. 6, each with a different value of MA. Solid curves show the quantity |β2|2|β1|2 over time for the same simulations. The agreement between the two quantities motivates the use of these two modes to characterize the system.

FIG. 12.

Dashed curves show the contribution to Reynolds stress τu made by fluctuations at kx=0.2 (left) and 0.4 (right) evaluated at z =0 vs time for three of the simulations shown in Fig. 6, each with a different value of MA. Solid curves show the quantity |β2|2|β1|2 over time for the same simulations. The agreement between the two quantities motivates the use of these two modes to characterize the system.

Close modal

From Fig. 12 alone, it is unclear whether the general decrease in ||β2|2|β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 |β1| and |β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.

FIG. 13.

Mode amplitudes |β2| (orange) and |β1| (blue) at kx=0.2 (left) and 0.4 (right) for MA=40 (solid), and 15 (dashed lines). As MA decreases, both mode amplitudes are reduced while also trending closer toward equipartition.

FIG. 13.

Mode amplitudes |β2| (orange) and |β1| (blue) at kx=0.2 (left) and 0.4 (right) for MA=40 (solid), and 15 (dashed lines). As MA decreases, both mode amplitudes are reduced while also trending closer toward equipartition.

Close modal

We have established that stable modes are directly responsible for reversals in the mid-layer Reynolds stress of low-kx fluctuations. To demonstrate that τu 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 τb. 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 τu and Maxwell stress τb, which are themselves separated into contributions from different kx. At all times, |τu| exceeds |τ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 τu and not τb, they are specifically due to τu 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.

FIG. 14.

Different contributions to momentum transport at four different times for the same MHD simulation as in Fig. 4. Solid lines show Reynolds stresses τu and dashed lines show Maxwell stresses τb (rescaled by a factor of 100 for improved visibility). Black lines are the total Reynolds and Maxwell stresses, while different colors indicate contributions from different kx. While τb is generally dominated by larger kx and is always negative, τu is dominated by smaller kx and changes sign over time.

FIG. 14.

Different contributions to momentum transport at four different times for the same MHD simulation as in Fig. 4. Solid lines show Reynolds stresses τu and dashed lines show Maxwell stresses τb (rescaled by a factor of 100 for improved visibility). Black lines are the total Reynolds and Maxwell stresses, while different colors indicate contributions from different kx. While τb is generally dominated by larger kx and is always negative, τu is dominated by smaller kx and changes sign over time.

Close modal

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 kx of magnetic energy. While the low-kx Reynolds stresses resemble the Reynolds stress of the eigenmodes τu[ϕj]—albeit broadened, consistent with the broadening of the flow profile and eigenmodes discussed in Sec. IV A—the low-kx Maxwell stresses bear less of a resemblance to τb[ψj], or to τb[ψj], particularly as time goes on. Furthermore, τb was never observed to change sign for the simulations considered here. When counter-gradient momentum transport occurs, it is always due to τu changing sign. This is consistent with the truncation error for the field, ϵT[ψ̂], generally exceeding the truncation error for the flow, ϵT[ϕ̂]: if stable and unstable modes alone captured the Maxwell stress as well as they capture the Reynolds stress, then τb at low kx would reverse sign whenever τu does. Instead, τb 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,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 τu. This is seen in Fig. 15, where τu and τb are shown for different MA, focusing on the first two counter-gradient transport phases. As field strength increases, |τ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-kx 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, τu 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 |τb(kx)| for large kx—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 MA5. For lower MA, the counter-gradient τu(kx=0.4) at this time does become partially reduced and then fully removed as MA decreases, and τb 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 τu amplitudes is similar to what happens during the down-gradient transport phases at t =40 and 85: |τu| becomes weaker as field strength increases at later times, but does not change with MA at earlier times.

FIG. 15.

The same breakdown of stresses as in Fig. 14, but for four simulations with Rm=500 and different MA. The two rows correspond to the first two instances of counter-gradient momentum transport (see, e.g., Fig. 6). Note that τb is rescaled by a factor of 10. As MA decreases, τb becomes more dominant, reducing net counter-gradient transport. At earlier times |τu| varies little with MA, but at later times it decreases with stronger fields.

FIG. 15.

The same breakdown of stresses as in Fig. 14, but for four simulations with Rm=500 and different MA. The two rows correspond to the first two instances of counter-gradient momentum transport (see, e.g., Fig. 6). Note that τb is rescaled by a factor of 10. As MA decreases, τb becomes more dominant, reducing net counter-gradient transport. At earlier times |τu| varies little with MA, but at later times it decreases with stronger fields.

Close modal

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.

The enhanced generation of small-scale fluctuations as magnetic field strength increases has been noted in terms of the distribution in kx 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 viscosity40and resistivity7 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.

FIG. 16.

Viscous (solid) and resistive (dashed) dissipation over time for simulations with the same initial conditions as those in Fig. 4, with different colors corresponding to different MA and Rm, and black corresponding to the hydrodynamic case. Even a weak magnetic field significantly increases total dissipation relative to the hydrodynamic case. Resistive dissipation increases with Rm, and peaks at a time that increases with Rm. Viscous dissipation is similar to the hydrodynamic case except when the amplified field increases beyond a threshold that depends on both MA and Rm.

FIG. 16.

Viscous (solid) and resistive (dashed) dissipation over time for simulations with the same initial conditions as those in Fig. 4, with different colors corresponding to different MA and Rm, and black corresponding to the hydrodynamic case. Even a weak magnetic field significantly increases total dissipation relative to the hydrodynamic case. Resistive dissipation increases with Rm, and peaks at a time that increases with Rm. Viscous dissipation is similar to the hydrodynamic case except when the amplified field increases beyond a threshold that depends on both MA and Rm.

Close modal

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 fixed21,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 Bx = 1, denoted by fj, and those obtained by linearizing about the instantaneous mean flow Ux and field Bxx at each time step, denoted by fj. The fj 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 fj modes do not correspond to the evolved system as directly as the fj 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 kx. Thus, the same eigenmode decompositions used in systems with a fixed28,29 or otherwise quasi-stationary32 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 MA = 100), with dissipation increasing as Rm increases or MA 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.

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 research68 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 work69 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 interiors70 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 hydrodynamic71,72 and MHD73 cases (similar to channel modes of the magnetorotational instability74–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 case78 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 work67 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 theories15 that were informed by studies of dominant nonlinear interactions30 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 

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.

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

1.
M.
Faganello
and
F.
Califano
, “
Magnetized Kelvin–Helmholtz instability: Theory and simulations in the Earth's magnetosphere context
,”
J. Plasma Phys.
83
,
535830601
(
2017
).
2.
F.
Rieger
, “
An introduction to particle acceleration in shearing flows
,”
Galaxies
7
,
78
(
2019
).
3.
K.
Kwak
,
R. L.
Shelton
, and
D. B.
Henley
, “
Si IV column densities predicted from non-equilibrium ionization simulations of turbulent mixing layers and high-velocity clouds
,”
Astrophys. J.
812
,
111
(
2015
).
4.
S.
Chandrasekhar
,
Hydrodynamic and Hydromagnetic Stability
(
Oxford University Press
,
1961
).
5.
P. G.
Drazin
and
W.
Reid
,
Hydrodynamic Stability
(
Cambridge University Press
,
1981
).
6.
D.
P. G
,
Introduction to Hydrodynamic Stability
(
Cambridge University Press
,
2002
).
7.
M. L.
Palotti
,
F.
Heitsch
,
E. G.
Zweibel
, and
Y.
Huang
, “
Evolution of unmagnetized and magnetized shear layers
,”
Astrophys. J.
678
,
234
244
(
2008
).
8.
R.
Keppens
,
G.
Tóth
,
R. H. J.
Westermann
, and
J. P.
Goedbloed
, “
Growth and saturation of the Kelvin–Helmholtz instability with parallel and antiparallel magnetic fields
,”
J. Plasma Phys.
61
,
1
19
(
1999
).
9.
T. W.
Jones
,
J. B.
Gaalaas
,
D.
Ryu
, and
A.
Frank
, “
The MHD Kelvin–Helmholtz instability. II. The roles of weak and oblique fields in planar flows
,”
Astrophys. J.
482
,
230
(
1997
).
10.
J.
Mak
,
S. D.
Griffiths
, and
D. W.
Hughes
, “
Vortex disruption by magnetohydrodynamic feedback
,”
Phys. Rev. Fluids
2
,
113701
(
2017
).
11.
H.
Jeong
,
D.
Ryu
,
T. W.
Jones
, and
A.
Frank
, “
The magnetohydrodynamic Kelvin–Helmholtz instability. III. The role of sheared magnetic field in planar flows
,”
Astrophys. J.
529
,
536
(
2000
).
12.
S. A.
Balbus
and
J. F.
Hawley
, “
Instability, turbulence, and enhanced transport in accretion disks
,”
Rev. Mod. Phys.
70
,
1
53
(
1998
).
13.
B.
Paxton
,
M.
Cantiello
,
P.
Arras
,
L.
Bildsten
,
E. F.
Brown
,
A.
Dotter
,
C.
Mankovich
,
M. H.
Montgomery
,
D.
Stello
,
F. X.
Timmes
 et al., “
Modules for experiments in stellar astrophysics (MESA): Planets, oscillations, rotation, and massive stars
,”
Astrophys. J. Suppl. Ser.
208
,
4
(
2013
).
14.
A.
Heger
,
N.
Langer
, and
S. E.
Woosley
, “
Presupernova evolution of rotating massive stars. I. Numerical method and evolution of the internal stellar structure
,”
Astrophys. J.
528
,
368
396
(
2000
).
15.
P. W.
Terry
,
B. J.
Faber
,
C. C.
Hegna
,
V. V.
Mirnov
,
M. J.
Pueschel
, and
G. G.
Whelan
, “
Saturation scalings of toroidal ion temperature gradient turbulence
,”
Phys. Plasmas
25
,
012308
(
2018
).
16.
C. C.
Hegna
,
P. W.
Terry
, and
B. J.
Faber
, “
Theory of ITG turbulent saturation in stellarators: Identifying mechanisms to reduce turbulent transport
,”
Phys. Plasmas
25
,
022511
(
2018
).
17.
G. G.
Whelan
,
M. J.
Pueschel
, and
P. W.
Terry
, “
Nonlinear electromagnetic stabilization of plasma microturbulence
,”
Phys. Rev. Lett.
120
,
175002
(
2018
).
18.
G. G.
Whelan
,
M. J.
Pueschel
,
P. W.
Terry
,
J.
Citrin
,
I. J.
McKinney
,
W.
Guttenfelder
, and
H.
Doerk
, “
Saturation and nonlinear electromagnetic stabilization of ITG turbulence
,”
Phys. Plasmas
26
,
082302
(
2019
).
19.
F.
Ebrahimi
,
S. C.
Prager
, and
D. D.
Schnack
, “
Saturation of magnetorotational instability through magnetic field generation
,”
Astrophys. J.
698
,
233
241
(
2009
).
20.
J.-H.
Xie
,
K.
Julien
, and
E.
Knobloch
, “
Subcritical saturation of the magnetorotational instability through mean magnetic field generation
,”
Mon. Not. R. Astron. Soc.
474
,
3451
3465
(
2018
).
21.
P. W.
Terry
,
D. A.
Baver
, and
S.
Gupta
, “
Role of stable eigenmodes in saturated local plasma turbulence
,”
Phys. Plasmas
13
,
022307
(
2006
).
22.
A. E.
Fraser
,
P. W.
Terry
,
E. G.
Zweibel
, and
M. J.
Pueschel
, “
Coupling of damped and growing modes in unstable shear flow
,”
Phys. Plasmas
24
,
062304
(
2017
).
23.
T. M.
Qian
and
M. E.
Mauel
, “
Observation of weakly damped modes using high resolution measurement of turbulence in a dipole confined plasma
,”
Phys. Plasmas
27
,
014501
(
2020
).
24.
C.
Ho
and
P.
Huerre
, “
Perturbed free shear layers
,”
Annu. Rev. Fluid Mech.
16
,
365
422
(
1984
).
25.
N. C.
Hurst
,
J. R.
Danielson
,
D. H. E.
Dubin
, and
C. M.
Surko
, “
Instability of an electron-plasma shear layer in an externally imposed strain flow
,”
Phys. Plasmas
27
,
042101
(
2020
).
26.
K.
Takamure
,
Y.
Ito
,
Y.
Sakai
,
K.
Iwano
, and
T.
Hayase
, “
Momentum transport process in the quasi self-similar region of free shear mixing layer
,”
Phys. Fluids
30
,
015109
(
2018
).
27.
K. D.
Makwana
,
P. W.
Terry
,
J.-H. H.
Kim
, and
D. R.
Hatch
, “
Damped eigenmode saturation in plasma fluid turbulence
,”
Phys. Plasmas
18
,
012302
(
2011
).
28.
D. R.
Hatch
,
P. W.
Terry
,
F.
Jenko
,
F.
Merz
,
M. J.
Pueschel
,
W. M.
Nevins
, and
E.
Wang
, “
Role of subdominant stable modes in plasma microturbulence
,”
Phys. Plasmas
18
,
055706
(
2011
).
29.
P. W.
Terry
,
K. D.
Makwana
,
M. J.
Pueschel
,
D. R.
Hatch
,
F.
Jenko
, and
F.
Merz
, “
Mode-space energy distribution in instability-driven plasma turbulence
,”
Phys. Plasmas
21
,
122303
(
2014
).
30.
K. D.
Makwana
,
P. W.
Terry
,
M. J.
Pueschel
, and
D. R.
Hatch
, “
Subdominant modes in zonal-flow-regulated turbulence
,”
Phys. Rev. Lett.
112
,
095002
(
2014
).
31.
A. E.
Fraser
, “
Role of stable eigenmodes in shear-flow instability saturation and turbulence
,” Ph.D. thesis (University of Wisconsin-Madison,
2020
): available at https://search.proquest.com/dissertations-theses/role-stable-eigenmodes-shear-flow-instability/docview/2445291638/se-2?accountid=14523.
32.
A. E.
Fraser
,
M. J.
Pueschel
,
P. W.
Terry
, and
E. G.
Zweibel
, “
Role of stable modes in driven shear-flow turbulence
,”
Phys. Plasmas
25
,
122303
(
2018
).
33.
M. J.
Pueschel
,
F.
Jenko
,
D.
Told
, and
J.
Büchner
, “
Gyrokinetic simulations of magnetic reconnection
,”
Phys. Plasmas
18
,
112102
(
2011
).
34.
M. J.
Pueschel
,
D.
Told
,
P. W.
Terry
,
F.
Jenko
,
E. G.
Zweibel
,
V.
Zhdankin
, and
H.
Lesch
, “
Magnetic reconnection turbulence in strong guide fields: Basic properties and application to coronal heating
,”
Astrophys. J., Suppl. Ser.
213
,
30
(
2014
).
35.
N.
Platt
,
L.
Sirovich
, and
N.
Fitzmaurice
, “
An investigation of chaotic Kolmogorov flows
,”
Phys. Fluids A
3
,
681
696
(
1991
).
36.
S.
Musacchio
and
G.
Boffetta
, “
Turbulent channel without boundaries: The periodic Kolmogorov flow
,”
Phys. Rev. E
89
,
023004
(
2014
).
37.
D.
Lucas
and
R.
Kerswell
, “
Spatiotemporal dynamics in two-dimensional Kolmogorov flow over large domains
,”
J. Fluid Mech.
750
,
518
554
(
2014
).
38.
K. M.
Case
, “
Stability of inviscid plane Couette flow
,”
Phys. Fluids
3
,
143
148
(
1960
).
39.
K. J.
Burns
,
G. M.
Vasil
,
J. S.
Oishi
,
D.
Lecoanet
, and
B. P.
Brown
, “
Dedalus: A flexible framework for numerical simulations with spectral methods
,”
Phys. Rev. Res.
2
,
23068
(
2020
).
40.
D.
Lecoanet
,
M.
McCourt
,
E.
Quataert
,
K. J.
Burns
,
G. M.
Vasil
,
J. S.
Oishi
,
B. P.
Brown
,
J. M.
Stone
, and
R. M.
O'Leary
, “
A validated non-linear Kelvin–Helmholtz benchmark for numerical hydrodynamics
,”
Mon. Not. R. Astron. Soc.
455
,
4274
4288
(
2016
).
41.
D. R.
Hatch
,
F.
Jenko
,
A. B.
Navarro
, and
V.
Bratanov
, “
Transition between saturation regimes of gyrokinetic turbulence
,”
Phys. Rev. Lett.
111
,
175001
(
2013
).
42.
D.
Biskamp
,
Magnetohydrodynamic Turbulence
(
Cambridge University Press
,
2003
).
43.
D. R.
Hatch
,
F.
Jenko
,
A. B.
Navarro
,
V.
Bratanov
,
P. W.
Terry
, and
M. J.
Pueschel
, “
Linear signatures in nonlinear gyrokinetics: Interpreting turbulence with pseudospectra
,”
New J. Phys.
18
,
075018
(
2016
).
44.
W.
Dong
,
E. W.
Tedford
,
M.
Rahmani
, and
G. A.
Lawrence
, “
Sensitivity of vortex pairing and mixing to initial perturbations in stratified shear flows
,”
Phys. Rev. Fluids
4
,
063902
(
2019
).
45.
C. R.
Harris
,
K. J.
Millman
,
S. J.
van der Walt
,
R.
Gommers
,
P.
Virtanen
,
D.
Cournapeau
,
E.
Wieser
,
J.
Taylor
,
S.
Berg
,
N. J.
Smith
,
R.
Kern
,
M.
Picus
,
S.
Hoyer
,
M. H.
van Kerkwijk
,
M.
Brett
,
A.
Haldane
,
J.
Fernández del Río
,
M.
Wiebe
,
P.
Peterson
,
P.
Gérard-Marchant
,
K.
Sheppard
,
T.
Reddy
,
W.
Weckesser
,
H.
Abbasi
,
C.
Gohlke
, and
T. E.
Oliphant
, “
Array programming with NumPy
,”
Nature
585
,
357
362
(
2020
).
46.
P.
Virtanen
,
R.
Gommers
,
T. E.
Oliphant
,
M.
Haberland
,
T.
Reddy
,
D.
Cournapeau
,
E.
Burovski
,
P.
Peterson
,
W.
Weckesser
,
J.
Bright
,
S. J.
van der Walt
,
M.
Brett
,
J.
Wilson
,
K. J.
Millman
,
N.
Mayorov
,
A. R. J.
Nelson
,
E.
Jones
,
R.
Kern
,
E.
Larson
,
C. J.
Carey
,
İ.
Polat
,
Y.
Feng
,
E. W.
Moore
,
J.
VanderPlas
,
D.
Laxalde
,
J.
Perktold
,
R.
Cimrman
,
I.
Henriksen
,
E. A.
Quintero
,
C. R.
Harris
,
A. M.
Archibald
,
A. H.
Ribeiro
,
F.
Pedregosa
, and
P.
van Mulbregt
, and
SciPy 1.0 Contributors
, “
SciPy 1.0: Fundamental algorithms for scientific computing in Python
,”
Nat. Methods
17
,
261
272
(
2020
).
47.
H.
Qin
,
Y.
Fu
,
A. S.
Glasser
, and
A.
Yahalom
, “
Spontaneous and explicit parity-time-symmetry breaking in drift wave instabilities
,” arXiv:2010.09620 (
2020
).
48.
P. G.
Baines
and
H.
Mitsudera
, “
On the mechanism of shear flow instabilities
,”
J. Fluid Mech.
276
,
327
342
(
1994
).
49.
E.
Heifetz
,
J.
Mak
,
J.
Nycander
, and
O. M.
Umurhan
, “
Interacting vorticity waves as an instability mechanism for magnetohydrodynamic shear instabilities
,”
J. Fluid Mech.
767
,
199
225
(
2015
).
50.
E.
Heifetz
and
A.
Guha
, “
Normal form of synchronization and resonance between vorticity waves in shear flow instability
,”
Phys. Rev. E
100
,
043105
(
2019
).
51.
J. R.
Carpenter
,
E. W.
Tedford
,
E.
Heifetz
, and
G. A.
Lawrence
, “
Instability in stratified shear flow: Review of a physical interpretation based on interacting waves
,”
Appl. Mech. Rev.
64
,
060801
(
2011
).
52.
S.
Kida
and
M.
Takaoka
, “
Vortex reconnection
,”
Annu. Rev. Fluid Mech.
26
,
169
177
(
1994
).
53.
E. G.
Zweibel
, “
Fast reconnection of weak magnetic fields
,”
Phys. Plasmas
5
,
247
251
(
1998
).
54.
Y. B.
Zel'Dovich
,
A. A.
Ruzmaikin
,
S. A.
Molchanov
, and
D. D.
Sokoloff
, “
Kinematic dynamo problem in a linear velocity field
,”
J. Fluid Mech.
144
,
1
11
(
1984
).
55.
D. G.
Dritschel
,
P. H.
Diamond
, and
S. M.
Tobias
, “
Circulation conservation and vortex breakup in magnetohydrodynamics at low magnetic Prandtl number
,”
J. Fluid Mech.
857
,
38
60
(
2018
).
56.
W.
N. O
, “
The expulsion of magnetic flux by eddies
,”
Proc. R. Soc. London, Ser. A
293
,
310
328
(
1966
).
57.
M.
Gaster
,
E.
Kit
, and
I.
Wygnanski
, “
Large-scale structures in a forced turbulent mixing layer
,”
J. Fluid Mech.
150
,
23
39
(
1985
).
58.
X.
Wu
and
X.
Zhuang
, “
Nonlinear dynamics of large-scale coherent structures in turbulent free shear layers
,”
J. Fluid Mech.
787
,
396
439
(
2016
).
59.
T.
Tatsuno
and
W.
Dorland
, “
Magneto-flow instability in symmetric field profiles
,”
Phys. Plasmas
13
,
092107
(
2006
).
60.
D.
Hughes
and
S.
Tobias
, “
On the instability of magnetohydrodynamic shear flows
,”
Proc. R. Soc. London, Ser. A
457
,
1365
1384
(
2001
).
61.
K. J.
Burns
, “
Flexible spectral algorithms for simulating astrophysical and geophysical flows
,” Ph.D. thesis (
Massachusetts Institute of Technology
,
2018
).
62.
P. J.
Schmid
, “
Nonmodal stability theory
,”
Annu. Rev. Fluid Mech.
39
,
129
162
(
2007
).
63.
S. C.
Reddy
,
P. J.
Schmid
, and
D. S.
Henningson
, “
Pseudospectra of the Orr-Sommerfeld operator
,”
SIAM J. Appl. Math.
53
,
15
47
(
1993
).
64.
J.
Squire
and
A.
Bhattacharjee
, “
Nonmodal growth of the magnetorotational instability
,”
Phys. Rev. Lett.
113
,
025006
(
2014
).
65.
J.
Squire
and
A.
Bhattacharjee
, “
Magnetorotational instability: Nonmodal Growth and the relationship of global modes to the shearing box
,”
Astrophys. J.
797
,
67
(
2014
).
66.
J.
Boyd
,
Chebyshev and Fourier Spectral Methods: Second Revised Edition
, Dover Books on Mathematics (
Dover Publications
,
2001
).
67.
J. B.
Parker
and
N. C.
Constantinou
, “
Magnetic eddy viscosity of mean shear flows in two-dimensional magnetohydrodynamics
,”
Phys. Rev. Fluids
4
,
083701
(
2019
).
68.
S.
Veilleux
,
R.
Maiolino
,
A. D.
Bolatto
, and
S.
Aalto
, “
Cool outflows in galaxies and their implications
,”
Astron. Astrophys. Rev.
28
,
2
(
2020
).
69.
A.
Grønnow
,
T.
Tepper-García
, and
J.
Bland-Hawthorn
, “
Magnetic fields in the galactic halo restrict fountain-driven recycling and accretion
,”
Astrophys. J.
865
,
64
(
2018
).
70.
P.
Garaud
, “
Double-diffusive convection at low Prandtl number
,”
Annu. Rev. Fluid Mech.
50
,
275
298
(
2018
).
71.
J. M.
Brown
,
P.
Garaud
, and
S.
Stellmach
, “
Chemical transport and spontaneous layer formation in fingering convection in astrophysics
,”
Astrophys. J.
768
,
34
(
2013
).
72.
T.
Radko
and
D. P.
Smith
, “
Equilibrium transport in double-diffusive convection
,”
J. Fluid Mech.
692
,
5
27
(
2012
).
73.
P. Z.
Harrington
and
P.
Garaud
, “
Enhanced mixing in magnetized fingering convection, and implications for red giant branch stars
,”
Astrophys. J.
870
,
L5
(
2019
).
74.
J.
Goodman
and
G.
Xu
, “
Parasitic instabilities in magnetized, differentially rotating disks
,”
Astrophys. J.
432
,
213
(
1994
).
75.
H. N.
Latter
,
P.
Lesaffre
, and
S. A.
Balbus
, “
MRI channel flows and their parasites
,”
Mon. Not. R. Astron. Soc.
394
,
715
729
(
2009
).
76.
H. N.
Latter
,
S.
Fromang
, and
O.
Gressel
, “
MRI channel flows in vertically stratified models of accretion discs
,”
Mon. Not. R. Astron. Soc.
406
,
862
(
2010
).
77.
M. E.
Pessah
, “
Angular momentum transport in protoplanetary and black hole accretion disks: the role of parasitic modes in the saturation of MHD turbulence
,”
Astrophys. J.
716
,
1012
1027
(
2010
).
78.
P.
Garaud
,
A.
Kumar
, and
J.
Sridhar
, “
The interaction between shear and fingering (thermohaline) convection
,”
Astrophys. J.
879
,
60
(
2019
).