1 Introduction

A long-standing debate into the nature of equilibrium magnetic fields in astrophysical plasmas was triggered by Parker (1972). He proposed that if an equilibrium magnetic field is subjected to an arbitrary, small perturbation, then the resulting magnetic field will in general not relax towards a smooth equilibrium, if the plasma in which the magnetic field is embedded is ideal so that magnetic reconnection is prohibited (the absence of additional perturbations and the presence of some damping mechanism are also assumed). He proposed, rather, that the field should relax towards a state containing tangential magnetic field discontinuities—thus leading to the onset of rapid magnetic reconnection and energy dissipation. Parker termed this process “topological dissipation”, since its onset is governed by the field line topology. Parker’s hypothesis—stated more precisely in Sect. 2—has been referred to variously as the “Parker problem” and “The Parker Magnetostatic Theorem”, amongst other names. It was originally put forward to explain the absence of observed large-magnitude, small-scale magnetic fields in turbulent plasmas such as the solar photosphere.

These days, Parker’s topological dissipation mechanism is best known as a key ingredient in his proposed mechanism for explaining the heating of the solar corona by “nanoflares” (Parker 1979, 1988). The nanoflare mechanism was proposed as an alternative to various wave-damping processes, and remains one of the most promising proposed solutions to the long-standing problem of coronal heating, first appreciated in the early 1940s. (For a history of the discovery of the hot corona see the discussions and references in Peter and Dwivedi (2014) and  Russell (2018), and for reviews of our current understanding of the topic see De Moortel and Browning (2015) and references therein.) There is now little doubt that the coronal plasma is heated due to the turbulent motions in the convection zone—this energy being transmitted up to the corona by the magnetic field before being dissipated. Parker argued that, because the photospheric flows tend to be slow (\(v\sim 1{-}10\) km/s) compared to the communication speed along coronal magnetic field lines (\(v_A\sim 1000\) km/s), the coronal field should evolve through a sequence of (approximate) equilibria. Due to the small plasma pressure, these equilibria are approximate force-free fields satisfying \((\nabla \times {\mathbf{B}})\times {\mathbf{B}}\approx \mathbf{0}\). Now, the idea is that as the field line footpoints are shuffled around by the photospheric flows, the field topology will periodically reach a state in which no smooth equilibrium exists. Rather, tangential discontinuities of the field ‘spontaneously’ form, leading to the onset of energy dissipation. Each such energy dissipation event is termed a “nanoflare”, with the X-ray corona being formed by a superposition of many such nanoflares occurring throughout the Sun’s atmosphere. Parker (1988) suggested that a typical nanoflare should release energy of order \(10^{24}\)ergs or less, with the largest having size less than \(10^{27}\)ergs. Since Parker’s original idea the term “nanoflare” has come to refer to any small, (observationally) unresolved heating event within a coronal loop, and substantial efforts have been made to infer the spectrum of nanoflare energies required to explain the observations (e.g., Reale 2014; Cargill et al. 2015).

In this review we discuss the debate surrounding the topological dissipation mechanism since it was first proposed by Parker (1979), and the current state of our understanding. It is worth noting that analogous problems to the existence of smooth force-free fields exist in other fields of research such as fluid dynamics, with vector fields satisfying \((\nabla \times {\mathbf{B}})\times {\mathbf{B}}= \mathbf{0}\) often referred to as “Beltrami” fields (see Sect. 4.1). We also discuss more broadly modelling approaches to understanding current sheet formation in response to random footpoint motions, and implications for coronal heating. The article is organised as follows. In Sect. 2 we state “The Parker Problem”. In Sect. 3 we discuss the energy injection process by footpoint motions. In Sects. 4 and 5 we present theoretical and computational work on the existence of force-free fields of arbitrary topology, and various approaches to demonstrating the formation (or otherwise) of tangential discontinuities in the geometry of the Parker problem. We go on in Sect. 6 to discuss the energy release process when magnetic reconnection is permitted to occur in the domain, while in Sect. 7 we consider formation of current sheets in the complex magnetic field topology of the corona—when magnetic field structures excluded in Parker’s original analysis are included. We finish in Sect. 8 with a summary.

2 Statement of the Parker hypothesis

The idea that led to the genesis of “The Parker Problem” is most strongly and precisely stated by Parker (1994), as the “magnetostatic theorem” (see also the discussion of Low 2010b). To begin, consider a plasma in force-balance between pressure forces and magnetic stresses, per

$$\begin{aligned} {\nabla p} = \frac{1}{\mu _0}(\nabla \times \mathbf{B})\times \mathbf{B}. \end{aligned}$$
(1)

Now, consider such a field in an infinite volume, \(\mathcal {V}\), between two parallel, perfectly-conducting planes, in which all magnetic field lines connect from one plane to the other. (The normal magnetic field is positive-definite on both planes and no closed field lines exist in \(\mathcal {V}\).) The plasma within \(\mathcal {V}\) is sufficiently highly-conducting that it can be treated as a perfectly conducting plasma that evolves according to

$$\begin{aligned} \frac{\partial \mathbf{B}}{\partial t}=\nabla \times (\mathbf{v}\times \mathbf{B}). \end{aligned}$$
(2)
Fig. 1
figure 1

Sketch by E. N. Parker showing the “arbitrary winding of the field lines” in the geometry of the Parker Problem, beginning with a the uniform field \(B_0\), with the field lines becoming “mixed and interlaced” after a time t in b. Image reproduced with permission from Parker (1994), copyright by OUP

The initial equilibrium satisfying (1) is usually taken to be the homogeneous field \(\mathbf{B}=B_0\mathbf{e}_z\) and the perfectly conducting boundaries are planes of constant z. We now consider the perturbation of this equilibrium field by tangential motions on one or both boundaries. These motions are taken to be smooth and differentiable (in order to represent, say, flows on the solar photosphere). Since the boundaries are perfectly conducting, the footpoints of the magnetic field lines are frozen to plasma elements there. Thus, if the boundary motions are disordered or ‘random’, the magnetic field lines in \(\mathcal {V}\) become tangled (see Fig. 1). Some time \(t=\tau \) later, suppose that the boundary flows stop. The magnetic field in \(\mathcal {V}\) is now out of equilibrium: we now allow the field to relax towards a static equilibrium, while holding the magnetic field line footpoints on the boundaries fixed. Since the plasma is perfectly conducting there is no magnetic reconnection in \(\mathcal {V}\). Thus the topology of the field has to remain the same (i.e., the magnetic field line connectivity between boundary points is preserved, as is the relative winding of any pair of field lines). If, however, no smooth equilibrium of a given topology exists then the process of relaxation must lead to some form of non-smoothness in the field. Typically in a plasma this non-smoothness corresponds to a tangential discontinuity, that is a jump of the tangential component of the magnetic field across a surface. This corresponds to a singular current along the surface and hence the name current sheet is used interchangeably. However, care has to be taken as for many authors a current sheet could also be of finite thickness (current layer).

Definition 1

(Parker hypothesis, general form) In the above setting the magnetostatic theorem (Parker 1994) states that in almost all cases (i.e., for almost all possible boundary flows) the magnetic field develops tangential discontinuities in \(\mathcal {V}\) during the relaxation to static equilibrium, due to the absence of a static equilibrium that is smooth.

This conjecture by Parker has never been proven or disproven, leading to the so-called “Parker Problem”. One of the principal mathematical difficulties arises from the need to precisely specify and maintain the magnetic field topology (field line tangling), as pointed out, for example, by Low (2010b) and Janse et al. (2010). Note that due to the typically small value of the plasma-\(\beta \) in the solar corona, the pressure gradient is often omitted from Eq. (1), meaning that the static equilibrium corresponds to a force-free magnetic field. This gives rise to a stronger version of the above definition:

Definition 2

(Parker hypothesis, force-free form) For almost all possible boundary flows the magnetic field develops tangential discontinuities in \(\mathcal {V}\) during relaxation to a force-free equilibrium.

For almost all practical applications we are unable to distinguish between the development of a singular current sheet (tangential discontinuity) and a sufficiently narrow and strong current layer. Indeed an exponentially growing current would provide a similarly good argument for the onset of magnetic reconnection in a plasma as the absence of a smooth force-free equilibrium. Thus one could weaken the above statement to

Definition 3

(Parker hypothesis, weak form) For almost all possible boundary flows the magnetic field develops current layers in \(\mathcal {V}\) during relaxation to a (force-free) equilibrium. The width of the current layers decreases and its strength increases exponentially with the complexity of the deformation.

This definition requires that we specify what precisely is meant by the complexity of the deformation. We will make this more precise in Sect. 5.3. For the moment we just assume that there is a measure of complexity of the boundary motion defining the deformation, which determines how strong and narrow the current layers become. The main difference to Definition 2 is that they never become singular.

The remainder of this article discusses various aspects of Parker’s mechanism for coronal heating by magnetic braiding. Note that for the majority of this review we confine our discussions to the geometry originally considered by Parker (as described above), in which all magnetic field lines connect from one perfectly-conducting plane to another. In spite of its apparent simplicity, this geometry has proved the principal focus of studies in the literature, as discussed below. We emphasise that in this geometry the existence of magnetic nulls, bald patches, and their associated separatrix surfaces is strictly excluded. The question of current sheet formation in magnetic fields that include such additional topological complexity is discussed briefly in Sect. 7.

In Sects. 4 and 5 we look at the different ways in which the above three hypotheses have been tested. First, in the next section we discuss the background in the context of the coronal heating problem.

3 Energy injection

In order to maintain the hot corona, sufficient energy must be supplied to balance the losses due to conduction and radiation. These losses are estimated at around \(100\mathrm {\ Wm}^{-2}\) for the quiet Sun, up to around a few times \(10^4\mathrm {\ Wm}^{-2}\) in active regions (Withbroe and Noyes 1977; Kano and Tsuneta 1996; Schonfeld et al. 2017). It is generally accepted that the energy source that maintains the corona against these losses is the solar convection, although how exactly that energy is transported to the chromosphere and corona and subsequently converted to heat remains a matter of debate. A clue that the magnetic field plays an important role comes from the observed correlation between magnetic field strength and temperature within coronal loops (Fisher et al. 1998; Pevtsov et al. 2003). What’s more, the hottest loops in active regions appear to be those that have the highest field strengths (being anchored at one end in a sunspot umbra) while also being subject to perturbations from convective motions (being anchored at the other end in a penumbral or plage region), as shown by Tiwari et al. (2017).

As described above, Parker (1972) proposed that the convective motions on the solar surface tangle the field lines in the corona around one another, increasing the magnetic energy which could then subsequently be dissipated as heat. The free magnetic energy that builds up in the corona is associated with magnetic field components (\(B_\perp \)) perpendicular to the ‘axial’ direction along the coronal loop. Parker (1988) considered a shear flow applied to a coronal loop, taking typical values for the magnetic field strength, flow velocity, and loop length of 100 G and 0.5 km s\(^{-1}\), and 100 Mm, respectively. He argued using energy balance that the average coronal value of the \(B_\perp \) field component should be of the order 25% of the ‘axial’ field along the loop, giving the ‘Parker angle’ of field lines of around \(14^\circ \) to the vertical. It should be noted that this estimate is based on an extremely simplified flow geometry, and can be expected at best to give a domain-averaged estimate of \(B_\perp \).

The flow of magnetic energy into (or out of) the solar atmosphere is quantified by the Poynting flux through the solar surface, obtained by integrating the flux of the Poynting vector

$$\begin{aligned} {\mathbf{F}}=\frac{1}{\mu _0}{\mathbf{E}}\times {\mathbf{B}}=-\frac{1}{\mu _0}({\mathbf{v}}\times {\mathbf{B}})\times {\mathbf{B}}=-\frac{({\mathbf{v}}\cdot {\mathbf{B}}){\mathbf{B}}}{\mu _0}+\frac{B^2\,{\mathbf{v}}}{\mu _0} \end{aligned}$$
(3)

through the surface, the above equation assuming an ideal evolution on the surface. The second term on the right-hand side is associated with flux emergence/submergence, which is excluded from Parker’s model; assuming zero flow through the solar surface S then the Poynting flux is

$$\begin{aligned} F_P=\int _S {\mathbf{F}}\cdot d\mathbf{A}=-\int _S\frac{1}{\mu _0}({\mathbf{v}}_t\cdot {\mathbf{B}}_t)\, B_n\, dA , \end{aligned}$$
(4)

where \(B_n\) is the normal component of \({\mathbf{B}}\) to the surface and \({\mathbf{v}}_t,{\mathbf{B}}_t\) are the tangential components of the flow and magnetic field, respectively, and the normal vector to the solar surface is taken to be upwards. Clearly, this flux may be either positive or negative depending on the relative orientations of \({\mathbf{v}}_t\) and \({\mathbf{B}}_t\).

One of the inherent problems in coronal physics is that the magnetic field and flow velocity cannot actually be measured at the coronal boundary, but rather are measured at the photosphere, below the intervening chromosphere; the nature of the energy transmission through the chromosphere remains a topic of important study. This complication notwithstanding, estimates of the Poynting flux through the photosphere are difficult since measuring the quantities in Eq. (4) with sufficient accuracy, spatial resolution and cadence remains challenging. Various techniques do exist that involve reconstructing either the horizontal flows or electric field from magnetogram and doppler data, that are now becoming more robust (Welsch et al. 2007; Kazachenko et al. 2014).

