EDP Sciences
Free Access
Issue
A&A
Volume 550, February 2013
Article Number A40
Number of page(s) 14
Section The Sun
DOI https://doi.org/10.1051/0004-6361/201219725
Published online 23 January 2013

© 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 re-enter the solar surface at another. The convective motions at these photospheric boundaries (or footpoints) are thought to reconfigure coronal fields and thereby cause free magnetic energy to accumulate, making the fields more susceptible to instability. Essentially, the kinetic energy of the convection zone is transported, via a Poynting flux, through the photosphere and stored within the coronal loop. In the case of a single loop, currents are created by those convective motions that cause the footpoints to be twisted. The approximate time scale for the twisting motions is long compared to the Alfvén time, and so the coronal loop (not necessarily illuminated) transitions through a series of force-free equilibria, which can be expressed as ∇ × B = α(r)B; where r is a position vector, and α = (μ0j·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 loop-like magnetic field, driven by continual convective motions, reaches the threshold for kink instability (Hood 1992; Browning & Van der Linden 2003; Haynes & Arber 2007; Srivastava et al. 2010). Coronal magnetic fields cover many thousands of kilometers (L≈ 50   Mm) and exist within a highly conductive environment: therefore, in the absence of instability, magnetic fields diffuse slowly (td ≈ 103  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 three-dimensional (3D) magnetohydrodynamic (MHD) models (Baty & Heyvaerts 1996; Velli et al. 1997; Arber et al. 1999; Baty 2000; Browning et al. 2008; Hood et al. 2009).

The coronal heating problem is most pronounced within active regions, where the heating requirement is approximately 107  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 large-scale solar flares (~1034  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 computationally-intensive process for more than a few examples. Hence, Bareford et al. (2010, 2011), used relaxation theory to identify the relaxed states for all of the threshold (i.e., marginally unstable) configurations.

An unstable field obeys relaxation theory if it relaxes towards the lowest energy state that conserves total magnetic axial flux and global magnetic helicity (Taylor 1974, 1986). This minimum energy state is a linear force-free field; i.e., α(r) is invariant and ∇ × B = αB. The helicity (K) indicates how intertwined a magnetic field is with itself (Berger 1999). Coronal loops have a non-zero normal flux at the footpoints; hence, the gauge-invariance of relative helicity (Berger & Field 1984; Finn & Antonsen 1985) makes it the more useful property: (1)Where A is the magnetic potential, B′ is the potential field with the same boundary conditions and A′ is the corresponding vector potential. Helicity-conserving relaxation has been seen in many laboratory experiments (Heidbrink & Dang 2000; Taylor 1986). It must be noted that helicity is subject to global resistive diffusion. Nevertheless, if localised dissipation occurs on small spatial scales (i.e., across shock fronts or within thin current sheets), the reduction in helicity will be negligible compared to the decrease in magnetic energy (Browning 1988; Browning et al. 2008). The original intention of relaxation theory was to explain laboratory plasma phenomena; but latterly, it has been frequently applied to the solar corona (Heyvaerts & Priest 1984; Browning et al. 1986; Vekstein et al. 1993; Zhang & Low 2003; Priest et al. 2005).

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

A long-standing problem has been how to apply relaxation theory in astrophysical contexts without the presence of conducting walls: simplistically, the relaxation should extend out to infinity and lead always to potential fields. Browning (1988) and Dixon et al. (1989) showed that relaxation theory could apply to volumes with free boundaries, but did not give a prediction for the spatial extent of the relaxed state. We find that a modification to Taylor relaxation theory is required before it can be used to estimate the energy released by a kink-unstable loop. In contrast to previous work, we calculate the helicity directly at specific times during each simulation (Browning et al. 2008 integrated the time differential of the helicity in order to show the change in δ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 second-order accurate predictor-corrector step that also incorporates artificial viscosity, ensuring shocks are captured accurately. Van Leer (1997) gradient limiters are used at the remap step in order to preserve monotonicity. The divergence-free condition (∇·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 cf 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 Rankine-Hugoniot jump conditions and the values of the coefficients have been chosen such that numerical oscillations behind shock fronts are prevented. Note, the force 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, Rb is the loop radius, ρ0 the initial mass density, and B1 the initial axial field at r = 0. The other variables are expressed in a similar manner; where L = 20  Rb is the loop length, tA = Rb/vA 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 Rb, ρ0 and B1: Values appropriate for a coronal active region can be obtained by setting Rb = 1  Mm, ρ0 = 1.6726 × 10-13  kg   m-3 and B1 = 50  G. Hence, the length becomes 20  Mm, vA ≈ 10  Mm   s-1 and η0 ≈ 4π × 106  Ω   m.

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

The computational domain is a 3D staggered grid: physical variables are not calculated at the same place for each cell in the domain. This approach improves numerical stability and allows conservation laws to be included in the computation. The domain size is Lx = Ly = 6 (− 3:+3) and Lz = 20 (−10:+10). Initially, the loop axis follows the z-axis and the loop radius is r = 1; therefore, the simulated loops all have an aspect ratio of 20. 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 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: 1282 × 256 (low) and 2562 × 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 line-tied coronal loop, see Gerrard et al. (2002), Gerrard & Hood (2003), Browning et al. (2008) and Hood et al. (2009).

thumbnail 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 (R1:R2:Rb:RB = 0.5:0.9:1:3, where RB is the distance from the initial loop axis to the edge of the envelope). The loop aspect ratio (L/ Rb) 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 current-neutralising layer so as to ensure the loop has (at least initially) zero net current: this improves on the model used by Browning & Van der Linden (2003) and Bareford et al. (2010), which allowed loops to have net current (i.e., an azimuthal field was usually present in the potential envelope). All currents are now created by convective motions local to the loop footpoints. Hence, a current neutralisation layer is introduced here, defined such that the azimuthal field (Bθ) always falls to zero at the loop boundary (Rb); therefore, Bθ is zero in the potential envelope. The loop’s radial α-profile is approximated by a piecewise-constant 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 RB = 3 (three times the loop radius); this is far enough away that the boundary conditions do not influence the plasma evolution.

Table 1

α-profiles, axial fluxes and field coefficients for the simulated loops.

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

where (i = 1,2,3) represents the sign of αi. The fields must be continuous at the inner radial boundaries, R1, R2 and Rb. (The positions are R1 = 0.5, R2 = 0.9 and Rb = 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 R2 and Rb.) Therefore, the coefficients Bj and Cj (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 B3θ(Rb) = 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 B1 = 1.

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

thumbnail Fig. 2

Linear kink instability thresholds for L/ Rb = 20. These thresholds have been cropped by a pair of dashed lines that indicate where Bz(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 close-fitting conducting shell (Browning & Van der Linden 2003), since in the case of zero net current, perturbations quickly fall to zero beyond the loop radius. A consequence of the Bessel function model is that the axial field changes sign if the values of α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 Bz(r) of uniform sign. The new stability space is closed by excluding those field profiles that have a Bz(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.

thumbnail Fig. 3

Initial analytically-determined Bz (solid) and Bθ (dashed) profiles for Luni and Lmix. 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 Luni (uniform twist) and Loop D will be denoted by Lmix (mixed twist).

Each of the loops indicated in Fig. 2 is subjected to a disturbance of the form, (17)where vr 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 c1 = 0.01 reduces the amplitude so that the perturbation is initially linear. This disturbance initiates the kink instability.

Both Luni and Lmix 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 Luni (α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 Luni only is much less (Wb ≈ 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 non-dimensional axial flux over the radius RB (Table 1, Col. 5). For loop Luni, this gives a multiplier of 4.8 × 1019 (assuming a radius of 1   Mm and a mean axial field strength of 50  G): thus, the dimensionalised initial loop energy is roughly 1021  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 (Emax ≈ 0.2). The maximum current, Jmax, also rises just before the decrease in W, indicating the formation of a helical current sheet. The nonlinear instability starts at around t = 50tA and within the next 50tA, approximately 70% of the total energy release has been achieved. Magnetic energy reduces more slowly after t = 100tA. The release of magnetic energy is of a similar size for both resolutions, and significantly, larger currents are recorded at high resolution, which is expected; otherwise, spatially-confined changes in current are missed and anomalous resistivity is reduced. The fact that the maximum current increases with resolution is indicative of current sheet formation. These structures have (possibly) infinite current density, so higher resolutions should reveal larger values of the maximum.

In the top two cases of Fig. 4, ideal MHD and zero background resisitivity, it can be seen that there is no detectable change in energy until around 40tA – 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.

thumbnail Fig. 4

Loop Luni: 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 (Jcrit = 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 non-zero 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 Luni; 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.

thumbnail Fig. 5

Loop Lmix: 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 300tA; this reduction is three orders of magnitude smaller than the energy release caused by the kink instability.

Loop Lmix (α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 Luni. Note, the initial magnetic energy is higher than the value given for the other loop. This is a consequence of setting B1 = 1: Lmix 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 2562 × 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 current-dependent resistivity and zero background resistivity – this is more compatible with the coronal environment. The breakdown of energy conservation associated with this parameter choice can be ignored since we are only interested in how the magnetic field behaves during and after the instability. However, this issue will have to be addressed for detailed studies of the thermodynamic evolution.

3.2. Magnetic field and current magnitude

Now we examine the magnetic field (and critical current distribution) at specific times during the simulations. Figure 6 shows how field lines, originally located within the core, become kinked as the instability takes hold. The dark grey field lines are drawn from the bottom boundary (or footpoint) and the light grey field lines are drawn from similar locations at the upper boundary. Loops Luni and Lmix 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   cf + ν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 Luni than for Lmix. In terms of azimuthal field, Lmix is the weaker of the two (Fig. 3), however, some importance should be attached to the fact that Lmix 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 Luni (e.g., Loop E), the energy release increases only slightly (δW ≈ 2.6).

Figure 7 shows cross sections of Luni at z = 0 (halfway along the loop) and Lmix at z = −2, which is roughly the centre of the only patch of significant viscous heating for Lmix at t = 60tA, 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 ≈ 50tA) 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 x-y plane, which are consistent with a cylindrical configuration, bounded by a current-neutralising layer: the arrows follow each other and the arrow sizes initially increase away from the axis, and then diminish before the loop edge.

thumbnail 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 = 60tA for Luni (top) and Lmix (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

thumbnail 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 Luni (top) and at z = −2 for Lmix (bottom).

Open with DEXTER

thumbnail Fig. 8

Spatial variation of current magnitude across the loop cross section at the final time of the simulation for Luni (left) and Lmix (right). All currents are now well below the critical value. The arrows are magnetic field vectors.

Open with DEXTER

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

3.3. Helicity conservation

DeVore (2000) showed how to calculate the magnetic helicity over an entire coronal volume above a photospheric bounding surface. The first step is to work out the magnetic vector potential for a current-free field that has the same distribution of vertical magnetic flux at the lower boundary. DeVore begins by deriving an expression for the scalar potential, using Green’s function for Laplace’s equation as the integration kernel, (18)where . The grid domain used by LARE3D has a Cartesian geometry: the coronal loop is initially represented as a straight cylinder within a rectangular box. The 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 gauge-invariant vector potential can be specified as, (21)by subtracting the helicity due to the potential field, and Eq. (21) can be re-expressed by expanding the cross product of the second term, (22)Finally, the gauge-invariant magnetic helicity is (23)The geometry used by DeVore differs significantly from that used here (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 five-point Newton-Cotes 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.

Table 2

Helicity at three times during the simulations (ηb = 0) for all five kink-unstable 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 (Kinitial = 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 ΔKW ≈ 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 (RB). 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: Rb ≤ r ≤ RB. Some form of partial relaxation is more likely to be relevant to the zero-net-current case; it is known from simulation results, presented here and in Hood et al. (2009), that reconnection is of limited extent, leaving much of the external field undisturbed. This contrasts with previous results for loops carrying net current (Browning et al. 2008), in which the disturbances in the nonlinear phase of the kink generally extended to the boundary of the simulation.

Bareford et al. (2011) maintained zero net current by fixing a neutralising loop surface at a specified radius (r = Rl). 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 Rl and RB. Naively, it might have been expected that the partially-relaxed 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 partially-relaxed state to include a region of azimuthal field extending to infinity, which is the consequence of allowing net current, since in the initial state, the azimuthal field is zero outside the loop. The existence of current sheets (here broadened into a narrow current-neutralising layer) bounding localised relaxed states has also been noted by Gimblett et al. (2006). The axial flux is conserved such that ψl of the threshold state is equal to ψl of the relaxed state. The subscript l denotes the relaxation radius, Rl; 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 = Rl (similarly, KB is the helicity from r = 0 to r = RB, 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 (Rl) 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   Rb (as for Luni), then the axial field would drop by a factor of 1.82, 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.)

thumbnail Fig. 9

Unstable and relaxed loop. K/ψ2 is conserved over the region enclosed by the dashed circle, which is the outer boundary (Rl) 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 ≤ Rl of the relaxed loop (Fig. 9, right). Also, helicity is absent from the potential envelope surrounding a zero-net-current loop; hence, Kb = KB (unlike magnetic energy, Wb ≠ WB). 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 Wl(αi1,αi2) is the energy of the threshold state and Wl(αl) is the energy of the relaxed state.

thumbnail Fig. 10

Variation with relaxation radius (Rl) of relaxed alpha (αl, left) and of dimensionless energy release (δW, right) for Luni (solid line) and for Lmix (dashed line).

Open with DEXTER

The relaxation radius (Rl) 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 Rl, however, this relationship is not linear; beyond a moderate expansion of 50% (Rl = 1.5) the energy release is ~99% of its maximum for Lmix and ~80% for Luni. 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 Rl ≳ 1.5.

4.1.1. Numerical relaxed states

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

thumbnail Fig. 11

Loop Luni: a comparison between the Bx (top row), By (middle row) and Bz (bottom row) magnetic field profiles obtained numerically (red line) and analytically (black line). The latter is calculated from the αl and Rl that best fit the numerical plot, which is taken from the final frame (t = 400tA) 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 300tA; 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|)   B1J1(| αl|r) and r ≤ Rl. We generate relaxed configurations for every value of Rl between 0.9 and 3.0 in increments of 0.01. Then, for each field component, it is checked which of the 211 possibilities has the lowest chi-squared value when compared with the numerical values for the same component. These field comparisons (Bx, By and Bz) 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 z-dependent 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 analytical-numerical comparisons for Luni at y = 0.

In general, the analytical plots, such as the ones in Fig. 11, show a good agreement with the numerics; however, the Rl and αl that best fit the red lines are independently chosen for each of the forty-five subplots (there are fifteen y - z positions and three field components). Figure 12 (left) shows the variation in these best-fit parameters: each symbol, by virtue of its position, gives the Rl − α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 x-axes (one for Rl, the other αl) means that all forty-five symbols can be clearly distinguished (the y-axis has no meaning beyond a simple sorting of the data points).

thumbnail Fig. 12

Loop Luni: the best-fit Rl − αl pairs when ηb = 0 (left) and when ηb = 10-4 (right). Black circles are for those best fits determined from Bz profiles, red plus signs are for Bx and blue crosses By.

Open with DEXTER

The derived averages are ⟨ Rl ⟩  = 1.76 ± 0.29 and ⟨ αl ⟩  = 0.28 ± 0.31; these two properties have a one-to-one 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(Rl) is small when Rl ≈ 1.8 (Fig. 10, right). Figure 12 suggests that the final numerically-calculated state is tending towards a relaxed loop that can be characterised by a localised invariant α-profile. (Note, all points would lie on a single vertical line if there were a perfect match between the analytical and numerical models.) The outlying points on the left have αl values that are further from ⟨ αl ⟩ than those on the right, since the relationship between αl and Rl 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 re-running the simulation with background resistivity; once again, ηb = 10-4. Figure 12 (right) yields ⟨ Rl ⟩  = 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 low-level noise and restricts the current, and thereby reduces the deviation associated with the analytical fit.

Table 3

Analytical and numerical comparison for the kink-unstable loops A–E.

thumbnail Fig. 13

Value of α (left) and plasma-β (right) over the x-y plane for Luni (ηb = 0). The grey circle approximates the loop cross section at z = 0 and t = 400tA.

Open with DEXTER

thumbnail Fig. 14

Left, the cosine of the angle between the current and the magnetic field over the x-y 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 Luni, also extends to other positions along the instability threshold, see Table 3. This table also confirms that the analytically-calculated helicities of Loops A–E are in approximate agreement with the numerical values, which were derived according to the procedure discussed in 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 locally-constant α profile. We show the value of α computed over the midplane (z = 0) of Luni, 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 anti-parallel (black) to the field for a high percentage of the cross sectional area (positions that are far from parallel, 20° < jB < 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 force-free state (although Taylor relaxation would predict a uniform pressure distribution). However, it seems that the loop is heading towards an approximate balance of forces (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 anti-parallel 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 zero-net-current cylindrical loops. These loop configurations have been identified by a linear analysis as being marginally kink unstable (Bareford et al. 2011). The simulations show that the instability quickly enters a nonlinear phase and magnetic energy declines sharply before leveling off. Furthermore, the amount of energy released matches the amount predicted by Taylor relaxation (Table 3), taking account of the fact that the relaxation is localised. Evidence for helicity conservation was presented, and the change in helicity was shown to be much smaller than the drop in magnetic energy. The implication of this result is that energy diffusion is occuring on small scales compared to the global length scale; i.e., within shocks associated with magnetic reconnection. The low values of kinetic energy during the unstable phase (compared to the internal energy) imply that these shocks occur near to reconnection sites. (We find that dissipation within current sheets only becomes significant when η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: slow-mode shocks are expected since only shocks of this type can reduce magnetic field strength.

Relaxation theory also predicts that the final relaxed state should have a constant α-profile. Although the final numerical α-profile still retains much fine structure, the final magnetic fields are well-modelled by a (localised) constant-α profile with some fluctuations superposed – this suggests the fine-scale structure is self-cancelling (i.e., it integrates out). Furthermore, the energy of the final numerical state is very well matched by the energy of the same constant-α state (the field energy is insensitive to the spikeness of the numerical data). The property of zero net current is retained after the instability. Typically, the loop expands radially, the field reconnecting with that present in the potential envelope. These results justify the choices made by Bareford et al. (2011) regarding the details of the relaxation process. Nevertheless, we were concerned that the thinness of the current neutralisation layer (Fig. 1) might have influenced the results. Hence, we also ran a resistive MHD simulation (with ηb = 0) for a zero-net-current 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 Lmix. Numerically, the energy release was 1.2, which was again consistent with the analytically-determined value.

It appears that the assumption of a Taylor-relaxed state, subsequent to a kink instability, has been verified by the work presented here. The relaxation does not extend over the full numerical volume, but over a region of smaller extent (out to a radius Rl, which is less than the full radius, RB). 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 Rl can be precisely determined from the field configuration at instability onset. However, the analytical work has revealed that for marginally-unstable loops, the energy release varies little with relaxation radius once Rl ≥ 1.5 (Fig. 10); hence, a calculation of the energy release does not necessarily require a precise prediction for Rl.

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

An issue for further research concerns the interaction of convective driving with the relaxation process. The LARE3D code could be used to help resolve this issue. It should be possible to choose a loop configuration (i.e., a set of α-parameters) that is just inside the threshold for linear kink instability and then, trigger the instability by applying a pre-determined velocity profile (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 re-simulated within LARE3D. This would reveal what effect, if any, curvature has on the energy release precipitated by kink instability. Of course, this procedure could also be applied to other improvements; e.g., gravity (with ρ(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

All Tables

Table 1

α-profiles, axial fluxes and field coefficients for the simulated loops.

Table 2

Helicity at three times during the simulations (ηb = 0) for all five kink-unstable loops.

Table 3

Analytical and numerical comparison for the kink-unstable loops A–E.

All Figures

thumbnail 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 (R1:R2:Rb:RB = 0.5:0.9:1:3, where RB is the distance from the initial loop axis to the edge of the envelope). The loop aspect ratio (L/ Rb) in this figure is 20.

Open with DEXTER
In the text
thumbnail Fig. 2

Linear kink instability thresholds for L/ Rb = 20. These thresholds have been cropped by a pair of dashed lines that indicate where Bz(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
thumbnail Fig. 3

Initial analytically-determined Bz (solid) and Bθ (dashed) profiles for Luni and Lmix. Field profiles for the stable loop (last row of Table 1) are also shown.

Open with DEXTER
In the text
thumbnail Fig. 4

Loop Luni: 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 (Jcrit = 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
thumbnail Fig. 5

Loop Lmix: 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
thumbnail 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 = 60tA for Luni (top) and Lmix (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
thumbnail 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 Luni (top) and at z = −2 for Lmix (bottom).

Open with DEXTER
In the text
thumbnail Fig. 8

Spatial variation of current magnitude across the loop cross section at the final time of the simulation for Luni (left) and Lmix (right). All currents are now well below the critical value. The arrows are magnetic field vectors.

Open with DEXTER
In the text
thumbnail Fig. 9

Unstable and relaxed loop. K/ψ2 is conserved over the region enclosed by the dashed circle, which is the outer boundary (Rl) of the relaxed loop.

Open with DEXTER
In the text
thumbnail Fig. 10

Variation with relaxation radius (Rl) of relaxed alpha (αl, left) and of dimensionless energy release (δW, right) for Luni (solid line) and for Lmix (dashed line).

Open with DEXTER
In the text
thumbnail Fig. 11

Loop Luni: a comparison between the Bx (top row), By (middle row) and Bz (bottom row) magnetic field profiles obtained numerically (red line) and analytically (black line). The latter is calculated from the αl and Rl that best fit the numerical plot, which is taken from the final frame (t = 400tA) 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
thumbnail Fig. 12

Loop Luni: the best-fit Rl − αl pairs when ηb = 0 (left) and when ηb = 10-4 (right). Black circles are for those best fits determined from Bz profiles, red plus signs are for Bx and blue crosses By.

Open with DEXTER
In the text
thumbnail Fig. 13

Value of α (left) and plasma-β (right) over the x-y plane for Luni (ηb = 0). The grey circle approximates the loop cross section at z = 0 and t = 400tA.

Open with DEXTER
In the text
thumbnail Fig. 14

Left, the cosine of the angle between the current and the magnetic field over the x-y 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

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.