The Effect of Anisotropic Viscosity on the Nonlinear Kink Instability

The kink instability of magnetohydrodynamics is believed to be fundamental to many aspects of the dynamic activity of the solar atmosphere, such as the initiation of flares and the heating of the solar corona. In this work, we investigate the importance of viscosity on the kink instability. In particular, we focus on two forms of viscosity; isotropic viscosity (independent of the magnetic field) and anisotropic viscosity (with a preferred direction following the magnetic field). Through the detailed analysis of magnetohydrodynamic simulations of the kink instability with both types of viscosity, we show that the form of viscosity has a significant effect on the nonlinear dynamics of the instability. The different viscosities allow for different flow and current structures to develop, thus affecting the behaviour of magnetic relaxation, the formation of secondary instabilities and the Ohmic and viscous heating produced. Our results have important consequences for the interpretation of solar observations of the kink instability.


Introduction
The solar corona is heated to millions of degrees, whereas the surface of the Sun beneath it is only heated to the order of thousands of degrees [1]. It is generally accepted that this heating is derived mainly from the extraction of energy from the stressed magnetic field in the solar atmosphere, although exactly how this process occurs is still an area of active research [2]. There has been much study of the direct transfer of magnetic energy to heat by Ohmic heating and the transfer through kinetic energy losses by viscous heating [2,3]. The relative importance and efficacy of these processes is found to depend on the models of plasma viscosity employed in the analyses. The aim of this work is to understand the effect of isotropic and anisotropic viscosity on the kink instability, a key phenomenon in coronal plasma dynamics.
In many situations in the solar corona, the strength of viscosity can be greater than that of resistivity by many orders of magnitude, even when anomalous resistivity is considered. Due to this, viscous heating can outperform Ohmic heating in certain coronal situations, particularly those involving reconnection [3,4,5,6]. These findings are dependent on how viscosity in the solar atmosphere is modelled. For a plasma in the presence of a strong magnetic field, viscosity is anisotropic and momentum transfer occurs preferentially in the direction of the magnetic field [7]. A general description of anisotropic viscosity in the solar atmosphere is given in [8,6]. More recent studies have demonstrated the importance of anisotropic viscosity for heating in investigations of three-dimensional (3D) null points [4], current sheet merging [5] and flux pile-up [9]. There is further evidence of the importance of anisotropic viscosity in other astrophysical applications, including the intracluster medium [10,11] and the solar wind [12]. In other solar applications, viscosity has a role to play in the damping of coronal instabilities [13] and waves [14,15,16], though not all these cases have been fully explored using an anisotropic model of viscosity.
Implementing the full Braginskii viscosity tensor [7] into existing magnetohydrodynamic (MHD) codes is not trivial. In the solar atmosphere, there are regions where the magnetic field is weak or vanishes (called null points) and at these locations the viscosity must transition smoothly from strongly anisotropic (in the direction of the magnetic field) to isotropic. Although the Braginskii tensor theoretically achieves this smooth transition, in practice the implementation of the full tensor in finite-difference MHD codes leads to numerical errors in the viscosity at locations where the magnetic field strength is very weak, due to a lack of sufficient resolution. Further exploration of this numerical issue can be found in [17]. One approach to solving this problem is to design a numerical scheme specifically to treat the full Braginskii tensor. Another approach is to devise a model that captures the main physics of the Braginskii tensor suitable for the solar corona and can be implemented in existing MHD codes [8]. The second option was taken up by MacTaggart, Vergori and Quinn [17], who developed a phenomenological model of anisotropic viscosity in the solar corona that captures the main physics for viscosity in the corona, namely parallel viscosity in regions of strong field strength and isotropic viscosity in regions of very weak or zero field strength. In this paper, we will refer to this model of viscosity as the switching model, for short. The model interpolates between isotropic and parallel viscosity based on how the local magnetic field strength affects the distribution of momentum transport. The interpolation itself can be adjusted, effectively changing the size of isotropic (i.e. weak field) regions, to compensate for resolution issues in simulations. In numerical tests of stressed null points, comparing the isotropic and anisotropic models showed that the use of isotropic viscosity overestimates the viscous heating by an order of magnitude [17].
As a first step of investigating the effects of anisotropic viscosity on coronal dynamics, we focus on the kink instability [18,19], believed to be a trigger for flares [20] and an important mechanism in the theory of coronal heating through nanoflares [21]. The instability has also been studied using shock viscosity [19,22] but a detailed investigation of the effects of Newtonian and Braginskii viscosity has not, to the best of our knowledge, been performed. In particular, the main aim of this paper is to provide insight into the effect of the choice of viscosity model on the nonlinear dynamics and relaxation of a twisted coronal loop, where the kink instability converts magnetic energy to heat through Ohmic heating generated via current structures and through viscous heating generated via flow structures. We aim to give an estimate of how well viscous heating (using both isotropic and anisotropic models) performs when compared with Ohmic heating. This study extends previous work [19] which also considers the kink instability in a zero-current loop (details given below). However, in contrast to [19], we consider only background resistivity and viscosity as the two mechanisms of heat generation (that is, we do not consider shock viscosity or anomalous resistivity). In performing this investigation, we also provide further validation of the switching model in a simpler topology to that used by MacTaggart et al. [17], in that there are no null points present in the field at any time.
The layout of the paper is as follows. The coronal loop model, MHD equations and viscosity models are described in Section 2. Simulation details and quantities used in the analysis of the simulations are described in Section 3. Detailed numerical results of a typical case of a kink instability are presented in Section 4 with a particular focus on how the different viscosity models affect its nonlinear evolution. The results of the typical case are confirmed and generalised by a parameter study in Section 5 where the dependences of the Ohmic and the viscous heating on the resistivity and the dynamic viscosity are presented. Our conclusions are summarised in Section 6.