One approach to the problem of estimating the Poynting flux has been to consider ‘typical’ motions, of either ‘twisting’ or ‘braiding’ type (Berger 1993; Zirker and Cleveland 1993). From such studies it is argued that certain types of motion are most efficient at injecting energy into the corona. On the other hand estimates exist for the Poynting flux based directly on observations. Yeates et al. (2014) and Welsch (2015) used Fourier Local Correlation Tracking (FLCT) to obtain photospheric velocities in plage regions, and found net-upward Poynting fluxes of around \(5\times 10^4\mathrm {\ Wm}^{-2}\), more than sufficient to balance the energy losses as per the discussion above (see Fig. 2a–c). Welsch (2015) also found that the per-pixel Poynting flux increased with the measured local magnetic field strength. Yeates et al. (2014) have compared their estimates with those from convection simulations, and the fluxes implied by the simulations are found to be slightly higher, but comparable. Shelyag et al. (2011, 2012) also studied vertical Poynting fluxes in simulations of the convection zone (see Fig. 2d–f), and identified horizontal motions of plasma in strong intergranular magnetic flux concentrations as providing the dominant contribution.

Fig. 2
figure 2

ac Respectively, the large scale magnetogram, a zoom in to the region of interest (with units of G), and the inferred Poynting flux density (in erg cm\(^{-2}\) s\(^{-1}\)). df Respectively, the modulus of the horizontal velocity, the vertical magnetic field, and the Poynting flux density at the photospheric level in magneto-convection simulations. Images reproduced with permission ac from Yeates et al. (2014), copyright by ESO; and df from Shelyag et al. (2011), copyright by the authors

The question whether the stored coronal magnetic energy should in general be expected to increase in response to random motions (corresponding to net positive Poynting flux) can also be approached from a topological viewpoint. This stems from the fact that the tangling of magnetic field lines puts a lower bound on the magnetic energy (Moffatt 1985). Now, by analogy with the statistics of random braids, it can be argued that if one applies a random set of motions to the ends of the field lines, it is more likely that the complexity of the tangling increases (Nechaev 1999). Intuitively, this is because there are are a greater number of paths that lead to further entanglement than to disentanglement. In other words the probability of increasing the tangling—and therefore the stored energy—is greater than that of decreasing the tangling. Since the field lines tend to become progressively more tangled, then it is to be expected that the free energy of the coronal field should steadily increase until some threshold for energy release is reached (see later).

All of the above suggests that the Poynting flux through the photosphere is sufficient to provide the energy to balance the radiative losses of the corona. However, the mechanism by which the energy conversion takes place, and where this occurs in the atmosphere, remain open questions. One question that naturally arises is the extent to which the photospheric motions actually braid the coronal field lines. Candelaresi et al. (2018) have made a first attempt at addressing this issue (extending earlier work by Yeates et al. 2012). They used flows extracted using FLCT of plage-region Hinode/SOT observations (see Yeates et al. 2012, and references therein), as well as flows taken from the top layer of magneto-convection simulations. They concluded that the flows from the magneto-convection simulations would induce substantial braiding on a timescale of a few minutes, while for the flows extracted from the magnetogram series the timescale was around 3 h. It should be noted, however, that the FLCT may under-estimate the helicity injection by the flow by as much as an order of magnitude (Welsch et al. 2007). Thus, more studies of this type with refined data and flow extraction methods are required to assess the rate at which field line footpoints are tangled by the photospheric flows.

4 Existence and structure of magnetohydrostatic equilibria: theory

4.1 Force-free equilibria: general properties

The original hypothesis by Parker (Sect. 2) was based on an argument against the existence of equilibria for fields of arbitrary topology. In the following we outline what is known on the existence and stability of force-free magnetic fields—often considered due to the small plasma-\(\beta \) in the corona. We begin with some simple fundamental arguments about the geometry of solutions to the force-free field equation.

First let us remark that the problem of finding solutions to the equation

$$\begin{aligned} {\mathbf{B}}\times \nabla \times {\mathbf{B}}=0 \end{aligned}$$
(5)

has a long history. It first appeared in fluid dynamics as a condition for a class of stationary solution to the Euler equation in a paper by Beltrami (1889). Correspondingly, the solutions are also known as Beltrami fields and the particular case of \(\nabla \times {\mathbf{B}}=\alpha {\mathbf{B}}\) with \(\alpha \) constant is sometimes also called a Trakalian field or linear force-free field: these are the Eigenfields of the curl operator. One should also note that the equation

$$\begin{aligned} \nabla \times {\mathbf{B}}(x) =\alpha (x) {\mathbf{B}}(x) \end{aligned}$$
(6)

is often used as an equivalent version to (5) although this is strictly speaking not correct, in that \({\mathbf{B}}\times \nabla \times {\mathbf{B}}=0\) allows solutions where, for instance, \(\nabla \times {\mathbf{B}}\) has a non-vanishing value (or even a singularity) at points where \(\mathbf{B}=0\) and hence \(\alpha (x)\) does not exist. The original current sheet solution by Syrovatskii (1971) is of this type. However, as we are interested foremost in the case where \(\mathbf{B}\) is non-vanishing in our domain we can work with (6) without restriction of generality. This equation looks on first sight deceptively simple, but since \(\alpha (x)\) is a function of the solution it is a tricky, non-linear partial differential equation. This can be seen by taking the divergence of (6) to see that the parameter \(\alpha \) has to be constant along field lines:

$$\begin{aligned} \nabla \times {\mathbf{B}}=\alpha {\mathbf{B}},\quad \Rightarrow {\mathbf{B}}\cdot \nabla \alpha =0. \end{aligned}$$
(7)

This means that, whenever \(\nabla \alpha \ne 0\), field lines lie on flux surfaces of \(\alpha =\rm{const.}\). Although the function \(\alpha \), and hence the geometry of these flux surfaces, are part of the solution of the force-free condition, their mere existence can be a severe restriction on the topology of the field. In particular, this is the case when we have recurrent field lines as for instance for a magnetic field in a torus. A simple topological argument can be used to show that the only closed \(\alpha =const.\) surfaces for fields with \(\mathbf{B}\ne 0\) are tori (Arnold 1986). This argument is based on the fact that the torus is the only closed compact surface in \(\mathbb {R}^3\) with Euler characteristic 0. The argument extends to the case of magnetohydrostatic equilibria, where field lines lie on surfaces of constant pressure, which are again constrained to be tori, for \(\nabla p\ne 0\). The typical solution of a force-free field in a toroidal domain (\(\mathbf{B}\ne 0\)) hence has to consist of regions where the field lines lie on surfaces which are topologically equivalent to nested tori. These regions can be separated by surfaces which are not tori, for which \(\nabla \alpha = 0\). In this geometry any field line is bound to a level surface of \(\alpha \). This excludes regions of volume filling ergodic field lines, which are a common phenomenon in fields which have no particular symmetry. Hence, regions with ergodic field lines can be force-free only with \(\alpha =\rm{const.}\) (Moffatt 1985; Vainshtein 1992). While there are indeed linear force-free fields which show such ergodic regions (for example the so-called ABC fields), there exist for instance no such fields in the torus with \(\mathbf{B}\ne 0\) and zero total helicity (Pontin et al. 2016).

The above considerations provide us with examples of how to prove non-existence of force-free fields for large classes of fields in a torus, or equivalently a cylinder with periodic boundaries in the z-direction (see also Tsinganos et al. 1984). However, this is not the setting Parker addressed. He was considering straightened out coronal loops, that is a magnetic field \(\mathbf{B}\ne 0\) with field lines connecting two opposite boundaries of the domain. In this case the surfaces \(\alpha = \rm{const.}\) end on the boundaries which is a less restrictive condition than periodicity, and hence the topological arguments for non-existence given above do not directly apply.

4.2 Linear perturbative analyses

Parker (1972) originally presented his hypothesis by considering small perturbations of a magnetic field about an equilibrium. Specifically, he took \({\mathbf{B}}=B_0{\mathbf{e}}_z+{\mathbf{b}}(x,y,z)\), with \(B_0\) constant, and made a series expansion of \({\mathbf{b}}\) and the pressure p in terms of a small parameter, \(\epsilon \), which is of order \(|{\mathbf{b}}|/B_0\):

$$\begin{aligned} p=\sum _{n=0}^\infty \epsilon ^n p_n, \quad {\mathbf{b}}=\sum _{{n=1}}^\infty \epsilon ^n {\mathbf{b}}_n , \end{aligned}$$
(8)

where \(\nabla \cdot {\mathbf{b}}_n=0\) and \(|{\mathbf{b}}_n|\sim B_0\) for all n. That \(\epsilon \) is small is a result of the assumption that the footpoint displacements (of order \(l_\perp \)) are small compared to the length of the loop (length of the domain in the z-direction, \(l_z\)). Inserting into the equilibrium equation

$$\begin{aligned} \frac{1}{4\pi }(\nabla \times {\mathbf{B}})\times {\mathbf{B}}=\nabla p \end{aligned}$$
(9)

and collecting terms in powers of \(\epsilon \) yields

$$\begin{aligned}&\nabla p_0+\left[ \nabla \left( p_1+\frac{B_0b_{1z}}{4\pi }\right) -\frac{B_0}{4\pi }\frac{{\partial {\mathbf{b}}_1}}{\partial z}\right] \epsilon \nonumber \\&\quad +\left[ \nabla \left( p_2+\frac{B_0b_{2z}}{4\pi }+\frac{b_1^2}{8\pi }\right) -\frac{B_0}{4\pi }\frac{{\partial {\mathbf{b}}_2}}{\partial z}-\frac{{\mathbf{b}}_1\cdot \nabla {\mathbf{b}}_1}{4\pi }\right] \epsilon ^2+\cdots =0. \end{aligned}$$
(10)

(Note that in this section we use Parker’s original notation in Gaussian (cgs) units to facilitate comparison with the original paper, as opposed to the SI (mks) units used above—hence the difference between Eqs. (1) and (9).) In Parker’s analysis, coefficients of powers of \(\epsilon \) are set to zero, with the first order term giving

$$\begin{aligned} \nabla \left( p_1+\frac{B_0b_{1z}}{4\pi }\right) =\frac{B_0}{4\pi }\frac{{\partial {\mathbf{b}}_1}}{\partial z}. \end{aligned}$$
(11)

Together with \(\nabla \cdot {\mathbf{b}}_1=0\), Eq. (11) implies \(\nabla ^2\left( p_1+B_0b_{1z}/4\pi \right) =0\). Now, since the length-scale of variations of \(\left( p_1+B_0b_{1z}/4\pi \right) \) on the boundaries is much smaller than the loop length, it is argued that the only admissible solution is that \(\left( p_1+B_0b_{1z}/4\pi \right) \) is constant throughout the domain (except, perhaps, in a boundary layer of thickness \(l_\perp /l_z\)). Inserting back into Eq. (11) shows that \(\partial \mathbf{b_1}/\partial z=0\). Parker then goes on to show that \(\partial \mathbf{b}_{n}/\partial z=0\) for all n. On the basis of this analysis, Parker (1972) concluded that the existence of a smooth equilibrium requires \(\partial \mathbf{b}/\partial z=0\), where \(\mathbf{b}\) is the magnetic field perturbation. The natural conclusion is then that any boundary displacement that is different on the two end-plates, which violates this condition, does not admit a smooth equilibrium. Similar conclusions were drawn by Tsinganos (1982).

As first realised by van Ballegooijen (1985), there is an error in the perturbation expansion outlined above, and used by Parker (1972) to obtain the condition \(\partial \mathbf{b}/\partial z=0\). In particular, van Ballegooijen (1985) argued that the term on the right-hand side of Eq. (11) is actually of order \(\epsilon \). It should therefore show up in the second order term in the expansion (10) (corresponding terms move from 2nd to 3rd order, etc.), meaning that Eq. (11) for the first-order term is modified to

$$\begin{aligned} \nabla \left( p_1+\frac{B_0b_{1z}}{4\pi }\right) =0. \end{aligned}$$
(12)

While we still arrive at the conclusion that \(\left( p_1+B_0b_{1z}/4\pi \right) \) must be constant, it no longer follows that \(\partial \mathbf{b}_{1}/\partial z= 0\). Thus the set of equilibrium solutions is generalised to include fields that depend on z. Indeed, the analysis by van Ballegooijen (1985) indicates that the components of the magnetic field perturbation \(\mathbf{b}\) can be expressed in a straightforward manner as a function of the boundary displacements. As such—it is argued—so long as the boundary motions are smooth, a corresponding smooth equilibrium can be found.

A similar analysis to that of van Ballegooijen (1985) was made by Zweibel and Li (1987), who considered small perturbations to both a uniform field and a weakly sheared field (their analysis builds on that of Sakurai and Levine 1981). They use a Lagrangian approach, writing the magnetic field in terms of the fluid displacement. This allows the frozen-in condition to be built into the equilibrium equation—overcoming the major challenge of specifying the topology of the perturbed field. The method contrasts with the Eulerian approach of van Ballegooijen (1985), however their conclusion is the same; that the magnetic field resulting from the class of small perturbations considered is explicitly (and uniquely) determined as a function of this displacement, and is therefore smooth (the boundary displacement is assumed to be smooth on physical grounds). Solutions are found to be approximately independent of z in the majority of the domain, with the notable exception of thin layers close to the line-tied boundaries, in which much of the stress in the field is taken up (see also Sakurai and Levine 1981). This conclusion was supported and expanded upon by Craig and Sneyd (2005), who developed a theory that includes arbitrary small footpoint displacements using a Fourier expansion. It should be noted that the interpretation of these solutions has been called into question by Low (2010a), though see also Craig (2010). Further examples of smooth equilibria deriving from small perturbations are presented by Rosner and Knobloch (1982) and Arendt and Schindler (1988).

One study in which a perturbative analysis does predict the formation of singular current sheets is the study of Bobrova and Syrovatskii (1979). They considered perturbations of the one-dimensional force-free magnetic field \({\mathbf{B}}=\left( \cos (\alpha z),\sin (\alpha z), 0\right) \), and demonstrated that the current becomes singular on certain surfaces in response to rather general boundary displacements in the absence of plasma pressure. However, it should be noted that the geometry they considered is very different to the one addressed by Parker’s hypothesis—in particular the domain is infinite along the direction of the magnetic field, with the perfectly conducting z-boundaries being tangent to the magnetic field lines.

