###### Abstract

It has been proposed that the million degree temperature of the corona is due to the combined effect of barely-detectable energy releases, so called nanoflares, that occur throughout the solar atmosphere. Alas, the nanoflare density and brightness implied by this hypothesis means that conclusive verification is beyond present observational abilities. Nevertheless, we investigate the plausibility of the nanoflare hypothesis by constructing a magnetohydrodynamic (MHD) model that can derive the energy of a nanoflare from the nature of an ideal kink instability. The set of energy-releasing instabilities is captured by an instability threshold for linear kink modes. Each point on the threshold is associated with a unique energy release and so we can predict a distribution of nanoflare energies. When the linear instability threshold is crossed, the instability enters a nonlinear phase as it is driven by current sheet reconnection. As the ensuing flare erupts and declines, the field transitions to a lower energy state, which is modelled by relaxation theory, i.e., helicity is conserved and the ratio of current to field becomes invariant within the loop. We apply the model so that all the loops within an ensemble achieve instability followed by energy-releasing relaxation. The result is a nanoflare energy distribution. Furthermore, we produce different distributions by varying the loop aspect ratio, the nature of the path to instability taken by each loop and also the level of radial expansion that may accompany loop relaxation. The heating rate obtained is just sufficient for coronal heating. In addition, we also show that kink instability cannot be associated with a critical magnetic twist value for every point along the instability threshold.

## Section 1 Introduction

The theory of nanoflare coronal heating (Parker, ?) postulates that small flares are sufficiently numerous to maintain a coronal temperature of millions of degrees. Observational studies have so far been unable to show conclusively whether nanoflares occur frequently enough to heat the corona (Krucker and Benz, ?; Parnell and Jupp, ?; Aschwanden and Parnell, ?; Parnell, ?): the smaller a flare the harder it is to distinguish from the coronal background. Nevertheless, we propose a coronal loop model, the purpose of which is to test the viability of the nanoflare theory with respect to coronal heating.

Random convective motions at the photosphere increase the energy of the magnetic fields that define coronal loops (plasma beta, 0.01). The magnetic field is repeatedly twisted by such motions; coronal conductivities are so large that the magnetic flux and plasma can be regarded as being frozen together. We assume that the kinetic energy imparted by the photosphere is dissipated via Direct Current heating: the timescale for photospheric turbulence is long compared to the Alfvén time (Klimchuk, ?). Hence, the loop can be represented as moving through a series of force-free states: , where / is the ratio of current density to magnetic field and is a position vector (Woltjer, ?). Magnetic stresses build within the loop until an instability is reached or dissipation is otherwise triggered.

Magnetic reconnection is thought to be the mechanism by which the excess magnetic energy is converted into heat. Observations of large-scale flares have revealed stong evidence for such a process (Fletcher, ?; Qiu, ?). In addition, several 3D MHD models have shown that current sheets form during the nonlinear phase of an ideal kink instability (Baty and Heyvaerts, ?; Velli, Lionello, and Einaudi, ?; Arber, Longbottom, and Van der Linden, ?; Baty, ?). The kink instability gives rise to helical current sheets that become the site of Ohmic dissipation, i.e., a heating event. These models are restricted to a narrow range of initial loop configurations. Expanding this range so that one could determine the relationship between say, the initial -profile and the resulting energy release would be too computationally expensive. However, the energy release can be calculated without recourse to following the complex dynamics of the reconnection process (Heyvaerts and Priest, ?).

When a magnetic field becomes unstable, relaxation theory predicts the field will transition to the lowest energy state that conserves total magnetic axial flux and global magnetic helicity (Taylor, ?, ?). The minimum energy (or relaxed) state is the well-known constant or linear force-free field, . The helicity (K) measures the degree to which the magnetic field is linked with itself (Berger, ?). However, the relative helicity (Berger and Field, ?; Finn and Antonsen, ?) is more useful since it is gauge invariant:

(\theequation) |

where is the magnetic potential, is the potential field with the same boundary conditions and is the corresponding vector potential.

Helicity conservation is not absolute. During relaxation, helicity is still subject to global resistive diffusion, but the change in helicity is negligible when compared to the drop in magnetic energy, so long as dissipation predominantly occurs within thin current sheets. The rates of dissipation for helicity and magnetic energy (W) are

(\theequation) |

where j = is the current density, l is the length scale of magnetic variation (i.e., current sheet thickness), L is the global length scale and is the resistivity (Browning, ?). Using and , the ratio of the dissipation rates reduces to . Hence, if , which is expected to be well satisfied for reconnecting current sheets within the highly conductive corona, where global resistive diffusion of helicity and energy are negligible. The relative sizes of the dissipation rates have been confirmed by MHD simulations, despite the coarseness of numerical grids (the difference between dissipation rates becomes more pronounced as the resistivity becomes smaller and falls below numerical precision). Browning et al. (?) showed that during the relaxation of a marginally (kink) unstable loop K/K 10 and W/W 10. Detailed estimates of coronal helicity dissipation are given by Berger (?); further justification for helicity-conserving relaxation is provided by laboratory experiments (Heidbrink and Dang, ?; Taylor, ?)

