Twisted versus braided magnetic flux ropes in coronal geometry
I. Construction and relaxation
Department of Mathematical SciencesDurham University,
Durham
DH1 3LE,
UK
email:
christopher.prior@durham.ac.uk
Received: 23 August 2015
Accepted: 21 November 2015
We introduce a technique for generating tubular magnetic fields with arbitrary axial geometry and internal topology. As an initial application, this technique is used to construct two magnetic flux ropes that have the same sigmoidal tubular shape, but have different internal structures. One is twisted, the other has a more complex braided magnetic field. The flux ropes are embedded above the photospheric neutral line in a quadrupolar linear forcefree background. Using resistivemagnetohydrodynamic simulations, we show that both fields can relax to stable forcefree equilibria whilst maintaining their tubular structure. Both end states are nonlinear forcefree; the twisted field contains a single sign of alpha (the forcefree parameter), indicating a twisted flux rope of a single dominant chirality, the braided field contains both signs of alpha, indicating a flux rope whose internal twisting has both positive and negative chirality. The electric current structures in these final states differ significantly between the braided field, which has a diffuse structure, and the twisted field, which displays a clear sigmoid. This difference might be observable.
Key words: magnetic fields / magnetic reconnection / magnetohydrodynamics (MHD) / Sun: corona / Sun: magnetic fields
© ESO, 2016
1. Introduction
Magnetic flux ropes are observed in the Sun’s coronal region over all activity cycles. They are believed to play a critical role in active region phenomena, such as sigmoid formation, coronal eruptions and coronal mass ejections (e.g. Mandrini et al. 2005; Okamoto et al. 2008; Cheng et al. 2011; Inoue et al. 2015). Additionally flux ropes can be observed in non active regions (e.g. Gibson et al. 2006b; Su et al. 2015). Through modelling it has been argued that flux ropes can either enter the corona by rising through the Sun’s convection zone (e.g. Archontis & Török 2008; Fan 2009; Manchester IV et al. 2004; MacTaggart 2011) or, alternatively, be formed in the coronal region due to gradual shearing (van Ballegooijen & Martens 1989) or strong shearing motions (e.g. Aulanier et al. 2010; Longcope & Beveridge 2007).
Often, however, authors address the subsequent behaviour of a toroidal flux rope placed at the heart of a background field. This approach is particularly prevalent in activeregion modelling where flux ropes with an internal twisted structure are forced into instability through various mechanisms, often leading to ejection (e.g. Fan & Gibson 2007; Kliem et al. 2004, 2010, 2012; Leake et al. 2014). This “kink” instability drives the flux rope to adopt a kinked structure through rotation about its apex (see Fig. 1). Such geometries appear as either S or Z shapes in projection, and are commonly referred to a “sigmoidal” in reference to the assumption that they explain the S and Zshaped sigmoid current structures viewed in softXray emission (e.g. Titov & Démoulin 1999; Gibson et al. 2006a; Green et al. 2007). It should be noted that the simulations mentioned above often yield S or Z shapes with significant symmetry, whilst in general there is a much larger variety of sigmoidal morphologies observed in the corona (Prior & Berger 2012).
Fig. 1 Idealisation of the process by which an initially toroidal tube a) rotates about its apex to form a kinked structure b). This tube in b) appears as an S shape when projected down. The morphology of b) is characteristic of a sigmoidal geometry. 

Open with DEXTER 
In addition to activeregion modelling, a second area of interest in coronal magnetic field modelling concerns the internal structure of tubular fields. A rough interpretation of Parker’s hypothesis (Parker 1972; Janse et al. 2010) is that the unusually high temperature of the coronal region can be explained by the creation of significantly complex field structure due to turbulent photospheric motions. This would lead to the formation of dense regions of current within the field, because of either field discontinuities or sharp field gradients. These sheets provide regions in which reconnection is promoted, enabling magnetic energy to be converted to heat in the surrounding region.
Geometrically, the existence of a large number of small current sheets manifests itself in complex entanglement of the magnetic field lines (see e.g. Janse et al. 2010). Field lines are curves f associated with B through solutions of the O.D.E (1)With this in mind a number of authors have addressed the question: what effect do complex (braided) patterns have on a magnetic field’s evolution? A significant number of studies have considered the problem of developing magnetic braiding through boundary motion. Magnetic fields defined in a Cartesian box, with the initial field structure a set of vertical lines, are subjected to a variety of boundary motions such as shearing flows of various scales. Various relaxation techniques are used to calculate the evolving field and current structure; often these techniques are quasistatic. Typically the evolution leads to the formation of thin layers of current and much interest is focused on the width of these layers (see WilmotSmith 2015). In addition a number of studies have attempted to estimate the amount of heating that braided reconnection could produce, in order to test Parker’s hypothesis (e.g. Galsgaard & Nordlund 1996; Craig & Sneyd 2005; Berger & AsgariTarghi 2009; Ng et al. 2012; Rappazzo & Parker 2013; Yeates et al. 2014; Pontin et al. 2011).
Of particular relevance to this study, van Ballegooijen et al. (2014) explored the relationship between footpoint motions and the possible formation of braided structures, finding braiding to be a dynamic rather than quasistatic phenomenon. This raises the possibility that complex braided structures could be built before the field relaxes to equilibrium. Of most relevance to this study are a series of numerical experiments which compare the effect of preexisting braided structures in a cylindrical domain on the field’s eventual evolution, by comparison to twisted cylindrical fields (e.g. Yeates et al. 2010, 2015; WilmotSmith et al. 2011; Wyper & Pontin 2014). One notable finding was that the braided fields typically released more energy than their twisted counterparts (for the same given starting energy). It was also found that the final end state differed, with the braided field splitting into two separate force free flux tubes of opposing chirality and the twisted field remaining as one tube with a single chirality. To the best of our knowledge there has been no systematic attempt to study the effect of significantly entangled field line configurations in flux ropes with realistic sigmoidal geometry (WilmotSmith 2015). The extra degree of freedom afforded by the (possibly) changing morphology of the rope might have a significant effect on the field’s evolution.
Fig. 2 Visualising the field creation process. a) Depicts a set of braided curves (in the form of the pigtail braid) on the left, and a set of twisted curves on the right. b) Depicts the pigtail braid embedded within a sigmoidal tube shape. c) Depicts the embedding of this tube in a Cartesian coordinate system, over the neutral line of a background quadrupolar flux distribution defined at the plane z = 0 (which plays the role of the photosphere in what follows). 

Open with DEXTER 
It is currently not possible to resolve the precise internal structure of coronal flux ropes using existing observations (Reale 2010). However, it may be possible to infer aspects of the internal structure indirectly if observable properties of flux ropes – such as their stability or overall shape – depend in a clear manner on the internal structure. For example the helical kink instability of twisted flux tubes can be used to link the eruption of toroidal magnetic flux rope to the presence of a twisted internal field structure (e.g. Linton et al. 1996, Titov & Démoulin 1999, Kliem et al. 2004). To determine whether such a clear dependency exists in coronal geometries is the overall goal of this work. There is no clear evidence either way of whether flux ropes with internally braided structure either exist or persist for any length of time in the coronal region, although recent highresolution observations by Cirtain et al. (2013) have been interpreted as the apparent braiding of coronal loop strands.
There is good reason then to consider modelling flux tubes with both realistic sigmoidal shape and complex internal structure, to see whether the observed evolution of a flux rope can reliably yield any information about its internal structure. With this in mind we detail here a technique for generating magnetic fields with arbitrarily complex tubular shape and internal structure. An example is shown in Fig. 2, where a braid is superimposed onto a sigmoidal tube, (a) to (b), a process described in Sect. 2. This tube is then embedded along the neutral line of a boundary flux distribution, (b) to (c), a process described in Sect. 3. Our technique gives complete control over both the global and internal structure of the flux rope, and hence provides a means by which we can explore the effect of internal flux rope topology on its eventual evolution. In Sect. 4 we compare tubes with the same sigmoidal morphology but both braided and twisted internal structures.
2. Creating a tubular field
In this section we describe the mathematical technique for constructing magnetic fields inside tubular domains of arbitrary axis shapes. We consider fields embedded in a Cartesian coordinate system with coordinates (x,y,z). Each field is defined by choosing (i) the tube axis curve; (ii) the tube crosssectional radius as a function of length along the axis; and (iii) a twist parameter describing the internal structure of the tube. In this paper, “twisted” flux ropes consist of a single one of these tubular domains, while “braided” ropes are generated by combining several of these tubular domains.
2.1. Tubular domains
Fig. 3 Tubular domains and curves. a) Depicts the tubular coordinate system T defined in Eq. (2). b) Depicts a tubular domain with u_{3}(s) = 0, ∀s. The curves drawn in this domain are curves of constant R and θ values. We see that the curves follow the shape of the tube. b) Depicts a tubular domain with the same axis curve as a), but with u_{3}(s) ≠ 0 so that the curves of constant (R,θ) value wind internally in the tube. 