In summary, various authors have developed solutions demonstrating that the original argument put forward by Parker (1972) was flawed. These solutions demonstrate the existence of smooth equilibria for broad classes of small perturbations to initial force-free fields, in the line-tied geometry envisaged by Parker (see Sect. 2). However, there remains no universal consensus that this rules out current sheet formation for generic small perturbations (e.g., Parker 2000; Low 2010b). It is also worth noting that while many of the linear analyses mentioned above demonstrate the existence of neighbouring smooth equilibria, they fail to address either the stability of these equilibria or their dynamic accessibility from the initial (unperturbed) equilibrium. As such they do not preclude the formation of current sheets during a dynamic evolution when the associated perturbations are applied. Assessing these issues requires a different approach to those described above, as discussed in the following sections.

Moreover, referring to the studies of van Ballegooijen (1985), Zweibel and Li (1987) and others, Parker (1994) notes that “the known continuous equilibrium solutions involve either only weak deformation of the field from a uniform state ...or a symmetry, degeneracy, or invariance (ignorable coordinate) of some form”. This is suggestive of a reformulation of the problem; that tangential discontinuities are expected to form when the footpoint displacements are of sufficient amplitude and complexity. As noted by Craig and Sneyd (2005), “this caveat means that the magnetostatic theorem no longer admits a precise mathematical statement”, making it difficult to unequivocally prove/disprove. In the following sections we turn to the case of finite-amplitude footpoint displacements.

4.3 Arguments against the existence of smooth equilibria

4.3.1 Parker’s optical analogy

A persuasive intuitive argument for the formation of tangential discontinuities was put forward by Parker, by means of the optical analogy (Parker 1989a, b, 1991, 1994). This optical analogy is based on the notion that a field line in a potential field takes the same path \({\mathbf{x}}(s)\) as an optical ray in a medium with index of refraction equivalent to \(|{\mathbf{B}}|\):

$$\begin{aligned} \frac{dx_i}{ds}=\frac{1}{n}\frac{\partial \psi }{\partial x_i} , \end{aligned}$$
(13)

where in the magnetic field case \({\mathbf{B}}=\nabla \psi \), \(n=|{\mathbf{B}}|\), in the optical case \(\psi \) is the phase of the electromagnetic wave and n is the index of refraction, and in both cases s is the arc-length along the path. This optical analogy is outlined below for the case of a force-free field. In such a force-free field

$$\begin{aligned} \nabla \times {\mathbf{B}}=\alpha {\mathbf{B}},\quad {\mathbf{B}}\cdot \nabla \alpha =0 , \end{aligned}$$
(14)

and flux surfaces of the magnetic field and the current density are coincident. Now, consider a flux surface of \(\nabla \times {\mathbf{B}}\). On such a flux surface \(\nabla \times {\mathbf{B}}\cdot {\mathbf{n}}=0\), where \({\mathbf{n}}\) is the normal to the flux surface, and therefore by Stokes’ theorem \(\oint {\mathbf{B}}\cdot d\mathbf{l}=0\) for any closed curve lying in the flux surface. This in turn means that in this surface \({\mathbf{B}}\) can be represented by a (two-dimensional) potential. Note that Parker supposes that the optical analogy applies only locally to flux surfaces—the global existence of such surfaces is far from assured, particularly in periodic geometries. Thus the path of the field lines in the flux surface is determined by Eq. (13).

Now, in optics, Fermat’s principle states that the optical length of a path followed by a light ray between two points is an extremum, where the optical length is the integral of the refractive index along the path. By analogy, it is argued, magnetic field lines in the flux surface should take the path that minimises

$$\begin{aligned} \int |{\mathbf{B}}(s)|\, ds. \end{aligned}$$
(15)

In order to do this, field lines should tend to avoid regions of high magnetic field strength. If the field strength is sufficiently high in some location on the flux surface, it is expected that field lines can be completely excluded from this region, effectively forming a ‘hole’ in the flux surface. This hole allows two bundles of magnetic flux that are elsewhere separated by the flux surface to come into contact (across the hole), and presumably the directions of the field lines in these bundles may not match, leading to a tangential discontinuity—see Fig. 3.

Fig. 3
figure 3

Illustration of the “optical analogy”, after Parker (1994). a, b Field lines are refracted away from a local maximum in the magnetic field strength (ellipse centred on the y-axis). In b the maximum is sufficiently strong to form a hole or exclusion zone for field lines in the surface. c ‘Hole’ formed in a stack of flux surfaces, through which field lines on either side of the stack may meet to form a tangential discontinuity

Parker has developed various example of holes or exclusion zones in flux surfaces based on the use of Fermat’s principle, for imposed local maxima of \(|{\mathbf{B}}|\) with various different geometries (Parker 1989b, 1994). However, while Parker argues persuasively that tangential discontinuities should form on the basis of these considerations, the fact that the field strength is specified a priori in the equilibrium with a particular distribution makes any direct verification of the proposition difficult to test. One concrete example is presented in Parker (1990); however, it does not conform to the geometry of the Parker problem, since lateral boundaries (perpendicular to the line-tied z-boundaries and parallel to \({\mathbf{B}}\)) are assumed to provide an external pressure force (i.e., these lateral boundaries are “squeezed in” in some range of z between the line-tied boundaries in order to provide a magnetic field enhancement).

4.3.2 Work of Low et al. on “topologically untwisted” fields

In a series of papers, Low (2006), Low and Flyer (2007), Janse and Low (2009) and Janse et al. (2010) considered what they describe as “topologically untwisted” fields (the former two papers dealing with the field outside a sphere and the latter two the field in a cylinder as per Parker’s hypothesis). These are fields for which the net circulation on any flux tube cross-section is zero, and thus they are either potential fields or fields in which the current is confined to sheets of zero thickness (tangential discontinuities). Janse and Low (2009) focus on the field inside a cylinder with perfectly-conducting boundaries and zero flux through the curved side boundary. They begin with a potential field in the volume and perturb the system by shrinking the length of the cylinder while maintaining the boundary distribution of \(\mathbf{B}\cdot \mathbf{n}\). They then crucially argue that since the field is topologically untwisted, and since no twist has been applied via the boundary perturbation, then the field must remain topologically untwisted. They proceed to calculate the new potential field in the compressed volume for specific examples and demonstrate that its topology is not the same as that of the initial field. Thus, they argue, if a relaxation ensues in which the topology is preserved, current sheets must form. However, Aly and Amari (2010) have questioned the assumption in this argument that the field must remain topologically untwisted. Indeed two counter-examples have been presented that suggest that the field can in general relax to a smooth equilibrium that is not topologically untwisted, but rather contains a distributed current—the first of these is based on a linear expansion scheme in 2D by Huang et al. (2009) (which, as pointed out by Janse et al. (2010) does not strictly contradict the claim of Janse and Low (2009) for fully 3D fields) and the second for the fully 3D field of Janse and Low (2009) by Pontin and Huang (2012) using a magneto-frictional simulation (see Sect. 5.2).

4.3.3 Uniqueness of equilibria in reduced MHD

Another argument in support of Parker’s hypothesis (general form, Definition 1) was put forward by Ng and Bhattacharjee (1998) who made use of the Reduced MHD (RMHD) approximation (see Sect. 6.3.1). They exploited the parallel between the RMHD equations and the 2D Euler equation in hydrodynamics to prove the theorem: “For any given footpoint mapping connected smoothly with the identity mapping ...there is at most one smooth equilibrium”. They point out the important implications for a smooth but unstable equilibrium: specifically that following loss of stability the system must relax towards a non-smooth state, the unique smooth state being the unstable one. They also demonstrate that this relaxation towards a non-smooth state—in which the field is not it equilibrium in the vicinity of the current sheets (regions of non-smoothness)—nevertheless can in principle minimise the magnetic energy. It is worth noting that the proof of the above theorem relies on the approximations of RMHD (see Sect. 6.3.1) and triply-periodic boundary conditions, and therefore the implications for the full MHD system with line-tied boundaries are unclear. Indeed, often equilibria obtained in the full MHD system exhibit boundary layers that are incompatible with the RMHD assumptions (see Sect. 4.2). While the uniqueness of smooth equilibria in a system governed by the full MHD equations has not been demonstrated in general, what has been shown for certain boundary conditions by Aly (2005) is that there is no smooth, force-free field topologically equivalent to the uniform field, other than the uniform field itself.

4.4 Arguments in favour of smooth equilibria

The studies described above in Sect. 4.3 address existence of equilibria in certain classes of field, but a general statement regarding the conditions for the absence of an equilibrium remains elusive. In this section we turn to arguments against the general form of Parker’s hypothesis.

4.4.1 Rigorous, general results on the existence of force-free fields

There are very few general results concerning the existence of force-free fields. The strongest statement so far made is due to Bineau (1972). In that work, an iterative scheme is constructed for obtaining a magnetic field \({\mathbf{B}}\) in a domain D satisfying

$$\begin{aligned}&\nabla \cdot {\mathbf{B}}=0,\quad \nabla \times {\mathbf{B}}=\beta \sigma {\mathbf{B}}, \quad {\mathbf{B}}\cdot \nabla \sigma =0, \quad \mathrm{in} \quad D({\mathbf{B}},S_1),\\&\nabla \cdot {\mathbf{B}}=0,\quad \nabla \times {\mathbf{B}}=0 , \quad \mathrm{in} \quad D-D({\mathbf{B}},S_1),\\&{\mathbf{n}}\cdot {\mathbf{B}}=\tilde{B}_n, \quad \mathrm{on}\quad S,\\&\sigma =\tilde{\sigma }, \quad \mathrm{on}\quad S_1, \end{aligned}$$

where S is a surface with normal vector \({\mathbf{n}}\) that is transverse to the field lines, \(S_1\) is a connected part of S on which the normal component of \({\mathbf{B}}\) is strictly greater than zero, \(\beta \) is a dimensionless constant, \(\tilde{B}_n\) and \(\tilde{\sigma }\) are given functions, and \(D({\mathbf{B}},S_1)\) is the sub-volume of D that is permeated by the magnetic flux passing through \(S_1\). Following construction of the iterative scheme, the question of existence of a force-free field as per the above equations becomes a question of whether the iteration converges. The iterative scheme involves adding successive perturbations to a potential field \({\mathbf{B}}^0\), and the convergence criterion then turns out to involve parameters including the strength of this field, the modulus of \(\sigma \) and its derivatives, and the maximum length of field lines within \(D({\mathbf{B}},S_1)\). The upshot is that the iteration converges for sufficiently small \(\beta \). The conclusion is that for any potential field \({\mathbf{B}}^0\) in D and any given \(\tilde{\sigma }\) on \(S_1\) “there exists a family of force-free fields in D that depends analytically on the parameter \(\beta \) in the neighborhood of \(\beta =0\)” (Bineau 1972). In one sense this is a powerful result: it means that force-free fields of arbitrary topology satisfying \(\nabla \times {\mathbf{B}}=\alpha {\mathbf{B}}\) exist for sufficiently small \(\alpha \) in a geometry consistent with the Parker problem. Note that since \(\alpha \) is not dimensionless we can always scale the domain (in the z-direction in Parker’s geometry) to make, for any given fieldline mapping, \(\alpha \) as small as we like. However, what we still do not know is how close to the potential field we must be. Or in other words, how large is the range of \(\alpha \), in the vicinity of zero, within which force-free fields are guaranteed to exist, and how does this depend on the field topology or other factors?

There do exist some general existence and stability results for the restricted case of linear (or constant-\(\alpha \)) force-free fields in a half-space (Aly 1992). However, such linear force-free fields are not relevant to the complex field of the solar corona. Aly (1990) derives sufficient conditions for linear ideal-MHD stability of general force-free fields, again in a half-space—loosely speaking that the product of the force-free parameter \(\alpha \) and the characteristic length scale of the field should be order-1 or smaller. Examples of exact force-free equilibria with complex 3D topology have been presented by Bogoyavlenskij (2000a, 2000b). However, as argued by Parker (2000), this does not itself imply that such smooth equilibria exist in general for arbitrary field topologies.

One of the principal reasons why a general proof regarding the existence or otherwise of equilibria in the Parker problem remains so challenging arises from the need to precisely specify and maintain the magnetic field topology while solving the force-free equations (see, e.g., Low 2010b). Indeed, Parker (1986) argued that the requirement to simultaneously satisfy \({\mathbf{J}}\times {\mathbf{B}}=0\) and the requirements for a given topology (field line connectivity) in general leaves the mathematical problem over-specified—or in other words that the two requirements are in general mutually incompatible. However, this claim was refuted by Antiochos (1987), who considered the relative winding of field lines about one another in a half-space above a line-tied photosphere, using an Euler-potential representation of the coronal field: \({\mathbf{B}}=\nabla \alpha \times \nabla \beta \) (note that while this geometry is closer to that in the corona, it is not the same as the standard Parker problem as defined in Sect. 2). He showed that the net winding of all field lines about a given coronal field line is entirely specified by the boundary conditions (boundary distributions of \(\alpha \) and \(\beta \)). Thus, he argued that the topological constraints on the coronal field are included in the boundary conditions, and do not provide an extra constraint on the solution of the force-free equations, meaning that the problem of finding a solution is not a priori over-specified (as had been argued by Parker 1986). The existence of a solution to those equations, however, is not directly addressed. In the following section, we describe another argument against the formation of tangential discontinuities in the volume in response to smooth boundary motions.

4.4.2 Permissibility of tangential discontinuities

