## Abstract

An analysis of the importance of shock heating within coronal magnetic fields has hitherto been a neglected area of study. We present new results obtained from nonlinear magnetohydrodynamic simulations of straight coronal loops. This work shows how the energy released from the magnetic field, following an ideal instability, can be converted into thermal energy, thereby heating the solar corona. Fast dissipation of magnetic energy is necessary for coronal heating and this requirement is compatible with the time scales associated with ideal instabilities. Therefore, we choose an initial loop configuration that is susceptible to the fast-growing kink, an instability that is likely to be created by convectively driven vortices, occurring where the loop field intersects the photosphere (i.e. the loop footpoints). The large-scale deformation of the field caused by the kinking creates the conditions for the formation of strong current sheets and magnetic reconnection, which have previously been considered as sites of heating, under the assumption of an enhanced resistivity. However, our simulations indicate that slow mode shocks are the primary heating mechanism, since, as well as creating current sheets, magnetic reconnection also generates plasma flows that are faster than the slow magnetoacoustic wave speed.

## 1. Introduction

The energy that heats the solar corona is thought to be released from magnetic fields that cover thousands of kilometres rising from and closing back down in the photosphere. As a consequence of these megametre scales, and as a result of the near-perfect conductivity within the solar atmosphere, coronal magnetic fields diffuse slowly (*t*_{d}≈10^{3} years). This leads to the important property that the magnetic field is *frozen* into the plasma, thereby allowing convective motions at the photospheric boundaries to continually reconfigure coronal magnetic fields. Hence, the work done by convective flows is stored within the corona as free magnetic energy. As the time scale associated with photospheric convection is long compared to the Alfvén time, the coronal magnetic field can regain equilibrium each time more energy is stored. This means that when an instability does eventually occur, the amount of energy accumulated is sufficient to heat the surrounding plasma to million degree temperatures. The radiation associated with this heating reveals the general topology of coronal magnetic fields, which are commonly observed as bundles of discrete coronal loops above solar active regions.

This paper is concerned with how an ideal kink instability [1–4]—that we assume is the eventual consequence of continual photospheric driving—can lead to the heating of a coronal loop. Many ideal instabilities generate current sheets, which when accompanied by resistive effects allow a change in magnetic topology, increasing the amount of energy that can be released. We investigate the ideal kink [5] as its growth rate is fast enough to match flare-like time scales. In addition, there is a substantial body of work [6–11], based on numerical simulations, that demonstrates how the kink instability creates sheets of high current density, facilitating magnetic reconnection.

The mechanism whereby magnetic energy is converted into heat has often been attributed to Ohmic dissipation [10–16]; usually, this is implemented by *switching on* an increased resistivity whenever the current exceeds some critical value. However, recent work involving coronal loop simulations [17] has revealed a strong connection between plasma heating and the shocks resulting from a kink instability. This is perhaps not surprising, since for several decades now, slow mode shocks [18] have been known to be an important feature of Petschek (or fast) reconnection [19], where the shortening of the current sheet allows shocks to form in the outflow region, enhancing the diffusion of magnetic energy. Much research has been done to investigate the critical parameters responsible for fast reconnection [20–22]. In addition, other workers have observed slow mode shocks within large-scale numerical simulations of the solar environment [23–25], and even observational data have provided strong evidence for shocks located at flare sites [26]. A common feature of numerical reconnection experiments is the use of a driven (or perturbed) current sheet, and, for models of the solar atmosphere, reconnection results from allowing magnetic flux to rise into the corona. In this work, we show that the conditions for reconnection and slow mode shocks can also occur when a coronal loop undergoes a kink instability. As coronal fields are line-tied at the photospheric boundary, vortical motions associated with convection will repeatedly twist the loop until it becomes kink unstable [27,28].

The paper is structured as follows. Section §2 describes the loop model in terms of its initial magnetic field. Next, we present the numerical code and focus on how well it performs with regard to shock capturing. We also discuss an improved implementation of anomalous resistivity, where, instead of an arbitrary constant, we use a variable critical current based on local plasma conditions. Our aim here is to rule out the possibility of Ohmic heating being underestimated in numerical experiments. Section §4 discusses the simulation results: specifically, when and where energy is released from the magnetic field, and how this energy is thermalized. Lastly, the results are summarized and our conclusions are given. The instructions, scripts and source code required to reproduce the results presented in this paper are held in an online GitHub repository [29].

## 2. Loop model

