# Coronal heating by the partial relaxation of twisted loops

###### Key Words.:

Instabilities, Magnetic fields, Magnetic reconnection, Magnetohydrodynamics (MHD), Plasmas, Sun: corona^{†}

^{†}offprints: M. R. Bareford

###### Abstract

Context:Relaxation theory offers a straightforward method for estimating the energy that is released when continual convective driving causes a magnetic field to become unstable. Thus, an upper limit to the heating caused by ensembles of nanoflaring coronal loops can be calculated and checked against the level of heating required to maintain observed coronal temperatures ().

Aims:We present new results obtained from nonlinear magnetohydrodynamic (MHD) simulations of idealised coronal loops. All of the initial loop configurations discussed are known to be linearly kink unstable. The purpose of this work is to determine whether or not the simulation results agree with Taylor relaxation, which will require a modified version of relaxation theory applicable to unbounded field configurations. In addition, we show for two cases how the relaxation process unfolds.

Methods:A three-dimensional (3D) MHD Lagrangian-remap code is used to simulate the evolution of a line-tied cylindrical coronal loop model. This model comprises three concentric layers surrounded by a potential envelope; hence, being twisted locally, each loop configuration is distinguished by a piecewise-constant current profile, featuring three parameters. Initially, all configurations carry zero-net-current fields and are in ideally unstable equilibrium. The simulation results are compared with the predictions of helicity-conserving relaxation theory.

Results:For all simulations, the change in helicity is no more than 2% of the initial value; also, the numerical helicities match the analytically-determined values. Magnetic energy dissipation predominantly occurs via shock heating associated with magnetic reconnection in distributed current sheets. The energy release and final field profiles produced by the numerical simulations are in agreement with the predictions given by a new model of partial relaxation theory: the relaxed field is close to a linear force free state; however, the extent of the relaxation region is limited, while the loop undergoes some radial expansion.

Conclusions:The results presented here support the use of partial relaxation theory, specifically, when calculating the heating-event distributions produced by ensembles of kink-unstable loops. The energy release increases with relaxation radius; but, once the loop has expanded by more than , further expansion yields little more energy. We conclude that the relaxation methodology may be used for coronal heating studies.

## 1 Introduction

The energy required to heat the solar corona is thought to originate from the magnetic fields that permeate the Sun’s atmosphere. The geometry of these fields is revealed by coronal loops, where the emitting plasma is constrained to follow the magnetic field: the plasma beta is extremely low (). Coronal loops are closed structures that emerge from the photosphere at one location and re-enter the solar surface at another. The convective motions at these photospheric boundaries (or footpoints) are thought to reconfigure coronal fields and thereby cause free magnetic energy to accumulate, making the fields more susceptible to instability. Essentially, the kinetic energy of the convection zone is transported, via a Poynting flux, through the photosphere and stored within the coronal loop. In the case of a single loop, currents are created by those convective motions that cause the footpoints to be twisted. The approximate time scale for the twisting motions is long compared to the Alfvén time, and so the coronal loop (not necessarily illuminated) transitions through a series of force-free equilibria, which can be expressed as ; where is a position vector, and / is related to the parallel electric current density (i.e., is a measure of the twist).

Energy release might be triggered when a loop-like magnetic field, driven by continual convective motions, reaches the threshold for kink instability (Hood 1992; Browning & Van der Linden 2003, Haynes & Arber 2007; Srivastava et al. 2010). Coronal magnetic fields cover many thousands of kilometers () and exist within a highly conductive environment: therefore, in the absence of instability, magnetic fields diffuse slowly (). Nevertheless, through an ideal kink instability, a coronal loop may be deformed such that magnetic flux surfaces are brought together. As the deformation continues, areas of high current are produced, allowing magnetic reconnection to take place. Hence, energy is released from the magnetic field during the nonlinear phase of an ideal kink instability, as shown by many three-dimensional (3D) magnetohydrodynamic (MHD) models (Baty & Heyvaerts 1996; Velli et al. 1997; Arber et al. 1999; Baty 2000; Browning et al. 2008; Hood et al. 2009).

The coronal heating problem is most pronounced within active regions, where the heating requirement is approximately , significantly higher than that necessary for the quiet Sun (Withbroe & Noyes 1977; Aschwanden & Acton 2001). These regions are often observed as dense thickets of transient coronal loops, the composition of which changes over the course of hours to days. Sudden reconfigurations are coincident with large-scale solar flares (). The signal strength of such phenomena means these events are more readily observed — detailed investigations have revealed strong evidence for magnetic reconnection (Fletcher 2009; Qiu 2009). Miniscule (nano)flares, far weaker than ones commonly observed, could maintain coronal temperatures if these less dramatic events occur with sufficient frequency (Hudson 1991). This type of flaring might be the sort initiated by the kink instability; however, the amount of energy released (i.e., the difference between the energy at instability onset and that of the relaxed field) depends on how much the magnetic field has altered before it relaxes. Bareford et al. (2010, 2011) identified the threshold for linear kink instability with respect to an idealised coronal loop model (both with and without net current). This work determined the subset of field configurations accessible via convective driving that are linearly kink unstable. One approach to calculating the energy release, is to represent each of these configurations within a nonlinear MHD code, allow the instability to take place and follow the reconnection dynamics until a relaxed state emerges. However, it is simply not feasible to run such a computationally-intensive process for more than a few examples. Hence, Bareford et al. (2010, 2011), used relaxation theory to identify the relaxed states for all of the threshold (i.e., marginally unstable) configurations.

An unstable field obeys relaxation theory if it relaxes towards the lowest energy state that conserves total magnetic axial flux and global magnetic helicity (Taylor 1974, 1986). This minimum energy state is a linear force-free field; i.e., is invariant and . The helicity () indicates how intertwined a magnetic field is with itself (Berger 1999). Coronal loops have a non-zero normal flux at the footpoints; hence, the gauge-invariance of relative helicity (Berger & Field 1984; Finn & Antonsen 1985) makes it the more useful property:

(1) |

