Contents

• The evolution of twisting coronal magnetic flux tubes
• 1 Introduction
• 2 Numerical model
• 3 Quasi-static and dynamic evolution of twisting flux tubes
• 4 Relaxation
• 5 Sigmoidal shape
• 6 Conclusions
• Appendix A: The modified two-step Lax-Wendroff method
• Appendix B: Treatment of the numerical parameters
• References

A&A 406, 1043-1059 (2003)
DOI: 10.1051/0004-6361:20030692

The evolution of twisting coronal magnetic flux tubes

T. Török^1,2 - B. Kliem²

1 - Astrophysical Institute Potsdam, 14482 Potsdam, Germany
2 - School of Mathematics and Statistics, University of St. Andrews, St Andrews, Fife KY16 9SS, UK

Received 23 December 2002 / Accepted 8 May 2003

Abstract
We simulate the twisting of an initially potential coronal flux tube by photospheric vortex motions, centred at two photospheric flux concentrations, using the compressible zero-beta ideal MHD equations. A twisted flux tube is formed, surrounded by much less twisted and sheared outer flux. Under the action of continuous slow driving, the flux tube starts to evolve quasi-statically along a sequence of force-free equilibria, which rise slowly with increasing twist and possess helical shape. The flux bundle that extends from the location of peak photospheric current density (slightly displaced from the vortex centre) shows a sigmoidal shape in agreement with observations of sigmoidal soft X-ray loops. There exists a critical twist, above which no equilibrium can be found in the simulation and the flux tube ascends rapidly. Then either stable equilibrium ceases to exist or the character of the sequence changes such that neighbouring stable equilibria rise by enormous amounts for only modest additions of twist. A comparison with the scalings of the rise of flux in axisymmetric geometry by Sturrock et al. (1995) suggests the former. Both cases would be observed as an eruption. The critical end-to-end twist, for a particular set of parameters describing the initial potential field, is found to lie in the range $2.5\pi<\Phi_{\rm c}<2.75\pi$. There are some indications for the growth of helical perturbations at supercritical twist. Depending on the radial profiles of the photospheric flux concentration and vortex velocity, the outer part or all of the twisted flux expands from the central field line of the flux tube. This effect is particularly efficient in the dynamic phase, provided the density is modeled realistically, falling off sufficiently rapidly with height. It is expected to lead to the formation of a cavity in which the twisted flux tube is embedded, analogous to the typical structure of coronal mass ejections.

Key words: MHD - Sun: activity - Sun: corona - Sun: coronal mass ejections (CMEs) - Sun: flares - Sun: magnetic fields

   
1 Introduction

The slow evolution of magnetic flux in the solar corona is governed by the condition of frozen-in plasma (Lundquist number $S\gg1$), by the strong variation of the plasma beta (the ratio of thermal pressure to magnetic pressure) from $\beta>1$ in the dense convection zone to $\beta\sim10^{-3}\dots10^{-2}$ in the inner corona, and by the related huge difference between the sound and Alfvén velocities in the convection zone on the one hand and the coronal Alfvén velocity on the other. Together these conditions imply that, to a good approximation, the coronal field evolves in a quasi-static manner through a series of nearly force-free equilibria while it is slowly driven by motions in a bottom layer in which the change of the magnitude of $\beta$ occurs; the latter is often simply taken to nearly coincide with the photosphere.

Departures from equilibrium at large spatial scales and the development of dynamics at the corresponding temporal scales characteristic of the corona occur sporadically in coronal mass ejections (CMEs), eruptive prominences, and flares. Many of these events lead to the opening of a closed magnetic configuration into interplanetary space. One way to model the occurrence of such eruptive events is to consider sequences of stable, magnetically closed, ideal MHD equilibria for continuously varying photospheric boundary conditions, searching for conditions that lead to a magnetically open equilibrium, an unstable equilibrium, or to complete absence of a neighbouring equilibrium at some point in the sequence. Dipolar or arcade-like potential magnetic fields are typically chosen as initial configurations, and shearing, twisting, or converging motions in the area of flux concentrations have been considered as driving photospheric motions (e.g., Sakurai 1979; for recent reviews see, e.g., Forbes 2000; Klimchuk 2001).

This paper presents a numerical study of loop-shaped flux tube evolution and stability within the framework of the ideal MHD equations in the limit $\beta\!=\!0$. The simulations are directed at the question whether a flux tube experiences a transition to eruptive behaviour if a sufficient amount of twist is accumulated by vortex motions at its photospheric footpoints. Observationally, this is motivated by two properties of CMEs (Hundhausen 1999). First, they appear to be ejected from parts of the atmosphere that are dominated by closed magnetic structures. Second, many CMEs include an eruptive prominence which often has the appearance of a strongly twisted, filamentary flux bundle. Theoretical investigations have also stimulated interest in single twisted flux tubes, since their eruption can be energetically favourable compared to the eruption of a whole sheared magnetic arcade. The eruption of a whole arcade requires the lifting of the whole overlying flux, which is energetically unfavourable (Sturrock 1991; Aly 1991) and requires an unrealistically large amount of shear (e.g., Amari et al. 1996a; Mikic & Linker 1994), or a special magnetic configuration, which enables the energetically favourable so-called magnetic breakout of the arcade by removing the overlying flux through magnetic reconnection (Antiochos et al. 1999). In the case of a single flux tube, part of the overlying flux may bend sidewards when twist is imposed at the footpoint area of the flux tube causing it to ascend (Amari et al. 1996b). The latter process requires less energy input than the lifting of the flux overlying an arcade of magnetic loops. It is not clear whether such loop twisting may lead to a real instability in the sense that magnetic energy is liberated; nevertheless, it is of interest as a building block for more sophisticated models of solar eruptions.

In our simulations twisting motions are applied at the photospheric footpoints of an initially potential coronal magnetic flux bundle. These motions serve both to create a twisted flux tube and to push it into the region of possible instability. Observational evidence for twisting photospheric motions can be found, e.g., in Kempf (1910), Hofmann et al. (1987), Brandt et al. (1988), Nightingale et al. (2002). Such observations are not very numerous, and a direct relationship between twisting photospheric motions and eruptive events is difficult to establish observationally. Part of this problem is clearly due to the difficulty in detecting such motions in observations of the photosphere, part is due to the fact that magnetic flux may generally emerge in an already twisted state (e.g., Lites et al. 1995; Leka et al. 1996) so that only moderate additional twist is required after emergence to destabilize it. There is also the possibility that eruptions arise from the creation of unstable twisted flux by magnetic reconnection (Pevtsov et al. 1996; van Driel-Gesztelyi et al. 2000; Moore et al. 2001) rather than from continuous photospheric twisting. Again, we regard our model as a building block and reference system for more sophisticated models of solar eruptions.

Continuous twisting of a flux tube of finite length inevitably leads to kink instability. The threshold for instability depends on the ratio of the azimuthal and axial field components of the flux tube, $B_\varphi/B_z$, the length-to-width ratio of the tube, L/r, and on the radial profile of the flux tube. It is usually expressed in terms of the flux tube twist

$$\Phi = \frac{LB_\varphi(r_0)}{r_0B_z(r_0)},$$ (1)

where r₀ is a characteristic radius, typically one at which the current density has fallen significantly compared with its peak value within the flux tube, or in terms of a radially averaged flux tube twist. For a straight, cylindrically symmetric flux tube with fixed (line-tied) ends and uniform twist (the Gold-Hoyle force-free equilibrium), the threshold was numerically found to be $\Phi_{\rm GH}\approx2.49\pi$(Einaudi & van Hoven 1983; Hood & Priest 1981). For force-free equilibria having other radial profiles, similar values for the critical average twist are found if the axial field component does not change sign (see, e.g., reviews by Einaudi 1990; Hood 1992 and references therein). Other values of the critical average twist are found for equilibria that include pressure gradients or a reversal of B_z at some radius, but these cases are not relevant to our study. For loop-shaped flux tubes with fixed ends, a threshold value is not known. It is not even clear whether an instability exists at all or whether the flux tube avoids the instability by expansion (in length and width) to ever larger equilibrium configurations as more twist is applied at its footpoints. There are, however, several indications that a critical amount of twist, above which a qualitative change of behaviour occurs, exists in this system as well (Amari et al. 1996b; Klimchuk et al. 2000). The nature of the change of behaviour has not yet been systematically investigated. Either a true instability is then reached, or successive equilibrium configurations have very rapidly increasing loop size for small additions of twist above this value. It is difficult to distinguish between these cases with finite-size simulation boxes. Also, on the Sun both cases would appear as an eruption of the flux tube if the apex height of successive equilibria quickly reaches the solar wind acceleration region.

Previous treatments of slow (quasi-static) twisting of a loop-shaped flux tube may be categorized according to whether stable force-free equilibria were calculated (Van Hoven et al. 1995; Klimchuk et al. 2000) or whether the twisting was enforced well into a region where the system became dynamic (Tokman & Bellan 2002; Amari et al. 1996b), which was then followed. In all investigations the twisted flux tube was found to inflate, to develop an S shape (though not necessarily of the sense suggested by the observations), and to remain stable for twists $\Phi\le2\pi$. Van Hoven et al. (1995) obtained a stable equilibrium for an applied end-to-end twist of $2.82\pi$, a value which presumably reflects the stabilizing influence of an overlying current-free magnetic arcade. Klimchuk et al. (2000) calculated stable equilibria for $\Phi=\pi$and $2\pi$, noting that "the relaxation becomes prohibitively slow for twists much larger than $2\pi$''. This might indicate a qualitative change of behaviour above a critical twist. In view of the typical size relations between coronal magnetic loops and sunspots, they applied the twist only to a sub-area of the photospheric flux concentration, displaced from the peak, which presumably leads to a more stable flux tube than the twisting of the whole centre of the flux concentration that was used in all other papers. Amari et al. (1996b) followed the evolution of a continuously and relatively slowly twisted flux tube using a reduced set of zero-beta ideal MHD equations. The density was determined from the prescribed condition of uniform Alfvén velocity, implying a substantial decrease of the density with height as is characteristic of the corona. These authors found a strong indication for eruptive behaviour above a critical twist, beyond which the quasi-static evolution of the system ceased and an accelerated rise of the flux tube commenced with up to near-Alfvénic velocities. Also, relaxed states could not be found for twists exceeding the critical value which was on the order of $2.5\pi$ (Amari, personal communication 2002). An analogous system was studied by Tokman & Bellan (2002). They included resistivity but prescribed uniform density. Dynamical behaviour was also found but only after magnetic reconnection between different parts of the twisted flux commenced. At this point, a twist of "a few turns'' up to a "highly twisted'' field was reached, depending on the vortex velocity (quantitative measures of the twist were not given). These findings, indicative of ideal MHD stability up to substantial amounts of twist, may be markedly influenced by the high inertia of the simulation plasma in the upper parts of the box.

