Braided magnetic fields: equilibria, relaxation and heating

The main purpose of this paper is to investigate the importance of magnetic braiding for plasma dynamics, and the first issue we address is the existence of smooth stable equilibria for braided magnetic topologies. A key ingredient in Parker's magnetic braiding mechanism for coronal heating is the hypothesis that for sufficiently complex magnetic braids such equilibria do not exist. The coronal field is expected to evolve through a sequence of such equilibria, since the Alfvén travel time along a coronal magnetic field line is much faster than the timescale of the photospheric driving. Plasma heating, mediated by reconnection, is then postulated when the field complexity reaches a state where no (continuous) equilibrium exists. Due to the low coronal plasma-β, such equilibria approximate Beltrami or force-free fields (f.f.f.):

$$\begin{eqnarray}&&\nabla \times \mathbf{B}=\alpha \mathbf{B}.\end{eqnarray} \tag{ 4 }$$

Little has so far been unequivocally proven regarding the existence of such f.f.f.s in non-symmetric configurations. One exception is the study of Bineau (1972), in which it was proved that in the absence of magnetic null points, f.f.f.s do exist for sufficiently small α (though no bounds are provided on α). There exist a number of different approaches to the problem, and—while there is mounting evidence in favour of thin but finite current layers as opposed to the tangential discontinuities proposed by Parker—there is presently no clear consensus, see Craig and Sneyd (2005), Pontin and Hornig (2015) and references therein. We refer the reader here also to the work on equilibria that include flow and plasma pressure, and can be modelled by so-called double-Beltrami fields (e.g. Yoshida et al 2001).

We consider in this paper a model braided magnetic field

$$\begin{eqnarray}&&\mathbf{B}=\underset{i=1}{\overset{6}{\sum}}\,\kappa {{\left(-1\right)}^{i}}\exp \,\left(-\frac{{{\left(x-{{x}_{i}}\right)}^{2}}+{{y}^{2}}}{2}-\frac{{{\left(z-{{z}_{i}}\right)}^{2}}}{4}\right)\\ &&\times \left(-y\,{{\mathbf{e}}_{x}}+\left(x-{{x}_{i}}\right)\,{{\mathbf{e}}_{y}}\right)~+~1{{\mathbf{e}}_{z}}\end{eqnarray} \tag{ 5 }$$

where ${{x}_{i}}={{\left(-1\right)}^{i+1}}$ and z_1..6  =  {−20, −12, −4, 4, 12, 20}. When the parameter $\kappa =1$, particular field lines in the domain (figure 1(a)) have a 'pigtail braid' structure (figure 1(b)). Such braiding could be generated in an initially uniform magnetic field via a sequence of rotational plasma motions on the boundary in a 'blinking vortex' pattern (Aref 1984). However, the particular pattern of this braid is not crucial to the qualitative results presented. Rather, the salient feature is the generic property of tangling of magnetic field lines that is well established in laboratory plasmas and has recently begun to be quantified in solar coronal observations (Yeates et al 2012, 2014).

Zoom In Zoom Out Reset image size

Here, we demonstrate for the first time that the topological entropy of this field increases as the parameter κ is increased, confirming that fields with higher κ are 'more braided' in a rigorous sense. The calculation of topological entropy is performed by directly evaluating the growth of a material line under successive iterations of the field line mapping (Newhouse and Pignataro 1993), using an adaptive procedure in which additional points are included along the line in regions of high curvature, as required to resolve this length. The result is plotted for different values of κ in figure 2. Fitting a straight line to the plot, we obtain a scaling

$$\begin{eqnarray}&&h\sim \kappa \,\left(2.76\pm 0.09\right)+0.39\pm 0.11\end{eqnarray} \tag{ 6 }$$

(the errors being the 95% confidence interval). That is, increasing κ corresponds to increasing the overall field complexity as quantified by h. We return to consider this scaling later.

Zoom In Zoom Out Reset image size

