## Abstract

Predicting the final state of turbulent plasma relaxation is an important challenge, both in astro-physical plasmas such as the Sun's corona and in controlled thermonuclear fusion. Recent numerical simulations of plasma relaxation with braided magnetic fields identified the possibility of a novel constraint, arising from the topological degree of the magnetic field-line mapping. This constraint implies that the final relaxed state is drastically different for an initial configuration with topological degree 1 (which allows a Taylor relaxation) and one with degree 2 (which does not reach a Taylor state). Here, we test this transition in numerical resistive-magnetohydrodynamic simulations, by embedding a braided magnetic field in a linear force-free background. Varying the background force-free field parameter generates a sequence of initial conditions with a transition between topological degree 1 and 2. For degree 1, the relaxation produces a single twisted flux tube, whereas for degree 2 we obtain two flux tubes. For predicting the exact point of transition, it is not the topological degree of the whole domain that is relevant, but only that of the turbulent region.

## 1. Introduction

Self-organization of turbulently relaxing plasma to a predictable minimum-energy state has been observed in laboratory confinement devices including the reversed-field pinch and the spheromak [1–4]. The so-called Taylor relaxation hypothesis assumes that the only relevant constraints on the dissipation of magnetic energy are the total magnetic flux and the total magnetic helicity. The latter is not an exact invariant in the presence of resistivity, but is known to be well preserved on typical timescales of relaxation processes. In order that all other ideal invariants are destroyed (such as helicity in subregions of the plasma [2], or other helicity moments [5]), the evolution must be sufficiently turbulent that magnetic reconnection is able to occur throughout the volume.

It has also been proposed that this Taylor relaxation theory might be applied to predict the energy released by rapid heating events in the solar corona [6], where magnetic energy is believed to be released through relaxation to a lower energy equilibrium. In this context, numerical magnetohydrodynamic (MHD) simulations have modelled the dynamic relaxation of various initially unstable equilibria, such as kink-unstable twisted magnetic flux ropes [7–11], or a magnetic field with a braided structure [12,13]. Our work has been motivated by the latter simulations, which showed that certain initial configurations self-organized into final equilibria whose magnetic topology was more complicated than predicted by Taylor theory, despite the occurrence of efficient reconnection. We identified the presence of an additional constraint beyond the total magnetic flux and helicity: the topological degree of the field line mapping [14,15].

The topological degree (defined in §2) is conserved provided that the degree of the boundary does not change. The latter can be ensured by having turbulent dynamics that are localized in the interior of the domain and do not affect the boundary. It is our goal in this paper to show, for a sequence of initial conditions of degree 1 which approach degree 2, how the final state suddenly switches from a single flux tube to a pair of flux tubes.

The assumption of localization is an important one for relaxation events in the solar corona. Unlike the reversed-field pinch, there are no conducting walls to define the relaxation volume [16]. Typically, coronal energy releases–for example, in solar flares–are highly localized in space. The extent of the relaxation region is determined by the connectivity of the magnetic field configuration, requiring either unstable configurations or very small-scale gradients to initiate the energy release. Dixon *et al.* [17] showed that Taylor theory may be applied to regions with a free boundary, although they did not specify where the boundary should be placed in any particular magnetic field. More recently, Bareford *et al.* [11] have shown that Taylor theory can give reasonable predictions of relaxed states in numerical solutions of kink-unstable magnetic flux tubes, provided that it is applied within the appropriate subregion.

Localized Taylor relaxation has also been applied in the context of tokamaks. In these devices, global Taylor relaxation does not describe the magnetic configurations that are observed. However, Hudson *et al.* [18] have developed a partial relaxation model where Taylor relaxation occurs in sub-volumes. These sub-volumes are separated by a discrete set of irrational flux surfaces that survive even in the presence of the chaotic field lines typical of non-axisymmetric magnetic fields. In another application, Gimblett *et al.* [19] have developed a model for edge-localized modes based on localized Taylor relaxation within only the outer region of the plasma.

In this paper, we consider a one-parameter family of initial magnetic configurations in a periodic (topologically toroidal) domain. These configurations, described in §3a, are chosen to have a ‘background field’ of gradually varying structure. This complements the particular configurations where this constraint was demonstrated previously [14], which had vanishing magnetic helicity.

## 2. Topological degree constraint

