Flute and kink instabilities in a dynamically twisted flux tube with anisotropic plasma viscosity

ABSTRACT

Magnetic flux tubes such as those in the solar corona are subject to a number of instabilities. Important among them is the kink instability that plays a central part in the nanoflare theory of coronal heating, and for this reason in numerical simulations, it is usually induced by tightly controlled perturbations and studied in isolation. In contrast, we find that fluting modes of instability are readily excited when disturbances are introduced in our magnetohydrodynamic flux tube simulations by dynamic twisting of the flow at the boundaries. We also find that the flute instability, which has been theorized but rarely observed in the coronal context, is strongly enhanced when plasma viscosity is assumed anisotropic. We proceed to investigate the co-existence and competition between flute and kink instabilities for a range of values of the resistivity and of the parameters of the anisotropic and isotropic models of viscosity. We conclude that while the flute instability cannot prevent the kink from ultimately dominating, it can significantly delay its development especially at strong viscous anisotropy induced by intense magnetic fields.

instabilities, MHD, plasmas, Sun: corona, Sun: magnetic fields

1 INTRODUCTION

The helical kink instability is a form of ideal magnetohydrodynamic (MHD) instability that occurs in highly twisted magnetic flux tubes such as those making up much of the solar corona (Reale 2014) and has been well studied in the coronal context (Hood & Priest 1979; Hood, Browning & Van der Linden 2009; Browning & Van der Linden 2003; Török & Kliem 2003; Török, Kliem & Titov 2004; Török & Kliem 2005; Bareford & Hood 2015; Quinn, MacTaggart & Simitev 2020c). Given its energetic non-linear development, it is considered a potential mechanism for heating the solar corona through the theory of nanoflares (Klimchuk 2006; Browning 1991) and a key mechanism in the production of solar flares (Hood & Priest 1979). Our previous work investigated a twisted magnetic flux tube already linearly unstable to the helical kink instability, focussing specifically on the effect of anisotropic viscosity on the non-linear dynamics (Quinn et al. 2020c). There and in most other investigations of the kink instability e.g. that of Hood et al. (2009), a perturbation is applied to an already significantly twisted flux tube. An alternative way to excite the instability (and the way employed here) is to start with an initially straight field and apply twisting motions at the boundaries to form a twisted flux tube which eventually becomes unstable. This kind of dynamic excitation of the kink instability represents more closely the actual evolution of magnetic flux tubes and the associated instabilities in the solar corona. In our simulations, the dynamic twisting of the flux tube reveals an additional instability, the flute instability, which has been theorized, for example by Priest (2013). While oscillations resembling flute perturbations have been found in simulations of coronal loops (Terradas, Magyar & Van Doorsselaere 2018), to our knowledge, this is the first time the flute instability has been investigated computationally in a coronal context.

The flute instability arises in magnetized plasmas where the plasma pressure gradient is oriented in the same direction as the field line curvature, that is the pressure and magnetic tension forces compete. This is similar to the competition between pressure and gravitational forces which gives rise to the Rayleigh–Taylor instability (RTI). In MHD terminology, the RTI is a typical example of an ideal interchange instability, where magnetic field lines are minimally bent and are, instead, exchanged during the evolution of the instability. The ideal flute instability is another example of an ideal interchange instability but confined to a cylindrical geometry, the term ‘flute instability’ referring to its likeness to a fluted column. In a twisted flux tube like a simple, unbraided coronal loop, the magnetic curvature is always directed towards the axis so the tube may be unstable to fluting when the pressure decreases outwards from a high-pressure core. Such a pressure distribution is generated in the flux tubes studied here as a result of the driving. The appearance of the flute instability is illustrated by, for example, the pressure contours in Fig. 5, where the perturbations follow the pitch of the twisted field.

Figure 1.

Radial velocity profile u_r(r) and acceleration profile u_t(t) of the driver (13) for parameters u₀ = 0.15, r_d = 5, and t_r = 2.

Figure 2.

Kinetic energy (a,b) and cumulative heating (c) as a function of time showing the development and measured growth rates γ and λ of the flute and kink instabilities, respectively. Resistivity value is η = 10⁻⁴ and Fig. 2(b) is an enlarged version of Fig. 2(a). The cumulative heating 〈Q_*〉, where * is either ν for viscous heating or η for Ohmic heating, is the respective heating term integrated both over space and from the initial moment up to the moment t in time. The viscous heating associated with the flute instability (that generated before t ≈ 28) is negligible compared to that associated with the kink instability (generated after t ≈ 28). While the isotropic model permits greater viscous heating (on the order of two orders of magnitude), the anisotropic model enhances Ohmic heating.

Figure 3.

Gradients in the current density generate gas pressure gradients through Ohmic heating. Axial current density (a), Ohmic heating (b), and pressure (c) as functions of the radial distance from the tube axis. All plots are from anisotropic cases when t = 20 through the plane z = 0. Note the sign of the axial current density j_z has been flipped for comparison and the Ohmic heating is given by Q_η = ηj². The gas pressure profile of an additional test-case where η = 0 is also shown. Line types are indicated in the legend.

Figure 4.

Gas pressure profiles during the development of the flute and kink instabilities. Shown are density plots of gas pressure at z = 0 with η = 10⁻⁴ and for the anisotropic viscosity model. Note the difference in colour scale in Fig. 4(c). The development in the isotropic case is similar.

Figure 5.

Simultaneous development of flute and kink instabilities in the isotropic and anisotropic cases illustrated by field lines and gas pressure contours. Field lines and contours of gas pressure (where p = 0.3) are plotted at t = 28. Also shown is the velocity driver |$u_r(\sqrt{x^2+y^2})$| at z = 2. The flute instability grows in both cases, though faster in the anisotropic case. The initial stages of the kink instability can also be observed in the field lines of the isotropic case in Fig. 5(a).

