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 behavior 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én speed, the plasma-*β*, and the velocity parameters, revealing surprising nuance in whether the velocity acts stabilizing or destabilizing. 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 behavior 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 unraveling 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 initialization in current sheets is intrinsically linked to tearing instabilities, which occur due to magnetic shear, and create the necessary reconnection points (X-points). 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 generalized 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. Following this work, many more studies have honed in on the growth rate modification of resistive tearing in the Hall regime.^{10–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–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,22} Finally, 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 generalized 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.^{26} Hence, 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,29} In this endeavor, the influence of equilibrium flow on the resistive tearing mode has already been studied extensively using both analytic^{30–34} and numerical^{28,29,35–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.^{39} Furthermore, 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}

A 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–48} and non-linear simulations^{49–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 stabilizing 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 Legolas (Refs. 15, 52, and 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

*x*-dependence only) and assume Fourier solutions

*ρ*,

**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 = \rho T$ denotes the pressure, $ J = \u2207 \xd7 B$ the current density, and $ \eta = 10 \u2212 4$, unless specified otherwise, and $ \gamma = 5 / 3$ are the resistivity and adiabatic index, respectively. Finally, $ k = k y \u2009 e \u0302 y + k z \u2009 e \u0302 z$ and

*ω*are the wave vector and frequency. The resulting algebraic problem is solved for the frequency and corresponding perturbed quantities by Legolas.

### A. Configurations

Now consider a semi-infinite plasma confined in the *x* direction 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.

#### 1. Harris current sheet

^{54}i.e.,

*x*= 0, so does

**k**, which we choose to be along the

*y*axis, $ k = k y \u2009 e \u0302 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 = \xb1 15 \u2009 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 | \u2273 10 \u2009 a B$). Due to the sharp transitions in the equilibrium profiles near the origin and approximately constant behavior away from the center, the problem is solved using a non-uniform grid concentrated near the origin, as described in the Appendix.

#### 2. Force-free magnetic field

*α*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 / \rho c$. The equilibrium is visualized in Fig. 1(b).

For this equilibrium, we choose **k** proportional to $ e \u0302 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 = \xb1 0.5$. Note that additional magnetic nullplanes appear in this interval for sufficiently large *α*.

### B. Conventions and technicalities

*1. ψ*-regimes

In the analytic literature surrounding the resistive tearing instability (notably Refs. 7 and 34), the mode is usually classified by the behavior of the normalized, magnetic *x*-perturbation amplitude $ B \u0302 1 x$, called *ψ* in Ref. 7 and subsequent literature, in a resistive layer $ [ x 0 \u2212 \delta , x 0 + \delta ]$ 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 $ \u223c \eta 3 / 5$ for constant-*ψ* modes, whereas for nonconstant-*ψ* modes, the growth rate scales as $ \u223c \eta 1 / 3$.^{7,34,40}

To visualize 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).

*ψ*and the resistive layer halfwidth

*δ*, the literature quantifies the matching quantity

^{7}

^{,}

*ψ*is obtained by solving the linearized, 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 artifact of this approach, the resulting solution has a discontinuity in $ \psi \u2032$. To eliminate discontinuities from the calculations, the matching quantity $ \Delta \u2032$ is introduced in the analytic approach. Consequently, this quantity appears in the analytic scaling laws and imposes the condition $ \delta | \Delta \u2032 | < 1$ for the constant-

*ψ*approximation to hold.

^{7,34,40}

*k*= 0.5 in Fig. 2(b). As shown in Ref. 59, defining

*δ*as the distance from the nullplane to the nearest inflexion point of $ J \u0302 z$, i.e., $ J \u0302 z \u2033 ( \delta ) = 0$, as shown in Fig. 2(c), is consistent with boundary layer theory. Though $ J \u0302 z$ is symmetric around the nullplane, we take

*δ*based on the now odd function $ Re ( J \u0302 \u2032 z )$.

To distinguish between analytic and numerical matching quantities and resistive layer halfwidths, we will refer to them with subscripts $A$ and $N$, respectively.