Open with DEXTER 
Our tubular domain takes the form of a tube of circular crosssection centred on an axis curve r(s): [ 0,L ] → R^{3}, whose start and end points lie in the z = 0 plane. An example is shown in Fig. 3a. We allow the tube radius R(s) to vary with s; this possibility might be crucial to try to capture the rapid change in tube radius which occurs in the photospherecorona transition region (e.g. van Ballegooijen et al. 2014).
To define a coordinate system in our tubular domain T, we choose an orthonormal framing (d_{1},d_{2},d_{3}) for r using its unit tangent vector d_{3} = r^{′} (total derivatives with respect to s will be denoted with a prime), and a vector field d_{1} which lies in the normal plane of d_{3} (d_{3}·d_{1} = 0,∀s). This basis is completed with the vector product d_{2} = d_{3} × d_{1}. This frame is then extended to form a curvilinear coordinate system by defining a map T(s,ρ,θ): [ 0,L ] × [ 0,1 ] × S^{1} → R^{3} as (2)
2.1.1. Coordinate frame
To fully specify the coordinate system T for a given axis curve r, we must specify d_{1} as a function of s. Intuitively, this controls how much the coordinate frame rotates around the axis as we move along the tube. As the basis {d_{i}} is orthonormal, it is convenient to write its evolution with s in terms of a triplet of scalar functions (u_{1}(s),u_{2}(s),u_{3}(s)) (see e.g. Antman 2005, Chp. 7) where (3)The functions u_{1}(s) and u_{2}(s) define rotations of the coordinate frame due to bending of the axis curve (about two independent directions), and in our case are fixed by the chosen axis curve r. More precisely, given r we define the functions u_{1}(s) and u_{2}(s) in terms of curvature κ and torsion τ of r, with (4)These specific forms for u_{1} and u_{2} are not the only possible choice (which amounts to a choice of the vector d_{1}), but they are the only forms for which the value of u_{3} is independent of the specific form of r (Bishop 1975).
The function u_{3}(s) then identifies the relative twisting of the frame about d_{3}, (with ). It represents the degree of freedom we have in choosing the direction of d_{1} as a function of s. For example, Fig. 3b shows an untwisted frame (u_{3}(s) = 0,∀s) with indicated lines of constant (ρ,θ). By comparison Fig. 3c depicts the same coordinate curves of a twisted coordinate frame (u_{3} ≠ 0). Both curves have the same curvatures (u_{1},u_{2}). The latter, twisted, frame is actually the Frénet frame for this example. To avoid introducing unnecessary twist through our choice of coordinate frame, we do not use the Frénet frame but rather an untwisted frame u_{3} = 0, meaning that γ ≠ 0 in general.
2.1.2. Metric geometry
Much of the analytic calculations which follow are performed in our curvilinear coordinate system T. To employ differential operators in this coordinate system we need to know the metric structure of T. Under our choice that u_{3} = 0, using the partial derivatives, (5)we find the metric tensor g_{ij} = (∂T/∂x_{i})(∂T/∂x_{j}), with x_{1} = s,x_{2} = ρ,x_{3} = θ for this coordinate system, to be (6)For constant tube crosssection R this simplifies to (7)In both cases the Jacobian can be found to be g = ρR^{2}(1 + ρR(u_{1}sinθ − u_{2}cosθ)) = ρR^{2}(1 + ρRκsin(θ − γ)) (using Eq. (4)). Thus we require that R(s) < 1 /κ(s) so that a tubular domain will not overlap itself locally, i.e. the coordinate system is locally well defined (see e.g. Frankel 2011). This does not preclude the tube from overlapping itself nonlocally, though for the purpose of this study we will be choosing curves such that this is not the case.
2.2. Field lines and the unit vector field
We define the magnetic field inside our tubular domain in two stages: first we choose a unit vector field N that determines the field line structure, then we create a divergencefree field B whose field lines are the same as those of N.
We define the unit vector field N in our tubular domain T by specifying its integral curves f [ ρ_{0},θ_{0} ] (s), s ∈ [ 0,L ], which will become the field lines of B. A specific curve is identified by its starting point (ρ_{0},θ_{0}) on the end s = 0 of the tube. Each curve is chosen to lie on a tubular surface of constant ρ, but to twist around this surface as we move along the tube, so that (8)where ψ(s) is a function (the same for all curves) that controls the twisting of the field, relative to the coordinate frame. Thus if ψ(s) ≡ 0 (and R(s) ≡ 0), the curves have the same morphology as the tube axis (the curves in Fig. 3b are such a case). If ψ(s) ≠ 0, then the curves twist around the tube axis. We do not require ψ(s) to be constant or even positivedefinite (the curves in Fig. 3c have ψ(s) = τ(s)). The set F = ^{{}f [ ρ,θ_{0} ]  ρ ∈ [ 0,1 ] ,θ_{0} ∈ S^{1}^{}} will cover the whole tubular domain as s varies from 0 to L. This specification means that the curves satisfy df/ ds> 0 i.e. the curves only pierce each disc of constant s once. Whilst it is perfectly possible to define curves which turn back on themselves along the direction of the tube axis, it is a harder task to define a set of curves which do this and also cover the whole domain. For our purpose the set of field line topologies covered by the definition F allows for sufficient complexity to be of interest, whilst not complicating matters.
Another possible further consideration would be to consider fields for which the functions R and ψ also depend on the ρ and θ coordinates in order to obtain more complex configurations. In what follows we prefer to allow for such further complexity by defining fields composed of multiple tubes (Sect. 2.4); this approach greatly simplifies the analytic considerations of what follows.
2.3. Divergencefree field
Our unit vector field N is turned into a divergencefree magnetic field B = φN by specifying its magnitude φ ≡ B at each point. This is determined entirely by choosing the distribution of φ on the boundary s = 0 of the tube. To see this, note that the divergencefree condition ∇·B = 0 implies that By integrating along a field line f [ ρ_{0},θ_{0} ] (l) (where l is the arclength parameter with N = df/ dl), we find that (11)where φ_{0} is the value of φ at the field line’s base point f [ ρ_{0},θ_{0} ] (0) = (ρ_{0},θ_{0},0).
To obtain an explicit expression for N for the set F, we take the derivative of a field line with respect to s, the arclength parameter of the central curve r. Thus (12)The factor λ accounts for the fact that s will differ in general from l, coinciding only when κ = 0 and u_{3} = 0. The derivative is (13)We may then find ∇·N, for which an explicit expression is derived in Appendix A. Note that if ψ′ = 0 and R′ = 0 then N is already divergencefree and φ is simply conserved along field lines.
2.4. Braided fields
Fig. 4 Braiding geometry. Panel a) depicts a pigtail braid whose centrelines are curves r_{i} embedded in a tube T_{b}. Panel b) illustrates a stirring boundary motion which could be used to create the pigtail structure in a). 

