Zonostrophic instabilities in magnetohydrodynamic Kolmogorov flow

This paper concerns the stability of Kolmogorov flow u = (0, sin x) in the infinite (x,y)-plane. A mean magnetic field of strength B0 is introduced and the MHD linear stability problem studied for modes with wave-number k in the y-direction, and Bloch wavenumber l in the x-direction. The parameters governing the problem are Reynolds number 1/nu, magnetic Prandtl number P, and dimensionless magnetic field strength B0. The mean magnetic field can be taken to have an arbitrary direction in the (x,y)-plane and a mean x-directed flow U0 can be incorporated. First the paper considers Kolmogorov flow with y-directed mean magnetic field, referred to as vertical. Taking l=0, the suppression of the pure hydrodynamic instability is observed with increasing field strength B0. A branch of strong-field instabilities occurs for magnetic Prandtl number P less than unity, as found by A.E. Fraser, I.G. Cresser and P. Garaud (J. Fluid Mech. 949, A43, 2022). Analytical results using eigenvalue perturbation theory in the limit k->0 support the numerics for both weak- and strong-field instabilities, and originate in the coupling of large-scale modes with x-wavenumber n=0, to smaller-scale modes. The paper considers the case of horizontal or x-directed mean magnetic field. The unperturbed state consists of steady, wavey magnetic field lines. As the magnetic field is increased, the purely hydrodynamic instability is suppressed again, but for stronger fields a new branch of instabilities appears. Allowing a non-zero Bloch wavenumber l allows further instability, and in some circumstances when the system is hydrodynamically stable, arbitrarily weak magnetic fields can give growing modes. Numerical results are presented together with eigenvalue perturbation theory in the limits k,l->0. The theory gives analytical approximations for growth rates and thresholds in good agreement with those computed.