This section describes the results of a numerical relaxation method for investigating the existence of braided equilibria. The philosophy is to define a magnetic field with a given topology, and then relax towards an equilibrium while exactly preserving that topology. In order to guarantee topology preservation a Lagrangian numerical mesh is employed, as described in Craig and Sneyd (1986) and Candelaresi et al (2014). A magnetofrictional (MF) evolution enforces a monotonic decay of the magnetic energy in the system, as a f.f.f. is approached. The critical question is whether this equilibrium is smooth or contains tangential discontinuities, i.e. current sheets. If it is discontinuous then the finite dissipation in a real plasma would lead to reconnection, topology change, and energy release. A number of studies using this approach have concluded that there is no evidence of formation of tangential discontinuities except in the presence of magnetic null points that we do not treat here (Craig and Sneyd 2005, Candelaresi et al 2015). In particular the magnetic field of equation (5) has been used as an initial condition for such ideal MF relaxations, first by Wilmot-Smith et al (2009) and subsequently by Candelaresi et al (2015) using an improved numerical scheme. In neither case was any tendency towards the formation of tangential discontinuities observed. Of course no numerical relaxation can reach an exact equilibrium. Thus, here we describe an extension of the MF relaxation. Specifically, we take the approximate equilibrium obtained after MF relaxation, and further relax this field using an MHD code (Galsgaard and Nordlund 1996) with a large viscosity applied. In this MHD code the resistivity is set explicitly to zero—due to the sixth-order spatial differencing numerical dissipation (and thus topology change) is minimised, so long as length scales of variations in $\mathbf{B}$ remain well above the grid-scale.

We have performed a systematic study varying the parameter κ in equation (5) using this approach, with both line-tied (perfectly conducting) and periodic boundary conditions in z (for spatial domain $x,y\in \left[-6,6\right]$ and $z\in \left[-24,24\right]$). The results for the line-tied case were presented by Pontin and Hornig (2015), who obtained a sequence of ideal equilibria for $\kappa <0.75$. The residual $\mathbf{J}\times \mathbf{B}$ forces present in the approximate equilibria could be balanced by plasma pressure for a plasma-β of order 10⁻³—thus the equilibria are as force-free as one would expect in the corona or a typical laboratory plasma. Figure 3 shows a sequence of the equilibria obtained, corresponding to different values of κ. The complexity of these magnetic braids increases for increasing κ, as shown by the topological entropy calculated above, which does not change during the ideal relaxation. These equilibria—that we emphasise are for the line-tied case—contain thin but well-resolved current layers (figure 3, middle panels). The equilibria share the same topology as the magnetic field (5) for values of the parameter κ up to $\kappa =0.75$.

Zoom In Zoom Out Reset image size

These results should be contrasted with the newly considered case where the z-boundaries are periodic. The simulations that we have performed using the same method but with periodic boundaries show no ideally accessible equilibrium for the range $0.5\leqslant \kappa \leqslant 0.75$. This non-existence of an equilibrium manifests through the generation of current layers that continue to thin until they reach the grid scale for all accessible numerical resolutions—this leads naturally to reconnection via numerical dissipation. The reason for this absence of an equilibrium for the periodic magnetic braids is explained in the following section.

In summary, for periodic magnetic braids of the form in equation (5), no equilibrium exists. For line-tied magnetic braids, a sequence of braided magnetic fields that are continuous equilibria does exist. However, these equilibria contain thin current layers, whose thickness is directly related to the field complexity, becoming exponentially thinner as the parameter κ is increased. We anticipate that this result for the line-tied case is a general one, not limited to the class of braids considered here, for the following reasons. First, the results are consistent with a number of other studies in which a magnetic field between two line-tied plates was subjected to sequences of shears, and an exponential thinning of current layers was observed as the number of shears (and the resulting field complexity in the volume) increased (van Ballegooijen 1988a, 1988b, Mikić et al 1989, Candelaresi et al 2015). Second, it can be argued that braided equilibria containing thin layers in the field line mapping must exhibit thin current layers on the same (or smaller) scale. This is explicitly demonstrated in the next section.

3.4.1. Exact equilibria.

