© ESO, 2017

1. Introduction

With the aid of improved spatial and temporal resolution in contemporary observing instruments, many recent studies have highlighted the abundance of wave energy in the solar corona (Aschwanden et al. 1999; Okamoto et al. 2007; Tomczyk et al. 2007; Parnell & De Moortel 2012). Whether this energy is dissipated sufficiently quickly to significantly contribute to the coronal heating budget remains a contentious issue (e.g. Arregui 2015; Cargill et al. 2016). However, the potential for waves to drive turbulent flows may have consequences for heating the corona as the small length scales that develop can lead to an increased rate of wave energy dissipation (see e.g. recent work by van Ballegooijen et al. 2011; Woolsey & Cranmer 2015).

The rapid damping of transverse kink oscillations in coronal loops has been widely detected in recent years. A number of studies have reported this phenomenon in both standing (e.g. Aschwanden et al. 2002; Pascoe et al. 2016) and propagating waves (e.g. Verth et al. 2010; Morton et al. 2014), and find decay rates far in excess of those expected from typical coronal dissipation. It is widely anticipated that the cause of this enhanced damping is the process of resonant absorption (or mode coupling in the case of propagating waves) through which energy is efficiently transferred from the kink mode to localised Alfvénic waves (Ionson 1978). This transfer is caused by the existence of a resonant layer in the loop on which the local Alfvén frequency corresponds to the frequency of the transverse wave mode (see Goossens & Ruderman 2011, for a comprehensive review).

This process may be of relevance to the study of coronal heating since the resultant Alfvénic waves are subject to dissipation through the process of phase mixing (Pagano & De Moortel 2017). The azimuthal waves are contained within the loop’s boundary layer and are associated with large gradients in both the velocity and magnetic fields. The large velocity shear is of particular interest for this paper as it may drive the formation of the Kelvin-Helmholtz instability (e.g. Heyvaerts & Priest 1983). Previous analytical studies (e.g. Browning & Priest 1984; Hollweg et al. 1990; Uchimoto et al. 1991) and numerical models (e.g. Terradas et al. 2008a; Antolin et al. 2014; Magyar & Van Doorsselaere 2016; Howson et al. 2017) have demonstrated that the KHI may be triggered in oscillating coronal loops. This may have further significance for coronal heating as the formation of the instability generates the possibility of an enhanced rate of wave energy dissipation. In particular, as the KHI develops, a regime of magnetohydrodynamic (MHD) turbulence develops, leading to a cascade of energy to smaller spatial scales and, ultimately, heating.

Whether the KHI is common in the solar atmosphere has not yet been conclusively established. Indeed, the Alfvénic waves which may trigger the instability are difficult to detect due to their weak compressibility. Furthermore, they are associated with length scales that reach the spatial resolving power of current observational instruments. Despite this, in the past few years, Foullon et al. (2011), Ofman & Thompson (2011) claim that the appearance of vortices in observations of coronal mass ejections may be evidence of the KHI. In these studies, however, the triggering mechanism is the presence of large plasma velocities parallel to the magnetic field, rather than azimuthal waves. Whilst no direct observations of the KHI forming in transversely oscillating coronal loops have been detected, Okamoto et al. (2015), Antolin et al. (2015) interpret plasma behaviour in prominence observations as evidence of resonant absorption and a turbulent regime induced by the KHI.

Several studies have demonstrated that resonant absorption remains important in various coronal loop configurations. For example, Van Doorsselaere et al. (2004), Terradas et al. (2006) find that the process is robust in the context of curved loops and Ruderman (2003) considers the process in the case of loops with elliptical cross-sections. Furthermore, studies considering general loop shapes (Pascoe et al. 2011) and multi-stranded structuring (Terradas et al. 2008b) highlight the insensitivity of resonant absorption (mode coupling) to these conditions. Additionally, studies find that loop stratification due to gravitational effects (Andries et al. 2005) and loops containing magnetic twist (Karami & Barin 2009; Ebrahimi & Karami 2016) show that these configurations do not prevent the transfer of energy to azimuthal Alfvénic modes. Comprehensive reviews of resonant absorption in various geometries are presented in Ruderman & Erdélyi (2009) and Goossens & Ruderman (2011).

Many studies have reported the presence of magnetic field twist in coronal structures that can be observed in prominences (and their subsequent eruption, e.g. Rust & Kumar 1996), in active regions (Démoulin et al. 2002) or predicted through numerical simulations of photospheric motions (Török & Kliem 2003) and flux emergence (e.g. Magara & Longcope 2003; Fan 2009). Indeed, shear flows at the photosphere are expected to introduce some twist into the corona. As a result, investigations into the effects of magnetic twist on the decay of kink mode oscillations are very relevant. It has been suggested that the presence of twist in coronal loops will stabilise the velocity shear through the magnetic tension force and so prevent the formation of the KHI (Soler et al. 2010). On the other hand, weak values of twist may still permit some development of the instability (Terradas, in prep.) and Zaqarashvili et al. (2015) find that standing kink and torsional Alfvén waves are always unstable to the KHI close to the antinode at the loop apex.

We aim to quantify the effect of twist on the non-linear evolution of kink modes by conducting a parameter study on the azimuthal component of the magnetic field. We introduce the numerical model used in Sect. 2 and describe our results in Sect. 3. A discussion and some conclusions are then presented in Sect. 4.

2. Numerical method

2.1. Initial configuration

We followed the work of Antolin et al. (2014), Howson et al. (2017) and modelled a coronal loop as a straight flux tube that was denser than the surrounding plasma. The loop was uniform along its length and had a circular cross-section with a dense, inner region (the core) and a boundary layer (the shell) over which the density increased from the external density, ρ_e = 8.4 × 10^-13 kg m^-3 to the internal density ρ_i = 3ρ_e. At all heights along the loop, the density was defined as $\mathit{\rho }\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}{\mathit{\rho }}_{\mathrm{e}}\mathrm{+}\frac{\mathrm{(}{\mathit{\rho }}_{\mathrm{i}}\mathrm{-}{\mathit{\rho }}_{\mathrm{e}}\mathrm{)}}{\mathrm{2}}\left(\mathrm{1}\mathrm{-}\tanh \left\{\frac{\mathit{R}\mathrm{-}{\mathit{r}}_{\mathit{a}}}{{\mathit{r}}_{\mathit{b}}}\right\}\right)\mathit{,}$(1)where r_a = 1 Mm is the loop radius measured to the centre of the boundary layer and ${\mathit{r}}_{\mathit{b}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{16}}$ Mm. This ensures the width of the boundary layer is approximately 0.4 Mm. The corresponding density profile across the loop is shown in Fig. 2. Measured to the centre of the boundary layer, the radius of the loop was 1 Mm and the length of the loop was 200 Mm. The initial set-up is shown in Fig. 1.