Introduction
The Kolmogorov flow, a periodic flow forced at a single wavenumber, is a classic flow to study due to its simplicity and its wide applications to the study of different geophysical and astrophysical systems.Its stability to linear disturbances is a classic problem first posed by Kolmogorov and studied by Meshalkin and Sinai (1961).These authors made use of continued fraction expansions to establish properties of the growth rate p(k), where k is a wavenumber in the y-direction, and determined a critical Reynolds number of Re c = √ 2. Close to onset of instability, for Re slightly larger than √ 2, it is large scale modes in y that are destabilised; more precisely, for Re = √ 2(1 + 3k 2 + • • • ) the most unstable mode has wave number k 1.This property allows the development of amplitude equations governing the flow on large space and time scales, obtained by Nepomniashchii (1976) and Sivashinsky (1985).Numerical simulations by She (1987) showed evolution from the most unstable scale to larger scales via an inverse cascade of vortex pairings, for a large scale allowed only in the y-direction.For large scales in both x-and y-directions, Sivashinsky (1985) showed evolution to a large-scale flow with chaotic temporal fluctuations, further explored in Lucas andKerswell (2014, 2015).
The stability problem posed by Kolmogorov is such a fundamental building block in stability theory that it has been elaborated in several studies by incorporating further physical phenomena.Frisch et al. (1996) included a β-effect, giving the gradient of a background planetary vorticity distribution; the gradient is oriented along the x direction (again following our conventions rather than those of the original paper), so that it does not interact directly with the transverse, basic state Kolmogorov flow u 0 .These authors derived an amplitude equation near to onset for a large scale in y, which they called the β-Cahn-Hilliard equation.For β = 0 this reduces to the PDE of Sivashinsky (1985) and simulations show that the inverse cascade of structures to large scales in y is arrested by the β-effect.These authors characterise the fundamental instability of the Kolmogorov flow as due to a negative eddy viscosity, in other words that the large-scale y-dependent modes have growth rate p = −ν E k 2 + • • • where the eddy viscosity ν E changes sign from positive below the threshold Re c = √ 2, to negative above, so destabilising the flow on large scales with the fastest growing modes determined by the next terms in this series (Dubrulle and Frisch 1991).
In terms of the geophysical motivation for these stability problems, any orientation of the background vorticity gradient, parameterised by β, with respect to the Kolmgorov flow is of interest.Manfroi and Young (2002) allow an arbitrary angle α between flow and gradient in a study of linear stability and nonlinear evolution using amplitude equations generalising those of Sivashinsky (1985) and Frisch et al. (1996).They find that the linear problem shows a delicate dependence of critical Reynolds number on angle α when unstable modes are allowed to adopt arbitrarily large scales in x and y.Another effect of geophysical relevance that may be included is stratification.Balmforth and Young (2002) considered the sinusoidal flow in the (x, z) plane with gravity in the −z direction and the flow directed in x, sinusoidal in z.These authors determined the behaviour of linear instabilities, depending on Reynolds, Richardson and Prandtl numbers, and derived an amplitude equation generalising that of Sivashinsky (1985).Simulations show that the inverse cascade of She (1987) is arrested by the presence of stratification over a wide range of parameters.
Relevant to the present paper, in astrophysical applications it is natural to introduce a magnetic field and study the coupled MHD system; as general motivation we note, for example, that the interaction between magnetic field, shear, and convection remains poorly understood in the solar tachocline (Hughes et al. 2007).Boffetta et al. (2000) considered the case in which a sinusoidal magnetic field (maintained by a source term in the induction equation) sits in a motionless fluid.This magnetic Kolmogorov system shows instabilities and an amplitude equation gives an inverse cascade to large scales.The recent paper Fraser et al. (2022) considers a background uniform magnetic field B 0 = (0, B 0 ) that is aligned with the Kolmogorov flow; this has no effect on the basic state flow but the elasticity of field lines affects perturbations depending on y, through the Lorentz force.These authors observe magnetic suppression of the instability first discussed by Meshalkin and Sinai (1961), as one might intuitively expect, but also two new families of unstable modes which only exist in the presence of magnetic field.One family exists for magnetic Prandtl number P < 1, for arbitrarily strong magnetic fields, provided the Reynolds number is above a threshold depending on P .This is studied numerically and growth rates obtained through asymptotic approximations for k 1; these authors refer to these modes as Alfvén Dubrulle-Frisch modes, as the instability can be linked to a change of sign of the eddy viscosity (Dubrulle and Frisch 1991).
Study of Kolmogorov flow instabilities is relevant to the formation of zonal flows in forced fluid systems, so-called 'zonostrophic instability' (Galperin et al. 2006).This process of jet formation has now been observed in many simulations, observations and experiments; see the representative studies, Vallis and Maltrud (1993), Read et al. (2007), Farrell and Ioannou (2008), Scott and Dritschel (2012), Srinivasan and Young (2012), Bouchet et al. (2013), andParker andKrommes (2014).Related to our work, Tobias et al. (2007) incorporated a magnetic field aligned with the x-direction of a planar fluid system with a β-effect present, a vorticity gradient in y.The system was driven by a body force with a given characteristic spatial scale.These authors observed the formation of jets in the x-direction for zero magnetic field, but then the suppression of jets, even at quite weak field strengths B 0 .For fixed nondimensional β, forcing and viscosity ν = 10 −4 , this process was explored by means of a series of runs with varying magnetic field strength B 0 and magnetic diffusivity η, and evidence for a threshold scaling law of B 2 0 ∼ η was observed (Tobias et al. 2011).Constantinou and Parker (2018) analysed Kelvin-Orr shearing wave dynamics for Rossby/Alfvén waves and the interplay between Reynolds and Maxwell stresses, providing evidence for this B 2 0 ∼ η threshold for jet formation.Durston and Gilbert (2016) focused on the couplings between large-scale zonal flow and zonal field in the presence of waves, calculating an effective viscosity and effective magnetic diffusivity, plus an effective cross transport term in which current gradients can drive the zonal vorticity; this and other transport effects are discussed in Chechkin (1999), Kim (2007), and Leprovost and Kim (2009).Parker and Constantinou (2019) interpret the presence or otherwise of jets in terms of the competition between a positive magnetic eddy viscosity term and a negative, purely hydrodynamic eddy viscosity.Note that while these studies of zonostrophic instability have many qualitative features in common with the topic of Kolmogorov flow instability, there are key difference that make any direct comparison difficult, even of scaling laws.The reason is that the studies referred to in this paragraph use a forcing which has a given spatial scale but is random in time, and the forcing is kept fixed while other parameters such as the viscosity, magnetic diffusivity, magnetic field or β are varied.In nondimensional terms, the input parameter is a Grashof number (formed from forcing strength, length scale and viscosity) in these systems (Durston and Gilbert 2016); the Reynolds number is then a diagnostic parameter, and can vary considerably in different regimes depending on the dominant balances in the Navier-Stokes equation between the forcing term, inertial term, viscous term and Lorentz force.However for stability of Kolmogorov flow, the basic state is fixed while the forcing is adjusted to maintain this: the input parameter is a Reynolds number, instead.
In the present paper we return to the classic set-up of steady, planar, Kolmogorov flow u 0 = (0, sin x) and consider the effect on its stability from magnetic field in the x-and ydirections.We find it convenient to refer to magnetic field in the y-direction, parallel to the flow as 'vertical field' and magnetic field in the x-direction, aligned with possible jet formation, as 'horizontal field' (even though gravity/stratification are not involved in our study).In §2 we set up the equations solved for linear perturbations with vertical magnetic field and in section 3 present numerical and analytical results, showing growth rates, thresholds and unstable mode structure.This section has common elements with the recent paper Fraser et al. (2022) (published while the present paper was in preparation); however we find it useful to set out the numerical results to compare with the later horizontal field case, and we present new analytical results in §3.1 for the 'weak vertical field branch'.The 'strong vertical field branch' in §3.2 is a primary focus for Fraser et al. (2022), and we give an alternative, matrix-based derivation of the asymptotic growth rate they obtain.Section 4 sets up the equations for horizontal magnetic field, with numerical results supported by analytical approximations in the limit k → 0 given in §5, together with the case of non-zero Bloch wavenumber , and §6 offers concluding discussion.Further analytical and numerical results will appear in Algatheem (2023).To keep the main body of the paper compact, we have developed analytical theory in appendices, building up in order of complexity rather than in the order in which the results are used.The method employed is perturbation theory for eigenvalues and eigenvectors of a matrix; naturally this is equivalent to methods used by other authors, but we find that is a systematic way of handling problems of increasing complexity, and gives insight into how couplings between individual flow and field modes can drive an instability.