#### 2. Shear ratio

*F*in Eq. (10), for the equilibrium flow, we define the angle-modulated Alfvén Mach number (like Ref. 34)

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

#### 3. 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 summarized in Table I as they were obtained by Chen and Morrison.^{34}

. | Constant-ψ
. | Nonconstant-ψ
. |
---|---|---|

$ R 0 \u226a 1$ | $ \gamma \u223c \eta 3 / 5$ | $ \gamma \u223c \eta 1 / 3$ |

$ R 0 \u2272 1$ | $ \gamma \u223c \eta 1 / 2$ | $ \gamma \u223c \eta 1 / 3$ |

$ R 0 > 1$ | Stabilized | Stabilized |

. | Constant-ψ
. | Nonconstant-ψ
. |
---|---|---|

$ R 0 \u226a 1$ | $ \gamma \u223c \eta 3 / 5$ | $ \gamma \u223c \eta 1 / 3$ |

$ R 0 \u2272 1$ | $ \gamma \u223c \eta 1 / 2$ | $ \gamma \u223c \eta 1 / 3$ |

$ R 0 > 1$ | Stabilized | Stabilized |

#### 4. Relative quantities

*γ*, which we define as

*γ*ranges from –1, which means the tearing instability is completely stabilized, to $ + \u221e$. If

*γ*= 0, background flow does not alter the growth rate. Similarly, we also define the relative numerical matching quantity

#### 5. 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 visualizations 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

^{34}by varying the velocity parameters.

#### 1. 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 \u2009 e \u0302 y , \u2009 \rho 0 = 1 , \u2009 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 the Appendix with parameters $ p 1 = 0.2 , \u2009 p 2 = 0 , \u2009 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) and 4(b) show Legolas's solutions for the $ B \u0302 x$ and $ B \u0302 y$ 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 ( \rho 1 ) = 0.01 \u2009 \rho 0$, at *t* = 0. Similarly, Figs. 4(d) and 4(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, while 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 $ \u2212 0.5 \lambda $ to $ 0.5 \lambda ( \u2243 26.18 )$, with $ \lambda = 2 \pi / | k |$ the wavelength, results in the linear perturbation of the magnetic field, presented as field lines and with a color map of the $ B 1 y$ component, in Figs. 4(c) and 4(f), for the case without and with flow, respectively. Note that since the amplitudes are normalized 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 behavior near the sheet. As expected from the schematic representation in Fig. 3, the visualization in Fig. 4(c) reveals a magnetic island in the center, 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 behavior 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 $ v \u0302 x$ and $ v \u0302 y$ perturbation amplitudes from Eq. (5) for the flowless case in panels (a) and (b) and for the case with flow in panels (d) and (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) and 5(f) display a 2D visualization of the magnetic field lines for the same *y*-interval and time as Figs. 4(c) and 4(f). Here, however, the panels are colored by the magnitude of the *y*-component of the velocity perturbation, $ v 1 y$. 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 color 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 toward the center of the magnetic island along the nullplane. Away from the sheet's center, the stream vectors cross the magnetic field lines, moving outward from the island's center, though with a much smaller speed than the inflow speed. The addition of background flow in Fig. 5(f) changes the picture significantly. While the flow speed is still much higher in the center 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.

#### 2. Matching quantity and resistive layer halfwidth

*a*= 1 ( $ v 0 = 0$) in this section. For this static Harris sheet configuration, the analytic matching quantity is given exactly by

_{B}^{7,40}

*ψ*modes, the analytic resistive layer halfwidth is given by

^{40}

*x*= 0). Note that this satisfies the constant-

*ψ*condition $ \delta A | \Delta A \u2032 | \u2243 0.124 < 1$.