A coronal loop is modelled as a straight cylinder with a loop axis that is initially coincident with the *z*-axis. Assuming the magnitude of the magnetic field on the axis, *B*_{0}, as a typical field strength and the radius of the cylinder, *R*_{0}, as a typical length, we use the magnetic field equations given by Hood *et al.* [13] to define a dimensionless initial loop equilibrium that has *zero* net current
2.1
and
2.2
where *r* is in the range 0≤*r*≤1, and *λ* is a constant. The somewhat complicated expression for *B*_{z} is actually derived from and (∇×** B**)×

**=0 (i.e. a force-free equilibrium is assumed). Thus, an expression for the ratio of current density to magnetic field,**

*B**α*=(

**⋅**

*j***)/(|**

*B***|**

*B*^{2}), which is continuous throughout the loop volume, can be determined, 2.3

These expressions, valid for 0≤*r*≤1, are illustrated by figure 1. Outside the loop, and *B*_{θ}=*α*=0. The footpoints of the loop are located at the *z*-limits (|*z*|=10) and so the loop apex is at *z*=0 (the azimuthal plane is synonymous with the *x*–*y* plane). The radial boundary is defined by the extent of the azimuthal field; this extent is itself an outcome of the spatial scale of the convective motions that generate the currents (and thereby azimuthal field). Essentially, localized footpoint motions will always define a loop that has zero net current. This loop configuration is linearly kink unstable when *λ*≥1.62. We use *λ*=1.8, as it is fully unstable, and allows us to compare our results with previous work [13,30].

At the start of the simulation, a weak perturbation is applied to the loop,
2.4
where *v*_{r} is the radial component of the perturbed velocity, *L*=20 is the loop length, *k*=0.3*π* is the wavenumber and . This disturbance initiates the kink instability and its form is constructed such that the amplitude of the perturbation falls off quickly with radius.

## 3. Numerical code

The numerical simulations are performed using a three-dimensional magnetohydrodynamic (MHD) code called LARE3D [31]. This code follows a **La**grangian **re**map scheme, hence its name. Essentially, the Lagrangian part, which is done using a second-order accurate predictor–corrector method, deforms the grid such that it moves with the plasma. The advantage of this technique is that additional physics (e.g. shock capturing) can easily be incorporated into the code. The remap stage involves the mapping of the plasma properties back to the original Cartesian grid; monotonicity is preserved through the use of van Leer [32] gradient limiters.

LARE3D solves the following resistive MHD equations:
3.1
3.2
3.3
3.4
with specific energy density,
3.5
where *ρ* is the mass density, ** v** is the plasma velocity,

**is the magnetic field,**

*B**P*is the thermal pressure,

*η*is the resistivity,

*J*is the current density,

*γ*=5/3 is the ratio of specific heats and

*μ*

_{0}is the magnetic permeability. The last terms of equations (3.2) and (3.4) feature tensors and are required to simulate the effect of hydrodynamic shocks; these are defined in §3a.

We normalize the variables in the MHD equations using dimensional values suitable for a coronal active region
3.6
where asterisk subscripts denote the normalized MHD variables. For example, *R*_{0}=1 Mm is the loop radius. The density and temperature are initially uniform (*ρ*_{0}=1.67×10^{−12} kg m^{−3} and *T*_{0}=20 000 K), and *B*_{0}=0.005 T is the initial axial field strength. Thus, at *t*=0, the plasma is dominated by the magnetic field: plasma-*β* is of order 10^{−5}. Other variables are expressed in a similar manner
3.7
where *L*=20 *R*_{0} is the loop length, *t*_{A}=*R*_{0}/*v*_{A} the radial Alfvén transit time, the Alfvén speed and the magnetic pressure. Hence, the length becomes 20 Mm, *t*_{A}≈0.3 s and *v*_{A}≈3.4 Mm s^{−1}.

The computational domain is a three-dimensional staggered grid: physical variables are not calculated at the same place for each cell in the domain, which improves numerical stability and allows conservation laws to be included in the computation. The domain size is *L*_{x}=*L*_{y}=4 (−2:+2) and *L*_{z}=20 (−10:+10). The loop is line-tied at *z*=−10, +10, which means, at those *z*-coordinates, the velocity components are set to zero. The velocity components are also zero at the boundaries for the other two directions. In addition, the normal derivatives of magnetic field, energy and density are zero at all boundaries. The simulations were run at three grid resolutions, 128^{2}×256 (low), 256^{2}×512 (medium) and 512^{2}×1024 (high). Unless stated otherwise figures 2–10 are taken from the medium resolution run. Using three different resolutions allows us to investigate whether the shock heating is real or merely an artefact of inadequate spatial resolution.

### (a) Shock resolution

The dominance of shock heating within simulations of kink-unstable loops [12,13,17] requires us to explain how this important mechanism is implemented. LARE3D uses shock viscosity [33] to capture the heating effect of shocks, which is implemented by adding an extra term to the energy equation. This term is expressed as the product of the rate of strain tensor,
3.8
and the shock tensor,
3.9
where ** v** is the velocity, is the magnetosonic speed,

*l*

_{s}is the distance across a grid cell in the direction normal to the shock front,

*s*is a similarly localized strain rate (the subscripts

*i*and

*j*denote the different spatial coordinates) and the shock viscosity coefficients

*ν*

_{1}=0.1 and

*ν*

_{2}=0.5 are constants. The terms within the rightmost set of brackets in equation (3.9) evaluate to zero when a cell undergoes uniform compression/expansion or rigid rotation. If we approximate these terms (and |

*s*|) by ∇

**/**