Twists in the range (3-15)$\pi$ were estimated in a sample of prominences with helical patterns near the point of their eruption (Vrsnak et al. 1991).

In the present paper we study a system that is similar to the one in Amari et al. (1996b) in several respects and we obtain similar indications for eruptive behaviour. We check for the existence of a critical twist beyond which no stable equilibrium exists by attempting to relax the system for different values of applied twist. Deviations from a force-free configuration are monitored as the twist is increased. We also investigate the influence of the inertia and the driving velocity on the dynamics by studying several initial density profiles, by including the continuity equation in the integration, and by reducing the driving velocity below commonly used values. Finally, by varying the vortex distance, we confirm that a low-lying highly sheared volume between the two vortices has negligible influence on the results.

We leave a study of the stabilizing influence of overlying current-free closed flux and the tendency of coronal loops to approach a constant coronal cross section for a subsequent investigation.

Following a description of the numerical method in Sect. 2, we study the quasi-static and dynamic phases of the evolution of a twisted flux tube, including the influence of the footpoint distance, initial density profile, and vortex velocity, in Sect. 3. The relaxation runs are presented in Sect. 4, the relationship to sigmoidal soft X-ray loops is discussed in Sect. 5, and our conclusions are given in Sect. 6.

   
2 Numerical model

   
2.1 Basic equations

We integrate the compressible ideal MHD equations employing the simplifying assumption of vanishing plasma beta:

   

$$\partial_t\rho -\nabla\cdot(\rho~\vec{u}),$$ = (2)

$$\rho~\partial_{t}\vec{u} -\rho~(~\vec{u}\cdot\nabla)~\vec{u} +\vec{j}\mbox{\boldmath$\times$ }\vec{B} +\nabla~\cdot\tens{T},$$ = (3)

$$\partial_{t}\vec{B} \nabla\mbox{\boldmath$\times$ }(~\vec{u}\mbox{\boldmath$\times$ } \vec{B}),$$ = (4)

where $\vec{B}$, $\vec{u}$, and $\rho$ are the magnetic field, the velocity, and the mass density, respectively. The current density is given by $\vec{j}=\mu_0^{~-1}~\nabla\mbox{\boldmath$\times$ }\vec{B}$. The last term in Eq. (3) is the viscous force, where $\tens{T}$ denotes the viscous stress tensor

   $$\tens{T}_{ij}=\rho~\nu~\left(\frac{\partial u_{i}}{\partial x... ...artial x_i} -\frac{2}{3}~\delta_{ij}~\nabla\cdot\vec{u}\right)$$

and $\nu$ is the kinematic viscosity. The bulk viscosity is neglected (since the viscous term is included only to facilitate relaxation toward equilibrium states and to improve numerical stability).

Since the plasma-beta is very small in the inner solar corona, we use the assumption $\beta=0$, which decouples the energy equation from the system. This guarantees that relaxed states, if found, are indeed force free, but on the other hand the dynamics is of course not fully described and our corresponding results must be interpreted with some caution. In particular, super-Alfvénic velocities are possible because the restoring force, $-\nabla p$, is not present if the fluid is accelerated by the Lorentz force $\vec{j}\mbox{\boldmath$\times$ }\vec{B}$. Nevertheless, the direction of dynamic evolutions that start from (near) force-free equilibrium states is determined by the Lorentz force and hence correctly described by our reduced system of equations.

It is also clear that force-free equilibria are independent of the density, which does not enter the force-free equation

$$\nabla\mbox{\boldmath$\times$ }\vec{B}=\alpha(\vec{x})\vec{B}.$$ (5)

If more than one equilibrium exists for the boundary conditions at t>0, then the evolutionary path, which depends on the distribution of density, determines which equilibrium is found. We do not expect this to be the case, because the initial potential field is uniquely determined by the boundary conditions and the applied vortex motions do not introduce knots or null points into the frozen-in field. This is supported by the simulations, which did not yield any indication for the existence of multiple equlibria at the same boundary condition. However, since the dynamics depends strongly on the inertia, we include the continuity equation as a step toward a more realistic description of the dynamics of arched flux tubes.

The equations are normalized by quantities which are derived from a characteristic length, $l_{\rm d}$, of the initial configuration, $\tilde{B}_{0z_{\rm max}}$ (the maximum normal component of the initial magnetic field in the bottom plane $\{z\!=\!0\}$), $\tilde{\rho}_0$, and $\tilde{v}_{\rm a0} =\tilde{B}_{0z_{\rm max}}/(\mu_{0}\tilde{\rho}_0)^{1/2}$ (the density and Alfvén velocity at the position of $\tilde{B}_{0z_{\rm max}}$). These are $l_{\rm d}$, $\tilde{v}_{\rm a0}$, $\tau_{\rm a}=l_{\rm d}/\tilde{v}_{\rm a0}$, $\tilde{\rho}_0$, $\tilde{B}_{0z_{\rm max}}$, and $\tilde{B}_{0z_{\rm max}}/\mu_{0}l_{\rm d}$ for $\vec{x}$, $\vec{u}$, t, $\rho$, $\vec{B}$, and $\vec{j}$, respectively. The normalization does not change the form of Eqs. (2)-(4), and Ampere's law becomes

$$\vec{j} = \nabla\mbox{\boldmath$\times$ }\vec{B}$$ (6)

after normalization. The normalized variables will be used henceforth.

   
2.2 Initial configuration

As initial condition we use a potential magnetic field created by two dipoles with vertical and oppositely directed moments $n\hat{\vec{z}}$ and $-n\hat{\vec{z}}$, located at (0,y₀,-z₀) and (0,-y₀,-z₀), respectively, with the plasma at rest:

 

$$\vec{B}_0 = \vec{B}_{0+}+\vec{B}_{0-},$$ (7)

  $$\begin{eqnarray*}\vec{B}_{0\pm} = \frac{\pm n} {[(x^{2}+(y \mp y_{0})^{2}+(z+z_... ... \mp y_{0})^{2}+2(z+z_0)^{2} \end{array} \right), \nonumber \\ \end{eqnarray*}$$

 

$$\vec{u}_{0} = 0.$$ (8)

The dipole depth is set to z₀=1.5 throughout this paper. The value of n follows from the condition $B_{0z_{\rm max}}=1$.

This field permits variation of the distance between the vortices without changing their profiles much, which will be used below to separate the influence of the shear that occurs for small dipole distances from the influence of the twist. Normalized dipole half distances y₀=0.5, 1, and 2 will be used. The field also yields a stronger concentration of the flux in the photosphere than the often used horizontal subphotospheric dipole. It is therefore probably a better model of the flux distribution in a sunspot group, but on the other hand, all of the overlying closed flux extends from the area of the flux concentration and becomes influenced by the photospheric vortices (Fig. 1). Although high overlying flux becomes only weakly twisted, there is, strictly speaking, no overlying flux that stays current-free. The stabilizing influence of such flux must be rather small or completely absent in our system, making it a useful reference case. We also note that the outer sign reversal of B_z on the bottom plane of the simulation is located sufficiently far from the vortex centres so that the velocities are nearly vanishing, and no current enhancements are induced, at these locations.

Figure 1: Field lines starting in the outer vortex area for a) the double-dipole field used in the present paper with dipole position y0=1 and  b) the single-dipole field used by Amari et al. (1996b). The plot shows field lines that start on the y axis in the range between the vortex centre and the outer point where the velocity has fallen to 0.1 of the peak value.

Three different initial density distributions are considered in this paper:

$$\rho_{0}=\left\{ \begin{array}{c} B_{0}^2 \\ \vert\vec{B}_0\vert \\ 1 \end{array} \right\}.$$ (9)

The first distribution is equivalent to uniform Alfvén velocity, $v_{\rm a}=1$, at the beginning of the calculation. This condition was also used by Amari et al. (1996b), but here it implies a lower density in the outer parts of the box, since the far field of a pair of antiparallel dipoles decays more rapidly with distance than that of a single dipole. Uniform density was used by Van Hoven et al. (1995) and Tokman & Bellan (2002). The distribution $\rho_0=\vert\vec{B}_0\vert$ is an arbitrarily chosen intermediate case.

We use different $\rho_0$ profiles for two reasons. First, they provide us with a check of the numerically obtained force-free equilibria, which are expected to be independent of the density. Second, the dynamics is strongly influenced by the inertia, with the low densities in the case $\rho _0=B_0^2$ enabling us to study the expansion of the flux from the flux tube axis, which may be relevant for the formation of a cavity in CMEs.

The calculations are started with a preparation phase in order to reduce the small spurious Lorentz forces that result from the discretization errors in the calculation of the initial current density. A relaxation of the initial configuration is performed with the vortices switched off. It is terminated when the maximum Lorentz force densities are ${\la}10^{-6}$. The time is then reset to zero.

   
2.3 Imposed vortex motions

Figure 2: a) Imposed velocity field for y0=1. b) Paths of fluid elements in the upper vortex for v0=10-2 after $75\tau _{\rm a}$. The vortex centre is located at $(0,1.03,-\Delta z)$. The thick arrows indicate the positions where two of the magnetic field line bundles in Fig. 3 start from. c) Profiles of the vortex $u_x(0,y,-\Delta z)$ ( thick line) and of the normal field component B0z(0,y,0) ( thin line).

A twisting velocity field is prescribed in the base plane throughout the driving phase of the simulation (actually it is applied in the first layer of grid points below $\{z\!=\!0\}$). It is chosen such that the velocity always points along the contours of B_0z(x,y,0), which leaves B_z invariant in this layer. This is achieved by

  

$$u_{x,y}(x,y,-\Delta z,t)=v_0f(t)\nabla^\perp~\{\zeta[B_{0z}(x,y,0,0)]\},$$ (10)

$$u_{z}(x,y,-\Delta z,t) = 0$$ (11)

