Dynamics of braided coronal loops
I. Onset of magnetic reconnection
A. L. WilmotSmith  D. I. Pontin  G. Hornig
Division of Mathematics, University of Dundee, Dundee, DD1 4HN, UK
Received 11 January 2010 / Accepted 2 April 2010
Abstract
Aims. The response of the solar coronal magnetic field to
smallscale photospheric boundary motions including the possible
formation of current sheets via the Parker scenario is one of open
questions of solar physics. Here we address the problem via a numerical
simulation.
Methods. The threedimensional evolution of a braided magnetic
field which is initially close to a forcefree state is followed using
a resistive MHD code.
Results. A longwavelength instability takes place and leads to
the formation of two thin current layers. Magnetic reconnection occurs
across the current sheets with threedimensional features shown,
including an elliptic magnetic field structure about the reconnection
site, and results in an untwisting of the global field structure.
Key words: magnetohydrodynamics (MHD)  plasmas  Sun: corona  magnetic fields
1 Introduction
Parker's notion of the ``topological dissipation'' of coronal magnetic fields (Parker 1972) continues to generate much debate. Simply put, Parker's suggestion is that following boundary motions of sufficient complexity the magnetic field of a coronal loop will be unable to ideally relax to a smooth forcefree equilibrium and instead tangential discontinuities in the field, corresponding to current sheets, will develop. In general terms the possible outcomes of relaxation are a development of singular current sheets (e.g. Ng & Bhattacharjee 1998; Janse & Low 2009), development of thin but nonsingular current layers (e.g. Longcope & Strauss 1994; Galsgaard et al. 2003), and a smooth equilibrium with largescale current features (e.g. van Ballegooijen 1985; Craig & Sneyd 2005; WilmotSmith et al. 2009a). It should also be noted that the distinction between the first two cases may be difficult to determine numerically.
In a previous a paper (WilmotSmith et al. 2009a) we considered the ideal relaxation of a braided magnetic field towards a forcefree equilibrium. The magnetic field configuration was based on the pigtail braid and imposed as an initial condition (rather than being built up through boundary motions). This particular braid was chosen since it is the simplest nontrivial braid with no net twist and, accordingly, is the most conservative realistic case to examine. More complex braided fields will, in general, contain components of the pigtailtype. There is ample motivation for modelling loops as having braided components. Photospheric turbulence subjects loop footpoints to a random walk, with motions of the fragments about each other acting to braid (or unbraid) the overlying loop. However, while there is some evidence for the existence of braided loop configurations (see the lefthand image of Fig. 1 for example), most coronal loops appear to be close to a potential field (e.g. Fig. 1, righthand image). This may be an effect of the large aspect ratios typical to loops, since a winding of a field line around another field line is almost undetectable when smoothed out over the length of the loop (see also Berger & AsgariTarghi 2009). Another reason could be that reconnection is very efficient in maintaining a low degree of topological complexity in loops. The present work is designed in part to test whether this is the case.
Figure 1: TRACE coronal loops. ( Left) A largescale tangled configuration and ( right) apparently smooth loops. 