Several authors have asserted that the formation of tangential discontinuities as predicted in Parker’s hypothesis is inadmissible in the line-tied geometry in question. The argument, which we present in the following, appears in slightly different forms in the papers by van Ballegooijen (1985, 1988a); Longcope and Strauss (1994), and Cowley et al. (1997). We follow most closely the latter. We proceed by assuming that a current sheet forms, and showing that this leads to a contradiction. This current sheet is represented in Fig. 4. Consider the two magnetic field lines, along the paths labelled \(\mathcal {C}_1\) and \(\mathcal {C}_2\), that are immediately on either side of this sheet (for concreteness, field line ‘1’ is in front of the sheet from the perspective shown, field line ‘2’ is behind). These field lines are chosen such that they are infinitesimally separated on the \(z=0\) line-tied boundary. Now, since the boundary flows are assumed to be smooth, these two field line footpoints must also have been infinitesimally separated at \(t=0\) (i.e., in the initial equilibrium). But since the initial equilibrium has a simple topology the footpoints of these field lines at \(z=L\) are infinitesimally separated at \(t=0\). Finally, invoking again the smooth nature of the boundary flows (this time at \(z=L\)), the field line footpoints remain infinitesimally separated at \(z=L\), like at \(z=0\), for all times.

Fig. 4
figure 4

Schematic of a current sheet (light shaded) in a magnetic field between line-tied perfectly-conducting boundary plates at \(z=0,L\) (dark shaded). On the front side of the current sheet the magnetic field is \({\mathbf{B}}_1\) and a representative magnetic field line follows the path \(\mathcal {C}_1\). Behind the current sheet, \({\mathbf{B}}={\mathbf{B}}_2\), and a field line following the path \(\mathcal {C}_2\) is displayed

The next step of the argument relies on the result, presented by Cowley et al. (1997), that in a force-free field the current flows entirely along the current sheet and never leaves it, and therefore the current path on a finite, simply-connected sheet must close. Thus,

$$\begin{aligned} \oint {\mathbf{B}}\cdot d\mathbf{l}=0 \end{aligned}$$
(16)

for any closed loop on either surface of the current sheet—either in front or behind in the perspective of Fig. 4 (true so long as the loop does not cross over the sheet). To construct such a closed loop, define \(\mathcal {C}_2^\prime \) to be a path that runs along \(\mathcal {C}_2\) but on the same side of the sheet as \(\mathcal {C}_1\). Similarly, define \(\mathcal {C}_1^\prime \) to run along \(\mathcal {C}_1\), but on the same side of the sheet as \(\mathcal {C}_2\). Now consider the loop integral of \({\mathbf{B}}\) along the path comprised of \(\mathcal {C}_1\) and \(\mathcal {C}_2^\prime \). This path remains everywhere in front of the sheet, and is parallel to \({\mathbf{B}}\) along \(\mathcal {C}_1\), but is not parallel to \({\mathbf{B}}\) along \(\mathcal {C}_2^\prime \). Therefore

$$\begin{aligned} \oint {\mathbf{B}}\cdot d\mathbf{l} = \int _{\mathcal {C}_1}{\mathbf{B}}\cdot d\mathbf{l} - \int _{\mathcal {C}_2^\prime }{\mathbf{B}}\cdot d\mathbf{l} = \int _{\mathcal {C}_1}B_1(l) dl - \int _{\mathcal {C}_2^\prime }{\mathbf{B}}_1(l)\cdot \hat{{\mathbf{B}}}_2(l) dl =0 , \end{aligned}$$
(17)

the final equality holding due to condition (16), and the ‘hat’ denoting a unit vector. Similarly, considering a loop behind the current sheet,

$$\begin{aligned} \oint {\mathbf{B}}\cdot d\mathbf{l} = \int _{\mathcal {C}_2}{\mathbf{B}}\cdot d\mathbf{l} - \int _{\mathcal {C}_1^\prime }{\mathbf{B}}\cdot d\mathbf{l} = \int _{\mathcal {C}_2}B_2(l) dl - \int _{\mathcal {C}_1^\prime }{\mathbf{B}}_2(l)\cdot \hat{{\mathbf{B}}}_1(l) dl =0. \end{aligned}$$
(18)

Now, across a tangential discontinuity \(|{\mathbf{B}}|\) does not change (Landau and Lifshitz 1960), so that \(|{\mathbf{B}}_1(l)|=|{\mathbf{B}}_2(l)|=B(l)\). Then, adding Eqs. (17) and (18) we have

$$\begin{aligned} \oint _{\mathcal {C}_1+\mathcal {C}_2} B(l)\left( 1-\cos \theta (l)\right) \, dl=0 , \end{aligned}$$
(19)

where \(\theta (l)\) is the angle between \({\mathbf{B}}_1\) and \({\mathbf{B}}_2\). The integrand is everywhere greater than zero, so the only way this equation can hold is if \(\cos \theta =1\) all along the field line, i.e., \(\theta =0\). This contradicts the presence of a tangential discontinuity, and we conclude that no such current sheet can exist.

The above argument appears to rule out the formation of tangential discontinuities in response to smooth boundary driving. While they are based on force-free fields, the argument was extended to the case of finite pressure by Field (1990). However, Ng and Bhattacharjee (1998) have pointed out that this argument applies only to current sheets with simple structure. They present an example of a current sheet that exhibits a ‘branching’ geometry, for which the above argument cannot be applied (in a cross-section through the sheet at constant z it exhibits a ‘Y’-type geometry). Cowley et al. (1997) admit that the above argument does not hold for current sheets of arbitrary topology, but attest that current sheets of “arbitrarily complicated geometries” are “pathological”, and should not be expected in practice. This claim is refuted by Ng and Bhattacharjee (1998), who argue that confining considerations to only current sheets of simple topology is “over-restrictive”. Chapter 8 of Parker (1994) discusses the formation of current sheets with non-trivial topologies by means of the optical analogy (see Sect. 4.3.1).

In conclusion, it appears that the simplest, topologically-planar, tangential discontinuities are prohibited from forming in response to smooth boundary motions in an ideal plasma. However, current sheets with more complicated topologies—e.g., branching—may perhaps arise. It is worth noting that the above proof assumes the absence of certain structures in the equilibrium magnetic field. In particular, we assumed that field lines whose footpoints are infinitesimally separated at one end must have footpoints infinitesimally separated at the other end as well. This implies no discontinuities are present in the field line mapping, i.e., there are no null points or separatrix surfaces in the volume. Similarly, strong gradients could be present in the mapping, indicating the presence of quasi-separatrix layers. These features are explicitly excluded in Parker’s original hypothesis, though we will consider their importance further in Sect. 7.

4.4.3 Ideal instabilities in the presence of line-tying

One way to form current sheets in a 3D geometry is through an ideal instability. Two prime candidates are the coalescence instability (with a guide field) of flux tubes (Finn and Kaw 1977) and the ideal kink instability (Kruskal et al. 1958; Hood and Priest 1979). In the absence of line-tying in both cases (i.e., when the system is periodic or invariant in the ‘guide field’ direction), singular current sheets form; at X-points of the in-plane field for the coalescence instability and on resonant surfaces for the ideal kink. Reasoning along the lines described in the previous section, Longcope and Strauss (1994) considered the coalescence instability, arguing that the introduction of line-tying in the guide field direction (thus placing the system within the Parker problem geometry) inhibits the formation of singular current sheets. They showed that the singular current sheets are replaced by intense but finite current concentrations. These tend to form with a perpendicular length-scale set by the field line mapping—see discussion and references in Craig and Pontin (2014).

Similar conclusions have been drawn regarding the effect of line-tying on the kink instability, at least in the linear phase. Specifically, in the line-tied system a current layer still forms at the radius that corresponds to the resonant surface in the infinite/periodic system. However, the thickness of the current layer scales inversely with the distance between the line-tied boundaries, and the growth rate of the instability is reduced—see Huang et al. (2010) and references therein. This is shown to be because—in the absence of a resonant surface—while the internal layer associated with the fastest-growing eigenfunction persists, it is smoother than its periodic counterpart. As pointed out by Huang et al. (2010), the issue of current sheet formation in the non-linear phase of the instability is not fully resolved.

5 Existence and structure of magnetohydrostatic equilibria: computational approaches

5.1 Ideal relaxation

In the foregoing sections, mathematical results on the existence and structure of magnetostatic equilibria have been discussed. In this section we turn to a computational approach for assessing this problem. In particular, the methods discussed in this section ask the following question: For a given magnetic topology, does a smooth equilibrium exist? This is approached by means of an ideal “relaxation” in which at \(t=0\) a non-equilibrium magnetic field with a chosen topology is taken, and the forces or free energy are then minimised in search of a corresponding equilibrium, while preserving this topology. These forces may be minimised by following various different evolutions. Considering the problem as a ‘thought experiment’ the details of this evolution are unimportant: the question at hand is the existence and structure (e.g., smoothness) of the final equilibrium for the given topology—and perhaps also its stability properties.

One such approach is described in Chodura and Schlüter (1981) and Moffatt (1985). The method is based on an ideal evolution of the magnetic field

$$\begin{aligned} \frac{\partial {\mathbf{B}}}{\partial t} - \nabla \times \left( {\mathbf{v}}\times {\mathbf{B}}\right) =0 , \end{aligned}$$
(20)

where the velocity is obtained from

$$\begin{aligned} {\mathbf{v}}= \gamma (\nabla \times {\mathbf{B}})\times {\mathbf{B}}. \end{aligned}$$
(21)

Here \(\gamma \) is a scaling factor that can be absorbed into the artificial time of the evolution. Since the change of magnetic energy W in the volume V is given by

$$\begin{aligned} \frac{dW}{dt}=-\int _V {\mathbf{F}}\cdot {\mathbf{v}}\, d^3 x = -\gamma \int _V \left( (\nabla \times {\mathbf{B}})\times {\mathbf{B}}\right) ^2 \, d^3 x , \end{aligned}$$
(22)

it is obvious that the magnetic energy decreases monotonically and this process comes to a halt only when \({\mathbf{v}}= \gamma (\nabla \times {\mathbf{B}})\times {\mathbf{B}}=0\), that is when a force-free field has been reached. Thus, the relaxation should produce a force-free field if the limit of \({\mathbf{B}}\) for \(t \rightarrow \infty \) exists. However, no proof of convergence exists, and the counterexample from above of an ergodic field with zero total helicity shows that in some cases the method can fail to converge to a smooth (i.e., sufficiently differentiable) field. This failure to converge can occur in various ways. It could be that no stationary limit is obtained or that the limit is not a smooth field. The latter, in turn, could include the possibility of current sheets. Various implementations of this method exist and are discussed in Sect. 5.2. The main problem of this method is that it requires a numerical solution. However, while there are recognised criteria for identifying a smooth solution, it is much more difficult, or practically impossible, to prove the absence of a smooth solution.

5.2 Computational approaches: ideal relaxation

5.2.1 The computational challenges

Direct computational assessment of Parker’s hypothesis brings a number of challenges. Foremost among these is the requirement to precisely maintain the magnetic topology during the simulated evolution, i.e., precisely maintain the magnetic field line mapping between the two line-tied boundaries. Typical computational approaches to solving the MHD equations employ finite differences on an Eulerian mesh. In such approaches the numerical diffusion is never identically zero, so that there is a finite ‘slippage’ in the magnetic connections between boundary points. This is exacerbated when structures form on the grid scale, as by definition they must if current sheets (or indeed thin, finite, current layers) form.

In the following sections two methods are described which seek to mitigate against these difficulties. However, in all cases the representation of current singularities remains problematic. Typically, in order to test the fidelity of a numerical simulation one would look for convergence with increasing numerical resolution. However, this convergence cannot occur in the case of a current sheet: instead the maximum current density represented on the grid will diverge. The rate of this divergence may be determined by the geometrical structure of the singularity, its strength, and the grid geometry. This is discussed further below.

5.2.2 Lagrangian magnetofrictional approaches

Under ideal MHD the vector \(\mathbf{B}/\rho \) evolves according to the equation

$$\begin{aligned} \frac{D}{Dt} \frac{\mathbf{B}}{\rho } = \left( \frac{\mathbf{B}}{\rho } \cdot \nabla \right) \mathbf{v} \end{aligned}$$
(23)

where D/Dt is the material derivative, \(\rho \) is the plasma density and \(\mathbf{v}\) the plasma velocity (this follows from the ideal induction equation (2) and mass conservation). This is identical to the evolution equation of a line element \(\mathbf{\delta x}\) in a flow (see, e.g., Moffat 1978). Thus, in order to tackle the problem of exactly preserving the magnetic topology, it is most convenient to use a Lagrangian description of the fluid, in which computational grid points are fluid elements. Now, suppose that the position vectors of these fluid elements (which are assumed to move in time) are represented by \({\mathbf{x}}({\mathbf{X}},t)\), where \({\mathbf{x}}({\mathbf{X}},0) = {\mathbf{X}}\). The modified ideal induction equation (23) implies that the magnetic fields at time \(t=0\) and at \(t>0\) are related by the pull-back under the mapping \({\mathbf{x}}({\mathbf{X}},t)\):

$$\begin{aligned} {\mathbf{x}}^{*}({\mathbf{B}}({\mathbf{x}},t),t) = {\mathbf{B}}({\mathbf{X}},0) \end{aligned}$$
(24)

(Candelaresi et al. 2014). This is a modern formulation of Alfvén’s Theorem (Alfvén 1943). Using the definition of the pull-back on Eq. (24), the magnetic field components \(B_i({\mathbf{x}},t)\) at \(t>0\) can be expressed in terms of the magnetic field components at \(t=0\) and the mesh deformation:

$$\begin{aligned} B_{i}({\mathbf{x}},t) = \frac{1}{\varDelta } \sum _{j=1}^3\frac{\partial x_{i}}{\partial X_{j}} B_{j} ({\mathbf{X}},0). \end{aligned}$$
(25)

Here \(\varDelta \) is the determinant of the Jacobian matrix \(\partial x_{i}/\partial X_{j}\) and measures the local fluid compression/expansion. Now, assuming an evolution equation for the positions of the fluid elements (grid points), these positions become the primary variables of the computation: the magnetic field is subsequently calculated via Eq. (25). This avoids the requirement for an explicit equation for the magnetic field evolution, the solution of which leads to the accumulation of ‘errors’ in the magnetic topology conservation (i.e., leads to ‘field line slippage’). The Lagrangian method thus preserves the topology of field lines exactly [through Eq. (24) or (25)] as well as the magnetic flux through each surface element of the grid. Additionally, \(\nabla \cdot {\mathbf{B}}= 0\) is automatically preserved. The aim then is to employ a suitable method to update the positions of the grid points in such a way as to redistribute the magnetic stresses and thus obtain an equilibrium (see Fig. 5).