where is the magnetic potential, is the potential field with the same boundary conditions and is the corresponding vector potential. Helicity-conserving relaxation has been seen in many laboratory experiments (Heidbrink & Dang 2000; Taylor 1986). It must be noted that helicity is subject to global resistive diffusion. Nevertheless, if localised dissipation occurs on small spatial scales (i.e., across shock fronts or within thin current sheets), the reduction in helicity will be negligible compared to the decrease in magnetic energy (Browning 1988, Browning et al. 2008). The original intention of relaxation theory was to explain laboratory plasma phenomena; but latterly, it has been frequently applied to the solar corona (Heyvaerts & Priest 1984; Browning et al. 1986; Vekstein et al. 1993; Zhang & Low 2003; Priest et al. 2005).

The relative ease with which relaxation theory can be applied meant that Bareford et al. (2010, 2011) were able to generate heating-event distributions from ensembles of idealised coronal loops, representing, albeit crudely, the population of coronal loops that exist within an active region. Each energy release is determined by where a loop crosses the instability threshold; this location is the outcome of a defined stochastic process. The distributions lead to an estimate for the heating rate that is just sufficient for coronal heating. However, the assumptions of this work regarding instability and relaxation theory have yet to be tested by a nonlinear 3D MHD code. The purpose of this paper is to elucidate further (Browning et al. 2008; Hood et al. 2009) the relaxation process and to understand how it can be applied to coronal loops that lack a conducting wall. We simulate a set of zero-net-current coronal loops that sample the linear instability threshold calculated by Bareford et al. (2011). This is a larger set than the one investigated by Hood et al. (2009), and futhermore, the loops represented here feature a current-neutralisation layer that maintains zero net current even if the currents inside the loop are predominantly single signed. (Loops that carry zero net current are preferred since the convective motions that twist the loop and thereby create azimuthal field are spatially localised; the field outside the loop is unaffected by motions within the loop cross-section and therefore remains purely axial.)

A long-standing problem has been how to apply relaxation theory in astrophysical contexts without the presence of conducting walls: simplistically, the relaxation should extend out to infinity and lead always to potential fields. Browning (1988) and Dixon et al. (1989) showed that relaxation theory could apply to volumes with free boundaries, but did not give a prediction for the spatial extent of the relaxed state. We find that a modification to Taylor relaxation theory is required before it can be used to estimate the energy released by a kink-unstable loop. In contrast to previous work, we calculate the helicity directly at specific times during each simulation (Browning et al. 2008 integrated the time differential of the helicity in order to show the change in ); thus, we are able to verify the extent of helicity conservation. We also examine the performance of the MHD code with regard to energy conservation. Any numerical dissipation will have implications for the modelling of plasma processes associated with heating, such as radiation. However, the plasma- is sufficiently low that any artificial resistivity should not influence the energy released by an unstable magnetic field, nor should it affect the evolution of a relaxing field.

The paper is structured in the following manner. Section 2 describes the numerical code used for the nonlinear MHD simulations, along with the loop model and the equations used to calculate the magnetic field. The corresponding instability threshold for the linear kink mode is introduced, as are the positions of the simulated loop configurations. Section 3 presents the results: specifically, how the different forms of energy vary over the course of the simulation, when the loop goes unstable, and how these results are affected by changes in the code parameters, such as spatial resolution and background resistivity. Following, the evolution of the loop is presented as regards magnetic field and current magnitude. Section 4 discusses how well the results fit a modified Taylor relaxation theory. Finally, in the last section, the results are summarised and our conclusions are given.

## 2 Numerical code

The nonlinear simulations are conducted using a 3D MHD Lagrangian Remap Cartesian code, called LARE3D (Arber et al. 2001). It is written in Fortran 90 and uses the Message Passing Interface (MPI) to achieve parallelisation. The Lagrangian step uses a second-order accurate predictor-corrector step that also incorporates artificial viscosity, ensuring shocks are captured accurately. Van Leer (1997) gradient limiters are used at the remap step in order to preserve monotonicity. The divergence-free condition () is maintained to machine precision by Evans & Hawley (1988) constrained transport.

LARE3D solves the resistive MHD equations given by

(5) |

with specific energy density,

(6) |

where is the mass density, is the plasma velocity, the magnetic field, the thermal pressure, is the resistivity (not magnetic diffusivity), is the current density, is the ratio of specific heats, and is the magnetic permeability. Viscous heating is represented by the last term of equation (5), which also incorporates artificial viscosity (Wilkins 1980) in order to capture the heating effect of shocks. This heating term is expressed as the product of the rate of strain tensor,

(7) |

and the shock tensor,

(8) |

where is the fast magnetoacoustic speed, is the distance across a grid cell in the direction normal to the shock front, is a similarly localised strain rate (the subscripts and denote the different spatial coordinates), and the artificial viscosity coefficients and are constants. The form of equation (8) is derived from the Rankine-Hugoniot jump conditions and the values of the coefficients have been chosen such that numerical oscillations behind shock fronts are prevented. Note, the force equation (2) also acquires a viscous term. Gravitational effects are ignored in this study, as are thermal conduction and radiation. The simulations are concerned with how the magnetic field changes in response to the kink instability; specifically, how much magnetic energy is released and how the field subsequently evolves. Conduction becomes important some time after the energy release and later, radiation is the dominant process. Note, numerical studies have shown that conduction can act on MHD time scales (Botha et al. 2011): the amount of energy released from the field is unaffected, but the kinetic energy parallel to the field is much reduced.

The MHD equations are made dimensionless by replacing the variables with dimensionless equivalents. For example,

where asterisks denote dimensional variables, is the loop radius, the initial mass density, and the initial axial field at . The other variables are expressed in a similar manner;

where is the loop length, is the radial Alfvén transit time, the Alfvén speed, and the magnetic pressure. The specific energy density, current density and resistivity (, and ) also have reference variables that can be expressed in terms of , and :

Values appropriate for a coronal active region can be obtained by setting , and . Hence, the length becomes , and .

The resistivity is taken to be non-uniform in these simulations,

where is the background resistivity (normally set to zero, since actual coronal resistivities are approximately ) and is the anomalous resistivity, which is only switched on when the current reaches or exceeds . The value of is set so that it is significantly higher than the maximum current at the start of the simulation. Super-critical currents appear when, during the nonlinear evolution of the kink instability, current sheets begin to form and decrease in thickness. The anomalous resistivity is intended to capture the dissipation occurring at scales below the grid resolution: at this scale, resistivity is enhanced by small-scale plasma instabilities.