Governing equations: vertical field
Our starting point is the system of equations for incompressible, constant density MHD, written in the form (2.1) Here the magnetic field is represented in velocity units, and f is an externally imposed rotational body force, used to maintain a given basic state for the system.In dimensional units this basic state is the Kolmogorov flow, namely the unidirectional flow in the (x, y)-plane specified by u 0 = U(0, sin(x/L)). (2.4) We use the length L and velocity U as the basis of our non-dimensionalisation, which results in the same equations (2.1-2.3)above but having ν now identified as an inverse Reynolds number, with Re = ν −1 , and η as an inverse magnetic Reynolds number, with Rm = η −1 .The non-dimensional Kolmogorov flow is then u 0 = (0, sin x). (2.5) Note that we drop the z-components of vectors where we can.A key quantity we will need is the magnetic Prandtl number (2.6) We will consider flows u and magnetic fields b lying in the (x, y)-plane, independent of z.For this we use a stream function ψ and vector potential a defined by (2.7) The governing equations may then be writtten in terms of the evolution of a scalar vorticity ω, and a: ) (2.10) Here J is the Jacobian of two functions in the plane and g is the z-component of the curl of the body force f .
We begin with the study of the stability of Kolmogorov flow in the presence of vertical magnetic field (that is, y-directed field) of strength B 0 (Fraser et al. 2022).Aiming for the most general set-up we also include a uniform horizontal flow (that is, an x-directed flow) of strength U 0 .We therefore adopt the following steady solution of the equations as our basic state, (2.11) or in our scalar-based formulation (2.12) The basic state magnetic field is shown in figure 1(a).The stability problem is parameterised by the four quantities {ν, B 0 , P, U 0 }.Note that while the mean horizontal flow specified by the parameter U 0 could be removed by a Galilean transformation, the Kolmogorov flow would then become a travelling wave.Thus, given we always take a steady Kolmogorov flow in the form (2.5), the effect of U 0 cannot be eliminated by this means.
To study the stability of this basic state we linearise, replacing and, droppping the subscript 1, we deduce the linear system (2.16) We now expand the fields in Fourier modes in x as (ψ, ω, a, j) = e pt+i x+iky n (F n , G n , H n , J n ) e inx + c.c. (2.17) Here p is the complex growth rate of the mode, k is the wavenumber in the y-direction with k > 0 without loss of generality, and is a Bloch or Floquet wavenumber in the x-direction Substituting these series into the linear equations (2.14-2.16)results in an infinite system of equations.For = 0 these may be written in the form: and for = 0 we simply replace n → n + wherever it appears (except as a subscript).This provides an eigenvalue problem giving a discrete set of eigenvalues with a dependence p(k, , ν, B 0 , P, U 0 ) in general, and the real part of the growth rate is unchanged on the replacement (k, ) → (−k, − ).
For a numerical solution we restrict −N ≤ n ≤ N for some integer N (typical values are 16, 32) and solve a discrete matrix problem written in the pentadiagonal form where ⊗ denotes the only non-zero entries.At a specified truncation N the (4N +2)×(4N +2) matrix is set up in Matlab, and an eigenvalue p with maximum real part is calculated.For a given parameter set {ν, B 0 , P, U 0 } the maximum real growth rate is defined as In what follows we will start by taking U 0 = 0, = 0 and only vary the vertical wavenumber k.The maximisation is then taken over a finite range of k-values, typically 100 values in the range 0 ≤ k ≤ 1, and any complex eigenvalues appear in complex conjugate pairs.We let k max (ν, B 0 , P ) be the corresponding maximising wave number, and we attach the appropriate (zero or positive) imaginary part to give p max (ν, B 0 , P ) as the (maximum) complex instability growth rate.It is then often instructive to plot Re p max , Im p max and k max .

Numerical and analytical results: vertical field
We use the numerical code as described above to produce eigenvalues so that we can explore the dependence on parameters.Our starting point is to investigate the effect of increasing the vertical magnetic field strength B 0 on the classic hydrodynamic instability of Kolmogorov flow.