Fig. 5
figure 5

Illustration of the Lagrangian relaxation computational approach, showing some representative grid lines (grey) and magnetic field lines both prior to relaxation (left), and in the numerically obtained equilibrium (right). The grid deforms to allow field lines to equilibrate the twist/stress along their length. Based on the simulations of Wilmot-Smith et al. (2009a)

It remains to decide upon an appropriate evolution equation for the mesh points (equivalently for the fluid velocity). One approach that has been used extensively is to employ a “magneto-frictional” evolution as discussed above to move towards an equilibrium. The first implementations of the method were made by Craig and Sneyd (1986) and Longbottom et al. (1998), while a similar approach was taken by Yang et al. (1986), who applied a magneto-frictional approach to find equilibria for magnetic fields represented in terms of Euler potentials. In the magneto-frictional model the fluid elements are moved in the direction of the force,

$$\begin{aligned} \frac{D{\mathbf{x}}}{Dt}=\gamma {\mathbf{F}}= \gamma (\nabla \times {\mathbf{B}})\times {\mathbf{B}}-\nabla p , \end{aligned}$$
(26)

where \(\gamma \) is a scaling factor that can be absorbed into the artificial time of the evolution. The integrated forces within the domain are then reduced until they fall below some chosen threshold. Craig and Sneyd (1986) argue that the stability properties of the resulting (numerical) equilibrium are not affected by the fictitious evolution.

Implementation of Eq. (26) requires calculation of \(\nabla \times {\mathbf{B}}\) on the numerical mesh, which is not rectangular after the initial time. This quantity can be expressed as a function of the initial magnetic field and grid deformation, analogous to Eq. (25) (Craig and Sneyd 1986). However, as first realised by Pontin et al. (2009) small discretisation errors can combine in such a way in this calculation to produce a current \({\mathbf{J}}=\nabla \times {\mathbf{B}}\) that is far from divergence-free, in locations where the numerical mesh becomes highly distorted. When this happens the relaxation is compromised: while the numerical value of \(\int _V{\mathbf{J}}\times {\mathbf{B}}\,d^3x\) can be minimised indefinitely, the true value of \(\int _V{\mathbf{J}}\times {\mathbf{B}}\,d^3x\) saturates at some finite level (or indeed increases). This finding motivated Candelaresi et al. (2014) to develop a new implementation of the numerical scheme, wherein \(\nabla \times {\mathbf{B}}\) is calculated using a so-called mimetic derivative that has the advantage of being (numerically) divergence-free (Hyman and Shashkov 1997). This alternative implementationFootnote 1 has been shown to substantially improve the accuracy of the force-free equilibrium produced, however is still adversely affected by high mesh distortion, losing numerical stability when grid cells become concave.

The above limitations notwithstanding, the implementations of Craig and Sneyd (1986) and Candelaresi et al. (2014) have been used to study the equilibrium properties of a number of different magnetic topologies relevant to the Parker problem. These include magnetic fields between perfectly conducting plates that are subject to coherent shear (Longbottom et al. 1998; Candelaresi et al. 2015) as well as more complex ‘braiding’ patterns (Craig and Sneyd 2005; Wilmot-Smith et al. 2009a; Pontin and Hornig 2015). These studies concluded that there is no indication of the formation of singular current layers during the relaxation, despite the range and complexity of perturbations considered. Instead, smooth current distributions are present in the equilibria obtained (the peak current converging to a finite value in convergence studies with the numerical resolution). Interestingly, in many cases the current density is seen to localise in ever thinner layers as the field perturbation increases (e.g., Pontin and Hornig 2015), as discussed in more detail in Sect. 5.3.

It is worth discussing briefly here the numerical representation of tangential discontinuities, on which these claims are based. As mentioned above, with increasing numerical resolution the current density can be expected to diverge in the presence of a tangential discontinuity. If the numerically represented current density is \(J\sim \delta B/\delta x\) then roughly speaking the jump \(\delta B\) should be constant with increasing resolution, while \(\delta x\) decreases. For a fixed Eulerian mesh one therefore expects \(J_{\max }\propto (\delta x)^{-1}\), whereas for a Lagrangian mesh a compression of mesh points towards the ‘sides’ of the current sheet can be expected, and flux conservation arguments show that a power-law dependence on resolution is expected; \(J_{\max }\propto (\delta x)^{-m}\) (Ali and Sneyd 2001). This type of peak current divergence with resolution has been reproduced in simulations of null point collapse with the Lagrangian magneto-frictional methods under discussion (Craig and Litvinenko 2005; Pontin and Craig 2005; Candelaresi et al. 2015), but is absent in configurations consistent with the Parker problem.

It should be noted that Low (2013) has raised concerns over the validity of the magneto-frictional approach, and thus the results derived from these ‘relaxation experiments’. He argued that the absence of fluid inertia leads to an incompatibility with the imposed boundary conditions. To test these concerns, Candelaresi et al. (2015) implemented various different evolution equations in their code as alternatives to Eq. (26), both for a fluid without inertia (the magneto-frictional case) and a fluid with inertia. They did not find any indication that the final state of the relaxation was dependent on the choice of evolution equation, for the range of configurations considered.

5.2.3 A “variational integrator” method

An alternative numerical Lagrangian mesh approach to seeking force-free equilibria exists. The mathematical basis of the method is similar to that described above. In particular, by working with Lagrangian labelling, the frozen-flux condition is built into the method via Eq. (25). The basis for the equations solved is the Lagrangian for ideal MHD presented by Newcomb (1962), which in Eulerian labelling is

$$\begin{aligned} L({\mathbf{v}},\rho ,p,{\mathbf{B}})=\int \left( \frac{1}{2}\rho v^2-\frac{p}{\gamma -1}-\frac{B^2}{2}\right) d^3 x , \end{aligned}$$
(27)

where \(\gamma \) is the ratio of specific heats and an adiabatic fluid is assumed such that \(p/\rho ^\gamma \) is constant. Writing this in Lagrangian labelling allows conservation laws for mass and magnetic flux to be built in, and the Euler–Lagrange equation for this Lagrangian (in the Lagrangian labelling) turns out to be the momentum equation (Newcomb 1962)

$$\begin{aligned} \rho _0 {\ddot{x}}_i - B_{0j}\frac{\partial }{\partial x_{0j}}\left( \frac{x_{ik}B_{0k}}{\varDelta }\right) + \frac{\partial \varDelta }{\partial x_{ij}} \frac{\partial }{\partial x_{0j}} \left( \frac{p_0}{\varDelta ^\gamma }+\frac{x_{kl}x_{km}B_{0l}B_{0m}}{2\varDelta ^2}\right) = 0 , \end{aligned}$$
(28)

where \({\mathbf{x}}({\mathbf{X}},0)={\mathbf{X}}\), \(\rho _0=\rho({\bf X},0)\), \(p_0=p({\bf X},0)\), \({\mathbf{B}}_0={\mathbf{B}}({\mathbf{X}},0)\), \(\varDelta =\det x_{ij}\) as above, the ‘dot’ denotes a time derivative, and summation over repeated indices is assumed. Note that in addition to Eq. (25) we have \(\rho =\rho _0/\varDelta \) (mass conservation) and \(p=p_0/\varDelta ^\gamma \). Zhou et al. (2014) describe a numerical method for solving this equation. They use discrete exterior calculus (Desbrun et al. 2005, and references therein) to first discretise the Lagrangian. They then use variational integrators for the time-discretisation on the Lagrangian directly, with the update scheme for \({\mathbf{x}}\) then determined from the discrete Euler–Lagrange equation (it is argued that this discretisation more faithfully preserves the energy than a direct discretisation of Eq. (28)). Of course, eventually during the relaxation the energy must be damped, and the authors introduce a friction term \(-\nu \rho \dot{{\mathbf{x}}}\) into the equation of motion for this purpose.

This variational integrator method has been used to demonstrate current sheet formation in the 2D configurations of coalescence instability (Zhou et al. 2014) and the Taylor problem (Zhou et al. 2016). Zhou et al. (2018) have used the method to study a 3D line-tied plasma in the geometry of the Parker problem. The field they studied is an extension of the Taylor problem (Hahm and Kulsrud 1985) to a 3D line-tied geometry, the linear solution to which was given by Zweibel and Li (1987). This linear solution predicts a current layer of finite width, though Zhou et al. (2018) point out that with finite amplitudes the solution becomes pathological for a system length (in the z-direction) exceeding a critical value. They confirm that the current layer is finite and smooth at short system lengths, but are unable to study systems beyond the critical length owing to a numerical instability induced by extreme mesh distortion. They nonetheless speculate that the current layer may become singular for lengths exceeding the critical value. It should be noted that they restrict the numerical method to the reduced MHD equations (see Sect. 6.3).

5.2.4 Summary: computational ideal relaxation

A different approach again is taken by van Ballegooijen (1988a, 1988b), once more based on a Lagrangian approach. A discussion of that study is deferred to the following section. In summary, computational methods exist for analysis of equilibrium states for given magnetic topologies via ideal relaxation. However, each has its limitations, and in particular Lagrangian mesh methods are susceptible to numerical instability when the mesh distortion becomes too high. While these methods have been shown to be capable of capturing singular current features in general, they are yet to find singularity formation in the Parker problem geometry.

5.3 Computational approaches: the case for progressively thinner current layers

As demonstrated in the preceding sections, there remains no consensus on the formation of tangential magnetic field discontinuities in the Parker problem. However, in this section we describe studies which have established that any smooth equilibrium obtained in the Parker problem geometry will tend to contain current layers, and moreover that as the magnetic field lines become more tangled, these current layers become thinner and more intense. Thus, regardless of the existence of smooth equilibria, Parker’s initial idea—that field line tangling should lead to intense currents—is conceptually correct. In the following we make this statement more precise.

One of the first studies to demonstrate the progressive thinning of current layers was performed by van Ballegooijen (1988a, b), using a representation of the magnetic field in terms of Euler potentials; \({\mathbf{B}}=\nabla \beta \times \nabla \gamma \). In van Ballegooijen (1988a) it is described how application of a sequence of shear flows with random strength and direction on one boundary leads to an exponential decrease in length scales in the field line mapping between the two z-boundaries (length scales of the functions \(\beta \) and \(\gamma \) in this case). Then in van Ballegooijen (1988b) the equilibrium field corresponding to the state attained following each successive shear is calculated numerically by an iterative energy minimisation procedure (the topology being fixed since \(\beta \) and \(\gamma \) are held constant on the boundaries). The equilibria thus obtained are found to contain smooth current layers, though these current layers become progressively (exponentially) thinner as the number of shears increases. Similar results were obtained more recently in the reduced MHD framework by Rappazzo and Parker (2013).

Fig. 6
figure 6

Contours of \(\alpha =J/B\) in the midplane of the domain in the quasi-equilibrium attained after (left-right) 4, 8, and 12 iterations of a simulation. Images reproduced with permission from Mikić et al. (1989), copyright by AAS

Mikić et al. (1989) expanded upon the studies of van Ballegooijen (1988a, b) by employing a dynamical approach, using a simplified version of the MHD equations. In particular they ignored pressure and density fluctuations as well as the inertial term in the equation of motion, took resistivity to be zero, and imposed a large, spatially uniform viscosity. They applied a sequence of random shear flows on one line-tied boundary (the other being held fixed) of a cubic domain discretised with a mesh of \(64^3\) grid points. Consistent with van Ballegooijen (1988a, b) they found that the current developed filamentary structures, that these structures became exponentially thinner with time, and that the peak current density grew exponentially with time (see Fig. 6).

Fig. 7
figure 7

Top: isosurface of the current density modulus \(|{\mathbf{J}}|\) at 60% of maximum. Middle: \(|{\mathbf{J}}|\) in the plane \(z=0\). Below: \(\log _{10}(Q)\) in the plane \(z=-24\). For twist parameter given by a \(k=0.5\), b \(k=0.6\), c \(k=0.7\). Images modified from Pontin and Hornig (2015)

Van Ballegooijen (1988b) found that the thickness of current filaments in the equilibrium field was determined by the field line mapping. The reason why this must generally be the case was demonstrated by Pontin and Hornig (2015). Their argument goes as follows. In a force-free field satisfying \({\bf J}=\frac{1}{\mu_0}\nabla\times{\bf B}=\alpha {\mathbf{B}}\), we have \(\alpha \) constant along field lines (\({\mathbf{B}}\cdot \nabla \alpha =0\)). Thus, the length scales of \(\alpha \) must correspond to the length scales of the field line mapping (so long as \(\alpha \) is not constant). Now, \(\alpha ={\mathbf{J}}\cdot {\mathbf{B}}=J_\Vert /|{\mathbf{B}}|\), or \(J_\Vert =\alpha |{\mathbf{B}}|\). Therefore \(J_\Vert \) inherits the length scales of \(\alpha \) (therefore of the field line mapping), perhaps with a multiplying factor of the field strength. (However, in the Parker problem geometry the magnetic field strength is approximately uniform.) The argument can also be extended to fields close to force-free equilibrium by considering the correlation length of \(\alpha \) along field lines—see Pontin et al. (2016). Pontin and Hornig (2015) demonstrated this theory by considering the ideal relaxation of a set of braided fields constructed by superimposing localised regions of twist to a homogeneous background field (following the approach first taken by Wilmot-Smith et al. 2009a). By choosing different values of the twist parameter (k, see Fig. 7) they were able to consider progressively “more braided” fields. This can be quantified by examining quantitatively the “squashing factor” of the magnetic field, Q (Titov et al. 2002; Titov 2007), which measures stretching and squeezing under the field line mapping. As the field becomes more braided, the squashing factor exhibits progressively larger numbers of thinner layers with higher values of Q (“quasi-separatrix layers”), see also Wilmot-Smith et al. (2009b). Pontin and Hornig (2015) showed that the equilibria obtained following ideal relaxation contained both quasi-separatrix layers and current layers whose thickness scaled exponentially with the twist parameter (see Fig. 7). What’s more, the scaling exponents were found to match within the fitting uncertainty, reinforcing the theoretical arguments of a causal relation.

