Coronal heating by the partial relaxation of twisted loops
^{1}
School of Mathematics and Statistics, University of St.
Andrews,
North Haugh, St. Andrews, KY16 9
SS Fife
UK
email:
michaelb@mcs.stand.ac.uk
^{2}
Jodrell Bank Centre for Astrophysics, Alan Turing Building, School
of Physics and Astronomy, University of Manchester, Oxford Road, M13 9 PL Manchester, UK
Received:
31
May
2012
Accepted:
24
November
2012
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 (T ≳ 10^{6} K).
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 threedimensional (3D) MHD Lagrangianremap code is used to simulate the evolution of a linetied 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 piecewiseconstant current profile, featuring three parameters. Initially, all configurations carry zeronetcurrent fields and are in ideally unstable equilibrium. The simulation results are compared with the predictions of helicityconserving 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 analyticallydetermined 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 heatingevent distributions produced by ensembles of kinkunstable loops. The energy release increases with relaxation radius; but, once the loop has expanded by more than 50%, further expansion yields little more energy. We conclude that the relaxation methodology may be used for coronal heating studies.
Key words: instabilities / magnetic fields / magnetic reconnection / magnetohydrodynamics (MHD) / plasmas / Sun: corona
© ESO, 2013
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 (β ≈ 0.01). Coronal loops are closed structures that emerge from the photosphere at one location and reenter 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 forcefree equilibria, which can be expressed as ∇ × B = α(r)B; where r is a position vector, and α = (μ_{0}j·B)/ B^{2} is related to the parallel electric current density (i.e., α is a measure of the twist).
Energy release might be triggered when a looplike 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 (L≈ 50 Mm) and exist within a highly conductive environment: therefore, in the absence of instability, magnetic fields diffuse slowly (t_{d} ≈ 10^{3} yr). 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 threedimensional (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 10^{7} erg cm^{2} s^{1}, 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 largescale solar flares (~10^{34} erg). 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 computationallyintensive 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 forcefree field; i.e., α(r) is invariant and ∇ × B = αB. The helicity (K) indicates how intertwined a magnetic field is with itself (Berger 1999). Coronal loops have a nonzero normal flux at the footpoints; hence, the gaugeinvariance of relative helicity (Berger & Field 1984; Finn & Antonsen 1985) makes it the more useful property: (1)Where A is the magnetic potential, B′ is the potential field with the same boundary conditions and A′ is the corresponding vector potential. Helicityconserving 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 heatingevent 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 zeronetcurrent 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 currentneutralisation 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 crosssection and therefore remains purely axial.)
A longstanding 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 kinkunstable 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 δK); 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 secondorder accurate predictorcorrector 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 divergencefree condition (∇·B = 0) is maintained to machine precision by Evans & Hawley (1988) constrained transport.
LARE3D solves the resistive MHD equations given by
with specific energy density, (6)where ρ is the mass density, v is the plasma velocity, B the magnetic field, P the thermal pressure, η is the resistivity (not magnetic diffusivity), J is the current density, γ = 5/3 is the ratio of specific heats, and μ_{0} is the magnetic permeability. Viscous heating is represented by the last term of Eq. (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 c_{f} is the fast magnetoacoustic speed, l is the distance across a grid cell in the direction normal to the shock front, s is a similarly localised strain rate (the subscripts i and j denote the different spatial coordinates), and the artificial viscosity coefficients ν_{1} = 0.1 and ν_{2} = 0.5 are constants. The form of Eq. (8) is derived from the RankineHugoniot jump conditions and the values of the coefficients have been chosen such that numerical oscillations behind shock fronts are prevented. Note, the force Eq. (3) 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, R_{b} is the loop radius, ρ_{0} the initial mass density, and B_{1} the initial axial field at r = 0. The other variables are expressed in a similar manner; where L^{∗} = 20 R_{b} is the loop length, t_{A} = R_{b}/v_{A} is the radial Alfvén transit time, the Alfvén speed, and the magnetic pressure. The specific energy density, current density and resistivity (ϵ, J and η) also have reference variables that can be expressed in terms of R_{b}, ρ_{0} and B_{1}: Values appropriate for a coronal active region can be obtained by setting R_{b} = 1 Mm, ρ_{0} = 1.6726 × 10^{13} kg m^{3} and B_{1} = 50 G. Hence, the length becomes 20 Mm, v_{A} ≈ 10 Mm s^{1} and η_{0} ≈ 4π × 10^{6} Ω m.
The resistivity is taken to be nonuniform in these simulations, where η_{b} is the background resistivity (normally set to zero, since actual coronal resistivities are approximately μ_{0} m^{2} s^{1}) and η_{c} = 0.001 is the anomalous resistivity, which is only switched on when the current reaches or exceeds J_{crit} = 15. The value of J_{crit} is set so that it is significantly higher than the maximum current at the start of the simulation. Supercritical 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 smallscale 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 L_{x} = L_{y} = 6 (− 3:+3) and L_{z} = 20 (−10:+10). Initially, the loop axis follows the zaxis and the loop radius is r = 1; therefore, the simulated loops all have an aspect ratio of 20. The loop is linetied at z = −10, +10, which means, at those zcoordinates, 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: 128^{2} × 256 (low) and 256^{2} × 512 (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 linetied coronal loop, see Gerrard et al. (2002), Gerrard & Hood (2003), Browning et al. (2008) and Hood et al. (2009).
Fig. 1 Schematic of a straightened coronal loop in the x  z plane (left) and in the x  y plane (right). The loop, comprises a core (dark grey), an outer layer (light grey) and a current neutralisation layer (blue); the whole loop is embedded in a rectangular potential envelope. The core radius is half the loop radius (R_{1}:R_{2}:R_{b}:R_{B} = 0.5:0.9:1:3, where R_{B} is the distance from the initial loop axis to the edge of the envelope). The loop aspect ratio (L/ R_{b}) in this figure is 20. 

Open with DEXTER 
The initial equilibrium model used in the latter two studies was extended by Bareford et al. (2011) to include an outer currentneutralising 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 (B_{θ}) always falls to zero at the loop boundary (R_{b}); therefore, B_{θ} is zero in the potential envelope. The loop’s radial αprofile is approximated by a piecewiseconstant function featuring three parameters (Fig. 1): the ratio of current to magnetic field is α_{1} in the core, α_{2} in the outer layer, α_{3} in the neutralisation layer and zero in the potential envelope. The free parameters are α_{1} and α_{2}, whereas α_{3} 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 R_{B} = 3 (three times the loop radius); this is far enough away that the boundary conditions do not influence the plasma evolution.
αprofiles, axial fluxes and field coefficients for the simulated loops.
The fields are expressed in terms of the wellknown 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:
where (i = 1,2,3) represents the sign of α_{i}. The fields must be continuous at the inner radial boundaries, R_{1}, R_{2} and R_{b}. (The positions are R_{1} = 0.5, R_{2} = 0.9 and R_{b} = 1, 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 R_{2} and R_{b}.) Therefore, the coefficients B_{j} and C_{j} (j = 2,3,4) are determined by the requirement of continuity of the magnetic field at all interfaces (Bareford et al. 2011). The value of α_{3} (the neutralisation layer current) is found, for a given (α_{1}, α_{2}), by numerical solution of B_{3θ}(R_{b}) = 0, ensuring that the net current is zero and that the azimuthal field vanishes outside the loop, see Eq. (14). From the nondimensionalisation of the magnetic field, the field coefficient at the core is B_{1} = 1.
The linear kink instability threshold for this currentneutralised 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 linetied αconfigurations.
Fig. 2 Linear kink instability thresholds for L/ R_{b} = 20. These thresholds have been cropped by a pair of dashed lines that indicate where B_{z}(r) starts to acquire a mixed polarity. The right threshold is sampled by a selection of five marginally unstable configurations (black circles). 

Open with DEXTER 
In contrast to the stability space for a loop of net current (Bareford et al. 2010), the instability threshold is open, see Fig. 2. This is very similar to the threshold found for a closefitting 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 α_{1} and α_{2} are too large. This reversal is unphysical and hence, introduces a restriction to the stable region of Fig. 2; namely, all stable configurations should have B_{z}(r) of uniform sign. The new stability space is closed by excluding those field profiles that have a B_{z}(r) of mixed polarity, since this could not be achieved directly by footpoint motions of an initially unidirectional field. The filled circles of Fig. 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 Fig. 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 δα_{1} ≈ δα_{2} – this is the reason why the part of the threshold curve where α_{2} > 0 is more finely sampled compared to α_{2} < 0. Note, the section of threshold on the left of Fig. 2 is merely the negative of the section on the right; hence, it is not necessary to sample both threshold sections.
Fig. 3 Initial analyticallydetermined B_{z} (solid) and B_{θ} (dashed) profiles for L_{uni} and L_{mix}. Field profiles for the stable loop (last row of Table 1) are also shown. 

Open with DEXTER 
The configurations listed in Table 1 encompass two types of loop; one class where α_{1} and α_{2} 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 L_{uni} (uniform twist) and Loop D will be denoted by L_{mix} (mixed twist).
Each of the loops indicated in Fig. 2 is subjected to a disturbance of the form, (17)where v_{r} is the radial component of the perturbed velocity, L is the loop length, is the radial coordinate, k = 1.1 is the wave number, θ = arctan(x/y), and the constant c_{1} = 0.01 reduces the amplitude so that the perturbation is initially linear. This disturbance initiates the kink instability.
Both L_{uni} and L_{mix} 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 L_{uni} (α_{1} = 2.25, α_{2} = 1.5) is linearly kink unstable, and the numerical simulation (Fig. 4, middle row) shows that it is also nonlinearly unstable. The magnetic energy, W, has an initial (dimensionless) value of 155.3 when integrated over the entire simulation domain, the internal (U) and kinetic (E) energies are also volume integrals, The initial magnetic energy integrated over the loop L_{uni} only is much less (W_{b} ≈ 20). All magnetic energies can be dimensionalised by making a simple correction to Eq. (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 nondimensional axial flux over the radius R_{B} (Table 1, Col. 5). For loop L_{uni}, this gives a multiplier of 4.8 × 10^{19} (assuming a radius of 1 Mm and a mean axial field strength of 50 G): thus, the dimensionalised initial loop energy is roughly 10^{21} J.
At the onset of instability, W undergoes a decrease, coincident with a rise in U and with a much more modest increase in kinetic energy (E_{max} ≈ 0.2). The maximum current, J_{max}, also rises just before the decrease in W, indicating the formation of a helical current sheet. The nonlinear instability starts at around t = 50t_{A} and within the next 50t_{A}, approximately 70% of the total energy release has been achieved. Magnetic energy reduces more slowly after t = 100t_{A}. 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, spatiallyconfined 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 Fig. 4, ideal MHD and zero background resisitivity, it can be seen that there is no detectable change in energy until around 40t_{A} – 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 Fig. 4 (η_{b} = 10^{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.
Fig. 4 Loop L_{uni}: the temporal variation in energy (magnetic, internal and kinetic), heating (Ohmic and viscous) and maximum current for ideal MHD (top row), for resistive MHD with η_{b} = 0 (middle row), and for η_{b} = 10^{4} (bottom row). The initial magnetic energy has been subtracted from the magnetic energy plots (solid lines, left column). The critical current (J_{crit} = 15) is indicated by the grey dashed horizontal lines (right column). For the η_{b} = 0 case, the maximum current plots are from the high (black) and low (grey) resolution simulations. 

Open with DEXTER 
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, Fig. 4 also shows that energy is not conserved for the two cases mentioned (ideal MHD and η_{b} = 0). At high resolution, only ~68% 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 10%. 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 nonzero background resistivity (η_{b} = 10^{4}) 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 η_{b} = 0. 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 η_{b} > 0: the internal energy increase is now approximately 96% of the magnetic energy decrease (the final kinetic energy is less than 0.2% of δW). A background resistivity of 10^{4} 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 η_{b} is halved, the percentage of the magnetic energy release lost to the simulation increases to around 10%. Further tests have shown that energy conservation is not improved by a lowering the critical current and thereby causing the anomalous resistivity (η_{c}) 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 L_{uni}; but η_{b} = 10^{4} is not realistic by coronal standards, and Fig. 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 η_{b} = 0: approximately 32% 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.
Fig. 5 Loop L_{mix}: the temporal variation in energy (magnetic, internal and kinetic), heating (Ohmic and viscous) and maximum current (low and high resolution) for resistive MHD with η_{b} = 0. The different plot lines follow the same scheme as that used for Fig. 4. 

Open with DEXTER 
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 (α_{1} = 0.5 and α_{2} = 0.1, see Table 1, bottom row) that is well within the stable region (as shown in Fig. 2). In the absence of an instability (and any applied resistivity), the energy declines by around one thousandth of one percent over 300t_{A}; this reduction is three orders of magnitude smaller than the energy release caused by the kink instability.
Loop L_{mix} (α_{1} = 2.54, α_{2} = −1.0) 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 L_{uni}. Note, the initial magnetic energy is higher than the value given for the other loop. This is a consequence of setting B_{1} = 1: L_{mix} 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 256^{2} × 512 simulations. The results for the other loop configurations on the threshold (Fig. 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 currentdependent 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 L_{uni} and L_{mix} 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 (Sect. 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 ρ L (ν_{1} c_{f} + ν_{2} L s ) – σ^{∗} 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 L_{uni} than for L_{mix}. In terms of azimuthal field, L_{mix} is the weaker of the two (Fig. 3), however, some importance should be attached to the fact that L_{mix} 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 α_{1} and α_{2} but keep the opposite signs, such that the total azimuthal field strength is comparable to L_{uni} (e.g., Loop E), the energy release increases only slightly (δW ≈ 2.6).
Figure 7 shows cross sections of L_{uni} at z = 0 (halfway along the loop) and L_{mix} at z = −2, which is roughly the centre of the only patch of significant viscous heating for L_{mix} at t = 60t_{A}, see bottom row of Fig. 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 Fig. 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 (Δt ≈ 50t_{A}) 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 Fig. 8; note the change in scale for the x and y 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 Fig. 8, the final state calculated by the numerical simulation is overlaid with magnetic field vectors in the xy plane, which are consistent with a cylindrical configuration, bounded by a currentneutralising layer: the arrows follow each other and the arrow sizes initially increase away from the axis, and then diminish before the loop edge.
Fig. 6 Magnetic field lines originating from the bottom left footpoint (dark grey) and from the upper right footpoint (light grey) are shown at t = 60t_{A} for L_{uni} (top) and L_{mix} (bottom). At the onset of instability, two plots are shown: one with isosurfaces of current (left) and the other with isosurfaces of  σ^{∗} , a proxy for viscous heating (right). 

Open with DEXTER 
Fig. 7 Spatial variation of current magnitude (left) and a viscous heating proxy (right) across the loop cross section at the apex (i.e., where z = 0) for L_{uni} (top) and at z = −2 for L_{mix} (bottom). 

Open with DEXTER 
Fig. 8 Spatial variation of current magnitude across the loop cross section at the final time of the simulation for L_{uni} (left) and L_{mix} (right). All currents are now well below the critical value. The arrows are magnetic field vectors. 

Open with DEXTER 
Loop L_{mix} expands less than L_{uni}, which is possibly due to the fact that the latter releases more energy. A single weaklytwisted 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 currentfree 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 x and y axes extend between − 3 and +3; hence, the integral limits given above. The photospheric boundaries are located at the limits of the z axis (z = −10, +10) and z = 0 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 gaugeinvariant vector potential can be specified as, (21)by subtracting the helicity due to the potential field, and Eq. (21) can be reexpressed by expanding the cross product of the second term, (22)Finally, the gaugeinvariant magnetic helicity is (23)The geometry used by DeVore differs significantly from that used here (Fig. 1), which features two separate photospheric boundaries at the limits of the z 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 z coordinate can simply be set to − 10.
Equations (18)–(23) have been implemented, using the fivepoint NewtonCotes integration formula. The marginally unstable loops (Fig. 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.
Helicity at three times during the simulations (η_{b} = 0) for all five kinkunstable loops.
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 ~2%.
Next, we calculate the ratio of the helicity variation to the change in magnetic energy, (24)both ΔK and ΔW are weighted by initial value. For four of the five loops, ΔK is around one order of magnitude lower than ΔW (Table 2, fifth column); these results are comparable with Browning et al. (2008). The relationship ΔK ≪ ΔW implies that magnetic energy dissipation is taking place on small spatial scales (i.e., shock fronts). One of the loops (E) has ΔK > ΔW; however, it is no coincidence that this loop also has the lowest helicity (K_{initial} = 1.16). The coarseness of the grid sampling used to calculate Eq. (23) – the x and y dimensions are sampled at every fourth cell and the z dimension at every other cell – is too great for such a small helicity. Hence, for this last loop, the grid sampling is increased from 4 × 4 × 2 to 2 × 2 × 2, which gives ΔK/ΔW ≈ 0.44. 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 r = 0 and r = 3 (R_{B}). The relaxed alpha was identified by assuming that ψ (axial flux) and K/ψ^{2} (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: R_{b} ≤ r ≤ R_{B}. Some form of partial relaxation is more likely to be relevant to the zeronetcurrent 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 (r = R_{l}). 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 R_{l} and R_{B}. Naively, it might have been expected that the partiallyrelaxed state would consist of a constant α field embedded directly in a potential field (α = 0) with no layer of reversed current. However, such an arrangement does not match the simulation results, see Sect. 4.1.1. Furthermore, it would be unphysical for the partiallyrelaxed 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 azimuthal field is zero outside the loop. The existence of current sheets (here broadened into a narrow currentneutralising layer) bounding localised relaxed states has also been noted by Gimblett et al. (2006). The axial flux is conserved such that ψ_{l} of the threshold state is equal to ψ_{l} of the relaxed state. The subscript l denotes the relaxation radius, R_{l}; it is the radial upper bound over which the associated property is calculated – the lower bound being the axis. For example, ψ_{l} is the axial flux from r = 0 to r = R_{l} (similarly, K_{B} is the helicity from r = 0 to r = R_{B}, the outer edge of the potential envelope). This method of conservation is illustrated by Fig. 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 (ψ_{l}) matches the initial loop flux (ψ_{b}), or it matches the initial flux within the radius (R_{l}) 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 Sect. 4.1.1). If the loop did not reconnect with its surroundings and expand to a radius of 1.8 R_{b} (as for L_{uni}), then the axial field would drop by a factor of 1.8^{2}, which is not observed. (The slight mismatch between the numerical and analytical axial field profiles seen in Fig. 11 is probably due to some limited free expansion of the loop.)
Fig. 9 Unstable and relaxed loop. K/ψ^{2} is conserved over the region enclosed by the dashed circle, which is the outer boundary (R_{l}) of the relaxed loop. 

Open with DEXTER 
The relaxation alpha (α_{l}) is determined by ensuring that the axial flux and helicity integrated over the dashed circle in Fig. 9 (left), match the values obtained when these same properties are integrated over 0 ≤ r ≤ R_{l} of the relaxed loop (Fig. 9, right). Also, helicity is absent from the potential envelope surrounding a zeronetcurrent loop; hence, K_{b} = K_{B} (unlike magnetic energy, W_{b} ≠ W_{B}). Thus, we do not need to choose which value of helicity to use. Once α_{l} is known, the energy release can be determined analytically; (25)where W_{l}(α_{i1},α_{i2}) is the energy of the threshold state and W_{l}(α_{l}) is the energy of the relaxed state.
Fig. 10 Variation with relaxation radius (R_{l}) of relaxed alpha (α_{l}, left) and of dimensionless energy release (δW, right) for L_{uni} (solid line) and for L_{mix} (dashed line). 

Open with DEXTER 
The relaxation radius (R_{l}) remains as a free parameter, although sensible values can be inferred from the numerical results. Figure 10 shows how α_{l} and δW vary with relaxation radius. As expected, these figures show an inverse relationship between α_{l} and δW. In general, δW increases with R_{l}, however, this relationship is not linear; beyond a moderate expansion of 50% (R_{l} = 1.5) the energy release is ~99% of its maximum for L_{mix} and ~80% for L_{uni}. This property of diminishing returns is also true for the other three loops indicated in Fig. 2. This means that the energy release is insensitive to the choice of relaxation radius, so long as it is assumed that R_{l} ≳ 1.5.
4.1.1. Numerical relaxed states
The final numericallydetermined state is merely the last snapshot provided by the simulation; it is expected to be close to the analyticallydetermined relaxed state.
Fig. 11 Loop L_{uni}: a comparison between the B_{x} (top row), B_{y} (middle row) and B_{z} (bottom row) magnetic field profiles obtained numerically (red line) and analytically (black line). The latter is calculated from the α_{l} and R_{l} that best fit the numerical plot, which is taken from the final frame (t = 400t_{A}) of the high resolution LARE3D simulation (η_{b} = 0). The comparisons are done at y = 0 for different z coordinates, z = −5 (left column), 0 (middle column) and 5 (right column). 

Open with DEXTER 
Each loop is simulated for at least 300t_{A}; 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, where B_{θ}(r) = (α_{l}/α_{l}) B_{1}J_{1}( α_{l}r) and r ≤ R_{l}. We generate relaxed configurations for every value of R_{l} 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 chisquared value when compared with the numerical values for the same component. These field comparisons (B_{x}, B_{y} and B_{z}) are performed over the x dimension for a selection of fifteen y  z coordinate pairs (y ∈ {−1, −0.5,0, +0.5, +1 }, z ∈ {−5,0, +5 }) – the horizontal dashed lines of Fig. 13 (left) indicate the y positions of the field plots.
As a result of the kink instability, the loop axis will be shifted from the origin in the x  y plane. The new axis position will be zdependent and can be recalculated by assuming that the relaxed loop is current neutralised (Fig. 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 analyticalnumerical comparisons for L_{uni} at y = 0.
In general, the analytical plots, such as the ones in Fig. 11, show a good agreement with the numerics; however, the R_{l} and α_{l} that best fit the red lines are independently chosen for each of the fortyfive subplots (there are fifteen y  z positions and three field components). Figure 12 (left) shows the variation in these bestfit parameters: each symbol, by virtue of its position, gives the R_{l} − α_{l} pair that best fits the numerical data extracted at specific y and z coordinates within the simulation domain – the symbol type denotes which field component is being matched. The use of two xaxes (one for R_{l}, the other α_{l}) means that all fortyfive symbols can be clearly distinguished (the yaxis has no meaning beyond a simple sorting of the data points).
Fig. 12 Loop L_{uni}: the bestfit R_{l} − α_{l} pairs when η_{b} = 0 (left) and when η_{b} = 10^{4} (right). Black circles are for those best fits determined from B_{z} profiles, red plus signs are for B_{x} and blue crosses B_{y}. 

Open with DEXTER 
The derived averages are ⟨ R_{l} ⟩ = 1.76 ± 0.29 and ⟨ α_{l} ⟩ = 0.28 ± 0.31; these two properties have a onetoone mapping. The dimensionless energy release is 5.87 ± 1.11. Despite the large scatter for α_{l}, the deviation for δW is comparatively modest, this is because d(δW)/d(R_{l}) is small when R_{l} ≈ 1.8 (Fig. 10, right). Figure 12 suggests that the final numericallycalculated 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 α_{l} values that are further from ⟨ α_{l} ⟩ than those on the right, since the relationship between α_{l} and R_{l} is not linear (Fig. 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 α_{l} 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 rerunning the simulation with background resistivity; once again, η_{b} = 10^{4}. Figure 12 (right) yields ⟨ R_{l} ⟩ = 1.74 ± 0.11 and ⟨ α_{l} ⟩ = 0.23 ± 0.09. The mean dimensionless energy release becomes ⟨ δW ⟩ = 6.06 ± 0.32. The resistivity smooths out lowlevel noise and restricts the current, and thereby reduces the deviation associated with the analytical fit.
Analytical and numerical comparison for the kinkunstable loops A–E.
Fig. 13 Value of α (left) and plasmaβ (right) over the xy plane for L_{uni} (η_{b} = 0). The grey circle approximates the loop cross section at z = 0 and t = 400t_{A}. 

Open with DEXTER 
Fig. 14 Left, the cosine of the angle between the current and the magnetic field over the xy plane. Right, the cosine of the angle between the magnetic field and the gradient of the pressure over the same area. Positions outside the grey circle, representing the loop cross section, have not been assigned a colour. The background resistivity is zero. 

Open with DEXTER 
The overall impressive level of agreement demonstrated for L_{uni}, also extends to other positions along the instability threshold, see Table 3. This table also confirms that the analyticallycalculated helicities of Loops A–E are in approximate agreement with the numerical values, which were derived according to the procedure discussed in Sect. 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 (Sect. 3.1 and Fig. 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 Fig. 11 clearly indicates that there are deviations from a simple locallyconstant α profile. We show the value of α computed over the midplane (z = 0) of L_{uni}, see Fig. 13. A grey circle based on the current magnitude plot of Fig. 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, α_{l} = 0.28. The influence of the initial α_{3} (the αvalue for the current neutralisation layer) can be seen in the limits of the colour bar. In addition, plasmaβ (Fig. 13, right) has, compared to initial values of 10^{3}, 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 Fig. 14; note, positions outside the loop are left unplotted. The current is either parallel (white) or antiparallel (black) to the field for a high percentage of the cross sectional area (positions that are far from parallel, 20° < j∠B < 160°, 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 forcefree state (although Taylor relaxation would predict a uniform pressure distribution). However, it seems that the loop is heading towards an approximate balance of forces (j × B = ∇P): Fig. 14 (right) shows that the angle between the field and the pressure gradient is on average close to zero. Interestingly, the parallel and antiparallel areas are often found next to each other and therefore might be expected to cancel at later simulation times. The plots of Figs. 13 and 14 are qualitatively similar to those taken at other z coordinates that are not too close to the footpoints (e.g., − 5 < z < 5).
5. Summary and conclusions
A nonlinear 3D MHD code has been used to simulate the evolution of a set of zeronetcurrent 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 η_{b} = 10^{4}.) 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: slowmode 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 wellmodelled by a (localised) constantα profile with some fluctuations superposed – this suggests the finescale structure is selfcancelling (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 (Fig. 1) might have influenced the results. Hence, we also ran a resistive MHD simulation (with η_{b} = 0) for a zeronetcurrent equilibrium that possessed a smooth αprofile, see Case 3 of Hood et al. (2009). The constant parameter, λ = 1.62, was set such that the equilibrium was on the threshold of instability. We found that the simulation results were similar to those produced by L_{mix}. Numerically, the energy release was 1.2, which was again consistent with the analyticallydetermined value.
It appears that the assumption of a Taylorrelaxed 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 R_{l}, which is less than the full radius, R_{B}). 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 R_{l} can be precisely determined from the field configuration at instability onset. However, the analytical work has revealed that for marginallyunstable loops, the energy release varies little with relaxation radius once R_{l} ≥ 1.5 (Fig. 10); hence, a calculation of the energy release does not necessarily require a precise prediction for R_{l}.
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 Taylorrelaxed 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 predetermined velocity profile (v_{θ}) 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 resimulated 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 ρ(z)), 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.
Acknowledgments
We thank the referee for helpful comments and M.R.B. acknowledges financial support from STFC.
References
 Arber, T. D., Longbottom, A. W., & Van der Linden, R. A. M. 1999, ApJ, 517, 990 [NASA ADS] [CrossRef] (In the text)
 Arber, T. D., Longbottom, A. W., Gerrard, C. L., & Milne, A. M. 2001, J. Comput. Phys., 171, 151 [NASA ADS] [CrossRef] (In the text)
 Aschwanden, M. J., & Acton, L. W. 2001, ApJ, 550, 475 [NASA ADS] [CrossRef] (In the text)
 Baty, H. 2000, A&A, 360, 345 [NASA ADS] (In the text)
 Baty, H., & Heyvaerts, J. 1996, A&A, 308, 935 [NASA ADS] (In the text)
 Bareford, M. R., Browning, P. K., & Van der Linden, R. A. M. 2010, A&A, 521, A70 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Bareford, M. R., Browning, P. K., & Van der Linden, R. A. M. 2011, Sol. Phys., 273, 93 [NASA ADS] [CrossRef] (In the text)
 Berger, M. A. 1999, Plasma Phys. Control. Fusion, 41, B167 [NASA ADS] [CrossRef] (In the text)
 Berger, M. A., & Field, G. B. 1984, J. Fluid. Mech., 147, 133 [NASA ADS] [CrossRef] (In the text)
 Botha, G. J. J., Arber, T. D., & Hood, A. W. 2011, A&A, 525, A96 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Browning, P. K. 1988, J. Plasma Phys., 40, 263 [NASA ADS] [CrossRef] (In the text)
 Browning, P. K., & Van der Linden, R. A. M. 2003, A&A, 400, 355 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Browning, P. K., Sakurai, T., & Priest, E. R. 1986, A&A, 158, 217 [NASA ADS] (In the text)
 Browning, P. K., Gerrard, C., Hood, A. W., Kevis, R., & Van der Linden, R. A. M. 2008, A&A, 485, 837 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 DeVore, C. R. 2000, ApJ, 539, 949 (In the text)
 Dixon, A. M., Berger, M. A., Browning, P. K., & Priest, E. R. 1989, A&A, 225, 156 [NASA ADS] (In the text)
 Evans, C. R., & Hawley, J. F. 1988, ApJ, 332, 659 [NASA ADS] [CrossRef] (In the text)
 Finn, J. M., & Antonsen, T. M. 1985, Commun. Plasma Phys. Controll. Fusion, 9, 111 (In the text)
 Fletcher, L. 2009, A&A, 493, 241 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Gerrard, C. L., & Hood, A. W. 2003, Sol. Phys., 214, 151 [NASA ADS] [CrossRef] (In the text)
 Gerrard, C. L., Arber, T. D., & Hood, A. W. 2002, A&A, 387, 687 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Gimblett, C. G., Hastie, R. J., & Helander, P. 2006, Plasma Phys. Control. Fusion, 48, 1531 [NASA ADS] [CrossRef] (In the text)
 Haynes, M., & Arber, T. D. 2007, A&A, 467, 327 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Heidbrink, W. W., & Dang, T. H. 2000, Plasma Phys. Control. Fusion, 42, L31 [NASA ADS] [CrossRef] (In the text)
 Heyvaerts, J., & Priest, E. R. 1984, A&A, 137, 63 [NASA ADS] (In the text)
 Hood, A. W. 1992, Plasma Phys. Controll. Fusion, 34, 411 [NASA ADS] [CrossRef] (In the text)
 Hood, A. W., Browning, P. K., & Van der Linden, R. A. M. 2009, A&A, 506, 913 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Hudson, H. S. 1991, Sol. Phys., 133, 357 [NASA ADS] [CrossRef] (In the text)
 Melrose, D. B., Nicholls, J., & Broderick, N. G. 1994, J. Plasma Phys., 51, 163 [NASA ADS] [CrossRef] (In the text)
 Priest, E. W., Longcope, D. W., & Heyvaerts, J. 2005, ApJ, 624, 1057 [NASA ADS] [CrossRef] (In the text)
 Qiu, J. 2009, ApJ, 692, 1110 [NASA ADS] [CrossRef] (In the text)
 Srivastava, A. K., Zaqarashvili, T. V., Kumar, P., & Khodachenko, M. L. 2010, ApJ, 715, 292 [NASA ADS] [CrossRef] (In the text)
 Taylor, J. B. 1974, Phys. Rev. Lett., 33, 1139 [NASA ADS] [CrossRef] (In the text)
 Taylor, J. B. 1986, Rev. Mod. Phys., 58, 741 [NASA ADS] [CrossRef] (In the text)
 Van der Linden, R. A. M. 1991, Ph.D. Thesis (Katholieke Universiteit Leuven) (In the text)
 Van der Linden, R. A. M., & Hood, A. W. 1999, A&A, 346, 303 [NASA ADS] (In the text)
 Van Leer, B. 1997, J. Comput. Phys., 135, 229 [NASA ADS] [CrossRef] (In the text)
 Vekstein, G. E., Priest, E. R., & Steele, C. D. C. 1993, ApJ, 417, 781 [NASA ADS] [CrossRef] (In the text)
 Velli, M., Lionello, R., & Einaudi, G. 1997, Sol. Phys., 172, 257 [NASA ADS] [CrossRef] (In the text)
 Wilkins, M. L. 1980, J. Comput. Phys., 36, 281 [NASA ADS] [CrossRef] (In the text)
 Withbroe, G. L., & Noyes, R. W. 1977, ARA&A, 15, 363 [NASA ADS] [CrossRef] (In the text)
 Yeates, A. R., Hornig, G., & WilmotSmith, A. L. 2010, Phys. Rev. Lett., 105, 8 [NASA ADS] [CrossRef] (In the text)
 Zhang, M., & Low, B. C. 2003, ApJ, 584, 479 [NASA ADS] [CrossRef] (In the text)
All Tables
Helicity at three times during the simulations (η_{b} = 0) for all five kinkunstable loops.
All Figures
Fig. 1 Schematic of a straightened coronal loop in the x  z plane (left) and in the x  y plane (right). The loop, comprises a core (dark grey), an outer layer (light grey) and a current neutralisation layer (blue); the whole loop is embedded in a rectangular potential envelope. The core radius is half the loop radius (R_{1}:R_{2}:R_{b}:R_{B} = 0.5:0.9:1:3, where R_{B} is the distance from the initial loop axis to the edge of the envelope). The loop aspect ratio (L/ R_{b}) in this figure is 20. 

Open with DEXTER  
In the text 
Fig. 2 Linear kink instability thresholds for L/ R_{b} = 20. These thresholds have been cropped by a pair of dashed lines that indicate where B_{z}(r) starts to acquire a mixed polarity. The right threshold is sampled by a selection of five marginally unstable configurations (black circles). 

Open with DEXTER  
In the text 
Fig. 3 Initial analyticallydetermined B_{z} (solid) and B_{θ} (dashed) profiles for L_{uni} and L_{mix}. Field profiles for the stable loop (last row of Table 1) are also shown. 

Open with DEXTER  
In the text 
Fig. 4 Loop L_{uni}: the temporal variation in energy (magnetic, internal and kinetic), heating (Ohmic and viscous) and maximum current for ideal MHD (top row), for resistive MHD with η_{b} = 0 (middle row), and for η_{b} = 10^{4} (bottom row). The initial magnetic energy has been subtracted from the magnetic energy plots (solid lines, left column). The critical current (J_{crit} = 15) is indicated by the grey dashed horizontal lines (right column). For the η_{b} = 0 case, the maximum current plots are from the high (black) and low (grey) resolution simulations. 

Open with DEXTER  
In the text 
Fig. 5 Loop L_{mix}: the temporal variation in energy (magnetic, internal and kinetic), heating (Ohmic and viscous) and maximum current (low and high resolution) for resistive MHD with η_{b} = 0. The different plot lines follow the same scheme as that used for Fig. 4. 

Open with DEXTER  
In the text 
Fig. 6 Magnetic field lines originating from the bottom left footpoint (dark grey) and from the upper right footpoint (light grey) are shown at t = 60t_{A} for L_{uni} (top) and L_{mix} (bottom). At the onset of instability, two plots are shown: one with isosurfaces of current (left) and the other with isosurfaces of  σ^{∗} , a proxy for viscous heating (right). 

Open with DEXTER  
In the text 
Fig. 7 Spatial variation of current magnitude (left) and a viscous heating proxy (right) across the loop cross section at the apex (i.e., where z = 0) for L_{uni} (top) and at z = −2 for L_{mix} (bottom). 

Open with DEXTER  
In the text 
Fig. 8 Spatial variation of current magnitude across the loop cross section at the final time of the simulation for L_{uni} (left) and L_{mix} (right). All currents are now well below the critical value. The arrows are magnetic field vectors. 

Open with DEXTER  
In the text 
Fig. 9 Unstable and relaxed loop. K/ψ^{2} is conserved over the region enclosed by the dashed circle, which is the outer boundary (R_{l}) of the relaxed loop. 

Open with DEXTER  
In the text 
Fig. 10 Variation with relaxation radius (R_{l}) of relaxed alpha (α_{l}, left) and of dimensionless energy release (δW, right) for L_{uni} (solid line) and for L_{mix} (dashed line). 

Open with DEXTER  
In the text 
Fig. 11 Loop L_{uni}: a comparison between the B_{x} (top row), B_{y} (middle row) and B_{z} (bottom row) magnetic field profiles obtained numerically (red line) and analytically (black line). The latter is calculated from the α_{l} and R_{l} that best fit the numerical plot, which is taken from the final frame (t = 400t_{A}) of the high resolution LARE3D simulation (η_{b} = 0). The comparisons are done at y = 0 for different z coordinates, z = −5 (left column), 0 (middle column) and 5 (right column). 

Open with DEXTER  
In the text 
Fig. 12 Loop L_{uni}: the bestfit R_{l} − α_{l} pairs when η_{b} = 0 (left) and when η_{b} = 10^{4} (right). Black circles are for those best fits determined from B_{z} profiles, red plus signs are for B_{x} and blue crosses B_{y}. 

Open with DEXTER  
In the text 
Fig. 13 Value of α (left) and plasmaβ (right) over the xy plane for L_{uni} (η_{b} = 0). The grey circle approximates the loop cross section at z = 0 and t = 400t_{A}. 

Open with DEXTER  
In the text 
Fig. 14 Left, the cosine of the angle between the current and the magnetic field over the xy plane. Right, the cosine of the angle between the magnetic field and the gradient of the pressure over the same area. Positions outside the grey circle, representing the loop cross section, have not been assigned a colour. The background resistivity is zero. 

Open with DEXTER  
In the text 