In other solar contexts, interchange instabilities can be found in the form of ballooning modes in arcades (Hood 1986), as the instability which forms tubes of specific size in the photosphere (Bunte 1993), and in the buoyancy of flux tubes (Schuessler 1984). However, the flute instability specifically is more commonly studied in fusion contexts (Wesson 1978; Mikhailovskii 1998; Zheng 2015). In fusion, the focus is generally on understanding how a particular plasma device may be stabilized to the instability in particular geometries such as that of the mirror machine (Jungwirth & Seidl 1965) or in toroidal geometries such as the tokamak (Shafranov 1968). The resistive flute instability (also known as the resistive interchange instability) can be excited even when the ideal flute instability is stabilized. As a result, this has been given significantly more attention (Johnson & Greene 1967; Correa-Restrepo 1983). While this body of research is useful and applicable in solar contexts, it is mostly limited to the study of the stability and linear development of the flute instability, the non-linear development being of secondary importance in the investigation of fusion devices. More detailed investigations of its non-linear development is required to understand its importance in the context of coronal dynamics and coronal heating. The development of the flute instability and its interaction with the simultaneously growing kink instability is the main focus of this work and the experiments described here represent an initial exploration into the non-linear flute instability in the solar corona.

In addition to our main goal, of particular interest here is the effect of anisotropic plasma viscosity, which in the following is found to strongly influence the growth of the flute instability. It is well known that viscosity in magnetized plasmas (such as those which make up the solar corona) is anisotropic and strongly dependent on the strength and direction of the local magnetic field (Hollweg 1986, 1985; Braginskii 1965). To take this into account, MacTaggart, Vergori & Quinn (2017) developed a phenomenological model of anisotropic viscosity that captures the main physics of viscosity in the solar corona as outlined in the analysis of Braginskii (1965), namely parallel viscosity in regions of strong field strength and isotropic viscosity in regions of very weak or zero field strength. For brevity, we will refer to this model of viscosity as ‘the switching model’. In Quinn et al. (2020c) and Quinn, MacTaggart & Simitev (2021), we implemented the switching model as a module for the widely used general MHD code Lare3d (Arber et al. 2001), and demonstrated significant effects of anisotropic viscosity on the development of the non-linear MHD kink instability and the Kelvin–Helmholtz instability. More generally, the interest in anisotropic viscosity stems from the open question of which heating mechanism (viscous or Ohmic) is dominant in the solar corona (Klimchuk 2006), an important facet of solving the coronal heating problem. Using scaling laws, it has been suggested that viscous heating (generated through anisotropic viscosity) can dwarf that of Ohmic heating (Litvinenko 2005; Craig & Litvinenko 2009). However, due to computational and observational limitations, this cannot be directly tested, and so the influence of other factors such as small-scale instabilities and turbulence is relatively unknown (Klimchuk 2006). In addition to directly heating the plasma, viscosity plays a part in the damping of instabilities and waves (Ruderman et al. 2000). It is this effect, we are most interested in here, and it shall be reported that the use of anisotropic viscosity permits the growth of the flute instability, which is otherwise strongly damped by isotropic viscosity.

The value of plasma resistivity also affects the development of the flute instability because pressure gradients generated through Ohmic heating substantially contribute to its growth. Ideally, our simulations would be performed using a realistic coronal resistivity values of approximately 10⁻⁸ (Craig & Litvinenko 2009), however, due to the dissipative nature of numerical schemes (particularly when using shock capturing techniques) this is computationally infeasible. To overcome this limitation, we perform and compare simulations at two computationally accessible resistivity values, 10⁻³ and 10⁻⁴, in an attempt to extrapolate results towards more realistic values. This comparison runs as an additional theme of the paper, if not a primary aim.

Our article is organised as follows. Section 2 introduces the flute instability and recalls relevant background on its linear stability analysis. Section 3 describes the governing equations, the coronal loop model, and the numerical parameters used. Section 4 presents the overall development of the flute instability before proceeding to compare simulations in the cases of various viscous anisotropy and various resistivity values. Here, results are organized by resistivity values as this allows to contrast isotropic and anisotropic viscosity cases more directly. Section 5 discusses the limitations of the simulations, with suggestions for future work, and Section 6 presents our conclusions in the wider context of coronal heating.

2 THE FLUTE INSTABILITY

In general, the stability of a cylindrical twisted magnetic flux tube is analysed using perturbations of the form

$$\begin{eqnarray*} \xi (r, \theta , z) = \xi (r) {\rm e}^{i(m\theta + kz -\omega t)}, \end{eqnarray*}$$

(1)

where ω is the oscillation frequency in time t, m, and k are the wavenumbers in the azimuthal and axial directions, θ and z, respectively, and r is the radial coordinate in cylindrical polars. The helical kink instability occurs for perturbations where m = 1, k ≠ 0 and is the only instability of this form which is a body instability; that is, it moves the entire body of the flux tube. Perturbations where m > 1 are termed flute or interchange instabilities.

When the magnetic field is sheared, as in a twisted magnetic flux tube, an interchange instability (such as the flute instability) is confined to a surface where the peaks and troughs follow the shear of the field. That is, the instability is confined to the surface where the perturbation wavevector (0, m/r, k) is perpendicular to the direction of the field, known as the ‘resonance surface’. In an axisymmetric twisted flux tube the resonance surface is located at a radius r specified by

$$\begin{eqnarray*} \frac{m}{r} B_{\theta }(r) + kB_z(r) \approx 0. \end{eqnarray*}$$

(2)

The stability of an infinite cylindrical flux tube to perturbations of the form (1) is given by the classical Suydam’s criterion (Suydam 1958)

$$\begin{eqnarray*} \frac{B_z^2 S^2}{4} + 2 r p^{\prime } \gt 0, \end{eqnarray*}$$

(3)

where S = rq′/q is a measure of the shear, q = 2|$\pi$|rB_z/LB_θ is the safety factor for a flux tube of length L and a prime denotes differentiation with respect to r (Mikhailovskii 1998). This applies to both flute and kink instabilities, although many additional effects such as line-tying are not incorporated into the corresponding linear analysis. The effect of line-tying on the kink instability is investigated by Hood & Priest (1979). Where (3) is not satisfied, the flux tube may be unstable to perturbations of the form (1). When m > 1, the perturbations remain local to resonant surfaces given by (2). When Suydam’s criterion is satisfied and the flux tube is linearly stable, it may still be unstable to non-local perturbations, where the shear and pressure are small enough that interchange perturbations do not need to satisfy (2). Additionally, the inclusion of resistivity generally reduces the stabilizing effect of the shear, permitting growth of a resistive interchange mode, albeit at a significantly slower rate than that of the ideal instability (Mikhailovskii 1998). At the values of resistivity studied here, the resistive growth rate is expected to be approximately two orders of magnitude less than the corresponding ideal rate. Furthermore, it will be found that the ideal linear analysis of Mikhailovskii (1998) is sufficient for understanding the flute instabilities investigated here since the associated flux tubes substantially fail the criterion (3). For these reasons, we consider only the ideal flute instability.