The computational domain is a 3D staggered grid: physical variables are not calculated at the same place for each cell in the domain. This approach improves numerical stability and allows conservation laws to be included in the computation. The domain size is (-3:+3) and (-10:+10). Initially, the loop axis follows the -axis and the loop radius is ; therefore, the simulated loops all have an aspect ratio of 20. The loop is line-tied at , which means, at those -coordinates, the velocity components are set to zero. The velocity components are zero at the boundaries for all directions. The normal derivatives of magnetic field, energy and density are zero at all boundaries. The simulations are run with two grid resolutions: (low) and (high). It is assumed that a result is not a numerical artefact if it is consistent across both resolutions.

### 2.1 Initial configuration

Some previous studies have used LARE3D to simulate the application of kink perturbations to a straightened line-tied coronal loop, see Gerrard et al. (2002), Gerrard & Hood (2003), Browning et al. (2008) and Hood et al. (2009).

The initial equilibrium model used in the latter two studies was extended by Bareford et al. (2011) to include an outer current-neutralising layer so as to ensure the loop has (at least initially) zero net current: this improves on the model used by Browning & Van der Linden (2003) and Bareford et al. (2010), which allowed loops to have net current (i.e., an azimuthal field was usually present in the potential envelope). All currents are now created by convective motions local to the loop footpoints. Hence, a current neutralisation layer is introduced here, defined such that the azimuthal field () always falls to zero at the loop boundary (); therefore, is zero in the potential envelope. The loop’s radial -profile is approximated by a piecewise-constant function featuring three parameters (Figure 1): the ratio of current to magnetic field is in the core, in the outer layer, in the neutralisation layer and zero in the potential envelope. The free parameters are and , whereas is dependent on the first two and is determined by the requirement of zero net current. The magnetic field is continuous everywhere, whereas the current has discontinuities, and the outer surface of the potential envelope, representing the background corona, is placed at (three times the loop radius); this is far enough away that the boundary conditions do not influence the plasma evolution.

The fields are expressed in terms of the well-known Bessel function model, generalised to the concentric layer geometry (Melrose et al. 1994; Browning & Van der Linden 2003; Browning et al. 2008). The field equations for the four regions (core, outer layer, neutralisation layer and envelope) are as follows:

(9) | |||||

(11) | |||||

(13) | |||||

(15) | |||||

(16) |

where () represents the sign of . The fields must be continuous at the inner radial boundaries, , and . (The positions are , and , so that most of the loop is similar to the one described by Bareford et al. (2010), but with a thin current neutralisation layer between and .) Therefore, the coefficients and () are determined by the requirement of continuity of the magnetic field at all interfaces (Bareford et al. 2011). The value of (the neutralisation layer current) is found, for a given (), by numerical solution of , ensuring that the net current is zero and that the azimuthal field vanishes outside the loop, see Eq. (2.1). From the nondimensionalisation of the magnetic field, the field coefficient at the core is .

The linear kink instability threshold for this current-neutralised loop was determined by the CILTS code (Van der Linden 1991; Browning & Van der Linden 2003) — it uses a bicubic Hermite finite element method to calculate the growth rates and eigenfunctions for specific line-tied -configurations.

In contrast to the stability space for a loop of net current (Bareford et al. 2010), the instability threshold is open, see Figure 2. This is very similar to the threshold found for a close-fitting conducting shell (Browning & Van der Linden 2003), since in the case of zero net current, perturbations quickly fall to zero beyond the loop radius. A consequence of the Bessel function model is that the axial field changes sign if the values of and are too large. This reversal is unphysical and hence, introduces a restriction to the stable region of Figure 2; namely, all stable configurations should have of uniform sign. The new stability space is closed by excluding those field profiles that have a of mixed polarity, since this could not be achieved directly by footpoint motions of an initially unidirectional field. The filled circles of Figure 2 identify the loop configurations (see also Table 1) that will be simulated by the LARE3D code (the initial field profiles for Loops B and D are given in Figure 3). All of these configurations are unstable to the ideal kink instability.

Loops of uniform twist (loops A–C) are a more likely result of correlated convective driving (Bareford et al. 2011), where the threshold is approached via a series of steps with — this is the reason why the part of the threshold curve where is more finely sampled compared to . Note, the section of threshold on the left of Figure 2 is merely the negative of the section on the right; hence, it is not necessary to sample both threshold sections.

Loop | |||||||||
---|---|---|---|---|---|---|---|---|---|

A | 2.42 | 2.4 | -13.08 | 15.3 | 0.1 | -0.0079 | 2.42 | -0.59 | 0.55 |

B () | 2.25 | 1.5 | -8.71 | 18.1 | 0.79 | -0.2 | -0.27 | 2.37 | 0.64 |

C | 2.15 | 0.53 | -4.95 | 20.8 | 0.61 | -0.15 | -0.94 | -1.85 | 0.74 |

D () | 2.54 | -1.0 | -0.84 | 21.0 | 0.91 | 0.50 | 0.92 | 0.38 | 0.74 |

E | 2.8 | -2.7 | 3.82 | 20.4 | 0.26 | 1.32 | -1.79 | 0.026 | 0.72 |

Stable | 0.5 | 0.1 | -0.96 | 27.8 | 0.97 | -0.0078 | 1.21 | 0.63 | 0.98 |

The configurations listed in Table 1 encompass two types of loop; one class where and have the same sign (A–C) and the other where these parameters have opposing signs (D and E). Loop B has been chosen as representative of the first loop type and Loop D of the second. Henceforth, Loop B will be referred to as (uniform twist) and Loop D will be denoted by (mixed twist).

Each of the loops indicated in Figure 2 is subjected to a disturbance of the form,

(17) |

where is the radial component of the perturbed velocity, is the loop length, is the radial coordinate, is the wave number, , and the constant reduces the amplitude so that the perturbation is initially linear. This disturbance initiates the kink instability.

Both and are simulated for low and high resolutions; only the high resolution is used for the other loop configurations. Each simulation runs until the magnetic field appears to have settled into a lower energy state.