Helicity conservation combined with the invariant nature of the relaxed -profile imply that helicity has simply become more evenly distributed within the loop. Thus, the relaxed can be calculated, as can the amount of energy liberated from the magnetic field during relaxation. The latter quantity can be interpreted as an upper limit to the heating event energy, since complete relaxation may not be attained. Issues concerning plasma response are outside the scope of the model presented in this paper.

Repeated episodes of slow photospheric driving may trigger an ideal MHD instability. Ideal (not resistive) instabilities are required in order for the time scale of the instability to match observations of the highly conducting corona, where resistive time scales are very long. Browning and Van der Linden (?) describe how a dynamic heating event is caused whenever the field exceeds the threshold for a linear kink instability in a cylinderical coronal loop. Extensive numerical simulations (Galsgaard and Nordlund, ?; Velli, Lionello, and Einaudi, ?; Lionello et al., ?; Baty, ?; Gerrard et al., ?) have demonstrated how energy release occurs during the nonlinear phase of the instability. Relaxation of the field towards a constant- state has also been demonstrated (Browning et al., ?; Hood, Browning, and Van der Linden, ?) alongwith a measurement of the energy release that is in good agreement with relaxation theory.

In summary, the model presented here is based on the idea that coronal loops are moved by photospheric motions through a series of force-free equilibria that will eventually culminate in a heating event whenever the ideal instability threshold is crossed. Bareford, Browning and Van der Linden (?, hereafter BBV?) calculated the heating event distributions for an ensemble of loops that possessed net current. This was done using a simple cylinderical field model in which the current profile ((r)) of the stressed field is represented by a two parameter family. (Theoretically, only two states on the closed instability threshold for two-parameter loops could yield zero-net-current configurations.) The resulting distributions were compatible with the energies required for coronal heating. Following BBV?, we use a simple cylindrical field model, and represent the nonlinear force-free field using a three-parameter family of current profiles in which is a piecewise constant. We now attempt to improve the realism of the model by including a mechanism that requires each loop to have zero net current. Thus, we consider the effects of twisting motions localised within the loop cross section (Hood, Browning, and Van der Linden, ?). Previously, the random nature of photospheric motions was represented by allowing different parts of the loop interior to vary independently. It is more reasonable however to assume some level of correlation, since it is likely that the whole of a loop footpoint will be subjected to the same convective eddy. Furthermore, we also consider the consequence of varying the aspect ratio of the loop (L /R). This paper will also check to see if there is a twist-based parameter that can be used as a simple diagnostic for loop instability, e.g., Hood and Priest (?).

The composition of the loop model and the equations used to express the magnetic field are given in the next section. Section 3 describes how the loop instability threshold is calculated and presents the equations used to determine the energy release associated with each relaxation. This section also analyses the loop configurations that follow the instability threshold. In addition, the loop’s path to instability is explained. The nanoflare energy distributions produced by the model are presented in Section 4. Finally, in the last section, the results are summarised and discussed.

## Section 2 Model

A loop is considered to evolve through equilibria as it is driven by slow photospheric footpoint motions. An idealised model of a straight cylindrical loop is used with the photosphere represented by two planes at z = 0, L; however, the essential physics should apply to more complex geometries. The stressed field is line-tied with (in general) a non-uniform (r) (where /). This is represented by a three-parameter family of piecewise-constant- profiles.

The model discussed here (Figure 1) is an improved version of the one used by BBV?, which allowed loops to have net current (i.e., an azimuthal field was present in the potential envelope). A net current contradicts the idea that the magnetic field is twisted by convective motions local to the loop footpoints. If the twisting motions are confined to some localised region, the field surrounding the loop should be purely axial. The currents generated by the twisting of the fields within the loop should close locally such that the loop carries zero net current. Hence, a current neutralisation layer is introduced here, defined such that the azimuthal field () always falls to zero at the loop boundary (R). Note that Hood, Browning and Van der Linden (?) undertook 3D numerical simulations with initial fields taken as twisted states with zero net current. There are, however, only two such fields on the marginal instability curve of BBV?, for which the current due to cancels that due to (that is 2.48, 0.95). Here, we construct a whole family of current-neutralised fields by adding the extra layer.

Following on from BBV? and Browning and Van der Linden (?), it is proposed that a relaxation event is triggered when the loop’s field becomes linearly unstable. The energy released, due to fast magnetic reconnection during the nonlinear development of the instability, can then be calculated using relaxation theory.

The loop’s radial -profile is approximated by a piecewise-constant function featuring three parameters. The ratio of current to magnetic field is in the core, in the outer layer, in the neutralisation layer and zero in the potential envelope. The free parameters are and , whereas is dependent on the first two and is determined by the requirement of zero net current. Note that the magnetic field is continuous everywhere, whereas the current has discontinuities. Recent work indicates that a discontinuous current has little discernable effect on the dynamics when compared to similar but continuous -profiles (Hood, Browning, and Van der Linden, ?). Following the investigations of Browning et al. (?), the outer surface of the potential envelope, representing the background corona, is placed at R = 3 (thrice the loop radius).

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

(\theequation) | |||||

(\theequation) | |||||

(\theequation) | |||||

(\theequation) | |||||

(\theequation) |