While Suydam’s condition gives an indication of the stability of a flux tube to a given perturbation, the linear growth rate of the ideal flute instability γ, defined as the imaginary part of ω in ansatz (1), can be determined via a stability analysis analogous to that of the Rayleigh–Taylor instability (Goldston 2020). The fastest growing mode in the r-direction is found to be the longest wavelength mode, while the fastest growing modes in the θ-direction are found to be those with the shortest wavelengths, i.e. large values of m in the notation of equation (1). In particular, for all modes with wavelengths in the θ-direction that are shorter than both the pressure-gradient scale-length and the radial height of the plasma, the growth rate λ tends to the limit

$$\begin{eqnarray*} \gamma ^2 = \frac{2|\nabla p|}{\rho R_c}, \end{eqnarray*}$$

(4)

where R_c is the radius of curvature of the magnetic field. We find that this expression gives a good estimate of the growth rate of the flute instability in our numerical simulations even at moderate values of m. Equation (4) only applies when the pressure gradient and radius of curvature vector are in the same direction; that is, the plasma is constrained by a concave magnetic field such that the pressure forces and magnetic tension forces are in competition. In a cylindrical, twisted flux tube, the field is always concave towards the central axis of the tube, so any inwardly directed pressure gradient is potentially unstable to fluting.

Throughout this paper, the twisted flux tube generated by the drivers has a pressure profile which is approximately axisymmetric, and independent of z away from the boundaries at z = ±2, and has a negative gradient, hence |∇p| may be written as −dp/dr. Similarly, away from the boundaries, the magnetic field has a negligible r component and little dependence on θ and z, allowing the field to be approximated as |${\boldsymbol B} = (0, B_{\theta }(r), B_z(r))^{\text{T}}$|⁠, in cylindrical coordinates (r, θ, z). For a twisted field of this form, the radius of curvature is given by

$$\begin{eqnarray*} R_c = \frac{1}{|({\boldsymbol b}\cdot \nabla) {\boldsymbol b}|} = \frac{r}{b_{\theta }^2}, \end{eqnarray*}$$

(5)

where |${\boldsymbol b} = {\boldsymbol B}/|{\boldsymbol B}|$| is the unit vector in the direction of the magnetic field and b_θ is the component of |${\boldsymbol b}$| in the azimuthal direction. These approximations allow the growth rate to be written as

$$\begin{eqnarray*} \gamma _{\text{ideal}}^2 = \frac{-2p^{\prime }}{\rho R_c}. \end{eqnarray*}$$

(6)

This approximation for the growth rate continues to hold while the flux tube remains relatively axisymmetric, that is while the kink instability remains in its linear phase.

The stability criterion (3) and the linear growth rate approximation (6) are useful only as a guide and for approximate analysis of the numerical simulations presented in this work. The precise form of the equilibrium state and the perturbations needed for the validity of equations (3) and (6) were used by Quinn et al. (2020c). In contrast, in the experiments reported in the following the system is driven and instabilities occur spontaneously due to random perturbations. As a result of the driving, the flux tube is also not in static equilibrium initially.

3 MATHEMATICAL FORMULATION AND NUMERICAL SET-UP

We consider the MHD equations for the density ρ, plasma velocity |${\boldsymbol u}$|⁠, pressure p, magnetic field |${\boldsymbol B}$|⁠, and internal energy ε, in their non-dimensionalized visco-resistive form

$$\begin{eqnarray*} \frac{D\rho }{Dt} = - \rho {\boldsymbol \nabla } \cdot {\boldsymbol u}, \end{eqnarray*}$$

(7a)

$$\begin{eqnarray*} \rho \frac{D{\boldsymbol u}}{Dt} = -{\boldsymbol \nabla } p + {\boldsymbol \jmath } \times {\boldsymbol B} + {\boldsymbol \nabla } \cdot {\mathsf \sigma }, \end{eqnarray*}$$

(7b)

$$\begin{eqnarray*} \frac{D{\boldsymbol B}}{Dt} = ({\boldsymbol B} \cdot {\boldsymbol \nabla }){\boldsymbol u} - ({\boldsymbol \nabla } \cdot {\boldsymbol u}){\boldsymbol B} + \eta \nabla ^2 {\boldsymbol B}, \end{eqnarray*}$$

(7c)

$$\begin{eqnarray*} \rho \frac{D\varepsilon }{Dt} = -p {\boldsymbol \nabla } \cdot {\boldsymbol u} + {Q}_{\nu } + {Q}_{\eta }, \end{eqnarray*}$$

(7d)

where η is the non-dimensionalized resistivity, |${\boldsymbol \jmath } = \nabla \times {\boldsymbol B}$| is the current density, and the terms |${Q}_{\nu } = {\mathsf \sigma } : {\boldsymbol \nabla }{\boldsymbol u}$| and |${Q}_{\eta } = \eta | {\boldsymbol \jmath } |^2$| are viscous heating and Ohmic heating, respectively. The system is closed by the inclusion of the equation of state for an ideal gas

$$\begin{eqnarray*} \varepsilon = \frac{p}{\rho (\gamma - 1)}, \end{eqnarray*}$$

(8)

with the specific heat ratio is given by γ = 5/3.

Two different models for the viscosity stress tensor |${\mathsf \sigma }$| will be compared and contrasted in this study. The first model is the conventional isotropic (or Newtonian) viscosity stress tensor used in the vast majority of the existing literature, so that,

$$\begin{eqnarray*} {\mathsf \sigma = }{\mathsf \sigma }_{\text{iso}} = \nu \mathsf{W}, \end{eqnarray*}$$