To further explore lengthscales in f.f.f.s, let us consider the case of an exact f.f.f. $\nabla \times \mathbf{B}=\alpha \mathbf{B}$. For the periodic case any magnetic field containing ergodic field lines is consistent only with a linear f.f.f., or piecewise-linear in separate ergodic regions. For our model field (5), if we identify the top and bottom z-boundaries with one another to make the field periodic, then for $\kappa >0.5$ the region $\sqrt{{{x}^{2}}+{{y}^{2}}}\lesssim 2$ is highly mixing and Poincare maps show a volume-filling ergodic region. Thus, in this region α must be constant. Now, for a constant-α f.f.f. $\mathbf{J}=\alpha \mathbf{B}$ and the current helicity is ${{H}_{J}}={\int}_{V}\mathbf{J}\centerdot \mathbf{B}\,{{\text{d}}^{3}}x=\alpha {\int}_{V}|\mathbf{B}{{|}^{2}}\,{{\text{d}}^{3}}x$. In our case by symmetry we have H_J  =  0, which implies that $\alpha =0$. However, this is a contradiction, since it is inconsistent with the presence of ergodic field lines, and we conclude that an exact equilibrium for this field topology for the periodic case does not exist.

Consider now the requirements for an exact equilibrium in the case where $z=\pm {{z}_{0}}$ are line-tied, perfectly conducting boundaries. Then α is constant along field lines that run between these boundaries. To understand the nature of the distribution of α perpendicular to $\mathbf{B}$, consider the following. Let us suppose that α has some global length scale, $\mathcal{O}\left(\ell \right)$, on some plane, say z  =  −z₀. Then, to find the length scales associated with the α distribution at, say, z  =  +z₀, we simply map α along field lines. We can relate the resulting length scales in the mapped quantity to the topological entropy, h, as follows. Recall that h may be interpreted as the limit of the stretching rate of material lines under repeated application of the mapping, the lengths of these lines L satisfing $L\sim \exp (hn)$ (n being the iteration number).

Consider a contour of α as such a material curve that is mapped in the 'flow' of $\mathbf{B}$. Then by definition the contour is stretched by a factor $\approx \exp (h)$ between z  =  −z₀ and z  =  z₀. Assuming constant area along the length of the braided flux tube, there must be a squeezing on average by a factor $\exp \left(-h\right)$ in the orthogonal direction. Thus, a typical contour of α on z  =  +z₀ will form an elongated structure with characteristic thickness $~{{L}_{\alpha}}\sim \ell \exp \left(-h\right)$. If the flux tube expands/contracts along its length then the area changes by a factor $B_{n}^{-}/B_{n}^{+}$, the ratio of the normal field components at the two boundaries for the elemental flux tube under consideration. Unless the braid has a high contraction/expansion factor this ratio will be $\mathcal{O}\,(1)$. Assuming an isotropic expansion, the length scales of α are modified to

$$\begin{eqnarray}&&{{L}_{\alpha}}\sim \ell \,{{\text{e}}^{-h}}{{\left(\frac{B_{n}^{+}}{B_{n}^{-}}\right)}^{-1}}.\end{eqnarray} \tag{ 7 }$$

Now, in a f.f.f., ${{J}_{\parallel}}\propto \alpha$, and thus the equilibrium must exhibit current layers that are also on the length scale ${{L}_{\alpha}}$ (or smaller). Note that Pontin and Hornig (2015) made a similar argument based on the squashing factor Q (see equation (1)), which shows that ${{L}_{\alpha}}\sim {{\left(Q\,B_{n}^{+}/B_{n}^{-}\right)}^{-1/2}}$.

Pontin and Hornig (2015) showed that the current layer thickness w_J in the equilibria for line-tied boundaries follows the scaling

$$\begin{eqnarray}&&{{w}_{J}}=\exp \left[\kappa \left(-5.59\pm 0.25\right)+2.53\pm 0.18\right].\end{eqnarray} \tag{ 8 }$$