with $\nabla^\perp:=(\partial_y,-\partial_x)$. A smooth function

   $$\zeta=B_{0z}^2\exp((B_{0z}^2-B_{0z_{\rm max}}^2)/\delta B^2),$$

chosen as in Amari et al. (1996b), determines the vortex profile. It is normalized such that the maximum of $\nabla^\perp\zeta$ equals unity. Its width is fixed by $\delta B=1$ in this paper, a value which yields a vortex width comparable to the width of the flux concentration. The velocities decrease rapidly for larger distances. These two vortices of the same sense are approximately centred at the dipole axes; their centre locations are at $\pm y_1$ with y₁=0.69, 1.03, 2.0 for y₀=0.5, 1.0, 2.0, respectively (Fig. 2).

In order to simulate a quasi-static evolution, the peak velocity, v₀, must be chosen much smaller than the normalized peak Alfvén velocity of the initial equilibrium (which is unity, or very close to unity, for all density profiles specified in Eq. (9)). Values v₀=10^-2and 10^-3 are employed (the former was most often used in the literature).

The function f(t) describes the temporal profile of the imposed twisting. The driving phase starts with a linear ramp ( $0\le t\le t_{\rm r1}=10$) to a value $f(t_{\rm r1})=1$, which is then held fixed. If a final relaxation phase is added, f(t) is linearly reduced to zero in an interval $t_{\rm r2}=1\mbox{--}5$, after which it is again held fixed.

   
2.4 Numerical implementation

As the physical domain we consider a finite cube of size $[-L_x,L_x] \mbox{\boldmath$\times$ } [-L_y,L_y] \mbox{\boldmath$\times$ } [0,L_z]$. Since both the initial configuration, Eqs. (7)-(9), and the twisting velocity field, Eqs. (10)-(11), possess line symmetry with respect to the z axis and since this symmetry is preserved throughout the evolution, the simulation domain can be restricted to the "half cube'' $[-L_x,L_x] \mbox{\boldmath$~\times~$ } [0,L_y] \mbox{\boldmath$~\times~$ } [0,L_z]$. Inside this area, the Eqs. (2)-(4) are discretized on a nonuniform Cartesian grid. We take $L_{x}\!=\!L_{y}\!=\!100$ and $L_z\!=\!200$. The grid spacing increases exponentially from a minimum $\Delta x_{\rm min}\!=\!\Delta y_{\rm min}\!=\!\Delta z_{\rm min}\!=\!0.03$in the vicinity of the origin to maximum values $\Delta x_{\rm max}\!=\!\Delta y_{\rm max}\!=\!9$ and $\Delta z_{\rm max}\!=\!12$ at the outer boundaries. In order to implement the boundary conditions, two layers of "ghost points'' are added to the box in each direction. Including these, the grid size is $195\times165\times150$.

Equations (2)-(4) are rewritten in a flux conserving form,

$$\partial_{t}\vec{U}+\partial_{x}\vec{F}(\vec{U}) +\partial_{y}\vec{G}(\vec{U}) +\partial_{z}\vec{H}(\vec{U}) =0,$$ (12)

where $\vec{U}\!=\!(\rho,\rho~u_{x},\rho~u_{y},\rho~u_{z},B_x,B_y,B_z)$is the vector of the integration variables and $\vec{F}$, $\vec{G}$, and $\vec{H}$ are nonlinear flux terms (e.g., Ugai & Tsuda 1977). A modified two-step Lax-Wendroff-method is employed for the numerical integration. It was adapted from that used by Ugai & Tsuda (1977) to our problem (see Appendix A). A variable time step is used according to the CFL criterion.

At the bottom boundary, the velocity is given by Eqs. (10)-(11). The tangential components of the magnetic field, B_x,y, are obtained by extrapolation (accurate to second-order). The normal component, B_z, is determined from $\nabla\cdot\vec{B}(z\!=\!0)=0$; it stays very close to $B_z(t\!=\!0)$throughout the simulation. For simplicity, the density is kept fixed at its initial value, $\rho(x,y,-\Delta z,t)=\rho_0$. Since the velocity points along the contours of B_z, while $\rho_0={\rm const.}$ at the contours of $\vert\vec{B}\vert$, this choice introduces an error for $\rho$ in the bottom plane if those contours disagree. This error is tiny for $y_0\ge1$and still relatively small for y₀=0.5, and it is irrelevant as long as the evolution stays quasi-static (see Sect. 2.1). We expect its influence on the dynamical evolution to be smaller than the effect of neglecting the pressure gradient force for all considered values of y₀. Closed upper and lateral outer boundaries are implemented by specifying $\rho=\rho_0$, $\vec{u}=0$, and $\vec{B}=\vec{B}_0$ (with $\vec{u}=0$ set also at the first layer of grid points inside these boundaries). The line-symmetric boundary conditions at the front plane $\{y\!=\!0\}$ are realized by even mirroring of $\rho$, u_z, B_x, and B_y at the zaxis (i.e., f(x,-y,z)=f(-x,y,z)) and by odd mirroring of u_x, u_y, and B_z (i.e., f(x,-y,z)=-f(-x,y,z)).

The huge variation of the density for the initial conditions $\rho _0=B_0^2$ and $\rho_0=\vert\vec{B}_0\vert$ in the chosen large box presents a challenge for any numerical scheme. Good stability properties of the scheme are also required by the stresses that develop between the twisted and the surrounding flux (especially if the latter is rigid by itself or embedded in fluid with high inertia) and by the strong dynamics that develop for strong twist. Although the Lax-Wendroff scheme is rather diffusive and although viscosity is included, further stabilization is needed for many runs. This is achieved by applying "artificial smoothing'' as introduced by Sato & Hayashi (1979) to $\rho$ and $\rho\vec{u}$after each integration step in addition to the stabilization by viscosity (see Appendix B for details). Both viscosity and artificial smoothing can be prescribed nonuniformly. We adopted the philosophy of supplying as much numerical diffusion as required by the different runs, while at the same time trying to keep identical numerical diffusion for runs that have similar initial conditions and are directly compared with each other in subsequent sections. There are two groups. Most runs with continuous driving in the base plane have a large ellipsoidal "inner region'' with uniform viscosity and smoothing and outward increasing numerical diffusion in the remainder of the box (where the densities are very low) and at the bottom plane (where the stresses are applied). All other runs have uniform viscosity and smoothing parameters in the whole box. Details are given in Appendix B and Tables B.1 and B.2. All analysis presented in this paper refers to states with the twisted flux tube completely inside the inner region or to states obtained with uniform viscosity and smoothing in the whole box.

The numerical diffusion gradually reduces the strong density gradients of $\rho _0=B_0^2$ and $\rho_0=\vert\vec{B}_0\vert$, leading to a slow increase of the density in the outer part of the box during the calculations. The Alfvén velocity then starts to decrease with growing distance from the dipoles also for $\rho _0=B_0^2$, but the drop of its height profile remains similar in magnitude to the conditions in the corona (e.g., Dulk & McLean 1978) for several $10^2\tau_{\rm a}$, i.e., during the whole driving phase of the runs presented in this paper as well as during short (successful) relaxations. This modified density distribution appears to be closer to coronal conditions than the initial $\rho _0=B_0^2$. The decline of the Alfvén velocity becomes less realistic during relaxation runs of ${\sim}10^3\tau_{\rm a}$ duration, which are considered unsuccessful anyway (Sect. 4).

   
3 Quasi-static and dynamic evolution of twisting flux tubes

In this section we describe and analyse the evolution of a continuously twisted system. First a typical set of parameters y₀, $\rho_0(\vec{x})$, and v₀ is chosen and characteristics of the rise of the twisted flux, the evolution of the magnetic energy and current density during the rise, and the deviation from a force-free configuration are discussed. We then proceed to a parametric study. Varying y₀ enables us to separate the influence of twist and shear, using different assumptions for $\rho_0(\vec{x})$ permits studying the influence of inertia on the evolution of twisted low-density coronal flux, and a smaller v₀ keeps the system closer to a force-free configuration. The effect of lowering the inertia is found to be particularly significant; it points to a possible explanation of the typical three-part structure of CMEs. See Table B.1 for an overview of the parameters used in the various runs.

   
3.1 Flux tube expansion

We choose a simulation with y₀=1, $\rho _0=B_0^2$, and v₀=10^-2 as our reference case. It is described in the following two subsections. All field lines of the initial potential field are symmetrical with respect to the $\{y\!=\!0\}$ plane. Most of the field lines that start in the region where the twisting velocity field is imposed are closed and connect the half planes $\{x,y>0,0\}$ and $\{x,y<0,0\}$. The central field line connecting the vortex centres has an initial apex height h₀=1.22 and length L₀=3.50. A flux-normalized volume can be assigned to each field line by the expression

   $$F(x,y,t) = \int_{C} \vert\vec{B}\vert^{-1} {\rm d}s,$$

where the integral is to be evaluated along the field line C starting at (x,y,0). It is a useful diagnostic of the central field line (Amari et al. 1996b). Its initial value in our reference run is F₀=F(0,y₁,0)=14.5.

Figure 3: Different evolution of the three groups of magnetic field lines that are directly influenced by the vortex motions: "central field lines'' (starting in the vicinity of a vortex centre) form the rising and expanding twisted flux tube; "inner field lines'' (initially starting between the vortex centres) form a weakly expanding inverse-S-shaped structure; and "outer field lines'' (initially starting beyond a vortex centre) show the strongest rise and expansion and a forward S shape. The reference run (y0=1, $\rho _0=B_0^2$, and v0=10-2) is shown at a) t=0 and b) t=75. The start regions of the selected inner and outer field lines are marked by the heavy arrows in Fig. 2b.