(9)

where ν is the viscous transport parameter, generally referred to as the viscosity,

$$\begin{eqnarray*} {\mathsf W} = \nabla {\boldsymbol u} + (\nabla {\boldsymbol u})^T - \frac{2}{3}(\nabla \cdot {\boldsymbol u}){\mathsf I}, \end{eqnarray*}$$

(10)

is the rate of strain tensor, and |${\boldsymbol I}$| is the 3 × 3 identity tensor. The second model, which is the one of actual interest, is the anisotropic viscosity stress tensor given by

$$\begin{eqnarray*} {\mathsf \sigma } = {\mathsf \sigma }_\text{aniso} = \nu \left[\frac{3}{2}({\boldsymbol W}{\boldsymbol b}\cdot {\boldsymbol b}) \left({\boldsymbol b} \otimes {\boldsymbol b} - \frac{1}{3}{\boldsymbol I} \right)\right], \end{eqnarray*}$$

(11)

where |${\boldsymbol b}$| is the unit vector in the direction of the magnetic field.

Expression (11) is identical to the strong field approximation of the general anisotropic viscosity tensor derived by Braginskii (1965). Expressions (9) and (11) arise as asymptotic limits of the more general switching model used in our earlier works (MacTaggart et al. 2017; Quinn et al. 2020c, 2021) which includes both isotropic and anisotropic contributions and can switch gradually between them depending on the strength of the magnetic field at a given spatio-temporal location. For example, in the vicinity of a null point where the magnetic field becomes weak the isotropic viscosity contribution becomes dominant in the switching model. Switching between the two limit cases is not relevant in this study, where the variations in the magnetic field are not significantly large.

The non-dimensionalization of equations (7) is identical to that used in the earlier works by Quinn et al. (2020c, 2021). A typical magnetic field strength B₀, density ρ₀, and length scale L₀ are chosen and the other variables non-dimensionalized appropriately. Velocity and time are non-dimensionalized using the Alfvén speed |$u_A = B_0 / \sqrt{\rho _0 \mu _0}$| and Alfvén crossing time t_A = L₀/u_A, respectively. Temperature is non-dimensionalized via |$T_0 = u_A^2 \bar{m} / k_B$|⁠, where k_B is the Boltzmann constant and |$\bar{m}$| is the average mass of ions, here taken to be |$\bar{m} = 1.2m_p$| (a mass typical for the solar corona) where m_p is the proton mass. Dimensional quantities can be recovered by multiplying the non-dimensional variables by their respective reference value (e.g. |${\boldsymbol B}_{\dim } = B_0 {\boldsymbol B}$|⁠). The reference values used here are B₀ = 5 × 10⁻³ T, L₀ = 1 Mm, and ρ₀ = 1.67 × 10⁻¹² kg m⁻³, giving reference values for the Alfvén speed u_A = 3.45 M ms⁻¹, Alfvén time t_A = 0.29 s, and temperature T₀ = 1.73 × 10⁹ K.

The following initial and boundary conditions are used to form a magnetic flux tube and excite instabilities by dynamic twisting. The magnetic field is prescribed as initially straight and uniform,

$$\begin{eqnarray*} {\boldsymbol B} = (0, 0, 1)^{\text{T}}, \end{eqnarray*}$$

(12)

in a cube of size [ − 2, 2]³, with further test simulations run using an elongated domain of size 4 × 4 × 20. Initially, the velocity is set everywhere to |${\boldsymbol u} = {\boldsymbol 0}$|⁠, the density to ρ = 1, and the internal energy to ε = 8.67 × 10⁻⁴. This corresponds to a typical coronal temperature of 10⁶ K and a plasma beta of 1.11 × 10⁻⁴. At the boundaries, the magnetic field, velocity, density, and energy are fixed to their initial values and their derivatives normal to the boundaries are set to zero except where twisting velocity ‘driver’, described below, is prescribed.

The flux rope is formed by prescribing a slowly accelerating, rotating flow at the upper z-boundary as

$$\begin{eqnarray*} {\boldsymbol u} = u_0 u_r(r) u_t(t) (-y, x, 0)^T, \end{eqnarray*}$$

(13)

where u_r(r) describes the radial profile of the twisting motion in terms of the radius r² = x² + y²,

$$\begin{eqnarray*} u_r(r) = u_{r0}(1 + \tanh (1 - r_d r^2)), \end{eqnarray*}$$

(14)

where r_d controls the radial extent of the driver, u_r0 is a normalizing factor, and u_t(t) describes the imposed acceleration of the twisting motion,

$$\begin{eqnarray*} u_t(t) = \tanh ^2(t/t_r), \end{eqnarray*}$$

(15)

where the parameter t_r controls the time taken to reach the final driver velocity u₀. The functions u_r(r) and u_t(t) are plotted in Fig. 1. At the lower boundary, the flow is in the opposite direction. This form of driver allows the system to be accelerated slowly enough that the production of disruptive shocks and fast waves is minimal. It is unavoidable that some waves are produced during the boundary acceleration; however, these usefully provide a source of noise which eventually forms a perturbation.

The driver velocity is set to u₀ = 0.15, the normalizing factor is u_r0 = 2.08, and setting r_d = 5 corresponds to a driver constrained to r < 1 and with a peak velocity at r ≈ 0.38. The ramping time is set to t_r = 2, resulting in an acceleration from 0 to u₀ over approximately five Alfvén times. These driver parameters correspond to a peak rotational period of T_R = 15.92, the length of time taken for one full turn to be injected by a single driver. Both drivers result in twist being added at a rate of 2|$\pi$| every 7.96 Alfvén times. The twist profile across the entire flux tube develops in such a way that by t ≈ 20, it is qualitatively similar to those studied by Quinn et al. (2020c), Hood et al. (2009), and Bareford & Hood (2015); however, the length of the flux tubes differs significantly. This configuration produces a z-directed tube of increasingly twisted magnetic field that eventually becomes unstable to both the flute instability and the helical kink instability.