Open with DEXTER 
In WilmotSmith et al. (2009a) an ideal Lagrangian relaxation scheme (Craig & Sneyd 1986; Pontin et al. 2009) was used to ideally evolve the braided field towards a forcefree state. A smooth near forcefree equilibrium was attained with largescale current features i.e. without any tangential discontinuities or even strong current concentrations. (This situation may be compared with Parker (1994) where the same pigtail braid situation is considered as a thought experiment and concluded to inevitably lead to tangential discontinuities. It appears that Parker's assumption of in the domain fails; we indeed have , where S is a crosssectional surface through the braid, but varies significantly between field lines).
Although the local current density in the ideal equilibrium is of largescale, a global quantity, the integrated parallel current, , has small scales (here the parallel indicates parallel to the magnetic field and the integral is taken along magnetic field lines). By small scales we mean that varies significantly between neighbouring field lines as the field lines pass though a particular plane, such as the lower boundary. In WilmotSmith et al. (2009a) it was suggested that these small scales in could lead to the development of a resistive instability and so a loss of equilibrium of the field. This mechanism would be distinct to that of Parker's topological dissipation but have the same consequence.
The consideration that for sufficiently small scales in the integrated parallel current then a forcefree magnetic field will be resistively unstable (the instability being dependent, for a given scale in the integrated parallel current, on the value of the resistivity) motivated us to put the end, near forcefree, state of the Lagrangian relaxation into a resistive MHD code and test whether or not the state is stable (for various values of resistivity allowed by numerical limitations). It turns out that the braided field as implemented in the resistive code is not stable. This finding, together with the subsequent evolution of the field, allows us to address a number of key questions related to MHD and the behaviour of coronal magnetic fields. The results are described in this series of papers. Here we describe the details relating to the numerical setup of the problem and the early evolution of the system. A subsequent paper addresses the long term evolution of the system.
2 Numerical scheme and simulation setup
2.1 Numerical scheme
The numerical scheme employed in the simulations that follow is described
briefly below (further details may be found in Nordlund & Galsgaard (1997) and at
http://www.astro.ku.dk/ kg). We solve the threedimensional resistive MHD
equations in the form
where is the magnetic field, the electric field, the plasma velocity, the resistivity, the electric current, the density, the viscous stress tensor, P the pressure, e the internal energy, the viscous dissipation and the Joule dissipation. An ideal gas is assumed, and hence . These equations have been nondimensionalised by setting the magnetic permeability , and the gas constant ( ) equal to the mean molecular weight (M). The result is that for a volume in which , time units are such that an Alfvén wave would travel one space unit in one unit of time.
The Eqs. (1)(6) are solved on staggered meshes; with respect to a mesh on which and e are defined at the body centre of the cell, and P are defined at face centres and and at edge centres. In this way the required MHD conservation laws are automatically satisfied. Derivatives are calculated using sixthorder finite differences that return a value which is shifted half a gridpoint up or down with respect to the input values. When the staggered mesh means that some quantity must be interpolated, data values are calculated using a fifthorder interpolation method at the relevant position. A thirdorder predictorcorrector method is employed for timestepping.
In all simulation runs we employ a spatially uniform resistivity model. Viscosity is calculated using a combined secondorder and fourthorder method (sometimes termed ``hyperviscosity''), which is capable of providing sufficient localised dissipation where necessary to handle the development of numerical instabilities (Nordlund & Galsgaard 1997). The effect is to ``switch on'' the viscosity where very short length scales develop, while maintaining a minimal amount of viscous dissipation where the velocity field is smooth.
2.2 Creating the initial condition
As discussed in Sect. 1, the initial state for the resistive MHD simulations is drawn from the final state of the Lagrangian relaxation experiment of WilmotSmith et al. (2009a). The quantities previously known from the Lagrangian code (see Craig & Sneyd 1986) are the magnetic field and the current and in the near forcefree relaxed state these are known on a highly distorted mesh. We describe below the process of constructing the field on the rectangular grid required for the resistive MHD simulations.
In order to ensure that the initial magnetic field is divergencefree, we work first with the vector potential
for .
In the relaxation scheme, the calculation of
requires only that we know the initial vector potential before
relaxation and the mesh deformation. In terms of the initial mesh
and the final ``relaxed mesh'' ,
the ith component of the final vector potential
is given (see Appendix) in terms of the initial vector potential by
To create the input magnetic field for the MHD simulations we then interpolate this vector potential onto a rectangular grid. Since the magnetic field components are facecentred on the staggered grid, the vector potential components are interpolated onto locations corresponding to edge centres of the desired grid. We then obtain the magnetic field by taking the curl of using the sixthorder finite differences described above, which yields magnetic field components at facecentres as required. An interpolation scheme using biharmonic spline radial basis functions was applied to , the particular scheme chosen to maximize the smoothness of the corresponding current density , which involves second derivatives of . To further improve this smoothness a simple fivepoint smoothing algorithm was finally applied to , before taking the curl.
The result of the above is that the initial braided magnetic field for our MHD simulations is divergencefree to accuracies on the order of truncation errors of the sixthorder finite differences (with typical maximum within the domain). The topology of the magnetic field turns out to be well conserved by the process, another important consideration for the experiment. However, a drawback is that the quality of the forcefree approximation is not perfectly maintained; the initial state is further from forcefree than the relaxed field of the Lagrangian experiment. Details and implications are discussed in Sects. 3 and 5.
Figure 2: Initial state of the simulation: ( left) Isosurface of current at 25% of the domain maximum and ( right) some particular magnetic field lines illustrating the braided nature of the field. 

Open with DEXTER 
2.3 Initial state
The initial magnetic field is given on a domain of size [24,24] in the vertical (z) direction and [6,6] in both the horizontal (x,y) directions. Covering this domain we take a uniform mesh of 320^{3} cells and employ closed boundary conditions in all three directions. The magnetic field is linetied, and can be very closely approximated by at the boundaries, i.e. it is directed perpendicular to the z boundaries and parallel to the x and y boundaries. To achieve the perpendicular condition the Lagrangian relaxation experiment was rerun on the larger horizontal domain (now [6,6] ^{2} compared with [4,4]^{2} in WilmotSmith et al. 2009a). The braided field is centered in the middle of the domain and the field is close to uniform in the region external to the braid. Accordingly the results presented here are shown for only the subsection of the full domain in which the important dynamics occur, specifically .
An isosurface of current in the initial state is given the lefthand image of Fig. 2. The current has largescales (see also the upperleft hand image of Fig. 5 where contours of current in a horizontal crosssection are shown) with two fingers of current extending vertically through the domain. Some sample field lines further illustrating the nature of the field are given in the righthand image of the same figure. Although nontrivial, the initial state has little magnetic energy in excess of the homogeneous field ( in ). The aspect ratio employed in the model is 1:8. Although this is larger than that of many previous simulations it is smaller than a realistic ratio for a coronal loop (1:50, say). In the configuration the poloidal field components are small compared with the toroidal components so that the field lines look almost straight. This level of braiding would be observationally difficult to distinguish from a potential field.
To initialise the simulation the dimensionless plasma density (Sect. 2.1) is set at throughout the domain and the internal energy as e=0.1. The result is a plasma at t=0 that lies in the range . For the results described in the main section of this paper (Sect. 3) we consider the early evolution of the system (up to t=14) with time measured in units of the Alfvén time. A uniform resistivity of has been taken and the effect of changing the resistivity is discussed at various points in the following text.
Figure 3: Maximum absolute value of the current a); velocity b); Lorentz force c); and total kinetic energy d) in the domain with time for a sequence of decreasing uniform resistivies as indicated in the figures. 