Thus, the fact that the thickness of current layers within the domain in an equilibrium is determined by length scales present in the field line mapping has been demonstrated with different theoretical arguments, and with computations in a range of magnetic configurations (Longcope and Strauss 1994; van Ballegooijen 1988a, b; Pontin and Hornig 2015). It naturally follows that, as the photospheric flows continually perturb the field line footpoints, and field lines become progressively more tangled, it is inevitable that the current layers within the domain will eventually become sufficiently thin and intense that rapid energy dissipation will occur. These arguments provide strong support for the weak form of Parker’s hypothesis as stated in Definition 3.

6 Reconnection onset and energy release

6.1 Flux braiding simulations

Parker’s idea that the corona is heated by dissipation of current layers formed in response to field line braiding sparked studies employing a wide range of approaches. One of these, not discussed so far, involves direct numerical simulation of the problem. We refer to these computations as “flux braiding simulations”. These studies are summarised in the excellent review by Wilmot-Smith (2015).

The philosophy of the majority of flux braiding simulations is fundamentally different from the ideal relaxation experiments described in the previous section. In particular, the focus shifts from the existence or otherwise of a smooth equilibrium, becoming rather whether formation of current layers (irrespective of whether they are finite or singular) and plasma heating ensues when time-dependent footpoint motions are applied. Instead of seeking an equilibrium structure and investigating its properties, the boundary driving is applied continually with the aim to reach a statistically-steady state whose properties are examined. In order to mimic the behaviour in the corona, characteristic timescales of the boundary flows should still be long compared with the Alfvén crossing time of the ‘coronal’ simulation domain—however it is worth noting that for computational expediency these timescales are typically not as well separated in the simulations as they are in reality.

There are various differences in the approach taken for the different flux braiding simulations, including between the equations solved (full 3D MHD with either zero or finite explicit resistivity, or reduced MHD, or variations of these), the form of the boundary motions applied (shear, rotational, etc.), the domain aspect ratio (from 1:1 to around 1:20), and time dependence of the driver.

6.2 Continually-driven systems: full MHD

6.2.1 Qualitative description of the energy dissipation

All of the flux braiding simulations described in this section share a common evolution: following the initiation of the boundary driving, the peak current grows within the domain, organising into thin layers that form and dissipate at locations throughout the domain, but that geometrically are always elongated in the z-direction along the loop, see Fig. 8. Reconnection within these current layers due to the finite resistivity (either explicit, hyperFootnote 2, or numerical) leads to dissipation of the magnetic energy, and a statistically steady state is generally reached, in which the injection of magnetic energy through the driven boundaries is balanced by the dissipation within the domain. We begin by discussing studies in which the full set of resistive MHD equations is solved (to be contrasted with reduced MHD studies in the following section). Notable variations of the numerical approach are mentioned during the exposition.

Fig. 8
figure 8

Illustrations of the statistically steady state attained in selected flux braiding simulations. a Current isosurfaces and magnetic field lines in the simulation of Galsgaard and Nordlund (1996) (for a unit cube). b Current isosurfaces in the simulation of Ng et al. (2012). c Current isosurfaces (left) and magnetic field lines (right) in the simulation of Rappazzo et al. (2007, 2008). Images reproduced with permission, copyright (a) by AGU; (b, c) by AAS

One of the first full-MHD flux-braiding studies was performed by Galsgaard and Nordlund (1996) who considered shear boundary flows of random phase, duration and velocity (on both boundaries, between 2 and 40% of the Alfvén speed; see illustration in Fig. 9a). Most of their simulations were conducted in a unit cube, although they also considered simulations where the box length in the z-direction was increased by a factor of 10. The authors used a hyper-resistivity coefficient that depended of the grid resolution, and they considered resolutions between \(24^3\) and \(136^3\). Some typical current distributions from these simulations are shown in Fig. 8a.

Fig. 9
figure 9

Illustrations of different photospheric drivers in flux braiding simulations. a Schematic of the effect of one (left) and two (right) shear velocity perturbations on the field lines in the domain. b Streamfunction contours for one time instant of the time-dependent driving flow employed by Longcope and Sudan (1994) and Ng et al. (2012), modified from Longcope (1993). c Contours of the streamfunction for the driving flow of Rappazzo et al. (2007, 2008)

One of the main findings of these simulations is that the free magnetic energy in the domain (energy in excess of the unperturbed, uniform field) depends strongly on the mean driving velocity: when this is 2% of the background Alfvén speed \(v_A\) the free energy is only 1.5% of the background, while it is 45% when the mean driving speed is 20% of \(v_A\). By contrast, the energy dissipation is found to be essentially independent of the resistivity (recall that this is a hyper-resistivity determined by the grid resolution). Based on a scaling law derived for the energy dissipation as a function of the driving timescale the authors conclude that the energy injection/dissipation is sufficient to heat the corona. It is worth noting that these simulations were run with a comparatively high plasma-\(\beta \) that was initially 0.5, and progressively increased in time due to the absence of radiative losses in the energy equation.

6.2.2 Turbulent nature of the statistically steady state

Another flux braiding simulation was performed by Hendrix and Van Hoven (1996)—they simplified the MHD equations by assuming \(\rho =p=\mathrm{constant}\) throughout, and employed time-dependent boundary flows in rotational cells, with a box aspect ratio (length:width) of 4:1 and grid resolutions up to \(256^2\times 31\) (the x- and y-directions being solved using a spectral method). Their main result was the finding that the energy dissipation process in the volume is due to a turbulent cascade to small scales. Based on magnetic energy spectra at the highest resolutions they identified an inertial range that flattened with time, approaching (but always remaining steeper than) the classical slope of \(k^{-3/2}\), k being the horizontal wavenumber (Kraichnan 1965).

Many other papers have studied the properties of the turbulence induced in the domain within flux braiding simulations, however the majority of these have solved only the reduced MHD equations (see the following section). One exception is the study of Dahlburg et al. (2012), who included substantially more physics in the system than Hendrix and Van Hoven (1996) by allowing \(\rho \) and p to vary, and including both thermal conduction and radiative losses in the energy equation. The grid resolution is \(128^3\) where again the x- and y-directions are solved using a spectral method, and the box aspect ratio is 5. With this additional physics included, they were able to show for the first time that the thermal energy of the plasma, like the magnetic and kinetic energies, attains a state that is statistically steady, albeit highly fluctuating in time and space. They emphasised that the temperature of the plasma in the domain (after the statistically steady state is reached) is highly structured, showing high spatial and temporal intermittency.

6.2.3 Dependence on driving flows

Magnetic flux braiding simulations have used a wide variety of boundary driving flows patterns. These usually take the form of sequences of rotations or shears, randomised within certain parameter ranges. These flow patterns are chosen largely for computational tractability—they are divergence-free and therefore do not lead to additional complications of plasma compression/rarefaction at the boundaries (moreover for the reduced MHD simulations, the flow is divergence-free by construction—see the following section). In studies that considered different driving patterns, the overall behaviour of the system is typically found to be much the same irrespective of the particular driver, and so it is tempting to conclude that the spatial pattern of the driving is largely unimportant in determining the coronal evolution.

However, Ritchie et al. (2016) showed by considering a more systematic driving that the coronal evolution can depend strongly on the driving pattern (supporting earlier ideas of Wilmot-Smith et al. 2011). In particular, they considered a driving flow given by a pair of vortices that ‘blink’ on and off in time in anti-phase with one another. When the vortices have the same sense of rotation the helicity injection has a preferential sign, while if they have opposite sense the net helicity injection (over many cycles) is zero. They showed that magnetic energy dissipation was very different between these two cases. When the net helicity injection is zero the magnetic energy in the domain shows only small fluctuations about the steady state, while for non-zero net helicity injection the energy tends to build slowly before exhibiting much larger sudden drops (that one might hypothesise are related to the onset of some large-scale kink-type instability in the domain). Furthermore, for the case with zero net helicity injection, Ritchie et al. (2016) also considered different spatial positioning of the vortical flows, and found that this affected the long-time-averaged magnetic energy in the domain. In particular, driving patterns that induce a higher degree of field line tangling (per unit time, say, as measured by the topological entropy, see Thiffeault 2010) result in higher stored energy. In a similar study, Knizhnik et al. (2019) examined the effect of helicity injection in a system with many (61) driving vortices imposed on the line-tied z-boundaries of a simulation domain of size \(3\times 3\times 1\). Like Ritchie et al. (2016) they found that the field in the domain developed large-scale twist in a shell surrounding the portion of the field that is driven, this being gradually produced by many small-scale reconnections within the volume. The evolution was not followed for as many driver periods as Ritchie et al. (2016), so no kink-type instability was observed. They found that the Poynting flux into the domain was greater for higher helicity injection, with similar overall heating in both cases but large-scale twist retained within the domain for high helicity injection rate. These results together indicate that certain characteristics of the driving flow may have a substantial impact on the resulting structure of the coronal field, and the distribution of energy release events.

6.3 Continually-driven systems: reduced MHD

6.3.1 Reduced MHD

Owing to the substantial reduction in complexity of the equations, many studies of the Parker problem, and many flux braiding simulations in particular, have made use of the reduced MHD equations (hereafter RMHD, see Kadomtsev and Pogutse 1974; Strauss 1976). This reduction in complexity allows analytical models to be developed, and also substantially reduces the computational expense, the latter being particularly beneficial for early studies when computational power was limited. The justification for this simplification is that the length of a typical coronal loop is much larger than its diameter. One therefore makes the following ansatz for the ordering of terms:

$$\begin{aligned} B_z\sim \nabla _\perp \sim 1, \quad |{\mathbf{B}}_\perp |\sim \frac{{\partial }}{\partial z}\sim \epsilon , \end{aligned}$$
(29)

where \(\epsilon \ll 1\). Following this ordering through the equations one finds that \(v_z\), variations of \(B_z\), and the pressure p are all of order \(\epsilon ^2\), so that p and \(v_z\) are absent from the leading order equations to be considered, and \(B_z=B_{z0}\), constant. Furthermore, \(\nabla \cdot {\mathbf{v}}=0\), so one can assume a constant density. Thus

$$\begin{aligned} {\mathbf{B}}=B_{z0}{\mathbf{e}}_z+{\mathbf{e}}_z\times \nabla _\perp \psi , \quad {\mathbf{v}}={\mathbf{v}}_\perp ={\mathbf{e}}_z\times \nabla _\perp \phi \end{aligned}$$
(30)

and to leading order the induction and momentum equations are

$$\begin{aligned} \rho _0\left( \frac{{\partial \varOmega }}{\partial t}+{\mathbf{v}}_\perp \cdot \nabla \varOmega \right)= & {} {\mathbf{B}}\cdot \nabla J_z+\nu \nabla _\perp ^2\varOmega , \end{aligned}$$
(31)
$$\begin{aligned} \frac{{\partial \psi }}{\partial t}+{\mathbf{v}}_\perp \cdot \nabla \psi= & {} -B_{z0}\frac{{\partial \phi }}{\partial z}+\eta \nabla _\perp ^2\psi , \end{aligned}$$
(32)

where \(\varOmega =\nabla _\perp ^2\phi \) is the z-component of the vorticity, and \(J_z=\nabla _\perp ^2\psi \) (see, e.g., Biskamp 1993). Due to the absence of pressure and density fluctuations, the only wave modes in the system are shear Alfvén waves, that propagate in the z-direction. Moreover, there is no coupling to an energy equation, so that any energy dissipated by viscous or resistive effects is immediately lost.

6.3.2 Dependence of the statistically steady state on \(\eta \)

The first direct simulations of flux braiding in the RMHD framework were carried out by Longcope and Sudan (1994). They applied a driving velocity composed of a ‘random’, time-dependent distribution of vortices intended to mimic the incompressible component of photospheric convection (see Fig. 9b), discretising the domain with up to \(48\times 48\) Fourier modes in xy, and 10 finite difference intervals in z. They confirmed the formation of intense, filamentary current concentrations in response to the applied driving. They examined the dependence of the resulting statistically steady state on the (spatially uniform) diffusivity \(\eta \), over a range \(\eta =10^{-3}\)\(6\times 10^{-3}\). This they accompanied with analysis of a model for the statistically steady state based on an ensemble of Sweet–Parker current sheets (Sweet 1958; Parker 1957). They concluded, both from the model and the simulations, that both the typical perpendicular field strength in the loop, \(\overline{B}_\perp \), and the energy dissipation rate, scale as \(\eta ^{-1/3}\). They noted, however, that extrapolation to coronal parameters would imply a free magnetic energy (\(\sim \overline{B}_\perp ^{\; 2}\)) in excess of the potential field energy. This being physically implausible, they noted that the rather strong inverse scaling with \(\eta \) must break down for some smaller value of \(\eta \) than those considered. Motivated by these findings, Ng et al. (2012) repeated the study of Longcope and Sudan (1994), but with advances in computing power were able to extend the range of diffusivity down to \(\eta =1.6\times 10^{-4}\) (this requiring grid resolution \(1024^2\times 128\)). A typical current distribution from their simulations is shown in Fig. 8b. They confirmed in their simulations that in the statistically steady state the total magnetic energy, kinetic energy and dissipation in the domain scale inversely with \(\eta \), see Fig. 10a, b. They also found that, as hypothesised by Longcope and Sudan (1994) the scalings of both the energy dissipation and perpendicular field strength exhibit a break, such that for \(\eta \le 10^{-3}\), a much shallower scaling around \(\eta ^{-1/5}\) is obtained. They argue that the reason for this break in the scaling is that for \(\eta \approx 10^{-3}\) the energy dissipation time, \(\tau _E\), is comparable to the driver correlation time, \(\tau _c\). They further argue that for \(\tau _E>\tau _c\) the energy dissipation rate should become independent of \(\eta \), using an argument based on random walk of field line footpoints (see also Ng and Bhattacharjee 2008). One note of caution related to the simulation scaling results of Ng et al. (2012) is that for \(\eta \) values smaller than the break point in the scaling, \(|{\mathbf{B}}_\perp |\approx 0.5B_z\). Thus \(|{\mathbf{B}}_\perp |\) is larger than would be expected for applicability of the RMHD approximation, and so extreme caution should be exercised in interpreting these results.