The problem formulated above is solved numerically using the staggered-grid, Lagrangian–Eulerian remap code for 3D MHD simulations Lare3D of Arber et al. (2001), where a new module for anisotropic viscosity has been included as detailed by Quinn et al. (2021). The resolution used in the current work is 512 grid points per dimension, comparable to the highest resolution kink instability studies of Hood et al. (2009) or medium resolution studies of Bareford & Hood (2015).

4 RESULTS

In an attempt to extrapolate to coronal resistivity values, we focus the attention on two selected pairs of simulations, one pair where the background resistivity is set to η = 10⁻³ and another where η = 10⁻⁴. As in the work of  Quinn et al. (2020c), only background resistivity is used. Each pair consists of one simulation using isotropic viscosity (9) and another one using the anisotropic model (11). The value of viscosity is set to ν = 10⁻⁴ in all cases.

The overall development of both the flute and the kink instabilities is broadly similar for the two values of resistivity and is described initially. Similar simulations were performed with a longer flux tube of length 20 instead of the tubes with length 4 shown here, and the results were found to be qualitatively similar. Focus is then placed on the detailed description of instabilities in the η = 10⁻⁴ cases, with the aim of comparing the effects of the two viscosity models. Then further features of the η = 10⁻³ cases are summarized.

4.1 Mechanism and general development of instability

Initially and in all cases computed, the twisting at the upper and lower boundaries gives rise to a pair of torsional Alfvén waves that proceed to travel along the tube from the upper and lower boundaries to their respective opposite boundaries. The interaction between these waves produces an oscillating pattern in the kinetic energy with a period of approximately four Alfvén times, equal to the time taken for an Alfvén wave to traverse the entire length of the domain as visible early in Fig. 2(a).

As the field continues to be twisted, currents form, due to the local shear in the magnetic field, and heat the plasma through Ohmic dissipation. Due to the radial form of the driver, the magnitude of the current density is greatest at the axis of the tube, then decreases radially outwards as seen in Fig. 3(a). The orientation of the twisting produces a current flowing in the negative z-direction for r⪅ 0.5. At r ≈ 0.5 (corresponding to the radius of peak driving velocity), the current switches orientation and is in the positive z-direction in a shell where 0.5 ⪅ r ⪅ 0.8. This form of a twisted field with an inner core of current in one direction surrounded by a shell of oppositely directed current is similar to the current configuration arising due to the field prescribed by Quinn et al. (2020c).

This current profile is reflected in the radial Ohmic heating profile (Fig. 3b) and, consequently, in the radial gas pressure profile (Fig. 3c). The highly pressurized core extends to r ≈ 0.2–0.4 (depending on the value of η) before increasing slightly around r ≈ 0.7. The secondary bump in gas pressure is due to the outer shell of current. The gas pressure gradient near the axis provides the outwardly directed gas pressure force which competes against the binding action of the magnetic tension and this provides the mechanism of flute instability excitation. The magnitude of the gas pressure gradient depends strongly on the value of resistivity η, with lower values producing shallower gradients which (as shall be seen) are more stable to the flute instability. Indeed, when η = 0, the radial gas pressure profile is nearly flat and the tube stable to the flute instability.

In all cases unstable to the flute instability, it occurs between t = 20 and t = 30. The continued driving at the boundaries eventually injects enough twist that the tube also becomes unstable to the kink instability. The kink initially develops linearly alongside or shortly after the flute instability and then erupts during its nonlinear phase, dominating the dynamics and disrupting the flute instability. The onset and the competition of the two instabilities is strongly affected by the value of η and the viscosity model used.

4.2 Instabilities at resistivity η = 10⁻⁴

We now describe the evolution and competition of flute and kink instabilities in case of resistivity η = 10⁻⁴. Fig. 4 shows the gas pressure profile of the anisotropic viscosity case (11) at time moments t = 26, 28, and 30 and at z = 0. At t = 26, the plasma begins to bulge out diagonally from the high-pressure core displaying an azimuthal wavenumber m = 4 as seen in Fig. 4(a) and indicating the presence of the flute instability. As the bulges move radially outwards into lower gas pressure regions they expand and accelerate, resulting in the entire gas pressure structure taking the shape of a four-leaf clover (Fig. 4b). By t = 30, the kink instability has disrupted the flute instability and is developing non-linearly as evident in Fig. 4(c). As is typical of non-linear kink development, the tube continues to release its stored potential energy in the form of kinetic energy and heat and the contained plasma becomes highly mixed. In the isotropic viscosity case which will not be illustrated by a separate figure, the flute instability is present but its growth is damped relative to the anisotropic case, and it is quickly outcompeted by the kink instability that dominates the dynamics.

Fig. 5 shows the effect the viscosity models have on the initial stages of the flute and kink instabilities in 3D. While the flute instability is observed in both cases, it is more pronounced in the anisotropic case, where it appears to disrupt the inner core of field lines and, as will be discussed further below, slows the growth of the kink instability. In the isotropic case, the growth of the flute instability is damped relative to the anisotropic case to the extent that the kink instability grows uninhibited and quickly disrupts the fluting.

Despite the flute instability appearing in the isotropic case (Fig. 5a), only in the anisotropic case can phases of linear growth of the flute and kink instabilities be seen in the kinetic energy profile as shown in Fig. 2(b). Here, the growth rates of the two instabilities are found to be γ = 0.69 for the flute and λ = 2.55 for the kink. The apparent phases of linear growth as measured from the kinetic energy time series, start at approximately t = 27 for the flute instability and t = 29.5 for the kink. In the isotropic case, the growth rate of the kink, λ = 2.97, is larger than in the anisotropic case, while the kinetic energy profile shows no evidence of flute instability growth.

The faster growth of the kink compared to that measured by Quinn et al. (2020c) is attributed to the relative aspect ratios of the flux tubes. The tube prescribed by Quinn et al. (2020c) has an aspect ratio of approximately 20 compared to the tube studied here which has an aspect ratio of approximately 4. While the total twist is similar in both tubes (after the drivers have injected twist up to t ≈ 20) the small aspect ratio results in more turns per unit length, leading to a faster growing instability.

