## Abstract

The eruption of multiple flux tubes in a magnetized plasma is proposed as a mechanism for explosive release of energy in plasmas. A significant fraction of the linearly stable isolated flux tubes are shown to be metastable in a box model magnetized atmosphere in which ends of the field lines are embedded in conducting walls. The energy released by destabilizing such field lines can be a large proportion of the gravitational energy stored in the system. This energy can be released in a fast dynamical time.

## 1. Introduction

The explosive release of energy from magnetically confined plasmas produces some of the most dramatic and destructive natural phenomena. In such events, a slowly evolving plasma suddenly erupts releasing a significant fraction of its stored magnetic, gravitational or pressure energy in a few tens of dynamical times (which is typically the Alfvén time or the free fall time). The stored energy is converted into some combination of heat, energetic particles, fast plasma flows and/or radiation. Tokamak disruptions [1], solar flares [2], coronal mass ejections [3], magnetospheric substorms [4] and edge localized modes in tokamaks [5–12] all exhibit this type of explosive behaviour. While there are a number of theories of these phenomena many of the central questions are still without quantitative answers. For example, What triggers the instability? What sets the timescale? How much energy is released? How much energy is given to energetic particles? Are there universal mechanisms? Here, we propose a specific explosive mechanism. Earlier weakly nonlinear analyses [13–20] demonstrated that near to the linear stability threshold *all* fine scale pressure and gravitational instabilities obey a generic (i.e. universal) equation that yields explosive dynamics. However, the subsequent evolution of these instabilities was not understood. In this paper, we argue that the fully nonlinear consequence of these instabilities is the eruption of multiple elliptical flux tubes on Alfvén timescales. As a specific example, we analyse the gravitational stability and fully nonlinear flux tube dynamics of a simple slab atmosphere that is line tied (i.e. one in which the ends of the field lines are embedded in conducting walls and therefore immovable). We suspect that the flux tube picture that emerges from our simple example captures some generic elements of the fully nonlinear evolution of *all* fine scale gravitational and pressure driven instabilities (e.g. ballooning modes [21–24]).

The dynamics of flux tubes in plasmas is an old subject—we cannot do justice to it all here. The stability of plasmas to the interchange of flux tubes was, we believe, first discussed by Edward Teller in 1954 in a classified meeting at Princeton on magnetically confined plasmas for fusion—see discussion in chapter 9 of [25]. This led to the classic analysis of Rosenbluth & Longmire in [26] (see also [27]). Parker made extensive use of a circular flux tube approximation in discussing dynamics in the solar convection zone; work that is reviewed in his famous monograph *cosmical magnetic fields* [28]. Spruit developed the flux tube theory including examining the stability of horizontal magnetic field in an atmosphere [29]. In fusion research, the theory of ‘blobs’ (isolated field aligned plasma structures that are similar to flux tubes but not necessarily frozen to the field) has been extensively researched (see [30,31] and references therein). Fan summarizes the work on gravitational stability and flux tubes for solar convection in reference [32], see also [33–35]. Much of this research on flux tubes concerns tubes in an otherwise unmagnetized plasma where circular tubes might be expected. We shall argue that tubes passing through a magnetized plasma are expected to be highly elliptical to minimize the stabilizing sideways perturbations of the surrounding field.

Just before an eruption, a magnetized plasma must be close to a linear instability boundary, because otherwise a large perturbation is needed to trigger instability. Ideal magnetohydrodynamic (MHD) instabilities just above the marginal stability point are of two distinct types: either *global* or *local* instabilities. The *kink mode* driven by current in a plasma cylinder is the archetype of the *global* (finite scale) MHD instability. At a critical current, the plasma becomes unstable to a single helical kink mode [36]. Early studies [37,38] on simple cylindrical cases showed that crossing the linear stability boundary resulted in bifurcation to nearby helical equilibria—see also [39]. Such helical states have also been predicted [40] and observed [41] to be the result of crossing a global resistive instability boundary. The dynamics associated with crossing a global instability boundary is not considered further in this paper. Just above, the linear instability boundary gravitational and pressure-driven ideal MHD instabilities are *local* instabilities (i.e. the eigenfunctions are localized in a narrow region). They also have a fine scale perpendicular to the field lines and are elongated along the field (e.g. the ideal ballooning mode [23]). *Local* modes driven by gravity become linearly unstable when the simple line tied equilibria addressed in this paper passes through the marginal linear instability threshold. In this transition, an infinity of local modes (all with infinitesimal perpendicular scale) becomes unstable as the profiles evolve through marginal stability, and the resultant threshold dynamics is complex.