Comparing this with the empirical scaling of ${{L}_{\alpha}}$ with h that would be obtained from equations (6) and (7), we conclude that equation (7) gives a rather weak upper bound on the current layer thickness—the current layers are in fact significantly thinner than the maximum permitted by the entropy argument. This difference could be down to several factors. First, h measures a global or average line stretching rate (and the argument above treats a global-scale α contour), while the current layers form in the vicinity of selected field lines—we assumed an equal squeezing of the α contour along its length to arrive at the factor $\exp \left(-h\right)$ while in practice this will vary between field lines. Second, h is determined by repeated iterations of the field line mapping, whereas in these line-tied simulations a material line or α contour is subject only to a single mapping iteration. The scaling of the current layer thickness in the domain can be shown to match very well with the scaling of the thinnest layers of Q present in the mapping (Pontin and Hornig 2015)—this quantity being derived from a single iteration of the mapping and being minimised over field lines (over x and y), thus removing the above two concerns.

3.4.2. Approximate equilibria.

We addressed above the nature of an exact equilibrium—however, any real (or simulated) magnetic field will only ever approximate such an exact f.f.f. Thus the question arises; what is the influence of a small departure from equilibrium on these results? To answer this we define

$$\begin{eqnarray}&&{{\alpha}^{\ast}}:=\frac{\mathbf{J}\centerdot \mathbf{B}}{{{B}^{2}}}=\frac{{{J}_{\parallel}}}{B};\quad {{\epsilon}^{\ast}}:=\frac{\parallel \mathbf{J}\times \mathbf{B}\parallel}{\parallel \mathbf{J}\parallel \parallel \mathbf{B}\parallel}=\frac{\parallel {{\mathbf{J}}_{\bot}}\parallel}{\parallel \mathbf{J}\parallel}.\end{eqnarray} \tag{ 9 }$$

The first quantity, ${{\alpha}^{\ast}}$, in general varies along field lines, but converges to the f.f.f. parameter α as we approach a force-free equilibrium. The second quantity is a dimensionless function that measures the proximity to a f.f.f. We have $0\leqslant {{\epsilon}^{\ast}}\leqslant 1$, where ${{\epsilon}^{\ast}}=0$ corresponds to an exact f.f.f. and ${{\epsilon}^{\ast}}=1$ to a maximal non-force-free state. The above argument for the non-existence of force-free equilibria in a periodic domain relied on the fact that α is constant along field lines. However, a similar argument holds for a near-force-free state if ${{\alpha}^{\ast}}$ has only small variations over the length of the periodic domain, or in other words if the correlation length of ${{\alpha}^{\ast}}$, l_corr is longer than the domain length L. We define the correlation length for a given field line by

$$\begin{eqnarray}&&\frac{1}{{{l}_{\text{corr}}}}:=\langle |\frac{1}{{{\alpha}^{\ast}}}\frac{\partial {{\alpha}^{\ast}}}{\partial s}|\rangle,\end{eqnarray} \tag{ 10 }$$

where s parameterises length along the field line and $\langle.\rangle$ is an average over the field line. Following a similar calculation as in Pontin et al (2009) we use

$$\begin{eqnarray}&&0=\nabla \centerdot \mathbf{J}=\nabla \centerdot \left({{J}_{\parallel}}{{\mathbf{e}}_{B}}\right)+\nabla \centerdot {{\mathbf{J}}_{\bot}}=\mathbf{B}\centerdot \nabla {{\alpha}^{\ast}}+{{\nabla}_{\bot}}\centerdot {{\mathbf{J}}_{\bot}},\end{eqnarray} \tag{ 11 }$$

which leads to

$$\begin{eqnarray}&&|\frac{1}{{{\alpha}^{\ast}}}\frac{\partial {{\alpha}^{\ast}}}{\partial s}|=|\frac{1}{{{J}_{\parallel}}/B}\frac{\mathbf{B}\centerdot \nabla {{\alpha}^{\ast}}}{B}|=|\frac{-{{\nabla}_{\bot}}\centerdot {{\mathbf{J}}_{\bot}}}{{{J}_{\parallel}}}|\approx \frac{{{\epsilon}^{\ast}}}{d},\end{eqnarray} \tag{ 12 }$$