*v**l*

_{s}, 3.10 where the linear term damps out the non-physical oscillations caused by finite differencing in the numerical code and also dissipates kinetic energy. More of the dissipation is handled by the quadratic term as the shock increases in strength. The

*form*of equation (3.10) can be derived from the Rankine–Hugoniot jump conditions; furthermore, if one neglects terms higher than second order, it is possible to obtain values of 1 and 2/3 for

*ν*

_{1}and

*ν*

_{2}, respectively. However, equation (3.9) is expressed in terms of tensors, which have been included in such a way so as to prevent non-physical results under MHD conditions. The energy dissipation associated with each grid cell depends on the level of non-uniform shearing. Note that there is no

*real*viscosity in the simulations presented here. The viscosity part of shock viscosity comes from the mathematical form used to capture shock heating. Shocks create steep velocity gradients and the heating implied by these gradients is captured by the viscous-like terms of equations (3.9) and (3.10).

Arber *et al.* [31] tested the shock-capturing capabilities of the LARE code (with *ν*_{1}=0.1 and *ν*_{2}=1) by running a series of numerical experiments based on the work of Brio & Wu [34]. A shock is initialized by having discontinuities in density, pressure and magnetic field along some boundary that divides a two-dimensional plane. We have repeated these tests using the values for the shock viscosity coefficients given here (*ν*_{1}=0.1, *ν*_{2}=0.5) at a resolution of 256×256. As with Arber *et al.*, we obtained results that agreed closely with those of Brio & Wu. Such comparisons are useful because Brio & Wu used a second-order upwind scheme, which avoids the problem of false oscillations behind the shock front.

### (b) Critical current density

In numerical simulations of coronal loops, the background resistivity is usually set to zero in order to represent the near perfect conductivity of the coronal environment. Ohmic heating only occurs when an anomalous (non-zero) resistivity is applied in response to the current density reaching a threshold value. Using an arbitrary value for the threshold is problematic when trying to determine the significance of Ohmic heating compared to other forms of heating. A high current density (at least double that of the maximum initial value) is assumed to be caused by a reduction in the local conductivity, which either arises from a small-scale turbulent process or is a consequence of collisionless reconnection. Thus, anomalous resistivity is meant to compensate for some physical process that cannot be modelled by LARE3D and this process may or may not be collisional. A better approach, considering the rarity of coronal plasma, is to set the threshold current density according to the current required for local conditions to be consistent with a collisionless plasma.

Gordovskyy *et al.* [28] showed how a spatially varying critical current density could be derived by matching the local electron drift speed with the local ion sound speed. We achieve the same result by following a slightly different derivation in order to make clear which combination of plasma properties is directly proportional to *j*_{crit}:
3.11
where *e* is the electron charge, *n*_{e} is the electron number density, *γ* is the ratio of specific heats, *Z* is the ion proton number, *k*_{B} is the Boltzmann constant, *T*_{e} is the electron temperature and *m*_{i} is the ion mass. When the two speeds are equal, the ions and electrons are about to separate, and this will lead to the higher (anomalous) resistivities associated with *collisionless* reconnection. For an ionized hydrogen plasma, *Z*=1 and *m*_{i}=*m*_{p}, and, since the plasma is still (although only just) a single fluid, *T*_{e}=*T*. Equation (3.11) therefore simplifies to
3.12
which can expressed in code normalized units,
3.13
through the use of the appropriate normalizing constants,
3.14
Note that the asterisk subscripts have been dropped from equation (3.13). Hence and using previously mentioned values, together with *γ*=5/3 and *m*_{f}=1, the proportionality constant is approximately 1.6×10^{8}.