Open with DEXTER 
Figure 4: Isosurfaces of current, , at 50% of the domain maximum for a sequence of increasing times (from left to right, t=4,6,8,10,12) showing the formation of the two initial current layers. 

Open with DEXTER 
3 Results
3.1 Basic properties
Figure 3 shows the maximum absolute values of the current, velocity and Lorentz force, and the kinetic energy ( ) for the time interval ( ) under consideration. The domain taken in all cases is the central section, , of the full box. The variation in quantities is shown for a sequence of uniform resistivities decreasing by over an order of magnitude, specifically and .
The initially high maximum Lorentz force, decreases rapidly over the first few time units. Both the high value and the decrease are artifacts of the method used to create the initial state. The interpolation required to transfer the state to the Eulerian grid (as described in Sect. 2.2) results in some noise in the initial magnetic field and current density. Some noise persists even with the application of a smoothing algorithm to the vector potential for the magnetic field and this is particularly noticeable in the Lorentz force rather than the magnetic field and current density alone. Thus whilst the final state of the Lagrangian relaxation experiment had a very small maximum Lorentz force (specifically ), the initial state here is further from forcefree. The decrease over then arises though a smoothing of the noise in the initial state.
Turning now to the remaining quantities shown in Fig. 3 two primary features are found. Firstly a growth in the kinetic energy and maximum velocity for occurs, and this growth is independent of resistivity, . For the remaining time considered there is no significant change in kinetic energy. Secondly, a linear growth in the current density occurs for . The rate at which the maximum current density increases is higher for lower resistivity, . At t=12 the maximum current density is achieved; this maximum is higher for lower resistivity but for all three resistivities the growth phase ends at the same time.
The lack of dependence of kinetic energy on resistivity may suggest an ideal instability is present. The subsequent linear growth of current would then be a nonlinear consequence of this instability rather than its initial appearance. This growth is clearly dependent on resistivity, being slower for higher values of . Little is known about the nonlinear phase of instabilities and such a dependence may still be consistent with an ideal instability with a nonlinear phase damped by . Strong conclusions are clearly difficult to draw at this stage. An additional consideration is that the implementation of the field on the new grid has resulted in significant Lorentz forces in the initial state. We return to these questions in Sect. 4 but now proceed to consider the the nature of the currents within the domain, now fixing .
3.2 Formation of current layers
Figure 4 shows isosurfaces of current at 50% of the domain maximum ( ) for a sequence of increasing times (for the initial state see Fig. 2). In the early stages ( ) the current diffuses slightly while maintaining its large scales in the perpendicular direction. A symmetric evolution follows and after the phase of current growth two current concentrations are present, centered at z=3.4 and z=3.6. We call these the two ``initial current layers''.
The formation of these two initial current layers is best illustrated by considering a horizontal cross section through the central plane (z=0). Figure 5 shows contours of the vertical component of current ( ) in that plane at t=0, 6, 12. The zcomponent is taken since it significantly dominates over the two horizontal components, as evident in the shape of the current sheets (Fig. 4). Note that in order to incorporate both sheets the crosssectional plane chosen, z=0, does not intersect the centre of either current sheet and so the magnitude of current in this plane is somewhat low in comparison to the domain maximum. The collapse of the two oppositely signed largescale fingers of current present in the initial state into two thin current sheets of correspondingly the same sign is clearly shown. Also evident is the formation of a weaker current envelope around the braided flux, separating it from the uniform background field. Crosssections of in the horizontal planes through the centres of the two current sheets are shown in the final two images of Fig. 7.
Figure 5: Contours of the vertical component of the current, J_{z}, in the central plane, z=0, at increasing times, illustrating the formation of the two initial current layers (first three images). The lower righthand image shows contours of integrated parallel current along field lines in the initial state at the central plane z=0. From this quantity field line tracers for the locations of initially high integrated parallel current have been determined and marked in crosses in the previous frames, as described in the main text. 