As the vortices are switched on, the volume permeated by field lines rooted in the area of the vortices starts to expand in all three spatial directions (as can be expected from the analogous behaviour of a twisted cylindrically symmetric configuration, in which the growing axial component of the field acts like additional plasma pressure, Sturrok 1994). These field lines now assume helical shapes and can be divided into three groups which evolve differently (Fig. 3; see also Tokman & Bellan 2002). Field lines starting in the vicinity of the vortex centres (central field lines) experience the strongest twist at the bottom plane (Fig. 2) while their footpoints do not move significantly. They form a central twisted flux tube which ascends and also expands laterally. Field lines initially starting in the region between the vortex centres (inner field lines) become sheared. They do not show a strong height expansion; instead they are streched, since their footpoint distance is initially increasing. Their projection onto the bottom plane shows a distinct inverse S shape, which is a direct consequence of the enforced motion of their footpoints. Finally, field lines initially starting beyond the vortex centres (outer field lines) exhibit the most rapid rise and lateral expansion while their footpoint distance is initially decreasing. The S-shaped projection of the outer field lines has a bending opposite to that of the inner field lines, due to the opposite displacement of the footpoints. The central flux tube starts to develop a weak inverse S shape (as the inner field lines), but reverses this early in the evolution (which is then still quasi-static) to show the same bending as the outer flux during the rest of the simulation. It should be noted that there is no sharp transition between these flux systems; it is not possible to find a clear boundary between the central flux tube and the outer or inner flux. The sigmoidal (S- or inverse-S) shape of the twisted flux is discussed further in Sect. 5.

The evolution of the central twisted flux tube can be divided into two markedly different phases, a quasi-static phase and a strongly dynamic phase. This can clearly be seen in the plots of its geometrical parameters and rise velocity in Fig. 4. The velocity of the flux tube, although permanently increasing, stays small ( $\vert\vec{u}\vert\ll v_{\rm a}$) up to $t\approx75$. The relaxation runs discussed in Sect. 4 show that the flux tube indeed stays nearly force free and, hence, runs through a sequence of near-equilibrium states in this phase. During $t\approx75\mbox{--}130$ the flux tube is rapidly accelerated to near Alfvén velocity, after which the expansion proceeds nearly linearly in length and width, i.e., nearly self-similarly (we find approximately $F\propto L^{2.2}$ in the quasi-static phase and $F\propto L^{2.9}$ in the dynamic phase). The rapid acceleration and the subsequent fast expansion both indicate that the system cannot settle to a new equilibrium after $t\approx75$. These findings are completely similar to those obtained by Amari et al. (1996b) for a different model of the initial magnetic field and density and support their conclusion that flux tubes erupt for sufficiently strong twist.

Figure 4: Evolution of height, length, volume, and apex rise velocity of the central field line of the twisted flux tube for the run shown in Fig. 3.

Figure 5 shows the temporal evolution of the velocity profile along the z axis (where u_x and u_y vanish). During the whole simulation, the outer flux expands even more rapidly than the central flux tube. (Super-Alfvénic velocities are reached because the pressure gradient is not included in the momentum equation; see Sect. 2.1. The expansion of the outer flux is slowed down only as it runs into the region of enhanced numerical diffusion; this was confirmed by a comparison run with a larger inner region of uniform numerical diffusion.) Obviously, the outer flux is not pushed upwards by the central flux tube in the present simulation but moves away from the central flux tube, primarily in response to its own footpoint twisting. This expansion of the outer flux is present already in the sequence of near-equilibria in the quasi-static phase. The main expansion occurs in the dynamic phase, provided the density decreases sufficiently rapidly with distance from the dipoles so that the effect of the inertia remains weak (particularly in our runs with $\rho _0=B_0^2$). The expansion includes the central flux in this run - opposite to the pinching expected from the nearly uniform twisting of the central flux. To see the origin of this apparent discrepancy, we consider the radial component of the Lorentz force density $\vec{f}=\vec{j}\mbox{\boldmath$\times$ }\vec{B}$ in a cylindrically symmetric (straight) flux tube, $\vec{B}=(0,B_\varphi(r),B_z(r))$, taken here as an approximation to our loop-shaped twisted flux tube:

$$f_r = -\left(\frac{1}{r}B_\varphi^2+\frac{1}{2}\partial_rB^2\right).$$ (13)

The first term in this expression always points inwards, but the second term may have either sign; it causes expansion of the flux if B²decreases sufficiently rapidly with distance from the axis. We have verified this to be the case in the bottom plane of our system and also in the cross-sectional plane at the apex of the central field line.

Figure 5: Velocity along the z axis at different times for the run shown in Fig. 3. The times correspond to the datapoints in Figs. 4a-c. Plus signs indicate the position of the respective apex heights of the central field line.

A pinching of the central flux towards the axis occurs for flatter functions B²(r) (e.g., Mikic et al. 1990). However, for a finite radial range of twist ( $B_\varphi\ne0$), the second term always leads to expansion of the flux at greater radii ( $B^2(r\!\to\!\infty)\to0$). In an ideal plasma, a volume of reduced density surrounding twisted core flux is necessarily formed. Pileup of flux and plasma must occur at the boundary to the unperturbed flux rooted beyond the photospheric vortex (although the resulting density profile may remain monotonously decreasing outwards if the initial profile is very steep and the vortex profile is smooth, as in our reference run). The resulting configuration bears striking resemblance to the three-part structure observed in many CMEs, which consists of a highly twisted prominence embedded in a cavity and a denser arc at the outer edge of the cavity (e.g., Hundhausen 1999). Further simulations including the pressure gradient term and initial magnetic field and density models that are better tuned to the conditions expected in the corona will be required to check for quantitative agreement of this effect with CME observations.

We note that the effects of pinching and expansion of flux may roughly cancel out if photospheric vortex motions are applied to a subarea of a photospheric flux concentration that is displaced from the point of maximum field strength, as in Klimchuk et al. (2000). The second term in Eq. (13) is then expected to change sign twice along all vortex revolutions closer to the centre than the point of maximum field strength. This may underlie the approximate constancy of the cross section of the flux tube for different amounts of twist found by Klimchuk et al. (2000).

Figure 6: Evolution of  a) magnetic energy;  b) current through the half plane $\{x,y>0,0\}$ ( plus signs: main current, asterisks: return current, diamonds: total current);  c) and  d) peak current density and peak Lorentz force density, respectively, in the whole box for the run shown in Fig. 3.

   
3.2 Evolution of energy, currents, Lorentz forces, and twist

As the magnetic field is changed due to the applied vortex motions, it becomes nonpotential, the total magnetic energy increases, and high current densities develop. The relative magnetic energy W(t)/W₀, where W₀ is the magnetic energy of the initial potential field, increases monotonously during the whole simulation (Fig. 6a). Only a weak tendency to saturate is seen after the onset of the dynamic phase. The energy, which is injected into the system by the vortex motions at its base, is then not only converted to magnetic energy but an increasing fraction is converted to kinetic energy of the plasma. At the end of the simulation, the kinetic energy amounts to about seven percent of the total energy increase of the system.

The strongest current concentration forms in the vicinity of the vortex centres due to the twist, however, also the shear between the vortices leads to buildup of currents in a central sheet that is formed cospatial with the inner field lines (Fig. 7). The shear and the resulting current increase with decreasing vortex distance y₀(Sect. 3.3). The total current through the bottom plane of our physical domain ( $\vert x\vert\le L_x$, $\vert y\vert\le L_y$) must vanish, since the applied vortices are compact (their velocities at the lateral boundaries are negligibly small). A net current through the bottom plane of the simulation domain ( $\vert x\vert\le L_x$, $0\le y\le L_y$) remains in our reference run, due to the shearing motions near the origin (Fig. 6b). This current nearly vanishes for larger vortex distance (it did not exceed 0.026 for y₀=2).

Figure 7: Isosurfaces $\vert\vec{j}\vert=0.2~j_{\rm max}$ (left) and $\vert\vec{j}\vert=0.5~j_{\rm max}$ (right) for the run shown in Fig. 3 at t=75.

Figure 8: $\vert\vec{j}\times\vec{B}\vert/\vert\vec{j}\vert\vert\vec{B}\vert$ along the central field line at t=54 ( solid), t=75 ( dotted), and t=102 ( dashed) for the run shown in Fig. 3.

Figure 9: Twist of the central flux tube for the run shown in Fig. 3.  a) Injected global (end-to-end) twist as a function of footpoint radius $r_{\rm foot}$ at t=75.  b) Evolution of the twist at the apex of the flux tube with $r_{\rm foot}=0.01$ ( asterisks) and $r_{\rm foot}=0.1$ ( diamonds). Also shown are the global twist ( dashed), $\alpha L/\Phi _{\rm a}$ at the apex ( triangles), and the (scaled) apex height, 0.1h/h0 ( solid line).

Despite the buildup of strong current densities, the Lorentz force densities remain rather small. Figures 6c, d show that the peak Lorentz force density on the grid is more than one order of magnitude smaller than the peak current density on the grid during the whole simulation. This is expected for the quasi-static phase but may be surprising for the dynamic phase, where it reflects the fact that the twisted flux can expand easily due to small inertia. These plots do not provide a quantitative measure of the deviation from a force-free state, since $\vert\vec{j}\vert$ and $\vert\vec{f}\vert$ do not necessarily peak at identical grid points. In Fig. 8 we plot $\vert\vec{j}\times\vec{B}\vert/\vert\vec{j}\vert\vert\vec{B}\vert$, a measure of the angle between $\vec{j}$ and $\vec{B}$, along the central field line for various times. This shows clearly that a deviation from a force-free state develops at the flux tube apex at the beginnig of the dynamic phase ($t\ga75$). The deviation continues to grow during the dynamic phase. Of course, during the whole driving phase, there must be Lorentz forces near the footpoints of the flux tube where the twisting is imposed, but these are hardly reflected in Fig. 8 because $\vert\vec{j}\vert$ and $\vert\vec{B}\vert$ are much larger at the footpoints so that a much smaller angle is sufficient to support the necessary Lorentz forces. Another measure of the deviation from a force-free configuration is given by the variation of the force-free parameter $\alpha=\vec{j}\cdot\vec{B}/B^2$ along the field lines, which is shown in Sect. 4.

The global twist of the central flux tube, $\Phi_{\rm g}$, or $2\pi$ times the number of enforced windings that the field lines make about its axis, is calculated from the imposed vortex velocity at the bottom boundary as

$$\Phi_{\rm g}=\frac{2}{r_{\rm foot}}~\bar{u}(r_{\rm foot})~t,$$ (14)

where $\bar{u}(r_{\rm foot})$is the average vortex velocity at a circle with radius $r_{\rm foot}$ centred at $(0,y_1,-\Delta z)$. This quantity is plotted in Fig. 9a. It is seen that the inner part of the flux tube ( $r_{\rm foot}\la0.1$) is nearly uniformly twisted (as is already clear from Fig. 2c). In addition we monitored the local twist at the apex of the central flux tube, $\Phi_{\rm a}$, which we define analogously to Eq. (1) as