Open with DEXTER 
As mentioned above, we construct “braided” fields by combining n tubular fields of the type defined above, inside a larger tubular domain T_{b}(s,ρ,θ). Thus our braid is defined by a set of n curves r_{i}(s) = T_{b}(s,ρ_{i}(s),θ_{i}(s)),i = 1,...n (an example is shown in Fig. 4a). These curves r_{i} can then be each individually identified as the axis curve r above and we can create a tubular field surrounding each r_{i}. It is this mechanism by which we add further complexity to the set of fields that we can generate. Various algorithms described in Berger (2001) or Bangert et al. (2002) could be used to construct the r_{i} for a particular braid; we give an explicit example in Sect. 3. We first construct this braid on a cylinder, providing functions ρ_{i}(s) and θ_{i}(s) in cylindrical coordinates for each r_{i}. We then use these functions as the tubular coordinates ρ_{i}(s), θ_{i}(s) in the overall domain T_{b}. A “stirring” boundary motion which could be used to create this braided field in an initially unbraided tube (with u_{3} = 0) is shown in Fig. 4b (a similar idea was used in Berger & AsgariTarghi 2009). There is a formal mathematical relationship between the stirring of a fluid in a twodimensional environment and the mathematical braids produced by creating threedimensional space curves which follow the time dependent stirring geometry. The pigtail pattern shown in Fig. 4b is the most complex possible stirring motion (Boyland et al. 2000).
We can additionally define the individual strands of the braid to have a net twist (a non trivial twisting function ψ_{i}(s)). This kind of field might mimic the effect of both static (unbraided) and moving vortex motions at the photospheric boundary. This idea will be explored further in part II of this study.
2.5. Practical matters
2.5.1. Interpolating the field
We intend to use these tubular fields (hereafter denoted B_{t}) in numerical simulations on a Cartesian grid. This means we need to interpolate our field from the tubular coordinate system (s,ρ,θ) onto a Cartesian grid. To do so we first take a Cartesian point C_{p} = (x,y,z) and check if it is in the tube by calculating its minimum distance to the tube axis. That is we find the solution s_{Cp} to (14)that minimises   (C_{p} − r(s)  . This can be done using any standard algorithm for root finding. If this minimum has C_{p} − r(s_{CP}) >R(s), then C_{P} does not lie inside the tube and we set B_{t} = (0,0,0). If C_{p} − r(s_{CP}) ≤ R it does lie inside the tube and we calculate its (s,r,θ) coordinates. We already have the s coordinate s_{CP} from our root finding algorithm; we get the radial coordinate from ρR = C_{p} − r(s_{CP}), finally we can solve (15)to obtain the principal value of θ_{0}. Once these coordinates have been obtained we know which field line f [ ρ_{0},θ_{0} ] the point C_{p} lies on and we can integrate Eq. (11) to obtain the components of B.
2.5.2. Ensuring the field is divergencefree
Whilst B is divergence free by construction, we found in practice that the interpolation error led to nonzero divergence when transferred to a finite Cartesian grid. Since maintaining the divergencefree condition to a highdegree of accuracy is important in MHD simulations (Tóth 2000), we added an additional divergence cleaning step. We take the field B produced by the interpolation routine described in Sect. 2.5.1, and find a vector potential A such that B = ∇ × A (this procedure is performed on a Cartesian grid) . This vector potential can then be curled numerically to give a field B that is numerically divergencefree to machine precision. The choice of A used in this procedure has an effect on how faithfully the magnetic field direction is maintained. We found the BiotSavart gauge of Prior & Yeates (2014) to work well for this task, as it ensured averaging of the errors over space rather than concentrating them along certain lines.
We found the process of transferring B_{t} to the Cartesian grid to be much more successful when the field strength function φ(ρ_{0},θ_{0}) tends smoothly to zero as ρ → R. For the examples in Sects. 3 and 4, we use the azimuthallysymmetric function φ_{0}(ρ_{0},θ_{0}) = 0.5 φ_{c}(cos(πρ_{0}) + 1), where φ_{c} is a constant. However, we should note that, other than differentiability, there are no strict conditions on the function’s radial decay; this function was merely a convenient choice. It avoids a discontinuous B at the edge of the tubular domain, and associated problems with interpolation. One could create distributions with far sharper radial decay, but care would be needed to ensure that the Cartesian grid used could resolve it sufficiently.
3. Example: twisted and braided flux ropes
This section describes the creation of some specific magnetic configurations designed to represent different possible internal structures (twisted and braided) of the same overall magnetic flux rope in the low corona. These fields are designed to be initial conditions for the comparative resistiveMHD simulations performed in Sect. 4 on each field.
3.1. Background field
Fig. 5 Field lines of twisted a) and braided b) flux ropes (primary colours) and the background linear force free field (light colouring). Field lines can be seen to peel of the ropes as a result of the influence of the background field B_{b}. These are the initial conditions for the simulations in Sect. 4. 