where (i = 1,2,3) represent the sign of . The fields must be continuous at the inner radial boundaries, R, R and R. (We choose R = 0.5, R = 0.9 and R = 1, so that most of the loop is similar to previous work, BBV?, with a thin current neutralisation layer between R and R) Therefore, the constants B and C (j = 2,3,4) can be expressed like so:

(\theequation) |

where

(\theequation) |

The value of (the neutralisation layer current) is found, for a given (, ), by numerical solution of B(R) = 0, ensuring that the net current is zero and that the azimuthal field vanishes outside the loop, see Equation 2. Thus, the equilibrium parameter space remains 2D (i.e., it is determined by , ) - the field profiles for a selection of loop configurations are given in Appendix B.

The magnetic flux through the loop and envelope is conserved:

(\theequation) | |||||

where the asterisks denote dimensionless quantities. Hence, in the model, is normalised to 1 and B is determined (noting that, in Equations 2-\theequation, B and C are functions of ). The normalised coronal loop radius (R = 1) is itself used to normalise the loop length (e.g. L = 20R), see Figure 1.

We assume that the loop evolves through a sequence of two-alpha fields (Equations \theequation-\theequation) as it is twisted by photospheric footpoint motions. It has been shown that the introduction of magnetic twist gives the coronal loop a circular cross section (Klimchuk, Antiochos, and Norton, ?). The loop model presented here has the same cross-sectional shape, but the loop radius (R) is held constant throughout the simulated photospheric driving. Purely azimuthal photospheric motions would cause a small expansion of the loop (Browning and Hood, ?) which we ignore; alternatively, small radial footpoint motions must be allowed in order to maintain constant loop radius. In any case, the sequence of loop equilibria explored by our model is clearly a small subset of the possible variation in field profiles that might arise from photospheric motions. As these random motions proceed, the loop continues to evolve through force-free equilibria until it becomes linearly unstable. We now discuss the calculation of the instability onset.

## Section 3 Linear kink instability thresholds

A coronal loop’s instability is constrained by the line-tying of the photospheric footpoints (Hood, ?). Hence, all perturbations are required to vanish at the loop ends (z = 0, L). Any linear perturbation can be decomposed as a sum, . We need only consider the m = 1 term however, since this azimuthal mode has been found to be the least stable (Van der Linden and Hood, ?) and is the dominant instability. The effect of such perturbations on the coronal loop are represented by the standard set of linearised ideal MHD equations. When the growth rate of a perturbation transitions from a negative value to a positive one, the loop has reached the threshold of an ideal linear instability. The instability threshold is a curve in 2-dimensional -space (, ). The properties of the loop (e.g., , and ) at these threshold points can be found by substituting the perturbation function into the linearised MHD equations, leading to an eigenvalue equation for the growth rates (Priest, ?). The growth rates and eigenfunctions of the most unstable modes are found numerically, for line-tied fields, with the CILTS code, described in Browning and Van der Linden (?) and Browning et al. (?). CILTS can be configured such that one of the loop’s free parameters is fixed whilst the other is incremented. The code terminates as soon as the real part of the eigenfunction falls below zero, i.e., the loop is no longer unstable to kink perturbations.

Figure 2 shows the instability threshold curves mapped by the CILTS code for three values of loop aspect ratio (L /R = 10, 20, 30). The longer the loop the smaller the -value required for instability, since, if the radius is held constant, longer loops are less affected by line-tying stabilisations (Priest and Hood, ?). The addition of a current neutralisation layer prevents the threshold curves from closing near the -space origin; this is unlike the net current case for which the threshold is a closed curve (BBV?, Figure 2). The open shape is indeed similar to the instability curve for loops with a conducting wall at R (Browning and Van der Linden, ?). This is because the eigenfunctions almost vanish at R, see Figures 4-6. Also, if is small, will be opposite in sign to , and the outer layer is stabilised by the neutralisation layer. However, the loop configurations become unrealistic as we increase the magnitude of and the axial field reverses. Positions outside the B reversal lines (Figure 2) represent loops that have axial fields of mixed polarity. These configurations cannot represent states attained by the twisting of a unipolar loop.

Before proceeding to calculate the energy release properties, it is first of interest to analyse the marginally unstable states. The threshold curves shown in Figure 2 have symmetry: a rotation of radians leaves the thresholds unchanged. Thus, it is sufficient to show how various properties (e.g., magnetic twist and energy release) vary along the parts of the threshold where 0 (Figure 3). For ease of plotting we can convert the threshold curves to a one dimensional form: the filled circles and associated numbers of Figure 3 (left) represent the tic marks and labels for the 1D threshold point axis, see Figures 7-11. Such figures will always be calculated using a loop of aspect ratio 20.

Firstly, we plot the linear eigenfunctions for a selection of marginally unstable -space points that follow the instability threshold. The location of these points can be determined from the threshold point number given at the top of each plot in Figures 4-6. The eigenfunctions of the marginally unstable modes are stongly radially confined; that is, there is almost no disturbance beyond the loop radius (R = 1). This contrasts sharply with the situation for loops with net current (Browning et al., ?; BBV?), in which the eigenfunction generally extends well into the potential envelope.

### 3.1 Instability threshold and critical twist