A previous series of papers [13–19] examined the weakly nonlinear evolution of gravitational and pressure-driven instabilities. Analysis showed that these instabilities exhibit generic explosive dynamics when passing slowly through the linear instability boundary. Specifically, narrowing fingers of plasma accelerate and push aside the surrounding magnetic field to release energy. We review this dynamics in §3. However, it is not clear from the weakly nonlinear analysis what happens to these fingers in the deeply nonlinear regime. Simulations have confirmed the nonlinear instability [42] and the formation of narrow fingers [17,42–44] but lose resolution before the instability reaches saturation. Here, we propose that the fingers evolve into eruptions of elliptical flux tubes.

In §4, we derive equations of motion for an isolated elliptical tube in the simple one-dimensional line tied equilibrium (with gravity). The tubes have dynamics that is analogous to a simple zero-dimensional *transcritical* bifurcation dynamics—see [45]. Specifically, they have nonlinear instability drive (with a quadratic amplitude dependence of the force), so that even linearly stable tubes can be metastable and will erupt given sufficient perturbation. Then, triggering metastable flux tubes by either a finite perturbation or by linear instability can yield an explosive release of energy on an Alfvén timescale. We show in §4 that flux tubes have one of three fates depending on the tube, the profiles and the plasma *β*. Some flux tubes have no lower energy states and will not erupt. Some lower *β* tubes erupt to a stable lower energy equilibrium with finite displacement. Finally, some higher *β* flux tubes erupt to a singular state with zero magnetic field (flux expulsion) along part of the tube, see §4b. The energy released by each tube in the evolution to these states can be a finite fraction of the energy stored in the tube.