Open with DEXTER 
3.3 Predictors of current layers
The field lines along which these two current sheets form turns out to be well predicted by the regions of high integrated parallel current in the initial state. In resistive MHD the integrated parallel current is related to the integrated parallel electric field via the relation (in the case of a uniform resistivity , as considered here). The integrated parallel electric field is a key quantity for 3D magnetic reconnection; for a localised reconnection region the maximum value of determines the rate of reconnection.
Shown in the lower righthand panel of Fig. 5 are contours of the integrated parallel current, , in the initial state, t=0, where the path of the integral is taken over magnetic field lines. Again the contour is shown in the central plane (z=0). To obtain this contour map, field lines have been integrated through 160^{2} grid of points covering the domain , z=0. Two peaks in the quantity are present and the structure is quite different to that of the current itself in the initial state. We now identify those field lines in this tracing procedure for which the value of is greater than or equal to 75% of the domain maximum, noting the locations where they intersect the lower boundary (where the locations of the field lines are held fixed). For the sequence of times of Fig. 5 we trace field lines starting from those locations on the lower boundary up through the domain and mark with a cross in that same figure their point of intersection with the z=0 plane. Here the difference in colours indicates field lines with positive (black crosses) and negative (white crosses) integrated parallel current, although this distinction is made only to facilitate identification of the locations. It is found that these field lines, traced from the initial locations of high integrated parallel current, are good indicators for the locations of formation of the two current layers. Since the flux on the lower boundary is held fixed these may be identified as the same field lines for as long as the evolution remains ideal. Whilst the evolution will be ideal only during the early stages of the simulation it is clear that the tracers do, nevertheless, provide a useful predictor for the locations of current sheet formation.
Locations of high integrated parallel current are not a commonly used indicator for current sheet formation. Indeed it is quasiseparatrix layers (QSLs), regions where the fieldline connectivity varies strongly (Priest & Démoulin 1995) that are widely thought of as likely sites of current sheet formation. To identify QSLs (as well as their intersections, hyperbolic flux tubes or HFTs) the squashing factor (Titov et al. 2002) is used. Usually denoted by Q, the squashing factor is an indictor of field line connectivity and takes on high values in regions where the field line mapping is strongly distorted. Regions of high Q outline QSLs. As discussed in WilmotSmith et al. (2009b), the braided magnetic field taken as the initial state here contains several QSLs. Contours of the squashing factor, Q, in the central plane (z=0) are shown in Fig. 6 at t=0 (upper) and t=12 (lower). For the calculation again 160^{2} points over the region have been used, a number comparable to the grid resolution.
Figure 6: Contours of the squashing factor, specifically , in the central plane, z = 0, for (upper) the initial state, t = 0 and (lower) the time t = 12. 