Prior to the onset of either instability, the flux tube is found to be linearly unstable to perturbations of the form (1) at t = 20 via Suydam’s criterion (3) as shown in Fig. 6(a). The criterion represents a balance between destabilizing pressure gradients and stabilizing magnetic shear, and in this case, the shear is so small and the pressure gradient so large that the tube is unstable over a wide range of radii, for 0.02 ⪅ r ⪅ 0.29. Using equation (6), the linear fluting growth rate γ is plotted as a function of r at t = 20 in Fig. 6(b). At any fixed moment, the radial dependence of the flute instability growth rate, γ(r), is a concave function and peaks at a certain radius that we denote by r_s in Fig. 6(b) and Table 1.

Figure 6.

Stability and linear growth rate of the flute instability. In panel 6(a), Suydam’s stability criterion (3) and its contributing terms are plotted and in panel 6(b) the predicted linear growth rates for the ideal (6) flute instability are plotted. Both plots are produced at t = 20 for η = 10⁻⁴ and using the anisotropic model. The location of the peak linear ideal growth rate r_s is also shown.

Table 1.

Quantitative differences in the observed perturbations between results with different viscosity models and resistivity values η. Measurement times are listed in the main text.

	η = 10−4		η = 10−3
	Anisotropic	Isotropic	Anisotropic	Isotropic
Theoretical average ideal linear growth rate of flute γ	0.88	0.88	1.73	1.73
Observed growth rate of flute γ	0.69	–	1.06	1.06
Observed growth rate of kink λ	2.55	2.97	0.15	0.15
Theoretical radius rs of peak ideal flute growth rate	0.125	0.125	0.163	0.163
Predicted axial wavenumber kflute	23.74	23.52	17.15	17.60
Observed axial wavenumber kflute	22.93	22.30	16.05	16.05
Observed axial wavenumber kkink	4.57	4.41	3.44	3.49
Cumulative viscous heat at t = 50	1.64 × 10−2	1.28	2.89 × 10−3	0.370
Cumulative Ohmic heat at t = 50	3.06	2.75	6.54	6.27

The location r_s of the peak fluting growth rate aligns well with the location of the resonant surface, where the observed perturbation appears to grow (Fig. 4a) and an estimate of the linear growth rate can be found by averaging γ over r, giving a theoretical growth rate of 0.88.

The flow and pressure profiles in the axial direction z at t = 26 are shown in Fig. 7(a), and at this moment, they assume the form of a superposition of the unstable modes with the largest amplitude. In particular, the fluting perturbation is most easily observed in the pressure profile and the kink instability is best revealed by either of the x- or the y-component of velocity which can serve as proxies for the radial velocity through the axis. Comparing the magnitudes of the profiles at this time suggests the flute instability is close to transitioning to its non-linear phase while the kink instability is still very much in its linear phase.

Figure 7.

Selected flow and pressure profiles and the resonance function defined by equation (2). (a) Gas pressure and velocity profiles in z at the fixed point (r, θ) = (0.101, 0), assuming the form of the most unstable modes. (b) The resonance surface mB_θ(r)/r + kB_z(r) as a function of r using the observed fluting perturbation wavenumbers. All plots are snapshots at t = 26 where η = 10⁻⁴ and the viscosity model is anisotropic.

The value of k for each instability mode is calculated as |$k = 2\pi /\tilde{\lambda }$|⁠, where |$\tilde{\lambda }$| is the wavelength of the perturbation, measured as the difference between the two peaks closest to z = 0 in Fig. 7(a) (thus minimizing the influence of line-tying on the measurement). This gives a value of k_flute = 22.93 and k_kink = 4.57 for the anisotropic model and k_flute = 22.30 and k_kink = 4.41. Using these values, it is observed that the fluting mode resonates with the field, that is mB_θ(r)/r + kB_z(r) ≈ 0, over a range of 0.15 ≲ r ≲ 0.225 (Fig. 7b). This is in close agreement with the predicted radius of peak linear flute growth r_s seen in Fig. 6(b).

Comparing the effect of the viscous models on the modes, in the isotropic case the fluting mode is damped, while in the anisotropic case the kink mode is diminished, explaining why the flute instability appears more readily in the anisotropic case (Fig. 2a).

4.3 Instabilities at resistivity η = 10⁻³

Fig. 8 shows a prolonged development of the flute instability and a slow nonlinear development of the kink instability at the higher resistivity value η = 10⁻³ in the case of anisotropic viscosity. Due to the enhanced Ohmic heating at η = 10⁻³, the gas pressure gradient is substantially stronger than at η = 10⁻⁴ and the flute instability is excited much earlier. Compared to the η = 10⁻⁴ cases, the instability occurs further from the axis, at r ≈ 0.16, and the larger gas pressure gradient drives the bulges in profile further from the axis during the non-linear phase (Fig. 8a). These bulges continue to extend outwards and mix the plasma as they develop. The kink instability can be observed displacing the axis of the tube diagonally upwards and to the right in Fig. 8(c). At this time in the η = 10⁻⁴ cases, the non-linear development of the kink was at a later stage of its development (Fig. 4c). The development of the kink then proceeds slowly, as it moves the axis of the tube further through the mixed region to eventually begin the reconnection process with the outer region of field that is typical of the instability in this kind of flux tube as was observed by Quinn et al. (2020c).

Figure 8.

Gas pressure profiles at z = 0 during the development of the flute and kink instabilities in the higher resistivity anisotropic case. The viscosity model is anisotropic and η = 10⁻³. In contrast to the case of η = 10⁻⁴, the non-linear development of the flute instability has time to mix the interior of the flux tube before the onset of the kink instability, the growth of which is likely affected by mixing of the plasma.

It is evident from the kinetic energy profile that the flute instability develops much earlier than in the η = 10⁻⁴ cases and grows at an increased rate of γ = 1.06 (Fig. 9b). The kink instability grows at a rate of λ ≈ 0.15, much slower than that observed in the η = 10⁻⁴ cases, and much lower than the flute instability. The time period between t ≈ 28 and t ≈ 32 is identified as the linear stage of the kink instability by inspecting the development shown in Fig. 8(a). This is broadly consistent with the tube surpassing a critical twist of 2.59π (Hood & Priest 1981; Török & Kliem 2003) between t = 30 and t = 35. One key observation is that, despite the early and disruptive growth of the flute instability, the kink instability still generates the bulk of the kinetic energy (Fig. 9a).