$$\Phi_{\rm a}=\frac{L~\vert B_\varphi(r_{\rm apex})\vert} {r_{\rm apex}B_{\hat{z}}(r_{\rm apex})},$$ (15)

using local cylindrical coordinates r, $\varphi$, and $\hat{z}$, oriented at the flux tube axis at the apex, and the length L of the central field line. Here $r_{\rm apex}$ is determined for every $r_{\rm foot}$from the condition of flux conservation; it is found to obey $r_{\rm apex}\approx (B_{0z_{\rm max}}/B_{\rm apex})^{1/2}r_{\rm foot}$. The implied assumption of circular cross section at both the feet and the apex of the flux tube is satisfied to a good approximation for all runs with $\rho _0=B_0^2$ and during the quasi-static phase of all other runs. The value of $\Phi_{\rm a}$ is actually obtained as an average along the circle of radius $r_{\rm apex}$ at the loop apex. Figure 9b shows that also $\Phi_{\rm a}$ is nearly uniform for $r_{\rm foot}\la0.1$. In the following, the twist will always be determined for $r_{\rm foot}=0.03$.

Since the simulation indicates nearly uniform twist for radii $r_{\rm apex}(r_{\rm foot}\le0.1)\ll L$, and rather small Lorentz forces, it is appropriate to compare the inner part of the twisted flux tube with the known solution for a uniformly twisted force-free cylindrical eqilibrium by Gold & Hoyle (1960). We plot the quantity $\alpha L/\Phi _{\rm a}$, with the force-free parameter $\alpha$ taken at the apex of the central field line, in Fig. 9b. There is a close agreement with the exact value $\alpha L/\Phi=2$ at the axis of the Gold-Hoyle equilibrium during the whole simulation, which indicates that the configuration stays relatively close to a force-free state in the vicinity of the apex of the central field line.

Figure 9b shows the temporal evolution of $\Phi_{\rm g}$ and $\Phi_{\rm a}$. They nearly agree during the quasi-static phase, indicating that the injected twist is then distributed nearly uniformly along the central flux tube. The onset of the dynamic phase occurs at $\Phi_{\rm g}\sim2.6\pi$ and $\Phi_{\rm a}\sim2.1\pi$ - values relatively close to the critical twist for onset of instability in the Gold-Hoyle equilibrium, $\Phi_{\rm GH}\approx2.49\pi$ (Hood & Priest 1981). A characteristic of the fast expansion of the central flux tube in the dynamic phase of the simulation is nearly constant apex twist, which stays close to the value at the onset of the dynamic phase while the global twist is linearly rising. The implied scaling $\alpha\propto L^{-1}$ is consistent with the nearly self-similar expansion noted above.

   
3.3 Influence of the central shear

Figure 10: Evolution of apex height and rise velocity of the central field line for $\rho _0=B_0^2$, v0=10-2, and different dipole half distances: y0=0.5 ( plus signs); y0=1 ( diamonds); y0=2 ( asterisks). h0=0.67, 1.22, 3.24 for y0=0.5, 1, 2, respectively.

Figure 11: Evolution of apex height and rise velocity of the central field line for y0=1, v0=10-3, and different initial density profiles: $\rho_0=\vert\vec{B}_0\vert$ ( plus signs); $\rho _0=B_0^2$ ( asterisks).

As a first step of our parametric study, we investigate to which extent the main properties of the dynamical evolution discussed above depend on the presence of shear motions between the vortices. Varying y₀influences mainly the amount of shear while the amount of twist stays nearly unchanged. The shear-driven current densities are negligibly small for y₀=2 and increase for decreasing y₀, reaching $\approx$60 percent of the peak current density induced by the twist for y₀=0.5(see Fig. 7 for the case y₀=1).

Studies of sheared magnetic arcades which are translationally invariant in the direction of shear have shown a rapid rise of the arcade in the sequence of ideal MHD equilibria if the footpoints on either side of the arcade axis are displaced from the initial potential configuration by more than about 2.5 footpoint distances (e.g., Amari et al. 1996a; Mikic & Linker 1994). Such displacements appear to be much larger than observed displacements in the photosphere but they do occur near the origin in our system for sufficiently long duration of the imposed vortex motions. (A similar situation existed in the previous investigations of flux tube twisting by Amari et al. 1996b and Tokman & Bellan 2002.) Although Fig. 5 does not contain indications that the rise of the central flux tube is driven from the region of the inner field lines, the question remains how strongly the general evolution is influenced by the shear. Amari et al. (1997) found that the eruption of the twisted flux tube occurred also for smaller shear than in Amari et al. (1996b), but they did provide neither an information on how the shear was varied nor quantitative relations between the characteristics of the shear and the eruption.

Figure 10 shows that the evolution of the twisted flux tube is qualitatively similar for all three values of y₀. A quasi-static and a dynamic phase can be distinguished in all cases, in agreement with the expectation (from Fig. 5 and Amari et al. 1997) that a sufficiently twisted flux tube erupts also in the absence of shear.

It is also seen that the system with the smallest amount of shear (and hence with the smallest total induced current in the box), y₀=2, becomes dynamic at the smallest amount of twist and reaches the highest expansion velocities. This surprising trend results from a secondary effect associated with the variation of y₀ in our initial field: the amount of outer flux that closes in our box decreases and its average driving velocity increases for increasing y₀.

   
3.4 Influence of the initial density distribution and the driving velocity

The distribution of density should have no influence on the evolution of our system in the quasi-static phase if this phase is indeed a non-bifurcating succession of neighbouring force-free equilibria and the driving velocity is small enough so that the inertia does not prevent reaching the neighbouring equilibria. It is, however, expected to have a strong influence on the expansion of the system in the dynamic phase, where an acceleration to near-Alfvénic velocities was observed in our reference run. The initial density distributions in Eq. (9) differ by several orders of magnitude at heights reached by the central and outer flux in the dynamic phase.

We have started this part of the parametric study using y₀=1, v₀=10^-2, and the three distributions $\rho_0(\vec{x})$ as given in Eq. (9) and found that the rise of the central flux tube differs by a small but noticable amount between the three cases already in the time interval in which our reference run was quasi-static. The rise shows a systematic lag, increasing with higher density, which indicates that those cases require a smaller v₀. (This is confirmed by the relaxation runs discussed in Sect. 4.) Huge differences occur for times $80\la t\la200$ (the dynamic phase of the reference run), in which the rise velocity of the central flux tube reaches only $\approx$0.2 for $\rho_0=\vert\vec{B}_0\vert$ and $\approx$0.02 for $\rho_0=1$.

In Fig. 11 we compare the cases $\rho _0=B_0^2$ and $\rho_0=\vert\vec{B}_0\vert$ for v₀=10^-3. In the quasi-static phase, the evolution of the central flux tube is now indeed nearly independent of the chosen $\rho_0$. In the dynamic phase, the rise velocities and the corresponding accelerations again show a strong dependence on the density, but the overall dynamic evolution is nevertheless qualitatively similar.

It is possible to create very strongly twisted flux tubes by choosing high densities in the outer parts of the box, e.g., $\rho_0=1$, and high driving velocities, v₀>10^-2 (Tokman & Bellan 2002). Such states are only seemingly quasi-static (the velocities are kept low by the inertia so that $\vert\vec{u}\vert\ll v_{\rm a}$, but the flows continue to rise for several $10^2~\tau_{\rm a}$) and cannot be taken as indication for ideal MHD stability of highly twisted flux tubes. Both of the above conditions cannot be met by the majority of coronal flux if it is twisted by photospheric motions. Prominences may be an exception, however. Strong overlying flux is expected to provide the prime contribution to the stability of twisted flux tubes in the corona.

Figure 12 shows that the outer flux for $\rho_0=\vert\vec{B}_0\vert$evolves similarly to the reference case (Fig. 5) but with a much weaker expansion. This difference is clearly due to the inertia, since it occurs only in the dynamic phase. Modeling the formation of a cavity around an erupting twisted flux tube thus requires a careful treatment of the density.

Figure 12: Velocity along the z axis at different times for the run with y0=1, v0=10-3, and $\rho_0=\vert\vec{B}_0\vert$. The plus signs indicate the position of the respective apex heights of the central field line.

The comparison of the runs with y₀=1, $\rho _0=B_0^2$, and different driving velocity in Figs. 4a, d and 11 (and also in Fig. 17 below) shows that the overall evolution is qualitatively identical and that the quasi-static phase and the point of transition to dynamic behaviour are also quantitatively very similar if the elapsed time is transformed to accumulated twist (equivalent to roughly a scale factor of 10 between the time axes). In the dynamic phase, far greater apex heights are reached by the slowlyer driven system, which has more time to evolve than the faster driven system for the same $\Phi_{\rm g}$. This substantiates the interpretation in Amari et al. (1996b) and in Sect. 3.1 that the flux tube erupts if its twist exceeds a critical value.

   
4 Relaxation

Figure 13: Apex heights of the continuously twisted or relaxed flux tube during the quasi-static phase of the runs with y0=1. Plus signs: continuously twisted, $\rho_0=\vert\vec{B}_0\vert$, v0=10-2; diamonds: continuously twisted, $\rho _0=B_0^2$, v0=10-2; asterisks: continuously twisted, $\rho _0=B_0^2$, v0=10-3; squares: relaxed states obtained from the run with $\rho_0=\vert\vec{B}_0\vert$, v0=10-2; and triangles: relaxed states from the run with $\rho _0=B_0^2$, v0=10-3. (The uppermost triangle does not represent a relaxed state, see text for details.)

The relaxation runs serve two purposes: first, to confirm that the twisted flux tube in our system evolves closely along a sequence of force-free equilibria in the quasi-static phase and second, to substantiate our interpretation that there is a critical amount of twist above which the flux tube cannot find a new equilibrium within the simulation box. To start a relaxation, all velocities in the bottom layer were reduced to zero within 1-5 Alfvén times and then kept at zero. The flux tube is assumed to have reached a relaxed force-free state if all of the following conditions are satisfied: (a) velocity of flux tube apex continuously falling for at least $120\tau_{\rm a}$ to values below $3\times10^{-5}$; (b) integrated Lorentz force density in box continuously falling in the same time interval; and (c) $\alpha\approx{\rm const.}$ along the central field line. Table B.2 gives an overview of the parameters used in the various relaxation runs.