Models of kink instability and viscosity
As a model of an idealised coronal loop, we consider a twisted flux tube in a Cartesian box of dimension [−2 Mm, 2 Mm] × [−2 Mm, 2 Mm] × [−10 Mm, 10 Mm] in the x, y and z-directions, respectively. The state of the plasma is typical of the corona, with density ρ initially 1.67 × 10 −12 kgm −3 , and with plasma pressure p such that the temperature of the plasma T is initially 2 × 10 4 K everywhere in the computational domain. The magnetic field B is constructed so that it is initially force-free and with zero axial current, line-tied at the boundaries, and twisted such that it is linearly unstable to the ideal kink instability. This configuration allows us to compare directly to previous studies that use identical magnetic field configurations [19,22,23].  The field outside the flux tube is straight and has a strength of 5 × 10 −3 T. Given this temperature and magnetic field strength, the plasma beta is initially β ≈ 10 −5 , a value realistic for the corona. The evolution of this flux tube is governed by the nonlinear MHD equations.

Model equations
Non-dimensionalised, the MHD equations take the form where u is the plasma velocity,  = ∇ × B is the current density, σ is the viscous stress tensor, η = 1/S is the normalised resistivity (equivalent to the inverse of the Lundquist number S), and use has been made of the material derivative, D/Dt = ∂/∂t + (u · ∇). The internal energy density is given by the equation of state for an ideal gas ε = p ρ(γ − 1) , where the specific heat ratio is given by γ = 5/3. The terms Q ν = σ : ∇u and Q η = η|| 2 model Ohmic heating and viscous heating, respectively. Using the nondimensionalisation scheme found in [24], reference values for the magnetic field B 0 , length L 0 and density ρ 0 are chosen to align with values typical for a coronal loop. The problem domain is nondimensionalised to [−2, 2] × [−2, 2] × [− 10,10] in the x, y and z-directions, respectively. Velocity and time are non-dimensionalised using the Alfvén speed u A = B 0 / √ ρ 0 µ 0 and Alfvén crossing time respectively. Temperature is non-dimensionalised via T 0 = u 2 Am /k B , where k B is the Boltzmann constant andm is the average mass of ions, here taken to bem = 1.2m p (a mass typical for the solar corona) where m p is the proton mass. The reference values of other variables are derived from these quantities and are listed in Table 1. Dimensional quantities can be recovered by multiplying the non-dimensional variables by their respective reference value (e.g. B dim = B 0 B). All further reference to variables will be to their non-dimensionalised values, unless stated otherwise.
We consider a force-free magnetic field for which the Lorentz force is zero so that (∇ × B) × B = 0. To satisfy this condition we choose, in cylindrical coordinates (r, θ, z), a magnetic field of the form B = α(r)∇×B, where α(r) is a function of a particular form that ensures the total axial current is zero. Aligning with previous work by Hood et al. [19], the smooth α(r) profile given as Case 3 in [19] is used. Using this profile, the equilibrium magnetic field B is written as  Figure 1: The initial field configuration. In (a) field lines are plotted corresponding to inner (red), outer (blue) and straight (yellow) regions of twist, with slices of α(r) shown at the footpoints. In the slices, red corresponds to α(r) > 0, blue to α(r) < 0 and white to α(r) = 0. In (b) the profiles of α(r) and the field components Bz and B θ across the flux tube are plotted.
for r ≤ 1 and B θ = 0 for r > 1, where λ is a parameter measuring the twist in the tube. The radial field throughout the domain is set to B r = 0. As is done in [19], we set λ = 1.8 to ensure the tube is unstable to the ideal kink instability. The equilibrium velocity for this magnetic field configuration is u = 0.
The form of α(r) in equations (6) and (7) splits the profile of the flux tube into three twist regions, the inner region of positive twist (r ≤ 0.5), the outer region of negative twist (0.5 < r < 1) and the straight-field region of zero twist (r ≥ 1) as shown in Figures 1(a) and (b). These figures also illustrate the equilibrium field. Since the inner region is more tightly twisted, this field configuration results in only the inner region becoming unstable to the kink instability, rather than the global instability seen in non-zerocurrent loops [18]. The regions of twist are used later to define a measure of reconnection.
Although we prescribe an initial temperature of T = 2 × 10 4 K, the equations simulated by the code are written using internal energy, thus we convert the temperature to internal energy using the non-dimensional relation ε = T /(1 − γ). Hence, the initial non-dimensionalised density and internal energy are uniformly given by ρ = 1, ε = 8.66 × 10 −4 , and have been non-dimensionalised using the reference values found in Table 1. The initial magnetic field and velocity are set to their equilibrium states, discussed above, with the addition of a small perturbation. In order to make a meaningful comparison of our results with those of [19], we use identical initial magnetic field and velocity perturbations, calculated via a linear stability analysis (in ideal MHD) applied to a similar flux tube that uses a constant, piecewise profile for α(r) [25,26,21].
At the boundaries, we satisfy the line-tied condition on the magnetic field by ensuring the field is constant and equal to its initial values given by equations (6) and (7). Similarly, on the boundaries we ensure that the density, internal energy and velocity u are constant and equal to their initial values. To close the system, the fluxes of all variables through each of the boundaries are set to zero. That is, on the x-boundary, and similarly, the y and z derivatives are zero on the y = ±2 and z = ±10 boundaries, respectively.