where we used in the last step that for approximately force-free states $|{{J}_{\parallel}}|\approx \parallel \mathbf{J}\parallel$ and introduced a typical length scale d for variations perpendicular to the field line: $\parallel {{\mathbf{J}}_{\bot}}\parallel /d\approx |-{{\nabla}_{\bot}}\centerdot {{\mathbf{J}}_{\bot}}|$. The result, ${{l}_{\text{corr}}}\approx \langle d/{{\epsilon}^{\ast}}\rangle$ shows that even for high-quality force-free approximations with ${{\epsilon}^{\ast}}\approx {{10}^{-3}}$ (as in our examples) the correlation length can be comparable with or smaller than the length of the domain (48 units in our case) if the perpendicular length scales become small as well ($d\approx 0.05$ for k  =  1 in equation (8)). One can turn this argument around to make a statement about the structure of approximate equilibria with ergodic field lines. Specifically, an approximate f.f.f. with ergodic field lines in a periodic domain of length L must have ${{l}_{\text{corr}}}\leqslant L$. This allows us to place an upper bound on the smallest perpendicular length scale for field lines in the ergodic domain ${{d}_{\text{min}}}\leqslant \langle d/{{\epsilon}^{\ast}}\rangle \langle {{\epsilon}^{\ast}}\rangle \leqslant L\langle {{\epsilon}^{\ast}}\rangle$, hence the field must contain current sheets of thickness $L\langle {{\epsilon}^{\ast}}\rangle$ or smaller. We note that this argument makes concrete the assertions of Taylor (1993) regarding filamentation of the current in the relaxed state of a tokamak plasma. Note further that Candelaresi et al (2015) have examined the correlation length in numerical relaxation simulations, and found that even when it becomes large on average, it may remain relatively small at areas of high complexity, hindering the numerical relaxation process.

The results described in the above section demonstrate that any sufficiently complex magnetic braid will be unstable in a resistive plasma owing to the development of thin current layers in the (approximate) equilibrium. Thus, in this section we describe the evolution following the onset of reconnection. This evolution was first followed in resistive MHD simulations with line-tied boundaries, starting from the magnetic field (5) (with $\kappa =1$), by Wilmot-Smith et al (2010) and Pontin et al (2011). It has subsequently been shown that repeating these simulations with periodic boundaries in the z-direction does not effect the qualitative nature of the relaxation although periodic boundaries do allow the field to reach a slightly lower energy end state (Yeates et al 2015). In this section we make a spectral analysis that demonstrates the turbulent properties of the relaxation for the first time and we briefly summarise the main features of the evolution and final state. The analysis that we present builds on the simulation of Pontin et al (2011), as well as two new simulations with the same setup and method, but with $\kappa =0.7$ and $\kappa =1.5$, thereby exploring for the first time the impact of the braid complexity. For concreteness we describe results of simulations with line-tied boundaries—however, the general results apply equally well to the periodic case.

The resistive relaxation of our magnetic braid involves initially the formation of two current layers (those that form under the ideal relaxation), that due to the finite resistivity undergo reconnection, leading to a cascade of further reconnection events. The global evolution of the system is characterised by an array of current layers that form and then dissipate, efficiently filling the volume of braided field lines—see figure 4(b). The dominant field component in the current layers is the 'guide' field along the tube (in the z-direction). Therefore the reconnection in these current layers takes the form of a continuous 'flipping' or 'slippage' of field lines (Priest and Démoulin 1995, Hornig and Priest 2003). The net result of these many localised reconnection events (around 30 when $\kappa =1$ and the magnetic Reynolds number ${{R}_{\text{m}}}={{10}^{-3}}$) is that the magnetic field 'unbraids', leading to a much simpler topological state (figure 4(c)).

Zoom In Zoom Out Reset image size

The resistive relaxation process described above has many of the properties of decaying turbulence. In figure 5 we plot the 2D magnetic energy spectrum over the xy-plane (averaged over z) for two of the simulation runs of Pontin et al (2011). We have an estimated spectral index close to  −5/3, which is within the range expected for decaying 2D MHD turbulence (Biskamp and Schwarz 2001). We also see that as R_m is increased we obtain a larger inertial range over which this scaling is obeyed, as expected. Together with this magnetic energy spectrum, a similar spectrum of kinetic energy is observed. The small-scale flows generated during the turbulent phase organise into decaying large scale flows in the final quasi-equilibrium (described below). This is consistent with the 'reverse dynamo' flow generation mechanism of Mahajan et al (2005).