In §4c, we solve the equations for flux tube motion for simple model profiles. Two specific cases are presented in §4c(i): case 1 illustrates a case with no flux expelled erupted states, and case 2 has a region of erupted flux tubes that are in a flux-expelled state. We choose parameters, so we change the field strength (and therefore the plasma *β*) keeping the density profile and *B*/*L* constant (where *L* is the length of the field lines). Both cases then have the same growth rate and quadratic nonlinearity profile and both are just above marginal stability in a very small section of the profile. The saturated states of the flux tubes in each case are, however, different. In case 1, the saturated state of those field lines with a lower energy equilibrium than the initial state are shown in figure 11. Flux tubes from outside the linearly unstable region in figure 11 are metastable and must be displaced finitely to be destabilized. The energy released by metastable flux tubes and the critical energy to excite the flux tube is shown for case 1 in figure 13. Flux tubes originating from lower down overtake the upper lines, and the maximum energy release is associated with metastable flux tubes (not linearly unstable flux tubes). The field strength in case 2 is about 5% less than in case 1. In this case, some flux tubes minimize their energy in a saturated displaced equilibrium and some evolve to the singular flux expelled state. The energy released by metastable flux tubes and the critical energy to excite the flux tube is shown in figure 14. The analysis of the one-dimensional dynamics of isolated elliptical flux tubes as developed in §4 misses the crucial issue of how the tubes interact and the fraction of the metastable tubes that erupt. If only the linearly unstable flux tubes (in the narrow localized unstable region) are destabilized, then clearly only a small energy release is possible. However, in §5, we estimate that a significant fraction of all the metastable tubes can be destabilized, so that a considerable fraction of the stored energy can be released in a few Alfvén times by this eruptive scenario (see §5). Triggering happens when the system gets close enough to marginal linear stability that noise is sufficiently large to destabilize barely stable field lines.

## 2. Equilibrium and equations

Here, we outline the one-dimensional model geometry and the equations to be solved.

### (a) Equilibrium

We consider a simple one-dimensional line tied magnetized atmosphere with magnetic field (*B*_{0}), gravitational acceleration (**g**), pressure (*p*_{0}) and density (*ρ*_{0}) given by
*B*_{0}(*x*), *p*_{0}(*x*) or *ρ*_{0}(*x*) (and a boundary condition if *ρ*_{0} is one of the two). Note that we use magnetic field units where *B*^{2} has units of energy density.

### (b) Perturbation: equations and boundary conditions

We adopt a simple ideal MHD system with scalar viscosity that captures the essential dynamics of fast explosive motion driven by the gravitational potential energy. The ideal motion drives the system to small scales thus dissipative processes are inevitably important. These processes are often not well modelled by a scalar viscosity—nonetheless, the basic picture of instability and metastability is, we believe, a robust generic feature of magnetically confined systems. The equation of motion is
*ν* is the small viscosity coefficient. The magnetic field obeys the familiar equation for frozen in field
**E**=−**v**×**B**. Density obeys the continuity equation:
*z*=0 and *z*=*L*. At the boundaries *z*=0 and *z*=*L*, the *x* and *y* components of the electric field are set to zero—i.e. *E*_{x}(*z*=0,*x*,*y*,*t*)=*E*_{x}(*z*=*L*,*x*,*y*,*t*)=*E*_{y}(*z*=0,*x*,*y*,*t*)=*E*_{y}(*z*=*L*,*x*,*y*,*t*)=0. We also set the pressure and density to be unperturbed at the boundaries—i.e. *p*(*z*=0,*x*,*y*,*t*)=*p*(*z*=*L*,*x*,*y*,*t*)=*p*_{0}(*x*) and *ρ*(*z*=0,*x*,*y*,*t*)=*ρ*(*z*=*L*,*x*,*y*,*t*)=*ρ*_{0}(*x*); thus, motion along the field through the boundary is allowed (i.e. **v**⋅**B**≠0 at *z*=0,*L*). These boundary conditions differ from those in reference [13] where all flow at the boundary is zero. The boundary conditions adopted in this paper are closely related to the behaviour on infinite field lines where the perturbed field line displacement, pressure and density must vanish as distance along the field line goes to infinity.

We will take the thermal conduction along the field line to be fast, so that the temperature is constant on the field line, i.e. *T*(*x*,*y*,*z*,*t*)=*T*_{0}(*x*_{0}), where *x*_{0} is the height of the field line at the walls (at *z*=0,*L*). Note that because the field is frozen into the stationary wall, and the moving plasma *x*_{0} is the original (Lagrangian) height of the field line. The pressure is then obtained from *p*(*x*,*y*,*z*,*t*)=*ρ*(*x*,*y*,*z*,*t*)*T*_{0}(*x*_{0})/*m*, where *m* is the ion mass. This is equivalent to the usual adiabatic equation for pressure with the ratio of specific heats equal to one (i.e. *p* obeys the same equation as *ρ*, equation (2.5)).

The slab system with the line tied boundary conditions (and no viscosity) will become linearly unstable above a critical length—see below. The passage through marginal stability may be effected by either lengthening the box or evolution of the (density, pressure or magnetic field) profiles by diffusion, heating or cooling. The first modes to become unstable have *k*_{y} is the wavenumber in the *y*-direction—a proof is given in appendix A of [46]. This proof is a simple extension for our chosen boundary conditions of the proof in reference [47] which used elements of references [48,49].

## 3. Small amplitude nonlinear motion

Here, we calculate the early, small amplitude, nonlinear evolution of fine *y* scale perturbations when the system is just above the linear marginal stability threshold. The treatment closely follows the development of Cowley & Artun in 1997 [13]—thus we omit considerable detail. We measure the distance above marginal stability by a dummy large parameter *n* where the growth rate of the most unstable perturbation is order *ρ*_{0}′=d*ρ*_{0}/d*x*. The typical *y* and *z* wavenumbers are *x* over a distance ** ξ** is of the form:

*x*-direction and the structure is elongated in the

*x*-direction compared with the

*y*-direction. The

*y*motion is small. This structure maximizes the motion in the direction of gravity thereby enhancing the release of potential energy. Expansion of the MHD equations in powers of

*n*

^{−1/2}first yields the form of displacement:

*ξ*

_{y}/∂

*y*

_{0}+∂

*ξ*

_{x}/∂

*x*

_{0}=0. The Cowley & Artun's [13] paper was formulated entirely in Lagrangian variables where the current position of a plasma element,

**r**, is related to the initial position of the same element,

**r**

_{0}, by

**r**=

**r**

_{0}+

**. Here, we have expressed the displacement in mixed Lagrangian–Eulerian variables—**

*ξ**z*is the current (Eulerian)

*z*position of a plasma element and

*x*

_{0}and

*y*

_{0}the initial (Lagrangian)

*x*- and

*y*-positions of the plasma element. This rather complicated representation is convenient because it makes the lowest-order displacement satisfy exactly the boundary conditions at

*z*=0,

*L*. In higher order, we obtain the evolution equation for

*ξ*:

*y*

_{0}average of the squared displacement,

*ξ*

^{2}and

*Γ*

_{0}is given by

*field line bending*, the second term is sometimes called the

*Parker instability drive*and the third term is the

*Rayleigh–Taylor instability drive*[47]. Note that we have assumed that

*C*

_{0},

*C*

_{2},

*C*

_{3},

*C*

_{4}are given by

*C*

_{2}) to the corresponding coefficients in reference [13]—but we have kept the same notation. From equation (3.3), we can derive an energy equation:

*V*=d

*x*

_{0}d

*y*

_{0}

*L*. The total energy

*C*

_{0}) arises from both the vertical acceleration and the acceleration along the tilted field lines. The stabilizing

*C*

_{2}linear term arises from the sideways,

*y*, displacement bending the field lines. The

*C*

_{3}stabilizing

*quasi-linear*term arises from the modifications of the mean pressure, density and magnetic field profiles by the perturbations (through

*C*

_{4}term arises from two effects: the change of the Rayleigh–Taylor drive at the displaced height (the

*C*

_{4}, because the field is convected with the flow to this order. The

*downdraft*outside rising fingers (or flux tubes) of plasma. The effect of the

*C*

_{4}explosive nonlinearity on dynamics is examined in the flux tube limit below. We presume that the system evolves through marginal stability—perhaps by a slow lengthening of the box

*L*(

*t*). Then, in some narrow region around the height

*Γ*

_{0}<0. In the unstable region, we can approximate

*C*

_{0},…,

*C*

_{4}by constants given by the functions evaluated at

*C*

_{0},…,

*C*

_{4}) see references [15,18]. The dynamical behaviour of equation (3.3) was examined in detail in references [13,14,16,17], see reference [17] for a detailed comparison of the asymptotic theory with direct simulations of MHD—we provide only a brief summary.

### (a) Linear eigenfunctions

For very small amplitude, *k*_{y} are given by
*x*_{0} and fine scale in *y*_{0}. When the equilibrium is evolving through marginal stability at a slow rate, we can write *Γ*_{A}*τ*_{E}≫1. Then, the characteristic linear growth rate for the first few exponentiations is

### (b) Nonlinear behaviour

The nonlinear terms in equation (3.3) become important when *C*_{4}>0, the explosive nonlinearity drives the motion in the upwards (*x* increasing) direction and stabilizes the downwards motion. We shall assume *C*_{4}>0 for discussion here, because the *C*_{4}<0 case can be obtained by reversing the sign of *ξ*. Numerical solution of equation (3.3) ([13,14,16–18]) reveals a generic scenario for the small amplitude nonlinear evolution. First, the instability grows linearly in the unstable region. Once the amplitude is large enough that the nonlinear terms (the *C*_{3} and *C*_{4} terms) begin to dominate, the upward moving plasma accelerates and narrows in *y*. The dynamics becomes characterized by fingers of explosively growing plasma of width Δ*y* (figure 1). The evolution approaches a finite time singularity—see references [13,14,16–18].

Here, we present an example numerical solution of equation (3.3) with coefficients derived from the model of §4c(i). This solution is calculated in the normalized tilde variables of §4c(i) where the coefficients (which are the same for both cases 1 and 2) are *C*_{0}=0.248, *C*_{2}=−0.352, *C*_{3}=0.044, *C*_{4}=0.216
*ν*=10^{−10} and initialize with the most unstable linear mode with *k*_{y}=5412—i.e. a wavelength of λ_{y}=0.00116. In figure 2, we plot the time behaviour of the energy terms of equation (3.7) at early and late times. The timesteps are d*t*=0.01. The spatial step size in the *x*-direction is around 6.84*10^{−4} (there are 200 points in a range of 0.136) and in the *y*-direction, the step size is 5.86*10^{−6} (there are 200 points in a wavelength). The simulation is terminated at a time of *t*=418 when the resolution becomes too poor.

We define the ‘finger’ width in *y*, Δ*y*_{0}, and the disturbed height, Δ*x*_{0} by
*t*_{0}=470.5. In this asymptotic regime, the contribution of the linear terms (*C*_{1} and *C*_{2} terms) is negligible, and the explosive nonlinearity (the *C*_{4} term) drives the growth of kinetic energy against viscous dissipation and the stabilizing quasi-linear term (the *C*_{3} term). In figure 2, it is clear that the nonlinear terms largely balance with the remaining drive providing the viscous heating and growth of kinetic energy. The explosive nonlinearity (*downdraft* (*y* (i.e. Δ*y*_{0} decreases). The stabilizing *C*_{3} nonlinearity is also reduced by narrowing the finger in *y*, thereby reducing the *y* average motion (i.e. *x* (i.e. increasing Δ*x*_{0}) to minimize the flattening of the mean profiles. The viscosity opposes the narrowing in *y* but does not halt it. We have not been able to find a simple analytic derivation of the scaling in equation (3.13)—i.e. with viscosity. However, the balance of the two nonlinear terms (*C*_{3} and *C*_{4}) yields (Δ*x*_{0})^{2}/Δ*y*_{0}∼*ξ* which is obeyed by the scaling in equation (3.13). In [16,17], the transition to explosive growth with finite larmor radius terms (rather than viscosity) opposing the narrowing in *y* was examined. Analytic expressions for the exponents of explosive growth were derived in this case. These are, of course, different from those with viscous growth given in equation (3.13). Fong also showed [17] that the dynamics of the asymptotic theory agreed with his full MHD simulations until resolution was lost. We note that the (quasi-linear) nonlinearity broadens the mode into the linearly stable region *detonation* [13,14,16–18] because of the analogy with chemical explosives. Linear instability is not necessary for explosive growth, because finite perturbations in equation (3.3) will destabilize linearly stable profiles and trigger detonation. The quasi-linear nonlinearity suppresses all but the largest amplitude fingers [17]. Thus, by the end of this small amplitude evolution, the dynamics consists of a few rapidly rising fingers of plasma—the flux in these fingers form elliptical flux tubes.

The treatment in this section cannot capture displacements as large as the width of the unstable region, because

## 4. Flux tube dynamics

Here, we examine the finite amplitude dynamics of single isolated narrow line tied flux tubes in our box equilibrium. The tubes have elliptical cross sections, elongated in the direction of motion (*x*) and narrower across (Δ*y*=*δ*_{1}≪Δ*x*=*δ*_{2}≪*L*) (figure 4). The exact cross-sectional shape of the tube is not important here—just that it is narrow and considerably elongated in the direction of motion (figure 4). This shape allows the tube to ‘knife’ through the plasma separating the surrounding field lines very little—indeed, we shall take *δ*_{1} to be sufficiently small that to lowest order the surrounding field lines are effectively unperturbed. We conjecture that such tubes are the later stage evolution of the fingers seen in the early-stage nonlinear development described above. The rising fingers in the small amplitude motion are, however, never independent (isolated). Thus, we must assume that as they evolve from fingers to a moving flux tube they become independent. The isolated tubes are presumed to move somewhat slower than the sound speed, because we are interested in the partially viscous behaviour near marginal stability and the saturated states of the flux tube.

Consider the field-aligned flux tube that is displaced through the plasma. The field inside the tube is denoted **B**_{in} and the field outside **B**_{out} (figure 5). The field in the tube is hardly bent in the *y*-direction, *x*=*x*(*x*_{0},*y*_{0},*z*,*t*) and *y*=*y*_{0}, where *x*_{0} is the undisplaced height of the field line and *y*_{0} is the undisplaced *y* position of the field line. Clearly, *x*_{0} and *y*_{0} are constant on field lines. Because the *y* dependence is not used further in this section, we omit the *y*_{0} dependence of *x* in subsequent expressions, i.e. we write *x*=*x*(*x*_{0},*z*,*t*). The field lines are tied to the wall therefore *x*(*x*_{0},0,*t*)=*x*(*x*_{0},*L*,*t*)=*x*_{0}. We can write
*B*_{z} is a function of *x* and *x*_{0} to be found and *z*- and *x*-direction, respectively. The force (per unit volume) on the plasma is
*p*+*B*^{2}/2 as the *total pressure* and **B**⋅**∇****B** as the *curvature force*. The forces across the narrow tube (in the *δ*_{1}/*C*_{fast} and equalize the total pressure inside and outside the tube. Thus, on the slow evolution time:
*x* and *z* along the tube; figure 4). We take the field and pressure outside the tube to be unperturbed so that
*B*_{y}∼*B*_{out}(*δ*_{1}/*δ*_{2}), where we have estimated that for a finitely displaced flux tube the external field line is bent a distance *δ*_{1} over a length of *δ*_{2}. The field line bending of the external field lines gives a sideways *y* force of order

The temperature is constant along the flux tube and at *z*=0, *L* is equal to *T*_{0}(*x*_{0}). Thus, *T*(*x*,*z*,*t*)=*T*_{0}(*x*_{0}), where *x*_{0}=*x*_{0}(*x*,*z*,*t*) is obtained from inverting *x*=*x*(*x*_{0},*z*,*t*). Because we are interested in stable displaced equilibrium states of the flux tube and slow eruptions from an unstable state, we set **F**⋅**B**=0. Thus,
*x*=*x*_{0} at *z*=0,*L*. From equations (4.6), (4.5), (4.4) and using the equilibrium relation equation (2.2), we obtain
*x*:
*the field line bending force and the buoyancy force from Archimedes' principle.* We have ignored the vertical forces from the corrections to the total pressure (

Thus (with

This relation is obeyed for a finitely displaced flux tube ((*x*−*x*_{0})∼*L*) with the cross-sectional shape of the weakly nonlinear regime (the left-hand side of equation (4.11) is *δ*_{2} increases, and *δ*_{1} decreases, so the inequality in equation (4.11) becomes even better satisfied. Thus, by the time the flux tube is finitely displaced, we expect the inequality to be well satisfied.^{1}

Because equation (4.9) only involves derivatives along the field line, we can consider dynamics of each field line separately. Using equations (4.7), (4.8) and (4.1) in equation (4.9), we obtain

Equilibrium position of the field lines satisfies *F*_{x}=0—a second-order nonlinear ordinary differential equation for *x*=*x*(*x*_{0},*z*,*t*) for each *x*_{0}. We integrate *F*_{x}=0 (using equation (4.12)) with the boundary conditions:
*x*(*x*_{0},*z*=*L*)=*x*_{0}. Therefore, by symmetry, the peak of the field line is then at *z*=*L*/2 with *x*=*x*_{peak} and (∂*x*/∂*z*)_{x0}=0. Thus, *x*_{peak}(*q*) must satisfy
*q* must satisfy the eigenvalue condition:
*q* until we find (∂*x*/∂*z*)_{x0}=0 at *z*=*L*/2. In section 4*a*, we use this procedure to find the equilibria for a simple model atmosphere.