Fig. 1 Configuration of the modelled loop.

Three different cases were considered for the nature of the magnetic field. In the following we describe the set up using cylindrical co-ordinates with the origin at the centre of the loop apex.

• Case 1: straight field – a uniform field of strength21 G aligned with the loop axis (y direction) wasincluded throughout the domain. In this case, in order to ensure thetotal pressure was uniform, the initial temperature was defined as$\mathit{T}\mathrm{=}\frac{{\mathit{P}}_{\mathrm{0}}}{\mathit{\rho }}\mathit{,}$(2)where P₀ was set to ensure the plasma-β = 0.05.

• Case 2: twisted field; | B_φ | greatest at core-shell boundary – as with the previous case, a magnetic field of strength 21 G was included in the domain. In this case the azimuthal component of the field was defined to be ${\mathit{B}}_{\mathit{\phi }}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ & \mathrm{if}\mathit{R}\mathrm{\le }\mathit{a,}\\ & \mathrm{if}\mathit{a}\mathit{<}\mathit{R}\mathrm{\le }\mathit{b,}\\ & \mathrm{if}\mathit{R}\mathit{>}\mathit{b}\mathit{.}\end{array}$(3)Here, τ is a parameter quantifying the level of twist, a is the radius of the core region and b is the radius of the loop (including the boundary region). The profile of B_φ is constant with height and we note that if τ = 0, the first case is recovered. Three cases of twist were considered, τ = {0.1,1.0,2.0} and are labelled as Case 2a, 2b and 2c, respectively. The amount of twist that these values produce is presented in Table 1.

  Table 1

  Twist parameter space.

  Fig. 2 Initial density profile for all cases through the cross-section of the loop. The density is normalised to the initial exterior density, ρe = 8.4 × 10-13 kg m-3.

  In order to maintain a constant field strength, we stipulated that the vertical magnetic field is ${\mathit{B}}_{\mathit{y}}^{\mathrm{2}}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{-}{\mathit{B}}_{\mathit{\phi }}^{\mathrm{2}}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathit{,}$(4)where B₀ = 21 G. In this geometry, the radially inwards tension force, T_R is given by ${\mathit{T}}_{\mathit{R}}\mathrm{=}\frac{{\mathit{B}}_{\mathit{\phi }}^{\mathrm{2}}}{\mathit{R}}\mathrm{·}$(5)Therefore, in order to balance the force associated with the magnetic twist, the gas pressure was defined as $\mathit{P}\mathrm{=}\mathrm{\int }\frac{{\mathit{B}}_{\mathit{\phi }}^{\mathrm{2}}}{\mathit{R}}\hspace{0.17em}\mathrm{d}\mathit{R,}$(6)and the constant of integration was selected to ensure the plasma-β = 0.05.

• Case 3: twisted field; | B_φ | greatest at shell-exterior boundary – once again, a magnetic field of strength 21 G was included throughout the domain. In this case the azimuthal component of the field was defined to be ${\mathit{B}}_{\mathit{\phi }}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ & \mathrm{if}\mathit{R}\mathrm{\le }\mathit{b,}\\ & \mathrm{if}\mathit{a}\mathit{<}\mathit{R}\mathrm{\le }\mathrm{2}\mathit{b,}\\ & \mathrm{if}\mathit{R}\mathit{>}\mathrm{2}\mathit{b}\mathit{.}\end{array}$(7)Here, the twist extends beyond the density-enhancement of the loop. This case was included in order to determine whether the location of twist within a loop affected the development of the KHI. Only ψ = 1.0 was considered and the corresponding level of twist is shown in Table 1. Once again, in order to maintain an initial equilibrium, the vertical magnetic field was given by Eq. (4) and the gas pressure was defined as in Eq. (6).

Since in each of the three cases, B_R = 0 and there is azimuthal and vertical invariance, we immediately see that ∇·B = 0 is satisfied.

Fig. 3 Magnitude of initial azimuthal magnetic field, | Bφ |, through the cross-section of the loop for each case. Here, Bφ is normalised to B0 = 21 G. We see that in all cases, | Bφ | ≪ | B0 |.

Fig. 4 Initial gas pressure profile through the cross-section of the loop for each case. Here, we have normalised the pressure to the external pressure, P0.

We quantified the amount of twist by determining the angle a field line moves through between the lower and upper boundaries of the domain. The cases considered are displayed in Table 1. The initial profiles of the azimuthal field strength and the pressure are shown in Figs. 3 and 4, respectively. We see that even in the most twisted configurations, the gas pressure inside the loop is always above 98% of the external gas pressure. In all cases the temperature profiles are very similar and decrease from an external temperature of 2.5 MK to an internal temperature of around 0.8 MK. Such a profile may be consistent with a loop in thermal non-equilibrium during its cooling phase.

Following Howson et al. (2017; which we will henceforth refer to as Paper I), in each case, the fundamental kink mode was excited using the following velocity profile $\begin{array}{ccc}& & {\mathit{v}}_{\mathit{x}}\mathrm{=}\mathit{A}\left(\mathit{\rho }\mathrm{-}{\mathit{\rho }}_{\mathrm{e}}\right)\cos \left(\frac{\mathit{\pi y}}{\mathrm{2}\mathit{L}}\right)\mathit{,}\\ & & {\mathit{v}}_{\mathit{y}}\mathrm{=}\mathrm{0}\mathit{,}\\ & & {\mathit{v}}_{\mathit{z}}\mathrm{=}\mathrm{0.}\end{array}$(8)Here, L = 200 Mm is the length of the modelled loop and the constant, A, is selected to produce a maximum velocity at the loop apex (y = 0) of 8.3 kms^-1. This oscillation corresponds to commonly observed low amplitude kink modes (e.g. Anfinogentov et al. 2015).

2.2. Numerical method