## 3 Numerical results

### 3.1 Energy and resistivity

Loop (, ) is linearly kink unstable, and the numerical simulation (Figure 4, middle row) shows that it is also nonlinearly unstable. The magnetic energy, , has an initial (dimensionless) value of 155.3 when integrated over the entire simulation domain,

the internal () and kinetic () energies are also volume integrals,

The initial magnetic energy integrated over the loop only is much less (). All magnetic energies can be dimensionalised by making a simple correction to equation (27) given in Bareford et al. (2010): since here the axial flux is not normalised to 1, the dimensionalising multiplier must first be divided by , which is the square of the non-dimensional axial flux over the radius (Table 1, column 5). For loop , this gives a multiplier of (assuming a radius of and a mean axial field strength of ): thus, the dimensionalised initial loop energy is roughly .

At the onset of instability, undergoes a decrease, coincident with a rise in and with a much more modest increase in kinetic energy (). The maximum current, , also rises just before the decrease in , indicating the formation of a helical current sheet. The nonlinear instability starts at around and within the next , approximately 70% of the total energy release has been achieved. Magnetic energy reduces more slowly after . The release of magnetic energy is of a similar size for both resolutions, and significantly, larger currents are recorded at high resolution, which is expected; otherwise, spatially-confined changes in current are missed and anomalous resistivity is reduced. The fact that the maximum current increases with resolution is indicative of current sheet formation. These structures have (possibly) infinite current density, so higher resolutions should reveal larger values of the maximum.

In the top two cases of Figure 4, ideal MHD and zero background resisitivity, it can be seen that there is no detectable change in energy until around — this is the time when a current sheet starts to form. The adjective ideal is italicised because the rate of magnetic diffusion in the corona is easily exceeded by numerical resistivity; therefore, any reduction in magnetic energy associated with a stable configuration is artificial. Hence, for the stable case, in which no current sheets ever form, there is no sudden onset of dissipation, only gradual dissipation as expected. In the bottom row of Figure 4 (), there is a continual dissipation of magnetic energy due to Ohmic effects; however, even here, there is a clear onset of enhanced dissipation at the point of current sheet formation.

This all agrees with previous work (Browning et al. 2008; Hood et al. 2009) and clearly indicates that energy dissipation (conversion from magnetic energy to internal energy) is closely associated with the existence of current sheets.

However, Figure 4 also shows that energy is not conserved for the two cases mentioned (ideal MHD and ). At high resolution, only of the energy released from the magnetic field is converted to internal energy, and by the end of the simulation, less than one percent of the energy release is in the form of kinetic energy. Halving the spatial resolution worsens the energy conservation by almost . This loss of energy is due to spatially unresolved current sheets, resulting in numerical diffusion. The results for Ideal MHD (i.e., zero magnetic diffusivity) are very similar to the ones produced for the resistive MHD case with no background resistivity. This shows that, for these cases, Ohmic dissipation does not contribute significantly to the rise in internal energy, which is instead mainly caused by viscous heating. The use of a non-zero background resistivity () changes the plots. The reduction in magnetic energy starts right away and is more drawn out and the maximum current is much less than when . These differences are all consistent with an increased resistivity. Crucially, Ohmic heating is now clearly present and, with the viscous heating, accounts for nearly all of the dissipated magnetic energy. In terms of energy conservation, LARE3D performs better when : the internal energy increase is now approximately of the magnetic energy decrease (the final kinetic energy is less than of ). A background resistivity of appears to be the smallest value that effectively minimises the effect of the numerical resistivity, such that the results are likely to be a reasonable description of the thermodynamics. If is halved, the percentage of the magnetic energy release lost to the simulation increases to around . Further tests have shown that energy conservation is not improved by a lowering the critical current and thereby causing the anomalous resistivity () to be applied earlier in the simulation.

A sufficiently large background resistivity does substantially mitigate the amount of energy lost by numerical resistivity, at least for ; but is not realistic by coronal standards, and Figure 4 (bottom left) reveals a significant drop in field energy before the kink instability takes effect. This initial decline is caused by global Ohmic diffusion rather than magnetic reconnection, since current sheets have not yet formed. The numerical resistivity comes about when current sheet widths begin to fall below the grid resolution; although magnetic energy continues to be dissipated, it is not converted to other forms of energy (e.g., internal or kinetic). However, shocks can still be resolved during the unstable phase, and these contribute to the internal energy via viscous heating. If the background resistivity is large enough it will limit currents and thereby prevent current sheets from becoming too thin. Numerical (that is to say artificial) resistivity will be a factor during the conversion process when : approximately of the reduction in magnetic field energy is not accounted for by the final internal energy — this is also true for the other loop configurations.

However, virtually all of this artificial resistivity is occurring during the energy conversion. The drop in magnetic energy is a robust result. To demonstrate further, we have used ideal MHD to simulate a loop ( and , see Table 1, bottom row) that is well within the stable region (as shown in Figure 2). In the absence of an instability (and any applied resistivity), the energy declines by around one thousandth of one percent over ; this reduction is three orders of magnitude smaller than the energy release caused by the kink instability.

Loop (, ) has a core that is oppositely twisted with respect to its outer layer. Figure 5 shows the same correspondence between the magnetic and internal energies that was seen for the previous loop. Again, the magnetic energy release is of the same size for both resolutions and, at the higher resolution, significantly larger currents are recorded; although the size of the release is around half that found for . Note, the initial magnetic energy is higher than the value given for the other loop. This is a consequence of setting : is less twisted and therefore has a higher axial flux, which results in a lower dimensionalised magnetic energy.

The general trends for magnetic, internal and kinetic energies are consistent between resolutions (for both loops), and most importantly, so is the size of the magnetic energy release. Therefore, simulations at higher resolution are not required — this paper will proceed with results taken only from the simulations. The results for the other loop configurations on the threshold (Figure 2) also suggest that linear instability evolves to a nonlinear stage, which gives rise to current sheets, magnetic reconnection and most importantly, shock heating. The following sections will present results based on a current-dependent resistivity and zero background resistivity — this is more compatible with the coronal environment. The breakdown of energy conservation associated with this parameter choice can be ignored since we are only interested in how the magnetic field behaves during and after the instability. However, this issue will have to be addressed for detailed studies of the thermodynamic evolution.