The magnetic energy in a narrow flux tube with flux d*ψ*=*B* d*A* is proportional to *l* integration is a line integral along the flux tube. We now show that the equilibria are stationary points of the magnetic energy. Let us define the energy functional
*z* integration is at fixed *x*_{0}. Varying *x*(*z*,*x*_{0},*t*) in *x*_{0} constant we get
*B*_{in}/∂*z*)_{x0}=(∂*x*/∂*z*)_{x0}(∂*B*_{in}/∂*x*)_{x0}. Thus, equilibrium, *F*_{x}=0, is a stationary state of the energy functional. To model the dynamics, we take a simple drag to balance the force, i.e.
*x*_{0} separately. The generalization of this section to flux tube motion in a general three-dimensional magnetic equilibrium is given in appendix B of [46]. Also, the details of how the boundary conditions at the conducting walls are met by a boundary layer is contained in appendix C of reference [46] see also reference [17].

### (a) Small amplitude flux tube motion

Here, we examine the small amplitude motion of the flux tube to connect the flux tube theory to the more general small amplitude theory given above in §3. Let us assume *x*(*x*_{0},*z*,*t*)=*x*_{0}+Δ*x* with Δ*x* much smaller than the typical scale height. Expanding equation (4.19) to second order in Δ*x*, we obtain
*z* from 0 to *L*, we get
*C*_{4} in equation (3.6). Clearly, the small amplitude flux tube dynamics captures the linear and explosive nonlinearity parts of equation (3.3)—the assumption of thin isolated flux tubes orders out the remaining terms in the forces of (3.3) and drag replaces the inertial and viscous terms. Again, we will take *C*_{4}>0, because the case *C*_{4}<0 is reproduced by changing the sign of *ξ*. Let *ξ*_{0} be the initial displacement, *γ*>0 displacements grow from infinitesimal amplitudes. For *ξ*_{0}>0, the growth accelerates (explosively) to infinity in a finite time, *ξ*_{0}<0, the growth saturates at an amplitude *ξ*=−*ξ*_{c}. In the linearly stable region *γ*<0, the displacement decays unless *ξ*>−*ξ*_{c} when it grows explosively. Thus, −*ξ*_{c} (which is positive for *γ*<0) is the *critical amplitude* to excite a metastable field line. Clearly, the small amplitude dynamics is of the transcritical form [45].