In each case, we numerically advanced the full, 3D, resistive MHD equations in normalised form using the Lagrangian-remap code, Lare3d (Arber et al. 2001). In order to avoid numerical instabilities, the code uses a staggered grid and implements an adaptive time step restriction so as not to violate the CFL condition whilst maintaining calculation efficiency. The equations are given by $\begin{array}{ccc}& & \frac{\mathrm{D}\mathit{\rho }}{\mathrm{D}\mathit{t}}\mathrm{=}\mathrm{-}\mathit{\rho }\mathrm{\nabla }\mathrm{·}{v}\mathit{,}\\ & & \\ & & \\ & & \end{array}$Here, variables have their usual meanings and the non-ideal contributions arise from the inclusion of viscous forces, F_ν, in the equation of motion (10), the associated heating term H_ν in the energy Eq. (11) and the resistivity (η)-dependent terms in (11) and the induction Eq. (12). The viscous force, F_ν, and viscous heating term, H_ν, are given by $\begin{array}{ccc}& & {{F}}_{\mathit{\nu }}\mathrm{=}\mathit{\nu }\left({\mathrm{\nabla }}^{\mathrm{2}}{v}\mathrm{+}\frac{\mathrm{1}}{\mathrm{3}}\mathrm{\nabla }\left(\mathrm{\nabla }\mathrm{·}{v}\right)\right)\mathit{,}\\ & & {\mathit{H}}_{\mathit{\nu }}\mathrm{=}\mathit{\nu }\left(\frac{\mathrm{1}}{\mathrm{4}}{\mathit{ϵ}}_{\mathit{i,j}}^{\mathrm{2}}\mathrm{-}\frac{\mathrm{2}}{\mathrm{3}}{\left(\mathrm{\nabla }\mathrm{·}{v}\right)}^{\mathrm{2}}\right)\mathit{,}\end{array}$where $\begin{array}{ccc}{\mathit{ϵ}}_{\mathit{i,j}}\mathrm{=}\frac{\mathit{\partial }{\mathit{v}}_{\mathit{i}}}{\mathit{\partial }{\mathit{x}}_{\mathit{j}}}\mathrm{+}\frac{\mathit{\partial }{\mathit{v}}_{\mathit{j}}}{\mathit{\partial }{\mathit{x}}_{\mathit{i}}}\mathrm{·}& & \end{array}$For comparison purposes, the domain used was identical to the high resolution case presented in Paper I. In all cases, the numerical domain had dimensions of 64 Mm × 200 Mm × 64 Mm and used 512 × 100 × 512 grid points. The loop was aligned with the y-direction and the initial velocity pertubation was aligned with the x-direction. The origin (x = y = z = 0) was located at the centre of the loop apex.

In order to ensure a standing mode was generated, the velocities were forced to zero on the y boundaries. All other variables have zero gradients on the edges of the domain. We note that since the magnetic twist is weak (B_φ ≪ B_y) in all of our simulations, this is sufficient to ensure a standing mode is generated. However, in cases with more significant twist, it may become necessary to use line-tying boundary conditions.

In order to minimise the boundary effects, a non-uniform grid was implemented, with decreasing spatial resolution in the x−z-plane away from the centre of the domain. To limit any effects of this non-uniformity, the grid was designed so that the oscillation and any subsequent fine scale dynamics occurred within a uniform high resolution region. The cell resolution in the centre of the domain corresponded to 15.9 km × 2000 km × 15.9 km. The y-axis had a uniform spatial resolution.

All simulations were essentially ideal since the values of resistivity, η = 10^-20 and viscosity, ν = 10^-20, are much smaller than the estimated numerical values. For the spatial resolution used, the numerical (and effective) values of η and ν are of the order 10^-5−10^-6 as discussed in Paper I.

3. Results

In all cases, the initial velocity profile induces the fundamental, standing kink mode with a maximum displacement at the loop apex of approximately 0.27 Mm (less than the loop radius of ≈1 Mm) and a period of approximately 280 s. The maximum displacement at each height is proportional to the maximum initial velocity at that height.

Fig. 5 Displacement of the loop apex through time for the straight case and a twisted case (Case 2b). The red lines are given by $\mathit{y}\mathrm{=}\mathrm{\pm }\mathrm{0.267}{\mathrm{e}}^{\mathrm{-}{\left(\frac{\mathit{t}}{\mathrm{1040}}\right)}^{\mathrm{2}}}$.

3.1. Resonant absorption

The displacement of the loop apex as a function of time for both the straight case and a twisted case (2b) is shown in Fig. 5. We observe that the damping rate is largely independent of the twist present and follows a Gaussian profile (red lines). This is in accordance with previous studies that have explained such a (initial) damping rate for both standing (Ruderman & Terradas 2013) and propagating (Hood et al. 2013) waves. This damping is an ideal process and is caused by the transfer of energy from the kink wave to azimuthal, Alfvénic modes contained within the loop’s boundary region through the well studied process of resonant absorption.

A resonance exists within the loop’s shell region at locations at which the kink frequency coincides with the local Alfvén frequency. In the straight case this is the location on which the density, $\mathit{\rho }\mathrm{=}\frac{{\mathit{\rho }}_{\mathrm{i}}\mathrm{+}{\mathit{\rho }}_{\mathrm{e}}}{\mathrm{2}}\mathrm{=}\mathrm{2}$. However, in the twisted cases, the increased length of the magnetic field lines causes the Alfvén frequency to be lower at this location. Consequently, the resonant layer moves to a new location further from the loop centre. However, the existence of any resonant layer is sufficient for the transfer of energy to continue (Terradas et al. 2008b; Pascoe et al. 2011) and so there is little change to the actual damping profile.

In Fig. 6 we highlight the concentration of flows in the boundary layer as the kink mode decays. We show the evolution of v² in the line x = xc,y = 0 where xc(t) is the location of the loop centre. This is a line perpendicular to the initial velocity in the horizontal cross-section at the loop apex. The observed oscillatory pattern arises from the periodic transfer of magnetic to kinetic energy associated with the kink wave. The nature of the initial velocity profile dictates that the majority of the kinetic energy is within the core region of the loop. However, as the kink mode couples with the azimuthal wave modes, energy is transferred into the loop’s shell region (approximately 0.8 Mm and 1.2 Mm from the loop centre). As this mode conversion proceeds, the oscillatory pattern survives; however, the length scales of the velocity field are now much smaller than previously. In non-ideal simulations, this would increase the rate of viscous dissipation but can also cause the loop to become Kelvin-Helmholtz unstable (see below). Although Fig. 6 displays the effects of resonant absorption for a twisted loop, a similar figure produced for the straight case is almost identical.

Fig. 6 Evolution of the square of the velocity along the line y = 0 (through the loop apex), z = 0. Here v2 is normalised to the maximum value obtained during the simulation. The figure is produced from Case 2b, a twisted field simulation.

Whilst the quantity of energy transferred from the kink mode to the Alfvénic modes is largely independent of the magnetic twist, the presence of an azimuthal component in the field excites additional harmonics of the azimuthal Alfvénic waves. In particular, the magnetic twist enhances the non-linearity of the system which is known to produce fluting modes (Ruderman et al. 2010; Ruderman & Goossens 2014; Terradas, in prep.). In Fig. 7, we show the time variation of v_φ along the line $\begin{array}{ccc}& & \mathit{x}\mathrm{=}{\mathit{c}}_{\mathit{x}}\mathrm{\left(}\mathit{y,t}\mathrm{\right)}\mathit{,}\\ & & \mathit{z}\mathrm{=}{\mathit{c}}_{\mathit{z}}\mathrm{-}\mathit{\zeta ,}\end{array}$where c_x(y,t) and c_z are the x and z coordinates of the loop centre, respectively and ζ = 1 Mm. This line moves with the global kink mode and is always contained within the loop boundary. It lies on the resonant layer in the straight case where the amplitude of the Alfvénic waves is expected to be greatest. For each height, we have normalised v_φ by the maximum initial velocity at that height.