There is the possibility that finite spatial resolution will underestimate the current density. We can correct for this by assuming that a typical current sheet thickness is equal to the Larmor proton radius,
3.15
which is the ratio of the thermal proton speed and the Larmor proton frequency (*B* is the local magnetic field strength). For a typical grid size of 4^{2}×20 and a resolution of 128^{2}×256, the grid cells have a minimum length of 4/128 Mm≈31.2 km: thus, the current can said to be underestimated by a factor (3.12×10^{4})/*r*_{lp}. The Larmor proton radius is proportional to (the constant of proportionality is roughly 7.2, again using quoted values). So, this extra factor required to account for finite resolution also depends on local conditions,
3.16
where *δx* is the minimum grid cell length in code normalized units (e.g. *δx*=4/128). Substituting in equation (3.15) gives
3.17
The scheme outlined here is intended to maximize the opportunities for Ohmic dissipation: we are making the strong assumption that the current density is always underesolved. Furthermore, our use of the ion and not the electron sound speed gives us a lower current threshold.

## 4. Numerical results

### (a) Energy and heating

In order to show that the shock heating is not a numerical artefact, we investigate the impact of increasing the number of grid cells—the colour of the plot line indicates the spatial resolution (figure 2). There is little difference in the amount of energy released by the magnetic field, which indicates that the energy release is entirely defined by the nature of the unstable equilibrium. The two most highly resolved simulations achieve instability about 20*t*_{A} earlier than the low-resolution case. Numerical dissipation (*E*_{dis}=*E*_{tot}(*t*)/*E*_{tot}(0)−1, where *E*_{tot}=*W*+*E*_{kin}+*E*_{int}) is lowest for the highest resolution case and is just under 0.5%. The amount of Ohmic heating reduces markedly with increasing resolution. Figure 3*a* shows the ratio of the total shock heating to the total Ohmic heating as a function of time. As *δx* is decreased, the critical current threshold rises (see equation (3.17)) and the opportunities for anomalous resistivity are much reduced. This relationship outweighs the fact that better resolved currents are typically higher (figure 3*b*).

What is significant about these results is that our attempts to maximize Ohmic dissipation by using a locally defined current threshold have instead confirmed the dominance of shock heating. Between the 128^{2}×256 and 256^{2}×512 simulations, the shock heating increases by about 10%; however, a further increase in spatial resolution leaves the level of shock heating more or less unchanged. The gradients in kinetic energy that LARE3D interprets as shock heating are not an artefact of inadequate resolution. Henceforth, figures and plots will be taken from the 256^{2}×512 simulation.

### (b) Shock heating

It is instructive to break down the shock heating into its constituent tensor components, where each component is identified by *ε*_{ij}−(1/3)*δ*_{ij}∇⋅** v** (equation (3.9))—the possible

*ij*pairs are

*xy*,

*xx*,

*yy*,

*xz*,

*yz*and

*zz*(figure 4). The shock heating is mostly due to the spatial variation of the

*x*and

*y*components of velocity in the azimuthal plane. Roughly, one-quarter of the heating is due to changes in

*v*

_{x}and

*v*

_{y}with respect to

*z*, and very little is contributed by the

*zz*component.

Figure 5 shows the distribution of shock heating at the beginning of the instability (*t*=95*t*_{A}). The shock heating is concentrated where the loop has kinked outwards against the external potential field. The curved bands of figure 5, which are indicative of Petschek reconnection, curve inwards towards the loop axis as one increases or decreases the *z* coordinate. Figure 5*a* is symmetrical about *y*=0; thus the top right quadrant is expanded so as to follow the shock heating in more detail. There are in fact two thin bands of shock heating that appear to trace the azimuthal magnetic field—between the two curved bands the shock heating is minimal. At this time, there are no super-critical currents around the midplane: Ohmic dissipation is therefore non-existent within figure 5.

If the dissipation of magnetic energy is responsible for the rise in shock heating then these two events must be connected by the Lorentz force. At every time step, the Lorentz term is calculated and used to update the velocity via the equation of motion. Shock heating might occur if such updates create sharp spatial gradients of velocity. Figure 5*b* shows a strong correlation between the Lorentz vectors (projected in the *x*–*y* plane) and the areas of strongest shock heating. The inner band of heating appears to be caused by Lorentz forces that are generally acting in a radially outward direction, whereas a radially inward force coincides with the outer band. In between the two bands, the Lorentz force is minimal and so is the shock heating. The reason for the shock heating *gap* must relate to the lack of a Lorentz force. Within the gap, current and magnetic field magnitudes are non-zero; hence, the only way for the Lorentz force (and the shock heating) to become negligible is for the current and magnetic field vectors to become aligned.

The Lorentz forces are dominated by the *x* and *y* components: moreover, it is striking to see two sets of anti-parallel Lorentz vectors (corresponding to the two bands of shock heating) (figure 5*b*). If these results are compatible with Petschek reconnection, there should be two shock front normals pointing outwards from the gap, which outlines the current sheet thickness. The sheet width would have a component out of the plane (parallel to the *z*-axis). Figure 6*a* allows us to examine the plasma flow pattern. Note that the magnitude of the in-plane velocity is given by the length of the arrow, and the longest arrows correspond to a velocity of 0.11, the limit of the *v*_{z} colour bar. Thus, at *x*=*y*=0.5, there is a plasma flow where the azimuthal and axial components are roughly equal (*v*_{z}≈0.7), and so figures 5 and 6 present us with an oblique line of sight into a plasma outflow region. The slow mode speed varies according to local conditions, as it is dependent on temperature, density, field strength and velocity, all of which are non-uniform. For this reason, we plot the logarithm of the slow mode speed Mach number (figure 6*b*) in order to confirm the presence of supersonic flows. Note also that the shock heating contours mostly surround regions coloured white: the slowing of shocked plasma to subsonic speeds correctly coincides with regions of heating.

Further insight is gained by comparing the distribution of Lorentz vectors with two other properties, specifically, plasma-*β* (figure 7*a*) and the axial current (figure 7*b*). Inside the shock heating gap, the plasma density and field strength are reduced, whereas temperature is increased; the plasma-*β* is therefore enhanced in this location. Turning our attention to figure 7*b*, the inner band of shock heating is just outside (radially speaking) a similarly curved band of positive axial current, and the outer band of shock heating coincides with the inward pointing Lorentz vectors, which have been generated by negative axial current. Where the axial current is strongest (i.e. most negative), there are no significant Lorentz vectors, since in that region ** j**×

**≈0. At the start of the simulation, there is a positive current flowing at the centre of the midplane, surrounded by a ring of weaker negative current, in accordance with the initial**

*B**B*

_{θ}field. When the instability develops, the kinking motion compresses any plasma in its path and thereby enhances the magnitudes of the axial currents. Figure 7

*b*reveals that the outwardly radial direction described by the kinking begins just inside the region of positive current. Hence, some positive

*j*

_{z}is amplified, creating the outward pointing Lorentz vectors: however, the larger region of negative axial current produces a stronger

*j*

_{z}, which creates the opposing Lorentz forces. Of course, these changes in the axial current profile also alter the azimuthal field; see the red (solid) plots of figure 8. As the kink instability gets underway, magnetic energy is converted into kinetic; this amplifies the kinking, giving rise to strong helical currents, which create a locally enhanced azimuthal field. Magnetic flux is conserved and so moves with the plasma; hence an enhanced azimuthal field must correspond to a diminished axial field—figure 8 shows this to be the case.

### (c) Global evolution

The plots so far have all been taken at the midplane (*z*=0). For this reason, the final two figures are intended to give a sense of the global changes in plasma temperature and magnetic field. Figure 9 shows cross sections of the temperature (scaled according to the values given in §3) at three times during the simulation. We find that the situation seen at the midplane is also applicable to other *z*-coordinates: at *t*=95*t*_{A} (the start of the instability), slow mode shocks are creating high temperatures at two other sections along the loop, near the *x*=−0.5 line. The hot (red) areas are where the loop is kinking outwards against the outer potential field. Subsequently, the plasma becomes more turbulent, before starting a slow relaxation to a uniform temperature.

The field configuration undergoes similar changes to those seen in other simulations of kink unstable loops. At instability onset, the magnetic field starts to undergo a helical kink; the following turbulent phase results in a loss of azimuthal field (*B*_{x,y}), which leaves a predominantly axial field by the end of the simulation (figure 10). It is notable that the kinking is more radially confined compared to past work (see Bareford *et al.* [17], fig. 6 and Hood *et al.* [13], fig. 10). We attribute this to two aspects of the initial equilibrium, a smooth *α*-profile and zero net current, that force *B*_{θ} to have a maximum at *r*≈0.4. Although the cited work also features loops of zero net current, *α*(*r*) has a piecewise-constant profile, resulting in high axial currents near the loop boundary at *r*=1. Another important difference is in the way Ohmic heating is handled. Bareford *et al.* [17] found that regions of Ohmic heating appeared more or less at the same locations (and times) as shock heating, although the latter was far stronger than the former. By contrast, the use of a variable critical current employed here produces small pockets of Ohmic heating scattered near the footpoints and around the midplane. In figure 10, we use arrows to point out two of these pockets. Ohmic heating ceases completely after *t*=170*t*_{A}, and at the end of the simulation, shock heating is an order of magnitude lower than what it was at 120*t*_{A}.

## 5. Summary and conclusion

We have taken an initial field configuration, representing a straightened coronal loop, and induced an instability with the intention of investigating how the magnetic energy released is converted into heat. The ideal kink instability converts mainly azimuthal field (i.e. *B*_{x} and *B*_{y}) into Lorentz forces, which increases kinetic energy, thereby amplifying the kink. The resulting compression of plasma generates high (compared with the initial state) helical currents; this increases the Lorentz forces further, which feeds back to the velocity and ultimately generates heating through shock viscosity. In other words, sharp gradients in velocity (created by Lorentz forces) shear and compress local grid cells. The entropy production implied by these sudden changes in velocity is captured by shock viscosity. The subsequent heating causes a drop in velocity via the viscous term (∇⋅** σ**) in the equation of motion. This velocity reduction then leads to field dissipation when those same velocities are used to remap field components back to the Eulerian grid via the induction equation. It is at this stage that numerical dissipation occurs, as the use of van Leer gradient limiters during the remap exaggerates the reduction in magnetic field strength. This is confirmed by the energy plots (figure 2)—too much magnetic energy is dissipated compared with the level of shock heating. Interestingly, this numerical dissipation is much reduced when the coronal loop is only marginally unstable (i.e.

*λ*=1.62), then the dissipation is only 0.36% of the initial total energy for a spatial resolution of 256

^{2}×512.

It seems that employing a more realistic scheme as regards the application of anomalous resistivity (§3b) reduces severely the opportunities for Ohmic dissipation (figure 3). In the simulations presented here, it is the shock heating that causes the magnetic field to lose energy and relax to a lower energy state. In theory, if there were no form of heating, energy would forever move back and forth between kinetic and magnetic. Instead, the shock heating weakens axial current flow and reduces *B*_{x} and *B*_{y}; in other words, the kinetic energy causing the shock is not simply returned to the field.

Figures 5–7 demonstrate that the plasma dynamics causing the heating are consistent with Petschek reconnection and slow mode shocks. The advantage of this dissipation mechanism is that it can operate within a collisonless regime such as the corona. One key assumption made here is that the net result of photospheric convective motions *twists* the magnetic fields that reach into the corona. Over time, these fields would become susceptible to the kink instabitiy, which, as this work shows, is a sufficient condition for coronal heating via slow mode shocks. We expect these findings will apply to more realistic loop models; ones that feature thermal condition, atmospheric stratification and loop curvature [27,28,30]. Shock heating may also occur for other types of instability (e.g. torus or ballooning), especially if it causes a coronal loop to rapidly expand into an overlying field.

## Funding statement

We gratefully acknowledge the compute time (STFC funded) granted for the DiRAC and UKMHD facilities. The computer hardware was funded by a BIS National e-Infrastructure capital grant ST/K00042X/1, DiRAC Operations grant ST/K003267/1 and Durham University—DiRAC is part of the National e-Infrastructure.

## Acknowledgements

We thank the referees for their helpful comments. The simulations were run on the UK MHD cluster at St Andrews University and on the COSMA Data Centric system at Durham University. The latter facility is operated by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk).

## Footnotes

One contribution of 13 to a Theo Murphy meeting issue ‘New approaches in coronal heating’.

- Accepted February 19, 2015.

- © 2015 The Author(s) Published by the Royal Society. All rights reserved.