Fig. 10
figure 10

a Total magnetic energy for the flux braiding simulations of Ng et al. (2012) with \(\eta =3\times 10^{-4}\) (blue) and \(\eta =5\times 10^{-3}\) (red), together with total kinetic energy for \(\eta =3\times 10^{-4}\) (green) and \(\eta =5\times 10^{-3}\) (orange). b Total joule dissipation for \(\eta =3\times 10^{-4}\) (blue) and \(\eta =5\times 10^{-3}\) (red). c Magnetic and kinetic energies for the flux braiding simulations of Rappazzo et al. (2007, 2008) with \(\eta =1.25\times 10^{-3}\) and \(v_A/u_{ph}=200\). d Total dissipation (ohmic plus viscous) for simulations with \(\eta =1/R\) as marked (based on standard second-order diffusion), together with a simulation using eighth-order dissipation terms, marked \(R_4\). Images reproduced with permission (a, b) from Ng et al. (2012), and (c, d) from Rappazzo et al. (2007, 2008), copyright by AAS

6.3.3 Intermittency in the statistically steady state

A substantial body of work addressing the Parker problem in the RMHD framework has focussed on the turbulent nature of the statistically steady state. The first such studies were performed by Einaudi et al. (1996) and Dmitruk and Gómez (1997). They in fact solved a 2D problem, replacing the \(\partial /\partial z\) terms in the RMHD equations with forcing functions. Einaudi et al. (1996) found that the average energy dissipation increases when \(\eta \) is reduced (the numerical resolution is increased), as do both the peak dissipation and degree of intermittency. Dmitruk and Gómez (1997) focus on this intermittency, and identify a scaling in the number N of energy release events with a given energy E within their simulation of the form

$$\begin{aligned} \frac{dN}{dE}\sim E^{-1.5} \end{aligned}$$
(33)

(individual energy release events being identified as peaks in a plot of volume integrated energy dissipation vs. time). They suggest that it is significant that this measured relation compares favourably with statistics of observed flare energies (e.g., Crosby et al. 1993; Wheatland 2008). These results were extended to fully-3D RMHD by Dmitruk and Gómez (1999) and Gómez et al. (2000).

6.3.4 Inertial range exponents in the statistically steady state

In order to probe the turbulent nature of the statistically steady state in the fully-3D RMHD problem, Dmitruk and Gómez (1999) analysed spectra of the magnetic energy \(E_m\) and kinetic energy \(E_k\) as functions of the in-plane wavenumber \(k_\perp =k_{xy}\). They found a magnetic energy spectrum \(E_m\sim k_\perp ^{-3/2}\), consistent with a Kraichnan-type spectrum, with the kinetic energy spectrum being a little shallower. They repeated their simulations for different parameter regimes, characterised by different ratios of \(\tau _A/\tau _p\), where \(\tau _A\) is the Alfvén travel time along the loop and \(\tau _p\) is the photospheric turnover time (analogous to the driver correlation time considered by Ng et al. 2012). Their results indicated that the total energy dissipation scales as \((\tau _A/\tau _p)^{3/2}\), consistent with the results of the 2D study of Dmitruk and Gómez (1997). Dmitruk and Gómez (1999) considered a range for \(\tau _A/\tau _p\), but with \(\tau _A/\tau _p\le 0.2\). This range was extended to larger values in a follow-up study by Dmitruk et al. (2003). They found for \(\tau _A/\tau _p=0.5\) that \(E_m\sim k_\perp ^{-2}\), and for \(\tau _A/\tau _p>0.5\) that \(E_m\sim k_\perp ^{-5/3}\). These results indicate the delicate nature of the problem at hand, though are probably not directly relevant to the heating of coronal loops, where it is expected that \(\tau _A/\tau _p\ll 1\).

The ideas discussed above were picked up by Rappazzo and co-workers. Their computational approach is similar to that of Dmitruk and co-workers, with the boundary driving involving a steady flow comprised of a superposition of vortices. However, they utilised both uniform resistivity/viscosity and hyper-resistive and hyper-viscous terms in the equations to extend the ranges of accessible (effective) Reynolds and magnetic Reynolds numbers. The domain aspect ratio is 1:10. (Note that since the driving flow is steady, these simulations are in the regime \(\tau _A/\tau _p=0\).) The first study of these simulations was reported by Rappazzo et al. (2007, 2008). They ran simulations with a ratio of Alfvén speed to typical photospheric driving speed, \(v_A/u_{ph}=\{50,200,400,1000\}\), finding that the magnetic energy spectrum steepens as the ratio increases, in agreement with Dmitruk et al. (2003). In analysing the energy dissipation rate, they argue that the relative sizes of the magnetic field strength and loop length are critical (in shorter loops the effect of the line-tying-enforced magnetic tension force is (on average) enhanced, weakening non-linear interactions). Based on dimensional analysis they assert that the energy dissipation should scale with the variable

$$\begin{aligned} f=\frac{l_\perp v_A}{L u_{ph}}. \end{aligned}$$
(34)

On this basis they propose a coronal heating scaling/dissipation rate

$$\begin{aligned} \epsilon \sim l_\perp ^2\rho v_A U_{ph}^2 \left( \frac{l_\perp v_A}{L u_{ph}} \right) ^\frac{\alpha }{\alpha +1} , \end{aligned}$$
(35)

where \(\alpha \) is the scaling index for the turbulence model. In their simulations they find strong turbulence (with Kolmogorov-like spectra with exponent \(-5/3\) in the perpendicular direction) for weak axial magnetic fields and long loops, while for strong axial magnetic fields and/or short loops a weaker turbulence regime is found (characterised by steeper spectral slopes of total energy). Similar to Galsgaard and Nordlund (1996) they argue that the energy dissipation should become independent of the magnetic Reynolds number for \(\eta \le 10^{-3}\). This assertion is made on the basis of the close comparison between the volume integrated dissipation for two of their simulations, in one of which they used second-order diffusion with \(\eta =1.25\times 10^{-3}\) and in the other of which they used eighth-order dissipation with a coefficient of \(10^{-19}\), see Fig. 10(d).

In a follow-up study, Rappazzo et al. (2010) highlighted that the properties of the turbulent dissipation in the loop are independent of the driving flow, by applying a continuous simple shear driving flow. In those simulations a single, monolithic current layer first forms and then breaks up in response to a tearing instability, which is followed by a transition to turbulence.

6.4 Summary: results of flux braiding simulations

While there are many variations in the approach taken in flux braiding simulations, the results build an overall picture of the system behaviour, as follows.

  • Following the initiation of the driving velocity, thin ribbons of current rapidly form and then dissipate throughout the domain, and a statistically steady state is reached.

  • Both reduced MHD and full MHD simulations indicate that the non-linear evolution can be described by a turbulent cascade (though see Klimchuk et al. 2010). The turbulence is anisotropic, and its properties (such as the scaling exponents in the inertial range) depend on the loop length and axial field strength, as well as potentially other properties associated with the driving velocity.

  • The free magnetic energy (in excess of the potential field energy) dominates over the kinetic energy, by a factor \(>10\) (see, e.g., Fig. 10a, c).

  • The free magnetic energy, or average \({\mathbf{B}}_\perp \), scales with driving speed. It also scales with the magnetic Reynolds number, at least for \(R_m\le 10^3\).

  • The average energy dissipation rate scales with driving speed, or with \(\tau _A/\tau _p\), but is (approximately) independent of the magnetic Reynolds number, at least for \(R_m\) sufficiently large \(\ge 10^3\). However, caution should be taken in interpreting these results, since there remains a huge disparity between the numerically accessible values of the magnetic Reynolds number and the true solar values.

  • The geometrical pattern of the driving flow has minimal effect on the energy dissipation. However, averaged quantities such as the helicity injection rate or correlation/turnover time may affect the characteristics of the energy release events in the domain.

6.5 Resistive relaxation simulations

The flux braiding simulations described so far in this section all achieve a (time-averaged) balance between energy input through continual boundary driving and dissipation in the volume. However, the discrepancy between the simulation parameters (particularly \(R_m\)) and those in the corona are large, meaning that implications for the corona must be inferred through often-uncertain extrapolation. This motivates the complementary approach taken in a series of papers by Wilmot-Smith et al. (2010, 2011) and Pontin et al. (2011), who consider a turbulent relaxation (in the absence of boundary driving) in a magnetic flux tube starting from an initially braided state, which can be viewed as generated by an ideal dynamics (\(R_m = \infty \)) over a certain time before reconnection sets in. This approach permits the details of the energy release process following onset of reconnection to be studied, which is relevant if—at coronal parameters—the relaxation/energy dissipation time is short compared with the timescale for field line tangling/energy injection.

The initial braided magnetic field used for this series of papers is taken as the final state of the ideal relaxation of Wilmot-Smith et al. (2009a) (see Sect. 5.2), and the authors considered a resistive MHD evolution with line-tying of \({\mathbf{B}}\) and zero flow on the boundaries. They found that following the onset of reconnection in a pair of initially forming current layers (those shown in Fig. 7), a cascade to smaller scales throughout the domain ensues. The overall result is an “unbraiding” of the field in the loop, such that the field line mapping after the resistive relaxation is much simpler than before, this unbraiding being achieved through reconnection in a large collection of discrete current layers (Pontin et al. 2011)—see Fig. 11. The presence of a turbulent inertial range, with a power-law slope of \(-5/3\) for the magnetic energy, is demonstrated by Pontin et al. (2016). The number of current layers formed during the relaxation, and time taken to unbraid the field, are both shown to increase with the magnetic Reynolds number. Even at the moderate value of \(R_m\) accessible in the simulations (100–5000), the energy release process lasts 5–8 loop Alfvén-crossing times (length-wise), with this likely to increase substantially at coronal values of \(R_m\). This energy release timescale is an important parameter when comparing with loop heating observations (Reale 2014). Another noteworthy result of these simulations is that not all of the energy stored in the initial state as twisting/braiding of field lines is released during the turbulent relaxation. Application of Taylor’s hypothesis in which the energy is minimised subject to the constraint of preserving the total helicity (Taylor 1974) would lead one to predict the final state of relaxation to be a uniform field. However, it turns out that only around 60% of this stored magnetic energy is released, with the rest remaining as large-scale twist within the domain (Pontin et al. 2011)—see the plots of \(\alpha ^*\) in Fig. 11. This led to the discovery of additional constraints on the relaxation—see Yeates et al. (2010).

Fig. 11
figure 11

Selected field lines (traced from fixed boundary footpoints), current isosurfaces, and contour maps of \(\alpha ^*=(1/48)\int _{z=-24}^{z=24} {\mathbf{J}\cdot \mathbf{B}}/{B^2}\,dl\), integrated along field lines and displayed on the plane \(z=-24\). The Alfvén speed in the non-dimensional simulations is 1, giving a loop Alfvén crossing time of 48 units. After Pontin et al. (2016), based on the simulations of Pontin et al. (2011)

Similar relaxation simulations are described by Rappazzo and Parker (2013)—they studied RMHD simulations starting from an initial non-equilibrium field constructed by superimposing a perturbation component \({\mathbf{B}}_\perp \) determined by a set of arbitrary Fourier modes onto a uniform field along the loop, \(B_z\). The boundaries are line-tied with zero flow. Like the above-mentioned studies, they observe a relaxation taking the form of decaying turbulence—this turbulent relaxation being triggered so long as the mean value of \(|{\mathbf{B}}_\perp |\) is greater than \(0.04B_z\). The relaxation again leads towards a final state containing regions of finite twist at large scales, rather than relaxing all the way to the homogeneous field.

6.6 Plasma response and consistency with observations

The primary focus of this article is the magnetic field evolution in the Parker problem. However, another crucial ingredient in determining what role (if any) Parker’s braiding mechanism plays in heating the corona is the plasma’s thermodynamic response to the energy deposited following the onset of reconnection. It is generally computationally infeasible to model the full plasma thermodynamics alongside the 3D magnetic field evolution. Instead, the plasma response is usually modelled in (ensembles of) 0D or 1D loops, with parameterised “nanoflare” energy deposition used to mimic the effect of heating events induced by magnetic reconnection (for example). The reader is directed to the comprehensive review by Reale (2014). The importance of combining theoretical work using different approaches for understanding coronal loop observations was recently highlighted by Klimchuk (2015).

Notable exceptions to the above—in which substantial thermodynamics are included alongside the 3D magnetic field evolution—are the flux braiding studies of Dahlburg et al. (2016, 2018); Rappazzo et al. (2018). Those authors have solved the full MHD equations including thermal conduction and optically thin radiation in the energy equation. They employed a driving flow composed of a superposition of vortical flows, similar to the approach of, e.g., Hendrix and Van Hoven (1996), see Fig. 9(c). In common with other flux braiding studies, Dahlburg et al. (2016) find a turbulent evolution ensues, characterised by current layers elongated along the loop. Magnetic energy dissipation and subsequent plasma heating in these current layers leads to a multi-thermal plasma distribution, with “strands” of hotter and cooler plasma distributed irregularly throughout the domain, on perpendicular length scales shorter than those presently observable (based on a dimensionalisation in which the simulated domain has size \(4\times 4\times 50\) Mm). The authors use simulated densities and temperatures to reconstruct observables such as the differential emission measure, which are argued to compare favourably with observations of active regions (see Fig. 12a). These results are extended by Dahlburg et al. (2018) to consider different loop lengths and magnetic field strengths. They find that the highest temperatures are reached for the shortest loops, and that stronger magnetic fields mean broader emission-measure distributions. Rappazzo et al. (2018) describe similar simulations, varying the ratio \(\tau _A/\tau _p\) (see Sect. 6.3). They conclude that for \(\tau _A/\tau _p\ll 1\) a broad distribution of temperatures is present in the loop, while for \(\tau _A/\tau _p= 1\) the temperature distribution is much narrower.