Fig. 7 Time evolution of the azimuthal component of the velocity vφ, along a line parallel to the loop axis within the shell region. The velocity is normalised to the maximum initial velocity at that height. The upper panel shows the simulation with a straight magnetic field (Case 1) and the lower panel shows Case 2c.

Fig. 8 Density deformation at the loop apex during the formation of the KHI. The density is normalised to the initial exterior density ρe.

In both cases we observe an oscillatory pattern in v_φ throughout the length of the loop that corresponds to the periodic reversal of the azimuthal waves. The period of these waves match the period of the global kink mode (≈280 s). Further, for both the straight and twisted simulations, we observe an increase in the amplitude of these waves through time as the resonant absorption injects energy into the Alfvénic modes. In the straight case, we notice that once the initial velocity variation has been accounted for, the energy transfer is roughly constant with height. However, the twisted case displays greater variation along the loop axis and, whilst the model is different, shows evidence of the helical fluting modes found in, for example, Ruderman et al. (2010). We note that whilst the vertical uniformity is violated in the twisted case, some symmetry is still present within the loop. In particular, if the plot is produced for the line x = c_x,z = c_z + ζ, then we reproduce the lower panel of Fig. 7 except with a reflection in y = 0 and a change of sign for v_φ.

3.2. Subsequent dynamics and the growth of the KHI

The large velocity shear that develops within the boundary region of the loop (see Fig. 6), becomes unstable to the KHI. The instability is readily observed in the density profile at all heights in the loop, but due to the greatest initial velocity, it is most apparent at the loop apex. In Fig. 8 we show the density cross-section at four different stages in the development of the KHI for all of the numerical simulations.

With the exception of Case 2c, in all of the experiments, Kelvin-Helmholtz vortices begin to form at around t = 600 s (t ≈ 680 s for Case 2c). For a while, three coherent vortices are observed on both sides of the loop (m = 3, see Paper I and Eq. (1) in Terradas et al. 2008a). Beyond this time, further vortices develop as the instability progresses and ultimately a turbulent regime with small spatial scales is generated. Prior to the onset of the KHI, the twisted cases exhibit only minor differences to the straight case. In particular, these differences would be extremely difficult to detect observationally. However, during the growth of the instability, the azimuthal field begins to have a more profound effect on the evolution of the loop.

Despite this, in terms of the density deformation, Case 2a displays almost identical behaviour to that of the straight case which may be expected due to the extremely weak twist present $\left(\frac{\mathrm{\left|}{\mathit{B}}_{\mathit{\phi }}\mathrm{\right|}}{\mathrm{\left|}{\mathit{B}}_{\mathit{y}}\mathrm{\right|}}\mathrm{\le }{\mathrm{10}}^{-3}\right)$ within the domain. Increasing the size of the azimuthal component of the magnetic field clearly suppresses the vortex formation. The cause of this inhibition is the magnetic tension force (Soler et al. 2010). Since resonant absorption proceeds in a similar fashion at all heights, the KHI tends to act in an homogeneous manner along the length of the loop. Certainly, observing the density cross-section a few grid points up (in the y-direction), looks similar to the panels shown in Fig. 8. Since the magnetic field is approximately frozen into the plasma, this ensures that straight field lines remain locally straight. On the other hand, since the density deformation is not azimuthally constant, twisted field lines experience different KHI effects at different heights. This causes field lines to bend and experience a tension force that acts to suppress further development of the instability.

In addition we observe that the location of twist within the loop can affect the growth of the instability. In particular, Case 3 has the same maximum twist as Case 2c but it is located further from the loop centre. Indeed, this case has a higher quantity of twist when integrated across the entire loop. Despite this, we see a lower level of vortex suppression in the Case 3 simulation than in Case 2c. However, once the first KH vortices have formed, there is evidence that the greater radial extent of the twist restricts the development of a turbulent regime at the edge of the shell region (compare the small scales observed in Row 1 Col. 3 with those in Row 5 Col. 3). Whilst the vortices that extend into the region of high twist experience an enhanced tension force that arrests the formation of smaller structures, there is little suppression of turbulence closer to the core of the loop.

Fig. 9 Greatest distance of plasma with density, ρ = 2 from the centre of the loop apex during the formation of the KHI. In the straight field case, initially this plasma is located on the resonance layer.

Fig. 10 Mean density along loop radii in the horizontal cross-section at the loop apex at t = 990 s. By this time, the KHI is well-developed in all cases. For comparison, the initial density profile is also included (black dashed line).

In Paper I, we established several metrics for determining the onset time of the KHI and measuring its subsequent growth. In order to quantify the extent of suppression caused by the magnetic twist, we examine these in Figs. 9, 10. A further comparison of the vorticity growth that was considered in Paper I is addressed in Sect. 3.4. In each figure we observe that the twist either postpones the onset of the KHI (Case 2c) or inhibits its subsequent growth. In the case of weak twist, it is only once the density deformation begins that the suppressive effects of the tension force arise.

In Fig. 9, we aim to summarise the panels of Fig. 8. Tracking the greatest distance of plasma of density, ρ = 2 (this corresponds to the location of the resonant layer in the straight case) from the centre of the loop allows the magnitude of the density deformation to be quantified. The suppressive effects of the magnetic twist are clear and again we note that Case 3 exhibits much lower suppression than Case 2c. Further, since the plasma being tracked is at the location of the resonant layer, we expect that any subsequent energy transfer into the Alfvénic modes will become less azimuthally uniform (Ruderman et al. 2010).

In Paper I, we observed that the KHI reduced the mean radial density gradient across the loop boundary (i.e. smoothed the density transition). In particular, simulations in which the KHI was partially or totally suppressed displayed a density profile which remained very similar to the initial density profile. In Fig. 10, we again observe this as the Case 2 simulations exhibit successively less smoothing as the twist was increased.

Fig. 11 Isosurface of density, ρ = 2 for the straight field case (left-hand panel) and Case 2c, a twisted field simulation (right-hand panel). The time shown is 990 s (three and a half wave periods) after the start of the simulation when the KHI is well-developed.

In Fig. 11, we highlight the density deformation caused by the KHI along the length of the loop. In the straight field case (left), the vortices observed are similar at all heights: since the velocity amplitude is not constant along the loop axis, the magnitude of the density deformation is largest at the antinode (loop apex), however, since the velocity shear has the same form (it is proportional to the maximum initial speed at that height), the KH-vortices are similar along the length of the loop. In contrast, the twisted magnetic field case (right) displays helical density deformation that is much smaller in magnitude than that observed in the straight case. This is in good agreement with Fig. 8.