### (b) Flux expulsion

In situations with substantial plasma beta (*B*_{in}≫(*δ*_{1}/*L*)*B*_{0}. But we can still apply our theory when the field line approaches the *B*^{2}_{in}=0 asymptotic limit closely if we consider

Suppose we have initially displaced the field line, so that *x*_{peak}=*x*_{zero}, where *B*_{in}(*x*_{zero},*x*_{0})=0 (figure 6). Then
*l*/d*x*≥1, where the equality holds when the field line is vertical and d*l*=d*x*. Thus, we can lower *z*=0, *x*=*x*_{0} to *z*=0, *x*=*x*_{zero}, horizontal with *x*=*x*_{peak}=*x*_{zero} for 0≤*z*≤*L* and vertical from *z*=*L*, *x*=*x*_{zero} to *z*=*L*, *x*=*x*_{0}. The horizontal segment makes no contribution to the energy, because *B*_{in}(*x*_{peak},*x*_{0})=0. We illustrate this minimum energy state with the continuous line in figure 6—note that the sharp corner in the field happens where *B*_{in}=0 and therefore has no *field line bending force*. The minimum energy state has

### (c) Flux tube motion in a model atmosphere

Here, we examine the flux tube motion in a simple model magnetized atmosphere. Let the unperturbed density and magnetic field be given by
*x*_{ρ} and *x*_{B} are the heights of the maximum density and minimum field, respectively. From equation (2.2), we find *mgL*_{ρ}/*T*_{0}(*x*_{0})=*gρ*(*x*_{0})*L*_{ρ}/*p*_{0}(*x*_{0})≪1. We define normalized variables as
*C*_{4} term drives upward motion for

#### (i) Numerical solutions for model

The model equilibrium is specified by five normalized parameters (ignoring *A*, *A*—the aspect ratio of the box. In figure 7, we plot the growth rate *A* is changed. Note how the system is just above marginal linear stability—i.e. the local growth rate is slightly positive in a narrow region around *A*=0.161604 and the second case (case 2) has *A*=0.1695. We choose these values, so that in case 1 all the field lines erupt without flux expulsion and in case 2 some of the field lines (those between

*Case* 1. *A*=0.161604. In this case, *t*=150, the field line has reached stable equilibria—the erupted saturated state for the initial conditions (i) and the unperturbed state for (ii). In figure 9, we show the peak amplitude of the field line (*q* until we find a solution that vanishes at *z*=0 and *z*=*L*). There are three equilibrium positions: the trivial

Using the shooting method, we have computed the saturated states and critical equilibrium states for the field lines from *x*_{0}∼0.2 and above *x*_{0}∼1.25, we found no equilibrium states other than the initial state. The saturated field line shapes are shown in figure 11. Note how the lower lines overtake the upper lines. In figure 12, we plot the saturated height and critical displacement as a function of the original field line position i.e.

*Case* 2. *A*=0.1695. In this case, *x*_{0}=0.295 and *x*_{0}=0.62 have a minimum value of *flux expulsion* takes place as treated in §4b. All the field lines that would have a minimum *x*_{0}<0.62) instead minimize their energy

## 5. Discussion and conclusion

The calculations presented in this paper support a model of explosive release of energy in magnetized plasmas by the destabilization of multiple metastable flux tubes. This eruption model is far from complete; indeed, a number of questions remain. Nonetheless, some results are clear. In §4, we demonstrated the metastability of some isolated thin elliptical flux tubes in a magnetized atmosphere. We show that such tubes can erupt on Alfvénic timescales when they cross the linear stability boundary, or when they are displaced by an amplitude greater than the critical amplitude. With viscous (or drag) dissipation, the flux tubes will relax to finitely displaced (saturated) equilibrium states releasing a significant fraction of their stored energy. The energy needed to destabilize all the metastable tubes is considerably less than the energy released (figures 13 and 14). In some high-pressure cases, the saturated equilibrium state is singular and the flux tube swells to infinite thickness thus reducing the field in the tube to zero—see §4b. We have also shown (see §3) that the weakly nonlinear behaviour near marginal stability yields growth in a narrow unstable region with erupting fingers pushing into and progressively destabilizing the metastable region—a process we have called *detonation*. More detail of this mechanism is given in references [16–19]. We have conjectured that these fingers evolve into flux tube eruptions. In the rest of this discussion section, we address unresolved issues for our eruption model qualitatively.

The size, number and shape of the flux tubes in an eruption must depend to some extent on the noise that creates the perturbation. We distinguish between slow evolution of the equilibrium and *noise* perturbations that temporarily move the system out of equilibrium. A large perturbation of a linearly stable plasma could trigger energy release but, at least in the early stage of eruption, the shape of the perturbation must determine the tubes that participate. However, noise levels are usually small in systems of interest. Thus, large perturbations are rare and would themselves require an explanation. In fusion experiments, the background *drift wave* turbulence provides a constant source of weak low-frequency noise with density perturbations of a few percentage at most. Triggering by other MHD instabilities is also seen. The noise in the solar corona comes from the convective motions that slowly perturb the foot points of the field lines. It seems likely, therefore, that eruptions begin in a region that is very close to being linearly unstable with perturbations (some part of the noise spectrum) that are close to the most unstable linear perturbations. The weakly nonlinear dynamics of marginally unstable atmospheres (see §3) shows that the dynamics evolves into a number of interacting explosively growing fingers. We argue that this is the beginning of the eruption of elliptical flux tubes—but not necessarily isolated elliptical flux tubes.

We have assumed that the erupting flux tubes are strongly elliptical in shape, i.e. *δ*_{1}≪*δ*_{2} in figure 5. The shape of the eigenmode in the linear regime has *ϵ*∼*γ*/*Γ*_{A}≪1 and *γ* is the growth rate, *Γ*_{A} the Alfven frequency and the typical equilibrium scale length is *L*_{γ} (see §3a) for definitions). The evolution in the weakly nonlinear regime depends on the dissipation—here the narrowing of the width of fingers in *y*_{0} depends on viscosity (in fusion finite larmor radius physics controls the narrowing, see [16,17]). If the resistivity is larger than the viscosity (in a small magnetic Prandtl number plasma), the eruption would be very different as the plasma would disconnect from the field. We will continue to assume that the resistive diffusion across the flux tube is slower than the eruption time. As a tube erupts, it will push the tubes in front of it causing some of them to become destabilized, it will also drag up tubes from below. This transfer of energy to metastable tubes is the *detonation* mechanism identified in the weakly nonlinear regime (see §3)—we do not have a good understanding of its efficiency in the fully nonlinear regime. If the detonation by each tube is efficient, then we expect the whole height of the metastable region to erupt—in §4c—this would lead to *δ*_{2}∼*L*_{ρ} (figure 11). The field lines surrounding the flux tube (coloured blue in figure 5) are bent sideways (in the *y*-direction). They will tend to flatten the tube with a force of order *y*-direction. This is however small compared with the forces in the *x*-direction and therefore, we expect it to make little difference to the final state.

If detonation was completely efficient, then we would expect the final state to be the lowest possible energy state. Such a state must pack in as many erupted tubes as possible. Let us now estimate the (*y*) distance between tubes (*D*) which maximizes the energy release. In figure 15, we illustrate two tubes spaced by *D* and displaced upwards by an average distance of *ξ*_{up}. Typically, *ξ*_{up} is of order *L*_{ρ} the gravitational scale height and as discussed above with efficient detonation, we expect *δ*_{2}∼*L*_{ρ}. The plasma between the tubes must be displaced down to make room for the flux tube rising. The plasma is compressible, but we will estimate the downwards displacement from rough incompressibility – *δ*_{1}*δ*_{2}∼*ξ*_{up}*δ*_{1}∼*ξ*_{down}*D*. We can estimate the energy needed to drive the downward motion as a fraction *A*_{1}<1 of the field line bending energy, i.e. *A*_{1}(*ξ*_{down}*B*_{0}/*L*)^{2}*Dδ*_{2}*L*=*A*_{1}(*ξ*_{up}*B*_{0}/*L*)^{2}(*δ*_{1}/*D*)^{2}*Dδ*_{2}*L*. The energy available from the upward going flux tube is some fraction *A*_{2}<1 of the gravitational energy, i.e. *A*_{2}(*ρ*_{0}*gξ*_{up})(*δ*_{1}*δ*_{2}*L*). The sideways motion of the field lines next to the flux tube gives another stabilizing energy (*δ*_{1}*B*_{0}/*L*)^{2}(*δ*_{1}*δ*_{2}*L*) which is always small for elliptical tubes. We release energy if the energy from the upwards moving flux tube exceeds the energy from the downwards motion, i.e. if

Note the factor *D* bigger than *δ*_{1} by a finite factor. The total energy release is then a finite fraction of the energy available if the magnetic field were absent—the gravitational energy ∼*ρ*_{0}*gL*_{ρ}*V*_{0}, where *V*_{0} is the volume of the metastable region.

Energy in the eruption is dissipated by viscosity or aerodynamic drag in our model. This would yield simple ion heating. If the motion does approach the speed of sound, we would expect the formation of shocks and possible acceleration of particles to non-thermal energies. The elliptical flux tubes in the saturated state have two current sheets, one on each side of the flux tube. These current sheets have opposite sign. Even if resistive diffusion is negligible during the eruption, it could act on a longer timescale on the saturated state reconnecting the field lines releasing more energy as heat. Such resistive diffusion would smooth out the magnetic field in the final state. We have ignored many secondary effects of the eruption; for example, we have not considered the possibility of secondary instabilities driven by the gradients across the flux tube or Kelvin–Helmholtz instabilities of the shear flow around the moving tube.

The discussion in this section is speculative. Clearly, high-resolution simulations that can follow the eruption to saturation would help resolve the many unanswered questions. Future work will pursue these questions and the extension of these ideas to more complex equilibria.

## Data accessibility

This work contains no experimental data. To obtain further information on the models underlying this paper please contact publicationsmanager{at}ccfe.ac.uk.

## Authors' contributions

S.C.C. acknowledges discussions on the issues presented in this paper with Omar Hurricane, Bryan Fong, Alex Schekochihin, Chris Ham, Jack Connor and Felix Parra. The mathematical analysis is due to S.C.C. and H.R.W. The paper was written by S.C.C. S.A.H. performed the computations in §3b and B.C. performed the computations in §4c(i).

## Competing interests

We declare we have no competing interests.

## Funding

This project has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement no. 633053 and from the RCUK Energy Programme (grant no. EP/I501045). S.A.H. is funded by the German Academic Exchange Service (DAAD-Stipendium für Doktoranden). H.R.W. is a Royal Society Wolfson Research Merit Award Holder. The views and opinions expressed herein do not necessarily reflect those of the European Commission.

## Footnotes

An invited Perspective to mark the election of the author to the fellowship of the Royal Society in 2014.

↵1 In the weakly nonlinear regime, the dominant vertical forces (those in equation (4.9)) cancel to the lowest order owing to the marginal stability. Therefore, the corrections to the force owing to the sideways (

*y*) motion*are*kept in equation (3.3).

- Received November 25, 2014.
- Accepted July 7, 2015.

- © 2015 The Author(s)

Published by the Royal Society. All rights reserved.