Anisotropic viscosity
Isotropic (Newtonian) viscosity is the most commonly used viscosity model for applications to the solar corona. Its viscous stress tensor is directly proportional to the rate-of-strain of a flow through the dynamic viscosity of the fluid, ν, where W = ∇u + (∇u) T − 2 3 (∇ · u)I is the traceless rate-of-strain tensor and I is the identity tensor. In the solar corona, apart from regions of very weak field strength, viscosity is dominated by the strong field limit of the full Braginskii tensor [6], where b = B/|B| is a unit vector in the direction of the magnetic field. The strong field tensor can be considered made up of two parts: the scalar strength of viscosity, given by the multiplication of the viscosity coefficient ν and the interaction of the strain rate tensor with the field (W b) · b; and the direction of action, involving only the direction of b. To illustrate how velocity gradients affect this tensor, consider a magnetic field where the frame of reference is such that the field is oriented in the z-direction, that is b = (0, 0, 1) T .
In this frame, and the full tensor in equation (11) can be written, in matrix form, as Consider a velocity only in the z-direction, u z , exhibiting a gradient in the same direction, ∂u z /∂z. In this case, equations (10) and (11) are equivalent, hence the anisotropic viscosity acts identically to isotropic viscosity in the direction of the field. In contrast, if that same velocity u z exhibits only a gradient perpendicular to the field, say ∂u z /∂x, there would be no strong field viscosity, since the shear gradient does not enter into expression (13), and the flow would be undamped. Similarly, if the flow consisted of only perpendicular velocities with gradients directed along the field, the viscosity would vanish and these flows would remain undamped. Indeed, a flow exhibiting only shear velocity gradients would remain undamped. For a further exploration of the behaviour of this tensor, we refer to the more comprehensive discussion of the terms of the full Braginksii tensor in [7] and to the physical interpretation found in [8]. As mentioned previously, the switching model of MacTaggart et al. [17] interpolates between equations (10) and (11), The interpolation function s(|B|) is derived in [17] by considering how the distribution of momentum transport depends on the magnetic field strength. The interpolating function is given by where a(|B|) is a constitutive function, chosen here to take its simplest, non-trivial form a(|B|) = a 0 |B| 2 . In practial terms, the parameter a 0 varies the size of the region where isotropic viscosity dominates in the numerical simulations. To align with earlier work in [17], we choose a 0 = 150. The application of the switching model in simulations of a stressed null point in [17] provides a fuller exploration of the isotropic feature of the model. Despite the switching itself being relatively unimportant here due to viscosity remaining fully anisotropic nearly everywhere, we still choose to use the switching model for two reasons: consistency with previous work [17], and as a tool to investigate the effect of prescribing anisotropy, which we will discuss in more detail later.
The switching viscosity is not only a physically realistic model for anisotropic viscosity in the corona, but can also be used to investigate, more deeply, the differences between isotropic and anisotropic viscosity. This is done by artificially setting the result of s in the code to fix the amount of anisotropy present. We will consider the balance between isotropic and anisotropic viscosity in more detail later.

Numerical setup and tools of analysis
Here we describe the computational code used to solve numerically the governing equations and we define the diagnostic tools used to analyse the simulation output results.

Numerical setup
We solve the governing MHD equations (1)-(4) numerically using the Lare3d code [24], which is widely used in the solar physics research community. Lare3d is based on a Lagrangian-Remap scheme with artificial viscosity and flux-limiters as shock-capturing devices. Since we are investigating the role of viscosity itself, the artificial viscosity (or shock viscosity) has been disabled. In order to compare features of our simulations with the results of Hood et al. [19] we have performed tests using shock viscosity, instead of either the switching or isotropic models. Using the default shock viscosity parameters present in the code, the behaviour closely mirrors that of isotropic viscosity with ν ≈ 5 × 10 −4 . When both switching and shock viscosity are enabled, the shock viscosity dominates and, again, the behaviour mirrors that of isotropic viscosity.
The simulations were run at a resolution of 350 × 350 × 700, with the exception of the parameter studies, which were run at a slightly higher resolution of 400 × 400 × 800. Since the switching viscosity only acts parallel to the magnetic field, in perpendicular directions numerical diffusion dominates. By running several simulations at resolutions of 250 × 250 × 500 up to 500 × 500 × 1000, it was found that the effect of resolution was negligibly small until around t = 150, well after the nonlinear phase of the instability. After this time there are some quantitative differences in outputs for different resolutions. However, the qualitative behaviour, which we describe later, does not strongly depend on the resolution.
We use an estimate ofν =η =ũ x L x /N 2 x for the numerical diffusion present in the simulations due to the finite difference scheme employed in Lare3d. Taking a typical velocity ofũ x = 1, i.e. the Alfvén velocity; N x = 350 as the number of grid-points in the x-direction; and L x = 4 as the length in the x direction, we estimate the numerical diffusion coefficient to beν =η ≈ 10 −5 . This provides a theoretical lower bound on simulating a physical viscosity or resistivity. In practice, however, we find setting the physical resistivity lower than η ≈ 5 × 10 −5 results in behaviour that does not converge with increasing resolution. This gives a practical lower bound for diffusion coefficients ofν =η ≈ 5 × 10 −5 . Thus, all results presented use physical diffusion coefficients (either viscosity or resistivity) greater than this lower bound.