3.3. Currents and field-aligned flows

One major effect of including magnetic twist within the oscillating loop is the formation of larger currents. In the model set-up, the azimuthal component of the magnetic field generates a field-aligned current which is in turn associated with a radially inwards tension force. This is initially balanced by adjusting the gas pressure accordingly. During the experiments, the resonant absorption process results in the formation of small scale currents within the loop’s boundary layer. In non-ideal simulations these are subject to dissipation which can heat the plasma (e.g. Pagano & De Moortel 2017, or see Paper I). Subsequently, the onset of the KHI generates even smaller scales in the magnetic field and hence larger currents, encouraging further heating (when larger dissipation coefficients are included).

Fig. 12 Magnitude of the current in the horizontal plane at the loop apex for the straight field case. Currents less than a fifth of the maximum size of the current in this plane are coloured in blue. By this stage in the simulation (t = 990 s), the KHI is well-developed.

In Fig. 12, we display the magnitude of the current in the loop apex (this is related to the Ohmic heating rate, ηj²) once the KHI is well-developed (at t = 990 s). Intricate currents form throughout the boundary layer of the loop as the magnetic field is distorted by the Kelvin-Helmholtz vortices. The small length scales in the magnetic field, and hence the current, are generated by the compressive flows caused by the instability (Lapenta & Knoll 2003). As such, the vortical structures produced (see Fig. 8) can also be observed in Fig. 12. Any Ohmic heating associated with these currents will be confined to the boundary region where the turbulent flows form in these simulations. However, during the development of the KHI in high amplitude regimes (see Magyar & Van Doorsselaere 2016), the loop can deform sufficiently to heat the core of the loop through mixing or directly through Ohmic heating as currents begin to develop throughout the cross-section of the loop. Consequently, unlike with resonant absorption induced phase mixing, the core of the loop can be heated in this manner.

The currents at the loop apex are dominated by the horizontal components j_x and j_z for the straight field case. However, this is not true for the twisted field cases. In Fig. 13, we display the maximum magnitude of the current (and each component) at the loop apex for each case as the simulations evolve. The currents in the top left-hand panel are normalised to the pre-KHI formation currents in each simulation to allow comparison of current growth factors. From the bottom left-hand panel, it is clear that the twisted cases generate much larger loop-aligned currents than the straight magnetic field case. Meanwhile, with the exception of | j_z | in Case 3 (discussed below), the horizontal currents remain roughly similar in all cases.

Fig. 13 Formation of large currents during the onset of the KHI. The top-left panel displays the maximum magnitude of the current in the horizontal cross-section at the loop apex. For each simulation, it is normalised to the maximum current in this plane at t = 560 s, immediately prior to the onset of the KHI, to show the relative growth of the currents during the development of the instability. The other three panels display the maximum magnitude of each component in the same cross-section. The units are dimensionless but are the same for each of the components.

Fig. 14 Horizontal components of the magnetic field (vectors) and the deformed density profile (coloured contour plot) at the loop apex once the KHI is well-developed. The result shown is produced from Case 2b, a twisted field simulation.

Since μj = ∇ × B, the y-component of the current is associated with radial gradients in B_φ and azimuthal gradients in B_R. In the straight case, at the loop apex, there is no horizontal field component (with the exception of small variations induced by the KHI) and so | j_y | is small. In the twisted cases, on the other hand, B_φ(R) is not constant and so there is a loop-aligned component present even prior to the development of the instability (t< 600 s; Fig. 13, bottom left-hand panel). Since the size of the radial gradient is dependent on the magnitude of the twist factor, we see that, for example, Case 2c has larger currents than Case 2b.

Following the onset of the KHI, the difference is enhanced by the vortical flows. In Fig. 14, we display the horizontal components of the magnetic field along with the deformed density profile for Case 2b. Considering a plasma element that is moving radially outwards on account of the KHI, we see that it advects field with a large azimuthal component into a region where there is lower twist. This process decreases the length scales on which the azimuthal component varies and so increases the currents. Further twisting motions associated with the KHI introduce radial components of the field and generates further currents.

Despite the suppressive effects of the azimuthal magnetic field on the density deformation, the enhanced currents observed will result in greater KHI-related heating for twisted regimes through Ohmic dissipation. Indeed, the third panel of Fig. 13 does demonstrate that the KHI in twisted cases is more energetic than in the straight case. However, as we see by comparing the Case 2 simulations in the top left hand panel of Fig. 13, increasing the strength of the magnetic twist produces smaller relative increases (from the pre-KHI level) in the magnitude of the current. We also note that despite the very small amount of magnetic twist present in Case 2a, there are still significant loop-aligned currents generated. This suggests that even though the field in Case 2a is locally almost straight, it may produce different heating signatures than the truly straight field case.

As mentioned above, Case 3 produces currents with a larger component in the z direction than the other cases. A later peak is also observed in j_x for Case 3 (top right-hand panel) which is larger than the x components of the currents generated in the other simulations. The horizontal components of the current reach much larger values in this case because the KH-vortices form in regions of stronger magnetic twist than in the Case 2 simulations. The transverse gradients produced by the KHI are therefore much greater than those that form in the other simulations.

One further consequence of the magnetic twist is the excitation of loop-aligned velocities. At all times in the straight case, any vertical flows present at the loop apex are extremely small. Since the gas and magnetic pressures are vertically uniform, there is no y-component of the magnetic tension and gravity is neglected, there are no forces to drive vertical flows at the loop apex. However, once twist in the magnetic field is introduced, the radial currents that form during the growth of the KHI, produce a vertical component of the Lorentz force. This is the case even for small amounts of magnetic field twist. Since these currents are highly spatially variable (some are orientated radially outwards and some are orientated radially inwards), both upwards and downward flows are generated.

In Fig. 15, we show the maximum value of v_y at the loop apex for all the simulations. We note that although Case 2a does not seem too different from the straight case, the velocities observed in the weakly twisted case are many orders of magnitude larger than those in the straight case. In all of the twisted cases, an increase in vertical flows is observed during the resonant absorption phase (t< 600 s), as currents forming in the shell region interact with the azimuthal magnetic field, generating a vertical Lorentz force. This effect becomes much more apparent as the KHI forms and larger currents are generated.

Case 3 displays significantly larger vertical flows then the Case 2 simulations due to the magnetic twist, and hence vertical Lorentz forces, being larger in the shell region where the largest currents form. Additionally, since the horizontal currents in the loop apex are greater in Case 3 (see upper right-hand panel of Fig. 13), this effect is further enhanced. Even prior to the formation of the KHI, we observe some loop-aligned flows being produced due to the interaction of the azimuthal magnetic field and the currents associated with the resonant absorption. Once the KHI forms, the suppressive effect of the magnetic twist ensures that for a short period the vertical flows are greater in Case 2b than Case 2c despite the stronger twist component present in the latter simulation.