Open with DEXTER 
In the coronal region flux ropes do not tend to exist in isolation. Typically the rope is embedded in and connected to a surrounding field. In this initial study, we embed our tubular field B_{t} in in a forcefree background B_{b}, ∇ × B_{b} = αB_{b}, an appropriate choice for the low corona where the magnetic pressure dominates the plasma pressure (e.g. Priest 2003). For simplicity here, we choose a linear forcefree field in the volume z> 0 specified by choosing a constant α and a boundary flux distribution for B_{z0} on z = 0. We define a quadrupolar distribution (found to be the distribution most commonly underlying solar filaments; Mackay et al. 2008), which we specify as a sum of Gaussians (16)where two of the A_{i} are positive and two negative. This produces a set of 2D Gaussian peaks centred at the coordinates (x_{i},y_{i}). In this particular study we use the parameters A_{1 / 2} = 1,A_{3 / 4} = −1,B_{i} = 5,i = 1,2,3,4 and positions (x_{1},y_{1}) = (−0.68π, − 0.3π), (x_{2},y_{2}) = (−0.2π,0.5π), (x_{3},y_{3}) = (0.3π, − 0.5π) and (x_{4},y_{4}) = (0.7π,0.4π), with the distribution restricted to a square −π ≤ x ≤ π, − π ≤ y ≤ π and an alpha value α = 0.5. The parameter A_{i} is chosen so that the maximum field strength of the background field is 1. This will allow us to set the relative strength of the embedded flux rope against 1. The parameter B_{i} chosen so that the flux sources are sufficiently narrow so as not to overlap. This flux distribution can be seen in Fig. 2c along with its neutral line which has has a zshaped morphology. We then generate the field on the domain z> 0 using the method detailed in Prior & Berger (2012). Field lines of the initial composite field B_{b} + B_{t}, for twisted and braided tubes, after interpolating to the Cartesian grid and divergence cleaning, are shown in Fig. 5. These examples are the initial fields whose subsequent evolution is simulated in Sect. 4. Field lines can be seen to peel off the tube structure, an effect which is more prevalent in the braided field. This is not a numerical artefact but is a result of the superposition of the two fields.
3.2. Tube morphology
To find an appropriate shape for the centreline r_{t} of the overall tube T_{b}, we use the general observation that solar filaments are located above neutral lines of the photospheric flux distribution (Mackay et al. 2010). The Zshaped neutral line, B_{z0}(x,y) = 0, of our background boundary flux distribution B_{b} forms a twodimensional curve (n_{x}(t),n_{y}(t)) t ∈ [ −1,1 ]. This can be seen in Fig. 2c. We define a threedimensional curve s to be (17)where t ∈ [ −f_{p1},f_{p2} ], 0 ≤ f_{p1}<f_{p2} ≤ 1 and h is the maximum height (zcoordinate) of r_{t}. Since the neutral line straddles the whole domain (Fig. 2c), choosing f_{p1} = 0 and f_{p2} = 1 would mean that the curve r_{t} had its footpoints on the box boundaries; this would not be helpful as part of the tube would then lie outside the domain. We found values h = 3, f_{p1} = 0.15 and f_{p2} = 0.85 gave a reasonable flux rope as depicted in Fig. 2c. This forms the axis curve of the overall tubular manifold T_{b}, within which we can either create a twisted field or a braided field (with several subtubes).
Care has to be taken as the parametrisation t will not be the arclength of the curve. Whilst arclength is a convenient parametrisation for analytic work, it is generally hard to use in practice. This means that all arclength derivatives detailed in Sect. 2 have to be calculated using (18)
3.2.1. Practical matters
The neutral lines in B_{z0} distributions of the type in Eq. (16) can be particularly sharp, while our field generation technique requires curves are at least differentiable. To deal with this, we smoothed the initial neutral line curve (n_{x},n_{y}) using a moving average and used a spline fitting routine in Mathematica 10 to obtain an analytic curve for the centreline.
Another feature of the 2D neutral line curve is that the sharp transition in its midsection can lead to curves r_{t} with high curvature values around their apex. The fact that the tubular coordinate system becomes illdefined when the product κR ≥ 1 would lead to a failure in the field generating algorithm of Sect. 2 (the field strength diverges). Thus we would have to limit the tube to be significantly thin (make R small). Alternatively we found that by either shuffling the positions of the Gaussian sources (changing the (x_{i},y_{i}) in Eq. (16)) or applying filtering techniques could reduce the curvature of r_{t} without significantly affecting the basic morphology of the tube. This allowed for thicker tubes allowing us to include more interesting topological structure within the tube.
We also found that if the product κR became near to 1 then the resolution used to discretize the integral in Eq. (11) had to be significantly high (> 50 000) in order to lead to a faithful recreation of the tubular field B_{t} after the cleaning process. This issue was overcome by using an adaptive integration grid based on the value of the product κR.
3.3. Internal structure
In this study we are not addressing the challenging process of generating braided structures by some timedependent physical process. Instead we compare two different internal structures within the same tube of radius R = 5 whose axis r_{t} was defined in Sect. 3.2 (Fig. 2b). The first is a twisted field with axis r_{t}, R = 5, and , where L is the arclength of r_{t} (so that is the total twisting). This is compared with a braided structure, made up of three untwisted (ψ_{i} = 0) internal strands of radius R_{i} = 0.2 which are linked in a pigtail braid (19)as illustrated in Fig. 2a. The curves are contorted to fit in the overall tube with axis r_{t}.
Finally, the distribution φ_{0}(ρ_{0},θ_{0}) must be specified, which, through Eq. (11), gives the magnetic flux on the boundary z = 0. The boundary flux of the background field B_{b} is given by Eq. (16), it has a maximum value of 1, which is also the maximum value of the whole field since it decays in z. We have chosen the azimuthallysymmetric form φ_{0}(ρ_{0},θ_{0}) = 0.5 φ_{c}(cos(πρ_{0}) + 1), where φ_{c} is a constant determining the maximum value of φ_{0}. This controls the relative strength of the tubular field with respect to the background field. In this study we choose φ_{c} = 1 to contrast and compare the behaviour of a flux tube of roughly equal strength to the background field, we also report briefly on results of φ_{c} = 3 and φ_{c} = 6, flux tubes significantly stronger than their background .
4. Numerical simulation
The superposition B_{t} + B_{b}, as described in Sect. 3, defines an initial configuration. This field will not generally be in equilibrium, unlike the toroidal twisted fields used in Fan & Gibson (2007), Kliem et al. (2004, 2010, 2012), Leake et al. (2014) for which an equilibrium within the background dipole field can be found (Titov & Démoulin 1999). The relative complexity of our background field, tube morphology and internal structure, make finding an analytical equilibrium unfeasible. That said, we note that there is evidence that braiding is a dynamic phenomenon (van Ballegooijen et al. 2014), and it is reasonable that we might see complex internal structures built up relatively quickly in the corona. In this paper, we investigate whether our twisted and braided configurations are physically feasible by studying their structural stability in resistiveMHD simulations. To keep things simple, we also avoid boundary motions. This allows us to focus on the questions of how the processes of ideal evolution and resistive reconnection alter the field, and whether an equilibrium that retains the magnetic structure can be found.
4.1. MHD equations
We evolve the initial field B_{t} + B_{b} by using the Lare3D Lagrangianremap code (Arber et al. 2001) to solve the resistiveMHD equations on a Cartesian box {− π ≤ x,y ≤ π,0 ≤ z ≤ 3π}, at resolution 264 × 264 × 396. The initial tube’s maximum height is less than 4 so this gives room for any reasonable vertical expansion. For simplicity, we apply line tied boundary conditions on all boundaries. This means that we must stop the simulations before reflecting waves from the boundary come to dominate. The code solves the following equations Here ρ_{d} is the mass density, v the plasma velocity, B the magnetic field, j the current density, p the plasma pressure, σ the stress tensor, ϵ the specific internal energy, η the resistivity, e the strain tensor, and γ = 5 / 3. The viscous term ∇·σ in Eq. (21) includes no background viscosity, but only shock viscosity to prevent unphysical oscillations and approximate the jump in entropy across shocks. The shock viscosity takes the form given in Bareford et al. (2013), and we use the same successful parameter values ν_{1} = 0.1, ν_{2} = 0.5. There is a corresponding heating term ϵ:σ in Eq. (23). We initially set ν = 1 and ϵ = 0.01 in nondimensional units. In these units, one unit of time is equal to the time taken by an Alfvén wave with B = ρ = 1 to move a unit distance in the box. The simulations presented here use a uniform resistivity of 5 × 10^{4}, a value found from experience to be just above the numerical resistivity at this resolution. This system and parameter set is the same as that used in Yeates et al. (2015). The initial velocity is assumed to vanish everywhere so that any motion will result from the field attempting to relax to equilibrium. The diffusion time of the rope/braid, the square of its radius divided by the resistivity, is approximately 500. In practice our simulations ended a long time before this, so any significant global evolution of the rope is likely due to dynamical behaviour.
4.2. Resulting relaxation
The results presented are for φ_{c} = 1, so that B_{t} is roughly of equal strength to its background B_{t}. The simulations were run until boundary effects became too significant and the timestep increased dramatically. In both cases we shall see the changing energy had begun to plateau at this point and in both cases had reached a forcefree state. The twisted simulation ran for just over twice as long (t ∈ [ 0,9.8 ]) as the braided simulation (t ∈ [ 0,4.2 ]).
Fig. 6 Snapshots of the twisted core at t = 0a), t = 4.9b) and t = 9.8c). The field is seen to expand and flatten slightly. There also appears to be a counterclockwise (from above) rotation of field lines at the apex of the rope. 

Open with DEXTER 
Fig. 7 Snapshots of the braided cores at t = 0a), t = 2.1b) and t = 4.2c). The curves are coloured based on their start points as belonging to one of the three initial braid footpoints. The interference of the background field is clear as an increasing number of field lines which would have originally belonged to the braid can be seen to have reconnected with the background field. 

Open with DEXTER 
Fig. 8 Snapshots of contours of constant current density   j   = 0.5 for the twisted rope at t = 0a), t = 4.9b) and t = 9.8c). In a) we see the flux rope’s initial morphology, in b) the outer layer of current expanding to expose the rope’s remaining core, and in c) that the core is clear and forms a tubular shape similar in radius to the initial tube but with a lower apex. 

Open with DEXTER 
Fig. 9 Snapshots of contours of constant current density   j   = 0.5 for the braided flux rope at t = 0a), t = 2.1b) and t = 4.2c). In a) the braid structure is clear. In b) the outer layer of current expands exposing a clear core seen in c); by comparison to the twisted case it is significantly fractured. 

Open with DEXTER 
Fig. 10 Snapshots of contours of the Lorentz force magnitude for the twisted rope, at t = 0a), t = 4.9b) and t = 9.8c). In a) the background is forcefree whilst the tube has an associated Lorentz force. In b) we see the force beginning to dissipate.c) the force has all but disappeared except on the boundary. 