To define the topological degree of a particular configuration, let *f*:*D*_{0}→*D*_{1}, where *f*=( *f*_{x},*f*_{y}), be the field line mapping from the lower boundary *D*_{0} to the upper boundary *D*_{1}. In other words, *f*(**x**_{0})∈*D*_{1} is the endpoint of the magnetic field line starting at **x**_{0}∈*D*_{0}. We assume that there is a strong enough guide field that all field lines pass from *D*_{0} to *D*_{1} without changing direction. We shall assume for simplicity that *D*_{1}=*D*_{0}, as in the periodic simulations presented in this paper. Field lines that satisfy *f*(**x**_{0})=**x**_{0} are known as fixed points of *f* (or periodic orbits in the case of periodic boundaries). The index of a fixed point describes the local structure of *f* around the fixed point and is defined as the local Brouwer degree of *f* (for more details, see [15]). Now let *D*⊂*D*_{0} be a subregion of *D*_{0}. The topological degree of *f* on *D*, denoted *T*(*D*), is defined to be the total (net) fixed point index, obtained by summing the indices of all isolated fixed points of *f* in *D*. One may express *T*(*D*) as the Kronecker integral
*D* [20]. As *T*(*D*) is an integer, the only way it can change under a continuous time-evolution of *f* is if one or more fixed points cross into or out of the boundary of *D*. So if *f* is fixed on the boundary of our turbulent region *D*, then *T*(*D*) must be preserved in time. In particular, this means that the relaxed state may be forced to contain more than one fixed point, implying certain magnetic substructure.

We use the convenient colour map technique introduced by Polymilis *et al.* [20] for visualizing fixed points of *f*, their indices and *T*(*D*). This is illustrated in figure 1 with the magnetic field
*x*,*y*) in *D*_{0}, according to the relative signs of *f*_{x}−*x* and *f*_{y}−*y*. Fixed points are readily identified as places where all four colours intersect. Furthermore, the topological degree *T*(*D*) of a region *D*⊂*D*_{0} may be identified by noting the anticlockwise sequence of colours around the boundary of *D*. In particular, the number of times that the full sequence of four colours (in the correct order) is repeated. For example, the degree of the full region shown in figure 1 is +1. Correspondingly, there is a net anticlockwise rotation of field lines around the boundary. Inside *D*, there are three fixed points: two ‘elliptic’ points with degree +1 and one ‘hyperbolic’ point in the centre with degree −1.

The topological degree relates the complexity of the field on the boundary of the domain to that of the interior field. This is similar to how Gauss' theorem relates the integrated electric field over a closed surface to the electric charge inside the surface. As for the topological degree, the surface integral over the electric field does not distinguish how many positive or negative electric charges are inside the domain: it only gives a net charge. For the topological degree, the analogue of the net charge is the sum of hyperbolic (degree −1) and elliptic (degree +1) periodic orbits. The simplest state (the smallest number of charges which give the correct net charge) is typically also the one with lowest energy. Thus, an efficient turbulent relaxation within an otherwise ideal plasma is expected to lead to the simplest force-free field consistent with the topological degree of the turbulent region.

## 3. Numerical set-up

### (a) Starting configurations

In this paper, we present resistive-MHD simulations for a family of initial magnetic configurations. Each is a superposition of two components **B**=**B**_{α}+**B**_{braid}, where **B**_{α} is a linear force-free field with constant *α*, and **B**_{braid} is a braiding magnetic field pattern consisting of six toroidal rings of magnetic flux. The field **B**_{braid} is orthogonal to **e**_{z} and vanishes on the boundaries of our domain. By contrast, the background field **B**_{α} is non-zero on all six boundaries of our domain. (For numerical convenience, we use a Cartesian domain.) By varying *α* and keeping **B**_{braid} fixed, we obtain a one-parameter family of initial configurations. For *α*=0, the configuration has degree 2, while for all positive values of *α*, it has degree 1. By decreasing the value of *α* towards 0, we can test when and how the transition affects the relaxed state. To initialize the other variables in our resistive-MHD simulations, we simply take zero initial velocity, constant density and constant pressure.

Note that the combined field **B** is not in equilibrium and leads to a dynamical evolution in the resistive-MHD equations. Previous simulations (in the *α*=0 case) have found consistent final-state topology whether or not the field is first subjected to an ideal relaxation before initiating the resistive-MHD evolution [14].

For **B**_{α}, we take the well-known axisymmetric constant-*α* magnetic field of Lundquist [21]. In cylindrical coordinates (*r*,*ϕ*,*z*), this takes the form
*J*_{0} and *J*_{1} are Bessel functions of the first kind. This field is readily shown to satisfy ∇×**B**_{α}=*α***B**_{α} for constant *α*. In the limit *α*→0, it reduces to a vertical, current-free magnetic field **B**_{0}=*B*_{0}**e**_{z}. In this paper, we fix *B*_{0}=1. The condition that *B*_{z}>0 everywhere in our domain puts an upper limit on the acceptable values of *α*. In particular, we require *α*<*α*_{r}, where *α*_{r}≈0.21 is the smallest root of **B**_{α} leads to a net electric current in the *z*-direction.

