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.
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 , solar flares , coronal mass ejections , magnetospheric substorms  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 . This led to the classic analysis of Rosenbluth & Longmire in  (see also ). 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 . Spruit developed the flux tube theory including examining the stability of horizontal magnetic field in an atmosphere . 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 , 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 . Early studies [37,38] on simple cylindrical cases showed that crossing the linear stability boundary resulted in bifurcation to nearby helical equilibria—see also . Such helical states have also been predicted  and observed  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 ). 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  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 . 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.
We consider a simple one-dimensional line tied magnetized atmosphere with magnetic field (B0), gravitational acceleration (g), pressure (p0) and density (ρ0) given by 2.1The equilibrium force balance is 2.2Thus, to fully specify the equilibrium, we must give two of the three functions of height, B0(x), p0(x) or ρ0(x) (and a boundary condition if ρ0 is one of the two). Note that we use magnetic field units where B2 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 2.3where ν is the small viscosity coefficient. The magnetic field obeys the familiar equation for frozen in field 2.4i.e. E=−v×B. Density obeys the continuity equation: 2.5The field lines are tied to walls at 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. Ex(z=0,x,y,t)=Ex(z=L,x,y,t)=Ey(z=0,x,y,t)=Ey(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)=p0(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  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)=T0(x0), where x0 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 x0 is the original (Lagrangian) height of the field line. The pressure is then obtained from p(x,y,z,t)=ρ(x,y,z,t)T0(x0)/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 where ky is the wavenumber in the y-direction—a proof is given in appendix A of . This proof is a simple extension for our chosen boundary conditions of the proof in reference  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 —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 , where ρ0′=dρ0/dx. The typical y and z wavenumbers are , and the mode is localized in x over a distance . The viscosity is treated as small . The Lagrangian displacement of the plasma, ξ is of the form: 3.1The perpendicular motion of field lines is predominately in the 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: 3.2and perpendicular incompressibility ∂ξy/∂y0+∂ξx/∂x0=0. The Cowley & Artun's  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, r0, by r=r0+ξ. Here, we have expressed the displacement in mixed Lagrangian–Eulerian variables—z is the current (Eulerian) z position of a plasma element and x0 and y0 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 ξ: 3.3where is the y0 average of the squared displacement, ξ2 and . The local linear drive Γ0 is given by 3.4where the first term is the stabilizing field line bending, the second term is sometimes called the Parker instability drive and the third term is the Rayleigh–Taylor instability drive . Note that we have assumed that , so that the terms in equation (3.4) cancel to dominant order. Using the equilibrium relation, equation (2.2), we can write in a form perhaps more familiar to some readers: 3.5This form demonstrates that instability is driven in a magnetized atmosphere by temperature decreasing upwards and magnetic field decreasing upwards (magnetic buoyancy). Note that in our case with thermal conduction along the field it is the temperature, not entropy gradient, that matters—see discussion in e.g. Balbus ). The coefficients C0,C2,C3,C4 are given by 3.6Note that these coefficients have a different normalization (and the sign of C2) to the corresponding coefficients in reference —but we have kept the same notation. From equation (3.3), we can derive an energy equation: 3.7where the integrals are taken over the whole plasma volume and dV =dx0 dy0L. The total energy is obviously of the form of kinetic energy and a nonlinear potential energy. Note that motion can only be driven by the terms that lower the potential energy. In this case, these are the linear drive in regions, where and the nonlinear drive, where . The inertial term (C0) arises from both the vertical acceleration and the acceleration along the tilted field lines. The stabilizing C2 linear term arises from the sideways, y, displacement bending the field lines. The C3 stabilizing quasi-linear term arises from the modifications of the mean pressure, density and magnetic field profiles by the perturbations (through ). The C4 term arises from two effects: the change of the Rayleigh–Taylor drive at the displaced height (the term) and the nonlinear Parker flows (the term). There is no nonlinear field line bending contribution to C4, because the field is convected with the flow to this order. The term in equation (3.3) guarantees incompressibility and drives a nonlinear downdraft outside rising fingers (or flux tubes) of plasma. The effect of the C4 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 , the system is just above local marginal stability and that everywhere else, it is locally stable Γ0<0. In the unstable region, we can approximate 3.8Thus, the unstable region has the width where . Because we expect the perturbation to be localized around the unstable region, we may replace the smoothly varying functions C0,…,C4 by constants given by the functions evaluated at . The dynamics close to marginal stability yields the generic form equation (3.3) for fine scale MHD instabilities (with complicated expressions for the constants C0,…,C4) see references [15,18]. The dynamical behaviour of equation (3.3) was examined in detail in references [13,14,16,17], see reference  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, , the motion is linear, and the most unstable eigenfunctions for a given ky are given by 3.9where the width and growth rate are given by 3.10The threshold for instability is . Above this threshold, the growth rate is positive (instability) for . For weak viscosity, , the peak growth rate is , and the wave number at the peak is . The eigenfunction and growth rate are consistent with the assumed scalings—i.e. , and , where the Alfvén frequency is , and the scale length is defined by . The linear modes, as expected, are localized in x0 and fine scale in y0. When the equilibrium is evolving through marginal stability at a slow rate, we can write with ΓAτE≫1. Then, the characteristic linear growth rate for the first few exponentiations is —much less than Alfvénic growth.
(b) Nonlinear behaviour
The nonlinear terms in equation (3.3) become important when for the explosive nonlinearity and for the quasi-linear nonlinearity. For C4>0, the explosive nonlinearity drives the motion in the upwards (x increasing) direction and stabilizes the downwards motion. We shall assume C4>0 for discussion here, because the C4<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 C3 and C4 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 C0=0.248, C2=−0.352, C3=0.044, C4=0.216 3.11and . Note that the system is weakly growing over a very narrow region (). We take ν=10−10 and initialize with the most unstable linear mode with ky=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 dt=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, Δy0, and the disturbed height, Δx0 by 3.12For the evolution in figures 2 and 3, we can fit the later stage nonlinear motion to 3.13with t0=470.5. In this asymptotic regime, the contribution of the linear terms (C1 and C2 terms) is negligible, and the explosive nonlinearity (the C4 term) drives the growth of kinetic energy against viscous dissipation and the stabilizing quasi-linear term (the C3 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 () drives the tip of the finger faster than the sides of the finger which are pushed down by the downdraft ()—thus the finger narrows in y (i.e. Δy0 decreases). The stabilizing C3 nonlinearity is also reduced by narrowing the finger in y, thereby reducing the y average motion (i.e. ) and by broadening the mode in x (i.e. increasing Δx0) 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 (C3 and C4) yields (Δx0)2/Δy0∼ξ 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  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 —field lines in this region are metastable, so that when knocked hard by the rising finger of plasma, they are destabilized. This mechanism of progressive destabilization was termed 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 . 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 . Thus, while the asymptotic regime is reached before the equations break down, the singularity itself is not. What then happens to the rising fingers of plasma? In section 4, we examine a scenario for the further evolution of exploding fingers of plasma.
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 Bin and the field outside Bout (figure 5). The field in the tube is hardly bent in the y-direction, . Thus, the equation for a field line in the tube is given (to leading order) by x=x(x0,y0,z,t) and y=y0, where x0 is the undisplaced height of the field line and y0 is the undisplaced y position of the field line. Clearly, x0 and y0 are constant on field lines. Because the y dependence is not used further in this section, we omit the y0 dependence of x in subsequent expressions, i.e. we write x=x(x0,z,t). The field lines are tied to the wall therefore x(x0,0,t)=x(x0,L,t)=x0. We can write 4.1where Bz is a function of x and x0 to be found and and are unit vectors in the z- and x-direction, respectively. The force (per unit volume) on the plasma is 4.2We shall refer to p+B2/2 as the total pressure and B⋅∇B as the curvature force. The forces across the narrow tube (in the -direction) are formally large, , and must balance to this order, i.e. 4.3Thus, fast waves propagate (at speed ) across the tube on a time δ1/Cfast and equalize the total pressure inside and outside the tube. Thus, on the slow evolution time: 4.4where ‘in’ refers to inside the tube and ‘out’ refers to just outside the tube (at the same x and z along the tube; figure 4). We take the field and pressure outside the tube to be unperturbed so that 4.5are known. This approximation is justified as follows. The correction to the external field (blue lines in figure 5) are of order By∼Bout(δ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 . Inserting these estimates in equation (4.3) would give a correction to the internal total pressure of order —small for highly elliptical flux tubes.
The temperature is constant along the flux tube and at z=0, L is equal to T0(x0). Thus, T(x,z,t)=T0(x0), where x0=x0(x,z,t) is obtained from inverting x=x(x0,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, 4.6where we have used the boundary condition that x=x0 at z=0,L. From equations (4.6), (4.5), (4.4) and using the equilibrium relation equation (2.2), we obtain 4.7Note that 4.8Therefore, the force in the tube in the direction of motion x: 4.9where we have used equations (4.4) and (4.5) and equilibrium force balance equations (2.2). Note that the force is simply 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 () in the second line of equation (4.9). This requires 4.10
Thus (with ), the ellipticity of the flux tube must exceed a critical value for the theory of this section to hold 4.11
This relation is obeyed for a finitely displaced flux tube ((x−x0)∼L) with the cross-sectional shape of the weakly nonlinear regime (the left-hand side of equation (4.11) is and the right-hand side is ). During detonation, δ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
Equilibrium position of the field lines satisfies Fx=0—a second-order nonlinear ordinary differential equation for x=x(x0,z,t) for each x0. We integrate Fx=0 (using equation (4.12)) with the boundary conditions: 4.13to obtain: 4.14The solution of this equation must satisfy the additional boundary condition that x(x0,z=L)=x0. Therefore, by symmetry, the peak of the field line is then at z=L/2 with x=xpeak and (∂x/∂z)x0=0. Thus, xpeak(q) must satisfy 4.15and q must satisfy the eigenvalue condition: 4.16In practice, it is simpler to solve equation (4.14) numerically by shooting and iterating q until we find (∂x/∂z)x0=0 at z=L/2. In section 4a, 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 dA is proportional to , where the dl 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 4.17where the z integration is at fixed x0. Varying x(z,x0,t) in keeping x0 constant we get 4.18where we have integrated by parts and used (∂Bin/∂z)x0=(∂x/∂z)x0(∂Bin/∂x)x0. Thus, equilibrium, Fx=0, is a stationary state of the energy functional. To model the dynamics, we take a simple drag to balance the force, i.e. 4.19This models much more complex viscous and aerodynamic drag on the moving flux tube. The form of the drag does not, of course, affect the final erupted state and the explosive nature of the eruption. Motion of the tube with drag reduces the energy (monotonically) because 4.20Thus, evolution takes the flux tube to an energy minimum. These are linearly stable positions of the flux tube—not necessarily a global minimum of the energy. We emphasize that in our approximation the dynamics of each field line is independent. We can therefore consider each x0 separately. The generalization of this section to flux tube motion in a general three-dimensional magnetic equilibrium is given in appendix B of . 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  see also reference .
(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(x0,z,t)=x0+Δx with Δx much smaller than the typical scale height. Expanding equation (4.19) to second order in Δx, we obtain 4.21where we have expanded the right-hand side of equation (4.8). Close to marginal stability, the first two terms on the right-hand side of equation (4.21) almost cancel for . Substituting this into equation (4.21) multiplying by and integrating over z from 0 to L, we get 4.22where is given in equation (3.4) and C4 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 C4>0, because the case C4<0 is reproduced by changing the sign of ξ. Let ξ0 be the initial displacement, and the linear growth rate. Then 4.23In the linearly unstable region, γ>0 displacements grow from infinitesimal amplitudes. For ξ0>0, the growth accelerates (explosively) to infinity in a finite time, ; for ξ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 .
(b) Flux expulsion
In situations with substantial plasma beta (), the field lines can erupt until one point or part of the flux tube has zero magnetic field—(i.e. ). The flux tube has then expanded to an infinite cross section to conserve the flux. Clearly, this breaks the assumptions of a thin isolated flux tube. Indeed, our theory must be restricted to cases with Bin≫(δ1/L)B0. But we can still apply our theory when the field line approaches the B2in=0 asymptotic limit closely if we consider . Here, we examine what happens in this limit.
Suppose we have initially displaced the field line, so that xpeak=xzero, where Bin(xzero,x0)=0 (figure 6). Then 4.24But dl/dx≥1, where the equality holds when the field line is vertical and dl=dx. Thus, we can lower by making the field line have three segments: vertical from z=0, x=x0 to z=0, x=xzero, horizontal with x=xpeak=xzero for 0≤z≤L and vertical from z=L, x=xzero to z=L, x=x0. The horizontal segment makes no contribution to the energy, because Bin(xpeak,x0)=0. We illustrate this minimum energy state with the continuous line in figure 6—note that the sharp corner in the field happens where Bin=0 and therefore has no field line bending force. The minimum energy state has 4.25This is clearly a local minimum in energy—not necessarily a global minimum. However, it is clear that if the field line peak reaches the point of zero field strength (the dotted line in figure 6) the motion will continue towards the minimum energy state. The dynamics of the section of the field line with zero, or at least very small, field strength is a topic for future work. It is clear, however, that the zero field section of the tube remains buoyant and will continue to rise.
(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 4.26where , and are constants and xρ and xB are the heights of the maximum density and minimum field, respectively. From equation (2.2), we find , where is a constant. We take large pressure to focus on the magnetized Rayleigh–Taylor instability (rather than on the Parker instability), so that mgLρ/T0(x0)=gρ(x0)Lρ/p0(x0)≪1. We define normalized variables as 4.27and the constants 4.28The normalized relative change in energy is . The weakly nonlinear behaviour (from equation (4.22)) in normalized variables is 4.29where , and . The nonlinear C4 term drives upward motion for and and downwards motion for .
(i) Numerical solutions for model
The model equilibrium is specified by five normalized parameters (ignoring in the limit of large pressure): , , A, and . We investigate two cases in which we fix , , and and vary A—the aspect ratio of the box. In figure 7, we plot the growth rate and the explosive nonlinearity, , for the numerical cases—these do not change as 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 . The first case (case 1) has 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 and ) erupt into a flux expelled state. While these cases are representative, we have made no attempt yet to survey all the possible cases. For example, we have not considered cases where the field lines erupt downwards—there could also be cases where some of the lines erupt downwards and some upwards.
Case 1. A=0.161604. In this case, and and the magnetic field energy is larger than the energy changes, so that for all and . In figure 8, we show the evolution of (i.e. the field line with ) from initial conditions: (i) just above the critical amplitude for nonlinear instability and (ii) just below the critical amplitude for instability. By the time 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 () note that the exploding field line reaches saturation at a value of . We show the relative energy change in figure 10. We also solved equation (4.14) for the equilibrium position of the field line by the shooting method (varying q until we find a solution that vanishes at z=0 and z=L). There are three equilibrium positions: the trivial ; the saturated state with and the critical unstable equilibrium . The saturated state is, of course, identical with the final state of (i) in figure 8. The minimum relative energy needed to excite the explosive behaviour is the energy of the critical equilibrium—for the field line this is . This is far less than the energy released going to the saturated state .
Using the shooting method, we have computed the saturated states and critical equilibrium states for the field lines from to (in steps of ). Below x0∼0.2 and above x0∼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. and in figure 13, we plot the relative energy of the saturated state and the critical energy to excite the nonlinear explosive motion.
Case 2. A=0.1695. In this case, and . The minimum value of (or equivalently ) for a given is at . Field lines between x0=0.295 and x0=0.62 have a minimum value of that is less than zero. Thus, this case has a region of flux tubes where flux expulsion takes place as treated in §4b. All the field lines that would have a minimum less than zero (i.e. field lines with 0.295<x0<0.62) instead minimize their energy by taking the limiting rectangular shape of §4b. Thus, for these field lines, we evaluate using (the appropriately normalized) equation (4.25) for all others, we use equation (4.17). The critical and saturated energy are plotted in figure 14.
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 and where ϵ∼γ/Γ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 y0 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 in the 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∼ξdownD. We can estimate the energy needed to drive the downward motion as a fraction A1<1 of the field line bending energy, i.e. A1(ξdownB0/L)2Dδ2L=A1(ξupB0/L)2(δ1/D)2Dδ2L. The energy available from the upward going flux tube is some fraction A2<1 of the gravitational energy, i.e. A2(ρ0gξup)(δ1δ2L). The sideways motion of the field lines next to the flux tube gives another stabilizing energy (δ1B0/L)2(δ1δ2L) 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 5.1
Note the factor is of order one for a profile near marginal stability. Therefore, we can release energy if we make 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 ∼ρ0gLρV0, where V0 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.
This work contains no experimental data. To obtain further information on the models underlying this paper please contact.
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).
We declare we have no competing interests.
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.
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.