Open with DEXTER 
Several regions of high Q are present in both snapshots. The two regions of highest Q in the initial state are associated with the two initial current layers at the later time, that is, the current sheets have formed along QSLs of the field. At the same time, several other regions of high Q are present (both at t=0 and t=12) that are not associated with any particular current features. For example, in the initial state there are eight distinct regions where is greater than 85% of its maximum value but only two of these regions correspond to particular current features at t=12. These results suggest that the integrated parallel electric current and the squashing factor could be used in conjunction for predicting current sheet formation more accurately than by using Q alone.
3.4 Plasma flows and reconnection
We move now to consider the nature of the plasma flows within the domain. First recall that the low values of the plasma beta (Sect. 2) imply the dynamics will be dominated by the Lorentz force rather than the gas pressure. With both the magnetic field and the current having stronger vertical than horizontal components we have that the Lorentz force is primarily in the horizontal direction and so, similarly, are the plasma flows. The Lorentz force drives a flow with a dipolar structure in the six regions of initially strongest current with direction dependent on the sign of the twist in that region. An illustration of such a flow is shown for the plane z=3.4 (a negative twist region corresponding to one of the sites of current sheet formation) in the first (upper left) image of Fig. 7.
The next four images in that figure show the development of the flow for a sequence of increasing times up to the point of maximum current. In the sequence the length of the arrows indicating flow direction have been normalized to each image. The background contours show the vertical component of current in that plane with the colour scale normalised to the current at t=12 (lowerlefthand image), given by the first colourbar. The sequence clearly shows the association of the stagnation part of the dipolar flow with the current intensification. The outflowing plasma from the location of stagnation sets up a counterflow to the initial dipolar structure leading to oppositely directed flows on either side of the weak enclosing current sheath. The final image shows plasma flows and current in the plane z=3.6 at t=12. This crosssection is across the second current sheet and shows the naturally expected inversion of the flow direction. The result is that the global flow structure is dominated by rotational components the direction of which varies both vertically along the structure and on either side (y>0, y<0) of the braided field.
Figure 7: Arrows indicate plasma velocities [V_{x},V_{y}] in horizontal crosssectional planes, superimposed on the vertical component of current, J_{z}. The first five images show structures about the upper current sheet (taking the plane z=3.4) while the final image is for the lower current sheet (taking the plane z=3.6). The images are given at various times, as indicated in the figures. 

Open with DEXTER 
As already indicated by Fig. 3, at t=12 the maximum magnitude of the plasma flows ( ) is a significant fraction of the Alfvén speed ( ). These strongest flows are associated with the global rotation and not with outflows from the two initial current layers (which have an associated ). However, we do expect magnetic reconnection to be taking place across these current sheets and so proceed to consider the nature of this reconnection. For this we concentrate again on the initial current layer centered at z=3.4 (a similar situation occurs about the other current layer) and focus on the structure of the field and flows in the plane perpendicular to the magnetic field at the location of maximum current magnitude, .
Figure 8: Here various quantities are considered in the plane perpendicular to the magnetic field at the region of maximum (negative) current. (Left) Arrows indicate plasma flows in that plane while the coloured contours show the outofplane component of vorticity. (Right) Arrows indicate magnetic field components in the plane superimposed on contours of the outofplane current. The colour table for each is blueyellow for negative  positive. 

Open with DEXTER 
In the lefthand image of Fig. 8 we indicate in more detail the nature of the stagnation flow about the current sheet. Superimposed on the background are contours of the outof plane component of vorticity, i.e. the component of vorticity in the direction of the magnetic field at the location of maximum current. The vorticity shows a quadrupolar configuration around the current sheet.
Perhaps more of a surprise is the structure of the magnetic field in the region. The righthand image of Fig. 8 shows the components of the magnetic field in the crosssectional plane under consideration. The magnetic field is shown to have an elliptic configuration about the current sheet, i.e. about the reconnection region. This finding contrasts with the twodimensional picture of magnetic reconnection under which the process can only occur at a hyperbolic (Xtype) nullpoint of the magnetic field. In threedimensions a much wider variety of possible reconnection sites exist. Reconnection may be associated with 3D nullpoints (e.g. Lau & Finn 1990; Priest & Pontin 2009), magnetic separators (which connect two nullpoints, e.g. Longcope & Cowley 1996; Pontin & Craig 2006; Haynes et al. 2007), or may occur in the absence of any such topological features (Schindler et al. 1988), the latter sometimes termed ``nonnull reconnection''. In particular, the local magnetic field structure need not be hyperbolic but may be elliptic (Hornig & Priest 2003), as recently found in some 3D numerical simulations of reconnection. For example, WilmotSmith & De Moortel (2007) considered reconnection occurring along a quasiseparator and found an elliptic field structure in perpendicular crosssections. The separator configuration of Parnell et al. (2010) showed an elliptic structure along a significant length of the separator under consideration. Parnell et al. (2010) also discussed the separator case theoretically, concluding an elliptic configuration would be a generic situation given a sufficiently strong current density along the separator. Our findings demonstrate that an elliptic field configuration may be present about reconnection sites in the nonnull case. Additionally, tracking these field lines back to the initial setup, an elliptic perpendicular field configuration is again found indicating that locally hyperbolic structures are not necessary for current intensification. As previously discussed, the squashing factor (Q) at the two reconnection sites is very high which demonstrates a further point; regions of highestQ within a domain may have a locally elliptic field configuration.
4 Nature of instability
Evidently the initial magnetic field configuration is not in a stable equilibrium. Since an exact equilibrium of the ideal relaxation code employed in WilmotSmith et al. (2009a) is known to be linearly ideally stable, the lack of stability could arise from one of a number of factors:
 1.
 A resistive instability. The relaxed state of WilmotSmith et al. (2009a) contains small scales in the integrated parallel current. Small enough scales (for a given resistivity) in this quantity are incompatible with a resistive equilibrium.
 2.
 A nonlinear ideal instability not previously found in the ideal evolution of WilmotSmith et al. (2009a) since that evolution only guaranteed linear stability.
 3.
 A lack of equilibrium in the initial state either since:
 The path to equilibrium of the ideal relaxation is fictitious and an exact equilibrium had not been reached. Numerical difficulties (described in Pontin et al. 2009) result in the final state of the ideal relaxation having , i.e. is not perfectly forcefree. Thus the final state of WilmotSmith et al. (2009a) need not be stable (or, indeed accessible via a real MHD relaxation dynamics).
 The technique used to transfer the initial state between codes has perturbed the magnetic field further from forcefree.