Figure 13 shows that the twisted flux tube in our runs with $\rho _0=B_0^2$ or $\vert\vec{B}_0\vert$ and v₀=10^-3 or 10^-2 indeed evolves closely along a sequence of force-free equilibria as long as the system remains quasi-static. The agreement improves for smaller densities in the outer part of the box ( $\rho _0=B_0^2$), which permit the central and outer flux to move more easily in response to the footpoint motions, and, of course, for smaller driving velocity. The apex heights of the relaxed states are nearly indistinguishable from the apex heights of the continuously twisted run with $\rho _0=B_0^2$ and v₀=10^-3 up to a global (end-to-end) twist $\Phi _{\rm g}\approx 2.5\pi$. This is true for both the relaxations started from this same run (triangles in Fig. 13) and the relaxations started from the run with $\rho_0=\vert\vec{B}_0\vert$ and v₀=10^-2 (squares in Fig. 13). Hence, the obtained equilibria are independent of the density distribution and vortex velocity, as expected.

Figure 14: Force-free parameter $\alpha$ along the central field line in the runs with y0=1 at a global twist of $\Phi _{\rm g}\approx 2.5\pi$. Dot-dashed: continuously twisted, $\rho_0=\vert\vec{B}_0\vert$, v0=10-2; solid: continuously twisted, $\rho _0=B_0^2$, v0=10-2; dotted: continuously twisted, $\rho _0=B_0^2$, v0=10-3; dashed: after relaxation from $\rho _0=B_0^2$, v0=10-3.

A sensitive diagnostic of the approach to a force-free state is given by the variation of the force-free parameter $\alpha$ along the field lines. In force-free states it is constant along each field line. In Fig. 14 we plot $\alpha$ along the central field line for the three runs shown in Fig. 13 and a relaxed state, all at nearly the same twist, $\Phi _{\rm g}\approx 2.5\pi$. Small deviations from a force-free state due to the driving at the footpoints are obvious for the two runs with v₀=10^-2, while the run with v₀=10^-3 and the relaxed state are nearly perfectly force free.

Figure 15: Relaxation runs near the critical twist, started from a system with y0=1 that was twisted using $\rho _0=B_0^2$ and v0=10-3. a) Apex height of twisted flux tube and b) total Lorentz force in the box during relaxation at $\Phi _{\rm g}=2.53\pi$ ( triangles), $2.75\pi$ ( squares), and $2.95\pi$ ( asterisks). The apex height for $\Phi _{\rm g}=2.95\pi$ is scaled to fit into the figure.

Figure 16: Lorentz force density at the z axis for relaxation runs withy0=1, $\rho _0=B_0^2$, v0=10-3. Thin lines are used for the run at $\Phi _{\rm g}=2.53\pi$ ( t=721, 801, 1041), thick lines are used for the runat $\Phi _{\rm g}=2.95\pi$ ( t=841, 1004, 1935). Solid lines show the system at the start of the relaxation, dotted lines show intermediate states, and dashed lines show the end of the runs. The Lorentz force density changes sign along the z axis during the relaxation at $\Phi _{\rm g}=2.53\pi$, so |fz(z)| is plotted in this case for $t\ge 801$. Crosses mark the respective apex heights of the twisted flux tube.

The following four figures illustrate our attempts to relax configurations with higher global twist, which have already entered the dynamic phase. In Fig. 15 we plot the apex height of the twisted flux tube and the Lorentz force density integrated over the box during three representative relaxation runs. These were started from that continuously twisted run with y₀=1 which remained closest to equilibrium during the quasi-static phase ( $\rho _0=B_0^2$ and v₀=10^-3) at global twists of $\Phi _{\rm g}=2.53\pi$, $2.75\pi$, and $2.95\pi$. The Lorentz force densities at the z axis are shown in Fig. 16 for $\Phi _{\rm g}=2.53\pi$ and $2.95\pi$. At $\Phi _{\rm g}=2.53\pi$, the flux tube reaches an equilibrium easily and quickly. At $\Phi _{\rm g}=2.95\pi$, it is impossible to relax the system. The apex height continues to rise and the total Lorentz force (which is dominated by the contributions from the loop footpoint areas) does not decay for more than $10^3\tau_{\rm a}$. The Lorentz forces at the z axis even increase during the relaxation attempt, primarily at the position of the flux tube. This shows that the gradual decrease of the rise velocity, apparent in Fig. 15a, is caused by the continuous slow rise of the density in the outer parts of the box, which is due to the numerical diffusion acting on the strong density gradient during the long integration. The relaxation run was terminated as the still ascending flux tube reached a height of $\approx$50 h₀. The relaxation at $\Phi_{\rm g}=2.75\pi$ is also considered unsuccessful. Using the diffusion parameters given in Table B.2, a continuous slow rise of the flux tube is found, while the total Lorentz force in the box stops decreasing at a value more than three times higher than the one reached in the relaxation at $\Phi _{\rm g}=2.53\pi$(Fig. 15). The Lorentz forces at the z axis first decrease and then start to increase, again primarily at the position of the flux tube. If the numerical diffusion is reduced to values applied during the twisting phase, the run becomes numerically unstable (as does the run at $\Phi _{\rm g}=2.95\pi$). We have therefore repeated the relaxation, applying an artificial gradual reduction of the density (up to a factor 50), to counteract the density increase due to the numerical diffusion. The flux tube then continues its rise to far greater heights than those shown in Figs. 15 and 13. This run was terminated at t=1350 and h/h₀=17. In contrast, the relaxed flux tube at $\Phi _{\rm g}=2.53\pi$ stays at its equilibrium position if a similar artificial density reduction is applied. Finally, an intermediate state of the system at $\Phi_{\rm g}=2.75\pi$ during the relaxation at a relatively low level of the Lorentz forces (t=900) was twisted further up to $\Phi _{\rm g}=2.95\pi$. The relaxation started from this stage gave the same results as the previous attempt at this amount of global twist.

Similar results were obtained from relaxation attempts using the system with $\rho_0=\vert\vec{B}_0\vert$ and v₀=10^-2 at twists $\Phi_{\rm g}>2.5\pi$. These runs begin with smaller velocities, due to the higher inertia, but again, the system did not relax at $\Phi_{\rm g}=2.91\pi$ during more than 10³ Alfvén times. At $\Phi_{\rm g}=2.65\pi$ a very slow evolution was found, in which the Lorentz forces at the z axis first decreased but started to increase again slowly after several 10²Alfvén times.

Figure 17: Apex height as a function of global twist for the runs with y0=1, $\rho _0=B_0^2$, and v0=10-2 ( dashed) and v0=10-3 ( solid). Triangles show apex heights at the end of relaxation runs for states calculated from the run with v0=10-3; the final two datapoints show still rising states, which were not relaxed at termination of the runs.

Figure 17 shows the apex height as a function of global twist for our reference run and the otherwise similar run with v₀=10^-3 in the full range of global twists covered by the simulations. Apex heights at the termination of relaxation runs (which represent only lower limits for $\Phi_{\rm g}\ge2.75\pi$) are added. The sharp transition above $2.5\pi$ is obvious and leads us to conclude that there is a critical twist, $\Phi_{\rm c}$, in our system above which either a stable equilibrium ceases to exist, or a sequence of stable equilibria exists but is characterized by huge increase of height for small increase of twist. The critical twist of the considered configuration lies in the range

   $$2.5\pi<\Phi_{\rm c}<2.75\pi.$$

Our numerical simulations cannot yield a final decision between the two possibilities, due to the finite box size. However, a strong indication is obtained from a comparison with related analytical results for the asymptotic behaviour of axially symmetric twisted force-free fields, which also show a rapid expansion for sufficiently large twist (or shear) angle. Sturrock et al. (1995) found that the flux rises exponentially with the square of the angle of footpoint displacement, which implies that a closed-field equilibrium exists for any finite angle. On the other hand, Uzdensky (2002) argued in favour of opening (loss of closed-field equilibrium) at a finite angle. The expression by Sturrock et al. (1995) fits the simulation data very well (Fig. 18) - but only for twists smaller than the above critical value. This indicates that loss of stable equilibrium and dynamic opening of the field occur in the loop-shaped configuration at $\Phi_{\rm g}=\Phi_{\rm c}$.

For practical purposes, the decision between the two possibilities is of limited importance anyway, since the final apex height of the relaxation run at $\Phi _{\rm g}=2.95\pi$ implies the rise of the flux tube by at least a factor of 50, probably at a speed comparable to the Alfvén velocity, which would observationally appear as an eruption. If we scale the footpoint distance in this simulation to a typical distance of the leading and the following polarity in an active region, $\sim$70 000 km, the minimum apex height of 50h₀ at $\Phi _{\rm g}=2.95\pi$ scales to $3~R_\odot$.

The critical twist in our system is only slightly larger than the critical twist of the uniformly twisted Gold-Hoyle equilibrium, $\Phi_{\rm GH}\approx2.49\pi$. This appears reasonable for the following two reasons. First, the inner part of our twisted flux tube is nearly uniformly twisted and, hence, similar to the Gold-Hoyle equilibrium, except for the curvature, which is relatively weak ( $r_{\rm curvature}/r_{\rm apex}(r_{\rm foot}\!=\!0.03, \Phi_{\rm g}\!=\!2.53\pi)\approx5$). Second, the surrounding untwisted flux and possibly also the less twisted outer flux in our system are expected to have a weak stabilizing effect. On the other hand, our system is less stable than the one studied by Van Hoven et al. (1995). We suppose that this difference is mainly due to the action of overlying arcade-like nearly potential flux in their system. An investigation of the stabilizing effect of such flux is left for future work.

Figure 18: Comparison of the apex heights of the twisted flux tube after (or at termination of) relaxation (symbols as in Figs. 13 and 17) with the rise characteristics $h\propto\exp({\rm const.}\cdot\Phi_{\rm g}^2)$ of axially symmetric configurations from Sturrock et al. (1995). A least-squares fit through the first eight datapoints is shown.

Figure 19: Field lines of the magnetic field ( solid) and current density ( dotted) in the central flux tube projected onto the bottom plane for y0=1, $\rho _0=B_0^2$, and v0=10-3. a) Driven phase, $\Phi _{\rm g}=2.53\pi$ (t=721). b) Driven phase, $\Phi _{\rm g}=2.95\pi$ (t=841; start of relaxation). c, d) Two states in relaxation run at $\Phi _{\rm g}=2.95\pi$ (t=1065 and t=1785).