3.4. Vorticity

Fig. 15 Maximum value of |vy| at the loop apex against time for all simulations.

In addition to the small scales that form in the magnetic field during the development of the KHI, small scales in the velocity field are also generated (e.g. Paper I). In non-ideal simulations, these gradients are susceptible to the effects of viscosity which, in turn, lead to the dissipation of wave energy. During resonant absorption, the excitation of azimuthal Alfvénic waves is associated with an increase in the loop-aligned component of vorticity. Following the subsequent growth of the KH vortices, the magnitude of the vorticity increases further.

In Fig. 16, we display the magnitude of the vorticity, | ω | in the horizontal cross-section at the loop apex. We show two different stages during the simulation for the straight field case (top row) and the twisted field Case 2c (bottom row). Prior to the formation of the instability, the profiles of | ω | are very similar. Once the KHI forms, differences become apparent as the characteristic vortices are suppressed in the twisted field case. However, the maximum size of the vorticity is very similar in both of these simulations, suggesting that despite the limited spatial extent of the vortices, the amount of kinetic energy is similar.

In the top row of Fig. 17, we display the three components of the vorticity integrated over the loop cross-section. In particular, at each time we calculate, $\begin{array}{ccc}& & \int \mathit{A}{\mathit{\omega }}_{\mathit{x}}\hspace{0.17em}\mathrm{d}\mathit{A}\mathrm{=}{\mathrm{\int }}_{\mathit{A}}\frac{\mathit{\partial }{\mathit{v}}_{\mathit{z}}}{\mathit{\partial y}}\mathrm{-}\frac{\mathit{\partial }{\mathit{v}}_{\mathit{y}}}{\mathit{\partial z}}\hspace{0.17em}\mathrm{d}\mathit{A,}\\ & & \int \mathit{A}{\mathit{\omega }}_{\mathit{y}}\hspace{0.17em}\mathrm{d}\mathit{A}\mathrm{=}{\mathrm{\int }}_{\mathit{A}}\frac{\mathit{\partial }{\mathit{v}}_{\mathit{x}}}{\mathit{\partial z}}\mathrm{-}\frac{\mathit{\partial }{\mathit{v}}_{\mathit{z}}}{\mathit{\partial x}}\hspace{0.17em}\mathrm{d}\mathit{A,}\\ & & {\mathrm{\int }}_{\mathit{A}}{\mathit{\omega }}_{\mathit{z}}\hspace{0.17em}\mathrm{d}\mathit{A}\mathrm{=}{\mathrm{\int }}_{\mathit{A}}\frac{\mathit{\partial }{\mathit{v}}_{\mathit{y}}}{\mathit{\partial x}}\mathrm{-}\frac{\mathit{\partial }{\mathit{v}}_{\mathit{x}}}{\mathit{\partial y}}\hspace{0.17em}\mathrm{d}\mathit{A,}\end{array}$for the panels from left to right. Here the area A is the horizontal cross-section at the loop apex. For all three components of ω, the straight field simulation (black lines) displays very little net vorticity when integrated across the loop cross-section. On the other hand, the inclusion of magnetic twist permits the non-linear generation of additional harmonics that are associated with the oscillating profiles in the panels in the top row of Fig. 17.

The most significant effect is observed in ω_y (second column). In the straight case, the azimuthal waves in the upper z half-plane have an equal and opposite contribution to those in the lower z half-plane. Meanwhile in the twisted case, the magnetic tension force associated with the azimuthal field has a positive contribution to one of the half-planes and has an inhibiting effect in the other. This causes the ω_y contribution of one set of azimuthal waves to be slightly larger than the set in the opposite half of the domain. The effect of the tension force reverses as the direction of motion changes, thus causing the observed oscillatory behaviour.

For the other two components of vorticity, the loop integrated values oscillate with twice the frequency of the kink mode. For both ω_x and ω_z, the associated panels in Fig. 17 are dominated by the contribution of loop-aligned gradients; $\frac{\mathit{\partial }{\mathit{v}}_{\mathit{z}}}{\mathit{\partial y}}$ and $\frac{\mathit{\partial }{\mathit{v}}_{\mathit{x}}}{\mathit{\partial y}}$, respectively. These gradients are associated with additional, higher harmonic azimuthal waves that are generated in the twisted field cases but are not observed in the straight field case.

Finally, the bottom row of Fig. 17 displays the maximum size of each component of the vorticity in the same horizontal plane. The central panel was used as a proxy for the growth of the KHI in Paper I. We see that despite being slightly suppressed when the KHI first forms, the maximum vorticity in the twisted case exceeds that observed in the straight case at arount t = 800 s. This supports the claim that the vortices in Case 2c have as much kinetic energy as the straight field vortices. There is a marked difference in the growth of the other two components of vorticity. In particular, the maximum size of | ω_x | and | ω_z | is around an order of magnitude larger in the twisted field case. This is caused by the azimuthal waves and the KH vortices having a significant helical component (see Fig. 11).

3.5. Loop-aligned numerical resolution

As shown above, one result of the inclusion of magnetic twist is that the KH vortices are helical in nature. Indeed, vertical structure in the density and velocity field becomes more profound in Cases 2 and 3 and, as more twist is included, this effect increases. Hence, both the horizontal and vertical resolution might affect the development of the KHI. Since, in our simulations, the loop-aligned spatial resolution is much coarser than the horizontal spatial resolution, we here investigate whether the vertical resolution is sufficient.

In all cases, the most unstable wave vector satisfies k·B = 0 such that the inhibiting effect of the right-hand side of Eq. (4.7) in Cowling (1976) is eliminated. Thus, an additional consequence of the inclusion of magnetic twist in our model is the modification of this wave vector, k. The azimuthal component of the magnetic field in Cases 2 and 3 ensures that k has a non-zero component parallel to the loop axis. In particular, it is no longer confined to the x−z-plane as it is in the straight field case.