### 3.2 Magnetic field and current magnitude

Now we examine the magnetic field (and critical current distribution) at specific times during the simulations. Figure 6 shows how field lines, originally located within the core, become kinked as the instability takes hold. The dark grey field lines are drawn from the bottom boundary (or footpoint) and the light grey field lines are drawn from similar locations at the upper boundary. Loops and follow the same course of events. Initially, the field lines are intertwined; then, during the growth of the instability, the currents become critical (indicated by the yellow, orange and red areas) and anomalous resitivity is applied. The dissipative effects of this anomalous resistivity have a minimal contribution to the internal energy (Section 3.1); instead magnetic energy is dissipated by the application of viscous heating at shock fronts. Hence, we also show a proxy for viscous heating, , which is equivalent to the shock tensor divided by — is represented by the cyan, blue and purple colours. In general, shocks are coincident with current sheets: i.e., the areas of high viscous heating, as indicated by , are cospatial with the largest currents. Note, the kink instability is more pronounced for than for . In terms of azimuthal field, is the weaker of the two (Figure 3), however, some importance should be attached to the fact that is twisted both ways. Opposing twists appear to mitigate the growth of the instability and thereby limit the energy release. If we increase the values of and but keep the opposite signs, such that the total azimuthal field strength is comparable to (e.g., Loop E), the energy release increases only slightly ().

Figure 7 shows cross sections of at (halfway along the loop) and at , which is roughly the centre of the only patch of significant viscous heating for at , see bottom row of Figure 6. Again, the colours represent current magnitudes (left column) and viscous heating (right column), and the plot times are the same as those used for Figure 6; i.e., shortly after the start of the kink instability. At this time, current sheets of narrow width start to form; furthermore, nearly all of these current sheets are associated with shock formation and viscous heating.

The unstable phase is over quickly () and by the end of the simulation the (reconnected) field lines have straightened considerably, indicative of a low constant- configuration. The areas where shocks have formed must be dispersed throughout the loop volume in order for helicity to be more evenly redistributed and thereby create a linear -profile. The reduction in the azimuthal components of the field lines should cause a radial expansion of the loop. At the initial equilibrium, the inward tension force of the azimuthal field is balanced by the outward magnetic pressure of the axial field; thus, if the tension decreases, the loop must expand before equilibrium can be regained. This behaviour is clearly demonstrated by Figure 8; note the change in scale for the and axes. The expansion of the loop is mostly associated with the reconnection of initially twisted field lines inside the loop with the ambient axial field. In Figure 8, the final state calculated by the numerical simulation is overlaid with magnetic field vectors in the - plane, which are consistent with a cylindrical configuration, bounded by a current-neutralising layer: the arrows follow each other and the arrow sizes initially increase away from the axis, and then diminish before the loop edge.

Loop expands less than , which is possibly due to the fact that the latter releases more energy. A single weakly-twisted flux tube best describes the final state for both loops. The following sections will show that the properties of these loops are consistent with relaxation theory.

### 3.3 Helicity conservation

DeVore (2000) showed how to calculate the magnetic helicity over an entire coronal volume above a photospheric bounding surface. The first step is to work out the magnetic vector potential for a current-free field that has the same distribution of vertical magnetic flux at the lower boundary. DeVore begins by deriving an expression for the scalar potential, using Green’s function for Laplace’s equation as the integration kernel,

(18) |

where . The grid domain used by LARE3D has a Cartesian geometry: the coronal loop is initially represented as a straight cylinder within a rectangular box. The and axes extend between and ; hence, the integral limits given above. The photospheric boundaries are located at the limits of the axis () and is the loop apex; Eq. (18) uses the first boundary position. The vector potential is constructed using,

(19) |

which becomes,

(20) | |||||

when the derivatives are moved inside the integral. Now the gauge-invariant vector potential can be specified as,

(21) |

by subtracting the helicity due to the potential field, and Eq. (21) can be re-expressed by expanding the cross product of the second term,

(22) | |||||

Finally, the gauge-invariant magnetic helicity is

(23) | |||||

The geometry used by DeVore differs significantly from that used here (Figure 1), which features two separate photospheric boundaries at the limits of the axis. Fortunately, the relative positions of the two boundaries mean that if the flux is cancelled for one it will be cancelled for the other, and so the lower bound coordinate can simply be set to .

Equations (18)–(23) have been implemented, using the five-point Newton-Cotes integration formula. The marginally unstable loops (Figure 2) have zero net current initially and should continue to do so during the simulation; thus, outside the loop the helicity is zero. This means the helicity, calculated using a straightforward cylindrical geometry, can be compared easily to that calculated for a Cartesian geometry, where the loop is enclosed within a rectangular box. The helicity is zero everywhere in the additional volume between the surface of the rectangular box and the outer edge of a cylindrical potential envelope.

Loop | Initial | Instability | Final | |
---|---|---|---|---|

A | 12.3 | 12.26 () | 12.22 | 0.09 |

B () | 12.29 | 12.27 () | 12.28 | 0.03 |

C | 10.47 | 10.46 () | 10.5 | 0.19 |

D () | 6.12 | 6.12 () | 6.11 | 0.13 |

E | 1.16 | 1.14 () | 1.18 | 1.32 |

Table 2 gives the helicities for each loop at three different times. The second time is the time of instability; i.e., when the loop is furthest from equilibrium. The helicity appears to be conserved: it varies little over the course of the simulation. For loops A–D, the helicity varies by less than one percent; whereas for Loop E the helicity increases by .

Next, we calculate the ratio of the helicity variation to the change in magnetic energy,

(24) |

both and are weighted by initial value. For four of the five loops, is around one order of magnitude lower than (Table 2, fifth column); these results are comparable with Browning et al. (2008). The relationship implies that magnetic energy dissipation is taking place on small spatial scales (i.e., shock fronts). One of the loops (E) has ; however, it is no coincidence that this loop also has the lowest helicity (). The coarseness of the grid sampling used to calculate Eq. (23) — the and dimensions are sampled at every fourth cell and the dimension at every other cell — is too great for such a small helicity. Hence, for this last loop, the grid sampling is increased from to , which gives . A finer sampling and/or higher resolution is required to bring this ratio down to below 0.1.

