Modification of the resistive tearing instability with Joule heating by shear flow

We investigate the influence of background shear flow on linear resistive tearing instabilities with Joule heating for two compressible plasma slab configurations: a Harris current sheet and a force-free, shearing magnetic field that varies its direction periodically throughout the slab, possibly resulting in multiple magnetic nullplanes. To do so, we exploit the latest version of the open-source, magnetohydrodynamic spectroscopy tool Legolas. Shear flow is shown to dramatically alter tearing behaviour in the presence of multiple magnetic nullplanes, where the modes become propagating due to the flow. Finally, the tearing growth rate is studied as a function of resistivity, showing where it deviates from analytic scaling laws, as well as the Alfv\'en speed, the plasma-beta, and the velocity parameters, revealing surprising nuance in whether the velocity acts stabilising or destabilising. We show how both slab setups can produce growth rate regimes which deviate from analytic scaling laws, such that systematic numerical spectroscopic studies are truly necessary, for a complete understanding of linear tearing behaviour in flowing plasmas.


I. INTRODUCTION
In many plasma configurations, from astrophysical to laboratory settings, the most dramatic and violent events are often driven by magnetic reconnection.During the reconnection process, magnetic field lines break and reconnect, drastically altering the magnetic topology, e.g. by unravelling a knotted field structure.Consequently, stored magnetic energy is converted into thermal and kinetic energy, accelerating particles and producing matter outflows.It is frequently observed in various phenomena, like coronal mass ejections, 1 the heliospheric current sheet, 2 and Earth's magnetosphere. 3 Due to its importance and ubiquity, reconnection has been a topic of great debate ever since Sweet 4 and Parker 5 first proposed their model, which was later followed by Petschek's model 6 to achieve faster reconnection rates in accordance with observations.Irrespective of the specific reconnection process, its initialisation in current sheets is intrinsically linked to tearing instabilities, which occur due to magnetic shear, and create the necessary reconnection points (Xpoints).These tearing instabilities come in collisional 7 (resistive) and collisionless 8 (mediated by the Hall effect and/or electron inertia) varieties, both of which have been studied extensively to understand their role in triggering fast magnetic reconnection.
To reach the fast reconnection regime, the influence of several effects on the resistive tearing mode has been studied for decades.For thin current sheets with a width on the order of the ion inertial length, the importance of the Hall term in the generalised Ohm's law was already shown forty years ago by Terasawa, 9 who showed that it enhances the growth rate of the resistive tearing mode.1][12][13][14][15] However, the role of the Hall term in magnetic reconnection is not limited to the enhancement of the resistive tearing growth rate.It has also been linked to finite Larmor radius effects, [16][17][18] it couples to the anisotropic electron pressure tensor, which enhances reconnection rates, 19,20 and it can induce a transition to whistler-mediated reconnection. 8,21,22Finally, Liu et al. 23 recently showed that the Hall term plays an important role in attaining the proper geometry for fast reconnection.
However, kinetic simulations suggest that the collisionless terms in the generalised Ohm's law do not dominate in collisional plasmas until the current sheet thins to the ion inertial scale, 24,25 though collisionless reconnection does occur at larger scales in plasmas with negligible collisionality due to electron inertia. 26Hence, for thicker, sufficiently collisional current sheets, the initial unstable perturbation is expected to be of the resistive tearing variety, with a negligible contribution due to collisionless terms.How these sheets then reduce to the ion inertial scale to establish fast reconnection is not fully understood, though processes like fractal reconnection 27 have been suggested.This then begs the question which other factors might influence the growth rate outside of the Hall regime.It is known that the growth rate is affected by various environmental factors, notably by the background flow. 28,29In this endeavour, the influence of equilibrium flow on the resistive tearing mode has already been studied extensively using both analytic [30][31][32][33][34] and numerical 28,29,[35][36][37] methods.Flow's influence is also not limited to the resistive tearing type, as its impact on collisionless tearing has been shown for flow parallel to the reconnecting field 38 and to a guide field. 39Furthermore, the role of flow, and particularly flow shear, extends well beyond the modification of tearing growth rates, as it may introduce the Kelvin-Helmholtz instability (KHI) into the system.This instability regularly becomes the dominant instability, 30,40 though it is also observed to co-exist with the plasmoid instability (secondary tearing) in current sheets, 41 at plasma interfaces, 42 and in turbulent reconnection. 43 profound understanding of instabilities and how to suppress them is also of vital importance in laboratory plasmas, and particularly for fusion research.In fusion devices like tokamaks, tearing instabilities lead to the formation of magnetic islands, which in turn disrupt plasma confinement.There have been many experiments, 44,45 linear studies, [46][47][48] and nonlinear simulations [49][50][51] studying the effects of both toroidal and poloidal flow.In these torus-like geometries, it is now generally accepted that the toroidal flow shear has a stabilising influence on the system by suppressing island formation.
In this work, we revisit the linear analysis of the resistive tearing mode in the presence of background flow using numerical means.Unlike earlier literature, we eliminate the need for approximations, like incompressibility, by employing the magnetohydrodynamic (MHD) spectroscopic code LEGO-LAS 15,52,53 (https://legolas.science).With this code we explore plasma stability parametrically for a selection of configurations and show that equilibrium flow can both enhance and suppress the growth rate of the resistive tearing mode depending on the parameter regime.The effect of equilibrium flow has been studied analytically by Chen and Morrison 34 and we connect our results to the power law predictions from their study.

II. SETUP AND CONVENTIONS
To study the linear properties of the resistive tearing instability, we linearise the dimensionless, compressible, resistive MHD equations around a one-dimensionally varying equilibrium (xdependence only) and assume Fourier solutions for perturbed quantities Here, ρ, v, T , and B are the usual quantities density, bulk velocity, temperature, and magnetic field, respectively, with subscripts 0 and 1 differentiating between equilibrium quantities and perturbations.Furthermore, p = ρT denotes the pressure, J = ∇ × B the current density, and η = 10 −4 , unless specified otherwise, and γ = 5/3 are the resistivity and adiabatic index, respectively.Finally, k = k y êy + k z êz and ω are the wave vector and frequency.The resulting algebraic problem is solved for the frequency and corresponding perturbed quantities by LEGO-LAS.
Note that our choice of energy equation, Eq. ( 3), including Joule heating, along with compressibility and a constant resistivity, differs from the isothermal closure used in the original derivation of the tearing mode. 7

A. Configurations
Now consider a semi-infinite plasma confined in the xdirection between two perfectly conducting plates, described in Cartesian coordinates.In such a plasma, we examine the resistive tearing instability for two separate equilibrium configurations: a Harris current sheet and a force-free magnetic field.

Harris current sheet
For the first configuration, we consider the popular Harris current sheet, as presented in Ref. 28.In their simulation setup they consider a typical Harris current sheet which they supplement with a similar velocity profile and a uniform density ρ 0 = 1.The temperature profile is simply obtained by demanding that the total equilibrium pressure (the sum of plasma pressure and magnetic pressure) is constant, 54 i.e.
with a solution for any constant p c , set to p c = 1 throughout this work, unless specified.The profiles are visualised in Fig. 1(a).Since B 0 (x) vanishes at x = 0, so does for any wave vector k, which we choose to be along the yaxis, k = k y êy .Hence, the magnetic nullplane (or resonant plane), defined as the location where F vanishes, 7 is always at x = 0 in this configuration, and magnetic reconnection will occur here.The locations of the perfectly conducting walls at x w = ±15 a B are chosen such that their effect on the tearing instability is negligible (according to Ref. 55 the effect of the conducting walls is negligible if they occur at a position |x w | ≳ 10 a B ). Due to the sharp transitions in the equilibrium profiles near the origin and approximately constant behaviour away from the center, the problem is solved using a non-uniform grid concentrated near the origin, as described in App. A. Equilibrium with the force-free magnetic field, Eqs.(11).Parameter values were chosen for visual clarity.

Force-free magnetic field
In the second configuration, the magnetic field strength is fixed whilst its direction varies continuously throughout the plasma slab, similar to the configuration in Refs.56 and 57.This force-free magnetic field profile is complemented with a constant density and temperature, and a linear velocity profile, as defined in Ref. 58, Sec.14.3.3, where ρ c , v c , and α are constant parameters.As the notation suggests, β 0 is the equilibrium plasma-β .Note that |B 0 | = 1 such that the dimensionless Alfvén speed c A equals c A = 1/ √ ρ c .The equilibrium is visualised in Fig. 1(b).For this equilibrium we choose k proportional to êy (and thus parallel to v 0 ), such that there is a magnetic nullplane at x 0 = 0, where F(x 0 ) = 0.The perfectly conducting boundaries are placed symmetrically around the nullplane, at x w = ±0.5.Note that additional magnetic nullplanes appear in this interval for sufficiently large α.

B. Conventions and technicalities
a. ψ-regimes.In the analytic literature surrounding the resistive tearing instability (notably Refs.7 and 34), the mode is usually classified by the behaviour of the normalised, magnetic x-perturbation amplitude B1x , called ψ in Ref. 7 and subsequent literature, in a resistive layer [x 0 − δ , x 0 + δ ] around the magnetic nullplane at x = x 0 , where F(x 0 ) = 0.If ψ is approximately constant across this resistive layer, the mode is called a constant-ψ mode.In general, this corresponds to short wavelengths (large k).On the other hand, for longer wavelengths (small k), the variation in ψ throughout the resistive layer is not negligible.In this case, the tearing mode is called a nonconstant-ψ mode.This distinction is important because, analytically, the growth rate of the tearing mode scales differently with resistivity η for the two regimes.In the static (and small v 0 ) case, the growth rate scales as ∼ η 3/5 for constant-ψ modes, whereas for nonconstant-ψ modes, the growth rate scales as ∼ η 1/3 . 7,34,40o visualise how ψ changes with k, ψ is shown for the tearing mode of the static Harris sheet from Sec. II A 1 (B c = 1, a B = 1, v c = 0) for a selection of wavenumbers in Fig. 2(a).
Alongside ψ and the resistive layer halfwidth δ , the literature quantifies the matching quantity 7 In the analytic approach, ψ is obtained by solving the linearised, ideal MHD equations outside of the resistive layer, because resistivity is negligible there, and matching them with the resistive solution inside the layer.However, as an artefact of this approach, the resulting solution has a discontinuity in ψ ′ .To eliminate discontinuities from the calculations, the matching quantity ∆ ′ is introduced in the analytic approach.Consequently, this quantity appears in the analytic scaling laws, and imposes the condition δ |∆ ′ | < 1 for the constantψ approximation to hold. 7,34,40umerically, ψ ′ has no true discontinuity.Instead, there is a steep but smooth reversal across the magnetic nullplane.Hence, we here define ∆ ′ as the difference in ψ ′ /ψ between the extrema on either side of the nullplane, as seen for the static Harris sheet with k = 0.5 in Fig. 2 symmetric around the nullplane, we take with Ĵ′ z (x M ) = max( Ĵ′ z ) and Ĵ′ z (x m ) = min( Ĵ′ z ), to limit the impact of the discretised grid and the numerical differentiation of Ĵz .However, if background flow is included, all perturbations become complex.To maintain consistency with the static configurations, the complex factor is chosen in such a way that Im( B1x ) is positive and symmetric around the nullplane.For ψ ′ /ψ this factor is eliminated and the real part is used to define ∆ ′ , but this allows us to define δ based on the now odd function Re( Ĵ′ z ).To distinguish between analytic and numerical matching quantities and resistive layer halfwidths, we will refer to them with subscripts A and N, respectively.
b. Shear ratio.Similarly to F in Eq. ( 10), for the equilibrium flow we define the angle-modulated Alfvén Mach number (like Ref. 34) Here, c A indicates the dimensionless Alfvén speed c A = |B 0 |/ √ ρ 0 .Following Ref. 34, the expression where the prime denotes the derivative with respect to x, acts as a diagnostic parameter to quantify the relative strength of the flow shear compared to the magnetic shear.We will refer to R 0 as the shear ratio.c.Scaling laws.Now that we have introduced the ψregimes and shear ratio R 0 , we can differentiate between the various conditions under which the growth rate scaling laws were derived as a function of the resistivity η.The scaling laws are summarised in Table I as they were obtained by Chen and Morrison. 34BLE I. Analytic scaling laws of the tearing growth rate with the resistivity η. 34 Constant-ψ Nonconstant-ψ d. Relative quantities.Throughout this article, we compare the tearing growth rate under the influence of equilibrium flow to the equivalent configuration without flow.In these cases, we opt to use the relative growth rate γ, which we define as Hence, γ ranges from −1, which means the tearing instability is completely stabilised, to +∞.If γ = 0, background flow does not alter the growth rate.Similarly, we also define the relative numerical matching quantity e. Solvers.Throughout this work, two main solvers are used in LEGOLAS, as described in Ref. 53.The first one is the default solver, QR-cholesky, which results in the full spectrum, but is quite slow.The second one is the shift-invert Arnoldi solver, which only computes a selection of modes and is thus faster, but requires a target to converge around.This latter method is preferred for parameter studies once we have an approximation of the growth rate to act as the target.

III. RESULTS
In this work, we numerically investigate the flow-sheared resistive tearing mode in two different configurations.In Sec.
III A we first present visualisations of the perturbed magnetic field and flow in a Harris sheet (see Sec. II A 1) due to the linear tearing mode, both in the absence and presence of equilibrium flow.Then, we systematically vary the parameters appearing in the shear ratio R 0 to evaluate how the growth rate is affected.In Sec.III B we repeat this parameter study for the force-free magnetic field configuration (see Sec. II A 2), where we also vary the plasma-β .

A. Harris current sheet
For the Harris current sheet setup (Sec.II A 1), the effect of shear flow on the resistive tearing mode was probed in Ref. 28 using non-linear, incompressible MHD simulations by computing the reconnection rate for a selection of test cases.Here, we compute the linear growth rate using the compressible equations for various parameter combinations, initially to look at how the growth rate scales with resistivity, and later density.Assuming k = k y êy , the shear ratio reduces to Hence, we also verify that when the flow shear exceeds the magnetic shear, i.e.R 0 > 1, the tearing instability is fully stabilised, 34 by varying the velocity parameters.

Linear tearing of the Harris sheet
To begin, we look at the influence of the resistive tearing mode on the magnetic field and flow.To do so, we consider the Harris sheet presented in Sec.II A 1, both with and without the background flow, Eq. ( 7).In both cases, we set our parameters to k = 0.12 êy , ρ 0 = 1, B c = 1, and a B = 1, and use v c = 0.25 and a v = 0.75 when equilibrium flow is included.This is solved on a symmetric, accumulated grid as described in App.A with parameters p 1 = 0.2, p 2 = 0, p 3 = 0.01, and p 4 = 5, resulting in 327 grid points.In this configuration, the sheet is expected to tear up into magnetic islands, as represented schematically in Fig. 3.For the flowless case, Figs.4(a,b) show LEGOLAS's solutions for the Bx and By perturbation amplitudes from Eq. ( 5), respectively, where the arbitrary factor was chosen in such a way that the maximal density perturbation equals 1%, i.e. max(ρ 1 ) = 0.01 ρ 0 , at t = 0. Similarly, Figs.4(d,e) show the amplitudes in the case with background flow.In either case, B z is not perturbed.Note that in the flowless case one component is purely real and one is purely imaginary whilst the presence of the background shear flow makes both components fully complex.
Substituting these perturbation amplitudes in Eq. ( 5) for t = 0 and y ranging from −0.5λ to 0.5λ (≃ 26.18), with λ = 2π/|k| the wavelength, results in the linear perturbation of the magnetic field, presented as field lines and with a colour map of the B 1y -component, in Figs.4(c) and (f), for the case without and with flow, respectively.Note that since the amplitudes are normalised and can be multiplied with any complex factor, in principle, the (dimensionless) time t presented here cannot be linked to a physical time without context.Further note that the range of x-coordinates was limited to focus on the behaviour near the sheet.As expected from the schematic representation in Fig. 3, the visualisation in Fig. 4(c) reveals a magnetic island in the centre, with dipping in the previously straight magnetic field lines on either side.If the system were to evolve linearly, however, due to the sharp peaks in the B y perturbation amplitude on either side of the magnetic nullplane (see panel b), the magnetic field lines inside the island would dip progressively deeper inwards at the nullplane (because they cannot cross), until a magnetic field line meets itself again at the origin, dividing the island into two smaller islands (closed magnetic fieldlines) on either side of the nullplane.Since this is not observed in non-linear simulations, the time when this behaviour starts to develop in the linear solution marks an upper bound on the transition time from the linear to the non-linear regime.The inclusion of background shear flow does not dramatically alter this magnetic field structure, aside from introducing a slight shear deformation of the field lines across the nullplane.For the chosen parameters, this cannot be seen clearly at t = 0, shown in Fig. 4(f), but the effect becomes more apparent as the perturbation grows.
Similarly, Fig. 5 presents the vx and vy perturbation amplitudes from Eq. ( 5) for the flowless case in panels (a,b), and for the case with flow in panels (d,e).Again, the z-component, v z , is not perturbed, and the flowless case has a purely real and purely imaginary v 1 -component.Identically to the magnetic field perturbation, the velocity perturbation amplitude also becomes complex when background flow is added.Contrary to the magnetic field perturbation amplitude though, where the newly introduced real/imaginary parts are an order of magnitude smaller than the original part, the real and imaginary parts of the velocity perturbation amplitudes are of a comparable order of magnitude.
Once more, Figs.5(c,f) display a 2D visualisation of the magnetic field lines for the same y-interval and time as Figs.4(c,f).Here, however, the panels are coloured by the magnitude of the y-component of the velocity perturbation, v 1y .In addition, the flow patterns are further highlighted with stream vectors up to a certain magnitude |v|, to show variations in the regions of smaller speed.From the colour map in Fig. 5(c it is immediately clear that the velocity in the flowless case is much higher at the magnetic nullplane than further away from it.In fact, the velocity indicates an inflow towards the centre of the magnetic island along the nullplane.Away from the sheet's centre, the stream vectors cross the magnetic field lines, moving outward from the island's centre, though with a much smaller speed than the inflow speed.The addition of background flow in Fig. 5(f) changes the picture significantly.Whilst the flow speed is still much higher in the centre of the sheet, the velocity now roughly follows the magnetic field at the island edges, resulting in a flow within the island.Again, the islands are too small at t = 0 to clearly distinguish the flow pattern inside, but it becomes apparent as the perturbation grows.

Matching quantity and resistive layer halfwidth
Before turning to the growth rate scaling, we evaluate the numerical matching quantity and resistive layer halfwidth, how they compare to analytic values, and how they are affected by introducing flow.Once again, we first consider the static Harris sheet with k = 0.5 êy , ρ 0 = 1, B c = 1, and a B = 1 (v 0 = 0) in this section.For this static Harris sheet configuration, the analytic matching quantity is given exactly by 7,40 Hence, for the chosen parameters we find ∆ ′ A = 3.Additionally, for constant-ψ modes, the analytic resistive layer halfwidth is given by 40 which evaluates to δ A ≃ 4.13 × 10 −2 here, where B ′ 0 (x) is evaluated at the nullplane (x = 0).Note that this satisfies the constant-ψ condition δ A |∆ ′ A | ≃ 0.124 < 1. Numerically (grid parameters p 1 = 0.2, p 2 = 0, p 3 = 0.01, p 4 = 5), as explained in Sec.II B, we find ∆ ′ N ≃ 2.56 and δ N ≃ 3.5 × 10 −2 (note the use of ≃ due to the x-discretisation in LEGOLAS).Again, the constant-ψ condition is satisfied, Hence, though the analytic and numerical values for ∆ ′ and δ are not in perfect agreement, they are of the same order and reasonably close.Additionally, though the wavenumber k = 0.5 is neither small nor large considering that the configuration is tearing unstable for k ≲ 1, the constant-ψ approximation holds, and thus we expect to recover a growth rate scaling proportional to η 3/5 in the static case.
Analytically, Hofmann 30 showed that if the velocity v 0 is proportional to the Alfvén velocity (and thus the magnetic field B 0 ) everywhere, the matching quantity is unaltered from the static case.For the Harris sheet with velocity profile Eq. ( 7), v 0 ∝ B 0 is only true if a v = a B .However, as shown in Fig. 6(a), the ψ ′ /ψ ratio is not identical for the static and stationary (v c = 0.5, a v = 1) case.This is further highlighted by the difference between both cases, shown in Fig. 6(b).Consequently, neither is the numerical matching quantity.Now, we find a numerical matching quantity of ∆ ′ N ≃ 3.77 and a resistive layer halfwidth of δ N ≃ 5.5 × 10 −2 .Hence, both the numerical matching quantity and resistive layer halfwidth depend on the velocity, even if v 0 ∝ B 0 .However, δ N |∆ ′ N | ≃ 0.21 < 1 still holds.Since ψ ′ /ψ is identical in our incompressible approximation, which also eliminates the Joule heating, 15 this deviation from Hofmann's result is presumably due to the use of a constant resistivity rather than a convectively perturbed resistivity.
The discrepancy between analytic and numerical ∆ ′ and δ , as well as the difference between the static and stationary Harris sheet, raises the question how the numerical matching quantity and resistive layer halfwidth depend on the various parameters.In Fig. 6(c), both the analytic and numerical matching quantities are shown side by side as a function of k for the static Harris sheet.For large wavenumbers, i.e. in the constant-ψ regime, the analytic and numerical ∆ ′ are in good agreement, but for small wavenumbers, the analytic matching quantity diverges to +∞, whereas ∆ ′ N is observed to decrease again.The maximal value of ∆ ′ N in our set of discrete k-values was achieved for k = 0.09, with δ N |∆ ′ N | ≃ 0.62, indicating a strong numerical deviation from the analytic results in the traditional nonconstant-ψ regime.Simultaneously, no maximum is observed in δ N , which continues to increase as k decreases.Note that the steplike behaviour of δ N is due to the discretisation of the x-coordinate in LEGOLAS (also in panel f).This deviation from analytic results in the nonconstant-ψ regime is not surprising.As k decreases, ψ(0) approaches zero and ∆ ′ N approaches an x −1 scaling.Due to the definition of ∆ ′ N , the evaluation of ψ ′ /ψ occurs near the edge of the resistive layer, and thus ∆ ′ N ∼ δ −1 N .Hence, ∆ ′ N decreases as δ N increases for k → 0. Consequently, ∆ ′ N is not a proper numerical equivalent of ∆ ′ A in the nonconstant-ψ regime, and thus δ N |∆ ′ N | cannot be used to distinguish between constant-ψ and nonconstantψ modes.
Similarly to Fig. 6(c), Fig. 6(d) presents the dependence of ∆ ′ N and δ N on the magnetic reversal halfwidth a B .Here, the numerical matching quantity appears to follow the analytic scaling reasonably well, though it is consistently larger than ∆ ′ A for our choice of k = 0.5.The numerical resistive layer halfwidth δ N , on the other hand, is observed to be constant as a function of a B , which is surprising considering that Furth et al. 7 derived an a −1 B dependence for the width of the region of discontinuity.
Finally, again adding the velocity profile Eq. ( 7), Figs.6(e,f) show the dependence of ∆ ′ N and δ N on the transition halfwidth a v of the velocity profile for various flow speeds v c .As expected, the effect of flow on the matching quantity increases with the flow speed v c .For sufficiently high v c and small a v , the numerical matching quantity is observed to increase initially with a v , whereas the resistive layer halfwidth decreases.As a v increases further, ∆ ′ N reaches a maximum and decreases again steadily.For the resistive layer halfwidth, we observe a monotone decrease with a v , which is more pronounced for larger v c .

Growth rate scaling with resistivity
In the seminal work by Furth, Killeen, and Rosenbluth 7 the authors derive power laws for the scaling of the incompressible tearing growth rate as a function of the resistivity η, for small values of η.Here, we introduce compressibility and consider a wide range of resistivity values.However, contrary to their work, we assume a constant resistivity without convective perturbations, and include Joule heating.From their derivations they conclude that the growth rate scales as a power law in η, Im(ω) ∼ η p , and that p depends on whether or not ψ is approximately constant across the magnetic nullplane.This ψ-classification remains equally important when shear flow is added. 34 the last section, we have established that the tearing mode of the Harris sheet from Sec. II A 1 falls into the constant-ψ regime for the parameters k = 0.5 êy , ρ 0 = 1, B c = 1, a B = 1, a v = 1, and v c sufficiently small.Now, we vary the resistivity η to compare to the literature's corresponding scaling laws.In Fig. 7(a), the flowless case is compared to the case with flow profile Eq. ( 7) for a selection of R 0values, which are expected to scale differently based on the analytic power laws presented in Ref. 34. (Note that the no flow case is not clearly visible since it almost coincides with the R 0 = 0.1 case.)All LEGOLAS runs were performed using grid parameters p 1 = 0.2, p 2 = 0, p 3 = 0.01, and p 4 = 5 (327 grid points).The desired value of R 0 was obtained by setting Since the analytic scaling laws are derived under the assumption that there is a clear temporal separation between the resistive diffusion time τ R = η −1 (in our dimensionless  7,59 comparison to these power laws is only meaningful for η < 10 −2 in this case.Therefore, to be on the safe side, all fitted power laws from Table I were limited to η < 10 −3 in Fig. 7(a).In accordance with the work of Chen and Morrison, 34 the scaling of the static Harris sheet's growth rate is found to be close to η 3/5 , as expected in the constant-ψ regime.For R 0 = 0.1 ≪ 1, the growth rate is slightly larger than the static case, but the same scaling law seems to hold.As R 0 increases to R 0 = 0.5, the scaling law changes to η 1/2 , which is also in line with the literature for R 0 ≲ 1.However, for R 0 = 0.8, none of the analytic power laws are a good fit.At small η, this case is observed to scale as ∼ η 0.44 .Note that this power law lies between the constant-ψ η 1/2 and nonconstant-ψ η 1/3 scaling laws, though the product δ N |∆ ′ N | remains significantly smaller than 1, once again confirming it is not a good metric to classify ψ.In addition, the growth rate deviates strongly from a power law well before η reaches 10 −2 .

units) and the Alfvén time τ
As expected, aside from the growth rate dropoff in the R 0 = 0.8 case, the other cases also deviate from the power laws as η crosses the 10 −2 threshold.Note that growth rates above 10 −2 should be interpreted with care, since diffusion of the equilibrium was neglected and might affect the dynamics significantly.

Density variation
To demonstrate that the growth rate does not care about the specific density value, but only about the Alfvén speed, first consider the static Harris sheet with k = 0.5 êy , B c = 1, a B = 1, and varying density (grid parameters: p 1 = 0.2, p 2 = 0, p 3 = 0.01, p 4 = 5).This is shown in Fig. 7(b), where a fit to the data shows that the growth rate scales approximately as ρ −0.232 0 .Comparing this to the analytic (incompressible) growth rate scaling ρ −1/5 0 (see e.g.Ref. 7, Eq. ( 54)), whilst reasonably close, the deviation is significant.This may be due to compressibility, especially considering that the plasma-β goes to ∞ at the nullplane in this configuration.
Since both the Alfvén speed c A and sound speed c s are proportional to ρ −1/2 (due to our choice of T 0 ∝ ρ −1 0 profile), the question becomes how the growth rate scales with either speed.To answer this, we now consider a second static Harris sheet with the same parameters as before, except we set B c = √ ρ 0 and p c = 10 6 .These choices ensure that the Alfvén speed is constant when varying the density whilst T 0 ≃ p c /ρ 0 , and thus c s ∼ ρ −1/2 0 .In this case we find a constant growth rate, shown in Fig. 7(c).Hence, we conclude that the growth rate does not depend on the specific density or sound speed, but scales with the Alfvén speed.
The effect of the inclusion of a non-zero flow of fixed size v c = 0.1 (and various transition halfwidths) depends on the case.In the latter case, the v c /c A -ratio is constant, and the growth rate modification does not vary with ρ 0 .In the first case however, the flow modifies the tearing growth rate variation as a function of the density ρ 0 , or better, because the ratio v c /c A changes.In fact, the flow introduces a critical density above which the tearing mode is fully damped as v c /c A → 1.This is highlighted in Fig. 7(b) by the dotted line representing where v c = c A .However, from the figure it is clear that the tearing mode is not always damped exactly when v c reaches c A , but depends on the flow transition halfwidth a v .Therefore, the Alfvén speed does not necessarily act as a transition value in the equilibrium speed with respect to tearing suppression.Though the tearing instability may vanish at a density lower than where the equilibrium speed equals the Alfvén speed, it appears the Alfvén speed still imposes an upper limit on the density above which the tearing mode vanishes, from the sharp dropoff there in the a v = 1.25 case.

Velocity variation
Since both parameters of the velocity profile (the maximal speed v c and halfwidth a v ) appear to play an important role, we vary both parameters simultaneously to identify the regions of stabilisation and further destabilisation of the tearing instability.Here, the maximal speed v c is kept sub-Alfvénic (v c < 1) because the flow-induced KHI dominates in the super-Alfvénic regime. 30The result is shown in Fig. 8(a) for fixed parameters k = 0.5 êy , ρ 0 = 1, B 0 = 1, and a B = 1, and was obtained with the shift-invert method on an accumulated grid with parameters p 1 = 0.2, p 2 = 0, p 3 = 0.01, and p 4 = 5 (327 grid points).Since this method requires an initial guess, the first guess was obtained from a run with the QR-cholesky solver for v c = 10 −2 and a v = 1, which was used to compute the growth rate for all values of a v and v c = 10 −2 with the shift-invert method.This array of growth rates was subsequently used as the initial guesses for the next value of v c , and so on.Visually speaking, each value (except those in the left column) in Fig. 8(a) was computed via shift-invert by providing the value to its immediate left as the initial guess.Note that whilst the tearing mode is fully damped in the top right corner of panel (a), the system is still unstable because the KHI appears here before v c reaches the Alfvén speed.This was again checked with the QR-cholesky solver in this region of the parameter space.(For the few parameter combinations where the shift-invert method failed to converge, the growth rate was calculated using the QR-cholesky solver.)Since the analytic scalings laws typically feature a scaling with the matching quantity, Fig. 8(b) shows the relative numerical matching quantity ∆ ′ rel for the parameter combinations where the tearing growth rate exceeds 10 −4 , and grey elsewhere.
In this figure, a few things stand out.First of all, the resistive tearing instability is clearly stabilised at the R 0 = 1 (orange dashed) line for all considered parameter combinations, but is already fully suppressed before this line is reached.Since the analytic growth rate scales with ∆ ′ A 34 and ∆ ′ A is only identical to the static case if v 0 ∝ B 0 , 30 the instability is only expected to vanish exactly at R 0 = 1 if a v = a B , so this is not surprising.However, even for a v = a B (yellow dotted line), the instability is already fully suppressed before R 0 = 1 is reached.This may be due to the earlier observation (see Sec. III A 2) that ∆ ′ N does differ between static and stationary Harris sheets here, even if v 0 ∝ B 0 .
Secondly, contrary to the non-linear observation in Ref. 28 that there exists a single critical a v value ∼ 0.35 where the transition from stabilising to destabilising occurs, we here observe that this critical a v depends on the maximal speed v c .Furthermore, their critical value lies in our stabilising region of the parameter regime across all velocities.Note though that the simulations in Ref. 28 are incompressible, include a nonzero viscosity, and lack Joule heating, any of which may affect this result.
Finally, comparing Figs.8(a) and (b), there is a discrepancy between the variations in the growth rate and the numerical matching quantity.At smaller a v , the growth rate is slightly stabilised whilst the matching quantity is slightly increased.Therefore, since we only varied the flow parameters, the influence of shear flow on the growth rate is not fully encapsulated in the numerical matching quantity, and the growth rate does not simply scale with ∆ ′ N .

B. Force-free magnetic field
Now, we turn to the setup from Sec. II A 2, where a magnetic field of fixed magnitude varies its direction periodically throughout the plasma slab.In this particular case, the shear ratio becomes Whilst there are three parameters in this expression, our parametric study will only focus on the variation of the equilibrium density ρ c and velocity coefficient v c , for reasons which will become clear after a demonstration of the role of α.

Multiple tearing modes
Since the parameter α regulates how fast the magnetic field's direction varies along x, it also determines how many magnetic nullplanes the system has in a certain x-interval for a given wave vector.Consequently, the number of tearing modes supported by the configuration depends on α.Additionally, if the equilibrium velocity is described by an odd function, like the linear profile in Eqs.(11), the spectrum is symmetric with respect to the imaginary axis.This is illustrated in Figs.9(a-d), where we varied the parameter α for parameters ρ c = 1, β 0 = 0.15, v c = 0.15, and k = 1.5 êy at 251 grid points.For (a) α = π/2, the magnetic shear is insufficient to induce a tearing instability.Increasing α without introducing an additional nullplane results in one non-propagating tearing mode (i.e.purely imaginary), visualised in (b) for α = 4.73884 (slightly more than 3π/2).In the presence of three nullplanes for α = 5π/2, (c) shows a pair of forwardbackward propagating instabilities and one non-propagating one.Finally, (d) contains only two pairs of forward-backward propagating tearing pairs for α = 4.1π, despite the presence of 5 nullplanes in the domain.If the equilibrium flow is removed, all tearing modes become non-propagating, i.e. purely imaginary, as demonstrated in Fig. 10(a) for the case with α = 4.1π.
For the flowless case in Fig. 10(a), the magnetic field perturbation amplitudes of the unstable modes are shown in panels (b) through (e).All positions of magnetic nullplanes are marked with a dash-dotted line.At the darker-coloured nullplanes we observe a dip in B1x and a sharp transition in B1y , similar to Figs. 4(a,b), indicative of tearing at this nullplane.Note that the dominant instability features tearing behaviour at the 3 central nullplanes, whereas the secondary instability only tears up the central nullplane whilst the tertiary mode is tearing up the same nullplanes as the dominant instability except for the central nullplane.Surprisingly, the least unstable mode does not show any signs of tearing.Additionally, the outer nullplanes do not show strong signs of tearing, presumably due to their proximity to the perfectly conducting boundaries, which exert a stabilising influence.
With the addition of flow, however, the situation changes.Now, the spectrum in Fig. 9(d) no longer has a single dominant instability, but a dominant pair and less unstable pair of forward-backward propagating instabilities.Furthermore, their perturbation amplitudes are now fully complex.For the dominant, forward-propagating (Re(ω) > 0) instability, the B1 perturbation amplitudes are shown in Figs.9(e,g,i), and similarly, those of the less unstable, forward-propagating instability are shown in Figs.9(f,h,j).The backwardpropagating counterparts are not shown, but are the mirror image with respect to x = 0 of the forward-propagating perturbation amplitudes.The magnetic nullplanes are again indicated by dash-dotted lines.
Contrary to the static cases, all four modes appear to tear all three central nullplanes, coloured in blue, though lighter blue lines indicate multiplication with a complex factor is needed to highlight the tearing behaviour.Additionally, the dominant, forward-propagating mode's largest amplitudes variations (i.e.strongest tearing) occur at the positive intermediate and central nullplanes, with smaller variations at the negative intermediate nullplane.The less unstable, forwardpropagating mode, on the other hand, has its largest amplitudes variations at the central nullplane, with smaller variations at both intermediate nullplanes.Hence, by introducing flow the tearing behaviour of all modes is altered significantly.From now on, the value of α is set to α = 4.73884 (i.e. the value used in Fig. 9(b) and Ref. 58).This ensures that the magnetic field makes between one-half and a full rotation in the considered domain, x ∈ [−0.5, 0.5], resulting in a single nullplane at x = 0 and a single tearing mode.Since there is only one nullplane, no multitearing occurs. 60,61A detailed look at multitearing is beyond the scope of this paper.

Matching quantity and resistive layer halfwidth
Like in the previous case of the Harris sheet, we evaluate the matching quantity and resistive layer halfwidth for the forcefree magnetic field configuration before moving on to the scaling of the growth rate (for now, we set aside our conclusion from that section regarding ∆ ′ N ).Again, we choose the parameters k = 1.5 êy , ρ c = 1, α = 4.73884, and β 0 = 0.15.When velocity is included, we set v c = 0.25.All runs were performed at 251 grid points.Numerically, the static case yields ∆ ′ N ≃ 5.67 and δ N ≃ 1.4 × 10 −2 , whereas the stationary case yields ∆ ′ N ≃ 5.84 and δ N ≃ 1.4 × 10 −2 .Hence, the constant-ψ condition is satisfied in both cases, with δ N |∆ ′ N | ≃ 0.079 < 1  and δ N |∆ ′ N | ≃ 0.082 < 1 for the static and stationary case, respectively.Therefore, we again expect a growth rate scaling proportional to η 3/5 in the next section.
Again, we look at the influence of the flow profile on the ψ ′ /ψ-ratio.In Fig. 11(a), the ψ ′ /ψ-ratio is shown for the static and stationary case.The difference between both cases is shown in Fig. 11(b).Clearly, the ψ ′ /ψ-ratio does not change wildly with the addition of this flow profile.However, the difference between both cases is not negligible either, explaining the difference in numerical matching quantity calculated above.Now the question arises how strongly the numerical matching quantity is impacted by a variation in wavenumber or speed, and whether the growth rate scales with ∆ ′ N .In Fig. 11(c), the numerical matching quantity and resistive layer halfwidth are shown as a function of k for the static case.Similarly to the Harris sheet, the numerical matching quantity initially increases with k before decreasing again.Interestingly though, the resistive layer halfwidth is decreasing as k increases, but the product δ N |∆ ′ N | never exceeds ∼ 0.20, which would imply that the constant-ψ approximation is valid across all k-values.However, as we concluded in Sec.III A 2, the definition of ∆ ′ N makes it unfit to compare to ∆ ′ A in the nonconstant-ψ regime, and ensures that δ N |∆ ′ N | < 1 is an insufficient criterion to classify a perturbation as a constant-ψ mode.Additionally, though the growth rate variation with the wavenumber appears to follow a similar trend as the numerical matching quantity, it appears they are not directly proportional.This is further highlighted by Fig. 11(d), where the numerical matching quantity and resistive layer halfwidth are shown as a function of v c for the stationary case.Here, the numerical matching quantity is observed to increase with v c initially, whereas the growth rate starts to decline sooner as v c increases.Here, the resistive layer halfwidth is observed to be more or less constant.Consequently, we henceforth abandon the study of the numerical matching quantity and its relation to the growth rate for this configuration, and solely focus on the modification of the growth rate by the flow profile.

Growth rate scaling with resistivity
Once again, we now turn to a comparison of the growth rate scaling with analytic predictions.The parameters are the same as in the previous section, except that β takes on various values and v c = 0.15 when flow is included.As demonstrated, for this equilibrium we are dealing with a constant-ψ mode and thus expect a scaling of Im(ω) ∼ η 3/5 in the absence of flow.For the case with flow, we expect an identical scaling since R 0 ≃ 0.03 ≪ 1. 34 This scaling is expected to hold for a/η > 10 2 , since the boundary layer approach used in analytic works is no longer appropriate as the resistivity increases, 59 and thus, substituting a quarter period for the transition halfwidth a = π/2α, for η < 3.3 × 10 −3 .Hence, all comparisons were limited to η < 10 −3 .Indeed, initially, the growth rate scales as η 3/5 , as shown in Fig. 12(a) for 301 grid points, independently of β .However, as η approaches 10 −3 , the growth rate starts to deviate from this scaling, with smaller β -values deviating sooner, and eventually decreases again.Adding velocity to this variation in resistivity steepens the growth rate dropoff for lower β -values, as evidenced by Fig. 12(b), going as far as eliminating the instability entirely.Of course, care should again be taken in the interpretation of growth rates above η ∼ 10 −2 , since diffusion of the equilibrium was neglected.

Density variation
Adopting the same approach as for the Harris sheet configuration, we here show that the growth rate of the forcefree magnetic field configuration also varies with the Alfvén speed.To do so, we first consider the flowless case with parameters k = 1.5 êy and α = 4.73884, for various values of β and varying the density ρ 0 (301 grid points).As can be seen in Fig. 13(a), a similar growth rate scaling, ρ −0.254 0 , to the Harris sheet is recovered for high plasma-β .For this configuration the deviation from the analytic ρ −1/5 0 scaling is even greater, though an explanation might be that this is due to the incompressible approximation breaking down at high β .However, as β decreases, the growth rate deviates even more from a power law, especially for larger densities (smaller Alfvén speeds).This may be due to the inclusion of Joule heating, whose contribution to the energy equation is more significant for small β (T 0 ∝ β ).Now adding a velocity with v c = 0.15, behaviour similar to the Harris sheet is observed in Fig. 13(b).Again, deviations from the flowless scaling are observed as the density approaches the critical value where the maximal flow in the domain equals the Alfvén speed (indicated with a dotted line).Above this threshold, the tearing mode is heavily damped and the growth rate goes to zero.
To show once again that this scaling is due to the change in Alfvén velocity, consider the above case but with B 0 the profile in Eqs.(11) multiplied with √ ρ 0 , such that the Alfvén speed is constant.The resulting growth rates are shown in Figs.13(c) and (d), for the static and stationary cases, respectively.In both cases, the growth rate is constant for all densities, with higher β resulting in a higher growth rate.Hence, the growth rate does not depend on the specific density, but on the Alfvén speed.Since β appears to affect the growth rate, its role is further investigated in conjunction with the velocity.

Velocity variation
Despite the simple velocity profile, the (v c , β )-parameter space reveals surprising complexity.After assuming a constant density ρ 0 = 1, the present parametric survey varied v c and β simultaneously for wave vector k = 1.5 êy , where the v c parameter was limited to the interval [0, 2], such that the equilibrium velocity remains sub-Alfvénic (|v 0 | ≤ c A ) on the entire domain, and β -values from 10 −2 to 10 2 were studied.The results are shown in Fig. 14, where all runs in panel (a) were performed at 301 grid points, whereas 201 grid points were used in panel (b).
As pointed out in Ref. 30, the introduction of flow in a system that is unstable to the resistive tearing mode can either stabilise or further destabilise the plasma.This is also immediately clear from Fig. 14(b), where blue indicates a stabilised system and red a strong increase in tearing growth rate.Whilst the plasma is mostly destabilised further by the presence of flow for large β , as clearly evidenced by Fig. 14, the destabilising effect does not scale monotonically with the velocity coefficient.Rather, the maximal destabilisation appears at some intermediate value between small speeds and the Alfvén speed.For small to intermediate β (≲ 1), on the other hand, both stabilising and destabilising influences are observed in significant fractions of the velocity space, with the strongest stabilising effect occurring at the Alfvén speed.Additionally, more than one stabilising-destabilising transition is observed in panel (b) along the speed axis for small β (≪ 1).This dependence on the plasma-β , which already appeared in a less pronounced way in Sec.III B 3, is not surprising.The ion sound Larmor radius ρ s is known to affect the growth rate, 59,62 and β enters in the radius through our definition of T 0 .Nevertheless, any dependencies on ρ s were derived using an isothermal closure, contrary to the equations used here.Since it is clear that the plasma-β affects the role of flow significantly, a more in-depth analysis of the influence of β on tearing modes in general offers interesting perspectives for future work.

IV. CONCLUSION
In this work we studied the linear regime of the resistive tearing mode in a compressible plasma with Joule heating in two different configurations featuring shear flow: a Harris current sheet and a force-free magnetic field varying its direction periodically throughout a plasma slab.
First, we visualised the magnetic field lines and flow patterns of the linearly perturbed Harris sheet in 2D, both in the absence and presence of a background flow.In either case, the magnetic field lines were pinched together periodically along the sheet to reconnect and form magnetic islands, as observed in simulations.If this perturbation is allowed to evolve linearly, the formation of two smaller sub-islands is observed inside the island.Since this behaviour does not occur in nonlinear simulations, the time where this formation is initiated in the linear evolution places an upper limit on the transition time from the linear to the non-linear regime.
For the velocity, the largest perturbation occurs at the magnetic nullplane for both cases with and without equilibrium flow, with plasma leaving the pinched regions and streaming towards the magnetic islands.The plasma further away from the nullplane has an almost negligible velocity compared to the plasma at the nullplane.However, for the flowless equilibrium the plasma was observed to move away from the magnetic islands here, whereas the plasma rotates inside the inner islands if a background flow is included.
Next, we introduced a numerical equivalent of the matching quantity ∆ ′ and resistive layer halfwidth δ .However, for both configurations the numerical quantities deviated from the analytic predictions.In the nonconstant-ψ regime in particular, the numerical matching quantity fails to capture the behaviour of the analytic matching quantity, due to the lack of sharp transition in ψ ′ /ψ near the nullplane for wider resistive layers.In addition, no clear scaling of the growth rate with the numerical matching quantity was observed.Therefore, we advocate against the use of the matching quantity as the sole quantifier of flow's influence on the tearing growth rate.
As a consequence of this deviation from analytic predictions in the matching quantity in the nonconstant-ψ regime, the product δ N |∆ ′ N | is not a good indicator of the validity of the constant-ψ approximation.Nevertheless, the literature's scaling laws hold in the constant-ψ regime, as far as we have observed.Indeed, the growth rates of both configurations were found to scale as η 3/5 in the absence of flow, as expected for a constant-ψ mode.Subsequently, when flow was introduced, the transition to the constant-ψ η 1/2 scaling was also observed for considerable flow speeds.However, in the case of the Harris sheet, a different scaling was found for a case with even larger flow speed.Of course, the growth rate scaling with resistivity deviated from the analytic scaling for stronger resistivities, as is to be expected, especially in the presence of flow.In both cases, flow was observed to introduce a cutoff resistivity above which the tearing mode is damped.In the case of the Harris sheet, this cutoff even lay inside the regime where the analytic scaling law is still expected to hold, though this was for the case that did not follow the analytic scaling law anyway.
Afterwards, we showed that the growth rate scales with the Alfvén speed, though the scaling was found to be closer to c A .Additionally, a significant deviation from this scaling was observed for larger densities, i.e. smaller Alfvén speeds, in the force-free field case with small plasma-β .The importance of the plasma-β was then further highlighted in this configuration, where an intricate interplay between the plasma-β and flow speed was observed.For large plasma-β , the flow was observed to destabilise the plasma further, with the strongest destabilisation occurring at intermediate flow speeds.For small plasma-β though, both stabilising and destabilising influences were observed, with the strongest stabilisation occurring at the Alfvén speed.Additionally, more than one stabilising-destabilising transition was observed along the flow speed axis for small plasma-β .Presumably, the importance of the plasma-β in this study is due to the inclusion of Joule heating, which is absent in preceding work.
Finally, it was shown that the addition of flow to an equilibrium with multiple magnetic nullplanes and tearing modes can modify the tearing behaviour significantly.The tearing behaviour of all modes was altered, with all modes tearing all nullplanes, albeit to different degrees.This is in stark contrast to the static case, where only the dominant mode tore all nullplanes, and the secondary and tertiary modes only tore the central and non-central nullplanes, respectively.
Looking ahead, future work could focus more on the role of the plasma-β , especially at small values.Additionally, LEGOLAS could be used to incorporate viscosity, 15 or to investigate transitions in instability dominance, from resistive tearing to the Kelvin-Helmholtz instability, at near-Alfvénic speeds.Also the Hall field 14 and electron inertia effects 15 are prime candidates for further investigation.Of course, a similar study to this one can be performed for the cylindrical tearing mode, 63 to compare the effects of axial and azimuthal flow.Finally, due to the growth rate depending on the specific flow profile, LEGOLAS could also be employed as a computationally inexpensive, diagnostic tool for concrete configurations, particularly for experiments and for comparison of linear theory to non-linear simulations.

FIG. 3 .
FIG. 3. Schematic representation of magnetic field lines altered by the tearing instability.

FIG. 6 .
FIG. 6.(a) ψ ′ /ψ for k = 0.5 in a static and stationary Harris sheet.(b) Stationary minus static difference in ψ ′ /ψ from panel (a).(c) Static ∆ ′ and δ as a function of k = k y êy .(d) Static ∆ ′ and δ as a function of a B .(e) ∆ ′ N and (f) δ N for k = 0.5 in a Harris sheet with flow of various transition halfwidths.

FIG. 7 .
FIG. 7. Resistive tearing growth rate of a Harris sheet, Eq. (6), as a function of (a) the resistivity η; (b) the density ρ 0 ; for parameters k = 0.5 êy , ρ 0 = 1, B c = 1, a B = 1, and (a) a v = 1; (b) v c = 0.1, if flow was included.In (b), the density where the maximal velocity of the equilibrium configuration coincides with the Alfvén speed is indicated with a dotted line.(c) Growth rate scaling with ρ 0 for constant Alfvén speed (B c = √ ρ 0 ) and p c = 10 6 .

FIG. 8 .
FIG. 8. (a) Relative tearing growth rate γ and (b) relative numerical matching quantity ∆ ′ rel with respect to the static case for combinations of the maximal speed v c and flow transition halfwidth a v .The yellow dotted line indicates the magnetic field transition halfwidth a B and the orange dashed line represents R 0 = 1.

FIG. 11 .
FIG. 11.(a) ψ ′ /ψ for k = 1.5 in the static and stationary force-free B 0 -configuration.(b) Stationary minus static difference in ψ ′ /ψ from panel (a).(c) Static δ N , ∆ ′ N , and tearing growth rate as a function of k = k y êy .(d) Stationary δ N , ∆ ′ N , and tearing growth rate as a function of v c .

FIG. 13 .
FIG. 13.Resistive tearing growth rate of a force-free magnetic field, Eqs.(11), for k = 1.5 êy , as a function of ρ 0 for various plasma-β (a) without flow and (b) for v c = 0.15.The dotted vertical line indicates where the maximal equilibrium speed equals the Alfvén speed.(c) and (d) show the growth rate for the same parameters as in (a) and (b), respectively, but with B 0 scaled with √ ρ 0 such that the Alfvén speed is constant.