We suppose that there is no complete loss of equilibrium as the twist is raised beyond $\Phi_{\rm c}$, but simply a loss of stability. The long-wavelength ideal kink mode is a probable candidate for the destabilization. If the twisted flux tube becomes kink-unstable, helical perturbations are expected to grow. Such perturbations are difficult to diagnose, since the flux tube develops a helical shape already in the quasi-static phase. Some indications for the growth of a kink perturbation can, however, be found. Projections of the twisted flux tube onto the bottom plane, plotted in Fig. 19 for the run with y₀=1, $\rho _0=B_0^2$, and v₀=10^-3, show two of them. First, the rotation of its upper part about the z axis, a measure of the helical distortion from the initially planar shape, accelerates at the onset of the dynamic phase and continues to increase during the relaxation at $\Phi _{\rm g}=2.95\pi$. Second, the angle between $\vec{j}$ and $\vec{B}$ increases in the upper part of the flux tube all through the relaxation run (as indicated already in Fig. 8). A similar, slightly weaker growth of the helical distortion of the flux tube occurred in the two-step relaxation run at $\Phi _{\rm g}=2.95\pi$described above. The magnetic energy in the box was found to decrease during these relaxation runs but only by a small amount ($\sim$0.1 percent), which is due to the low densities in the upper parts of the box and the only slightly supercritical value of the twist.

Figure 20: Contours of $\vert\vec{j}(z\!=\!0)\vert$ for the run with y0=1, $\rho _0=B_0^2$, v0=10-3 after relaxation at $\Phi _{\rm g}=2.06\pi$. Contours with $\alpha >0$ (jz>0) are shown solid, contours with $\alpha <0$ (jz<0) are shown dotted. (The contours of $\alpha$ in the bottom plane are similar to the contours of $\vert\vec{j}\vert$, see Fig. 3 in Török & Kliem 2002.) Bz>0 in the whole area shown. The vortex centre is marked by a plus sign, the location of peak current density is marked by an asterisk. These are the start points of the two flux tubes shown in Fig. 21.

   
5 Sigmoidal shape

As mentioned in Sect. 3.1 and apparent from Fig. 19, the projection of the central flux tube onto the bottom plane is forward-S shaped during most of the evolution. Since the force-free parameter $\alpha$ is negative in the central flux tube for our chosen sense of twist, this shape is opposite to the typical bending of sigmoidal soft X-ray loops in the solar corona (Pevtsov et al. 1997; Nightingale et al. 2002). It has been suggested by Magara & Longcope (2001) that such loops are distinguished from the neighbouring flux by a higher current density (supposing that the heating is related to the current density). Figure 20 shows that the highest current density in our system does not occur on the central field line. Although the displacement of the current density peak from the vortex centre is rather small, its effect on the shape of the associated flux bundle is significant. This flux bundle has an inverse S shape (Fig. 21), before and after relaxation, in agreement with the tendency of soft X-ray sigmoids that are rooted in $\alpha <0$ regions. It rises with increasing twist in a manner similar to the central flux tube but with a somewhat smaller apex height (typically by a factor $\approx$0.8). Figure 21 also shows that the type of the sigmoidal shape (forward or inverse S) of force-free fields is in general not only related to the sign of $\alpha$. Our simulations, which use a different method of twist injection into the corona than Magara & Longcope (2001), support their conjecture that current-related heating distinguishes sigmoidal soft X-ray loops from the surrounding coronal fields.

Figure 21: Projection onto the bottom plane of the flux tubes that extend from the vortex centre and from the area of maximum current density; their apex heights are 2.6 and 2.2, respectively. Field lines of $\vec{B}$ (solid) and $\vec{j}$ (dotted) are shown. Parameters are as in Fig. 20. The contours $\alpha =0$ (dashed) are included.

   
6 Conclusions

Our conclusions can be summarized as follows.

(1) A twisted magnetic flux tube that is created from a potential field and continuously driven by slow vortex motions at its footpoints starts to evolve quasi-statically through a sequence of stable force-free ideal-MHD equilibria to a very good approximation. The sequence of equilibria is characterized by a weak expansion with increasing twist, which agrees with the scaling obtained by Sturrock et al. (1995) for axisymmetric configurations, and is found to be independent of the distribution of density and the driving velocity.

The field lines become helical, with forward and inverse S shapes coexisting closely. The S shape of the flux tube that extends from the vortex centre is opposite to the S shape of the flux tube that extends from the (slightly displaced) region of maximum current density; the latter is consistent with the typically observed bending of S-shaped soft X-ray loops.

(2) There is a critical end-to-end twist, $\Phi_{\rm c}$, above which the flux tube does not find a new equilibrium in the simulation box but enters a dynamic phase, developing expansion velocities no longer negligible in comparison to the Alfvén velocity. Then there is either no stable equilibrium, or stable equilibria are characterized by huge increase of apex height for small increase of twist beyond $\Phi_{\rm c}$. A comparison of the simulation results with the scaling of flux expansion given by Sturrock et al. (1995) suggests that there is no stable equilibrium for $\Phi>\Phi_{\rm c}$. Both cases have the appearance of an eruption if the simulation is scaled to parameter values characteristic of the solar corona. This supports similar conclusions reached by Amari et al. (1996b) for a different initial potential field.

The quantitative properties in the dynamic phase depend strongly on the density, $\rho(\vec{x})$, and on the driving velocity in the base plane. There are indications that helical perturbations grow during this dynamic phase and also during relaxation runs started at supercritical twist.

(3) All these properties remain unchanged if the vortices are so close to each other that shear is created between them in addition to the twist, at least as long as the current density induced by the shear remains smaller than the current density induced by the twist.

(4) The critical twist in our system (with normalized footpoint half-distance $y_0\!=\!1$) lies in the range $2.5\pi<\Phi_{\rm c}<2.75\pi$.

(5) The outer flux that surrounds the twisted flux tube and is also rooted in the area of the vortices experiences an expansion from the central field line of the flux tube as the twist rises. The inner part of the flux tube may experience a pinching or an expansion, depending on the structure of the initial field, $\vec{B}_0(\vec{x})$, and vortex, u_x,y(x,y,0). The expansion can be discerned in a simulation only if the density is modeled sufficiently realistically such that the Alfvén velocity is not too rapidly decreasing with height. It leads to a volume of reduced density surrounding the twisted flux tube. In a more realistic model than the one studied here, it is expected to lead also to pileup of magnetic flux and density at the boundary to the untwisted and unsheared flux further out. Such a configuration is similar to the three-part structure, consisting of a prominence, a cavity, and an outer arc, which is often observed in CMEs.

Acknowledgements

This paper benefitted from very constructive comments by the anonymous referee, T. Amari, and N. Seehafer. We gratefully acknowledge also comments by, and discussions with, P. Bellan, P. Démoulin, K. Galsgaard, A. Hood, J. Klimchuk, T. Neukirch, K. Shibata, V. Titov, and L. van Driel-Gesztelyi. B.K. acknowledges the hospitality of the Solar Group at the Nobeyama Radio Observatory, where part of this work was completed. This investigation was supported by BMBF/DLR grants No. 50 OC 9901 2 and 01 OC 9706 4, and by EU grant No. HPRN-CT-2000-00153. The John von Neumann-Institut für Computing, Jülich granted Cray computer time.

   
Appendix A: The modified two-step Lax-Wendroff method

Starting from Eq. (12) the two-step Lax-Wendroff method is derived using a second order Taylor expansion of $\vec{U}$ in time, where the time derivatives of $\vec{U}$ are replaced by spatial derivatives of $\vec{F}$, $\vec{G}$, and $\vec{H}$ (e.g., Richtmyer & Morton 1967). One iteration of the scheme consists of two steps. We have used a modified version of the first step, which calculates auxiliary values at $t=(2n+1)\Delta t$ from the values at $t=2n\Delta t$,

 

                     U_i,j,kⁿ⁺¹ = U_i,j,kⁿ

$$-\frac{\Delta{t}}{2\Delta{x}} \left(F_{i+1,j,k}^n-F_{i-1,j,k}^n\right)$$

$$-\frac{\Delta{t}}{2\Delta{y}} \left(G_{i,j+1,k}^n-G_{i,j-1,k}^n\right)$$

$$-\frac{\Delta{t}}{2\Delta{z}} \left(H_{i,j,k+1}^n-H_{i,j,k-1}^n\right)$$ (A.1)

and a standard version of the second step, which calculates the physical values at $t=(2n+2)\Delta t$,

 

                     U_i,j,kⁿ⁺² = U_i,j,kⁿ

$$-\frac{\Delta{t}}{\Delta{x}} \left(F_{i+1,j,k}^{n+1}-F_{i-1,j,k}^{n+1}\right)$$

$$-\frac{\Delta{t}}{\Delta{y}} \left(G_{i,j+1,k}^{n+1}-G_{i,j-1,k}^{n+1}\right)$$

$$-\frac{\Delta{t}}{\Delta{z}} \left(H_{i,j,k+1}^{n+1}-H_{i,j,k-1}^{n+1}\right).$$ (A.2)

In standard versions of the Lax-Wendroff scheme, an averaging term

$$\begin{eqnarray*}\frac{1}{6}~\left(U_{i-1,j,k}^n+U_{i+1,j,k}^n +U_{i,j-1,k}^n+U_{i,j+1,k}^n+U_{i,j,k-1}^n+U_{i,j,k+1}^n\right) \end{eqnarray*}$$

is used in place of the first term on the r.h.s. of Eq. (A.1). This has the consequence that the effective first order spatial derivatives at the grid point (i,j,k) are calculated by averaging the derivatives obtained at the six nearest neighbours [(i-1,j,k), (i+1,j,k), (i,j-1,k), (i,j+1,k), (i,j,k-1), and (i,j,k+1)], thus enhancing the diffusion and improving the numerical stability, but it is not advantageous if one requires derivatives of quantities which exhibit a strong spatial variation, as for example the magnetic field components in our model in the regions of strong field. This is true even for good grid resolution. Therefore, we modified the scheme by replacing the averaging term with U_i,j,kⁿ. Now the spatial derivatives are calculated directly at the grid points and the numerical diffusion as well as the discretization errors of the current density are reduced. We performed test runs with both versions of the scheme and found similar dynamical evolution of our model systems.