## 4 Partial relaxation model

### 4.1 Analytical calculation

In Browning & Van der Linden (2003), Browning et al. (2008) and Bareford et al. (2010), a loop with net current was assumed to relax such that it expanded to fill the entire potential envelope: the -profile was invariant between and (). The relaxed alpha was identified by assuming that (axial flux) and (the normalised helicity) were conserved over the loop and envelope, in accordance with Taylor relaxation. The only limit to the relaxation is the position of an (unphysical) bounding wall. Hence, the relaxed state always represented a threefold radial expansion of the initial state. Later, Bareford et al. (2011), considered a range of relaxation radii: . Some form of partial relaxation is more likely to be relevant to the zero-net-current case; it is known from simulation results, presented here and in Hood et al. (2009), that reconnection is of limited extent, leaving much of the external field undisturbed. This contrasts with previous results for loops carrying net current (Browning et al. 2008), in which the disturbances in the nonlinear phase of the kink generally extended to the boundary of the simulation.

Bareford et al. (2011) maintained zero net current by fixing a neutralising loop surface at a specified radius (). The neutralising surface is imposed by fixing the field coefficients of the potential envelope, so that they do not change during relaxation. In the relaxed state, the envelope is the region between and . Naively, it might have been expected that the partially-relaxed state would consist of a constant field embedded directly in a potential field () with no layer of reversed current. However, such an arrangement does not match the simulation results, see Section 4.1.1. Furthermore, it would be unphysical for the partially-relaxed state to include a region of azimuthal field extending to infinity, which is the consequence of allowing net current, since in the initial state, the azimthal field is zero outside the loop. The existence of current sheets (here broadened into a narrow current-neutralising layer) bounding localised relaxed states has also been noted by Gimblett et al. (2006). The axial flux is conserved such that of the threshold state is equal to of the relaxed state. The subscript denotes the relaxation radius, ; it is the radial upper bound over which the associated property is calculated — the lower bound being the axis. For example, is the axial flux from to (similarly, is the helicity from to , the outer edge of the potential envelope). This method of conservation is illustrated by Figure 9. Left, is the marginally unstable threshold state; the radius that this loop will expand to is indicated by the dashed circle; right, is the expanded relaxed loop. We might expect that, either the relaxed loop flux () matches the initial loop flux (), or it matches the initial flux within the radius () eventually attained. The former is correct if the field freely expands into the surrounding field. The latter choice is made if the loop radially expands by reconnection; i.e., the loop eats into the surrounding axial field. This outcome agrees much better with the simulation results (again see Section 4.1.1). If the loop did not reconnect with its surroundings and expand to a radius of (as for ), then the axial field would drop by a factor of , which is not observed. (The slight mismatch between the numerical and analytical axial field profiles seen in Figure 11 is probably due to some limited free expansion of the loop.)

The relaxation alpha () is determined by ensuring that the axial flux and helicity integrated over the dashed circle in Figure 9 (left), match the values obtained when these same properties are integrated over of the relaxed loop (Figure 9, right). Also, helicity is absent from the potential envelope surrounding a zero-net-current loop; hence, (unlike magnetic energy, ). Thus, we do not need to choose which value of helicity to use. Once is known, the energy release can be determined analytically;

(25) |

where is the energy of the threshold state and is the energy of the relaxed state.

The relaxation radius () remains as a free parameter, although sensible values can be inferred from the numerical results. Figure 10 shows how and vary with relaxation radius. As expected, these figures show an inverse relationship between and . In general, increases with , however, this relationship is not linear; beyond a moderate expansion of 50% () the energy release is 99% of its maximum for and 80% for . This property of diminishing returns is also true for the other three loops indicated in Figure 2. This means that the energy release is insensitive to the choice of relaxation radius, so long as it is assumed that .

#### 4.1.1 Numerical Relaxed States

The final numerically-determined state is merely the last snapshot provided by the simulation; it is expected to be close to the analytically-determined relaxed state.

Each loop is simulated for at least ; so, the sooner the instability occurs, the closer the loop will be to a fully relaxed state by the end of the simulation. The Cartesian components of the analytical relaxed field are as follows,

(26) | |||||

(27) | |||||

(28) |

where and . We generate relaxed configurations for every value of between 0.9 and 3.0 in increments of 0.01. Then, for each field component, it is checked which of the 211 possibilities has the lowest chi-squared value when compared with the numerical values for the same component. These field comparisons (, and ) are performed over the dimension for a selection of fifteen coordinate pairs (, ) — the horizontal dashed lines of Figure 13 (left) indicate the positions of the field plots.

As a result of the kink instability, the loop axis will be shifted from the origin in the plane. The new axis position will be -dependent and can be recalculated by assuming that the relaxed loop is current neutralised (Figure 8): moving inwards from the edge of the envelope, the loop boundary is detected whenever the -value rises above some minimal value. Any axis shift is taken into account before the analytical results are compared with the numerical data. Figure 11 presents a subset of the analytical-numerical comparisons for at .

In general, the analytical plots, such as the ones in Figure 11, show a good agreement with the numerics; however, the and that best fit the red lines are independently chosen for each of the forty-five subplots (there are fifteen positions and three field components). Figure 12 (left) shows the variation in these best-fit parameters: each symbol, by virtue of its position, gives the pair that best fits the numerical data extracted at specific and coordinates within the simulation domain — the symbol type denotes which field component is being matched. The use of two -axes (one for , the other ) means that all forty-five symbols can be clearly distinguished (the -axis has no meaning beyond a simple sorting of the data points).