Open with DEXTER 
Fig. 11 Snapshots of contours of the Lorentz force magnitude for the braided rope, at t = 0a), t = 0.49b) and t = 4.2c). In a) the background is forcefree whilst the braid structure has an associated Lorentz force. The force has reduced significantly in b), leaving a fracture structure. c) The force has all but disappeared except on the boundary. 

Open with DEXTER 
Fig. 12 Plots of the magnetic, kinetic, internal and (normalised) total energy as a function of time for both the braided and twisted simulations. The, magnetic, kinetic and internal energy plots show the difference from their values at t = 0. In both cases the magnetic energy is converted into kinetic and internal energy. Two phases can be identified. The first is a rapid change as the magnetic is converted into kinetic energy and internal energy increase (t ∈ [ 0,0.1 ] for the braid and t ∈ [ 0,35 ] for the twisted field). In the second phase the kinetic energy levels off and then falls steadily, whilst the internal energy increases gradually. The loss of magnetic energy of the twisted field is larger than the braided field. 

Open with DEXTER 
Fig. 13 Viscous (unbroken lines) and ohmic (dashed lines) heating during the simulations. In a) we see, similar to the energy plots, the braided quantities increase at a quicker rate than the twisted quantities. In both cases the viscous dissipation is larger than the ohmic heating by roughly a factor of 5. b) Shows the same quantities with the braided quantities scaled by the ratio of the initial twisted to braided energy (1.65), at the point where the simulations end the viscous dissipation is roughly the same. 

Open with DEXTER 
Fig. 14 Snapshots of the vertically integrated current for the twisted field at t = 0a), t = 4.9b) and t = 9.8c). In a) we see the Zsigmoidal shape of the initial tube. b) Shows that the tube’s outer current layer has expanded. In c) the outer shell has further expanded leaving an inner rope structure which is far less kinked than in a). The rotation of the tube core from a) to c) is counterclockwise. 

Open with DEXTER 
Fig. 15 Snapshots of the vertically integrated current at t = 0a), t = 2.1b) and t = 4.2c). In a) we see the braided structure projected. b) Shows that the tube’s outer current layer has expanded. In c) the outer shell has further expanded leaving an inner shell with a significantly fractured set of current peaks. 

Open with DEXTER 
The evolution of the field lines composing the cores of the twisted and braided tubes is shown in Figs. 6 and 7 respectively. The field lines are drawn by choosing a set of curves anchored at locations in the initial footpoints of the twisted tube and in the individual strands of the braid. The twisted curves clearly expand during the evolution. It can also be seen they rotate (anticlockwise as seen from above). A number of field lines peel away from the structure (there is a significant change in their end points), and the increasing number of such field lines suggests that reconnection between the flux rope and the background field is occurring. The evolution of the braided field lines of Fig. 7 is less clear, other than a general radial expansion of the tube and a decrease in the height of the core’s apex. Again, there also appears to be substantial reconnection with the background field with an increasing number of field lines which begin at the footpoint of the original braid being directed away from the core of the rope/braid.
Contour plots of the current density (  j  ), shown in Figs. 8 (twisted) and 9 (braided), give more clarity. We choose a fixed level of 0.5, half of the maximum in the background field. The braided structure can be seen with much greater clarity than for the field line plots. In both cases the outer layer of current (between the tube and the surrounding background field) expands to leave a core of current. The twisted core is itself a tube, whilst the braided curve has a fractured structure. This is consistent with the results of WilmotSmith et al. (2011), who found fields with complex braiding to develop a more diffuse smallscale current structure. In our simulations the outer expanding front continues until it reaches the linetied boundary of the domain where it is reflected back and the simulation is stopped.
Plots of the Lorentz force magnitude (  j × B  ) are shown in Figs. 10 and 11; both fields relax such that the Lorentz force is significantly reduced (except at the line tied boundary). For the t = 9.8 braid configuration the maximum Lorentz force (away from the boundary) is 0.69 and the mean value 0.006, by comparison at t = 0 the maximum value is 20.6254 and the mean is 0.093. For the t = 4.2 braid configuration the maximum Lorentz force (away from the boundary) is 0.47 and the mean value 0.00195, by comparison at t = 0 the maximum value is 44.96 and the mean is 0.040. Also shown in Fig. 12 are plots of the change in magnetic, kinetic, and internal energy and their sum. The magnetic energy is converted into kinetic and internal energy, whilst the total sum is well conserved. The change in each quantity is larger for the twisted field, though we note that the initial magnetic energy of the twisted field over the braided field has a ratio ≈ 1.65, and taking this ratio into account the changes are comparable. Initially the magnetic energy is converted into mostly kinetic energy, after which this kinetic energy decreases steadily at a slow rate. With our linetied boundaries and in the absence of a uniform viscosity, small scale motions persist in our simulation. Comparing the braided and twisted plots we see that the change in kinetic energy occurs at a significantly faster rate for the braided field, in which the more complex field structure generates dissipative losses at a faster rate (Fig. 13a). After this, the internal energy increases significantly before levelling off, as a result of both ohmic heating and viscous heating (Fig. 13a); in common with previous simulations by Bareford et al. (2013). There is more viscous heating than ohmic heating, also the values of both quantities are larger for the twisted field. However, when the difference in initial energy is accounted for, we find both types of dissipation to be roughly similar for the two fields (Fig. 13b). So in both cases the fields have relaxed to fairly stable structures. That the initial kinetically driven change is insensitive to the resistivity was suggested by rerunning the simulation with a resistivity 5 × 10^{3} (an order of magnitude larger); barring a small increase in the amount of ohmic dissipation the plots remained almost unchanged.
Figures 14 and 15 show the distribution of a proxy of the emission which might be viewed by lineof sight imaging in extremeultraviolet or Xray wavelengths (Cheung & DeRosa 2012). We average the square of the current density (j·j) along all field lines in the domain, then integrate this quantity vertically to mimic the line of sight view. The initial sigmoidal shape can be seen in (a) for both the twisted and braided fields. The field lines which peel off the ropes are also present, these field lines have significantly contorted geometry and high current. In (b), for both fields, we see the expansion of the outer current layer. In Figs. 14b and c the outline of the emergent core flux rope seen in Fig. 8c is again visible; relative to the initial sigmoid (a) it is far less kinked. The indicated rotation of the tube is of the same chirality as the rotation of the twisted field lines seen in Fig. 6. By comparison the final braided shape in 15c is much more complex and diffuse.
To summarise: both twisted and braided fields relax to a stable state for which the tubular structure is still clearly present. In both cases there is an initial current expansion to leave a core structure with sigmoidal morphology. The final braided state is much more complex than the twisted state in that its current structure is diffuse and composed of small scale current islands, whilst the twisted case retains a more regular tubular morphology. The integrated current plots (Figs. 14 and 15) are of particular interest as they raise the possibility that one might be able to recognise the signature of a more complex internal tube morphology through lineofsight observation.
Whilst we can clearly observe large scale rotation of the twisted field, it is not clearly apparent what is happening internally in the (initially) braided tube. To investigate the internal topology, we consider an average local twisting rate along field lines.
Fig. 16 Snapshots of the local twisting distribution L_{f}(x,y) on a subdomain of the boundary plane which contains the initial tube footpoint, at the beginning a), middle b) and end c) of the twisted field simulation. The sign is largely negative except at the tube edge. 

Open with DEXTER 
Fig. 17 Snapshots of the local twisting distribution L_{f}(x,y) on a subdomain of the boundary plane which contains the initial tube footpoint, at the beginning a), middle b) and end c) of the braided field simulation. Both positive and negative values are present and the magnitude increases in b) and c). 