Figure 9.

Kinetic energy and heating as a function of time in the cases where η = 10⁻³. The results from both viscosity models are shown. The flute instability appears earlier than where η = 10⁻⁴ and the growth rate of the kink instability is decreased. The cumulative heating 〈Q_*〉, where * is either ν for viscous heating or η for Ohmic heating, is the respective heating term integrated both over space and from the initial moment up to the moment t in time. While the heat generated via viscous heating is orders of magnitude lower when using anisotropic viscosity, Ohmic heating is enhanced by the use of the anisotropic model.

Due to the influence of the drivers on the kinetic energy, the fluting growth rate is difficult to estimate from the kinetic energy profile as accurately as in the η = 10⁻⁴ cases. Since the kink instability occurs after the development of the fluting, its growth rate is similarly difficult to gauge. Nevertheless, it is clear that the flute instability grows at a rate of the same order as that in the η = 10⁻⁴ cases. It is also apparent that the kink instability grows much slower in the η = 10⁻³ cases.

Table 1 summarizes the quantitative differences between the results for the two models of viscosity and the two values of the resistivity. The theoretical average growth rate is computed as the mean across the radius of the ideal estimate (6) and is in good agreement with the observed rate in each case, particularly in the less resistive case η = 10⁻⁴ which better represents ideal plasma. The discrepancy between predicted and observed growth rates is due, in part, to the growth rate estimate (6) being derived under the assumption of asymptotically large values of m ≫ 1, while the observed mode has a finite value of m = 4. Despite this, the predicted growth rate is of a similar magnitude to the observed rate. The location r_s of the peak growth rate provides a prediction of where the instability will initially grow. This radius is used in conjunction with the resonance equation (2), with m = 4, to predict the axial wavenumber of the mode with the greatest linear growth, i.e. the mode most likely to be observed. Again, these are in good agreement with the observed fluting wavenumbers, which are measured at times just prior to the accelerated development of the flute instability, that is at t = 22 when η = 10⁻³ and t = 26 when η = 10⁻⁴. The kink wavenumber is measured at t = 26 in both cases. Overall, the agreement between predicted and observed growth rates and mode wavenumbers allows us to conclude that the observed instability is the ideal flute instability and that expression (6) can be effectively applied to estimate the growth rate of the flute instability within coronal loops.

Also listed are estimates for the cumulative heat generated via viscous and Ohmic heating during the simulations. As is also found in previous studies of viscous heating in kink instabilities (Quinn et al. 2020c), anisotropic viscous heating is approximately two orders of magnitude lower than isotropic and the use of anisotropic viscosity enhances Ohmic heating.

5 DISCUSSION

Perturbing a magnetic flux tube by dynamic twisting of the flow at the cylinder bases leads to excitation of the flute instability in addition to the well-studied kink instability. Our aim in performing the reported numerical experiments was to explore the flute instability in a dynamically twisted coronal flux tube, specifically focussing on the effect of anisotropic viscosity on the development of the instability. In addition, we wish to understand the effect the instability has on the proceeding kink instability, the effect on the overall heating generated through viscous and Ohmic dissipation, and the effect that varying resistivity has on the development of the flute. Our findings are discussed below.

We have found evidence of the flute instability using both models of viscosity; however, isotropic viscosity damps the initial growth such that it does not develop beyond its linear phase before the faster growing kink instability disrupts the flux tube and halts the development of the flute. Given that most numerical studies of the kink instability employ isotropic or shock viscosity, this is likely why the flute instability has not been previously reported.

Counter-intuitively, the growth rate of the kink instability is lower in the weakly dissipative anisotropic cases, compared to the strongly dissipative isotropic cases as one would expect the kink instability to grow more quickly (or at least be unaffected) when using anisotropic viscosity. Indeed, the simulations reported by Quinn et al. (2020c) display this behaviour, where the kink instability grows faster in the anisotropic cases. We speculate that it is the presence of the flute modes that negatively affects the growth of the kink instability. It seems unlikely that in the linear regime the flute and kink modes are able to directly couple, given that the kink instability generally presents at the axis of a flux tube and the flute at some resonant surface away from the axis. We believe that, instead of a direct coupling, the linear kink instability is disturbed by the more complex magnetic field configuration that arises due to the mixing caused by the nonlinear development of the flute modes. The complexity in the field can be seen by comparing Figs 5(a) and (b).

Beyond the effect on the growth of the two observed instabilities, the two viscosity models greatly affect both the viscous heating and Ohmic heating rates illustrated in Figs 2(c) and 9(c). Anisotropic viscosity is naturally less dissipative than isotropic viscosity and generates approximately two orders of magnitude less total viscous heat than isotropic viscosity. This is somewhat offset by anisotropic viscosity permitting greater kinetic energy release and enhanced mixing, which in turn enhances Ohmic dissipation of heat through the generation of strong localized current sheets. The overall effect is that more heat is generated when the viscosity is anisotropic, similar to what has been observed in previous work (Quinn et al. 2020c). In the context of coronal heating, this is encouraging: the use of a less dissipative viscosity model actually results in greater overall heating. How this finding generalizes to more realistic coronal resistivities, and whether it holds true for other coronal instabilities, should be the subject of further investigations.

It is difficult to distinguish the effect of the flute instability on the viscous or Ohmic heating from that of the viscosity itself, particularly since the flute instability is quickly disrupted by the kink instability. However, it can be concluded from the plots of cumulative heat (Figs 2c and 9c) that the bulk of the viscous and Ohmic heat is generated in the non-linear phase of the kink instability. There is additionally some non-negligible Ohmic heating generated prior to the onset of the kink; however, this is attributed to the large-scale currents associated with the twist in the field, rather than any currents created by the flute instability. This leads us to conclude that the flute instability itself has little direct impact on coronal heating, but instead can affect the heating rate by slowing the development of the kink instability.