The braiding field **B**_{braid} was introduced by Wilmot-Smith *et al*. [22]; its construction is based on the pigtail braid, with six toroidal rings of flux,
*x*_{i}=*k*_{i}=(1,−1,1,−1,1,−1), *z*_{i}=(−20,−12,−4,4,12,20). This pattern of flux is efficient at ‘mixing’ the field lines while having zero net helicity and leads to a demonstrably chaotic field line mapping in our periodic domain [13]. It is effectively this region of efficient mixing that generates small magnetic scales enabling current sheets to form, leading to magnetic reconnection. The extent of this region determines the region of turbulent relaxation in which the field is able to relax efficiently.

Figure 2 shows illustrative magnetic field lines for the combined states with *α*=0.001, 0.01, 0.05 and 0.1. Although the field line connectivity is significantly altered, **B**_{braid} is, energetically, a relatively small perturbation to the background field **B**_{α}. Denoting the magnetic energy by *E*_{mag}=〈*B*^{2}〉/(2*μ*_{0}), one has that *E*_{mag}(**B**_{0}+**B**_{braid})≈1.008*E*_{mag}(**B**_{0}), while *E*_{mag}(**B**_{0.1}+**B**_{braid})≈1.009*E*_{mag}(**B**_{0.1}). It should be noted that **B**_{α} is *not* the minimum-energy (Taylor) state for our configuration, except when *α*=0. This is because the magnetic helicity of the combined field **B** differs from that of **B**_{α}.

### (b) Numerical simulations

The Lare3D Lagrangian-remap code^{1} is used to solve the resistive-MHD equations in a Cartesian box {−8≤*x*≤8,−8≤*y*≤8,−24≤*z*≤24}, at resolution 320×320×240. We apply periodic boundary conditions in *z* and line-tied boundary conditions in *x* and *y*. The code solves the non-dimensionalized equations
*ρ* 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 density,

*η*the resistivity,

**the strain tensor and**

*ε**et al.*[23]. The viscous term ∇

**in (3.4) includes no background viscosity, but only a shock viscosity to prevent unphysical oscillations and approximate the jump in entropy across shocks. The shock viscosity takes the tensor form given in [11], and we use the same parameter values**

*σ**ν*

_{1}=0.1 and

*ν*

_{2}=0.5. There is a corresponding heating term

*ε***in (3.6). We initially set**

*σ**ρ*=1 and

*ϵ*=0.01 in non-dimensional 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 our box. The simulations presented here use a uniform resistivity of

*η*=5×10

^{−4}. Previous simulations of the

*α*=0 case found that the topology of the relaxed state is not sensitive to the choice of

*η*, although the details of the turbulent relaxation do change [12].

## 4. Results

For all values of *α*, there is an initial phase of turbulent relaxation until approximately *t*=100, followed by a more gradual resistive dissipation. This pattern is the same as the earlier simulations with *α*=0 [12] and was also seen for the relaxation of a kink-unstable loop [11]. Huang *et al.* [24] find a similar distinction between quasi-static resistive evolution and the onset of a dynamical phase, in resistive reduced-MHD simulations of a randomly structured field.

In the turbulent phase of our simulations, the dynamics consists of a cascade from initially large to smaller current sheets, which interact with one another to dissipate magnetic energy during the relaxation. Figure 3 shows the appearance of these current sheets at *t*=50 during the turbulent phase, in a cross section at the mid-plane *z*=0. In each case, there is a distinguished turbulent region outside which there are no significant currents or dynamics. The shape of the turbulent region is more circular for the run with *α*=0.1, owing to the influence of the different background field. Figure 4*a* shows the maximum current throughout the domain as a function of time. All four runs follow a bursty, intermittent pattern of maximum current in the turbulent phase, followed by a smooth evolution with lower maximum current during the gradual, resistive phase. The run with *α*=0.05 maintains a higher maximum current for longer than the others: this is due to the interaction of one of the resulting flux tubes with the background field, as will be discussed below.

