Flux Rope Formation Due to Shearing and Zipper Reconnection

Zipper reconnection has been proposed as a mechanism for creating most of the twist in the flux tubes that are present prior to eruptive flares and coronal mass ejections. We have conducted a first numerical experiment on this new regime of reconnection, where two initially untwisted parallel flux tubes are sheared and reconnected to form a large flux rope. We describe the properties of this experiment, including the linkage of magnetic flux between concentrated flux sources at the base of the simulation, the twist of the newly formed flux rope, and the conversion of mutual magnetic helicity in the sheared pre-reconnection state into the self-helicity of the newly formed flux rope.


Introduction
Eruptive solar flares and coronal mass ejections (CMEs) contain large highly twisted magnetic flux ropes, and a key question is how the twist is created, since the twist in the preeruptive state is much smaller than what is observed later in the solar wind in an interplanetary coronal mass ejection (ICME) or magnetic cloud (MC) (Webb, 2000;Démoulin, 2008;Vourlidas, 2014).
Prior to the onset of an eruption, the magnetic structure around a prominence is highly sheared and slowly evolves through a series of equilibria. A magnetic flux rope is often present, or is formed during the earliest stages of the eruption. Many studies of the buildup to the onset of flares and CMEs have been undertaken, with a flux rope being formed in several ways, including flux emergence (Archontis and Hood, 2008), cancellation of flux (van Ballegooijen and Martens, 1989), and reconnection at a quasi-separator  or a separator (Longcope and Beveridge, 2007). Our investigation will study the flux rope formation process itself, rather than the eruption onset either due to a loss of equilibrium (Priest and Forbes, 1990;Lin and Forbes, 2000) or an equilibrium instability (e.g. Priest and Forbes, 2000;Priest, 2014) often attributed to either the kink instability (Hood and Priest, 1979;Fan and Gibson, 2003;Török, Kliem, and Titov, 2004) or the torus instability (Török and Kliem, 2005;Kliem and Török, 2006), which controls most of the dynamic and explosive behaviour witnessed during a flare or CME.
It is generally accepted that one effect of the eruption process is to drive reconnection in the region below the erupting flux rope. Such underlying reconnection may either create a new flux rope or increase the flux and twist of a pre-existing weakly twisted rope (Gibson et al., 2004(Gibson et al., , 2006. The erupting prominence is invariably seen to be more highly twisted than before it erupted (e.g. Mackay et al., 2010;Mackay and Yeates, 2012). Furthermore, three-dimensional (3D) reconnection naturally tends to create twist (Berger and Field, 1984;Hornig and Priest, 2003;Priest, Longcope, and Janvier, 2016).
The characteristics of flux ropes observed in interplanetary space have also been compared with the properties of their source regions at the Sun. In such cases, the poloidal flux has been found to be roughly equal to the total magnetic flux that is reconnected at the Sun (Qiu et al., 2007), suggesting that flux ropes in the interplanetary medium are mainly created at the Sun during the initial phase of a CME by magnetic reconnection. An interesting feature of solar flares and CMEs (described by e.g. Fletcher, Pollock, and Potts, 2004) is that they have two distinct phases in which the nature of reconnection changes (Yang et al., 2009;Qiu et al., 2010). During the rise phase, a flare brightens in Hα, at one point on each flare ribbon, and there is a single flare loop connecting the two points in extreme ultra-violet (EUV) and hard X-rays. These chromospheric brightenings spread rapidly along the polarity inversion line (PIL) to form the flare ribbons, and at the same time, a whole arcade of flare loops is formed whose footpoints lie in the ribbons. The ribbons and arcade are created by what Priest and Longcope (2017) call "zipper reconnection". Then the main phase of the flare is created by "main-phase reconnection" and is characterised by the outwards motion of the ribbons in a direction perpendicular to the PIL and the rise of the flare arcade. Zipper reconnection explains the initial motion of flare brightenings parallel to the PIL. It builds up magnetic twist and flux in the core of the erupting flux rope during the rise phase. Then, during the main phase, reconnection adds more twist and flux to the central core of the flux rope. Concepts of magnetic helicity conservation and the initial configuration itself are used to quantify the build-up of both twist and magnetic flux in the erupting flux rope (Priest and Longcope, 2017).
Zipper reconnection has features in common with "flux cancellation" (Martin, Livi, and Wang, 1985;van Ballegooijen and Martens, 1989) and "tether cutting" (Moore et al., 2001;Amari et al., 2003a,b;Aulanier et al., 2010) eruption models, but two new factors in our analysis are the role of magnetic helicity conservation in calculating the resulting flux rope twist, and the progression of the reconnection along the PIL. Flux cancellation refers to the cancellation of magnetic flux at the photosphere, usually driven by photospheric motions towards the PIL and leading to the formation of a flux rope (and subsequently a prominence), whereas zipper reconnection can also be initiated after motions parallel to the PIL, can occur above the photosphere, and is involved in the rise phase of a two-ribbon flare. Tether cutting reconnection was proposed as a means of initiating a flare by reconnecting two sheared flux ropes under a coronal arcade to form a sigmoid, whereas zipper reconnection refers to the reconnection of a series of loops placed along and inclined to the polarity inversion line and is a response to the eruption.
In light of such models, several recent papers are particularly relevant to our work. Wang et al. (2017) combine solar and interplanetary observational data to give a detailed overview of the dynamic formation of a single magnetic flux rope in the solar corona during a tworibbon flare. Crucially, these authors derive estimates for how the poloidal and toroidal flux and the field line twist evolve over time. From a modelling perspective, Inoue et al. (2018) use cutting-edge techniques, combining non-linear force-free field (NLFFF) extrapolations (based upon magnetogram images of an active region that produced an M6.6-class flare) and magnetohydrodynamic (MHD) simulations of the formation and eruption of a magnetic flux rope due to tether-cutting reconnection. While much focus is placed on the eruptive phase of the flux tube development, these authors also evaluate the build-up of magnetic twist both before and after the tether-cutting reconnection forms the flux rope.
In this paper, we set up a numerical experiment to study zipper reconnection and determine whether it can indeed form a large flux rope. In Section 2 we describe our simple numerical model, where a pair of magnetic bipoles (initially linked by a pair of untwisted magnetic flux tubes) are sheared and undergo magnetic reconnection to form a twisted overlying magnetic flux rope. We study several properties of this resulting configuration and its evolution over time in Section 3; these properties include the magnetic flux and connectivity at the base of the model (in Section 3.1), the formation and evolution of the flux rope (in Section 3.2), and an analysis of the flux rope twist and self-helicity (in Section 3.3). We discuss our findings and outline our conclusions and future work in Section 4.

Computational Setup
We model a single "zippette" flux rope formation by starting with two untwisted flux tubes beside each other and allow them to relax to reach a near-potential state. Once the overall field has relaxed, one tube is then moved with respect to the other so that the whole configuration is skewed, as shown in Figure 1. Unlike other similar previous experiments, the shear phase alone is responsible for the configuration that evolves once reconnection begins. Previous experiments combine the shear phase with inflow towards the PIL in order to initiate reconnection. Such an inflow is not necessary in our experiment.
The evolution of the magnetic field is followed by solving the MHD equations where ρ is the plasma density, p the plasma pressure, v the velocity, B the magnetic field, j = ∇ × B/μ the current density, η the magnetic diffusivity, γ = 5/3 the ratio of specific heats, and μ = 4π × 10 −7 is the magnetic permeability. In addition, ∇ · B = 0. These are solved using the Lagrangian remap code, (Lare3D) described in Arber et al. (2001). Lare3D J. Threlfall et al. uses dimensionless variables, based on a magnetic field strength, B 0 = 10 G, a length, L = 10 7 m, and a mass density, ρ 0 = 1.67 × 10 12 kg m −3 . Thus, the typical Alfvén speed is V A = 690 km s −1 and the typical time is t = 14.5 s. The reference current density is j 0 = B 0 /(μL) = 8×10 −5 A m −2 and the reference magnetic diffusivity is η 0 = 6.9×10 12 m 2 s −1 .
The magnetic diffusivity consists of two parts, namely a background value, η b , and an anomalous value, η a , that is only non-zero when the current rises above a critical value. This allows the formation of strong currents and ensures that the reconnection remains spatially local. In dimensionless variables, and η b = 10 −4 and η a = 10 −3 . In the subsequent sections, time units are quoted in Alfvén times, τ A = L/V A , and lengths in terms of L. The computational domain considered is −8 ≤x ≤ 8, −8 ≤ȳ ≤ 8 and 0 ≤ z ≤ 20 (where barred quantities represent dimensionless variables in the numerical domain).

Initial Configuration
The initial (pre-shear) configuration can be seen in Figure 2. The vertical magnetic field on the photospheric boundary consists of four sources, two positive and two negative, and each has the same form, ±B 0 exp (−r 2 /r 2 0 ), and flux. r is a local radial coordinate centred on x = ±1 and y = ±2.5 and r 0 is a measure of the width of each source. The contours of vertical magnetic field strength in Figure 2a illustrate the magnetic configuration imposed on the base of the simulation, while Figure 2b shows the corresponding 3D magnetic field, illustrating that the configuration initially consists of two untwisted flux tubes. The field is left to relax until t = 1200, by which time it is very close to potential. To simplify the discussion, we label each of the strong polarity regions in Figure   regions of positive or negative polarity in the positive x domain (with R + initially centred on [1, 2.5] and R − initially centred on [1, −2.5]). Analysis of the initial connectivity linking these regions (seen in Figure 2b) shows that field lines that originate in L + (R + ) all terminate in L − (R − ); there is no field linking the left and right halves of the domain in the initial configuration.

Shearing Motions
From t = 1000 to t = 1600, we impose a slow, sub-Alfvénic shearing velocity on the photospheric base that switches on at t = 1200 and off at t = 1400. Thus, on z = 0 and for Both v x and v z remain zero on this boundary and all other variables have zero normal derivative. This shearing motion moves both sets of sources until L + is brought approximately opposite R − . The source motion is not completely ideal; all four sources are slightly spread out during their movement. As L + and R − move closer together, the current density above the photosphere increases. From t = 1400 until t = 1600, the driving has stopped, but the current above the photosphere now exceeds j crit and the loops begin to reconnect. Our experiment is primarily designed to study the reconnection of these super-critical "coronal" currents. Our inclusion of a smaller background resistivity means that in addition, we also anticipate a small amount of source diffusion and possibly flux cancellation when opposite sources are close together. While non-negligible, these effects will play a secondary role to the super-critical currents and the value of η a that drive zipper reconnection, and that dominate the behaviour revealed in this experiment.

Results
We now examine various aspects of the results of this experiment.

Evolution of Base Connectivity and B z Flux
We first assess how the connectivity of the sources evolves over time, both during and after the shearing phase of the experiment. In order to do so, we must first redefine the regions labelled L + , L − , R + and R − . The initial positive and negative sources (shown in Figure 2b) overlap. Distinct contours for sources in either half-plane can only be formed for values of |B z | > 0.37 (as shown in the original image); smaller values of B z are unable to distinguish between right and left sources of the same sign. For the time being, the positive flux source core regions are defined by the critical B z = 0.37 contours, while the negative source core regions are defined by B z = −0.37. It should be noted that such definitions mean that our "sources" already omit some of the total flux imposed at the base.
To estimate how the connectivity of our sources evolves, we calculate a Cartesian grid of initial positions within each positive or negative sources, before tracing a magnetic field line through the domain from each point. If the final position of the field line lies in the left (x < 0) half-plane, the initial position from which the field line is traced is coloured blue, while initial positions whose field lines terminate in the right (x > 0) half-plane are coloured red. For brevity, we only present the connectivity of field lines traced from the positive B z sources (L + and R + ) in Figure 3 at four different stages during and after the shearing phase of the experiment.
In Figure 3, we build up a picture of how the connectivity of each source changes with time. We can see that no connectivity changes are apparent during the shearing phase (t = 1300τ A , in Figure 3a) or at the end of the shearing phase (t = 1400τ A , in Figure 3b). However, by t = 1500τ A , Figure 3c shows the existence of a large portion of the R + region coloured blue, meaning that links have formed with the left half-plane. A very small red inclusion is also notable in L + ; R + has much greater connectivity changes than L + . A short time later, (t = 1600τ A , in Figure 3d), the blue inclusion in R + and the red inclusion in L + are much more similar in size. It is also apparent from Figure 3 that the area surrounded by contours of |B z | = 0.37 is decreasing with time. Indeed, by t = 2000τ A , no flux is visible on the base that satisfies |B z | > 0.37; the polarity regions broaden over time, reducing the peak value of B z as they broaden.
To account for this effect, and to continue to monitor the connectivity of each of the source regions, we reduce the critical value of B z used to define the core of each source region. After t = 1600τ A , each strong polarity region is instead defined using contours of |B z | = 0.2. At this value, at t = 1600τ A , our sources once again remain distinct (weaker B z values would not reveal two positive and two negative sources) and broad (stronger B z values reduce the area of each source).
The connectivity maps in Figure 4 show that the sources at the base continue to diffuse at and after t = 1600τ A , and that the critical contours of |B z | = 0.2 omit large regions of weak positive and negative flux. Regarding the connectivity, we once again see that the blue incursion into R + is larger than the red incursion in L + at t = 1600τ A and remains so throughout the remainder of the experiment. The connectivity map in R + also evolves somewhat differently to L + : as the blue inclusion increases in area, a red feature appears to "break off" from the primary red region from t = 2000τ A . By the end of the experiment, most of the flux in the source that remains connected to the right half-plane is only found on the extreme right of R + , while a separate "island" of red field lines has now formed near the centre of R + in Figure 4d.
Another quantity of interest is the amount of magnetic flux through each source, and how much flux links which sources. In order to estimate this, we evaluate the vertical magnetic field B z at each location on the grid (previously used to perform field line tracing). To estimate the total flux through each source, we sum B z in x and y directions. We also estimate the total flux that ultimately links to either right or left half-planes, only summing values in areas with field lines that link to x > 0 or to x < 0. We display the total flux (and its components, linked either to the right or left half-planes) in Figure 5, comparing the flux through |B z | = 0.37 contours in Figure 5a   There are several noteworthy aspects of Figure 5. As expected, the total flux through both R + and L + decreases with time. The total flux through R + is greater than the total flux through L + during the driving phase. Also during this phase (and also illustrated by Figures 3a-3b), the total flux in L + is entirely linked to the left half-plane, while the total flux in R + is entirely linked to the right half-plane. This changes almost immediately after the driving phase ends: flux from R + terminating at x < 0 rises from zero at approximately 1410τ A , matched by a corresponding drop in flux linking R + with x > 0. Evidence of changing connectivity in L + occurs much later, with Figure 5a demonstrating an increase from zero of L + flux linked to x > 0 at approximately 1510τ A , with a corresponding decrease in L + flux linkage with x < 0.
Turning to the evolution shown in Figure 5b, decreasing the critical value of |B z | used to define each source region widens the area covered by each source and therefore increases the flux values recovered (relative to Figure 5a, which uses |B z | = 0.37 to define each source). It can be seen that the dominant connectivity of the R + region changes at later times: prior to t ≈ 2100τ A , most of the total flux through R + terminates in the right half-plane. From  t ≈ 2100τ A until the end of the experiment, the majority of flux through R + now links to the left half-plane. Such an effect also occurs in the flux through L + , whose connectivity with the left half-plane always remains greater than the connectivity with the right half-plane. It is also worth noting that the connectivity of both sources appears to be tending towards equipartition, approaching equal levels of flux through each source linking to the left and right half-planes.
Finally, we present a comparison of the angles between each of the specific source regions, as defined in Figure 1, in order to estimate the magnetic helicity. Each angle is defined between the locations of peak positive and negative B z flux on the base of the simulation domain. These locations are not affected by any choice of critical B z threshold used to outline the sources. The evolution of the angles between these locations is seen in Figure 6.
The figure shows that θ 1 and θ 2 are identical until the velocity driver is switched on, while θ 3 and θ 4 remain identical for all time. When the velocity driver is switched on, θ 1 increases rapidly, rising above π/2 before reducing to a narrower angle than seen before the driving phase. The rise and fall can be attributed to the L + source region moving relative to R + and R − . The angle θ 1 is widest when L + lies at a y-value halfway between R + and R − (as approximately shown in Figure 3a). L + continues to move down below y = 0, while R + and R − move upwards in y, causing θ 1 to become much narrower by the end of the driving phase (shown e.g. in Figure 3d). Meanwhile, according to Figure 6, θ 2 (and indeed θ 3 and θ 4 ) reduces over the course of the driving phase, as the angle subtended by R + , L − , and R − closes and the sources are progressively sheared.

Overlying Flux Rope: Formation
One of the goals of this experiment was to establish if this process could lead to the formation of new flux ropes in the manner described in the zipper model (Priest and Longcope, 2017). By tracing selected magnetic field lines in the domain, it is easy to see that a new twisted flux rope does indeed form between sources R + and L − after the shearing phase of the experiment ends. As shown in Figure 7, this flux rope expands to fill the simulation domain over time, and overlies an arcade of magnetic field lines that link L + and R − .

Overlying Flux Rope: Twist and Helicity
The field lines that form the overlying flux rope in Figure 7 are twisted, forming a helical path around a central axis. A key question is what the twist of these field lines is (relative to the flux rope axis) and how the twist is distributed throughout the tube. To calculate the twist at any location in the tube, we again create a grid of initial positions, this time distributed in the xz-plane at y = 0. The grid is restricted to lie within a chosen contour of the perpendicular magnetic field strength (i.e. B 2 x + B 2 z ), which we use to define our flux rope cross-section. The grid is used to trace field lines forwards (to one negative flux source at the base) and backwards (to the positive source), with the recovered locations at the base then used to calculate an angle of rotation relative to the flux rope axis (defined where B 2 x + B 2 z = 0 in the [x, y = 0, z] plane). In Figure 8, we have colour-coded each position on each grid to represent the angle through which that field line rotates around the flux rope axis in three different snapshots as the flux rope rises and expands (and have included the perpendicular field strength in the background for context).
From Figure 8, we see that the majority of field lines within the flux rope are twisted by between 1.9 -2.5π . The highest values of twist form in two thin bands, which broaden over time. By the end of the experiment, the peak values of twist are focussed close to the axis of the flux rope, but small regions of weaker twist (< 1.7π ) have become apparent at the edges of the flux tube.
The twist of field lines that form the flux rope (about its axis) is a crucial component of magnetic helicity. The conservation of total magnetic helicity is a natural consequence of magnetic flux tube formation and eruptions, since it is conserved during 3D reconnection (see e.g. Priest, Longcope, and Janvier, 2016). The total magnetic helicity (H T ) can be expressed as the sum of the self-helicity (H s ) resulting from the twist of each flux tube in the system and the mutual helicity (H M ) of each tube relative to each other; for a system containing N tubes, this amounts to Hence, in a system containing two tubes labelled A and B, this becomes The self-helicity of any flux rope can be calculated as a product of the twist (φ) and the magnetic flux F passing through the tube (with a perpendicular cross-section S). For simplicity, we use the definition outlined in Priest, Longcope, and Janvier (2016), using the mean twist of the flux tubeφ, namely, In Figure 9 we compare the average twist and self-helicity of the flux rope over time, and include the kinetic and magnetic energies of the system to provide context. The average angle of footpoint rotation shown in Figure 9a is calculated using the method described earlier; creating a grid of initial positions throughout the midplane of the flux rope in each simulation snapshot, tracing field lines back to their photospheric sources, and evaluating the rotation of each field line at one source to its position at the other source relative to the flux rope axis. In order that we may better compare with other works that measure twist, we have also used the method of Berger and Prior (2006), in which a twist parameter (t w ) evaluates the integrated parallel current along a given field line, i.e.
where dl is an arc length along the field. From Figure 9a, both measures differ slightly in values recovered, but the overall trends agree. A direct comparison of the two measures and how (in certain cases) these values are related is given in Appendix A. The average field line twist of the flux rope has a maximum shortly after the shearing phase ends, after which it falls, tending to one turn (or a footpoint rotation close to 2π ) by the end of the experiment. This result is, in some sense, to be expected. The flux rope forms almost immediately after the shearing phase ends, before expanding to fill the simulation volume. The formation itself occurs through reconnection at a large current sheet above the photosphere created by the shearing motion. This reconnection imparts a large amount of twist (up to 3π along some field lines) into the newly formed flux rope. Following the flux rope formation, additional reconnection (of smaller, fragmented sub-critical currents below the flux rope) allows more field lines to be brought into the flux rope itself; this (and the lack of overlying field) causes the flux rope to expand, while much of the newly added field is twisted at or slightly below 2π , reducing the average twist over time.
The energies shown in Figure 9b also follow this pattern. As expected, our configuration contains much more magnetic energy than kinetic throughout the experiment. The background resistivity causes a slow decrease in magnetic energy over time, but the fastest rate of magnetic energy release (other than during the very early relaxation phase) can be found immediately and briefly after the shearing phase ends. At this point, anomalous resistivity has begun to act on the current above the critical value. The shear phase also briefly injects a small amount of kinetic energy into the system, but in general, the system continues to lose kinetic energy over time. Using the footpoint field line rotation as the measure of twist, Figure 9b suggests that the self-helicity of the flux rope steadily increases from zero over time. The increase in self-helicity shown in Figure 9b results from increasing amounts of magnetic flux in the tube. As the magnetic flux rope expands, magnetic reconnection below the tube continues to draw in new magnetic flux, contributing to the flux rope expansion and increasing levels of flux, while slowly reducing the average field line twist.
The mutual helicity can be found by using the angles θ 1 and θ 2 (for flux tubes A and B lying parallel to one another) or θ 3 and θ 4 (in the case of a flux tube C that overlies a second flux tube D), identified in Figure 1, namely that However, the total magnetic flux through a specific contour of B z at the base is not constant due to the broadening of magnetic flux sources. Using conservation of helicity, Priest and Longcope (2017) estimated that both the overlying and underlying flux tubes would have a twist of π , provided that the self-helicity is split equally between the two tubes. In our experiment, however, diffusion at the computational base implies that only an overlying tube forms. Conservation of total magnetic helicity implies that its average twist will therefore be 2π , as observed, since all the self-helicity goes into the overlying tube.
Finally, we broadly place our findings in the context of other recent investigations. In the event studied in Wang et al. (2017), two reconnection phases are recovered; however, the first phase (which, in the conceptual model of Priest and Longcope, 2017, would correspond to the zipper reconnection phase) displays characteristics of both zipper reconnection and main-phase reconnection. This rules out a direct comparison between the temporal evolution of magnetic twist in both cases. Our recovery of a magnetic flux rope with an average twist of 2π might initially appear to be much less than the typical values of twist found by Wang et al. (2017). However, we note that our model represents only one single component (a socalled "zippette") of the full zipper reconnection picture (which consists of many zippettes occurring one after the other). Now that we have shown that it is possible to form a magnetic flux rope through this mechanism, the next step would be to investigate how several zippettes might combine to form a full zipper reconnection event, and how each zippette event contributes to the evolution of magnetic twist in such a flux rope. Such a model would provide a better comparison for the observations of Wang et al. (2017). Our twist values are much more closely aligned to those found in Inoue et al. (2018). Indeed, in Figure 9a we evaluate the field line twist using the same method, and recover values that agree well with those found in their model (noting that we only evaluate the twist parameter of field lines that lie within the flux rope, while Inoue et al., 2018, calculate this for many field lines, and use this to infer that a twisted flux rope has indeed formed on some of those field lines, when t w increases above unity). However, we would also emphasise that a direct comparison between values of t w and the twist given by the angle of footpoint rotation relative to axis footpoints, even in relatively simple cases, can be complex and potentially misleading (as shown in Appendix A, and discussed further in e.g. Liu et al. (2016)).

Conclusions and Future Work
In conclusion, we have demonstrated that a helical flux rope can be formed directly by shearing and reconnection of two untwisted flux tubes (without the need for an inflow velocity). We have examined the properties of this system in detail, observing that the outermost flux sources in the system preferentially undergo reconnection. Although the sources overlap and diffuse slightly at the base, we find that the average level of twist in the flux rope tends to 2π by the end of the experiment. This twist can be explained simply by evoking total magnetic helicity conservation and so equating the final self-helicity of the flux rope to the initial mutual helicity of the system. This is only the first step in an investigation into the mechanism of zipper or zippette reconnection. Whether and how a sequence of such events can form a single erupting flux rope overlying a magnetic arcade of post-flare-loop-like structures remains to be seen. The magnetic topology present in this configuration remains relatively simple with the initial configuration consisting of two untwisted flux tubes side by side. In future, we plan to conduct experiments with much more realistic magnetic structures and to decrease the diffusion near the base, so as attempt to show how a single erupting magnetic flux rope overlying a magnetic arcade of flare loops may be produced.
Key questions in future include: Is the same reconnection achieved if the positive and negative sources do not overlap? What is the effect of modifying the form of the velocity driver applied at the base to move the flux sources as solid bodies? What is the effect of reducing the diffusion near the base so as to prevent the sources spreading out?
Acknowledgements The authors gratefully acknowledge the financial support of STFC through the Consolidated grant, ST/N000609/1, to the University of St Andrews. This work used the DIRAC 1, UKMHD Consortium machine at the University of St Andrews, the DiRAC Data Centric system at Durham University, operated by the Institute for Computational Cosmology, and the DiRAC Data Analytic system at the University of Cambridge, operated by the University of Cambridge High Performance Computing Service. These systems are operated on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment was funded by BIS National E-infrastructure capital grants (ST/K00042X/1 and ST/K001590/1), STFC capital grants (ST/K00087X/1, ST/H008861/1 and ST/H00887X/1) and DiRAC Operations grants (ST/K003267/1 and ST/K00333X/1). DiRAC is part of the National E-Infrastructure.

Disclosure of Potential Conflicts of Interest
The authors declare that they have no conflicts of interest.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. 1 + r 2 2 0 /L 2 .
(16) Figure 10 shows the variation of and t w with r for the twist profile defined above. It is clear that 2 0 and t w are almost identical for the constant twist field. However, there is a significant difference between the two measures when the non-uniform twist profile is used. In particular, while remains positive throughout the loop cross section, t w changes sign. In addition, if we calculate the average value of t w across the circular cross section of the loop for which r L, then the value is zero. This is not true for the twist. While it is computationally easier to calculate t w in a more general, non-cylindrical loop, we believe that (where possible) it is better and clearer to define the twist as the angle a field line rotates about the loop axis. Figure 10 Radial variation of footpoint twist ( ) and integrated parallel current (t w ), using the example twist profile defined by Hood et al. (2016), where L = 10, a = 1 and 0 = 1. Uniform twist profiles are shown as dashed lines, while radially varying profiles are thick solid curves (for key, see legend).