Connectivity
As part of the analysis in Section 4, we present a practical measure of magnetic reconnection -the mean change in field line connectivity ∆Φ c . We determine the connectivity Φ c for field lines by analysing the start and end points of a sample of magnetic field lines. Any field line that begins at one location in one of the footpoints will map to a corresponding point in the opposite footpoint. Following field lines from their starting point at z = −10 to their end point at z = 10, we label lines depending on the twist regions in which they start and end. Initially, the field lines within each distinct region map one-to-one to the same region. As the field reconnects radially during the instability, field lines begin to start and end in different twist zones. By tracking the number of field lines that have changed twist zones within a set time period we get a practical estimate for the rate of reconnection, as well as a visual guide of where that reconnection is occurring. It should be noted that this measure of reconnection does not take into account azimuthal reconnection (that is reconnection within the same twist zone). As such, it is only a partial measure of reconnection.
In practice, magnetic field output is saved from the code at intervals of ∆t = 5. We use the visualisation tool Mayavi [27] to compute the magnetic field lines over a grid of starting points (x i , y j ) at a given time n∆t, where n indexes the output files. This process gives a connectivity map Φ at time n∆t is found by subtracting one connectivity map from the previous and then taking the mean across all points (x i , y j ),

Parallel electric field
Another useful measure of magnetic reconnection is the maximum value of the integral of the electric field parallel to the magnetic field E = η( · B)/|B| along a magnetic field line [28,29,30], where C is a magnetic field line with start and end points within the footpoints at z ± 10.
Similarly to the calculation of connectivity, we employ Mayavi to compute magnetic field lines using a grid of field line starting points (x i , y j ) at a given time. We track the local value of the modulus of the parallel electric field, |E | = |η · B| along the magnetic field lines. The parallel electric field is then summed along each of the field lines to give a distribution Φ(x i , y j ) across the profile of the field. The maximum of this distribution gives a measure of the reconnection rate.
It can be argued that the reconnection rate calculated by taking the global maximum is only the rate for one region of magnetic diffusion, and the nonlinear phase of the kink instability creates multiple diffusion regions in its development. One way to calculate the reconnection rate for each region is via the algorithm described in [31], which dissects the distribution Φ(x i , y j ) into separate regions before finding the maxima corresponding to the reconnection rate per diffusion region. In practice, we find the current structures created by the kink instability to be simple enough that this extended analysis is unnecessary.

Other observables
To further characterise our results we use of the volume-integrated parallel and perpendicular kinetic energies, the magnetic energy, and the total Ohmic heating generated by time T , The time and volume-integrated viscous heating rate can be written in the form for the isotropic viscous stress tensor (10) and in the form for the switching viscous stress tensor (14), respectively. The linear growth of the instability ends around t = 35 and the inner field compresses into the outer field, creating a current sheet. Between t = 45 and 50 this current sheet enables reconnection between the two regions. The transition for the switching case is qualitatively similar. In all three plots, the viscosity and resistivity are ν = 10 −4 and η = 5 × 10 −4.5 , respectively.

Nonlinear evolution of a typical case
In this section, we present results from simulations with a specific choice of viscosity and resistivity. This provides an opportunity to analyse, in detail, the onset and evolution of the kink instability in a typical case. Parameter studies illustrating that the observed dynamics is typical are presented further below in Section 5. We compare the results from simulations of two cases. Isotropic viscosity is used in the first case and switching viscosity is used in the second case. The choice of viscosity tensor is the only difference between the two cases. In this section, we use the diffusion parameters ν = 10 −4 , η = 5 × 10 −4.5 , both small but suitably above the threshold of numerical diffusion discussed in Section 3. The chosen value of ν is well within the range of typical values found in the real corona, that is, in our nondimensionalisation, between 10 −8 and 10 −3 [16]. All other parameters are identical in both cases and are kept fixed to the values specified in Section 3. Due to the strength of the field and lack of null points, it is measured that s = 1 throughout the entire domain, thus the switching model reverts to the strong field approximation of the Braginskii tensor (11).