Zoom In Zoom Out Reset image size

A number of factors control how fully developed the turbulence becomes—one being the value of R_m as shown in figure 5. In addition, analysing the simulations with different values of κ (the initial field complexity), we find that the turbulence becomes more developed and the overall relaxation timescale increases as κ is increased. This is illustrated in figure 6 by isosurfaces of the current density for two resistive relaxation simulations with different values of κ. The plots are made at a time when a maximal number of discrete current layers appear in the volume—at t  =  55 for $\kappa =1.5$, and t  =  45 for $\kappa =0.7$. Note for reference that the Alfvén travel time along the length of the simulation domain in z is approximately 48 time units. We observe that the number of small scale current layers is greater for greater κ. The time taken to reach the final relaxed state (the nature of which is discussed in the next section) is estimated to be t  =  120 for $\kappa =0.7$, t  =  200 for $\kappa =1$, and $t\gtrsim 400$ for $\kappa =1.5$ (or approximately 2.5, 4 and 8 crossing times, respectively).

Zoom In Zoom Out Reset image size

Taylor's relaxation hypothesis posits that the only conserved topological quantity during any sufficiently turbulent relaxation is the total magnetic helicity, and a variational calculation shows that the final state should be the linear force-free field (Woltjer 1958, Taylor 1974). The theory was very successful in predicting the relaxed state in a reverse-field pinch and to a lesser extent in other devices (Ortolani and Schnack 1993), although more recently some departures of the final equilibrium from the 'Taylor state' have been noted experimentally (Cappello et al 2008). It has also been proposed that the hypothesis may apply in the Sun's corona (Heyvaerts and Priest 1984).

A major result of the study of Pontin et al (2011) was that the Taylor (linear f.f.f.) state is not reached after the turbulent relaxation of the magnetic braid (5). Since the net magnetic helicity is zero, Taylor's hypothesis would predict the uniform field $\mathbf{B}=1{{\mathbf{e}}_{z}}$ as the final state. However, what is found instead is a non-linear f.f.f. containing a pair of oppositely-twisted flux tubes, embedded within the uniform background field. This is shown by the plots of ${{\alpha}^{\ast}}$ (see equation (9)) in figure 4 (recall that for an exact f.f.f., this coincides with the force-free parameter α). The fact that the equilibrium is a non-linear f.f.f. implies the presence of additional constraints on the relaxation.

One such constraint, identified by Yeates et al (2010), is the topological degree (Polymilis et al 2003) of the field. This quantity is preserved during the relaxation so long as the topological degree of the boundary remains unchanged (as for example when the turbulent relaxation is confined to the interior) and it may render the Taylor state inaccessible. For the magnetic braid of equation (5) this topological degree—given by the sum of the indices of periodic orbits of the field—is 2, which is inconsistent with the predicted Taylor state. A further promising tool is the field line helicity, defined as in equation (3). Unlike H, its distribution uniquely characterises the topology of the magnetic field (Yeates and Hornig 2013), so understanding its evolution during turbulent relaxation could in principle allow prediction of the topology of the final state. Recently, Russell et al (2015) derived an evolution equation for the field line helicity, and analysed the importance of its different terms for slipping reconnection in a complex magnetic field. The reconnection primarily rearranges line helicity while having a relatively small impact on the total helicity, providing a new justification that relaxation approximately conserves total helicity while changing how the helicity density is distributed. However, the relaxation cannot exchange field line helicity arbitrarily, which means the Taylor state may be prohibited in some cases, even when the topological degree would permit it. This could potentially explain why the contour map of ${{\alpha}^{\ast}}$ shown in figure 4(c) is not piecewise-constant, contrary to what would be expected if locally Taylor states (linear f.f.f.s) were produced in each of the flux tubes.