With these considerations in mind we examine the evolution of only the middle section of the initial state of the braided magnetic field that is, we cutout the section of the field in the above described experiment that lies in , at t=0. This field is inserted as an initial condition in a new run, now keeping the flux fixed on , the new upper and lower boundaries of the domain. To maintain consistency we use the same resolution in the horizontal direction 320^{2} and a similar effective resolution in the vertical direction, 128 cells over .
Figure 9: (Upper) Vertical component of current in the central plane (z=0) at t=20 for the an MHD evolution on only the central section, of the full field, as described in the main text. (Lower) Maximum current (solid line) and total kinetic energy (dashed line) and maximum Lorentz force (dotdashed line) over time for the same field. 

Open with DEXTER 
In the evolution of this new ``middlethird'' field we find the system adjusts from its initial condition (with zero plasma velocity) to an approximate equilibrium in which the current structure is qualitatively similar to that of the initial state. To illustrate, contours of current in the central plane (z=0) are shown at t=20 in Fig. 9 (upper panel) which may be compared with Fig. 5 (upper left) where the corresponding contours in the initial state, t=0, are shown. The maximum current in the domain is shown as a solid line in the upper panel of Fig. 9 (lefthand axis), the kinetic energy integrated over , as a dashed line (lefthand axis) and the maximum Lorentz force within the domain as a dotdashed line (righthand axis). As in the evolution of the full field, the Lorentz force in this initial state is high as a result of the transfer between grids and decreases rapidly over the first few time units. However, the Lorentz force subsequently stabilises at a low value as the system readjusts to equilibrium. The result that the middle section of the braid alone is in an equilibrium state suggests that an instability is present in the full braided field (rather than a lack of equilibrium due to numerical artefacts). Furthermore the instability is of a longwavelength type with the full structure of the braided field being key.
Returning to the evidence of Figure 3, during the time where the instability is clear in the current growth there is only a very slight dependence of the maximum velocity within the domain on resistivity and no increase in kinetic energy. This effect may be due to the confining nature of the strong background field external to the braided field structure which results in a deflection of the outflowing plasma around the boundary of the braided field (see Fig. 7). Prior to no dependence of the flow on resistivity is seen suggesting an ideal dynamics where the system is adjusting to the distance from equilibrium. The current growth in does have a clear separation according to resistivity . However the increase is linear (rather than exponential) suggesting a dominant nonlinear phase and the growth is slower for higher resistivity suggesting the nonlinear phase is damped by resistive effects.
5 Conclusions
A previous paper (WilmotSmith et al. 2009a) considered the ideal relaxation of a braided magnetic field towards a forcefree equilibrium. Here we have taken the final state of that relaxation and used it as an initial condition for a full resistive MHD simulation.
The braided field is not in a stable equilibrium; two thin current sheets form after a short time (around a quarter of the time for an Alfvén wave to cross the numerical box in the vertical direction). The linear rate of increase of current density and the maximum strength of the current is found to increase with decreasing resistivity although the evolution of the total kinetic energy in the domain is independent of resolution. We conclude that the instability is possibly an ideal one although the details of how it occurs have not been determined. The wavelength of the instability was tested by considering the MHD evolution of the middle third section of the braid alone and this new field was found to be stable. Hence a longwavelength dependence is implied.
The initial configuration contains many regions of high squashing factor, Q, corresponding to quasiseparatrix layers. Plasma flow across QSLs is often thought likely to lead to current sheet formation. The two current sheets that form here do align with two of the highest regions of Q although the remaining regions of high Q do not correspond to any particular current features. The locations of the two current layers are well predicted by the peak values of the integrated parallel current in the initial state. In the central plane this quantity shows two clear peaks and it is at these regions that the two current layers later form.
These two current layers correspond to reconnection sites. Perpendicular to the layers the magnetic field has an elliptic structure, an admissible property of three dimensional reconnection that has only recently been found. The flow about the reconnection sites is of an asymmetric stagnation type, although largescale rotational flows dominate the global structure; the effect of reconnection is not a strong acceleration of the flow but a more subtle untwisting process leading to the change in magnetic field topology.
The longerterm evolution of the system will be considered in a future paper.
Appendix A: Appendix
Equation (7) can be derived from an evolution for the vector potential of a frozenin magnetic field:
This equation is equivalent to the Liederivative for a differential oneform, , associated with the vector field (Hornig 1997). Hence (A.1) can be written as
(A.2) 
and we can express (Abraham et al. 1988), or more conveniently , where maps the initial to the final coordinates and the star indicates the pullback operation. This last equation written in components of the vector potential is (7). Acknowledgements
The authors would like to thank Nasser AlSalti for helpful suggestions and Klaus Galsgaard for the use of his code. A.W.S. and G.H. acknowledge financial support from the UK's STFC. The simulations described were carried out on the STFC and SFC (SRIF) funded linux clusters of the UKMHD consortium at the University of St Andrews.
The Transition Region and Coronal Explorer, TRACE, is a mission of the StanfordLockheed Institute for Space Research (a joint program of the LockheedMartin Advanced Technology Centre's Solar and Astrophysics Laboratory and Stanford's Solar Observatories Group), and part of the NASA Small Explorer program.
References
 Abraham, R., Marsden, J. E., & Ratiu, T. 1988, Manifolds Tensor Analysis and Applications, Applied Mathematical Sciences (New York: SpringerVerlag), 75, 370 [Google Scholar]
 Berger, M. A., & AsgariTarghi, M. 2009, ApJ, 705, 347 [NASA ADS] [CrossRef] [Google Scholar]
 Craig, I. J. D., & Sneyd, A. D. 1986, ApJ, 311, 451 [NASA ADS] [CrossRef] [Google Scholar]
 Craig, I. J. D., & Sneyd, A. D. 2005, Sol. Phys., 232, 41 [NASA ADS] [CrossRef] [Google Scholar]
 Galsgaard, K., Titov, V. S., & Neukirch, T. 2003, ApJ, 595, 505 [NASA ADS] [CrossRef] [Google Scholar]
 Haynes, A. L., Parnell, C. E., Galsgaard, K., & Priest, E. R. 2007, Proc. Roy. Soc. A, 463, 1097 [NASA ADS] [CrossRef] [Google Scholar]
 Hornig, G. 1997, Phys. Plasmas, 4, 646 [NASA ADS] [CrossRef] [Google Scholar]
 Hornig, G., & Priest, E. R. 2003, Phys. Plasmas, 10, 2712 [NASA ADS] [CrossRef] [Google Scholar]
 Janse, Å. M., & Low, B. C. 2009, ApJ, 690, 1089 [NASA ADS] [CrossRef] [Google Scholar]
 Lau, Y.T., & Finn, J. M. 1990, ApJ, 350, 672 [NASA ADS] [CrossRef] [Google Scholar]
 Longcope, D. W., & Cowley, S. C. 1996, Phys. Plasmas, 3, 2885 [NASA ADS] [CrossRef] [Google Scholar]
 Longcope, D., & Strauss, H. R. 1994, ApJ, 437, 851 [NASA ADS] [CrossRef] [Google Scholar]
 Ng, C. S., & Bhattacharjee, A. 1998, Phys. Plasmas, 5, 4028 [NASA ADS] [CrossRef] [Google Scholar]
 Nordlund, A., & Galsgaard, K. 1997, A 3D MHD code for parallel computers. Technical report, Astronomical Observatory, Copenhagen University [Google Scholar]
 Parker, E. N. 1972, ApJ, 174, 499 [NASA ADS] [CrossRef] [Google Scholar]
 Parker, E. N. 1994, Spontaneous Current Sheets in Magnetic Fields With Applications to Stellar XRays. (New York: Oxford University Press, Inc.), 23 [Google Scholar]
 Parnell, C. E., Haynes, A. L., & Galsgaard, K. 2010, JGR, 115, A02102 [NASA ADS] [CrossRef] [Google Scholar]
 Pontin, D. I., & Craig, I. J. D. 2006, ApJ, 642, 568 [NASA ADS] [CrossRef] [Google Scholar]
 Pontin, D. I., Hornig, G., WilmotSmith, A. L., & Craig, I. J. D. 2009, ApJ, 700, 1449 [NASA ADS] [CrossRef] [Google Scholar]
 Priest, E. R., & Demoulin, P. 1995, JGR, 100 (A12), 23443 [NASA ADS] [CrossRef] [Google Scholar]
 Priest, E. R., & Pontin, D. I. 2009, Phys. Plasmas, 16, 122101 [NASA ADS] [CrossRef] [Google Scholar]
 Schindler, K., Hesse, M., & Birn, J. 1988, JGR, 93, 5547 [NASA ADS] [CrossRef] [Google Scholar]
 Titov, V. S., Hornig, G., & Démoulin, P. 2002, JGR, (A8) SSH 31, 1164 [Google Scholar]
 van Ballegooijen, A. A. 1985, ApJ, 298, 421 [NASA ADS] [CrossRef] [Google Scholar]
 WilmotSmith, A. L., & De Moortel, I. 2007, A&A, 473, 615 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 WilmotSmith, A. L., Hornig, G., & Pontin, D. I. 2009a, ApJ, 696, 1339 [NASA ADS] [CrossRef] [Google Scholar]
 WilmotSmith, A. L., Hornig, G., & Pontin, D. I. 2009b, ApJ, 704, 1288 [NASA ADS] [CrossRef] [Google Scholar]