The stability of the Lax-Wendroff scheme requires also that the viscous term be included in the second step of the iteration only (e.g., Roache 1972).

   
Appendix B: Treatment of the numerical parameters

Artificial smoothing is applied to selected integration variables $\psi=(\rho,\rho~u_x,\rho~u_y,\rho~u_z)$ after each full time step by the replacement (Sato & Hayashi 1979)

  $$\begin{eqnarray*}\vec{\psi}_{i,j,k}^{n}&\longrightarrow&\sigma\vec{\psi}_{i,j,k}... ...\left. \vec{\psi}_{i,j,k-1}^{n}+\vec{\psi}_{i,j,k+1}^{n}\right). \end{eqnarray*}$$

The smoothing is not applied to $\vec{B}$ to avoid unphysical currents and Lorentz forces. It would have been desirable to smooth the velocities only, but numerical stability required to smooth also the density.

 

 

Table B.1: Main parameters of the simulations with continuously applied twisting. $\sigma _{u}$ and $\sigma _{\rho }$ are the smoothing parameters of the variables $\rho~u_{x,y,z}$ and $\rho$, respectively. The subscripts "in'' and "out'' refer to the "inner region'' of the simulation domain with uniform $\nu$ and $\sigma$ and to the rest of the box, respectively. $\nu _{\rm in}=0.1$ in all runs.

v0	$\rho_0$	y0	$\nu_{\rm out}$	$1\!-\!\sigma_{u}$	$1\!-\!\sigma_{\rho,\rm {in}}$	$1\!-\!\sigma_{\rho,\rm {out}}$
0.01	B02	0.5	0.1-1	0.0005	0.0005	0.0005-0.1
0.01	B02	1.0	0.1-1	0.0005	0.0005	0.0005-0.1
0.01	B02	2.0	0.1-1	0.0005	0.0005	0.0005-0.1
0.01	$\vert\vec{B}_0\vert$	1.0	0.1-1	0.0005	0.0005	0.0005-0.1
0.01	1	1.0	0.1-1	0.0005	0.0005	0.0005-0.1
0.001	B02	1.0	0.1	0.0005	0.0005	0.0005-0.01
0.001	$\vert\vec{B}_0\vert$	1.0	0.1	0.0001	0.0001	0.0001

 

 

Table B.2: Main parameters of the relaxation runs. These runs had uniform $\nu =0.1$ and $\sigma$ (as given in the table). $t_{\rm relax}$ is the duration of the runs. $t_{\rm r2}=1$ for the runs with $\rho _0=B_0^2$, and $t_{\rm r2}=5$ for the runs with $\rho_0=\vert\vec{B}_0\vert$ (except the one at $\Phi_{\rm g}=2.50~\pi$, in which also an artificial reduction of $\rho$ was applied).

y0	$\rho_0$	v0	$(1-\sigma)$	$\Phi_{\rm g}$	$t_{\rm relax}$
1.0	B02	0.001	0.0004	2.06 $\pi$	120
			0.0004	2.30 $\pi$	120
			0.001	2.53 $\pi$	280
			0.002	2.75 $\pi$	753
			0.005	2.95 $\pi$	1095
1.0	$\vert\vec{B}_0\vert$	0.01	0.0001	0.39 $\pi$	120
			0.0001	0.86 $\pi$	200
			0.0001	1.33 $\pi$	200
			0.0001	2.28 $\pi$	600
			0.0005	2.50 $\pi$
			0.0001	2.65 $\pi$	1120
			0.0001	2.91 $\pi$	1358

In particular the very low densities in the case $\rho _0=B_0^2$ required a high numerical diffusion in the outer regions of the simulation domain. Therefore, we implemented both the artificial smoothing of the density and the viscosity (Eq. (2.1)) such that they could be prescribed nonuniformly. In an "inner region'', including about one third of the box size in each direction ( $\sqrt{x^2+y^2+(0.5z)^2}\le30$), but excluding a bottom layer described below, the parameters determining the diffusion, $\nu$ and $(1-\sigma)$, were prescribed uniformly. Outside this region, the viscosity and the density smoothing were either continued uniformly or permitted to increase exponentially towards the outer boundaries, as the stability of the individual runs required. The bottom layer, where the highest gradients and Lorentz forces occur, also required higher diffusion than the rest of the inner region for stability. Nonuniform smoothing of density and velocity, increasing from the values in the inner region to $(1-\sigma)=0.05$ in a layer between $\{z\!=\!z_1\}$ and $\{z\!=\!0\}$ was applied, with z₁=1.2 in the driving phase and z₁=0.3 during relaxation. The chosen values of the numerical parameters are compiled in Table B.1 for the continuously driven runs and in Table B.2 for the relaxation runs.

References

 
• Aly, J. J. 1991, ApJ, 375, L61 NASA ADS
   
• Amari, T., Luciani, J. F., Aly, J. J., & Tagger, M. 1996a, A&A, 306, 913 NASA ADS
   
• Amari, T., Luciani, J. F., Aly, J. J., & Tagger, M. 1996b, ApJ, 466, L39 NASA ADS
   
• Amari, T., Luciani, J. F., Aly, J. J., & Mikic, Z. 1997, in Coronal Mass Ejections, ed. N. Crooker, J. A. Joselyn, & J. Feynman (Washington: AGU), 101
   
• Antiochos, S. K., DeVore, C. R., & Klimchuk, J. A. 1999, ApJ, 510, 485 NASA ADS
   
• Brandt, P. N., Scharmer, G. B., Ferguson, S., et al. 1988, Nature, 335, 238 NASA ADS
   
• Dulk, G. A., & McLean, D. J. 1978, Sol. Phys., 57, 279 NASA ADS
   
• Einaudi, G. 1990, in Physics of Magnetic Flux Ropes, ed. C. T. Russell, E. R. Priest, & L. C. Lee (Washington: AGU), 43
   
• Einaudi, G., & van Hoven, G. 1983, Sol. Phys., 88, 163 NASA ADS
   
• Forbes, T. G. 2000, J. Geophys. Res., 105, 23 153
   
• Gold, T., & Hoyle, F. 1960, MNRAS, 120, 89 NASA ADS
   
• Hofmann, A., Rendtel, J., Aurass, H., & Kalman, B. 1987, Sol. Phys., 108, 151 NASA ADS
   
• Hood, A. W. 1992, Plasma Phys. & Contr. Fusion, 34, 411 NASA ADS
   
• Hood, A. W., & Priest, E. R. 1981, Geophys. Astrophys. Fluid Dynamics, 17, 297 NASA ADS
   
• Hundhausen, A. 1997, in Coronal Mass Ejections, ed. N. Crooker, J. A. Joselyn, & J. Feynman (Washington: AGU), 1
   
• Hundhausen, A. 1999, in The Many Faces of the Sun, ed. K. T. Strong, J. L. R. Saba, B. M. Haisch, & J. T. Schmelz (New York: Springer), 143
   
• Kempf, P. 1910, Astron. Nachr., 185, 197 NASA ADS
   
• Klimchuk, J. A. 2001, in Space Weather, ed. P. Song, H. J. Singer, & G. L. Siscoe (Washington: AGU), 143
   
• Klimchuk, J. A., Antiochos, S. K., & Norton, D. 2000, ApJ, 542, 504 NASA ADS
   
• Leka, K. D., Canfield, R. C., McClymont, A. N., & van Driel-Gesztelyi, L. 1996, ApJ, 462, 547 NASA ADS
   
• Lites, B. W., Low, B. C., Martinez Pillet, V., et al. 1995, ApJ, 446, 877 NASA ADS
   
• Magara, T., & Longcope, D. W. 2001, ApJ, 559, L55 NASA ADS
   
• Mikic, Z., & Linker, J. A. 1994, ApJ, 430, 898 NASA ADS
   
• Mikic, Z., Schnack, D. D., & van Hoven, G. 1990, ApJ, 361, 690 NASA ADS
   
• Moore, R. L., Sterling, A. C., Hudson, H. S., & Lemen, J. R. 2001, ApJ, 552, 833 NASA ADS
   
• Nightingale, R. W., Shine, R. A., Brown, D. S., et al. 2002, in Multi-Wavelength Observations of Coronal Structure and Dynamics - Yohkoh 10th Anniversary Meet., ed. P. C. H. Martens, & D. Cauffman (Amsterdam: Pergamon), 149
   
• Pevtsov, A. A., Canfield, R. C., & Zirin, H. 1996, ApJ, 473, 533 NASA ADS
   
• Pevtsov, A. A., Canfield, R. C., & McClymont, A. N. 1997, ApJ, 481, 973 NASA ADS
   
• Richtmyer, R. D., & Morton, K. W. 1967, Difference methods for Initial-Value Problems, 2nd ed. (New York: Wiley-Intersience), Ch. 13.4
   
• Roache, P. J. 1972, Computational Fluid Dynamics (Albuquerque: Hermosa), Ch. V-F
   
• Sakurai, T. 1979, PASJ, 31, 209 NASA ADS
   
• Sato, T., & Hayashi, T. 1979, Phys. Fluids, 22, 1189 NASA ADS
   
• Sturrock, P. A. 1991, ApJ, 380, 655 NASA ADS
   
• Sturrock, P. A. 1994, Plasma Physics (Cambridge University Press), Ch. 13.4-13.7
   
• Sturrock, P. A., Antiochos, S. K., & Roumeliotis, G. 1995, ApJ, 443, 804 NASA ADS
   
• Tokman, M., & Bellan, P. M. 2002, ApJ, 567, 1202 NASA ADS
   
• Török, T., & Kliem, B. 2002, ESA SP-506, 781
   
• Ugai, M., & Tsuda, T. 1977, J. Plasma Phys., 17, 337 NASA ADS
   
• Uzdensky, D. A. 2002, ApJ, 574, 1011 NASA ADS
   
• van Driel-Gesztelyi, L., Manoharan, P. K., Démoulin, P., et al. 2000, J. Atmospheric & Solar Terr. Phys., 62, 1437
   
• Van Hoven, G., Mok, Y., & Mikic, Z. 1995, ApJ, 440, L105 NASA ADS
   
• Vrsnak, B., Ruzdjak, V., & Rompolt, B. 1991, Sol. Phys., 136, 151 NASA ADS