The turbulent phase is also evident in the total energies shown in figures 4 and 5. For example, the total kinetic energy *E*_{kin}=〈*ρv*^{2}〉/2 is significant mainly during the turbulent phase and follows a quite similar pattern in all runs. The oscillations seen in *E*_{kin} and also in the magnetic energy *E*_{mag} have a period consistent with torsional Alfvén waves, launched from the initial flux ring locations and counter-propagating in *z*. Although these waves dominate the frequency spectrum, the dynamics are nonetheless turbulent in the sense that the chaotic field line mapping produces a cascade to smaller spatial scales throughout the braided region. This cascade and the resulting unpredictable bursts of reconnection are important for removing energy and restructuring the magnetic field. Previous simulations of the *α*=0 case have shown that consistent relaxed states are obtained whether the initial state contains discrete flux rings (as here), or is first subjected to an ideal relaxation [12], in which case there is a broader frequency spectrum.

In our resistive simulations, the dissipation of magnetic energy must be compared to that of the background **B**_{α} field under resistive diffusion alone. Figure 5 shows that the turbulent phase is characterized by a much faster dissipation of magnetic energy than would be expected from diffusion of **B**_{α} (dashed line). In these plots, the energy is normalized by *E*_{pot}, which is the energy of a uniform vertical field **B**=*B*_{0}**e**_{z} with *B*_{0} chosen to give the same magnetic flux as **B**_{α}. This is the minimum possible energy for each configuration in our periodic domain, ignoring all helicity constraints (and also the constraint of line-tying on the side boundaries). Some of the magnetic energy is lost by ohmic dissipation, but the majority of magnetic energy is dissipated by viscous heating at shock fronts, generated by the turbulent reconnection [11]. During the turbulent phase, the rate of ohmic heating is only 20–50% that of viscous heating. This is evident in figure 4*c*, which shows the cumulative viscous and ohmic heating in each run. Both heating rates level off after the end of the turbulent phase, although the asymptotic ohmic heating rates depend on *α*, reflecting the resistive decay of the background field (which decays like

In this paper, our main focus is on the magnetic topology of the end states. Here, ‘end state’ means the gradually decaying configuration that remains after the turbulent phase has ended. It is evident from the magnetic field lines at *t*=300 (figure 2*e*–*h*) that there is a difference between the end states for *α*=0.001 and *α*=0.01, as compared with *α*=0.05 and *α*=0.1. In the former two runs, there are two oppositely twisted flux tubes, whereas in the latter two runs there is only a single flux tube. We remark that Parker [25] proposed that static MHD equilibria cannot have more complicated topology than either a single flux tube or two parallel flux tubes of opposite twist, and our end states are in accord with this.

The separation into either one or two tubes is clearly seen in figure 6, which shows the average value **j**⋅**B**/*B*^{2} along each magnetic field line. The quantity λ is the current helicity density (we avoid the symbol *α* which refers specifically to the background field **B**_{α}). In a force-free equilibrium, which approximately holds after the turbulent relaxation, λ is constant along each field line. Note that the separation into two tubes is not merely a transient phenomenon: the two twisted tubes for *α*=0.001 or *α*=0.01 actually repel one another and will not eventually merge together. Rather, their currents (twist) will continue to individually decay on the resistive timescale.

The transition between end states with single and double flux tubes occurs at a critical *α* between 0.01 and 0.05. In our case, the region of turbulence coincides with the region of field line mixing, namely the kidney-shaped region best seen in the colour maps of figure 7. The transition in the final state is triggered when a particular hyperbolic (index −1) periodic orbit in the initial state moves inside the mixing region. This hyperbolic orbit is located well outside the mixing region at (*x*,*y*)≈(−4.5,0) when *α*=0.001 (figure 7*a*), moves closer [at (*x*,*y*)≈(−3.8,0.06)] for *α*=0.01 (figure 7*d*) and is eventually inside the mixing region for *α*=0.05 (figure 7*g*). This changes the topological degree of the turbulent region from 2 to 1.

Note that there is an asymmetry in the two tubes produced by the turbulent relaxation, and this asymmetry increases as *α* is increased. This is seen by comparing figure 6*e* and 6*f* with figure 6*b* and 6*c*. Firstly, the separating motion of the tubes in the *x*-direction is influenced by the background field. (If there were no background field, the tubes would simply move apart symmetrically about *x*=0.) Note that we have repeated the simulation with a larger domain in *x* with identical results at *t*=300, confirming that the background field causes the asymmetry, rather than the numerical boundary conditions. Secondly, the pattern of reversed-sign **B**_{α} with respect to the two tubes, there is a more significant current sheet outside the left-hand tube than outside the right-hand tube, seen clearly for *α*=0.01. For *α*=0.001, the background field is too weak to produce notable asymmetries.

As *α* is increased further beyond 0.01, the separation of the two tubes becomes so small that the left-hand tube is eventually engulfed by the right-hand tube. The run with *α*=0.05 is interesting because it is just past the transition point between double and single tube final states. In this run, the initial turbulent relaxation leaves a vestige of the second tube at *t*=100 (figure 6*h*), with a strong current sheet outside it. This current sheet is sharp enough that it undergoes resistive decay by time *t*=300, removing the second tube. From this, we see that the precise location of the transition point between asymptotic states with one and two tubes is likely to be dependent on the resistivity. On the other hand, the nature of the final state of the turbulent relaxation (e.g. at *t*=100) is conjectured to be independent of the resistivity.

## 5. Discussion

This numerical experiment shows that one must choose the boundary appropriately if one is to correctly predict the end-state topology based on the topological degree of the initial state. The practical application of such a prediction is therefore dependent on being able to predict the extent of the turbulent relaxation sufficiently accurately. In our case, the region of turbulent relaxation is largely determined by pre-existing mapping complexity in the initial magnetic field. Therefore, one makes the correct prediction by considering the chaotic mixing region of the colour maps in the initial states (figure 7).

In other situations, it may be difficult to predict the extent of the turbulent region before the onset of dynamical relaxation. For example, Bareford *et al.* [11] began with a laminar magnetic field structure not containing current sheets. Only once the kink instability had led to the onset of a turbulent relaxation, did it become clear that the extent of the turbulent region would be about 1.8 times the diameter of the initial loop. It was suggested by Bareford *et al.* [11] that, due to the presence of zero net vertical current, their relaxation region was more localized than previous simulations by Browning *et al.* [9] in which turbulence filled the whole domain. However, our simulations with a net vertical current still have a localized turbulent region (e.g. the *α*=0.1 case presented here, or the ‘*S*^{3}’ case described by Wilmot-Smith *et al.* [13]).

From a practical point of view, it is very desirable to predict not only the topology (e.g. number of flux tubes) of the end state, but also the amount of magnetic energy released. A possible approach is to apply Taylor theory—assuming conservation of total magnetic helicity—restricted to the turbulent region [11]. This would predict a linear force-free field within that region. For our *α*=0.1 simulation, the field does relax to a much smoother and symmetric spatial distribution of λ. But, according to the topological degree, the cases *α*=0.01 and *α*=0.001 cannot relax to the Taylor state, and indeed this is what our simulations show. We see the formation of two separate flux tubes of oppositely signed λ. However, even in the case where the topological degree is consistent with a Taylor state, we find that the resulting flux tube is surrounded by a region of oppositely signed λ, such that a field with constant (or piecewise-constant) λ is not clearly appropriate.

The physical nature of the degree constraint is nothing more or less than the freezing-in of the magnetic topology on the side boundaries of the turbulent (non-ideal) region. This constraint will exist whenever the turbulent region is localized within a wider ideal region. In our parameter study, the transition between final states with one and two flux tubes may be thought of as a change in the dominance of the contribution to the field line mapping from **B**_{braid} compared with **B**_{α}. But ultimately, it is the initial degree of the mapping restricted to the turbulent region that constrains the evolution.

Our assumption of a periodic domain is inessential. Although the results here are presented for the case of periodic *z*-boundaries, we have repeated the simulations for line-tied *z*-boundaries (**v**=0), as would be appropriate for the fast relaxation of coronal loops in the solar atmosphere. The qualitative finding of a transition between double and single tube final states as *α* is increased remains valid. The main difference is that the two tubes for *α*=0.001 and *α*=0.01 are restricted from moving apart by the line-tying of their magnetic footpoints.

Finally, we note that, although we have illustrated with resistive-MHD simulations, the degree constraint is purely a property of the global magnetic field. It is applicable more generally, relying neither on the fluid approximation nor any particular physics assumed within individual reconnection sites.

## Data accessibility

The Lare3D source code needed to replicate our calculations is available from the CCPForge repository http://ccpforge.cse.rl.ac.uk/gf/.

## Authors contributions

A.R.Y. drafted the manuscript. A.J.B.R. and A.R.Y. ran the numerical simulations and generated the output. All authors contributed to design of the study and analysis of the results, and gave final approval for publication.

## Competing interests

We declare we have no competing interests.

## Funding

The authors were supported by STFC consortium grant nos. ST/K000993/1 and ST/K001043/1 to the universities of Dundee and Durham.

## Acknowledgements

Numerical simulations used the UKMHD cluster at the University of St Andrews, funded by STFC and SRIF. We thank the referees for useful suggestions.

## Footnotes

↵1 Lare3D is available from http://ccpforge.cse.rl.ac.uk/gf/.

- Received January 8, 2015.
- Accepted April 30, 2015.

© 2015 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.