Following BBV?, we look for a single twist-related parameter that takes on a critical value whenever the loop reaches the threshold (Appendix B shows the twist profiles for a selection of loop configurations, stable and unstable). It has often been postulated that instability can be identified with a single critical twist value irrespective of the detailed field profiles. The average twist can be calculated in several ways;

(\theequation) |

Equation \theequation is the average twist weighted by area, while 3.1 and 3.1 have been used by Velli, Einaudi, and Hood (?) and Baty (?). Note, Equation 3.1 can be calculated analytically, see Appendix A. denotes the average twist, weighted by area, over the core, outer layer and current neutralisation layer. The tilde () and hat () symbols are used to indicate the other equations.

None of the twist averages (Figure 7) is invariant around the whole threshold curve, although 5 (Equation 3.1) for the majority of threshold points. This value is approximately in line with the oft-quoted result of 2.49, the critical twist for a loop of aspect ratio 10 (Hood and Priest, ?). Each threshold point has a radial twist profile; these profiles feature reversed twist until around point 60, where the profile becomes single signed. After this point, the three average-twist plots converge to values between 5 and 10. At higher threshold points, the plots diverge and for Equation 3.1 and \theequation the averages increase sharply.

Finally, we consider the proposal of Malanushenko et al. (?), that a critical value of normalised helicity (equivalent, in our terms, to the normalised loop helicity, , over the range 0 - R) indicates instability onset.

Figure 8 shows that the normalised helicity is certainly not the same for every threshold point, even if and have the same sign.

### 3.2 Path to instability

BBV? used a random walk process to simulate a loop being twisted by turbulent photospheric motions. In other words, a loop performed a sequence of fixed-length steps of random direction within -space until the instability threshold was crossed. We will follow this process for zero-net-current loops too, however, we will also employ spatially correlated random walks. This is to allow the correlation between the inner and outer parts of the loop to be varied. In particular, it is more likely that the twisting will be fairly uniform across the loop (i.e., the change in is similar to the change in ).

When a loop begins its random walk (i.e., when it emerges from beneath the photosphere) it is assigned a random starting position within the stable region of -space equilibria (i.e., the loop may have some initial twist). It is more likely however, that the initial twist will be small and that the initial value of is similar to (or correlated with) the initial value. Furthermore, the change in -coordinates that occur whenever the loop steps through -space, in response to photospheric driving, should also be correlated. The initial coordinate of the walk is chosen from a normal distribution centred on zero. A standard deviation is chosen such that the probability of the initial value representing an unstable configuration is negligible. Similarly, the initial coordinate is chosen such that the mean is the initial coordinate.

The step values, and , are determined by assuming a step size, , and . Hence, is also chosen from a normal distribution, but this time the mean is and is chosen such that the mean is .

As the standard deviation of the normal distribution used to select and is decreased, the range of threshold crossings narrows. In other words, the walks follow the = line more closely. A standard deviation of 0.1 will be used for the simulations presented in this paper, since this value restricts the threshold crossings to points where the twist is single-signed, see Figure 9. Correlated walks therefore are predisposed to maintaining the realism of loop configurations, since it is expected that in general photospheric motions do not create loops that have reversed twist.

Figure 9 shows that a loop might cross a B reversal line before it reaches the instability threshold. If this happens, the loop is discarded and the simulation resumes with a new loop that has a stable -configuration. Once a loop reaches the instability threshold, it becomes linearly unstable. At this point, the field releases energy and transitions to a lower-energy state defined by Taylor relaxation: helicity is conserved and the -profile relaxes to a single value.

### 3.3 Energy release calculation

We allow each loop of an ensemble of 10 loops to undergo a single relaxation (BBV?). Initially, a loop starts from an assigned stable state. The field profile then undergoes a random walk (which may or may not be correlated) until it crosses the instability threshold; whereupon, the loop relaxes and the profile transitions to the relaxation line ( = ). The relaxation () will, of course, vary depending on where the threshold was crossed; is found by helicity conservation (Taylor, ?; Heyvaerts and Priest, ?; Taylor, ?; Browning and Van der Linden, ?). In mathematical terms, we find the roots of the following equation:

(\theequation) |

where and are the coordinates of the instability threshold crossing (conservation of axial flux is assured through the normalisation = 1). The helicity can be expressed as follows:

(\theequation) |

where I(r) is the axial current and L is the loop length (Finn and Antonsen, ?). Axial flux is also conserved. The full expressions for helicity and energy are given in Appendix A. The energy difference between the unstable and relaxed states can be calculated thus:

(\theequation) |

This is the relaxation energy: the energy that is liberated from the magnetic field during the event. How much of this energy is converted to heat depends on the plasma response; thus, the relaxation energy represents the upper limit of the energy available for plasma heating.

In BBV?, a loop with net current was relaxed such that the -profile became invariant over the range 0-R. Hence, the relaxed state always represented a threefold radial expansion of the threshold state (i.e., from R to R), the relaxation encompassed both the loop and the potential envelope. Numerical simulations (Browning et al., ?) indicate that this is a good model for loops with net current. However, for loops with zero net current, the instability is more radially confined and the reconnection activity is correspondingly localised; it does not extend to the outer boundary (Hood, Browning, and Van der Linden, ?; Bareford et al., ?).

