Physical role of topological constraints in localised magnetic relaxation

Predicting the final state of turbulent plasma relaxation is an important challenge, both in astrophysical 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, while 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.


I. INTRODUCTION
Self-organization of turbulently relaxing plasma to a predictable minimum-energy state has been observed in laboratory confinement devices including the reversedfield pinch and the spheromak [1][2][3][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][8][9][10][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]. * anthony.yeates@durham.ac.uk The topological degree (defined in Section II) is conserved provided that the degree of the boundary does not change. The latter can be ensured by having a turbulent dynamics which is localized in the interior of the domain and does 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 reversedfield pinch, there are no conducting walls to define the relaxation volume [16]. Typically, coronal energy releasesfor 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 relax-ation 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 Section III A, 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.

II. 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 end-point 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 Yeates & Hornig [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 around the boundary of D [20]. Since 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 utilize the convenient color 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 B = ∇ × Ae z + e z , A = 0.6 sin 2 x cos(0.5y) + cos(0.3x) cos(0.3y). (2) The color map assigns one of four colors (in this paper, we use shades of gray) to each point (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 colors intersect. Furthermore, the topological degree T (D) of a region D ⊂ D 0 may be identified by noting the counterclockwise sequence of colors around the boundary of D. In particular, the number of times that the full sequence of four colors (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 counter-clockwise 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.

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 nonzero 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 initialise 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 where 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 J 0 (α r √ 128) = 0. (This is when the first field reversal occurs at the corners of the domain.) It should also be noted that B α leads to a net electric current in the z-direction.
The braiding field B braid was introduced by [22]; its construction is based on the pigtail braid, with six toroidal rings of flux, The parameters used are 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 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 [23] 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 Here ρ 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 γ = 5/3 the ratio of specific heats. Details of the numerical methods are given by Arber et al. [24]. The viscous term ∇σ in (6) includes no background viscosity, but only a shock viscosity to prevent unphysical oscillations (this takes the form given in [11]). There is a corresponding heating term εσ in (8). 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].

IV. 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. [25] 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(b), 4(c) and 5. For example, the total kinetic energy E kin = v /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 rings locations and counter-propagating in z.
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 normalised 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 during the turbulent phase the rate of ohmic heating is only 20−50% that of viscous heating. Thus the majority of magnetic energy is dissipated by viscous heating at shock fronts, generated by the turbulent reconnection [11]. Correspondingly, Figure 4(c) shows that the cumulative internal energy E int increases faster during the turbulent phase. Subsequently the rate is larger for the simulation with α = 0.1, owing primarily to the volume ohmic heating resulting from resistive decay of the background field.
In this paper, our main focus is on the magnetic topology of the end states. Here "end state" means the grad- tubes, while in the latter two runs there is only a single flux tube. The separation into either one or two tubes is clearly seen in Figure 6, which shows the average valueλ of λ = 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 kidneyshaped region best seen in the color 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 7a), moves closer [at (x, y) ≈ (−3.8, 0.06)] for α = 0.01 (Figure 7b), and is eventually inside the mixing region for α = 0.05 (Figure 7c). This changes the topological degree of the turbulent region from 2 to 1.
Notice 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 panels (f) and (j) in Figure 6 with panels (e) and (i). Firstly, the separating motion of the tubes in the xdirection 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λ around each tube is different. Owing to the direction of rotation of 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 noticeable 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 6g), 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.

V. 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 color 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 localised 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.