We estimate the required vertical resolution as follows. Since k·B = 0 and B_R ≈ 0 when the instability is triggered, we have $\frac{\begin{array}{c}\mathrm{|}\\ \mathrm{|}\\ \mathrm{|}\end{array}{\mathit{k}}_{\mathit{y}}\begin{array}{c}\mathrm{|}\\ \mathrm{|}\\ \mathrm{|}\end{array}}{\left|{\mathit{k}}_{\mathit{\phi }}\right|}\mathrm{\approx }\frac{\left|{\mathit{B}}_{\mathit{\phi }}\right|}{\begin{array}{c}\mathrm{|}\\ \mathrm{|}\\ \mathrm{|}\end{array}{\mathit{B}}_{\mathit{y}}\begin{array}{c}\mathrm{|}\\ \mathrm{|}\\ \mathrm{|}\end{array}}\mathrm{·}$(16)Meanwhile, in order to resolve the vertical component of k at least as well as we resolve the horizontal component, we require${\mathrm{\Delta }}_{\mathit{y}}\mathrm{\le }\frac{{\mathit{\lambda }}_{\mathit{y}}{\mathrm{\Delta }}_{\mathit{\phi }}}{{\mathit{\lambda }}_{\mathit{\phi }}}\mathit{,}$(17)where Δ_y and Δ_φ are the vertical and azimuthal numerical resolution (grid widths) and ${\mathit{\lambda }}_{\mathit{y}}\mathrm{=}\frac{\mathrm{2}\mathit{\pi }}{{\mathit{k}}_{\mathit{y}}}$ and ${\mathit{\lambda }}_{\mathit{\phi }}\mathrm{=}\frac{\mathrm{2}\mathit{\pi }}{{\mathit{k}}_{\mathit{\phi }}}$ are the vertical and azimuthal wave lengths, respectively. We note that Δ_φ is confined by the horizontal grid widths Δ_x and Δ_z such that ${\mathrm{\Delta }}_{\mathit{x}}\mathrm{\le }{\mathrm{\Delta }}_{\mathit{\phi }}\mathrm{\le }\sqrt{{\mathrm{\Delta }}_{\mathit{x}}^{\mathrm{2}}\mathrm{+}{\mathrm{\Delta }}_{\mathit{z}}^{\mathrm{2}}}\mathit{.}$(18)Consequently, since Eq. (17) must be satisfied throughout the instability forming region, we obtain ${\mathrm{\Delta }}_{\mathit{y}}\mathrm{\le }\frac{{\mathit{B}}_{\mathit{y}}{\mathrm{\Delta }}_{\mathit{x}}}{{\mathit{B}}_{\mathit{\phi }}}\mathrm{·}$(19)In particular, we immediately notice that the greater the magnetic twist, B_φ, the more refined the vertical grid has to be in order to ensure the simulation is well resolved. We also note this is only an upper bound on Δ_y, and in the limit B_φ → 0, more restrictive bounds are required in order to resolve the initial kink mode, for example.

In our model, for the most twisted cases, Eq. (19) yields Δ_y ≲ 1200 km, suggesting that the vertical component of the wave vector k is less well resolved than the horizontal component. In order to determine whether this produced significant numerical artefacts in our simulations, we re-used the set-up of Case 2c but with double the number of grid points in the y direction.

Fig. 16 Magnitude of the vorticity at two different times, t = 420 (left) and t = 990 s (right) in the horizontal cross-section through the loop apex. The top row corresponds to Case 1, the straight field simulation, and the second row corresponds to Case 2c. Here we have normalised the vorticity to the maximum observed in either simulation.

Fig. 17 Top row: the three components of the vorticity integrated over the horizontal plane at the loop apex. Bottom row: maximum size of each component of the vorticity in this plane. For both of the rows, the black line corresponds to the straight field simulation (Case 1) and the red line corresponds to Case 2c.

Fig. 18 Top row: density profile along the most under-resolved vertical line at t = 900 s for the high resolution (blue) and low resolution (red) simulations. In the first two panels, the position of grid points are denoted with a cross, and the third panel shows a comparison between the two simulations. Bottom row: comparison of the three velocity components along the same line for the two different resolutions.

In Fig. 18, we display the density and velocity profiles along the most under-resolved vertical line in the low resolution simulation at a time when the instability is well-developed. The line lies within the vortex-forming region in the boundary of the loop. We observe that for −45 Mm ≤ y ≤ −35 Mm, the density profile is approaching the limit of the spatial resolution, particularly in the low resolution (100 grid points) simulation. Despite this, there is little difference between the density profile (third panel) and the velocity field (bottom row) produced by the two experiments. We found that there was no significant difference in the growth rate of the KH vortices or indeed in any of the plots displayed previously. However, we note that these resolution effects are expected to become more significant in regimes with greater magnetic twist.

4. Discussion and conclusions

We have conducted an investigation into the effects of magnetic twist on the dynamics of a transversely oscillating coronal loop. In all cases considered, the level of twist was weak in the sense that B_φ ≪ B_y (where B_y is the loop-aligned component of the field) and the loop was not kink unstable.

The fundamental standing kink wave was induced by an initial velocity profile and was observed to decay due to resonant absorption. Initially, this phenomenon, through which energy is transferred from the standing kink mode to localised Alfvénic waves, proceeds largely independently of the twist present within the loop. The azimuthal component of the magnetic field does modify the Alfvénic modes slightly, but the global dynamics remain unchanged when compared to the straight field case (see Paper I). In particular, the twisted field permits an additional, higher harmonic, azimuthal wave that contributes an oscillating net vorticity at the loop apex.

The azimuthal waves are contained within the boundary region of the loop and are associated with a significant radial velocity shear. In dissipative media, these waves will deposit their energy as heat through the process of phase mixing; however, due to the approximately ideal regime presented, this was not observed. Once a critical shear has been reached, the Kelvin-Helmholtz instability is triggered in all cases. The density deformation that is associated with the KHI is partially suppressed in cases with twisted magnetic field. Further, for the extent of the parameter space considered here, more twist causes greater suppression. In contrast with the dissipation dependent suppression observed in Paper I, we note that the twist does not significantly delay the onset of the KHI, but instead limits its subsequent growth. This is in agreement with the velocity shear associated with the resonant absorption growing independently of the twist.

Previous studies have demonstrated that other parameters such as the amplitude of the initial perturbation, the width of the boundary layer (Magyar & Van Doorsselaere 2016; Terradas, in prep.), and the dissipation coefficients (Howson et al. 2017) can modify the growth rate of the instability. In the first of these articles, the authors find that decreasing the amplitude of the initial perturbation decreases the instability’s growth rate. Similarly, the second work finds that in both straight and twisted field regimes, thicker boundary layers restrict the growth rate of the KHI. Both of these effects, along with the stabilisation via enhanced dissipation coefficients, are caused by a reduced velocity shear associated with the Alfvén waves. In the case of twisted field, however, the shear is only weakly affected and instead the instability growth rate is inhibited by the magnetic tension force.

Despite the suppression of the density deformation, the KHI in twisted fields may have significance for the coronal heating problem because it generates the small length scales required for an enhanced dissipation rate (Antolin et al. 2014; Howson et al. 2017). Indeed, Paper I demonstrated that even with lower dissipation coefficients, the instability may result in earlier heating than is produced by phase mixing alone at high dissipation. Whilst the suppression of the KH vortices is readily observed in both the density cross-section and the horizontal velocity field, smaller spatial scales form in the magnetic field in the cases with greater twist. Since twisted cases have an azimuthal component that is radially non-uniform, the characteristic, radial movement of plasma, and hence magnetic field (since it is approximately frozen-in to the plasma), associated with the KHI can generate large currents. Indeed, if coronal loops are twisted in some way, then some of the density deformation caused by the KHI would be suppressed, but importantly stronger currents are generated and hence greater Ohmic heating is expected. Meanwhile, analysis of the vorticity at the loop apex reveals that although the density deformations produced in the twisted cases are smaller in size, they have as much kinetic energy as in the straight case.