Linear phase
The linear development of the kink instability lasts until t ≈ 35 as illustrated in Figure 2(a) and has a measured linear growth rate of σ = 0.13. Since the initial velocity perturbation is calculated from an ideal and inviscid MHD model with a piecewise constant α(r) in the equilibrium configuration, the perturbation does not necessarily represent the most unstable mode for the setup of the simulation. For this reason there is a brief transient period before the exponential rise of the instability at t ≈ 10, as shown in Figure 2(a). The isotropic model damps this initial velocity perturbation more than the switching model, leading to a small difference in kinetic energy during the growth of the linear instability, although the growth rate appears to be identical across the two models. The duration of the linear phase is also unaffected by the choice of viscosity model. Initially, the instability occurs in the inner region of twist, r < 0.5, where the magnetic field kinks helically. This section of the magnetic field compresses into the outer region, creating a current sheet along the length of the tube as shown in Figure 2(b). As the field continues to be compressed, it provides a magnetic pressure force that stalls the linear growth. The greater kinetic energy in the switching case leads to greater compression and thus a larger (though not notably stronger) current sheet. After this point, the growth of the kink instability is no longer in the linear phase.
During the transition from the linear to the nonlinear phase, field lines in the current sheet between the regions of inner and outer twist start to reconnect (Figure 2(b) and (c)). This happens sooner in the switching case, due to the larger compression.

Nonlinear phase
Although the choice of viscosity model has a small effect on the linear phase of the kink instability, it does play an important role in the development of the nonlinear phase. By examining the kinetic energies (KEs) in Figures 3(a) and (b), a pattern emerges in both cases that has similarities with the nonlinear behaviour of kink instabilities described in Hood et al. [19]. Shortly after the linear phase, at t ≈ 50, the KEs for both viscosity models exhibits a sharp rise, with the KEs associated with the switching model attaining higher amplitudes. At the same time, a sharp rise is also found in the maximum current as seen in Figure 3(c) and, leading on from this spike, the current magnitudes associated with the switching model are larger than those associated with the isotropic model. Returning to the KEs, the energies associated with the switching model are greater until t ≈ 175, after which, in only the isotropic case, we see a clear secondary spike in perpendicular kinetic energy, and a large increase in parallel kinetic energy, much greater than the energy shown in the switching case. It is difficult to detect this new phase in the maximum current (Figure 3(c)), but it is found in other quantities related to magnetic reconnection.   Figure 4 display similar trends to those found in the perpendicular kinetic energy plot, Figure 3(a). For both isotropic and anisotropic viscosity there appears to be two major peaks in the reconnection measures that align with peaks in the perpendicular kinetic energy. This is much more obvious in the isotropic case. We can conclude that both viscosity models allow for two phases of reconnection but the time at which they occur is significantly modified by the form of viscosity chosen. It is, therefore, clear that the form of viscosity is having a significant effect on the nonlinear evolution of the kink instability, both on the flow dynamics and the reconnection of the magnetic field. We will now examine, in more detail, the two important phases indicated by the isotropic time series, and how the switching case differs.

First phase: t ≈ 65-100
At t = 65, an intense current structure appears near the centre of the tube for both viscosity models, although it is much stronger in the switching case as illustrated in Figure 5. Since the viscous damping associated with parallel viscosity is much less than that of isotropic viscosity, the flows in the switching case are stronger than those in the isotropic case (Figures 3(a) and (b)). The faster flows drive stronger reconnection in the central current structure (see Figure 4) and the interaction of these processes leads to stronger outflows and finer-scale structures in the switching model case compared with the isotropic model case. Evidence of this behaviour can be seen by comparing the current and flow structures in Figure 5. The effects of this phase can also be seen in the magnetic energy evolution, shown in Figure 3(d). Between times t = 100 and 125, due to stronger reconnection in the switching case, the magnetic field relaxes marginally faster than that of the isotropic case, before the secondary instability begins in the isotropic case around t = 125.

Second phase: t ≈ 125-175
The contrast between fine-scale current and flow structures for the switching model, and the smoother, larger-scale structures of the isotropic model continues to be present at later times. Figure 6 shows the same data as Figure 5 but for the times t = 125, 150 and 175. Looking at the slices for t = 125, there is more fine-scale structure generated in the switching case compared to the isotropic case, as in the first phase described above. This second phase, however, marks the beginning of a significant change in behaviour in the isotropic model case. From Figure 3, the parallel KE for the isotropic model exhibits a rapid and large increase in kinetic energy, characteristic of a secondary instability. To a lesser extent, there is also growth in the perpendicular KE, and the two reconnection measures for the isotropic model. In the second phase, these three measures increase to eventually become greater than their corresponding values for the switching model, around t = 175. In order to understand this significant difference in the behaviour between the two models, we will first consider the slices in Figure 6. At t = 125 (panels (a) and (d) of Figure 6), the difference in behaviour between the models is similar to the first phase but some new features appear. The KE in the isotropic model begins to increase and, as mentioned before, appears to signify a secondary instability. In Figure 6(a), two new current sheets have formed at the top and bottom of the tube. A three-dimensional (3D) visualisation of these current sheets is shown in Figure 7(a). The outflow from the reconnection occurring within these current sheets then creates two new symmetric vortices on the right hand side of the tube, advecting the field into the centre of the tube. This behaviour can be seen clearly in Figure 6(b) where vortex motion compresses the magnetic field and forms a central region of enhanced current density. Later, as seen at t = 175 in Figure 6(c), the central current region becomes stronger due to continued compression and reconnection ensues, becoming stronger than the switching model case (see Figure 4). The outflows from this current region then feed into the vortical motions that drive the compression. In this way, a feedback loop is set up, and the reconnection within the current structure continuously drives the flow, resulting in an instability. This kind of interaction between multiple current sheets is also seen in [19]. Due to this secondary instability, magnetic relaxation now becomes faster for the isotropic case. The magnetic energy for this case now dips below that of the switching case, as shown in Figure 3(d).
During this phase, the kinetic energy in the switching model case also increases but to a much smaller extent compared to the isotropic model case. Although the current densities in Figures 6(d)  exhibit finer-scale structure compared to the isotropic case, the magnitude of the current density within the tube becomes weaker with a more uniform profile developing in time. The dominating current sheets are on the edge of the tube, as also indicated in Figure 7(b).