Weak vertical field branch
Figure 2(a) shows the real part of the growth rate p(k, ν, B 0 , P ) for ν = η = 0.4 and so P = 1, plotted against k for given values of the magnetic field strength B 0 .Here B 0 is increased from zero in steps of 0.05 as we read down the family of curves.The top curve relates to the purely hydrodynamic case.As we increase B 0 we note two effects: first the peak is reduced, in other words the magnetic field acts to suppress the instability, as found by Fraser et al. (2022).Secondly, for large scales, namely small k, a new branch of decaying modes appears, with growth rates largely independent of field strength.Figure 2(b) shows the imaginary part of the growth rate p.This is zero for the purely hydrodynamic case and remains zero for this branch as it is suppressed by the field.The new branch for low k has a non-zero imaginary part which becomes more prominent as B 0 is increased and we read up the curves.Some investigation shows that the new branch is in fact a damped Alfvén wave on the vertical magnetic field.For zero background flow u 0 , a vertical field supports Alfvén waves with and the real and imaginary parts are shown dotted in figure 2. Since the Alfvén waves are modified by the background Kolmogorov flow, the agreement is not exact, but this is clearly the origin of these small-k damped modes.
A typical unstable mode is shown in figure 3 for parameter values corresponding to the peak in the lowest curve in figure 2 hydrodynamic problem as we increase B 0 from small values, we refer to this as the weak vertical field branch, to be contrasted with a strong field branch we encounter shortly.
Having seen a particular example of how the magnetic field suppresses the hydrodynamic instability by plotting p(k, ν, B 0 , P ), we now show results where we maximise over k for each set of the parameters.Figure 4(a) shows numerical results for Re p max (ν, B 0 , P ) with P = 1 as a colour plot across the (ν, B 0 )-plane.The white curve shows the threshold for instability (actually set for Re p max having a small positive value).The horizontal axis B 0 = 0 is the hydrodynamic case, where the white curve crosses at ν c = 1/ √ 2. Instability occurs in the region below the white curve, and we can see that it is suppressed as B 0 increases, up to the point where B 0 0.7 and the instability is entirely eliminated.We do not show Im p max , which is zero within the region of instability.
We can develop perturbation theory (as in, for example, Frisch et al. 1996, Manfroi andYoung 2002) to calculate approximate growth rates valid for k → 0. In appendix C we give the details.One key result is the formula (C.19), giving the growth rate, and showing clearly how the effect of the magnetic field is to suppress the hydrodynamic B 0 = 0 instability (as seen in figure 2(a)).For unstable modes it is necessary that ν < 1/ √ 2 and in this case maximising the growth rate over values of k gives (3.24) Putting p max = 0 gives the threshold of instability as the straight line in the (ν, B 0 )-plane: Figure 4(b) shows the theoretical growth rate and the threshold marked by a white (straight) line, showing good agreement with the full numerics.The perturbation theory is developed about the point ν c = 1/ √ 2, B 0 = 0 of the onset of the hydrodynamic instability.Hence the agreement is particularly good near this point; elsewhere the theory gives results that are qualitatively correct only, for example the theoretical value of B 0 which suppresses instability for all ν is given by with B * 0.58 for P = 1 in panel (b), which is a little lower than the actual value B * 0.7 seen for the numerical results in panel (a).

Strong vertical field branch
Although the magnetic field acts to suppress the instability for magnetic Prandtl number P = 1, this is not the whole picture, and investigations for P < 1 show the presence of a strong vertical field branch, as found by Fraser et al. (2022).This strong vertical field branch is analysed in appendix B, using a scaling in which B 0 = O(k −1 ) as k → 0. The pertubation theory then involves a leading order undamped Alfvén wave with frequency p 0 = ikB 0 = O(1).The coupling of this wave with the flow field leads to potential instability with a growth rate given in (B.9) as: (3.28) For example if P = 1/2 then ν * = 1/ √ 12 0.28, in good agreement with the vertical white lines in figure 5. Thus the instability persists for arbitrarily large magnetic fields, provided the viscosity is below this Prandtl-number dependent threshold, in other words provided the Reynolds number Re = ν −1 > ν −1 * .
All the above results have been taken for Bloch wave number = 0. Introducing brings in an extra degree of freedom and allows the possibility of new instabilities.However in the vertical field case increasing | | from zero has only a stabilising effect (Fraser et al. 2022).For example figure 6 shows growth rates in the ( , k)-plane for weak and strong field cases in panels (a) and (b).The vertical axis gives = 0 and we see an island of unstable modes in each case.However to either side of the vertical axis, the growth rates are diminished and allowing = 0 has little impact.For this reason we will not consider further for the vertical field case, except to mention that the theory in appendix B may be extended to incorporate = 0, as detailed in Algatheem (2023).
4. Governing equations: horizontal field, with U 0 = 0 In this section we study the stability of Kolmogorov flow in the presence of horizontal field B 0 , and to reduce the complexity of the problem we restrict to the case of zero horizontal mean flow U 0 = 0. We thus adopt the following basic state, a steady solution of the equations (2.1-2.3): or in the scalar formulation (4.30)The basic state magnetic field is shown in figure 1(b), with horizontal field lines distorted by the background Kolmogorov flow, becoming increasingly extended in the limit of small η.Note, to pick up a comment in the introduction, that the body force required to maintain the Kolmogorov flow increases with field strength B 0 and with magnetic Reynolds number Rm = η −1 , unlike in many large-scale simulations of zonostrophic instability, where the magnitude of a random forcing is held fixed, while other parameters are varied.
The stability problem is parameterised by {ν, B 0 , P, U 0 = 0}.The corresponding linear system is where the fields represent the perturbation to the basic state.The resulting equations for the Fourier modes in x are ) for = 0 and, as elsewhere, for = 0 we replace n by n + .This infinite system of linear equations may then be truncated and set up as a matrix eigenvalue problem, analogously to that in (2.20) for vertical field; the matrix now takes a heptadiagonal form.

Numerical and analytical results: horizontal field
We have used Matlab to obtain growth rates p(k, , ν, B 0 , P ) (here U 0 = 0) and we focus first on the case = 0.

