1. Introduction
The tearing mode is a fundamental instability of resistive plasma current sheets. It provides a major mechanism for the rapid disruption and dissipation of current sheets on magnetohydrodynamic (MHD) scales. Understanding the tearing instability is therefore crucial for explaining instabilities in Tokamak plasmas (Rutherford Reference Rutherford1973; White et al. Reference White, Monticello, Rosenbluth and Waddell1977; Connor et al. Reference Connor, Cowley, Hastie, Hender, Hood and Martin1988), as well as explosive phenomena in space and astrophysical plasmas, including coronal mass ejections (Lin et al. Reference Lin, Murphy, Shen, Raymond, Reeves, Zhong, Wu and Li2015), solar flares (Lu et al. Reference Lu, Feng, Warmuth, Veronig, Huang, Liu, Gan, Ning, Ying and Gao2022) and nanoflares (Jafari, Vishniac & Xu Reference Jafari, Vishniac and Xu2021).
The tearing mode was first analysed by Furth, Killeen & Rosenbluth (Reference Furth, Killeen and Rosenbluth1963) and Coppi et al. (Reference Coppi, Galvao, Pellat, Rosenbluth and Rutherford1976). In an infinitely long current sheet, the growth rate of the most unstable tearing mode,
$\gamma _m$
, scales as
$\gamma _m \tau \sim S^{-1/2}$
, where
$\tau = a/V_A$
is the Alfvén crossing time,
$a$
is the current sheet thickness and
$V_A$
is the upstream Alfvén speed. The Lundquist number is defined as
$S = aV_A/\eta$
, where
$\eta$
is the magnetic diffusivity. In most space and laboratory plasma environments, the particle collision frequency is low and
$\eta$
is consequently very small. The corresponding Lundquist number is therefore typically very large (Ji & Daughton Reference Ji and Daughton2011), implying an extremely small tearing-mode growth rate.
Over the past two decades, it has been recognised that the tearing growth rate in a macroscopic current sheet depends not only on the Lundquist number, but also on the current-sheet aspect ratio. For a two-dimensional current sheet of finite length
$L$
, the maximum growth rate expressed in terms of the macroscopic Alfvén time is
$\gamma _m \tau _L \sim S_L^{-1/2} (a/L)^{-3/2}$
, where
$S_L= LV_A/\eta$
and
$\tau _L = L/V_A$
. Thus, the tearing growth rate is controlled jointly by the Lundquist number and the current sheet aspect ratio (Shibata & Tanuma Reference Shibata and Tanuma2001; Loureiro, Schekochihin & Cowley Reference Loureiro, Schekochihin and Cowley2007; Uzdensky, Loureiro & Schekochihin Reference Uzdensky, Loureiro and Schekochihin2010).
In particular, as a current sheet thins and its aspect ratio approaches
$a/L \sim S_L^{-1/3}$
, the tearing mode growth rate becomes independent of
$S_L$
and is of order
$\tau _L^{-1}$
. This regime is known as ‘ideal’ tearing (Pucci & Velli Reference Pucci and Velli2013; Del Zanna et al. Reference Del Zanna, Landi, Papini, Pucci and Velli2016; Papini et al. Reference Papini, Landi and Del Zanna2018). Numerical simulations have confirmed that the tearing instability can rapidly disrupt a long, thin current sheet into a large number of plasmoids (Bhattacharjee et al. Reference Bhattacharjee, Huang, Yang and Rogers2009; Huang & Bhattacharjee Reference Huang and Bhattacharjee2013; Tenerani et al. Reference Tenerani, Velli, Rappazzo and Pucci2015; Shi et al. Reference Shi, Tenerani, Velli and Lu2019), enhancing the reconnection rate.
More recently, tearing-mode theory has been extended to incorporate a variety of physical effects, including Hall physics (Pucci, Velli & Tenerani Reference Pucci, Velli and Tenerani2017; Papini et al. Reference Papini, Landi and Zanna2019), guide fields (Shi et al. Reference Shi, Velli, Pucci, Tenerani and Innocenti2020), different equilibrium magnetic field configurations (Pucci et al. Reference Pucci, Velli, Tenerani and Del Sarto2018; Shi et al. Reference Shi, Artemyev, Velli and Tenerani2021), neutral particles (Pucci et al. Reference Pucci, Singh, Tenerani and Velli2020) and plasma jets (Shi Reference Shi2022). Many of these theoretical models adopt the incompressible approximation. Nevertheless, the influence of compressibility on the tearing instability has also been considered in several earlier studies. Numerical calculations by Tachi, Steinolfson & Van Hoven (Reference Tachi, Steinolfson and Van Hoven1985) examined the compressible tearing mode using a linearised radiative-MHD code and found that compressibility slightly reduces the tearing growth rate relative to the incompressible case. Recent work by De Jonghe & Keppens (Reference De Jonghe and Keppens2024) analysed compressible linear tearing modes using the LEGOLAS eigenvalue solver (De Jonghe et al. Reference De Jonghe, Claes and Keppens2022) and found that the tearing mode growth rate increases with plasma-
$\beta$
, implying that incompressible tearing mode has larger growth rate than the compressible tearing mode. We note, however, that Tachi et al. (Reference Tachi, Steinolfson and Van Hoven1985) considered only the
$\beta =0.01$
case, whereas De Jonghe & Keppens (Reference De Jonghe and Keppens2024) examined the
$\beta$
-dependence at only a single wavenumber.
Two-dimensional compressible-MHD simulations have also shown that compressibility can modify the critical Lundquist number for the onset of plasmoid instability (Ni et al. Reference Ni, Ziegler, Huang, Lin and Mei2012; Baty Reference Baty2014). In particular, these simulations generally find a lower critical Lundquist number at higher plasma-
$\beta$
, suggesting enhanced tearing instability at higher
$\beta$
, in qualitative agreement with the trend reported by De Jonghe & Keppens (Reference De Jonghe, Claes and Keppens2024).
Despite these previous efforts, a systematic investigation of the plasma-
$\beta$
dependence of the linear compressible tearing mode over a wide wavenumber range is still lacking. An independent and systematic quantitative analysis of the linear eigenvalue problem therefore remains valuable. In this study, we derive the linearised compressible-MHD equation set for a force-free current sheet and calculate the linear growth rate of tearing mode by solving the resulting boundary-value-problem. This approach provides an independent determination of the linear eigenvalues and eigenfunctions, and enables us to quantify systematically the effects of compressibility under the assumptions considered here. We further compare the compressible results with their incompressible counterparts, and characterise the dependence of the growth rate and perturbation structure on plasma-
$\beta$
. The remainder of the paper is organised as follows. In § 2, we derive the governing equation set and boundary conditions. In § 3, we present the numerical results. In § 4, we summarise this study and discuss the implications of our results.
2. Derivation of equations to solve
2.1. Background fields
We start from the compressible resistive-MHD equation set with an adiabatic closure:
\begin{align} \frac {\partial \rho }{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }\left ( \rho \boldsymbol{U} \right ) = 0, \end{align}
\begin{align} \rho \left ( \frac {\partial \boldsymbol{U}}{\partial t} + \boldsymbol{U} \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{U} \right ) = - \boldsymbol{\nabla }\left ( P + \frac {B^2}{2 \mu _0} \right ) + \frac {1}{\mu _0} \boldsymbol{B} \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{B}, \end{align}
\begin{align} \frac {\partial \boldsymbol{B}}{\partial t} = \boldsymbol{\nabla }\times \left ( \boldsymbol{U} \times \boldsymbol{B} \right ) + \eta {\nabla} ^2 \boldsymbol{B} ,\end{align}
\begin{align} \frac {\partial P}{\partial t} = -\boldsymbol{U} \boldsymbol{\cdot }\boldsymbol{\nabla }P - \kappa P \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{U} ,\end{align}
where
$\rho$
,
$\boldsymbol{U}$
,
$\boldsymbol{B}$
and
$P$
denote density, velocity, magnetic field and thermal pressure, respectively;
$\mu _0$
is the vacuum magnetic permeability and
$\kappa$
is the adiabatic index.
Coordinate system adopted in this study.
We note that, in the solar corona and lower solar atmosphere, various thermodynamic processes, such as thermal conduction and radiative cooling, may become important and thus invalidate the adiabatic closure adopted here. A detailed investigation of these effects, however, is beyond the scope of the present study.
As illustrated in figure 1, we assume that the
$z$
-axis is aligned with the wave vector (
$\boldsymbol{k}$
) direction and that the background fields vary only along the
$x$
-direction. Here, we consider a force-free current sheet in which the background density
$\rho _0$
and pressure
$P_0$
are uniform, the background velocity vanishes everywhere, and the background magnetic field takes the form
\begin{align} \boldsymbol{B_0}(x) = B_0 \, \tanh {\left ( x/a \right )} \hat {e}_z + B_0 \, \mathrm{sech} \left ( x/a \right ) \hat {e}_y, \end{align}
so that
$\left | \boldsymbol{B_0} \right | = \text{Const.}$
everywhere. We note that the background field is in exact equilibrium only when
$\eta =0$
. For weakly collisional plasmas,
$\eta \rightarrow 0$
and the tearing mode growth time scale is much shorter than the diffusion time scale of the background field. We can therefore treat the background field as stationary. Strictly, the stationary-background approximation requires the tearing growth time to be shorter than the resistive diffusion time of the equilibrium. This condition becomes marginal or invalid at the smallest Lundquist numbers considered here. Results in this regime should therefore be regarded as properties of the corresponding stationary eigenvalue problem rather than of a self-consistently evolving resistive current sheet.
2.2. Equations for perturbations
We denote the perturbations in density, velocity and magnetic field by
$\rho _1$
,
$\boldsymbol{u}$
and
$\boldsymbol{b}$
, respectively. The pressure perturbation is simply
\begin{align} p_1 = C_s^2 \rho _1 ,\end{align}
where
$C_S = \sqrt {\kappa P_0 /\rho _0}$
is the sound speed. Linearisation and normalisation equation (2.1) gives
\begin{align} \gamma\kern-2pt \rho = - \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}, \end{align}
\begin{align} \gamma\kern-1pt \boldsymbol{u} = - \beta \boldsymbol{\nabla }\rho -\boldsymbol{\nabla }\left ( \boldsymbol{B} \boldsymbol{\cdot }\boldsymbol{b}\right ) + \boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{b} + \boldsymbol{b} \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{B}, \end{align}
\begin{align} \gamma\kern-1pt \boldsymbol{b} = -\boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{B} + \boldsymbol{B} \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u} - \boldsymbol{B} \left ( \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} \right ) + \frac {1}{S} {\nabla} ^2 \boldsymbol{b}, \end{align}
where
$\rho = \rho _1/\rho _0$
. The magnetic field is expressed in Alfvén-speed unit through
$B \rightarrow B/\sqrt {\mu _0 \rho _0}$
, all velocities are normalised by the upstream Alfvén speed
$V_A = B_0 / \sqrt {\mu _0 \rho _0}$
and all lengths are normalised by the current sheet thickness
$a$
. Here,
$\beta = C_s^2 / V_A^2$
is the plasma beta, which is
$\kappa /2$
times the conventional definition
$P_0/(B_0^2/2\mu _0)$
. Additionally,
$S = a V_A / \eta$
is the Lundquist number.
For the infinitely long current sheet in the
$z$
direction considered here, the perturbations take the form
\begin{align} \left ( \begin{array}{c} \rho (x,z,t) \\[2pt] \boldsymbol{u}(x,z,t) \\[2pt] \boldsymbol{b}(x,z,t) \end{array} \right ) = \left ( \begin{array}{c} \rho (x) \\[2pt] \boldsymbol{u}(x) \\[2pt] \boldsymbol{b}(x) \end{array} \right ) e^{\gamma t + i kz}. \end{align}
Substituting these expressions into (2.4), we obtain the following ordinary-differential-equation set:
\begin{align} \gamma\kern-1pt b_x = - {k\kern-1.3pt B}_z u_x + \frac {1}{S} \left ( b_x^{\prime \prime } - k^2 b_x \right ), \end{align}
\begin{align} \gamma\kern-1pt b_y = - B_y^\prime u_x + kB_z u_y - B_y \big(u_x^\prime + k u_z \big) + \frac {1}{S} \left ( b_y^{\prime \prime } - k^2 b_y \right ) ,\end{align}
\begin{align} \gamma\kern-1pt u_x = \frac {\beta }{\gamma }\left ( u_x^{\prime \prime } + k u_z^\prime \right ) - \left ( B_z \frac {b_x^\prime }{k} + B_y b_y \right )^\prime + k\kern-1.3pt B_z b_x, \end{align}
\begin{align} \gamma\kern-1pt u_y = - k\kern-1.3pt B_z b_y + B_y^\prime b_x, \end{align}
\begin{align} \gamma\kern-1pt u_z = - \frac {k\beta }{\gamma } \left ( u_x^\prime + k u_z \right ) +k\kern-1.3pt B_y b_y + B_z^\prime b_x ,\end{align}
where the prime denotes
${\rm d}/{\rm d}x$
. We have eliminated
$\rho$
using the continuity equation and eliminated
$b_z$
using the divergence-free condition
$\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{b} = 0$
. We also adopt the mappings:
$i\kern-1pt b_x \rightarrow b_x$
,
$i\kern-1pt u_y \rightarrow u_y$
and
$i\kern-1pt u_z \rightarrow u_z$
, so that the equation set becomes purely real. Accordingly, we restrict our analysis to the purely growing tearing branch and therefore take
$\gamma$
to be real, as in the incompressible case. We note that (2.6e
) can be used to eliminate the term
${\beta }/{\gamma }( u_x^{\prime \prime } + k u_z^\prime )$
in (2.6c).
In the limit
$\beta \rightarrow +\infty$
, (2.6e
) reduces to
$u_x^\prime + k u_z = 0$
, i.e. the incompressible condition
$\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} = 0$
. In this case, (2.6c
) can be simplified to
\begin{align} \frac {\gamma }{k} \left ( u_x^{\prime \prime } - k^2 u_x \right ) = B_z \left ( b_x^{\prime \prime } - k^2 b_x \right ) - B_z^{\prime \prime } b_x, \end{align}
which, together with (2.6a
), forms the two-equation set for the classical incompressible tearing instability (Pucci & Velli 2013). In the limit
$\beta \rightarrow 0$
, the order of the equation set is reduced because
$u_x^{\prime \prime }$
terms vanish.
2.3. Boundary condition
To numerically solve (2.6) and determine the eigenvalue
$\gamma$
, we need to specify the boundary conditions of the system. Far from the current sheet, the background magnetic field becomes nearly-uniform, with
$B_z \approx \text{Const.}$
and
$B_y \approx 0$
. The perturbation amplitudes should therefore decay with increasing distance from the centre of the current sheet. Assuming uniform background fields and writing the perturbed fields in the form
$f(x) = \bar {f} e^{\alpha x}$
, (2.6) gives
\begin{align} \left [ \begin{array}{c@{\quad}c} k\kern-1.3pt B & \gamma - \dfrac {1}{S} ( \alpha ^2 - k^2 ) \\[15pt] \gamma - \dfrac {1 + k F}{\gamma } \beta \alpha ^2 & \dfrac {B}{k} ( \alpha ^2 - k^2 ) \end{array} \right ] \left ( \begin{array}{c} \bar {u}_x \\[5pt] \bar {b}_x \end{array} \right ) = 0 ,\end{align}
where
$B$
is the asymptotic value of
$B_z$
and should be
$\pm 1$
,
\begin{align} F = - \frac {k\beta }{\gamma ^2 + k^2 \beta }, \end{align}
such that
$u_z = F u_x^\prime$
(2.6e
). Requiring the determinant of the matrix on the left-hand side to vanish gives the equation for
$\alpha (k,\gamma |\beta )$
:
\begin{align} \begin{aligned} p_4 \alpha ^4 + p_2 \alpha ^2 + p_0 = 0 \end{aligned} \end{align}
with
\begin{align} p_4 & = \frac {C \beta } {S \gamma }, \nonumber\\[3pt] p_2 & = - \left ( B^2 + \frac {\gamma }{S} \right ) - C \beta \left ( 1 + \frac {k^2}{S \gamma } \right ), \nonumber\\[6pt] p_0 & = \gamma ^2 + k^2 \left (B^2 + \frac {\gamma }{S} \right ), \end{align}
where
\begin{align} C = 1 + kF = \frac {\gamma ^2}{\gamma ^2 + k^2 \beta }. \end{align}
In the classical incompressible case,
$\alpha = \pm k$
, with the sign of
$\alpha$
chosen to ensure that the perturbations decay away from the current sheet. Equation (2.10) gives two solutions for
$\alpha ^2$
, one of which is close to
$k^2$
; we use this solution to specify the boundary condition.
3. Results
We use the boundary-value-problem solver implemented in the SciPy package (Virtanen et al. Reference Virtanen2020) to solve (2.6). The solver employs a fourth-order collocation algorithm with a damped Newton method, and can simultaneously solve for unknown parameters (
$\gamma$
in this case). The computational domain is set to
$x/a\in [-15,15]$
. The solver adaptively adjusts both the locations and the number of grid points, typically between 1000 and 4000 points, to achieve convergence for a prescribed tolerance. For a system
$y^\prime = f(x,y)$
, the tolerance is defined through the normalised residual criterion:
\begin{align} \varepsilon \geqslant \mathrm{r.m.s.}\left ( \frac {y^\prime - f(x,y)}{1 + |f(x,y)|} \right ). \end{align}
We adopt
$\varepsilon = 10^{-3}$
throughout this study. To assess the sensitivity of the results to this choice, we also performed calculations with
$\varepsilon = 5 \times 10^{-4}$
and
$10^{-2}$
for
$S=10^6$
and found that the resulting growth rates were insensitive to the tolerance over this range.
Figure 2 shows an example of the solved eigenfunctions
$u_x$
,
$b_x$
and
$\rho$
for
$S=10^4$
and
$k=0.127$
. Panel (b) is a zoom-in view of panel (a), as indicated by the grey shading. In each panel, we show three curves corresponding to three values of
$\beta$
: the incompressible limit (
$\beta =+\infty$
),
$\beta =1$
and
$\beta =0$
. The velocity and magnetic field perturbations are only weakly affected by variations in
$\beta$
. In contrast, the density perturbation depends strongly on
$\beta$
. As
$\beta$
decreases towards zero, a strong density perturbation develops within the inner layer. From the ranges of the vertical axes, one can see that, in the limit
$\beta \rightarrow 0$
, the amplitude of
$\rho$
becomes significantly larger than those of
$u_x$
and
$b_x$
, which is associated with the sharp variation of
$u_x$
across the inner layer.
Solved eigenfunctions for
$S=10^4$
and
$k=0.127$
(the most unstable mode for
$S=10^4$
). The right plot is a zoom-in view of the region marked by grey shading in the left plot. In each plot, the panels from top to bottom show
$u_x$
,
$b_x$
and
$\rho$
, respectively. Black lines correspond to the incompressible case (
$\beta \rightarrow +\infty$
), blue dashed lines to
$\beta =1$
and orange dotted lines to
$\beta =0$
.
In figure 3, we show the calculated dispersion relation
$\gamma -k$
for
$S=10^2$
(panel a) and
$S=10^4$
(panels b and c). Here panel (c) is a zoom-in view of the region enclosed by the black box in panel (b). In each panel, different curves correspond to different values of
$\beta$
: from dark to light,
$\beta$
increases from zero to infinity. In general, the growth rate increases as
$\beta$
decreases and the increase of
$\gamma$
is more pronounced for
$S=10^2$
than for
$S=10^4$
. For each combination
$(S,\beta )$
, we calculate the maximum growth rate
$\gamma _m$
and the corresponding wavenumber
$k_m$
. In figure 4, we plot
$\gamma _m-S$
(panel a) and
$k_m-S$
(panel b) for different values of
$\beta$
. For reference, we plot
$\gamma \propto S^{-1/2}$
and
$k \propto S^{-1/4}$
, which are the asymptotic scaling relations for the incompressible tearing mode, as black dashed lines. Clearly, for large
$S$
(
$S \geqslant 10^4$
), all curves converge, indicating that
$\beta$
has only a negligible influence on the maximum tearing mode growth rate. For smaller
$S$
, the curves diverge and the increase in
$\gamma _m$
associated with decreasing
$\beta$
becomes more significant. Here,
$k_m$
also varies with
$\beta$
: smaller
$\beta$
corresponds to smaller
$k_m$
, i.e. longer wavelength.
Dispersion relation
$\gamma -k$
for different
$\beta$
: (a)
$S=10^2$
; (b)
$S=10^4$
; (c) zoom-in of the region marked by the black box in panel (b).
Figure 3. Long description
Panel A: A line graph shows the dispersion relation for S = 10^2. The x-axis is labeled k and ranges from 0.10 to 0.50. The y-axis is labeled γ and ranges from 0.0325 to 0.0475. Different lines represent various values of β, including 0.00, 0.10, 1.00, 10.00, and +∞. Panel B: A line graph shows the dispersion relation for S = 10^4. The x-axis is labeled k and ranges from 0.00 to 0.30. The y-axis is labeled γ and ranges from 0.0002 to 0.0060. Different lines represent various values of β, including 0.00, 0.10, 1.00, 10.00, and +∞. Panel C: A line graph shows a zoom-in of the region marked by the black box in panel (b) for S = 10^4. The x-axis is labeled k and ranges from 0.10 to 0.15. The y-axis is labeled γ and ranges from 0.00596 to 0.00606. Different lines represent various values of β, including 0.00, 0.10, 1.00, 10.00, and +∞.
(a) Maximum growth rate as a function of
$S$
for different
$\beta$
. (b) Wavenumber for maximum growth rate as a function of
$S$
for different
$\beta$
.
We quantify the relative increase in
$\gamma _m$
as
$\beta$
decreases from
$+\infty$
to zero, i.e.
\begin{align} \varDelta (S) = \frac {\gamma _m(S,\beta =0) - \gamma _m(S,\beta =+\infty )}{\gamma _m(S,\beta =+\infty )}, \end{align}
and show the result in figure 5. Consistent with figure 4, the relative increase in
$\gamma _m$
is negligible for
$S\gt 10^4$
, remaining below 1 %, but becomes more significant at smaller
$S$
: for
$S=10^2$
,
$\varDelta$
is approximately 4 %, while for
$S=10$
,
$\varDelta$
can exceed 10 %.
Relative increase of
$\gamma _m$
as
$\beta$
decreases from
$+\infty$
to zero (3.2), as a function of
$S$
.
As already mentioned, the density perturbation can be very strong in the low-
$\beta$
regime. In figure 6, we show the ratio of the maximum amplitudes (i.e. the peaks of the eigenfunctions) of
$\rho$
and
$b_x$
for the most unstable mode as a function of
$S$
. The density perturbation decreases with increasing
$S$
, although the decrease is relatively weak. In contrast, the density perturbation is highly sensitive to
$\beta$
, with smaller
$\beta$
producing stronger density perturbations. Therefore, in the compressible tearing mode, strong density perturbations can develop near the centre of the current sheet. These density perturbations may influence the transition from the linear to the nonlinear stage of the instability.
Maximum amplitude of density fluctuation divided by maximum amplitude of
$b_x$
for the most unstable mode as a function of
$S$
. Different curves correspond to different
$\beta$
.
Figure 6. Long description
A line graph with the horizontal axis labeled S and the vertical axis labeled max(|rho|)/max(|bx|). The graph includes four data lines representing different values of beta: 0.0, 0.1, 1.0, and 10.0. Each line shows the relationship between S and the ratio of max(|rho|) to max(|bx|). The line for beta equals 0.0 starts at a higher value and decreases as S increases. The line for beta equals 0.1 follows a similar trend but starts at a lower value. The line for beta equals 1.0 starts even lower and shows a slight decrease. The line for beta equals 10.0 starts at the lowest value and remains relatively flat as S increases.
4. Conclusion
In this study, we derived the governing equations for the tearing mode from the fully compressible resistive-MHD equations and computed its linear growth rate in a Harris-type force-free current sheet using a boundary-value-problem solver. Our results show that compressibility generally increases the tearing mode growth rate. This increase is weak at large Lundquist numbers (
$S \gt 10^4$
), but becomes non-negligible at small Lundquist numbers (
$S\lt 10^2$
). Nevertheless, even though the tearing mode growth rate depends only weakly on plasma-
$\beta$
at large
$S$
, the associated density perturbation can be strong and may exceed the magnetic field perturbation in relative amplitude as
$\beta$
decreases. Thus, in low-
$\beta$
current sheets, such as those in the lower solar corona, the growth of tearing instability can generate strong density perturbations near the centre of the current sheet. These strong compressive fluctuations may influence the subsequent nonlinear evolution and thermodynamics of the current sheet.
The increase of tearing mode growth rate at small Lundquist numbers suggests, within the assumptions of the present linear model, that compressibility may lower the critical Lundquist number for the triggering of plasmoid instability relative to the incompressible case. This trend is opposite to that reported in previous compressible-MHD simulations (Ni et al. Reference Ni, Ziegler, Huang, Lin and Mei2012; Baty Reference Baty2014), which found a lower critical Lundquist number at higher plasma-
$\beta$
. We emphasise that these simulations do not employ the same equilibrium configuration as that considered here. In particular, thermal-pressure gradients contribute to maintaining force balance in those models, resulting in a non-uniform plasma-
$\beta$
across the current sheet, whereas the present force-free equilibrium has uniform thermal pressure. This difference in equilibrium structure may contribute to the contrasting plasma-
$\beta$
dependence.
Previous linear studies also do not provide a direct one-to-one comparison with the present model. Tachi et al. (Reference Tachi, Steinolfson and Van Hoven1985) found that compressibility slightly reduces the tearing-mode growth rate, but their calculation employed a radiative-MHD model that includes radiative cooling, Joule heating and temperature-dependent resistivity, rather than the adiabatic closure adopted here. De Jonghe & Keppens (Reference De Jonghe and Keppens2024) found a plasma-
$\beta$
dependence of the tearing-mode growth rate using the LEGOLAS eigenvalue solver, but their force-free equilibrium differs from ours and the
$\beta$
dependence was examined only at a fixed wavenumber, rather than by identifying the most unstable mode over the full dispersion relation for each
$S$
and
$\beta$
. These differences in equilibrium configuration, thermodynamic treatment and numerical methodology make a direct quantitative comparison difficult.
We therefore do not regard the discrepancy with previous results as fully resolved. A definitive comparison would require calculations using similar equilibrium configuration and physical assumptions. In particular, determining the linear tearing-mode growth rate for equilibria with non-uniform thermal pressure and conducting compressible-MHD simulations initialised with the same force-free equilibrium considered here would help clarify the origin of the discrepancy. These investigations, however, are beyond the scope of the present study and are left for future work.
Acknowledgements
Editor Francesco Califano thanks the referees for their advice in evaluating this article.
Funding
This work is supported by NSF SHINE #2548299 and NASA ECIP #80NSSC26K0321.
Declaration of interests
The authors report no conflict of interest.
Use of artificial intelligence (AI) tools
ChatGPT 5.6 Sol was used solely for language editing and polishing of the manuscript.
S=104
k=0.127
S=104
ux
bx
ρ
β→+∞
β=1
β=0
γ−k
β
S=102
S=104
S
β
S
β
γm
β
+∞
S
bx
S
β
Open access