Late-time states
For both cases, the asymptotic relaxed magnetic field is a linear force-free field. The route to this asymptotic state, however, depends on the viscosity model used. At the late time of t = 600, there remain clear differences in the field structure between the two models resulting from the different nonlinear evolutions, as can be seen in Figure 8. At t = 600, the magnetic field in the isotropic case (Figure 8(a)) appears straighter, indicative of more efficient magnetic relaxation. Indeed, Figure 3(d) shows that more energy has been extracted from the field in the isotropic case. At t = 600, the current density and energies (see Figure 3) are still non-zero, so further relaxation is expected. For coronal applications, however, these late times are not as important as the early phases, described above, when the initial and secondary instabilities develop.

Viscous and Ohmic heating
Over the lifetime of the entire instability the switching model allows for the generation of more Ohmic heating (Figure 9(a)). This is despite the long, secondary phase of reconnection produced in the isotropic  case. The greater heating in the switching case is due to two factors: the greater compression created by faster flows, creating stronger or larger current sheets and the more numerous current sheets created by more complex flows. However, isotropic viscous heating dominates that of the switching model by two orders of magnitude (Figure 9(b)) ultimately leading to greater overall heating in the isotropic case (Figure 9(c)). Physically, this is due to anisotropic viscosity only performing significant damping when velocity gradients align appropriately with the magnetic field (that is, when (W b) · b is non-zero). Comparing Ohmic and viscous heating (Figures 9(a) and (b)), Ohmic heating outperforms viscous heating in both cases, by an order of magnitude in the isotropic case and by three orders in the switching case. Even though we use similar values for the diffusion of the magnetic field η and the velocity ν, during the kink instability the current sheets produced are much stronger than the gradients in velocity, hence the Ohmic heating dissipates more energy than the viscous heating.
Due to the relationship between (W b) · b and Q ν (equation (22) with s ≈ 1), the small magnitude of Q ν in Figure 9(b) implies that (W b) · b is small everywhere. With the anisotropic viscous heating being heavily dependent on the magnetic field direction and since (W b) · b is small everywhere in the kink simulation, it follows that the anisotropic viscous heating is always lower in magnitude compared to the isotropic viscous heating, which is not bound by the diection of the magnetic field.

The effect of anisotropy on feedback reconnection
We have described the nonlinear evolution of the kink instability for the cases of only isotropic viscosity and (practically) only anisotropic viscosity. In order to determine how much of each extreme form of viscosity controls the evolution of the secondary instability, we can fix the interpolation between isotropic and anisotropic viscosity by prescribing the switching function s as a constant in equation (14), instead of letting s rely on the local field strength |B|. It should be noted that the simulations in which we fix s are no longer physically realistic, but we use them only to estimate the amount of anisotropy required in the viscosity to disrupt the secondary instability. Since the interpolation involves only s 2 instead of s, in practice we fix the value of s 2 .
By letting s 2 in equation (14) take values between 0 and 1, it is found that there is not a smooth transition between the two extremes of behaviour. Instead, we find a critical value of s 2 , between 0.5 and 0.6, below which (closer to isotropic) we see flows simple enough to create and sustain feedback reconnection, and above which (closer to anisotropic) we see flows complex enough to disrupt this feedback situation. This behaviour can be seen in how the kinetic energy time series changes with s 2 in Figure 10.

Parameter study
In order to confirm that the results of Section 4 are typical, and to further understand how they vary, two parameter studies are performed; one varying viscosity, keeping all other parameters constant; and one varying resistivity, again keeping all other parameters constant.
In the first study we vary the viscosity as ν = 5 × 10 −n , where the index n takes the values 4.75, 4.5, 4.25, 4 and 3.75, while keeping resistivity constant at η = 5 × 10 −4.5 . This range of viscosities represents values that are typically used in simulations, with a lower bound above numerical diffusion and an upper bound below physically unrealistic values for the corona.
In the second study we vary the resistivity as η = 5 × 10 −m , where the index m takes the values 4.75, 4.5, 4.25, 4, 3.75, and 3.5, while keeping the resistivity constant at ν = 5 × 10 −4.5 . Similar to the limits on viscosity, any lower resistivities become comparable to numerical diffusion. Higher resistivities diffuse the field so quickly that the instability does not have time to grow. Figure 11 shows the maximum kinetic energy produced by the two instabilities found in the isotropic case in Section 4. The maximum kinetic energy provides a useful measure of the efficacy of an instability,  Figure 11: Maximum kinetic energy corresponding to initial instability and secondary instability as functions of resistivity η and viscosity ν. (a) resistivity η is fixed at 5 × 10 −4.5 and (b) viscosity ν is fixed at 5 × 10 −4.5 . In both plots are shown the maximum kinetic energy produced by the initial instability (green, 1-marker) and the maximum kinetic energy produced by the secondary instability (purple, 2-marker). Only results using the isotropic viscosity are shown. particularly when comparing the relative magnitudes of the initial and secondary instabilities. Since we only find evidence of the secondary instability in the isotropic case, we do not show the results from the switching case.