The derived averages are and ; these two properties have a one-to-one mapping. The dimensionless energy release is . Despite the large scatter for , the deviation for is comparatively modest, this is because is small when (Figure 10, right). Figure 12 suggests that the final numerically-calculated state is tending towards a relaxed loop that can be characterised by a localised invariant -profile. (Note, all points would lie on a single vertical line if there were a perfect match between the analytical and numerical models.) The outlying points on the left have values that are further from than those on the right, since the relationship between and is not linear (Figure 10, left) - this explains the large deviation for the relaxation alpha. Low levels of magnetic field that fluctuate around the zero line tend to be fitted by low values; whereas values significantly greater than the mean imply that high currents still exist within the final numerical state. Both of these effects can be ameliorated by re-running the simulation with background resistivity; once again, . Figure 12 (right) yields and . The mean dimensionless energy release becomes . The resistivity smooths out low-level noise and restricts the current, and thereby reduces the deviation associated with the analytical fit.

Analytical | Numerical | |||||
---|---|---|---|---|---|---|

Simulation | ||||||

A | 1.9 | 0.25 | 12.89 | 12.19 | 8.31 | |

B () | 1.76 | 0.28 | 12.92 | 12.3 | 5.8 | |

C | 1.8 | 0.2 | 11.03 | 10.52 | 3.15 | |

D () | 1.59 | 0.19 | 6.45 | 6.14 | 2.15 | |

E | 1.18 | 0.19 | 1.22 | 1.15 | 2.58 |

The overall impressive level of agreement demonstrated for , also extends to other positions along the instability threshold, see Table 3. This table also confirms that the analytically-calculated helicities of Loops A–E are in approximate agreement with the numerical values, which were derived according to the procedure discussed in Section 3.3. The correspondence between the numerical and analytical energy releases is the most significant finding. There is evidence to suggest that this correlation persists even when different settings are used for the LARE3D parameters controlling resistive MHD (Section 3.1 and Figure 4). In addition, these results are consistent with previous work (Browning et al. 2008; Hood et al. 2009).

### 4.2 Final Current Distribution

The jaggedness of the numerical field profiles shown in Figure 11 clearly indicates that there are deviations from a simple locally-constant profile.

We show the value of computed over the midplane () of , see Figure 13. A grey circle based on the current magnitude plot of Figure 8 (left) has been used to approximate the shape of the cross section. There are many, albeit confined, areas that are far from the calculated mean, . The influence of the initial (the -value for the current neutralisation layer) can be seen in the limits of the colour bar. In addition, plasma- (Figure 13, right) has, compared to initial values of , undergone localised increases of up to two orders of magnitude.

The cosine of the angle between the current and the magnetic field is plotted in Figure 14; note, positions outside the loop are left unplotted. The current is either parallel (white) or anti-parallel (black) to the field for a high percentage of the cross sectional area (positions that are far from parallel, , account for approximately 30% of the loop cross section). Almost three quarters of the energy released from the field becomes internal energy (i.e., thermal pressure), which may be associated with departures from a force-free state (although Taylor relaxation would predict a uniform pressure distribution). However, it seems that the loop is heading towards an approximate balance of forces (): Figure 14 (right) shows that the angle between the field and the pressure gradient is on average close to zero. Interestingly, the parallel and anti-parallel areas are often found next to each other and therefore might be expected to cancel at later simulation times. The plots of Figures 13 and 14 are qualitatively similar to those taken at other coordinates that are not too close to the footpoints (e.g., ).

## 5 Summary and conclusions

A nonlinear 3D MHD code has been used to simulate the evolution of a set of zero-net-current cylindrical loops. These loop configurations have been identified by a linear analysis as being marginally kink unstable (Bareford et al. 2011). The simulations show that the instability quickly enters a nonlinear phase and magnetic energy declines sharply before leveling off. Furthermore, the amount of energy released matches the amount predicted by Taylor relaxation (Table 3), taking account of the fact that the relaxation is localised. Evidence for helicity conservation was presented, and the change in helicity was shown to be much smaller than the drop in magnetic energy. The implication of this result is that energy diffusion is occuring on small scales compared to the global length scale; i.e., within shocks associated with magnetic reconnection. The low values of kinetic energy during the unstable phase (compared to the internal energy) imply that these shocks occur near to reconnection sites. (We find that dissipation within current sheets only becomes significant when .) Further research could reveal exactly how magnetic reconnection spreads through the loop volume in response to a kink instability, which would also reveal where the loop is heated and when. It is this widespread dispersal of reconnection and shock heating that ensures helicity conservation. At present, we cannot confirm the nature of the shocks: slow-mode shocks are expected since only shocks of this type can reduce magnetic field strength.

Relaxation theory also predicts that the final relaxed state should have a constant -profile. Although the final numerical -profile still retains much fine structure, the final magnetic fields are well-modelled by a (localised) constant- profile with some fluctuations superposed — this suggests the fine-scale structure is self-cancelling (i.e., it integrates out). Furthermore, the energy of the final numerical state is very well matched by the energy of the same constant- state (the field energy is insensitive to the spikeness of the numerical data). The property of zero net current is retained after the instability. Typically, the loop expands radially, the field reconnecting with that present in the potential envelope. These results justify the choices made by Bareford et al. (2011) regarding the details of the relaxation process. Nevertheless, we were concerned that the thinness of the current neutralisation layer (Figure 1) might have influenced the results. Hence, we also ran a resistive MHD simulation (with ) for a zero-net-current equilibrium that possessed a smooth -profile, see Case 3 of Hood et al. (2009). The constant parameter, , was set such that the equilibrium was on the threshold of instability. We found that the simulation results were similar to those produced by . Numerically, the energy release was 1.2, which was again consistent with the analytically-determined value.

It appears that the assumption of a Taylor-relaxed state, subsequent to a kink instability, has been verified by the work presented here. The relaxation does not extend over the full numerical volume, but over a region of smaller extent (out to a radius , which is less than the full radius, ). In this sense, the relaxation is partial. A relaxed state can only be identified if the relaxation radius is known; at present, it is unclear how can be precisely determined from the field configuration at instability onset. However, the analytical work has revealed that for marginally-unstable loops, the energy release varies little with relaxation radius once (Figure 10); hence, a calculation of the energy release does not necessarily require a precise prediction for .