Open with DEXTER 
For a field line f(l) of arclength L_{f} whose footpoint coordinates are (x_{f},y_{f}), we define the integrated quantity (26)which represents the mean rotation of the local field lines around f(l). For a linear forcefree field j·B = αB·B and L_{f} is just the linear forcefree parameter α, which would be constant throughout the domain. This quantity was calculated for relaxing braided and twisted cylindrical fields in WilmotSmith et al. (2011), Yeates et al. (2015). These authors found that a twisted field relaxed to a forcefree state with a single sign of α within the tube, although the surrounding background field prevented it from reaching the spatially constant α predicted by the Taylor relaxation hypothesis. By contrast, they found braided tubes to relax to form two forcefree flux tubes with α values of equal and opposite sign, entirely contrary to the Taylor hypothesis. Here we consider the distribution L_{f}(x,y) of local twisting for each field line anchored at points (x,y). Since we are really interested in the flux rope core we concentrate on the domain −1.320365 ≤ x ≤ −0.010365, −2.3702 ≤ y ≤ −1.1702, which contains the flux rope’s starting footpoint.
Evolution of L_{f}(x,y) for the twisted field is shown in Fig. 16. It has one dominant (negative) sign consistent with the sign of twisting initially imposed in the tube. The magnitude of twisting decreases with time, consistent with the loss of twisting through rotation induced writhing of the tube, which was observed in Figs. 8 and 14. The twistwrithe decomposition of helicity (Berger & Prior 2006) tells us to expect a loss of twisting upon global rotation of the tube. The other feature is the region of positive L_{f} values on the tube edge. This corresponds to the set of field lines which connect to the background field, as seen in Fig. 6. The increase in size of this positive L_{f} region is consistent with the reconnection of the tube and the background field observed in Figs. 6b and c. Outside the tube there is a constant value of L_{f} = 0.5 from the background linear forcefree field.
The distribution of L_{f}(x,y) for the braided field is shown in Fig. 17. At all three times both signs of L_{f} are observed and present with significant magnitudes. Comparing (b) and (c) to (a) we see an increase in the area of which the magnitude of L(x,y) is greater than 2 (for both signs), indicating a move towards more inhomogeneous and extreme twisting geometries; this contrasts to the general decrease in magnitude of L_{f} in the tube’s interior, due to writhing of the twisted field. In (b) and (c) there are several observable regions of a particular sign of L_{f}, in particular there are two large regions with a dominant sign L_{f}, one positive one negative, though neither are homogeneous. Whilst it is not the clear splitting into two tubes found in WilmotSmith et al. (2011), Yeates et al. (2015), it is indicative that the tube has not relaxed to a linear forcefree state. Unlike the twisted field it is not clear what the effect of reconnection with the background field has had on the distribution.
We remark that the same set of simulations and analysis were performed for fields with φ_{c} = 3 and φ_{c} = 6, i.e. flux ropes of strength significantly greater than the background field B_{b}. The relative weakness of the background field meant that the field lines of both tubes expanded to greater radii. However, we did not observe any significant kinking reminiscent of the type of eruptive event seen in the toroidal twisted tube simulations (Fan & Gibson 2007; Kliem et al. 2004, 2010, 2012; Leake et al. 2014). The current structure was qualitatively similar to the results reported above (their magnitude aside) and in all cases the field relaxed to a forcefree state with the energy plots levelling off as in Fig. 12. Finally, accounting for magnitude, the L_{f} distributions were very similar. We do not present that data here as it does not add significantly to the conclusions. We remark that it is possible the line tied boundary conditions of the simulation could have affected these simulations more significantly than the φ_{c} = 1 case. Performing the simulations in a larger box might have allowed the flux tube to displace the background field in a more significant manner, allowing for more dramatic changes in morphology. This possible effect would be more marked for the higher φ_{c} simulations as they expand to a greater degree. In this study we wanted a reasonable resolution for which we could perform a number of simulations in order obtain an understanding of how such sigmoidal braided structures would respond. Increasing the domain size whilst maintaining the same resolution would lead to a significant increase in simulation time and this task is left for a future study.
5. Conclusions
A new technique for generating magnetic flux tubes of arbitrary axial geometry, varying radius and internal geometry, was introduced in Sect. 2. A procedure for embedding this field over the neutral line of a boundary flux distribution was also established (Sect. 3). The ultimate aim for introducing this technique is to explore the relationship between the internal structure and evolving global morphology of coronal flux ropes with realistic geometry. It is hoped that this will provide insight to complement the currently available observational data for the corona.
Sect. 4 reports the results of initial simulations of a sigmoidal flux tube embedded in a linear forcefree field, generated by a quadrupolar distribution. Tubes were created containing both twisted and braided internal fields. In both cases the tube was seen to relax to a state with significantly reduced Lorentz force, whilst its magnetic, kinetic and internal energies began to plateau, indicating the tube had stabilised itself. The current structure of the two fields differed significantly, with the braided field developing a more complex structure. The current density projected onto the photospheric boundary differed significantly for the two tubes in their relaxed state, with the twisted field exhibiting the usual sigmoidal morphology, whilst the braided field yielded a more complex small scale distribution. This suggests a possible means of observing evidence of flux tubes with significantly tangled internal structures.
The main advantage of this procedure is its flexibility. Apart from a certain level of smoothness, the axial geometry is not constrained. The control over the internal structure is also significant as one can use either multiple subtubes or overlapping fields to generate a significant variety of internal geometries. The fields produced are divergence free to a high degree of accuracy. By allowing for varying tube radii one might also be able to account for the rapidly varying tube radii seen in the photospherecorona transition layer. There are some technical issues that must be considered when creating the fields. For axial curves with significant curvature the tube radius must be limited in order to ensure the tubular coordinate system Eq. (2) is well defined (Sect. 2.1.2). Also, more complex tubular structures will require either large or adaptive grids in order to be resolved. One potential drawback is the fact that the composite field of the braided flux rope B_{t} and the background field B_{b} will not generally be in equilibrium. There is no way of knowing that the structure created will be stable. This may not be too much of a problem as the corona is a dynamic environment, so starting at a static equilibrium is not essential in principle. That said, there are significant advantages to beginning at an equilibrium state, such as the ability to perform parameter studies which can highlight the effect of a particular parameter, or change in state, in triggering the various MHD instabilities which flux ropes can undergo (see Schmieder et al. 2013, for a review). If a stable forcefree state is required one could apply an ideal relaxation (e.g. van Ballegooijen 2004; Pontin et al. 2011). The example simulations of Sect. 4 suggest that placing the tube over the photospheric neutral line could be a means of producing a forcefree tube with complex internal topology and specific sigmoidal structure. A limitation of the current method for specifying the internal field is the fact that it is not easy to include singular points in the field.
The simulations provide some tentative evidence that tangled flux ropes could persist for significant time periods in the corona. That said, there are a number of questions raised by the simulations which deserve attention in future studies. Would the flux tubes have remained stable if they had been contained in a larger domain which may have given more freedom for the tube to displace the containing background field? If the reflecting waves caused by the line tied boundary conditions had not forced the simulation to be terminated, could reconnection have eventually rendered the fields unstable, or alternatively reach a genuine linear force free state? What would have happened with background fields B_{b} of differing α, a potential field with the same boundary flux or a nonlinear forcefree field, and how important is the distribution of the boundary flux sources? Can some distributions stabilise the flux rope where others cannot? Finally, we chose a somewhat idealised initial pigtail braid distribution. There is a huge range of different internal distributions which differ significantly from this pigtail in both braided and twisted field structure; are there certain classes of internal topologies which behave in a similar manner?
Part II of this study will tackle some of these questions through a multitude of simulations in which the parameters mentioned in the previous paragraph are varied systematically. An additional study will recreate erupting twisted toroidal flux rope simulations and compare the same initial system with an internally braided flux rope replacing the twisted one. In this case we know that there is a guaranteed loss of stability of the twisted field, so one could test whether a braided distribution could stabilise the flux rope when a twisted geometry does not. Beyond this, the next step would be to include boundary motions. A good number of the studies summarised in WilmotSmith (2015) take initially straight fields in Cartesian box and attempt to braid them with small scale shearing motions. Using the techniques detailed here one could perform the a similar study in a sigmoidal flux rope whose initial internal structure is undistorted (a field specified by the map T with u_{3} = 0).
Acknowledgments
Chris Prior would like to acknowledge the support of an Addison Wheeler Fellowship. A.R.Y. thanks STFC for financial support. This work made use of the Hamilton supercomputer facility provided by the HPC Service of Durham University.
References
 Antman, S. S. 2005, Nonlinear problems of elasticity, Applied Mathematical Sciences (New York: Springer), vol. 107 [Google Scholar]
 Arber, T., Longbottom, A., Gerrard, C., & Milne, A. 2001, J. Comput. Phys., 171, 151 [NASA ADS] [CrossRef] [Google Scholar]
 Archontis, V., & Török, T. 2008, A&A, 492, L35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Aulanier, G., Török, T., Démoulin, P., & DeLuca, E. E. 2010, ApJ, 708, 314 [NASA ADS] [CrossRef] [Google Scholar]
 Bangert, P. D., Berger, M. A., & Prandi, R. 2002, J. Phys. A: Mathematical and General, 35, 43 [NASA ADS] [CrossRef] [Google Scholar]
 Bareford, M., Hood, A., & Browning, P. 2013, A&A, 550, A40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Berger, M. A. 2001, J. Phys. A: Mathematical and General, 34, 1363 [NASA ADS] [CrossRef] [Google Scholar]
 Berger, M. A., & AsgariTarghi, M. 2009, ApJ, 705, 347 [NASA ADS] [CrossRef] [Google Scholar]
 Berger, M. A., & Prior, C. 2006, J. Phys. A: Mathematical and General, 39, 8321 [NASA ADS] [CrossRef] [Google Scholar]
 Bishop, R. L. 1975, Amer. Math. Month., 82, 246 [CrossRef] [Google Scholar]
 Boyland, P. L., Aref, H., & Stremler, M. A. 2000, J. Fluid Mech., 403, 277 [NASA ADS] [CrossRef] [Google Scholar]
 Cheng, X., Zhang, J., Liu, Y., & Ding, M. 2011, ApJ, 732, L25 [NASA ADS] [CrossRef] [Google Scholar]
 Cheung, M. C., & DeRosa, M. L. 2012, ApJ, 757, 147 [NASA ADS] [CrossRef] [Google Scholar]
 Cirtain, J. W., Golub, L., Winebarger, A. R., et al. 2013, Nature, 493, 501 [NASA ADS] [CrossRef] [Google Scholar]
 Craig, I., & Sneyd, A. 2005, Sol. Phys., 232, 41 [NASA ADS] [CrossRef] [Google Scholar]
 Fan, Y. 2009, ApJ, 697, 1529 [NASA ADS] [CrossRef] [Google Scholar]
 Fan, Y., & Gibson, S. 2007, ApJ, 668, 1232 [NASA ADS] [CrossRef] [Google Scholar]
 Frankel, T. 2011, The geometry of physics: an introduction (Cambridge University Press) [Google Scholar]
 Galsgaard, K., & Nordlund, Å. 1996, J. Geophys. Res.: Space Phys., 101, 13445 [NASA ADS] [CrossRef] [Google Scholar]
 Gibson, S., Fan, Y., Török, T., & Kliem, B. 2006a, Space Sci. Rev., 124, 131 [NASA ADS] [CrossRef] [Google Scholar]
 Gibson, S. E., Foster, D., Burkepile, J., de Toma, G., & Stanger, A. 2006b, ApJ, 641, 590 [NASA ADS] [CrossRef] [Google Scholar]
 Green, L., Kliem, B., Török, T., van DrielGesztelyi, L., & Attrill, G. 2007, Sol. Phys., 246, 365 [NASA ADS] [CrossRef] [Google Scholar]
 Inoue, S., Hayashi, K., Magara, T., Choe, G., & Park, Y. 2015, ApJ, 803, 73 [NASA ADS] [CrossRef] [Google Scholar]
 Janse, Å., Low, B., & Parker, E. 2010, Phys. Plasmas, 17, 092901 [NASA ADS] [CrossRef] [Google Scholar]
 Kliem, B., Titov, V., & Török, T. 2004, A&A, 413, L23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kliem, B., Linton, M., Török, T., & Karlickỳ, M. 2010, Sol. Phys., 266, 91 [NASA ADS] [CrossRef] [Google Scholar]
 Kliem, B., Török, T., & Thompson, W. 2012, Sol. Phys., 281, 137 [NASA ADS] [Google Scholar]
 Leake, J. E., Linton, M. G., & Antiochos, S. K. 2014, ApJ, 787, 46 [NASA ADS] [CrossRef] [Google Scholar]
 Linton, M., Longcope, D., & Fisher, G. 1996, ApJ, 469, 954 [NASA ADS] [CrossRef] [Google Scholar]
 Longcope, D., & Beveridge, C. 2007, ApJ, 669, 621 [NASA ADS] [CrossRef] [Google Scholar]
 Mackay, D. H., Gaizauskas, V., & Yeates, A. R. 2008, Sol. Phys., 248, 51 [NASA ADS] [CrossRef] [Google Scholar]
 Mackay, D. H., Karpen, J. T., Ballester, J. L., Schmieder, B., & Aulanier, G. 2010, Space Sci. Rev., 151, 333 [NASA ADS] [CrossRef] [Google Scholar]
 MacTaggart, D. 2011, A&A, 531, A108 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Manchester IV, W., Gombosi, T., DeZeeuw, D., & Fan, Y. 2004, ApJ, 610, 588 [NASA ADS] [CrossRef] [Google Scholar]
 Mandrini, C., Pohjolainen, S., Dasso, S., et al. 2005, A&A, 434, 725 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ng, C., Lin, L., & Bhattacharjee, A. 2012, ApJ, 747, 109 [NASA ADS] [CrossRef] [Google Scholar]
 Okamoto, T. J., Tsuneta, S., Lites, B. W., et al. 2008, ApJ, 673, L215 [NASA ADS] [CrossRef] [Google Scholar]
 Parker, E. 1972, ApJ, 174, 499 [NASA ADS] [CrossRef] [Google Scholar]
 Pontin, D., WilmotSmith, A., Hornig, G., & Galsgaard, K. 2011, A&A, 525, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Priest, E. R. 2003, Solar magnetohydrodynamics (Springer), 1 [Google Scholar]
 Prior, C., & Berger, M. 2012, Sol. Phys., 278, 323 [NASA ADS] [CrossRef] [Google Scholar]
 Prior, C., & Yeates, A. 2014, ApJ, 787, 100 [NASA ADS] [CrossRef] [Google Scholar]
 Rappazzo, A., & Parker, E. 2013, ApJ, 773, L2 [NASA ADS] [CrossRef] [Google Scholar]
 Reale, F. 2010, Liv. Rev. Sol. Phys., 7, 5 [Google Scholar]
 Schmieder, B., Démoulin, P., & Aulanier, G. 2013, Adv. Space Res., 51, 1967 [NASA ADS] [CrossRef] [Google Scholar]
 Su, Y., van Ballegooijen, A., McCauley, P., et al. 2015, ApJ, 807, 144 [NASA ADS] [CrossRef] [Google Scholar]
 Titov, V., & Démoulin, P. 1999, A&A, 351, 707 [NASA ADS] [Google Scholar]
 Tóth, G. 2000, J. Comput. Phys., 161, 605 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 van Ballegooijen, A. A. 2004, ApJ, 612, 519 [NASA ADS] [CrossRef] [Google Scholar]
 van Ballegooijen, A. A., & Martens, P. C. H. 1989, ApJ, 343, 971 [NASA ADS] [CrossRef] [Google Scholar]
 van Ballegooijen, A., AsgariTarghi, M., & Berger, M. 2014, ApJ, 787, 87 [NASA ADS] [CrossRef] [Google Scholar]
 WilmotSmith, A. 2015, Philos. Trans. Roy. Soc. London A: Mathematical, Physical and Engineering Sciences, 373, 20140265 [NASA ADS] [CrossRef] [Google Scholar]
 WilmotSmith, A., Pontin, D., Yeates, A., & Hornig, G. 2011, A&A, 536, A67 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wyper, P., & Pontin, D. 2014, Phys. Plasmas, 21, 102102 [NASA ADS] [CrossRef] [Google Scholar]
 Yeates, A., Hornig, G., & WilmotSmith, A. 2010, Phys. Rev. Lett., 105, 085002 [NASA ADS] [CrossRef] [Google Scholar]
 Yeates, A., Bianchi, F., Welsch, B., & Bushby, P. 2014, A&A, 564, A131 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Yeates, A., Russell, A., & Hornig, G. 2015, Proc. Roy. Soc. A, 471, 20150012 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Explicit Expression for ∇·N