Effect on the secondary instability varying diffusion parameters
Looking at Figure 11(a), it is observed that increasing ν reduces the kinetic energy generated in both instabilities. For small values of ν we find the secondary instability causes more energy to be produced than the first, however as ν increases, this relationship reverses, with the initial instability causing more energy to be produced than the secondary one for large ν. This reversal suggests that the greater kinetic energy produced by the initial instability for low values of ν is causing a stronger current sheet to form, enhancing reconnection, and producing a stronger secondary instability.
The effect of resistivity η on the secondary instability is to suppress it entirely when η is large. Since the secondary instability is driven by reconnection outflows, it is not surprising that we can find values of η for which the reconnection outflows do not feedback to produce the instability. show the total heat generated by t = 400 via viscous Q ν and Ohmic Q η dissipation as we vary the strength of viscosity ν. It should be noted that, to allow the trend in the anisotropic viscous heating to be seen in the plot, it has been multiplied by a factor of 10. Before discussing the apparent trends in the heating when we vary viscosity, it is useful to note that, just as in the typical case described previously, for the range of ν shown, isotropic viscous heating remains approximately two orders . This is to capture the behaviour of only the initial nonlinear evolution of the instability, neglecting any further instabilities like the secondary instability found in Section 4. Note, the maxima do not necessarily occur at the same time and this particular parameter study has been performed fixing ν at a slightly different value to the previous parameter studies.
of magnitude greater than the anisotropic viscous heating, and the Ohmic heating is consistently higher when using anisotropic viscosity than when using isotropic.
Since viscous dissipation (equations (21) and (22)) has a functional dependence on ν and Ohmic dissipation (equation (20)) does not, we could naively assume that as we vary viscosity we should see some trend in the viscous dissipation for both models and no trend in the Ohmic dissipation. The trends that are observed broadly adhere to this but, unexpectedly we do find some trend in the Ohmic heating when using isotropic viscosity.
When employing the switching model, the Ohmic heating appears to be independent of ν, whereas when employing the isotropic model, we find a small trend of decreased Ohmic heating with increased ν. These trends can be explained by considering the effect of viscosity on compressive flows and current densities. During the kink instability, Ohmic heating, being proportional to the square of the local current density, is increased when an already sheared magnetic field is compressed by flows perpendicular to the field, increasing the local current density. Thus, as the speeds of perpendicular flows increase, so does the Ohmic heating. These perpendicular flows are effectively only damped by isotropic viscosity. Since we find the maximum kinetic energy (Figure 12(c)) decreases with ν in only the isotropic case, and remains constant in the switching case, it is appropriate that the Ohmic heating decreases with ν in the isotropic case and is negligibly dependent on ν in the switching case.
If varying ν does not change the dynamics in the switching case, the functional dependence of Q ν on ν (see equation (22)) suggests we should observe an increase in anisotropic viscous heating with ν. Figure 12(a) reveals precisely this.
The relationship between the isotropic viscous heating at ν appears non-trivial. Given the decrease in maximum kinetic energy (Figure 12(a)) with ν, it is expected the isotropic viscous heating should also decrease. However, this is not what is observed. Although there appears to be a slight decreasing trend in the isotropic viscous heating when ν is increased past 10 −4 , the left-most point is clearly an outlier. This suggests the secondary instability is having a significant and non-trivial effect on the heating. Indeed this is also suggested by the subtle change of gradient in the maximum kinetic energy on the left-hand side of the Figure 12(c). Due to this particular parameter study producing only five data points, we cannot discuss these trends with much confidence. A more detailed parameter study should be performed, investigating more values of ν within and beyond the range studied here.

Dependence of linear growth rate on viscosity
For each value of η we calculate the linear growth rate σ of the onset of the kink instability by plotting the logarithm of the kinetic energy against time and measuring the gradient during the period of linear growth (as is done in Figure 2). Figure 13(a) plots these growth rates against η, and Figure 13(b) shows the maximum kinetic energy calculated as the maximum prior to t = 125. For every η, this time is between the peaks of the kinetic energy corresponding to the first and secondary instabilities. Taking the maximum before this time allows us to capture only the behaviour of the initial instability, since this is the instability of interest in this section.
It can be seen from the relationship between the growth rate and ν for both viscosity models that isotropic viscosity appears to begin to suppress the kink instability, for larger ν, while the switching viscosity does not (Figure 13(a)). This is also apparent from the relationship between the maximum kinetic energy and ν for both models (Figure 13(b)). This difference between the viscosity models results from the anisotropic viscosity being so weak that the dynamics of the initial onset of the kink instability are not significantly affected by a significant increase in ν. Figure 14(a) and 14(b) show the total heating generated by t = 400 via viscous Q ν and Ohmic Q η dissipation as we vary the strength of resistivity η. Just as in Section 5.2, the anisotropic viscous heating is multiplied by 10. Again, it is useful to note that, across the entire range of η studied here, the isotropic viscous heating remains approximately two orders of magnitude greater than the anisotropic viscous heating and the Ohmic heating produced when using switching viscosity is consistently higher than that produced when using isotropic viscosity. This aligns with the results when varying viscosity, as discussed above.