Energy release could be limited if the unstable loop attains an equilibrium that is less than fully relaxed (i.e., the -profile remains nonlinear) and still conserves helicity. There is perhaps, for some field configurations, another constraint that decides the relaxed state, such as the topological degree of the field line mapping between the ends of the loop, as investigated by Yeates et al. (2010). They examined two braided magnetic field configurations (one based on the simple pigtail braid and the other more complex). Both configurations underwent turbulent relaxation, leading to a final state that conserved topological degree and was less relaxed than that predicted by Taylor theory — the final state for the pigtail braid featured two flux tubes of opposite twist. Nevertheless, it is possible for the Taylor-relaxed state and the state that preserves topological degree to coincide. In our case the invariants given by Yeates et al. do not provide any extra constraint, making our results consistent with their predictions. This would explain the level of agreement between the LARE3D simulations and Taylor relaxation.

An issue for further research concerns the interaction of convective driving with the relaxation process. The LARE3D code could be used to help resolve this issue. It should be possible to choose a loop configuration (i.e., a set of -parameters) that is just inside the threshold for linear kink instability and then, trigger the instability by applying a pre-determined velocity profile () at one of the footpoints. Loop curvature has also not been considered. The linear stability codes require an analytical form for the magnetic fields. If a loop is to retain its curvature, it can only be simulated numerically, which means choices have to be made concerning loop parameters (e.g., length, radius and -profile). Usefully, those straightened loop configurations that are kink unstable and are likely to be reached by convective driving have been identified. These configurations could be adapted to include curvature and re-simulated within LARE3D. This would reveal what effect, if any, curvature has on the energy release precipitated by kink instability. Of course, this procedure could also be applied to other improvements; e.g., gravity (with ), conduction and radiation. However, a feature that improves the realism of the loop model may not be important as regards kink instability and Taylor relaxation.

Acknowledgements. We thank the referee for helpful comments and M. R. B. acknowledges financial support from STFC.

## References

- (1) Arber, T. D., Longbottom, A. W., & Van der Linden, R. A. M. 1999, ApJ, 517, 990
- (2) Arber, T. D., Longbottom, A. W., Gerrard, C. L., & Milne, A. M. 2001, J. Comput. Phys., 171, 151
- (3) Aschwanden, M. J., & Acton, L. W. 2001, ApJ, 550, 475
- (4) Baty, H., & Heyvaerts, J. 1996, A&A, 308, 935
- (5) Baty, H. 2000, A&A, 360, 345
- (6) Bareford, M. R., Browning, P. K., & Van der Linden, R. A. M. 2010, A&A, 521, A70
- (7) Bareford, M. R., Browning, P. K., & Van der Linden, R. A. M. 2011, Sol. Phys., 273, 93
- (8) Berger, M. A., & Field, G. B. 1984, J. Fluid. Mech., 147, 133
- (9) Berger, M. A. 1999, Plasma Phys. Control. Fusion, 41, B167
- (10) Botha, G. J. J., Arber, T. D., & Hood, A. W. 2011, A&A, 525, A96+
- (11) Browning, P. K., Sakurai T., & Priest, E. R. 1986, A&A, 158, 217
- (12) Browning, P. K. 1988, J. Plasma Phys., 40, 263
- (13) Browning, P. K., & Van der Linden, R. A. M. 2003, A&A, 400, 355
- (14) Browning, P. K., Gerrard, C., Hood, A. W., Kevis, R., & Van der Linden, R. A. M. 2008, A&A, 485, 837
- (15) DeVore, C. R. 2000, ApJ, 539, 949
- (16) Dixon, A. M., Berger, M. A., Browning, P. K., & Priest, E. R. 1989, A&A, 225, 156
- (17) Evans, C. R., & Hawley, J. F. 1988, ApJ, 332, 659
- (18) Finn, J. M., & Antonsen, T. M. 1985, Commun. Plasma Phys. Controlled Fusion, 9, 111
- (19) Fletcher, L. 2009, A&A, 493, 241
- (20) Gerrard, C. L., Arber, T. D., & Hood, A. W. 2002, A&A, 387, 687
- (21) Gerrard, C. L., & Hood, A. W. 2003, Sol. Phys., 214, 151
- (22) Gimlett, C. G., Hastie, R. J., & Helander, P. 2006, Plasma Phys. Control. Fusion, 48, 1531
- (23) Haynes, M., & Arber, T. D. 2007, A&A, 467, 327
- (24) Heidbrink, W. W., & Dang, T. H. 2000, Plasma Phys. Control. Fusion, 42, L31
- (25) Heyvaerts, J., & Priest, E. R. 1984, A&A, 137, 63
- (26) Hood, A. W. 1992, Plasma Physics and Controlled Fusion, 34, 411
- (27) Hood, A. W., Browning, P. K., & Van der Linden, R. A. M. 2009, A&A, 506, 913
- (28) Hudson, H. S. 1991, Sol. Phys. 133, 357
- (29) Melrose, D. B., Nicholls, J., & Broderick, N. G. 1994, J. Plasma Phys., 51, 163
- (30) Priest, E. W., Longcope, D. W., & Heyvaerts, J. 2005, ApJ, 624, 1057
- (31) Qiu, J. 2009, ApJ, 692, 1110
- (32) Srivastava, A. K., Zaqarashvili, T. V., Kumar, P., & Khodachenko, M. L. 2010, ApJ, 715, 292
- (33) Taylor, J. B. 1974, Phys. Rev. Lett., 33, 1139
- (34) Taylor, J. B. 1986, Rev. Mod. Phys., 58, 741
- (35) Van der Linden, R. A. M. 1991, Ph.D. Thesis (Katholieke Universiteit Leuven)
- (36) Van der Linden, R. A. M., & Hood, A. W. 1999, A&A, 346, 303
- (37) Van Leer, B. 1997, J. Comput. Phys., 135, 229
- (38) Vekstein, G. E., Priest, E. R., & Steele, C. D. C. 1993, ApJ, 417, 781
- (39) Velli, M., Lionello, R., & Einaudi, G. 1997, Sol. Phys., 172, 257
- (40) Wilkins, M. L. 1980, J. Comput. Phys., 36, 281
- (41) Withbroe, G. L., & Noyes, R. W. 1977, Ann. Rev. Astr. Ap., 15, 363
- (42) Yeates, A. R., Hornig, G., & Wilmot-Smith, A. L. 2010, Phys. Rev. Lett. 105, 8
- (43) Zhang, M., & Low, B. C. 2003, ApJ, 584, 479