In order to obtain the magnetic field N from Eq. (11) we need the divergence of N. Since we can express this field generally in tubular coordinates (ρ,θ,s), it would be easier to calculate the divergence in this coordinate system. In an arbitrary coordinate system (x_{1},x_{2},x_{3}) the divergence takes the form (A.1)where N_{i} are the coordinates of N in tubular coordinates { ∂/∂x_{i} } (see e.g. Frankel 2011). Comparing Eqs. (13) to (5), we see that, in our tubular coordinates,(A.2)For our tubular coordinate system, we found in Sect. 2.1.2 that g = ρR^{2} + ρ^{2}R^{3}(u_{1}sinθ − u_{2}cosθ). With N_{s} = 1 /λ, N_{θ} = ψ′/λ, and N_{ρ} = R′ρ/ (Rλ), we obtain Note that λ has an additional dependency on s, as the function θ is not just the coordinate but varies through the relation θ = θ_{0} + ψ(s). This is not the case for the derivatives of g which relate to the coordinate map T not the field line f.
All Figures
Fig. 1 Idealisation of the process by which an initially toroidal tube a) rotates about its apex to form a kinked structure b). This tube in b) appears as an S shape when projected down. The morphology of b) is characteristic of a sigmoidal geometry. 

Open with DEXTER  
In the text 
Fig. 2 Visualising the field creation process. a) Depicts a set of braided curves (in the form of the pigtail braid) on the left, and a set of twisted curves on the right. b) Depicts the pigtail braid embedded within a sigmoidal tube shape. c) Depicts the embedding of this tube in a Cartesian coordinate system, over the neutral line of a background quadrupolar flux distribution defined at the plane z = 0 (which plays the role of the photosphere in what follows). 