All Figures
Figure 1: TRACE coronal loops. ( Left) A largescale tangled configuration and ( right) apparently smooth loops. 

Open with DEXTER  
In the text 
Figure 2: Initial state of the simulation: ( left) Isosurface of current at 25% of the domain maximum and ( right) some particular magnetic field lines illustrating the braided nature of the field. 

Open with DEXTER  
In the text 
Figure 3: Maximum absolute value of the current a); velocity b); Lorentz force c); and total kinetic energy d) in the domain with time for a sequence of decreasing uniform resistivies as indicated in the figures. 

Open with DEXTER  
In the text 
Figure 4: Isosurfaces of current, , at 50% of the domain maximum for a sequence of increasing times (from left to right, t=4,6,8,10,12) showing the formation of the two initial current layers. 

Open with DEXTER  
In the text 
Figure 5: Contours of the vertical component of the current, J_{z}, in the central plane, z=0, at increasing times, illustrating the formation of the two initial current layers (first three images). The lower righthand image shows contours of integrated parallel current along field lines in the initial state at the central plane z=0. From this quantity field line tracers for the locations of initially high integrated parallel current have been determined and marked in crosses in the previous frames, as described in the main text. 

Open with DEXTER  
In the text 
Figure 6: Contours of the squashing factor, specifically , in the central plane, z = 0, for (upper) the initial state, t = 0 and (lower) the time t = 12. 

Open with DEXTER  
In the text 
Figure 7: Arrows indicate plasma velocities [V_{x},V_{y}] in horizontal crosssectional planes, superimposed on the vertical component of current, J_{z}. The first five images show structures about the upper current sheet (taking the plane z=3.4) while the final image is for the lower current sheet (taking the plane z=3.6). The images are given at various times, as indicated in the figures. 

Open with DEXTER  
In the text 
Figure 8: Here various quantities are considered in the plane perpendicular to the magnetic field at the region of maximum (negative) current. (Left) Arrows indicate plasma flows in that plane while the coloured contours show the outofplane component of vorticity. (Right) Arrows indicate magnetic field components in the plane superimposed on contours of the outofplane current. The colour table for each is blueyellow for negative  positive. 

Open with DEXTER  
In the text 
Figure 9: (Upper) Vertical component of current in the central plane (z=0) at t=20 for the an MHD evolution on only the central section, of the full field, as described in the main text. (Lower) Maximum current (solid line) and total kinetic energy (dashed line) and maximum Lorentz force (dotdashed line) over time for the same field. 

Open with DEXTER  
In the text 
Copyright ESO 2010