Dependence of heating on resistivity
As in the parameter study varying ν, we also observe a non-trivial relationship between the viscous heating for both models and η (Figure 14(a)). The isotropic viscous heating reveals an decreasing trend over all values of η studies here, however there is a clear jump in heating between approximately 10 −3.9 and 10 −3.7 . Given that these are the values of η where we observe strong influence of the secondary instability on the kinetic energy output (see Section 5.1), these results suggest it is the kinetic energy produced by the secondary instability that is being damped at low values of η.
Just as in the parameter study varying ν, the anisotropic viscosity shows very little variability with η ( Figure 14(a)), even with the heating multiplied by a factor of 10 for plotting. Despite the dynamics significantly changing with η, the effect of the anisotropic viscosity is so small that we observe very little change in the heating. . This is to capture the behaviour of only the initial nonlinear evolution of the instability, neglecting any further instabilities like the secondary instability found in Section 4. Note, the maxima do not necessarily occur at the same time and this particular parameter study has been performed fixing ν at a slightly different value to the previous parameter studies.
We observe that Ohmic heating increases with increasing η (Figure 14(b)). This is to be expected given the functional dependence of Q η on η, however the actual scaling is not linear in η, as might be predicted from equation (20). Rather, we find that Q ν varies linearly with log 10 (η 1/2 ) for the range of η studied here. Without a more comprehensive parameter study covering more values of η, it is difficult to reason why the scaling takes this form. However, we are able to state with confidence that the use of anisotropic viscosity is consistently enhancing Ohmic heating across the range of η studied here. This is due to the kink instability producing more kinetic energy in the switching case, which better compresses the magnetic field, creating stronger current sheets and thus enhancing Ohmic heating.

Dependence of linear growth rate on resistivity
As is done in Section 5.2.2, we calculate the linear growth rates for each value of η. These and the maximum early time (t < 125) kinetic energy are shown in Figure 15. The plots show that the use of the switching model seems to consistently amplify the growth of the kink instability, shown both in the growth rate and in the kinetic energy. Beyond this, the two models of viscosity show similar trends with η.
Both plots in Figure 15 show that the kink instability is strongly inhibited for values of η greater than approximately 10 −2.5 . This can be explained by the initial diffusion of the magnetic field being so fast-acting for large values of η that the instability is totally suppressed. The increased suppression of the instability with strength of Ohmic diffusion can be seen in both plots as η increases past 10 −3 .

Summary and discussion
We have studied the linear and nonlinear development of the MHD kink instability with two different viscosity models. The first is isotropic (Newtonian) viscosity, which is the most commonly used viscosity model in coronal loop studies. The second is anisotropic viscosity, representing the strong-field limit of Bragkinskii viscosity with a preferred direction parallel to the magnetic field. The implementation of anisotropic viscosity is via the switching model [17] which is suitable for coronal applications.
By considering particular (low) values of the viscosity and resistivity, we find that the effect of the different viscosity models on the linear onset of the kink instability is marginal. The significant differences appear in the nonlinear phase. Two main phases of evolution can be identified which highlight the differences between the effects of the two viscosity models. The anisotropic (switching) case produces more kinetic energy at the onset of the nonlinear phase of the instability-the first phase. It also produces flows and current sheets with smaller length scales compared to the isotropic case and this allows the magnetic field to relax faster due to more efficient reconnection.
In the second phase, the isotropic case exhibits a secondary instability, which is not found in the anisotropic case. This new instability leads to enhanced reconnection and faster magnetic relaxation, compared to the anisotropic case. The simulations are run for 600 Alfvén times (a long time period for coronal applications) and the behaviour of the second phase continues for all of this time.
We have also run a series of parameter studies varying the strengths of viscosity and resistivity. We find the qualitative results of the two phases of the detailed investigation hold true over a range of viscosities and resistivities, including the existence of the secondary instability. Notably, over all parameters studied, viscous heating is consistently overestimated by the isotropic model, and Ohmic heating is consistently enhanced by use of the switching model.
Although there can be much variability in the nonlinear behaviour of the kink instability, our results reveal an important general finding. At the beginning of the nonlinear phase, anisotropic (parallel) viscosity allows for the development of smaller length scales (both flows and current sheets), compared to isotropic viscosity, leading to more efficient reconnection and faster magnetic relaxation, at least initially. As has been described in detail, isotropic viscosity can produce other effects later. However, for coronal applications, it is the initial nonlinear phase of the instability that is likely to be of most interest since, in reality, a coronal loop will interact with others on a longer time scale, thus affecting the nonlinear evolution. Since anisotropic viscosity is a more realistic model for viscosity in the corona, our results will be useful in the interpretation of observations of coronal loops that are kink unstable. This topic will be an interesting avenue of future research.