Instability for Bloch wave number = 0
With zero Bloch wave number , the instability has periodicity 2π in the x-direction and 2π/k in the y-direction.Figure 7 shows plots of the growth rate p(k, ν, B 0 , P ) against k for ν = η = 0.1 (P = 1) and B 0 increasing as detailed in the caption.Focusing on the real part of p in panel (a) we observe that the magnetic field initially suppresses the instability, going from the blue B 0 = 0 curve to the lower, red B 0 = 0.05 curve.However increasing B 0 further, the green B 0 = 0.1 curve reveals a double-peaked growth rate and then these two peaks increase as B 0 is increased, as indicated in the figure caption.For these stronger fields, the second peak is the lower of the two, and is associated with non-zero imaginary part Im p of the growth rate, as shown in panel (b), while the dominant instability of the first peak has Im p = 0.In fact for all our studies of instability for horizontal field with U 0 = 0 and = 0, we observe that the dominant instability is always direct or non-oscillatory, that is Im p = 0.
To give a more global picture of these results for horizontal field, we now show p max (ν, B 0 , P ) as a colour plot in the (ν, B 0 )-plane with P = 1 in figure 8(a), with the white line denoting the instability threshold Re p = 0.For modest magnetic fields we observe the suppression of the purely hydrodynamic instability as in the case of vertical field in figure 4. However as B 0 is increased another branch of instability emerges from the bottom left of the panel 8(a) and shows increasing growth rates, particularly for smaller viscosities ν.This new branch of instabilities is perhaps not surprising (Durston and Gilbert 2016), given that the basic state horizontal field in (4.29) and depicted in figure 1(b) has a wavey structure, and for P = 1 becomes increasingly convoluted as η = ν is decreased.Other values of P give plots that are similar in structure.
To gain analytical results and understanding, in appendix A we discuss perturbation theory for the horizontal field system, taking the limit k → 0 while retaining B 0 and other parameters of order unity.The resulting leading order equations involving G 0 , G ±1 , H 0 and H ±1 split into the two independent 3 × 3 matrix systems giving the two branches of instability evident in figure 8(a).We discuss them in turn.
The first system involves G 0 and not H 0 , in other words is dominated by a large-scale flow and not a large-scale field.We call this the G 0 or flow branch of the horizontal field instability.Analysis gives equation (A.31), which we reproduce here as This gives the leading growth rate as a function of the parameters times k 2 ; it represents an eddy viscosity ν E seen by large-scale modes and the instability is marked by this quantity becoming negative.While it is not possible to maximise this expression over k to gain a complete analysis of the instability, it does give the instability threshold Re p = 0, by setting (5.37) This formula is plotted as the blue curve in panel 8(b) and shows good agreement with the numerical results for the lower branch in panel 8(a).For B 0 = 0 we recover the hydrodynamic result ν c = 1/ √ 2, and this analysis tells us how this hydrodynamic instability, domimated by the large-scale flow in G 0 , is suppressed by interaction with the magnetic field.If we maximise B 2 0 as a function of ν in (5.37), we find that this occurs at (5.38) and putting this into B 2 0 gives an unwieldy expression for the threshold value B * above which the horizontal field suppresses the Kolmogorov instability.We do not present it here but give further discussion in §6.
Note also that from (5.37), (5.39) and so the instability emerges with this slope from the origin ν = 0, B 0 = 0 of the figure.A typical example of the instability is shown in figure 9: the stream function in panel (a) shows its nature as a zonostrophic instability giving horizontal jets, with modifications to the field in panel (b).
The second system arising from perturbation theory involves H 0 and not G 0 : it is dominated by a large-scale magnetic field and so we refer to this as the H 0 or field branch of the horizontal field instability.The result of perturbation theory gives (A.37), reproduced here as (5.40) The onset of instability again can be interpreted as a transport quantity becoming negative, here the eddy magnetic diffusivity η E (see Chechkin 1999).Note that this instability, being a directly growing instability does not connect to the strong-field branch of the vertical field (which is an over-stable wave).
The threshold for instability is given by [ (5.41) This curve is plotted on figure 8(b) in red and again shows good agreement with the numerical results for the field branch in panel (a).Note that for fixed P , B 0 → ∞ as ν → ν * with ν * = P 3/2, (5.42) and so a viscosity larger than this is enough to prevent the field branch instability no matter how strong the field.We also have for small fields and viscosities that the threshold (5.41) is given by and so for P = 1 both thresholds (5.39) and (5.43) emerge from the origin with the same slope, though for general P they are different.A typical example of the instability is given in figure 10.The perturbation flow now does not have a zonostrophic jet structure, but shows closed eddies in panel (a).Panel (b) however shows a banded structure in the magnetic field (showing the dominant role of H 0 ), and indicates a tendency for the background mean field to segregate into bands of stronger and weaker horizontal field.

Instability for Bloch wave number = 0
Finally, we consider horizontal field for the case of non-zero .This allows an instability to take up a scale 2π/k in the y-direction and 2π/ , as → 0, in the x-direction.It turns out that instabilities can occur for = 0 even when the system is stable for = 0, in the case of horizontal field (unlike the situation for vertical field).11(a).The corresponding growth rates are quite small, and so this is otherwise not immediately evident on these plots.
We therefore turn to theory for = 0, developed in appendix D, which gives a growth rate in (D.13) of This approximation reveals an instability that crucially relies on having a non-zero Bloch wavenumber, = 0, with and k both small.If we fix the parameters ν, P and B 0 we can consider the growth rate p as a function in the ( , k) plane.Setting the quantity inside the square root to zero to find a threshold, we see that the region of instability is demarcated by the pair of straight lines given by . (5.45) The growth rate p tells us about the instability growth rate as we increase k and from zero, but to find how this eventually decreases, we would need to go to next order in perturbation theory, which is impractical and unlikely to be informative.To give a qualitative feel for the growth rate we will add on the diffusive suppression term − 1 2 (ν +η)k 2 that is certainly present at next order, and look at To see how this works, figure 12(a) shows growth rates plotted in the ( , k)-plane for B 0 = 0.2 and ν = η = 0.75, parameters corresponding to stability for = 0.There is now a region of instability taking a 'butterfly' form, outlined by the white curve Re p = 0.The straight blue lines are given by (5.45) and are tangential to the white lines at the origin, confirming the theory.In panel (b) we show an analogous figure for the 'fixed up' growth rate in (5.46).The agreement between the results in the two panels (a) and (b) is excellent near the origin, but then further out the agreement is only qualitative, and quite rough, as we might expect.For an example of an unstable configuration, figure 13 shows the flow and field for the fastest growing mode in figure 12.We observe the structure of a large-scale jet-like structure but at an oblique angle to the y-direction, this being allowed by the non-zero Bloch wavenumber .
Returning to the bigger picture, for instability at a general points in the (ν, B 0 ) plane we need the quantity [• • • ] inside the square root in (5.44) to be positive for some values of k and .Equivalently, it corresponds to requiring that the straight lines in (5.45) have a finite slope.It can be checked that this gives a threshold for instability: (5.47) Values of field smaller than this give instability and so this formula gives the threshold curve in the (ν, B 0 ) plane for this family of = 0 instabilities.This threshold is shown as a dashed curve in figure 8(b).We thus see in this panel that allowing = 0 means that almost the whole of the parameter ranges shown give instability, all except for a small curved triangular region at the top right, and this is in agreement with the white curve obtained numerically in figure 11(a).The agreement is not perfect because the growth rates in the top right corner become small and the unstable region shrinks away in the ( , k)-plane, as the thresholds are approached, and so the precise location of the white curve becomes hard to resolve without further work.
Note that for P = 1, from (5.47) the instability is cut off at ν = √ 3, and in general the instability requires ν < ν * = P (2 + P ). (5.48) Numerical experimentation shows that for ν less than this value, the region of instability in the ( , k)-plane becomes vanishingly small as B 0 tends to zero or tends to the value given in (5.47).The maximum magnetic field for the = 0 instability found here is given by maximising B 0 over ν: the maximum occurs at and we pick this up in the final discussion section, §6.

Discussion
In this study we have explored instability of the classic Kolmogorov flow in the presence of magnetic field which is either vertical, aligned with the flow, or horizontal, aligned with possible jet formation.In the vertical field case we have obtained new analytical results for the maximum growth rate (3.24) and magnetic field threshold (3.25), that show the suppression of the original instability found by Meshalkin and Sinai (1961) by vertical magnetic field.In particular we have obtained a value B * in (3.26) as an estimate of the field magnitude required to suppress the Kolmogorov instability.This corresponds to a threshold that is For the strong field branch we have confirmed the numerical results of Fraser et al. (2022) and provided an alternative derivation of their growth rate formula (3.27) and threshold (3.28).We note that these authors also find a third branch of instabilities, which they term 'varicose Kelvin-Helmholtz' modes, which we have not observed at the Reynolds numbers we have used.
The case of horizontal field, broadly relevant to several studies of jet formation where the field is aligned with potential jets (Tobias et al. 2007, Durston and Gilbert 2016, Constantinou and Parker 2018), shows more complex structure, unsurprisingly given the wavey nature of the magnetic field in the basic state, seen in figure 1(b).We observe again the suppression of the purely hydrodynamic instability, the flow or G 0 branch, when the magnetic field strength is increased, with a threshold of B 0 = B * for complete suppression with B * given by substituting ν 2 from (5.38) into (5.37).For large P , we have ν 2 1/4 while for small P , ν 2 P/2.The value of B * for suppression then amounts to: Interestingly this bears comparison with Tobias et al. (2007) who have ν = 10 −4 fixed and 10 −1 ≤ η ≤ 10 −6 in their runs.For the greater values of η used, P 1 and so a threshold B 2 * ∼ η is indicated above, and found in these full numerical simulations.Also, note that at this threshold we have that the forcing magnitude is fixed in magnitude in (4.29) and so there is, at least roughly, a correspondence of working with fixed forcing amplitude as in their paper and the fixed Kolmogorov flow in ours.However we should remark that these authors use a non-zero value of β and we have zero, and so further work would need to be done to make a sound comparison.
A feature of the horizontal field problem is that for increasing magnetic field strengths a further branch of instabilities emerges, the field or H 0 branch, also seen in Durston and Gilbert (2016).Analytical formulae for the threshold of instability are given for each branch.The field branch exists provided the Reynolds number is above a threshold, Re > ν −1 * with ν * given in (5.42).Allowing a Bloch wavenumber = 0 in the x-direction, in addition to the wavenumber k in the y-direction, allows a new branch of instabilities.In particular a magnetic field, no matter how weak, can destabilise the Kolmogorov flow provided the Reynolds number Re = ν −1 > ν −1 * with ν * given in (5.48).For example at P = 1, the purely hydrodynamic instability is present for Re > √ 2 but an MHD instability is present for arbitrarily weak but non-zero horizontal magnetic field provided Re > 1/ √ 3.For sufficiently large magnetic field this instability is again suppressed, and making use of (5.49) (with ν 2 √ 2P 3 for small P and ν ( √ 2 − 1)P 2 for large P ) we find a threshold Note that the order of limits could be important: in our discussion in this paper we are fixing any value of P and then allowing k and to tend to zero.Other limits are possible and could be explored by appropriate scalings in our calculations.
Underlying our study is matrix eigenvalue perturbation theory as set out in the appendices, a flexible tool for these types of problems.We find it gives greater clarity than using a multiple scales formulation or applying perturbation theory to roots of a polynomial, even though all these methods are ultimately equivalent.Note that while many of the instabilities seen by us and by other authors can be characterised as involving a negative eddy viscosity term −ν E k 2 or a negative eddy magnetic diffusivity term −η E k 2 at large scales, the growth rate p(k, ) in the case of horizontal field shows a complicated dependence on k and in (5.44).Although this = 0 instability occurs at arbitrarily large scales, it cannot be categorised as involving a simple negative eddy transport effect.This arises because we are applying perturbation theory to a repeated eigenvalue of the limiting k → 0, → 0 problem.Looking to the future, it would be interesting to pursue further research on the Kolmogorov flow as an MHD system, particularly on the nonlinear evolution of instabilities and inverse cascades (Fraser et al. 2022, Algatheem 2023), and on the interaction of magnetic field with a β-effect and Rossby waves.
A.3.Field or H 0 branch In the second system (A.9-A.10), the large-scale field H 0 is present but no large-scale flow.We write the system as M v = pv with The matrix series for M now has .33)We have that the inverse of the 2×2 block of M 0 is given as in (A.26) with ν and η interchanged and the same ∆.The calculation proceeds as before.At leading order in the eigenvalue problem (A.14) we take the same solution as that given in (A.27).At first order, we have .34)This gives p 1 = 0 and we solve (A.16) for v 1 as At the next order (A.21) yields p 2 and so or, with P = ν/η, Further analysis commmences from equation (5.40).
Appendix B: Vertical strong field, with U 0 = 0, = 0 In this appendix we turn to the vertical field system.Here there are two types of instability and two analyses that we will set out in this appendix and the next one.The calculation in this appendix is designed to capture the properties of the strong field branch seen for η > ν in figure 5 and is equivalent to that set out in Fraser et al. (2022).Mathematically we need to consider the limit when B 0 → ∞ as k → 0, and we find that relating these via B 0 = O(k −1 ) is most informative.We will set the Bloch wavenumber = 0 and take no mean flow U 0 = 0.
If we write out the vertical field equations truncated to G 0 , H 0 , G ±1 and H ±1 , rewrite in terms of G ± and H ± defined in (A.5), then we obtain the equations (without any further approximation) in the form M v = pv with where G 0 = k 2 G 0 as usual.The fields G + and H + are decoupled (as we have U 0 = 0) and so not considered further.Before expanding M in powers of k, for strong vertical field we rescale B 0 = k −1 B 0 with B 0 fixed in the limit k → 0. Then expanding M gives the matrices For an approximate growth rate p we use the expansion (A.14) and solve order by order.At leading order (A.15) we focus on the eigenvalues p 0 = ±iB 0 of M 0 , corresponding to largescale undamped Alfvén waves.We will focus on the upper sign without loss of generality, and take Here w 0 is the left eigenvector as usual, with w 0 (M 0 − p 0 ) = 0 and w 0 v 0 = 2.
Moving to the first order, from (A.20) we rapidly find We now need to solve (A.16) for v 1 .To find a solution we clearly need only invert the 2 × 2 lower right block of M 0 − p 0 to calculate with the inverse determinant ∆ now defined as With this it is straightforward to calculate p 2 from (A.21) Taking the real part of p 2 , and then putting back B 0 = kB 0 and p = p This is taken up in the main body of the paper as (3.27).
Appendix C: Vertical weak field, with U 0 = 0, = 0 We now continue with our analysis of vertical field instabilities for = 0 and U 0 = 0. We studied the strong field branch with B 0 = O(k −1 ) in the previous appendix.In the present appendix we will take B 0 = O(k 2 ): this addresses the branch of vertical field instability as seen in figure 4 and allows us to resolve the question of how magnetic field suppresses the purely hydrodynamic instability onset and reduces the critical value of ν below ν c = 2 −1/2 .We will see that our results will be correct qualitatively even when B 0 is as large as order unity while k → 0. Nonetheless it is convenient to refer to this as the 'weak field branch' to contrast with the strong field branch in the previous appendix.
We start with the matrix system (B.1) for instability in the presence of vertical field: M v = pv.However before expanding M 0 in powers of k we first rescale B 0 = k 2 B 0 , and set Here we are going to develop perturbation theory around the critical point ν c (U 0 ) for the purely hydrodynamic problem, with ν c (0) = 2 −1/2 .Expanding in powers of k yields 3) It will be convenient to let In keeping with the hydrodynamic case we will be expanding the system about a zero eigenvalue p 0 = 0 of M 0 with a flow field specified in G 0 .However we should note that p 0 is a twice repeated eigenvalue and we have two left and two right eigenvectors which we distinguish with † and ‡: v T 0 † = w 0 † = (1, 0, 0, 0, 0, 0), v T 0 ‡ = w 0 ‡ = (0, 0, 0, 1, 0, 0).(C.5) In the perturbation theory for a general matrix M we would need to take v 0 as a general combination of v 0 † and v 0 ‡ , to be determined further in the expansion.However here, to avoid unnecessary algebra we will jump to the solution we need, and take for a flow-dominated eigenfunction.We verify that this works as we delve into the expansion.
Looking to the first order equation (A.16) we apply w 0 † and w 0 ‡ to the left-hand side, which only gives p 1 = 0 and we have However when we solve for v 1 we can add not only a multiple of v 0 to the solution (which would have no effect in the calculation) but also a multiple of the purely magnetic eigenvector v 0 ‡ .Thus we solve for v 1 in the form where b is an unknown constant, to be determined.
At the next order we will aim to take p 2 = 0 so as to push various the effects down the series in powers of k: this is achieved if we fix ∆ = 2. Thus at second order we set p 2 = 0, ∆ = 2, (C.9) and note that this then fixes This is the critical value for onset for the pure hydrodynamic case, with a mean flow but with zero magnetic field.A strong enough mean flow |U 0 | > 2 −1/2 is enough to suppress any instability and so from now on we take U 2 0 < 1 2 so that ν 0 is real and positive.
In the case of zero mean flow, U 0 = 0, we have ν c = 2 −1/2 , δ = 2P 2 , and this simplifies to We continue the discussion in the main body of the paper; see equations (3.24, 3.25).Note that going to fourth order here suggests that the modes G ±2 and H ±2 might also be needed to give a correct evaluation of p 4 ; this needed to be checked and was -we found that the couplings are too weak, in terms of the powers of k involved, to give a contribution to the growth rate to the order taken above.
Appendix D: Horizontal field, with U 0 = 0, = 0 Our final calculation brings in the Bloch wavenumber .This can be done in the case of vertical field (Algatheem 2023), but there increasing from zero seems to have only the effect of suppressing the = 0 instability, so we do not consider this further.We instead study horizontal field, where having = 0 can enhance instability.We take U 0 = 0 to keep the problem manageable.We omit straightforward but messy details and in fact will only go up to first order in perturbation theory.
Our starting point is the equations (4.34, 4.35) with n replaced by n + , and we consider only the modes G 0 , G ±1 , H 0 , H ±1 .We set, as in the original horizontal field problem (appendix A), Once we have = 0 in the problem, we have to ask how scales as k → 0. It turns out that the appropriate scaling to gain useful results is = k = O(k), (D.2) so we hold and G 0 constant while k → 0. We now follow the usual procedure of making these substitutions, expressing the equations in terms of G 0 , G ± , H 0 , H ± and writing the system as and We are now ready to calculate p.The matrix M 0 has now lost the attractive block structure present in the earlier expansions as a consequence of the scaling of .Nonetheless M 0 has a double zero eigenvalue p 0 = 0 with right eigenvectors v †0 = (1, 0, 0, 0, 0, 0) T , v ‡0 = (0, 0, 0, 1, 0, 0) T , (D.7) and left eigenvectors w 0 † = (1, 0, 0, 0, 2 ∆η −1 iB 0 (ν + η), 2 ∆η −1 (η 2 − B 2 0 )), w 0 ‡ = (0, 0, 0, 1, 0, 0).(D.8) In the previous case of a double-zero eigenvalue (in appendix C) we anticipated the structure of v 0 (as dominated by the hydrodynamic field G 0 ).Here we cannot do so and so we set v 0 = bv 0 † + cv 0 ‡ , (D.9) for some constants b and c.
Now, looking at the first order equation (A.16), namely p 1 v 0 = (M 0 − p 0 )v 1 + M 1 v 0 with p 0 = 0, we can apply either of the two vectors w 0 † and w ‡ on the left, to gain two equations, We have gained this equation by only going to the first order matrix M 1 , but it reveals an instability that crucially relies on having a non-zero Bloch wavenumber, = 0. We continue our discussion in the main body of the paper, from (5.44).
.21) with the maximum taken over a grid of k and values.