Fig. 12
figure 12

a Modified from Dahlburg et al. (2016). Synthesised emission from their continuously driven ‘D’ simulations at resolution \(128^2\times 144\). b Modified from Pontin et al. (2017). Synthesised emission from their relaxation simulations, resolution \(360^2\times 240\)

Another study to include thermal conduction and optically thin radiation, and to synthesise observables, was carried out by Pontin et al. (2017) (see Fig. 12b). They considered decaying turbulence during the relaxation of an initially braided magnetic flux tube, as in the simulations of Wilmot-Smith et al. (2010); Pontin et al. (2011) described in Sect. 6.5. They aimed to address the issue that, while some recent observations of coronal loops reveal distinguishable crossing bright strands (the previously assumed visible signature of magnetic braiding) the majority of observed loops do not (but rather appear to have either no substructure or to be composed of “well-combed” strands). They showed that clear crossing of heated strands appeared only in certain cases, in spite of the fact that underlying field was braided in all cases considered. The conclusion is then that an absence of a clear “braided” appearance in which bright strands cross one another does not exclude that the underlying field lines are braided. This is seen also in the simulations of Dahlburg et al. (2018), where no crossing strands are present in the synthesised emission (Fig. 12a). It was recently shown by Pontin et al. (2020) that the spectral properties of synthesised line profiles from simulations of braiding-induced turbulence are consistent with key properties of observed spectra from coronal emission lines, including the magnitude of the non-thermal broadening and its scaling with intensity as well as the non-Gaussian nature of the spectra.

The studies described in this section represent only a start in determining the consistency of the braiding mechanism with observations of the hot corona. In future, coupling with the lower atmosphere through the all-important transition region, perhaps by linking more closely with the one-dimensional loop models that include more realistic thermodynamics, will be a key step.

7 Implications of coronal magnetic complexity

7.1 Magnetic topology of the corona: formation of current sheets

Parker’s braiding hypothesis for coronal heating considers the simplest case in which the magnetic field extends from one positive polarity to one negative polarity on the photosphere. However, it is now well-established that the coronal magnetic field has a highly complex structure, with photospheric flux concentrations present at all observable scales, following a power-law distribution with flux (Parnell et al. 2009). What’s more, these photospheric flux fragments are continually emerging, cancelling, splitting, and merging, with estimates that all of the coronal flux is ‘recycled’ (reconnected) on average on timescales of hours (Close et al. 2005; Bellot Rubio and Orozco Suárez 2019).

The complex array of photospheric flux fragments means that the coronal magnetic field must be pervaded by a web of null points, their associated separatrix surfaces and separator lines, and quasi-separatrix layers (Schrijver and Title 2002; Longcope and Parnell 2009; Platten et al. 2014). These structures—which represent either discontinuities (null points, separatrices, separators) or strong gradients (quasi-separatrix layers, hyperbolic flux tubes) in the field line mapping—are all potential sites for accumulation of intense currents, and thus magnetic reconnection. Two-dimensional magnetic null points (X-points) are well known to ‘collapse’ in response to certain external perturbations (Dungey 1953; Imshennik and Syrovatsky 1967). In a line-tied geometry, slow boundary motions lead to a sequence of approximate equilibria which in the ideal limit contain current sheets at the null that are tangential discontinuities with singular current density (Green 1965; Syrovatskii 1971). The same is true of the three-dimensional magnetic null points present in the solar corona (Klapper et al. 1996; Pontin and Craig 2005; Fuentes-Fernández and Parnell 2013), as well as separator field lines (Longcope and Cowley 1996). Quasi-separatrix layers are additional sites where intense currents may accumulate, in response to certain boundary flows, though most likely these currents form in layers of finite thickness even in the ideal limit (see Démoulin 2006, for a review). All of the above has led Priest et al. (2002) to propose a modification to Parker’s braiding scenario, called “flux tube tectonics”, in which the current sheet formation—and subsequent reconnection and heating—take place at the web of null points, separatrices and separators threading the corona.

A number of authors—either explicitly or implicitly—made use of the fact that singular current sheets form in the ideal limit at magnetic nulls to support Parker’s hypothesis of spontaneous formation of tangential discontinuities. We have excluded these from the discussion of Sects. 4 and 5 since they do not conform to the geometry of the original Parker problem, in which all field lines connect from one boundary plate to the other and magnetic nulls are thus excluded. A number of authors considered magnetic field equilibria above a plane, in a configuration that begins with a bald-patch (Titov et al. 1993) at the photosphere (a null point being located beneath the photosphere), initially proposed by Moffatt (1987). A potential field extrapolation based on the final flux distribution would contain a null point in the corona. However, applying the topological constraints of an ideal evolution between the initial and final states a current sheet is instead present, extending from the photosphere up to the location at which the X-point would be located. A number of variants of this problem have been studied, including different boundary perturbations and extension to spherical geometry (Low 1987, 1989, 1991, 1992; Vekstein et al. 1991). These results lend additional weight to the notion of current sheet formation in the corona in general, but do not directly address the Parker problem as stated in Sect. 2. Further studies of this type, in which magnetic nulls and/or separatrices are present include the studies by (Kumar et al. 2014, 2015).

7.2 Large-scale 3D MHD simulations of the corona

Around the beginning of this century, computing power became sufficient to permit MHD simulations of substantial parts of the corona. The philosophy of these models is to reproduce as faithfully as possible the observed structure of the corona by including as much physics as is computationlly feasible. This requires realistic pressures and temperatures, and thus in addition to solving the induction equation together with conservation of mass and momentum, these models must solve an energy equation that accounts for thermal conduction and radiative losses. Some such models use a parameterised heating function, but in most the heating is directly by ohmic dissipation, consistent with Parker’s braiding mechanism. It is only possible to resolve the convective driving in smaller simulations. Larger-scale simulations use a velocity driver that mimics some key properties of observed photospheric flows. For a comprehensive review of the successes and caveats of these models, see Peter (2015).

Fig. 13
figure 13

Simulated TRACE 171 (top) and 195 (middle) images, together with a simulated magnetogram (bottom), showing three identified bright loops. Image reproduced with permission from Gudiksen and Nordlund (2005b), copyright by AAS

The first models of this type were presented by Gudiksen and Nordlund (2002, 2005a, b) (see Fig. 13). They demonstrated the formation of a loop-dominated corona in their simulations, concluding that a DC-type heating consistent with Parker’s braiding mechanism was responsible for creating these few-million-degree loops. They do note that the simulated region is rather smaller than a typical active region, meaning that they do not produce loops as long as in the corona, meaning for example that loop over-densities (with respect to the surroundings) are not as extreme as those observed. In that case the photospheric/boundary driving flow was based on a model that mimics some properties of the incompressible component of photospheric granulation.

Bingert and Peter (2013) made a detailed study of individual energy release events in one of these models. They found that the heating occurred in transient bursts whose energy content followed a power-law distribution. They also found that the energy flux into the coronal plasma peaks for event energies around \(10^{17}\) J, consistent with the energy of a nanoflare proposed by Parker (1988). Kanella and Gudiksen (2018) performed a similar study, employing an intensive and sophisticated algorithm to pick out individual (ohmic) heating events in space and time. They too found a power law distribution of heating events in both energy and duration (but not in event volume). In general, the height distribution of the energy deposition and temperature typically matches the observations quite well, with the energy input per particle (calculated by dividing the energy input by the mass density) peaking somewhere near the transition region (Hansteen et al. 2010; Bingert and Peter 2011).

Fig. 14
figure 14

Comparison of observed emission and synthesised emission in a large-scale, data-driven MHD simulation of an active region. Note that the peak count value in the observations (corresponding to 1.00) is 3500 DN/pixel, which is a factor of six higher than in the model. Image reproduced with permission from Warnecke and Peter (2019), copyright by the authors

Warnecke and Peter (2019) took this simulation approach one step further by employing a time sequence of observed magnetograms to update the lower boundary condition for the magnetic field. This is combined with a driving flow to mimic granulation of the type used by Gudiksen and Nordlund (2002, 2005a, b). The aim in this “data-driven” simulation is to see whether the observed coronal emissions are reproduced (previous studies having shown that the coronal-type emissions could be reproduced in a general sense). It turns out that the characteristics of the coronal emission are rather well reproduced (see Fig. 14), with visible loops in many of the same locations, super-imposed upon a diffuse background corona. The authors speculate that the discrepancies observed (e.g., lack of small-scale structure in the active region core, under-estimate of the overall shear in the loops) could be improved in future by increasing the resolution of the underlying magnetograms, as well as making use of the observed horizontal components of the magnetic field. The fact that the observed coronal emission is well reproduced by such simulations driven by footpoint shuffling is a strong indication that the braiding/flux-tube tectonics model is a key part of the heating process.

What is perhaps surprising is that these large-scale MHD simulations reproduce so well the energy fluxes and temperatures required to maintain a hot corona, despite the fact that they drastically under-resolve the braiding process itself. This has led to speculation that the plasma heating is relatively insensitive to the details of the mechanism that extracts the energy from the coronal field. In the larger-scale models, single current layers typically fill entire ‘loops’ in the coronal volume. Since small spatial and temporal scales are missing from the modelling, the heating is of a ‘DC’ type consistent with Parker’s original ideas—it is likely that AC or wave heating would make a larger contribution if these smaller spatio-temporal scales were included. It is worth noting that, while the overall heating rates are reproduced in these models, many more detailed features of the observations are not: for example the non-thermal broadening of synthesised spectral lines is much smaller than typically observed.

8 Summary: established results and open questions

Ever since it was first proposed by Parker (1972), the idea of “topological dissipation” has stimulated significant debate. While the mechanism is at the heart of one of the major coronal heating theories (nanoflare heating), there remain significant open questions. The problem of proving existence or non-existence of equilibria in the geometry of the Parker problem is so challenging mathematically not least due to difficulty in specifying the field topology analytically (see the discussion of Low 2010b). Computationally, eliminating numerical dissipation (equivalently preserving the field topology) as well as distinguishing finite current layers from tangential discontinuities on a discrete grid, provide significant difficulties. Thus, the existence of smooth equilibria for arbitrary topology—and accessibility from an initially smooth field via smooth boundary motions—remains an open problem.

In the scenario addressed by Parker’s hypothesis, we begin from a smooth magnetostatic equilibrium, and perform smooth boundary motions (Sect. 2). The original claim by Parker (1972) was that no smooth equilibrium could in general be found even when the boundary field line footpoint displacement was small. However, van Ballegooijen (1985) demonstrated an error in Parker’s calculation, further arguing that a smooth equilibrium should always exist in this case (Sect. 4.2). Assessment of the case of finite-amplitude footpoint perturbations typically requires a computational approach. While various arguments have been proposed in support of Parker’s hypothesis (Sect. 4.3), there has been no concrete demonstration of formation of tangential discontinuities. On the other hand, while certain classes of fields have been shown to exhibit smooth equilibria, there remains no general statement. Thus far, Parker’s hypothesis in the form of either Definition 1 or 2 remains to be proven or disproven.

In spite of the aforementioned lack of progress in determining the existence or otherwise of tangential discontinuities in Parker’s scenario (as per Definitions 1 and 2), there is now overwhelming evidence that Parker’s original idea was conceptually correct. Specifically, the weak form of the hypothesis (Definition 3) has been proven by a number of approaches: continual braiding of the coronal field lines by boundary motions leads to equilibria containing current layers that grow thinner and more intense at an exponential rate. This must inevitably lead to the onset of reconnection and energy dissipation. The larger the magnetic Reynolds number, the longer it will take to reach the dissipation scale, but it must be reached in the end.

On the basis of the above, it is clear that magnetic flux braiding can heat the solar corona. However, there remain many gaps in our understanding that must be filled before we can assess the efficiency of the mechanism, and therefore whether it provides a substantial contribution to heating the corona in practice. Flux braiding simulations in the geometry of the Parker problem show qualitative behaviour that is robust between the different approaches (as summarised in Sect. 6.4). However, there remain significant steps to be taken in applying these results to the corona: (i) the scaling of various quantities with parameters such as the magnetic Reynolds number remains uncertain, (ii) the relative timescales of energy storage (field line tangling) and energy release are crucial for the nature of the heating, but are not well known (the former due to observational limitations and the latter at least in part to the aforementioned magnetic Reynolds number dependence), (iii) the effect of the huge additional complexity of the chromosphere, in terms of magnetic structure, plasma structure, and additional physics beyond MHD. What we do know is that even in the corona there are many likely sites of current sheet formation, reconnection, and energy dissipation that are excluded from the problem as posed by Parker (Sect. 7.1).

Finally, global simulations of the corona show a surprisingly good agreement with observations (given the difference in parameter regimes). Plasma is heated in these models primarily by dissipation of current sheets induced by footpoint shuffling, and the resultant coronal loops successfully reproduce many properties of the observations (Sect. 7.2). This provides evidence that the braiding mechanism really does provide some appreciable contribution to heating the corona. Parker’s two linked ideas of topological dissipation and nanoflare heating of the corona by field line braiding have stimulated a wide range of studies that have in turn greatly advanced our knowledge of the Sun and of plasma physics processes more broadly. They promise to continue to do so for the foreseeable future.