Additionally, once the KHI currents begin to form, the azimuthal component of the magnetic field produces a vertical component of the Lorentz force that drives upward and downward flows. These are not observed in the completely straight case but are significant even in weakly twisted regimes. This offers the possibility of detecting magnetic twist in coronal loops by searching for loop-aligned flows at the apex in decaying kink oscillations.

One further effect of the inclusion of magnetic twist is that, contrary to the straight case, the development of the KHI triggered vertical asymmetry along the loop (as a consequence of the loop-aligned flows). This was observed both in the velocity field and the cross-sectional density profile. From a regime that was vertically symmetric about the loop apex, the instability generated significant non-uniformity between the two halves of the loop. When magnetic twist is included, these effects should be considered when vertical symmetry is used in order to model only one half of the loop and hence save computational time.

A major criticism of previous work considering the Kelvin-Helmholtz instability as a driver of MHD turbulence in the corona and hence a possible coronal heating mechanism is that even weak magnetic twist will suppress the development of the instability. Whilst we have demonstrated that, to some extent, this is indeed the case, the growth of much larger currents in the twisted field cases more than compensates for the level of suppression when considering energy dissipation. The ultimate effect of the twist is an increase in the difficulty of detecting the KHI vortices observationally rather than a decrease in its efficiency as a heating mechanism. Whilst the reduced density deformation may reduce the chance of observing the instability directly, the presence of field aligned flows may offer an alternative detection mechanism.

Acknowledgments

The authors would like to thank the referee for their helpful comments and suggestions. The authors would also like to thank Prof. P. J. Cargill for helpful discussions. The research leading to these results has received funding from the UK Science and Technology Facilities Council (consolidated grant ST/N000609/1) and the European Union Horizon 2020 research and innovation programme (grant agreement No. 647214). This work used the DIRAC 1, UKMHD Consortium machine at the University of St Andrews and the Darwin Data Analytic system at the University of Cambridge, operated by the University of Cambridge High Performance Computing Service on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment was funded by a BIS National E-infrastructure capital grant (ST/K001590/1), STFC capital grants ST/H008861/1 and ST/H00887X/1, and DiRAC Operations grant ST/K00333X/1. DiRAC is part of the National E-Infrastructure.

References

All Tables

All Figures

	Fig. 1 Configuration of the modelled loop.
In the text

	Fig. 2 Initial density profile for all cases through the cross-section of the loop. The density is normalised to the initial exterior density, ρe = 8.4 × 10-13 kg m-3.
In the text

	Fig. 3 Magnitude of initial azimuthal magnetic field, | Bφ |, through the cross-section of the loop for each case. Here, Bφ is normalised to B0 = 21 G. We see that in all cases, | Bφ | ≪ | B0 |.
In the text

	Fig. 4 Initial gas pressure profile through the cross-section of the loop for each case. Here, we have normalised the pressure to the external pressure, P0.
In the text

	Fig. 5 Displacement of the loop apex through time for the straight case and a twisted case (Case 2b). The red lines are given by $\mathit{y}\mathrm{=}\mathrm{\pm }\mathrm{0.267}{\mathrm{e}}^{\mathrm{-}{\left(\frac{\mathit{t}}{\mathrm{1040}}\right)}^{\mathrm{2}}}$.
In the text

	Fig. 6 Evolution of the square of the velocity along the line y = 0 (through the loop apex), z = 0. Here v2 is normalised to the maximum value obtained during the simulation. The figure is produced from Case 2b, a twisted field simulation.
In the text

	Fig. 7 Time evolution of the azimuthal component of the velocity vφ, along a line parallel to the loop axis within the shell region. The velocity is normalised to the maximum initial velocity at that height. The upper panel shows the simulation with a straight magnetic field (Case 1) and the lower panel shows Case 2c.
In the text

	Fig. 8 Density deformation at the loop apex during the formation of the KHI. The density is normalised to the initial exterior density ρe.
In the text

	Fig. 9 Greatest distance of plasma with density, ρ = 2 from the centre of the loop apex during the formation of the KHI. In the straight field case, initially this plasma is located on the resonance layer.
In the text

	Fig. 10 Mean density along loop radii in the horizontal cross-section at the loop apex at t = 990 s. By this time, the KHI is well-developed in all cases. For comparison, the initial density profile is also included (black dashed line).
In the text

	Fig. 11 Isosurface of density, ρ = 2 for the straight field case (left-hand panel) and Case 2c, a twisted field simulation (right-hand panel). The time shown is 990 s (three and a half wave periods) after the start of the simulation when the KHI is well-developed.
In the text

	Fig. 12 Magnitude of the current in the horizontal plane at the loop apex for the straight field case. Currents less than a fifth of the maximum size of the current in this plane are coloured in blue. By this stage in the simulation (t = 990 s), the KHI is well-developed.
In the text

Fig. 13 Formation of large currents during the onset of the KHI. The top-left panel displays the maximum magnitude of the current in the horizontal cross-section at the loop apex. For each simulation, it is normalised to the maximum current in this plane at t = 560 s, immediately prior to the onset of the KHI, to show the relative growth of the currents during the development of the instability. The other three panels display the maximum magnitude of each component in the same cross-section. The units are dimensionless but are the same for each of the components.

In the text

	Fig. 14 Horizontal components of the magnetic field (vectors) and the deformed density profile (coloured contour plot) at the loop apex once the KHI is well-developed. The result shown is produced from Case 2b, a twisted field simulation.
In the text

	Fig. 15 Maximum value of |vy| at the loop apex against time for all simulations.
In the text

	Fig. 16 Magnitude of the vorticity at two different times, t = 420 (left) and t = 990 s (right) in the horizontal cross-section through the loop apex. The top row corresponds to Case 1, the straight field simulation, and the second row corresponds to Case 2c. Here we have normalised the vorticity to the maximum observed in either simulation.
In the text

	Fig. 17 Top row: the three components of the vorticity integrated over the horizontal plane at the loop apex. Bottom row: maximum size of each component of the vorticity in this plane. For both of the rows, the black line corresponds to the straight field simulation (Case 1) and the red line corresponds to Case 2c.
In the text

Fig. 18 Top row: density profile along the most under-resolved vertical line at t = 900 s for the high resolution (blue) and low resolution (red) simulations. In the first two panels, the position of grid points are denoted with a cross, and the third panel shows a comparison between the two simulations. Bottom row: comparison of the three velocity components along the same line for the two different resolutions.

In the text