We therefore consider that the relaxation radius, R, can be anywhere in the range R(= 1) R R(= 3). If R = R, we have complete relaxation as previously considered (Browning and Van der Linden, ?; Browning et al., ?; BBV?); otherwise relaxation is localised. is constant between 0 and R and the fields in the remaining envelope (where = 0 and R r R) are fixed so that they do not change during relaxation; this maintains current neutralisation, albeit via an infinitely thin current-neutralising surface. Axial flux is conserved, such that (over 0-R) of the threshold state is equal to (over 0-R) of the relaxed state. K/ is conserved in an identical manner (in our previous work, conservation was always over 0-R and since the total axial flux was normalised to 1, conserving K/ was identical to conserving K). Likewise, the energy release is the energy of the threshold state over 0-R minus the energy of the relaxed state over the same radial range. In fact, the energy of the remaining potential envelope is unchanged, so that the energy release could also be taken over the entire volume (0-R); similarly, the envelope has zero helicity before and after relaxation.

## Section 4 Distribution of energies and coronal heating considerations

We now proceed to the main task of our work, which is to calculate the distribution of magnitudes of the sequence of heating events generated by random photospheric driving. First, we show how various properties vary along the instability threshold.

### 4.1 Helicity and energy

The left panel of Figure 10 plots total helicities of the threshold states. A total helicity (or flux) is one calculated over the range 0 - R, i.e., the loop and envelope.

None of the threshold states have sufficient helicity for the relaxed state to feature helical modes (Taylor, ?) and so all relaxed states are cylinderically symmetric.

The left panel of Figure 11 confirms that each threshold state corresponds to a relaxed state. Both of the graphs in Figure 11 feature two plots; the dashed line represents minimum relaxation (R = R) and the solid line represents full relaxation (R = R). The right panel shows that, in general, W is affected by R (although, there is one part of the threshold where the energy release is insensitive to relaxation radius); hence, the energy distributions that appear in the next section are calculated for minimum and maximum relaxation radii.

The energies shown in Figure 10 (right) and Figure 11 (right) are given as dimensionless quantities; BBV? derived the following expression for calculating a dimensional energy,

(\theequation) |

where R is the loop radius in the corona and B is the mean axial field in the corona. Assuming typical values (R = 1 and B = 0.01 ), we obtain dimensional energy values of . Thus, the top end of the W scale ( 0.05) for R = 3 is equivalent to . This is in the microflare range, but nanoflare energies will be obtained for weaker fields.

### 4.2 Flare energy distribution

An expression for the energy flux is derived by considering the loops in the ensemble as spatially separated but flaring simultaneously. All the energy input from the photosphere is dissipated, in a long-term time-average over many events, since the instability threshold limits the accumulation of stresses within the coronal magnetic field. The energy flux, F, is thus;

(\theequation) |

where N is the average number of steps taken to reach the threshold and is the time taken to complete each step in the random walk. Similar expressions exist in the literature based on random photospheric twisting (Sturrock and Uchida, ?; Berger, ?; Zirker and Cleveland, ?; Abramenko, Pevstov, and Romano, ?).