Figure 5 .
Figure 5. Shown are (a) instability growth rate Re pmax and (b) frequency Im pmax for vertical field plotted in the (ν, B 0 ) plane, with P = 1 2 .The results shown are obtained computationally and the white curve in panel (a) shows zero growth rate.
Figure 5 shows (a) the real part and (b) imaginary part of the growth p max (ν, B 0 , P ) for P = 1/2, that is η = 2ν.The threshold Re p max = 0 is shown as a white curve in both figures.Looking from the bottom of the panel 5(a) (increasing B 0 ) we see that the curving white line, showing the weak field branch in panel 4(a), turns to become a near-vertical line, demarcating a new branch with non-zero frequency Im p max evident in panel 5(b).

Figure 11
Figure 11.(a) Instability growth rate Re p and (b) imaginary part Im p shown for horizontal field, plotted in the (ν, B 0 ) with P = 1, U 0 = 0, and = 0.The maximising values of k and of are shown in panels (c, d) respectively.

Figure 11
Figure 11(a) shows a similar general structure to figure 8(a), though note that the white curve has all but disappeared.The purely hydrodynamic instabilities are suppressed by increasing B 0 , i.e. the G 0 or flow branch; this is not really visible now in panel (a) but the branch is evident in panel (c) showing k max , with max = 0 in panel (d) here.Then, the H 0 or field branch of instability appears clearly in panel (a) for larger B 0 .However looking at panels (b) and (d) now the field branch has regions with non-zero imaginary part to the growth rate, regions with increasing from zero, and regions with constant at = 0.5 indicating instabilities with 4π periodicity in x (further investigated in Algatheem 2023).Returning to figure 11(a), the reason for the loss of the white curve, compared with figure 8(a), is that there is now a new instability, reliant on having = 0 which means that almost all of the region of stability in the = 0 figure 8(a) is now unstable in figure11(a).The corresponding growth rates are quite small, and so this is otherwise not immediately evident on these plots.