It is likely that the m = 4 azimuthal mode is excited due to influences from the boundaries in the Cartesian box, for example through the interaction of reflected fast waves generated in part by the driver. Performing a similar experiment in a cylindrical numerical domain, or prescribing a variety of perturbations in the Cartesian domain may reveal other, faster growing modes. The modes may also be influenced by nonlinear coupling between the m > 1 and m = 1 modes, as is found in the study of kink and flute oscillations (Ruderman 2017; Terradas et al. 2018).

As the current distribution, which develops as the flux tube is twisted, is similar to that found in the initial flux tube configuration of Quinn et al. (2020c), the question arises why the fluting instability is not observed in the latter. Although the current distribution (and thus heating and pressure distributions) in the tubes of Quinn et al. (2020c) may support the flute instability, the tube is initially perturbed with a motion close to an unstable eigenmode of the kink instability, resulting in the instability growing immediately from t = 0. In contrast, in the tubes studied here, such a perturbation must grow from numerical noise, allowing a secondary, fluting modes perturbation to also develop and become significant enough to observe. As an explorative alternative to prescriptive perturbations, we recommend the use of numerical noise in the study of coronal instabilities.

Since the main driver of the flute instability is the gas pressure gradient generated through Ohmic heating, it is prudent to ask if the same gas pressure gradient could be generated using physical coronal values of the resistivity, which are estimated to be approximately η = 10⁻⁸ (Craig & Litvinenko 2009), and are thus much smaller than those studied here. Additionally, the simulations presented here do not incorporate radiation or thermal conduction, two processes which would remove energy (and hence reduce gas pressure) from high-pressure regions in a coronal loop and thus could prevent meeting the required conditions for the growth of a flute instability. Indeed, at η = 10⁻⁴ the flute instability was more quickly outcompeted by the kink instability and appeared to have little impact on the resultant dynamics, which mirror those of other kink instability studies (Hood et al. 2009). This suggests that even lower values of resistivity would result in flux tubes without any significant flute instability, at least for this form of driver and mechanism of gas pressure generation. Regardless, coronal loops with strong radial gas pressure gradients have been observed (Foukal 1975), and such loops may be unstable to the flute instability. Modelling of a prescribed flute-unstable flux tube, as opposed to the dynamically stressed loop investigated here, would provide a useful comparison to observations; however, it may be difficult to prescribe a tube which is not also susceptible to kinking. Linear stability analyses of this kind of flux tube (a dynamically created zero total axial current tube) focus on the kink instability (Browning & Van der Linden 2003) so do not provide much insight into the potential for fluting without a kink.

While our results show that a flux tube can be unstable to the flute instability and yet the faster growing kink instability can quickly dominate when the gas pressure gradient is small enough, the opposite case is also observed. A faster growing flute instability appears to slow the growth of the kink instability although, importantly, it does not fully disrupt the development of the kink. Understanding the balance between the non-linear growth rates of the two instabilities is important for prediction of whether the flute instability may be found at all in the real solar corona, or whether its realistic growth rate is too slow compared to that of the kink instability.

6 CONCLUSION

This paper details the non-linear development of two ideal instabilities, the kink and the flute instabilities, both of which develop naturally in the course of twisting an initially straight magnetic flux tube. This provides a different approach to that employed in the simulations performed in the earlier study by  Quinn et al. (2020c) in that the instabilities are not excited by any prescribed perturbations but, instead, the field is dynamically driven into an unstable state and the perturbations provided by noise in the system. Not only is the kink instability excited due to the twist in the field, but nearly simultaneously a pressure-driven flute instability can also be excited in unstable pressure gradients generated by Ohmic heating. Simulations were performed with two values of resistivity, η = 10⁻³ and 10⁻⁴, and for two forms of viscosity, isotropic, and anisotropic. The results prove an initial and important first step towards understanding non-linear flute instabilities in the solar corona.

It has been shown that the flute instability can be quickly dominated by the kink instability if the kink grows substantially faster than the flute. However, if the flute has time to develop non-linearly, it mixes the plasma within the flux tube, generating small-scale current sheets and releasing some magnetic energy. The overall effect of this mixing is to slow the growth of the kink instability. The slowed growth of the kink does not appear to significantly impact the kinetic energy released during its evolution, only the time over which it is released.

These numerical experiments have provided evidence that the flute instability can occur in twisted magnetic flux ropes and grow at similar rates to the kink instability. Further estimation of the relative growth rates in more realistic coronal loop set-ups is required to fully understand if the flute instability plays a pertinent role in coronal loop physics.

ACKNOWLEDGEMENTS

We would like to thank David MacTaggart for his input on the thesis chapter on which this paper is based. JQ was funded via an EPSRC studentship: EPSRC DTG EP/N509668/1.

Results were obtained using the ARCHIE-WeSt High Performance Computer (www.archie-west.ac.uk) based at the University of Strathclyde.

The authors acknowledge the use of the UCL Myriad High Performance Computing Facility (Myriad@UCL), and associated support services, in the completion of this work

This work used the DiRAC Extreme Scaling service at the University of Edinburgh, operated by the Edinburgh Parallel Computing Centre on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk), specifically Thematic Project allocation ACTP245. This equipment was funded by BEIS capital funding via STFC capital grant ST/R00238X/1 and STFC DiRAC Operations grant ST/R001006/1. DiRAC is part of the National e-Infrastructure.

APPENDIX: DATA AVAILABILITY

A custom version of Lare3d (Arber et al. 2001) has been developed, where a new module for anisotropic viscosity has been included. The version including the new module can be found at https://github.com/jamiejquinn/Lare3d, and has been archived by citeQuinn et al. (2020a). The version of Lare3d used in the production of the results presented here, including initial conditions, boundary conditions, control parameters and the anisotropic viscosity module, can be found in the repository of  citeQuinn et al. (2020b). Associated running scripts for generating, building and running simulations on a cluster is also provided Quinn (2022a). The data analysis and instructions for reproducing all results found in this report may be also found at https://github.com/JamieJQuinn/coronal-fluting-instability-analysis and has been archived (Quinn 2022b).

All simulations were performed on a single, multicore machine with 40 cores provided by Intel Xeon Gold 6138 Skylake processor running at 2 GHz and 192 GB of RAM, although this amount of RAM is much higher than was required; a conservative estimate of the memory used in the largest simulations is around 64 GB. Most simulations completed in under 2 d.