Open with DEXTER  
In the text 
Fig. 3 Tubular domains and curves. a) Depicts the tubular coordinate system T defined in Eq. (2). b) Depicts a tubular domain with u_{3}(s) = 0, ∀s. The curves drawn in this domain are curves of constant R and θ values. We see that the curves follow the shape of the tube. b) Depicts a tubular domain with the same axis curve as a), but with u_{3}(s) ≠ 0 so that the curves of constant (R,θ) value wind internally in the tube. 

Open with DEXTER  
In the text 
Fig. 4 Braiding geometry. Panel a) depicts a pigtail braid whose centrelines are curves r_{i} embedded in a tube T_{b}. Panel b) illustrates a stirring boundary motion which could be used to create the pigtail structure in a). 

Open with DEXTER  
In the text 
Fig. 5 Field lines of twisted a) and braided b) flux ropes (primary colours) and the background linear force free field (light colouring). Field lines can be seen to peel of the ropes as a result of the influence of the background field B_{b}. These are the initial conditions for the simulations in Sect. 4. 

Open with DEXTER  
In the text 
Fig. 6 Snapshots of the twisted core at t = 0a), t = 4.9b) and t = 9.8c). The field is seen to expand and flatten slightly. There also appears to be a counterclockwise (from above) rotation of field lines at the apex of the rope. 

Open with DEXTER  
In the text 
Fig. 7 Snapshots of the braided cores at t = 0a), t = 2.1b) and t = 4.2c). The curves are coloured based on their start points as belonging to one of the three initial braid footpoints. The interference of the background field is clear as an increasing number of field lines which would have originally belonged to the braid can be seen to have reconnected with the background field. 

Open with DEXTER  
In the text 
Fig. 8 Snapshots of contours of constant current density   j   = 0.5 for the twisted rope at t = 0a), t = 4.9b) and t = 9.8c). In a) we see the flux rope’s initial morphology, in b) the outer layer of current expanding to expose the rope’s remaining core, and in c) that the core is clear and forms a tubular shape similar in radius to the initial tube but with a lower apex. 

Open with DEXTER  
In the text 
Fig. 9 Snapshots of contours of constant current density   j   = 0.5 for the braided flux rope at t = 0a), t = 2.1b) and t = 4.2c). In a) the braid structure is clear. In b) the outer layer of current expands exposing a clear core seen in c); by comparison to the twisted case it is significantly fractured. 

Open with DEXTER  
In the text 
Fig. 10 Snapshots of contours of the Lorentz force magnitude for the twisted rope, at t = 0a), t = 4.9b) and t = 9.8c). In a) the background is forcefree whilst the tube has an associated Lorentz force. In b) we see the force beginning to dissipate.c) the force has all but disappeared except on the boundary. 

Open with DEXTER  
In the text 
Fig. 11 Snapshots of contours of the Lorentz force magnitude for the braided rope, at t = 0a), t = 0.49b) and t = 4.2c). In a) the background is forcefree whilst the braid structure has an associated Lorentz force. The force has reduced significantly in b), leaving a fracture structure. c) The force has all but disappeared except on the boundary. 

Open with DEXTER  
In the text 
Fig. 12 Plots of the magnetic, kinetic, internal and (normalised) total energy as a function of time for both the braided and twisted simulations. The, magnetic, kinetic and internal energy plots show the difference from their values at t = 0. In both cases the magnetic energy is converted into kinetic and internal energy. Two phases can be identified. The first is a rapid change as the magnetic is converted into kinetic energy and internal energy increase (t ∈ [ 0,0.1 ] for the braid and t ∈ [ 0,35 ] for the twisted field). In the second phase the kinetic energy levels off and then falls steadily, whilst the internal energy increases gradually. The loss of magnetic energy of the twisted field is larger than the braided field. 

Open with DEXTER  
In the text 
Fig. 13 Viscous (unbroken lines) and ohmic (dashed lines) heating during the simulations. In a) we see, similar to the energy plots, the braided quantities increase at a quicker rate than the twisted quantities. In both cases the viscous dissipation is larger than the ohmic heating by roughly a factor of 5. b) Shows the same quantities with the braided quantities scaled by the ratio of the initial twisted to braided energy (1.65), at the point where the simulations end the viscous dissipation is roughly the same. 

Open with DEXTER  
In the text 
Fig. 14 Snapshots of the vertically integrated current for the twisted field at t = 0a), t = 4.9b) and t = 9.8c). In a) we see the Zsigmoidal shape of the initial tube. b) Shows that the tube’s outer current layer has expanded. In c) the outer shell has further expanded leaving an inner rope structure which is far less kinked than in a). The rotation of the tube core from a) to c) is counterclockwise. 

Open with DEXTER  
In the text 
Fig. 15 Snapshots of the vertically integrated current at t = 0a), t = 2.1b) and t = 4.2c). In a) we see the braided structure projected. b) Shows that the tube’s outer current layer has expanded. In c) the outer shell has further expanded leaving an inner shell with a significantly fractured set of current peaks. 

Open with DEXTER  
In the text 
Fig. 16 Snapshots of the local twisting distribution L_{f}(x,y) on a subdomain of the boundary plane which contains the initial tube footpoint, at the beginning a), middle b) and end c) of the twisted field simulation. The sign is largely negative except at the tube edge. 

Open with DEXTER  
In the text 
Fig. 17 Snapshots of the local twisting distribution L_{f}(x,y) on a subdomain of the boundary plane which contains the initial tube footpoint, at the beginning a), middle b) and end c) of the braided field simulation. Both positive and negative values are present and the magnitude increases in b) and c). 

Open with DEXTER  
In the text 