Numerically (grid parameters $ p 1 = 0.2 , \u2009 p 2 = 0 , \u2009 p 3 = 0.01 , \u2009 p 4 = 5$), as explained in Sec. II B, we find $ \Delta N \u2032 \u2243 2.56$ and $ \delta N \u2243 3.5 \xd7 10 \u2212 2$ (note the use of $\u2243$ due to the *x*-discretization in Legolas). Again, the constant-*ψ* condition is satisfied, $ \delta N | \Delta N \u2032 | \u2243 0.090 < 1$. Hence, though the analytic and numerical values for $ \Delta \u2032$ 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 \u2272 1$, the constant-*ψ* approximation holds, and thus, we expect to recover a growth rate scaling proportional to $ \eta 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 \u221d B 0$ is only true if *a _{v}* =

*a*. However, as shown in Fig. 6(a), the $ \psi \u2032 / \psi $ ratio is not identical for the static and stationary ( $ v c = 0.5$,

_{B}*a*= 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 $ \Delta \u2032 N \u2243 3.77$ and a resistive layer halfwidth of $ \delta N \u2243 5.5 \xd7 10 \u2212 2$. Hence, both the numerical matching quantity and resistive layer halfwidth depend on the velocity, even if $ v 0 \u221d B 0$. However, $ \delta N | \Delta \u2032 N | \u2243 0.21 < 1$ still holds. Since $ \psi \u2032 / \psi $ is identical in our incompressible approximation, which also eliminates the Joule heating,

_{v}^{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 $ \Delta \u2032$ 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 $ \Delta \u2032$ are in good agreement, but for small wavenumbers, the analytic matching quantity diverges to $ + \u221e$, whereas $ \Delta \u2032 N$ is observed to decrease again. The maximal value of $ \Delta \u2032 N$ in our set of discrete *k*-values was achieved for *k* = 0.09, with $ \delta N | \Delta \u2032 N | \u2243 0.62$, indicating a strong numerical deviation from the analytic results in the traditional nonconstant-*ψ* regime. Simultaneously, no maximum is observed in $ \delta N$, which continues to increase as *k* decreases. Note that the steplike behavior of $ \delta N$ is due to the discretization 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, $ \psi ( 0 )$ approaches zero and $ \Delta \u2032 N$ approaches an $ x \u2212 1$ scaling. Due to the definition of $ \Delta \u2032 N$, the evaluation of $ \psi \u2032 / \psi $ occurs near the edge of the resistive layer, and thus, $ \Delta \u2032 N \u223c \delta N \u2212 1$. Hence, $ \Delta \u2032 N$ decreases as $ \delta N$ increases for $ k \u2192 0$. Consequently, $ \Delta \u2032 N$ is not a proper numerical equivalent of $ \Delta \u2032 A$ in the nonconstant-*ψ* regime, and thus, $ \delta N | \Delta \u2032 N |$ cannot be used to distinguish between constant-*ψ* and nonconstant-*ψ* modes.

Similarly to Fig. 6(c), Fig. 6(d) presents the dependence of $ \Delta \u2032 N$ and $ \delta 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 $ \Delta \u2032 A$ for our choice of

*k*= 0.5. The numerical resistive layer halfwidth $ \delta N$, on the other hand, is observed to be constant as a function of

*a*, which is surprising considering that Furth

_{B}*et al.*

^{7}derived an $ a B \u2212 1$ dependence for the width of the region of discontinuity.

Finally, again adding the velocity profile Eq. (7), Figs. 6(e) and 6(f) show the dependence of $ \Delta \u2032 N$ and $ \delta 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*, the numerical matching quantity is observed to increase initially with

_{v}*a*, whereas the resistive layer halfwidth decreases. As

_{v}*a*increases further, $ \Delta \u2032 N$ reaches a maximum and decreases again steadily. For the resistive layer halfwidth, we observe a monotone decrease with

_{v}*a*, which is more pronounced for larger $ v c$.

_{v}#### 3. 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 ( \omega ) \u223c \eta 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}

In Sec. III A 2, 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 \u2009 e \u0302 y , \u2009 \rho 0 = 1 , \u2009 B c = 1$, *a _{B}* = 1,

*a*= 1, and $ v c$ sufficiently small. Now, we vary the resistivity

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

_{0}values, 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 , \u2009 p 2 = 0 , \u2009 p 3 = 0.01$, and $ p 4 = 5$ (327 grid points). The desired value of

*R*

_{0}was obtained by setting $ v c = a v B c R 0 / a B$.

Since the analytic scaling laws are derived under the assumption that there is a clear temporal separation between the resistive diffusion time $ \tau R = \eta \u2212 1$ (in our dimensionless units) and the Alfvén time $ \tau A = a B \rho 0 \u2009 B c \u2212 1$,^{7,59} comparison to these power laws is only meaningful for $ \eta < 10 \u2212 2$ in this case. Therefore, to be on the safe side, all fitted power laws from Table I were limited to $ \eta < 10 \u2212 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 $ \eta 3 / 5$, as expected in the constant-*ψ* regime. For $ R 0 = 0.1 \u226a 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 $ \eta 1 / 2$, which is also in line with the literature for $ R 0 \u2272 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 $ \u223c \eta 0.44$. Note that this power law lies between the constant-*ψ* $ \eta 1 / 2$ and nonconstant-*ψ* $ \eta 1 / 3$ scaling laws, though the product $ \delta N | \Delta \u2032 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 \u2212 2$.

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 \u2212 2$ threshold. Note that growth rates above $ 10 \u2212 2$ should be interpreted with care, since diffusion of the equilibrium was neglected and might affect the dynamics significantly.

#### 4. 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 \u2009 e \u0302 y , \u2009 B c = 1$, *a _{B}* = 1, and varying density (grid parameters: $ p 1 = 0.2 , \u2009 p 2 = 0 , \u2009 p 3 = 0.01 , \u2009 p 4 = 5$). This is shown in Fig. 7(b), where a fit to the data shows that the growth rate scales approximately as $ \rho 0 \u2212 0.232$. Comparing this to the analytic (incompressible) growth rate scaling $ \rho 0 \u2212 1 / 5$ [see, e.g., Ref. 7, Eq. (54)], while reasonably close, the deviation is significant. This may be due to compressibility, especially considering that the plasma-

*β*goes to $\u221e$ at the nullplane in this configuration.

Since both the Alfvén speed $ c A$ and sound speed $ c s$ are proportional to $ \rho \u2212 1 / 2$ (due to our choice of $ T 0 \u221d \rho 0 \u2212 1$ 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 = \rho 0$ and $ p c = 10 6$. These choices ensure that the Alfvén speed is constant when varying the density while $ T 0 \u2243 p c / \rho 0$, and thus, $ c s \u223c \rho 0 \u2212 1 / 2$. 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 \u2192 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.

#### 5. 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 stabilization and further destabilization 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.

^{30}The result is shown in Fig. 8(a) for fixed parameters $ k = 0.5 \u2009 e \u0302 y , \u2009 \rho 0 = 1 , \u2009 B 0 = 1$, and

*a*= 1 and was obtained with the shift-invert method on an accumulated grid with parameters $ p 1 = 0.2 , \u2009 p 2 = 0 , \u2009 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 \u2212 2$ and

_{B}*a*= 1, which was used to compute the growth rate for all values of

_{v}*a*and $ v c = 10 \u2212 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 while 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 scaling laws typically feature a scaling with the matching quantity, Fig. 8(b) shows the relative numerical matching quantity $ \Delta \u2032 rel$ for the parameter combinations where the tearing growth rate exceeds $ 10 \u2212 4$, and gray elsewhere.

_{v}In this figure, a few things stand out. First of all, the resistive tearing instability is clearly stabilized 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 $ \Delta \u2032 A$^{34} and $ \Delta \u2032 A$ is only identical to the static case if $ v 0 \u221d B 0$,^{30} the instability is only expected to vanish exactly at $ R 0 = 1$ if *a _{v}* =

*a*, so this is not surprising. However, even for

_{B}*a*=

_{v}*a*(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 $ \Delta \u2032 N$ does differ between static and stationary Harris sheets here, even if $ v 0 \u221d B 0$.

_{B}Second, contrary to the non-linear observation in Ref. 28 that there exists a single critical *a _{v}* value $ \u223c 0.35$ where the transition from stabilizing to destabilizing occurs, we here observe that this critical

*a*depends on the maximal speed $ v c$. Furthermore, their critical value lies in our stabilizing region of the parameter regime across all velocities. Note though that the simulations in Ref. 28 are incompressible, include a non-zero viscosity, and lack Joule heating, any of which may affect this result.

_{v}Finally, comparing Figs. 8(a) and 8(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 stabilized while 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 $ \Delta \u2032 N$.

### B. Force-free magnetic field

While there are three parameters in this expression, our parametric study will only focus on the variation of the equilibrium density $ \rho c$ and velocity coefficient $ v c$, for reasons which will become clear after a demonstration of the role of *α*.

#### 1. 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)–9(d), where we varied the parameter *α* for parameters $ \rho c = 1 , \u2009 \beta 0 = 0.15 , \u2009 v c = 0.15$, and $ k = 1.5 \u2009 e \u0302 y$ at 251 grid points. For (a) $ \alpha = \pi / 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), visualized in (b) for $ \alpha = 4.738 \u2009 84$ (slightly more than $ 3 \pi / 2$). In the presence of three nullplanes for $ \alpha = 5 \pi / 2$, (c) shows a pair of forward–backward-propagating instabilities and one non-propagating one. Finally, (d) contains only two pairs of forward–backward-propagating tearing pairs for $ \alpha = 4.1 \pi $, despite the presence of five 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 $ \alpha = 4.1 \pi $.

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-colored nullplanes, we observe a dip in $ B \u0302 1 x$ and a sharp transition in $ B \u0302 1 y$, similar to Figs. 4(a) and 4(b), indicative of tearing at this nullplane. Note that the dominant instability features tearing behavior at the three central nullplanes, whereas the secondary instability only tears up the central nullplane, while 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 stabilizing 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 ( \omega ) > 0$] instability, the $ B \u0302 1$ perturbation amplitudes are shown in Figs. 9(e), 9(g), and 9(i), and similarly, those of the less unstable, forward-propagating instability are shown in Figs. 9(f), 9(h), and 9(j). The backward-propagating 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, colored in blue, though lighter blue lines indicate multiplication with a complex factor is needed to highlight the tearing behavior. 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, forward-propagating 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 behavior of all modes is altered significantly.

From now on, the value of *α* is set to $ \alpha = 4.738 \u2009 84$ (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 \u2208 [ \u2212 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,61} A detailed look at multitearing is beyond the scope of this paper.

#### 2. 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 force-free magnetic field configuration before moving on to the scaling of the growth rate (for now, we set aside our conclusion from that section regarding $ \Delta \u2032 N$). Again, we choose the parameters $ k = 1.5 \u2009 e \u0302 y , \u2009 \rho c = 1 , \u2009 \alpha = 4.738 \u2009 84$, and $ \beta 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 $ \Delta \u2032 N \u2243 5.67$ and $ \delta N \u2243 1.4 \xd7 10 \u2212 2$, whereas the stationary case yields $ \Delta \u2032 N \u2243 5.84$ and $ \delta N \u2243 1.4 \xd7 10 \u2212 2$. Hence, the constant-*ψ* condition is satisfied in both cases, with $ \delta N | \Delta \u2032 N | \u2243 0.079 < 1$ and $ \delta N | \Delta \u2032 N | \u2243 0.082 < 1$ for the static and stationary case, respectively. Therefore, we again expect a growth rate scaling proportional to $ \eta 3 / 5$ in Sec. III B 3.

Again, we look at the influence of the flow profile on the $ \psi \u2032 / \psi $ ratio. In Fig. 11(a), the $ \psi \u2032 / \psi $ ratio is shown for the static and stationary case. The difference between both cases is shown in Fig. 11(b). Clearly, the $ \psi \u2032 / \psi $ 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 $ \Delta \u2032 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 $ \delta N | \Delta \u2032 N |$ never exceeds $ \u223c 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 $ \Delta \u2032 N$ makes it unfit to compare to $ \Delta \u2032 A$ in the nonconstant-*ψ* regime, and ensures that $ \delta N | \Delta \u2032 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.

#### 3. 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 Sec. III B 2, 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 ( \omega ) \u223c \eta 3 / 5$ in the absence of flow. For the case with flow, we expect an identical scaling since $ R 0 \u2243 0.03 \u226a 1$.^{34} This scaling is expected to hold for $ a / \eta > 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 = \pi / 2 \alpha $, for $ \eta < 3.3 \xd7 10 \u2212 3$. Hence, all comparisons were limited to $ \eta < 10 \u2212 3$. Indeed, initially, the growth rate scales as $ \eta 3 / 5$, as shown in Fig. 12(a) for 301 grid points, independently of *β*. However, as *η* approaches $ 10 \u2212 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 $ \eta \u223c 10 \u2212 2$, since diffusion of the equilibrium was neglected.

#### 4. Density variation

Adopting the same approach as for the Harris sheet configuration, we here show that the growth rate of the force-free magnetic field configuration also varies with the Alfvén speed. To do so, we first consider the flowless case with parameters $ k = 1.5 \u2009 e \u0302 y$ and $ \alpha = 4.738 \u2009 84$, 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, $ \rho 0 \u2212 0.254$, to the Harris sheet is recovered for high plasma-*β*. For this configuration, the deviation from the analytic $ \rho 0 \u2212 1 / 5$ 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 \u221d \beta $).

Now adding a velocity with $ v c = 0.15$, behavior 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 $ \rho 0$, such that the Alfvén speed is constant. The resulting growth rates are shown in Figs. 13(c) and 13(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.

#### 5. Velocity variation

Despite the simple velocity profile, the $ ( v c , \beta )$-parameter space reveals surprising complexity. After assuming a constant density $ \rho 0 = 1$, the present parametric survey varied $ v c$ and *β* simultaneously for wave vector $ k = 1.5 \u2009 e \u0302 y$, where the $ v c$ parameter was limited to the interval $ [ 0 , 2 ]$, such that the equilibrium velocity remains sub-Alfvénic ( $ | v 0 | \u2264 c A$) on the entire domain, and *β*-values from $ 10 \u2212 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 stabilize or further destabilize the plasma. This is also immediately clear from Fig. 14(b), where blue indicates a stabilized system and red a strong increase in tearing growth rate. While the plasma is mostly destabilized further by the presence of flow for large *β*, as clearly evidenced by Fig. 14, the destabilizing effect does not scale monotonically with the velocity coefficient. Rather, the maximal destabilization appears at some intermediate value between small speeds and the Alfvén speed. For small to intermediate *β* ( $ \u2272 1$), on the other hand, both stabilizing and destabilizing influences are observed in significant fractions of the velocity space, with the strongest stabilizing effect occurring at the Alfvén speed. Additionally, more than one stabilizing–destabilizing transition is observed in panel (b) along the speed axis for small *β* ( $ \u226a 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 $ \rho 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 $ \rho 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 visualized 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 behavior does not occur in non-linear 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 toward 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 $ \Delta \u2032$ 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 behavior of the analytic matching quantity, due to the lack of sharp transition in $ \psi \u2032 / \psi $ 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 $ \delta N | \Delta \u2032 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 $ \eta 3 / 5$ in the absence of flow, as expected for a constant-*ψ* mode. Subsequently, when flow was introduced, the transition to the constant-*ψ* $ \eta 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.

Afterward, we showed that the growth rate scales with the Alfvén speed, though the scaling was found to be closer to $ c A 1 / 2$ than to the analytic $ c A 2 / 5$. 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 destabilize the plasma further, with the strongest destabilization occurring at intermediate flow speeds. For small plasma-*β* though, both stabilizing and destabilizing influences were observed, with the strongest stabilization occurring at the Alfvén speed. Additionally, more than one stabilizing–destabilizing 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 behavior significantly. The tearing behavior 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.

## ACKNOWLEDGMENTS

This work was supported by funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program, Grant Agreement No. 833251 PROMINENT ERC-ADG 2018. J.D.J. acknowledges further funding by the UK's Science and Technology Facilities Council (STFC) Consolidated Grant No. ST/W001195/1. R.K. is further supported by Internal Funds KU Leuven through the Project No. C14/19/089 TRACESpace and an FWO Project No. G0B4521N.

We thank the referees for their constructive and thought-provoking comments, which helped to improve the manuscript significantly.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Jordi De Jonghe:** Formal analysis (lead); Investigation (lead); Methodology (equal); Visualization (lead); Writing – original draft (lead); Writing – review & editing (lead). **Rony Keppens:** Funding acquisition (lead); Methodology (equal); Writing – original draft (supporting); Writing – review & editing (supporting).

## DATA AVAILABILITY

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

All functionalities to reproduce the results presented here are available in Legolas v2.1.1. The Legolas code is freely available under the GNU General Public License. For more information, visit https://legolas.science/.

### APPENDIX: ACCUMULATED GRID

Due to the heavily localized transitions in the equilibrium profiles of the Harris current sheet, an equidistant grid would require many more grid points than a centrally accumulated grid to properly resolve the region of steepest change in the middle. Therefore, this study opted for a grid constructed using the algorithm below^{64} for an interval $ x \u2208 [ a , b ]$ and function *f*.

Declare array $ auxGrid$ |

Declare array $ finalGrid$ |

$ auxGrid ( 1 ) \u2190 a$ |

$ i \u2190 1$ |

while $ auxGrid ( i ) < ( a + b ) / 2$ do |

$ auxGrid ( i + 1 ) \u2190 auxGrid ( i ) + f ( auxGrid ( i ) )$ |

$ i \u2190 i + 1$ |

end while |

$ finalGrid ( 1 ) \u2190 a$ |

$ finalGrid ( 2 i \u2212 1 ) \u2190 b$ |

$ \kappa \u2190 ( ( a + b ) / 2 \u2212 auxGrid ( i \u2212 1 ) ) / ( auxGrid ( i ) \u2212 auxGrid ( i \u2212 1 ) )$ |

for j from 1 to i − 1 do |

$ finalGrid ( j + 1 ) \u2190 auxGrid ( j ) + \kappa f ( auxGrid ( j ) )$ |

$ finalGrid ( 2 i \u2212 j \u2212 1 ) \u2190 a + b \u2212 finalGrid ( j + 1 )$ |

end for |

$ finalGrid ( i + 1 ) \u2190 ( a + b ) / 2$ |

Declare array $ auxGrid$ |

Declare array $ finalGrid$ |

$ auxGrid ( 1 ) \u2190 a$ |

$ i \u2190 1$ |

while $ auxGrid ( i ) < ( a + b ) / 2$ do |

$ auxGrid ( i + 1 ) \u2190 auxGrid ( i ) + f ( auxGrid ( i ) )$ |

$ i \u2190 i + 1$ |

end while |

$ finalGrid ( 1 ) \u2190 a$ |

$ finalGrid ( 2 i \u2212 1 ) \u2190 b$ |

$ \kappa \u2190 ( ( a + b ) / 2 \u2212 auxGrid ( i \u2212 1 ) ) / ( auxGrid ( i ) \u2212 auxGrid ( i \u2212 1 ) )$ |

for j from 1 to i − 1 do |

$ finalGrid ( j + 1 ) \u2190 auxGrid ( j ) + \kappa f ( auxGrid ( j ) )$ |

$ finalGrid ( 2 i \u2212 j \u2212 1 ) \u2190 a + b \u2212 finalGrid ( j + 1 )$ |

end for |

$ finalGrid ( i + 1 ) \u2190 ( a + b ) / 2$ |

## REFERENCES

*Electromagnetic Phenomena in Cosmical Physics*

*AAS-NASA Symposium on the Physics of Solar Flares*(NASA Special Publication,

*Magnetic Reconnection in Plasmas*

*Magnetohydrodynamics of Laboratory and Astrophysical Plasmas*

*Numerical Models for Differential Problems*