In other words, is the time taken for to change by /R; we may estimate, a timescale for this process as follows. Based on axial values, corresponds to a change in magnetic twist (L/2)( /R); taking L /R = 20 and = 1 gives 10. If this is caused by photospheric twisting motions of magnitude v for a time interval , we find ()R/v, where R is the footpoint radius (at the photosphere). With typical values of R = 200 and v = 1 , we obtain 2000 ; note that this is consistent with the expected correlation time for photospheric motions; granule lifetimes are of the order 10 (Zirker and Cleveland, ?). Hence has a linear relationship with the loop length, = 100(L []). Applying the previously used values for B and R, gives a dimensional flux of . This result is applicable to Active Regions (a value for the Quiet Sun can be obtained by setting B = 0.001 ; this simply lowers the multiplier (10) by 2 orders of magnitude.

The energy distributions given below are each derived from a loop ensemble, that is a collection of 10 loops flaring simultaneously. Since each loop relaxes only once we can sidestep the complications that come with allowing loops to survive many relaxations: for example, a loop may shrink or implode after flaring (Janse and Low, ?), and a different instability threshold would have to be applied should the aspect ratio be altered as a consequence.

#### 4.2.1 Distribution of ”nanoflares”

Examination of Figure 12 yields three key points. Firstly, the total energy released increases with aspect ratio, but the average step count, N, decreases. This is to be expected since loop volume increases with aspect ratio, whereas the size of the stability region shrinks, see Figure 2. Secondly, as indicated before (Figure 11, right), increasing the relaxation radius increases the energy released. And thirdly, correlated walks mean higher step counts, however, whether or not there is also an increased energy release depends on the relaxation radius.

If R = R(= 1), the energy release from correlated walks is reduced compared to the uncorrelated distributions; whereas, complete relaxation, R = R(= 3), leads to an increased energy release. This less-than-straightforward point is consistent with the plot that shows the variation in energy release along the threshold for both values of R, see Figure 11 (right). A correlated walk would favour crossings around threshold point 90; when R = R the energy release is almost at its lowest for this part of the threshold, whereas the opposite is the case when R = R. This is also true for the thresholds applicable to loops of aspect ratio 10 and 30.

For loops of aspect ratio 10, correlated walks produce distributions that have high-energy cut-offs - this feature is an artifact of the simple two- model. It is caused by the fact that when L /R = 10, the relaxation line intersects the B reversal line before the instability threshold.

When one calculates the dimensional heat fluxes (Equation \theequation) one finds that flux is weakly dependent on aspect ratio. Further examination reveals that any dependence on aspect ratio can only come from , which is determined by the coordinates of the instability threshold. incorporates a length factor in units of the loop radius, i.e., , where is the dimensionless energy release per unit of dimensionless length. Substituting the full expression for the step time (, Section 4.2) into Equation \theequation gives

(\theequation) |

again, applying previously used values, this simplifies to erg cm s. The length terms cancel and the ratio is effectively a constant.

For distributions derived from uncorrelated walks and minimal relaxation (R = R), F 3-410 . Using correlated walks instead, diminishes the fluxes to 0.9-210 . Increasing the relaxation radius to R will reverse this reduction and yield F 7-1010 . This last result is also true for distributions based on uncorrelated walks and full relaxations. When R = R correlated walks do lead to higher energy releases, however, these walks are longer and have higher step counts, which means the flux remains roughly constant.

Finally, in Figure 13 we show the logs of the flare energy distributions presented in Figure 12. The distributions calculated from uncorrelated walks give log plots that almost match the critical gradient for coronal heating. Although the log plots of the correlated (Gaussian-shaped) distributions do not follow power-laws, we include these results for completeness.

## Section 5 Summary and conclusions

We have investigated the distribution of energy releases in an ensemble of coronal loops driven by random photospheric footpoint motions, using Taylor relaxation theory. The twisting has been assumed to be localised within the loop cross section, so that the loop is always without net current (the azimuthal field vanishes at - and beyond - the loop boundary). This means we are genuinely studying individual loops, rather than (unrealistically) allowing the potential envelope outside the loop to be twisted as the loop evolves, as in previous work (BBV?).

A relaxation event is triggered whenever the loop becomes unstable to an ideal kink instability; during the nonlinear phase, current sheets form and subsequent rapid reconnection occurs. A distribution of events is built up by allowing loop equilibria to evolve through a random walk, representing the effects of turbulent photospheric footpoint motions, until the linear stability threshold is reached. The effectiveness of the relaxation approach is that energy release is easily calculated for a wide family of profiles, this is extremely difficult with 3D numerical simulations. Furthermore, relaxation theory becomes a better representation for very high conductivities, which cannot be accessed by present day simulations. This approach can, of course, be extended to more complex field models than the simple cylindrical coronal loop models used here. The energy fluxes obtained are sufficient for coronal heating; the fluxes agree with the oft-quoted estimates of Withbroe and Noyes (?). The heating flux is larger for photospheric motions with higher temporal correlation: within our model, this is represented by larger step sizes for the random walk (so that the loop is coherently twisted until it reaches instability, rather than randomly twisting and untwisting many times). We thus show that dissipation within loops could be sufficient for coronal heating, but in reality, this is likely to form only part of the coronal heating input. Topological complexity arising from, for example, the discrete nature of photospheric flux sources (Priest, Heyvaerts, and Title, ?) and braiding motions (Parker, ?), will also play a strong role.

A distribution of heating events, or nanoflares, is obtained, for a variety of conditions. For the case of spatially uncorrelated twisting motions, in which the motions may vary strongly across the loop cross section, a power law distribution of energy versus occurrence frequency is obtained, with a slope slightly steeper than the critical value of -2 required for nanoflare heating to be effective (Hudson, ?). For strongly correlated twist motions, in which the twist in the outer part of the loop is close to that in the inner core, a peaked energy distribution is obtained, with almost Gaussian shape. The former case reflects the distribution of available energies around the instability threshold, whereas as the latter is mainly determined by the allowable range of twist profiles. It should be noted that these distributions (Figures 12-13) are obtained for an ensemble of identical loops: in reality, much broader distributions will result due to variations in axial field strengths and photospheric driving. The true nanoflare distribution is a convolution over more than one parameter.

The effect of loop aspect ratio has been considered and been found to have little impact on energy flux. The higher volume of large aspect ratio loops is counteracted by the smaller stability region (instability occurs at lower values). As the aspect ratio is increased beyond 30, we expect the stability region to reduce by smaller and smaller amounts. In other words, the region will converge to a minimum area. This has been shown for constant- loops, see Figure 4 of Browning and Van der Linden (?). Hence, assuming this expectation is verified, the energy flux will be independent of aspect ratio (Equation \theequation), assuming that the same axial field strength is applied to all members of the loop ensemble. Presumably, there is a dependence between loop size and , so an ensemble that features some distribution of field strengths will still depend (albeit indirectly) on the aspect ratio.

Contrary to the B profiles of Appendix B, the axial field at the loop footpoints should not change during the random walk or during relaxation. The reason for this discrepancy is that preservation of the footpoint axial field introduces a z dependency - the field becomes two dimensional. However, if the length of the loop exceeds its radius, a 1D field approximation - such as the one used by the model presented here - still remains adequate for a substantial portion of the loop. Zweibel and Boozer (?) and Browning and Hood (?) show that the z dependence is confined to thin boundary layers near the footpoints. Hence, the difference in energies for loops represented by 1D and 2D fields is negligible especially if L /R 10 (see also Lothian and Browning, ?; Robertson, Hood, and Lothian, ?). Dalmasse, Browning and Bareford (?) have investigated a simpler loop, having just a core and outer layer (with a conducting wall at r = 1), by calculating the energy releases according to Taylor relaxation for a representative sample of threshold configurations. This was done using both 1D and 2D fields, with the latter maintaining the axial field at the footpoints. The resulting energy releases differ by less than 1% between the 1D and 2D cases.

The results presented here are based on a loop model that has a thin current-neutralising layer (this approximates to a current sheet), in which the fields discontinuously change at the loop edge. The main reason for this choice is so that the fields inside the loop are close to the previously-studied two- model (BBV?), and thus a comparison can be made with previous work. Also, such fields correspond to twisting within an isolated flux source, whilst the flux which surrounds the loop in the corona originates from untwisted separated sources. Interestingly, the ideal instability threshold in this case is very similar to that found for a close-fitting conducting wall at the loop edge, as originally used by Browning and Van der Linden (?). This is because the thin current layer forces unstable perturbations to vanish (almost) at the loop edge. In numerical simulations, the choice of a thin current layer has consequences in allowing resistive modes to be significant; although for realistic values of the resistivity (unattainable in simulations) the growth rate of such modes is extremely slow. Preliminary studies have also been undertaken with a thicker current-neutralising layer. In this case, a closed stability threshold curve can be obtained, and the results are more similar to previous work (BBV?).

One important consequence of considering loops with zero net current is that the reconnection activity tends to be localised near the loop and thus relaxation is likely to be incomplete (rather than including a large part of the surrounding potential field). We consider here two limiting cases: localised relaxation, in which only the loop volume relaxes to a minimum energy (constant-) state, and the surrounding potential envelope remains unaffected; and complete relaxation, in which the loop and the potential envelope relax out to the external boundary. The latter is clearly the true minimum energy state. Numerical simulations (Bareford et al., ?) indicate an intermediate situation, but somewhat closer to the completely localised relaxation. In fact, the loop reconnects with some of the surrounding axial field, but only to a limited extent. This is an important issue for understanding relaxation in the Sun, where the extent of relaxation is not defined by conducting walls: in contrast with laboratory plasmas (Taylor, ?). This is discussed more fully in a companion paper (Bareford et al., ?). In general, complete relaxation naturally gives larger energy releases, but the choice of relaxation radius does not strongly affect the distribution of heating events. Future work will use numerical simulations to explore the transition to instability, and the effects of continual driving.

Acknowledgements. The authors thank Jim Klimchuk for constructive comments and M. R. B. thanks STFC for studentship support.

## Appendix

## Section A Expressions for loop properties

Expressions for some key quantities (, K and W) are given here. For compactness, these are given only for 0 and 0, while special cases (e.g., = 0) must be dealt with separately. Expressions for constant- fields can be recovered by setting = , which gives more familiar formulae. The superscripts and subscripts that accompany each quantity term denote the upper and lower radial bounds over which the quantity is calculated. The vacuum permeability, , used in the magnetic energy expressions, is set to 1.

### a.1 Average magnetic twist

(\theequation) |

### a.2 Magnetic helicity

(\theequation) |

### a.3 Magnetic energy

(\theequation) |

## Section B Magnetic field profiles for a selection of -space points

The magnetic axial field, B, azimuthal field, B, and magnetic twist, (r), profiles are presented for a selection of stable and unstable loop configurations, see Fig. 14.

The empty circle on the = line in Fig. 14 is the relaxed state of the unstable configuration identified by the filled circle on the threshold. The radius of the relaxed state is 1.5.

## References

- 2006 Abramenko, V. I., Pevstov, A. A., and Romano, P.: 2006, Astrophys. J. Lett. L81, 646.
- 1999 Arber, T. D., Longbottom, A. W., and Van der Linden, R. A. M.: 1999, Astrophys. J. 517, 990.
- 2002 Aschwanden, M. J. and Parnell, C. E.: 2002, Astrophys. J. 572, 1048.
- 2010 Bareford, M. R., Browning, P. K., and Van der Linden, R. A. M.: 2010, Astron. Astrophys. 521 A70.
- 2011 Bareford, M. R., Browning, P. K., and Hood, A. W.: 2011, in preparation.
- 1996 Baty, H. and Heyvaerts, J.: 1996, Astron. Astrophys. 308, 935.
- 2000 Baty, H.: 2000, Astron. Astrophys. 360, 345.
- 2001 Baty, H.: 2001, Astron. Astrophys. 367, 321.
- 1984 Berger, M. A.: 1984, Geophys. Astrophys. Fluid Dynamics 30, 79.
- 1991 Berger, M. A.: 1991, Astron. Astrophys. 252, 369.
- 1999 Berger, M. A.: 1999, Plasma Phys. Control. Fusion 41, B167.
- 1984 Berger, M. A. and Field, G. B.: 1984, J. Fluid. Mech. 147, 133.
- 1986 Browning, P. K., Sakurai, T., and Priest, E. R.: 1986, Astron. Astrophys. 158, 217.
- 1988 Browning, P. K.: 1988, J. Plasma Phys., 40, 263.
- 1989 Browning, P. K. and Hood, A. W.: 1989, Solar Phys. 124, 271.
- 2003 Browning, P. K. and Van der Linden, R. A. M.: 2003, Astron. Astrophys. 400, 355.
- 2008 Browning, P. K., Gerrard, C., Hood, A. W., Kevis, R., and Van der Linden, R. A. M.: 2008, Astron. Astrophys. 485, 837.
- 2011 Dalmasse, K., Browning, P. K., and Bareford, M. R.: 2011, in preparation.
- 1985 Finn, J. M. and Antonsen, T. M.: 1985, Communications in Plasma Physics and Controlled Fusion 9, 111.
- 2009 Fletcher, L.: 2009, Astron. Astrophys. 493, 241.
- 1997 Galsgaard, K. and Nordlund, Å.: 1997, J. Geophys. Res. 102(A1), 219.
- 2001 Gerrard, C. L., Arber, T. D., Hood, A. W., and Van der Linden, R. A. M.: 2001, Astron. Astrophys. 373, 1089.
- 2000 Heidbrink, W. W. and Dang, T. H.: 2000, Plasma Phys. Control. Fusion 42, L31.
- 1984 Heyvaerts, J. and Priest, E. R.: 1984, Astron. Astrophys. 137, 63.
- 1992 Hood, A. W.: 1992, Plasma Phys. Control. Fusion 34, 411.
- 2009 Hood, A. W., Browning, P. K., Van der Linden, R. A. M.: 2009, Astron. Astrophys. 506, 913.
- 1979 Hood, A. W. and Priest, E. R.: 1979, Solar Phys. 64, 303.
- 1981 Hood, A. W. and Priest, E. R.: 1981, Geophys. Astrophys. Fluid Dynamics 17, 297.
- 1991 Hudson, H. S.: 1991, Solar Phys. 133, 357.
- 2007 Janse, Å. M. and Low, B. C.: 2007, Astron. Astrophys. 472, 957.
- 2000 Klimchuk, J. A., Antichos, S. K., and Norton, D.: 2000, Astrophys. J. 542, 504.
- 2006 Klimchuk, J. A.: 2006, Solar Phys. 234, 41.
- 1998 Krucker, S. and Benz, A. O.: 1998, Astrophys. J. 501, L213.
- 1998 Lionello, R., Velli, M., Einaudi, G., and Mikić, Z.: 1998, Astrophys. J. 494, 840.
- 2000 Lothian, R. M. and Browning, P. K.: 2000, Solar Phys. 194, 205.
- 2009 Malanushenko, A., Longcope, D. W., Fan, Y., and Gibson, S. E.: 2009, Astrophys. J. 702, 580.
- 1994 Melrose, D. B., Nicholls, J., and Broderick, N. G.: 1994, J. Plasma Phys. 51, 163.
- 1972 Parker, E. N.: 1972, Astrophys. J. 174, 499.
- 1988 Parker, E. N.: 1988, Astrophys. J. 330, 474.
- 2004 Parnell, C.E.: 2004, in: R. W. Walsh, J. Ireland, D. Danesy, and B. Fleck (eds.), Proc. SOHO 15 Workshop - Coronal Heating, ESA SP-575, 227.
- 2000 Parnell, C. E. and Jupp, P. E.: 2000, Astrophys. J. 529, 554.
- 1987 Priest, E. R.: 1987, Solar Magnetohydrodynamics, D. Reidel Publishing Company, Dordrecht (The Netherlands), Chap. 7.
- 1979 Priest, E. R. and Hood, A. W.: 1979, Solar Phys. 64, 303.
- 2002 Priest, E. R., Heyvaerts, J. F., and Title, A. M.: 2002, Astrophys. J. 576, 533.
- 2005 Priest, E. R., Longcope, D. W., and Heyvaerts, J.: 2005, Astrophys. J. 624, 1057.
- 2009 Qiu, J.: 2009, Astrophys. J. 692, 1110.
- 1992 Robertson, J. A., Hood, A. W., and Lothian, R. M.: 1992, Solar Phys. 137, 273.
- 1981 Sturrock, P. A. and Uchida, Y.: 1981, Astrophys. J. 246, 331.
- 1974 Taylor, J. B.: 1974, Phys. Rev. Lett. 33, 1139.
- 1986 Taylor, J. B.: 1986, Rev. Mod. Phys. 58, 741.
- 1999 Van der Linden, R. A. M., and Hood, A. W.: 1999, Astron. Astrophys. 346, 303.
- 1990 Velli, M., Einaudi, G., and Hood, A. W.: 1990, Astrophys. J. 350, 428.
- 1997 Velli, M., Lionello, R., and Einaudi, G.: 1997, Solar Phys. 172, 257.
- 1977 Withbroe, G. L. and Noyes, R. W.: 1977, Ann. Rev. Astron. Astrophys. 15, 363.
- 1958 Woltjer, L.: 1958, Proc. Natl. Acad. Sci. USA, 44, 489.
- 1993 Zirker, J. B. and Cleveland, F. M.: 1993, Solar Phys. 144, 341.
- 1985 Zweibel, E. G. and Boozer, A. H.: 1985, Astrophys J. 295, 642.