© ESO, 2013

1. Introduction

Coronal mass ejections (CMEs) are one of the few solar phenomena that directly affect the Earth, therefore they are one of the main drivers of space weather (Schwenn 2006). Understanding their origin and propagation is therefore key, not only for the Earth but also for understanding how the Sun loses both magnetic flux and magnetic helicity (Low 1996; Lynch et al. 2005). Over the years, a wide variety of models have been put forward to explain the origin of CMEs, and a discussion of these can be seen in the reviews of Forbes et al. (2006) and Chen (2011). While a variety of models exist, one of the leading models for explaining CMEs is the flux rope ejection model (Forbes & Isenberg 1991; Amari et al. 2000; Fan & Gibson 2007). The flux rope ejection model may itself be split into two categories. Flux ropes formed quasi-statically from perturbations of already existing coronal arcades (Aulanier et al. 2005; Mackay & van Ballegooijen 2006a) or those formed dynamically during the emergence of magnetic flux (Manchester et al. 2004; Archontis & Hood 2008). Within the present paper we focus on those formed quasi-statically from existing coronal arcades.

For this category of flux rope, the life span normally undergoes three separate stages of evolution: the formation, equilibrium and eruption phases. To begin with, the flux rope forms over time periods of days to weeks as a coronal arcade is perturbed by photospheric motions and flux cancellation at a PIL (van Ballegooijen & Martens 1989; Forbes 1991). An example can be seen in Cheng et al. (2010) where the flux rope formation was observed over three days, along with the signature of flux cancellation at the photosphere. Similar dynamics have been observed by Liu et al. (2010), where a sigmoidal flux rope was formed by the joining of two J-shaped loops which then quickly erupt. The formation phase is then followed by a period during which the flux rope lies in near equilibrium with its surroundings. Evidence for this exists in the form of sigmoids (McKenzie & Canfield 2008), solar filaments (Mackay et al. 2010), and coronal cavities (Hudson et al. 1999). During both the formation and equilibrium phases the magnetic structure does not change significantly over time scales much longer than the characteristic Alfvén time.

Finally, there is the ejection phase where the flux rope loses equilibrium and is ejected out of the solar corona due to magnetic forces. This occurs over time periods of a few hours, where flux rope ejections are believed to be one of the main progenitors of CMEs. It is possible to track the evolution of a flux rope in observed CMEs from the ejection phase (Cheng et al. 2011) to their propagation into interplanetary space (Wood et al. 1999). While many flux ropes follow the evolution described above, there are also situations where flux ropes are simultaneously formed and ejected (Cheng et al. 2011; Archontis et al. 2009).

It can be seen that the formation and ejection of magnetic flux ropes involve a wide variety of time scales. The formation phase happens very slowly (days to weeks or months). Under these circumstances the formation can be approximated using a quasi-static zero-β model where the magnetic field evolves through a series of equilibrium states (Forbes 1991; Mackay & van Ballegooijen 2006a,b; Yeates & Mackay 2009; Yeates et al. 2010). In contrast the ejection phase is rapid, where a flux rope can travel at the speed of several hundreds of km s^-1 and within a few hours the coronal magnetic field reconfigures to a state significantly different from the pre-eruption state. One important aspect of the ejection is that the coronal plasma is compressed and heated. As a result, the plasma may no longer be in the low β-regime and full MHD modelling must now be applied.

The wide variety of time scales involved, in the formation of the flux rope and then its eruption, poses a considerable challenge to theoretical models. Within the present paper our aim is to combine two models, which we believe accurately reproduce individual aspects of flux rope formation and ejection over their individual time scales. We will first use the global non-linear force-free field (GNLFFF) model of Mackay & van Ballegooijen (2006a) which neglects the role of the plasma to consider the formation of the flux rope in the zero-β, quasi-static approximation. Once the flux rope forms and is about to erupt, we apply a full MHD approach, through using the stressed magnetic field configuration from the GNLFFF model as an initial condition in MHD simulations using the AMRVAC code (Keppens et al. 2012). With this we can accurately describe the full life span of the evolution of a flux rope from its slow generation to its rapid ejection.

While we choose to apply these two methods, full MHD simulations may in principle be used to simulate both the formation and eruption of the flux rope. However, presently this is not practicable for two reasons. The first is limited computational resources. To model only the eruptive phase which corresponds to approximately two hours of physical time, we require over 5000 cpu hours. If we were to model the full three weeks of formation and ejection, assuming a linear scaling, we would need three months for a single run. Additionally due to the small MHD signal travel time, a large number of iterations would be required over the three months of computation. This would significantly increase round-off errors and numerical diffusion in the AMRVAC code.

While our approach is one way of considering the origin of CMEs, other authors have tackled the modelling of flux rope ejections using different approaches. An example is that of Reeves et al. (2010) who modelled a solar eruption as the outcome of a slow evolution preceding a fast and sudden eruption. In addition, a number of authors have considered alternative methods for forming flux ropes which do not use flux cancellation. For example, Fan (2009) and Archontis et al. (2009) explain the flux rope formation in the framework of flux emergence. Savcheva et al. (2012) explain the formation of an observed flux rope due to the rotation of the foot points of an active region in the presence of magnetic diffusion. While in Amari et al. (2003) the flux rope is formed (and then ejected) as a result of shearing followed by flux cancellation. In addition to these studies, the ejection of flux ropes may be explained by several mechanisms in numerical models: Török & Kliem (2005, 2007), Fan (2010) and Aulanier et al. (2010) use the Torus instability as the main driver for a full ejection. Archontis & Hood (2012) describes the eruption of a flux rope in the context of flux emergence with an horizontal ambient magnetic field antiparallel to the emerging magnetic field (see also Antiochos et al. 1999). Roussev et al. (2012) consider the solar eruption as a reorganization of the solar corona after magnetic flux and helicity injection. Finally, in Amari et al. (2011) the convergence of foot points causes the ejection.

While there are the many models for the formation and ejection of flux ropes we have decided to focus only on the flux cancellation model in this paper, as first described by van Ballegooijen & Martens (1989) and simulated by Mackay & van Ballegooijen (2006a). An important aspect of this work is that we form the flux rope out of a smoothly varying continuous magnetic field configuration, which is different from many other techniques. This means that the coronal magnetic field is not predisposed to allow reconnection, as happens in some models which use distinct flux configurations. The outline of the paper is as follows. The coupling of the GNLFFF and the MHD models is non-trivial and will be explained in detail in Sect. 2 with additional information in Appendices A and B. In Sect. 3 we describe the outcome of our MHD simulation which produces the ejection of the flux rope. In Sect. 4 we discuss the results and finally draw conclusions in Sect. 5. It should be noted that this is the first stage in this process where we have made a number of simplifications due to the new technique being employed. In future studies these assumptions will be relaxed.

2. Method

2.1. Previous simulations

In order to pursue the objective of modelling the full life span of a flux rope from the slow quasi-static built-up to the rapid ejection, we adopt two different numerical models for the two different phases in the life of a flux rope. We couple the two techniques through using the output from one model as the initial condition for the other.

Mackay & van Ballegooijen (2006a) apply a GNLFFF model that produces flux ropes by means of the relative motion and interaction of two magnetic bipoles. In the model, one bipole is slightly more tilted with respect to the equator than the other. The two bipoles are then advected under the effect of solar differential rotation, meridional flow and surface diffusion. Differential rotation has the effect of shearing the magnetic field lines of the bipoles and produces an axial magnetic field over the PIL. This axial magnetic field is then localized above the PIL due to surface diffusion. As a result of flux cancellation a flux rope forms above the PIL of the more tilted bipole after 19 days of evolution. An illustration of this 3D magnetic field configuration at the time of flux rope formation can be seen in Fig. 1. An important aspect of this 3D magnetic field with flux rope is that it is produced out of a continuous magnetic field and in the area surrounding the flux rope no antiparallel magnetic field lines exist.

We identify this as the time when a flux rope has clearly formed in the simulation. If the GNLFFF simulation is continued, the newly formed flux rope becomes unstable and is ejected out of the computational box (a full description is given in Mackay & van Ballegooijen 2006a,b). However, the simulations of Mackay & van Ballegooijen (2006a) neglect the effect of the plasma and use a magnetic relaxation technique, so the ejection of the flux rope does not follow the true dynamics of an erupting system. In particular the time scale of the ejection is too slow compared to an observed eruption. For the purposes of the present paper we assume the evolution modelled by Mackay & van Ballegooijen (2006a) to be accurate until day 19. We then adopt the stressed non-potential magnetic field that exists on this day as the initial magnetic field condition for numerical MHD simulations using ARMVAC which will follow the dynamics of the eruption.

The coupling of the two models is a non-trivial process and is a major technical challenge of the presented work. In order to accomplish this some assumptions and approximations are required (Sect. 2.2). From a mere technical point of view, the successful combination of these two modelling techniques is an achievement which produces a new tool for long term studies.

2.2. Coupling of the two models

To couple the two models we need to generate a full set of MHD variables (density ρ, velocity v, total energy density e and magnetic field B = ∇ × A) from the GNLFFF variables (vector potential A and magnetic field B). In both codes the variables are defined in the spherical domain (r, θ, φ). However, only the magnetic field, B, is common to both codes (see Fig. 2). The three components of the magnetic field are therefore directly transported from the GNLFFF code to the MHD code as described in Appendix A. The transport is designed to preserve the gradients of B, thus the forces in the system. However for the other MHD variables assumptions need to be made.

As the velocity, v, in the GNLFFF model, is a non-physical relaxation velocity we do not import it into the MHD simulation. We choose to initially set v = 0 in the MHD model. However for cases where the imported magnetic field is not in equilibrium a non-zero velocity develops within the MHD simulation after one time step. To complete the set of MHD variables, we assume that the density of the plasma, ρ, is given by (1)where we include a background uniform distribution to avoid extremely low density values occurring where the magnetic field is weak. Typical parameters are chosen to be ρ₀ = 10^-17 g/cm³, B₀ = 1 G, ρ_bg = 10^-21 g/cm³. We also assume a uniform thermal pressure, p, (2)where p₀ = 0.006 dyne/cm².

2.3. ARMVAC code

Once we have a full set of MHD variables from the coupling of the models, we perform a MHD simulation. An important feature of our MHD simulations is that no photospheric boundary motions are applied. This means that any dynamics occurring within the MHD simulation are purely due to the stress developed in the GNLFFF simulation before the data is transferred into the MHD simulation. We use the AMRVAC code mainly developed at the KU Leuven to run the simulations (Keppens et al. 2012). The code solves the MHD equations in their conservative form and we make use of the magnetic field splitting technique (Powell et al. 1999) and consider (3)where B₀ is a static, intrinsic magnetic field and B₁ is the difference (B − B₀). The following equations must hold for B₀The MHD equations we solve are

where the expression for the total energy density e is (11)and γ = 5/3 is the ratio of specific heats, and t is the time.

The use of the magnetic field splitting technique allows us to remove any errors that arise when transporting the variables from the GNLFFF code to the MHD code, due to round-off, discretization, and truncation. In effect, it minimizes spurious forces that are numerically generated. In addition, by neglegting the energy associated to the term in Eq. (11) we avoid the thermal pressure being a residual subject to round off errors where the plasma β is extremely low.

In our simulation, B₀ is chosen to be the potential field deduced from the photospheric magnetic field from the GNLFFF simulation, and B₁ is simply computed as the difference between B of the GNLFFF and the prescribed B₀. Given the specified photospheric magnetic field, the resulting coronal potential magnetic field is unique.

To check the stability of this technique, we have performed several simulations with different B₀ and therefore B₁ (corresponding to the same B). Different choices of B₀ corresponds to different choices of the photospheric magnetic field compared to that transported from the GNLFFF model. The smaller B₁ becomes, the less magnetic flux is free to numerically reconnect. The potential magnetic field B₀ derived from the photospheric magnetic field at Day 19 of the GNLFFF simulation is the most natural choice, in order to have the smallest B₁, however the potential magnetic field derived from any photospheric magnetic field of the GNLFFF is a suitable B₀.

In the simulation we present here, the numerical diffusion has no role in driving the early evolution of the system, which is solely due to the initial condition being out of magnetohydrostatic equilibrium, i.e. the initial evolution is the same and independent of the choice of B₀. On the other hand, after tenths of Alvén times (~30), if the system has highly departed from the initial condition, the exact evolution of the system does depend to some degree on how the magnetic field numerically reconnects and the initial choice of B₀ can play a role. It should however be noted that in all simulations the ejection progresses in a similar manner no matter the choice of B₀. The simulation we present here is for the smallest B₁ possible.

To evolve the system of Eqs. (7)−(10) forward in time we use a Total Variation Diminishing Lax-Friedrichs scheme (Tóth & Odstrcil 1996).

The size of the computational domain is the same as in the GNLFFF model where it extends over 1.5 R_⊙ in radial extent starting from r = R_⊙. The colatitude, θ, spans from θ = 30° to θ = 100° and the longitude, φ, spans over 90°. The boundary conditions are treated with a system of ghost cells. Open boundary conditions are imposed at the outer boundary, reflective boundary conditions are set at the θ boundaries and the φ boundaries are periodic. At the lower boundary we impose constant boundary conditions taken from the first four θ-φ planes of cells derived from the GNLFFF model. Therefore the actual computational domain is composed of 128 × 128 × 128 cells. The validity of the technique to combine the two models is tested in Appendix B through importing an equilibrium magnetic field configuration.

3. Simulations of an erupting flux rope

3.1. Properties of the initial condition

Figure 3 illustrates the initial condition for the magnetic field in the MHD simulation of the erupting flux rope which comes from day 19 of the simulations of Mackay & van Ballegooijen (2006a). The magnetic field lines of the flux rope connect between the two polarities of the left hand side (LHS) bipole (red lines, Fig. 3). Overlying the flux rope, magnetic field lines create an arcade like structure. Three additional arcade systems are also present in the initial condition: one above the Right Hand Side (RHS) bipole, another connecting between the two bipoles (both shown as blue lines in Fig. 3) and finally an arcade lies high in the corona connecting between the external polarities of the bipoles (green lines, Fig. 3). The magnetic field configuration shown in Fig. 3 is a single smoothly varying continuous system.

In the lower corona where the magnetic field strength is approximately 10 Gauss, the prescribed initial density is around 10^-16 g/cm³. The density drops by several orders of magnitude in the outer corona to a final value of 10^-21 g/cm³ at 2.5  R_⊙, where the magnetic field intensity is ~0.01 Gauss. The pressure is initially uniform (p = 0.006 dyne/cm²) and in the domain β ranges from 10^-3 near the centre of the bipoles to 10³ at the outer boundary. The initial Alfvén speed is of the order of 10⁸ cm/s. The time step of the AMRVAC simulation is around dt = 5 × 10^-3 s, and snapshots of the simulation are sampled every Δt ≈ 70 s, which corresponds to about half an Alfvén time (τ_Alf = 120 s). Hereafter we express the time in units of τ_Alf.

To understand how the pre-stressed initial condition of the magnetic field will drive the dynamics, Fig. 4a shows a map of the radial component of the Lorentz force in the r − φ plane along the centre of the bipoles. In all plots white/pink denote positive values, blue zero values and green/black negative values. It is clear from Fig. 4a that there is a strong positive (radially outward) component of the Lorentz force localized under and around the flux rope in the initial condition.

3.2. Early phase of the ejection

At t = 0 τ_Alf the positive Lorentz force present underneath the flux rope pushes the flux rope and the plasma contained within it upwards and it immediately starts to rise. Therefore the initial condition as imported from the GNLFFF model results in an erupting system. Simultaneously the foot points of the flux rope are subjected to a twisting motion due to the θ and φ components of the Lorentz force (not shown). As the flux rope moves up a horizontal Lorentz force acts below the flux rope to push oppositely oriented magnetic field lines on either side of the PIL towards the PIL. When these magnetic field lines encounter one-another they reconnect below the flux rope. This creates deeply dipped magnetic field lines which increase the outward component of the Lorentz force.

Figure 5 illustrates this process, where in Fig. 5a we see two individual magnetic field lines at t = 2.61 τ_Alf. These magnetic field lines are typical of the arcade magnetic field lines that connect the outer and central part of the LHS bipole. These magnetic field lines are pushed towards one another and reconnect below the flux rope to produce a low, deeply dipped, magnetic field line underneath the flux rope. In Fig. 5b at t = 2.67 τ_Alf we see the reconnected magnetic field line where the magnetic field line has the same outer edge starting points as before. Through this process, by t = 2.67 τ_Alf the radial component of the Lorentz force under the flux rope doubles in magnitude compared to its initial value.

In Fig. 6 the bulk plasma velocity at t = 11.6 τ_Alf can be seen. Figure 6a shows a map of the radial component of the velocity in the r − φ plane, along the centre of the bipoles, while Fig. 6b shows a map of the radial component of the velocity in the θ − φ plane zoomed in at the bipoles. In addition, in Fig. 6b contours of the radial component of the magnetic field are shown (solid contours positive flux, dashed negative flux) and the green vectors denote the horizontal velocity field. In Fig. 6b the velocity vectors are typically of the order of 2 × 10⁶ cm/s while velocity vectors of magnitude lower than 1.3 × 10⁶ cm/s are not shown. Figure 6a shows that the plasma above the LHS bipole is ejected outwards with a velocity higher than 8 × 10⁶ cm/s. This motion is initially due to the positive radial Lorentz force present under the flux rope in the initial condition and secondly due to the reconnection of the magnetic field lines below the flux rope. Figure 6b clearly shows that in the proximity of the PIL of the LHS bipole (above which the flux rope lies) the plasma moves upwards. Simultaneously horizontal plasma motions converge towards the PIL. The convergence is most prominent in the zone around θ =  [70°,75°] and φ =  [38°,42°].

To follow the ejection of the flux rope and the motion of the plasma, in Fig. 7 (left column) density distributions (maps) and magnetic field lines (green lines) can be seen at various times. In the initial stages of the flux rope ejection, the arcades above the bipoles remains largely unchanged (compare Fig. 7a and b). By t = 17.4 τ_Alf the flux rope has risen to 1.4  R_⊙, where it has expanded and has a velocity of ~ 50 km s^-1. In addition, it can be clearly seen that dense plasma is ejected outwards by the motion of the flux rope, (Figs. 7a and b).

In this early phase of the dynamics, the flux rope that is ejected moves as a complete structure enclosed in the space of the arcade that lies above the LHS bipole. At t = 23.2 τ_Alf, the map of the radial Lorentz force in the r − φ plane through the centre of the bipoles illustrates the forces acting on the ejection (Fig. 4b). The Lorentz force above the LHS bipole is strongest below r = 1.3 R_⊙ where the axis of the flux rope is approximately at r = 1.4 R_⊙. In addition to the force that ejects the flux rope, the motion of the ejecting flux rope generates additional Lorentz forces above it. Directly above the flux rope the ejection increases the curvature of the magnetic field lines and produces a downward directed magnetic tension (highlighted by the black dashed contour in Fig. 4b). However, slightly higher than this there is a compression of the magnetic field which leads to a larger magnetic pressure and correspondingly upward directed Lorentz force which is approximately twice as strong as the downward component (Fig. 4b). It should be noted that in all of the simulations considered no matter the choice of B₀, the early stages of the eruption occur in the same way.

3.3. Propagation of the flux rope

After the onset of the eruption, the flux rope moves upward, confined by the overlying arcade of the LHS bipole. As it does so, it expands in θ and φ directions. The propagation of the ejected magnetic field can be seen in Fig. 7 (right column) where different magnetic flux systems are denoted with different colours. In Fig.7e the blue and yellow magnetic field lines denote magnetic field lines that in the initial condition connect between the outer edges and central part of the LHS bipole. The blue magnetic field lines start from the positive polarity and the yellow magnetic field lines from the negative polarity. For these magnetic field lines the same starting points are used in all subsequent images. In addition to this there are some magnetic field lines (not shown) that lie along the entire length of the PIL and belong to the flux rope. At high heights magnetic field lines of the large-scale overlying arcade (green lines) can be seen. Between the times shown in Figs. 7e and f, the flux rope rises and the blue/yellow arcade magnetic field lines reconnect underneath the flux rope and add magnetic flux to the flux rope. This can be seen through the two sets of magnetic field lines changing their connectivity in the early stage of the ejection to have intermingled footpoints. In the early stage of the ejection the external magnetic field lines are not significantly affected by the motion of the flux rope (Fig.7f) and the flux rope moves upwards as a coherent structure.

After t = 17.4 τ_Alf (Fig. 7g), the ejected flux rope starts to interact with the overlying magnetic field. At this point the magnetic field configuration becomes significantly more complex and part of the magnetic field of the flux rope changes connectivity. A new phase of the ejection starts where some of the magnetic field lines connecting the outer polarity of the LHS bipole change their connectivity and connect to the outer polarity of the RHS bipole: denoted as yellow lines in Fig.7g. This leads to two systems of magnetic field lines moving upwards: one connecting the two polarities of the LHS bipole (blue lines) and the other connecting the external polarities of the LHS and RHS bipoles (yellow lines). The ejection is now no longer a simple structure and the evolution follows this dynamics until t = 34.8 τ_Alf.

As the flux rope continues to move up, its magnetic field interacts with the open magnetic field lines that connect to the source surface at 2.5  R_⊙ and a further phase of magnetic reconnection occurs. Some of the original footpoints of the flux rope now become connected to the outer boundary (yellow lines) however it is still possible to identify a coherent magnetic structure that connects the two polarities of the LHS bipole (Fig.7h, blue lines). The evolution of the magnetic field lines follows a similar evolution until the end of the simulation. At t = 69.6 τ_Alf the ejection reaches the upper boundary and the dynamics of the simulation are no longer followed.

3.4. Propagation of density fronts, role of thermal pressure, and conservation of energy

The plasma within and around the flux rope is initially dragged upwards by the flux rope motion. This includes the high density region (ρ ~ 10^-16 g/cm³) which is concentrated inside the flux rope. A coherent ejection of plasma lasts in the simulation as long as the flux rope moves as a single identity (before t = 34.80 τ_Alf). During this phase, the plasma distribution is clearly tied to the magnetic configuration of the flux rope, and has a cylindrical shape (see Figs. 7a, b).

After the flux rope reconnects with the open magnetic field lines and is no longer a clearly identified structure the coherence of the plasma also decreases. As is shown in Fig. 7d the density distribution splits into two main parts: there is one column of density above the LHS bipole and a wide bow shaped front overlying both bipoles. In between the two density structures the bulk plasma velocity is still directed upwards while the two density structures are pushed sideways.

In order to examine whether the evolution of the system leads to an ejection of plasma, along with the flux rope, we follow a density front of ρ = 5 × 10^-18 g/cm³ in Fig. 8. The front starts from a height of 1.3  R_⊙ and is accelerated by the ejecting magnetic flux rope. The average speed of the front is about 100 km s^-1. Since the evolution is adiabatic, as the flux rope moves up and the plasma is displaced, the thermal pressure behind the flux rope drops. At the same time, the ejecting flux rope pushes aside the overlying high-β plasma and carries upward the low-β plasma contained within it. As a result, the low-β region containing the flux rope expands and rises (see white lines in Figs. 7a−d).

The initial dynamics of the erupting flux rope is solely determined by the magnetic field as there are initially no pressure gradients. Once the eruption is under way, the key effect of the pressure gradient is to act as a counter force to the motion of the flux rope, through compression ahead and decompression behind. However, as long as the ejection of the flux rope takes place in a low β regime, the pressure gradients are not able to counter the dynamics. This does however mean that, the initial background configuration can have a significant impact on whether or not the ejection hits the outer boundary at 2.5 R_⊙. An eruption which hits 2.5 R_⊙ always takes place as long as the density front can move into an environment where β < 1. Thus the initial position of the β = 1 surface determines whether or not the eruption will be a full eruption, or a quenched one. We have performed additional simulations with a background pressure ten times larger than presented in the simulation here. It clearly shows that the eruption stops as soon as the flux rope travels out of the low-β region surrounding the bipoles. This is illustrated in Fig. 9 where we show maps of density with magnetic field lines superimposed along with the β = 1 contour. The state of the system is shown at the initial and final times: t = 0 τ_Alf and t = 69.6 τ_Alf. Although the initial magnetic field of the simulations are identical, as the background thermal pressure is ten times higher, the contour at β = 1 is lower in the atmosphere and now lies within the high density region. The flux rope and the dense plasma are ejected as before in the initial stages but now they reach a height of 1.6  R_⊙ after t = 40 τ_Alf where the ejection is quenched. With the higher background pressure the rise of the flux rope is slower and it stops once the dense structure crosses over the β = 1 zone. This shows that the thermal pressure is critical in determining whether a full or quenched eruption takes place.

Finally, in order to check the validity of our simulation where we do not apply any external boundary driving we plot the variation of total energy in the domain as function of time in units of the initial total energy. The total energy is computed as the integral over the whole domain of the total energy, e which includes thermal, magnetic and kinetic energy. We show in Fig. 10 that the total energy is conserved over the simulation and as a result the eruption is not driven and is only due to forces present in the initial condition. At the end of the simulation approximately ~1% of the energy is lost, but this can be attributed to the open boundary conditions imposed at 2.5  R_⊙ which allows energy to escape the domain.

4. Discussion

In the work presented here we have followed the full life span of a flux rope from formation to eruption. In doing so we have developed a new technique for coupling distinct forms of simulations, where the slow formation of the flux rope is modelled by a quasi-static approximation (Mackay & van Ballegooijen 2006a) and the eruption through a full MHD simulation. We now discuss a number of issues that have arisen in the present simulations and will be considered in future studies.

4.1. Magnetic configuration of the system

In our simulation, the flux rope forms and exists in a smoothly varying continuous magnetic configuration produced as a result of photospheric motions and flux cancellation. Within the magnetic field configuration the stress and forces applied have also been developed in a consistent way as may occur on the Sun. In the immediate surroundings above the flux rope there is the bipolar arcade and above this there is a large scale external arcade. The orientation of the magnetic field lines in the arcade are similar to that of the flux rope. This means that to some extend, the magnetic configuration is unfavorable for an ejection to escape.

4.2. Do we model the generation of a CME?

In the early stages the flux rope is initially ejected outwards as a single structure and it is possible to clearly follow its evolution. Within the flux rope trapped plasma is expelled and can be clearly identified. However, as soon as the flux rope moves beyond the symmetrical overlying arcade of the LHS bipole the dynamics become significantly more complex due to the interaction of the flux rope with the external magnetic field. After magnetic reconnection occurs high up it is no longer possible to clearly identify the flux rope. However, once the flux rope is no longer clearly identified the plasma does maintain its upward motion. Through following a dense front of plasma it reaches the top of our domain at 2.5  R_⊙ travelling at a speed of around 100 km s^-1. While such height reached by the dense front is promising for producing a CME, more investigations need to be performed to unambiguously answer the question as to whether or not we have modelled the generation of a CME.

In our simulation the characteristic three-component structure of CMEs is only hinted at, as we have just an indication of dense and void regions, but not the typical light-bulb shape. However this shape as seen in LASCO images is normally seen at radii much larger than that considered in the present simulations. This may also be partly due to the fact that the flux rope exists in a smoothly varying continuous magnetic field configuration and there is not a distinct flux rope cavity in the overlying magnetic field. Although the velocity of the expelled front is enough to let the plasma reach 2.5  R_⊙ it only matches the slowest tail of the CME population. Nevertheless, with this initial study, based on an initial condition produced by slow quasi-static build up of stress using observed motion on the Sun, we have successfully modelled a plausible ejection for a CME.

4.3. What initiates the ejection?

The Lorentz force is the only force that exists in the initial condition, and it determines the initial evolution of the system. The initial condition of our simulation is in non-equilibrium and the radial component of the Lorentz force underneath and in the newly formed flux rope pushes it upwards. This force is a result of the formation of the flux rope such that it is too large to be held down by the overlying arcade. Therefore, the initial evolution and ejection of the flux rope is intrinsic to its existence in this scenario.

However, more complex dynamics occur later on in the simulation, when magnetic reconnection dominates the evolution of the flux rope. It occurs in two ways. First, the upward motion of the flux rope is enhanced by reconnection of magnetic field lines underneath the flux rope as it rises. Secondly, the magnetic field lines of the flux rope reconnect with the overlying open and closed magnetic field lines. This does not influence the onset of the dynamics of the ejection, but it is relevant for the geometry and shape of the ejection at later stages.

In summary, in the simulation there are different locations of reconnection. However magnetic reconnection only plays a major role for the ejection after the onset of the rise has occurred. Magnetic reconnection does not initiate the ejection, which starts because of a situation of non-equilibrium after the formation of the flux rope.

4.4. Is the role of thermal pressure important?

In the simulation, the thermal pressure is initially uniform. Though simple to implement it is not realistic for the solar corona. We however choose this as it does not add any additional forces in the gravity-less model. Within the presented simulation the key effect of the thermal pressure is to oppose the motion of the flux rope, through compression in front and rarefactions behind the flux rope. However, despite of the opposition from the thermal pressure the flux rope still reaches a height of 2.5  R_⊙. Additional simulations show that if the background pressure is increased to be 10 times larger then the eruption is quenched in the corona.

An important aspect to consider is that our simulation does not include gravity. The gravitational force exerted on the plasma may prevent it and the flux rope from escaping from the sun. If we assume a stratified atmosphere in the initial conditions, the gravity force would be initially balanced by an upward directed pressure gradient. If we model the gravity in a properly stratified atmosphere, the confining power of the gravity might be negated by the upward directed pressure gradient. Additional simulations show that in the majority of cases gravity does not prevent the ejection, but only slows down the ejection. These results will be the subject of a future paper.

4.5. Differences with the GNLFFF model

Our work aims to bridge two different phases of the flux rope life span: the built up and the ejection. It starts from day 19 of the GNLFFF model after the flux rope has formed. Our MHD simulation then models the eruption of the flux rope over a period of a few hours after its construction and the onset of non-equilibrium.

In the simulations of Mackay & van Ballegooijen (2006a) after day 19, the flux rope also starts to rise in the corona. However it only reaches the higher corona after a number of days as the GNLFFF model is designed to damp evolution. In contrast the MHD code presented here properly describes the dynamic evolution of the system as it departs from the initial condition: in the presented simulation the ejection of the flux rope is complete in slightly more than two hours.

An important aspect of the simulations is that during the ejection above the LHS bipole the RHS bipole does not change significantly. Therefore if the erupting coronal magnetic field were to be allowed to settle back down to a force-free equilibrium, a slow quasi-static evolution could then be considered once again to follow the formation of a second flux rope above the RHS bipole. In future we will consider taking the post-eruption relaxed MHD state and inserting it back into the GNLFFF model to determine if another eruption can be produced at a different location.

An alternative approach to that considered here is to couple the GNLFFF with the MHD simulation at an earlier stage of the evolution, before the flux rope is fully formed. Then, by applying evolving boundary conditions in the MHD simulation identical to the ones of the GNLFF, the formation of the flux rope may be followed in the MHD framework. However, such attempt answers different questions than those addressed in the present paper and poses significant computational problems, relating to flux cancellation boundary condition in the MHD code and computational requirements.

5. Conclusions

In the work presented here, we investigate whether a flux rope slowly formed as result of flux cancellation and photospheric shearing motions in a smoothly varying continuous magnetic field distribution can lead to a fast ejection of plasma from the solar corona. In order to consider the full life span of the flux rope, we have coupled two models: the first is a GNLFFF model that describes the quasi-static formation of the flux rope over time periods of days to weeks and the second is a full MHD simulation that follows the evolution of the system once the flux rope is formed and becomes unstable. In our model no additional shearing or stress of the magnetic field is carried out in the MHD simulation, where the resulting dynamics are solely a result of the stress built up between the two bipoles during the slow evolution.

In our simulation the flux rope initially rises and then interacts with the overlying magnetic field resulting in significant magnetic reconnection above and below. The coronal plasma follows the magnetic field evolution and is expelled out to 2.5  R_⊙. During this process the flux rope is deformed, so in future studies we need to investigate how the flux rope can conserve its shape during the eruption and how the plasma can be further accelerated, since the stress currently accumulated in the GNLFFF model is only enough to produce a modest ejection.

In many descriptions of the onset of CMEs, magnetic reconnection is used to produce the flux rope in a configuration ready to erupt (e.g. Aulanier et al. 2010) and we agree with this vision. In our scenario magnetic reconnection is used to produce the flux rope, but is not necessary in the early stages of the ejection. which is caused by a non-equilibrium. This is significantly different from other scenarios where magnetic reconnection does not play any role, like the Torus instability mechanism (which does not account for the flux rope formation).

We now discuss the role of thermal pressure in the modelling of our eruption. The thermal pressure plays no role in the initiation of the flux rope ejection, but can play a role in the propagation phase. In the present paper we have used a simple atmosphere as we focus more on the coupling of the two models rather than on the atmosphere profile. It is however important to include the role of thermal pressure if we want to move towards a realistic model of a solar eruption from generation to propagation into the interplanetary space. In the present simulations, thermal pressure inhibits the propagation of the flux rope and can quench the ejection in some scenarios (e.g. it obstructs the development of Torus instability Kliem & Török 2006). However, if the thermal pressure decays with height as in a stratified atmosphere its inhibition effect can be neglected. Therefore we will further investigate the effect of thermal pressure and of the atmospheric profile in future studies through included gravity in our model.

Another aspect of this study is to analyse additional magnetic configurations which have been subjected to the accumulation of more stress. In the GNLFFF model a flux rope can grow beyond the size of the flux rope described here. If the eruption is not released, the foot point motion and the flux cancellation can build up a stronger force underneath and the flux rope may then be expelled more violently.

An important achievement of the present work is the development of a technique that successfully couples two different models of the solar corona evolution. We develop such a technique because many different solar phenomena evolve over drastically different time scales. It is currently not possible to consider the full life span of flux ropes produced out of observed surface motions with only one model as computing power is currently insufficient to describe the formation and the ejection processes. In such a framework our efforts allow us to adopt the distinct models, i.e. different assumptions and approximations, during different stages of the coronal active region evolution. This opens up the possibility of many future studies of the Sun’s global magnetic field evolution which would not be possible otherwise.

Acknowledgments

D.H.M. would like to thank STFC, the Leverhulme Trust and the European Commission’s Seventh Framework Programme (FP7/2007−2013) for their financial support. P.P. would like to thank the European Commission’s Seventh Framework Programme (FP7/2007−2013) under grant agreement SWIFF (project 263340, www.swiff.eu) for financial support. These results were obtained in the framework of the projects GOA/2009-009 (KU Leuven), G.0729.11 (FWO-Vlaanderen) and C 90347 (ESA Prodex 9). The research leading to these results has also received funding from the European Commission’s Seventh Framework Programme (FP7/2007-2013) under the grant agreements SOLSPANET (project n 269299, www.solspanet.eu), SPACECAST (project n 262468, fp7-spacecast.eu), eHeroes (project n 284461, www.eheroes.eu). The computational work for this paper was carried out on the joint STFC and SFC (SRIF) funded cluster at the University of St Andrews (Scotland, UK).

References

Appendix A: Interpolation of B from the GNLFFF code to ARMVAC

To import B from the GNLFFF model to the MHD code, we perform a 3D second order interpolation in spherical coordinates. In Fig. A.1 we see an illustration of the r-φ plane on a flat surface. The black grid, represents the GNLFFF grid, where individual magnetic field components (B^GNLFFF) are defined at cell ribs. The red grid represents the MHD grid, where all components of B^VAC are defined at the cell centres. The two grids have a different number of points and different spacing: the GNLFFF grid is composed of 106 × 172 × 181 cells, where the cell size is non-uniform (see Mackay & van Ballegooijen 2006a). In contrast the AMRVAC grid is composed of 132 × 128 × 128 cells uniformly spaced. As the black and red cells in Fig. A.1 are of a different size they do not directly match. Due to this, the relative position and values of the components of B^GNLFFF and B^VAC cannot be determined by simple symmetry conditions. To interpolate from the GNLFFF code to the MHD code, we identify the cell of the GNLFFF grid into which each point of the MHD grid falls. We then use the values of the adjacent GNLFFF cells for the computation of the interpolated magnetic field:

(A.1)where k_r,j_θ,i_φ are the indexes identifying the position in the MHD grid and are the values of the components of B^GNLFFF closest to the coordinates of the point (k_r,j_θ,i_φ) and Δs = (dr,rdθ,rsinθdφ). The adjacent cells are used to compute the gradients of the components of the magnetic field, and the Hessian matrix | H^GNLFFF(r₀,θ₀,φ₀) | through finite differences. The second order interpolation is performed as it guarantees a correct transporting of the spatial derivatives of B between the codes and thus the Lorentz Force.

In Fig. A.2 we compare the three components of B in the GNLFFF model and the MHD code along the φ direction for a typical case. The profiles of each component of the magnetic field in the GNLFFF model (continuous lines) and the magnetic field interpolated in the MHD code (crosses) match well. At the lower boundary, where the gradients of the magnetic field components are the steepest, the interpolated magnetic field departs by at most 1% from the original magnetic field.

Appendix B: Equilibrium test

To check the validity of our approach, before we consider an erupting system we perform a test where a magnetic field configuration in force balance (i.e. (∇ × B) × B = 0) is transported from the GNLFFF model into ARMVAC. Through this we verify that the Lorentz force is properly transported between the codes and that the MHD numerical modelling does not generate artificial forces. This test is performed using an initial condition similar to the case we wish to investigate, consisting of two magnetic bipoles that lie close to one-other but are not overlapping (see Figs. B.2a and b). The magnetic field configuration used is a potential magnetic feld computed from the photospheric magnetic field at Day 21 in Mackay & van Ballegooijen (2006a). The B₀ magnetic field is also a potential magnetic field but computed from the photospheric magnetic field at Day 19. Therefore, in such set up, although both B and B₀ are potential magnetic fields, they are not identical, as they correspond to a different photosphetic magnetic field. This results in a non-zero B₁ that is allowed to evolve in time, according to the MHD equations presented in Sect. 2.3.

As the magnetic field is potential, we expect this configuration to be in stable equilibrium and we should not observe any significant departure from the initial distribution of magnetic field, density, and pressure. Also the assumptions we make about the spatial distribution of ρ and p should not introduce any dynamics. Our assumption in Eq. (1) and the typical values of parameters we use lead to a uniform Alfven speed of ~10⁸ cm/s at t = 0 s. With these parameters we obtain an atmospheric profile in the radial direction for the density and plasma beta, β = 8πp/B² taken from the centre point of the bipoles, as shown in Fig. B.1. The plasma beta varies from 10^-3 at the photosphere to over 10³ at the top of the box, while the coronal density drops by over five orders of magnitude.

With these values the Alfvén crossing time, i.e. the time for an Alfvén wave to travel from one side of the domain to the other, is about 1000 s. While this is one of the time scales that can be determined, for our study a better time scale is the time it takes an Alfvén wave to travel across a single bipole. This Alfvén travel time is around τ_Alf = 120 s so we define this as the characteristic time in the simulation.

In order to test the stability of the system we let the system evolve for a large number of Alfvén times. Figure B.2 shows the initial configuration of the magnetic field at t = 0 (Figs. B.2a and b) and after t ~ 140 τ_Alf (Figs. B.2c and d), where the green lines represent magnetic field lines plotted from the same starting point at both times. Over the 80 Alfvén times, the magnetic field does not show any significant evolution. At the end of the simulation, changes to the magnetic field are of the order of 1% of the initial value. Thus importing the magnetic field from the GNLFFF code has not generated any additional forces.

Note that the plasma density distribution does not show any signifcation evolution, even though small scale changes occur. At most, the bulk velocity is of the order of 10⁴ cm/s, but only appear in regions external to the bipoles, where the density is very low. This shows that the transporting between the GNLFFF and MHD models is possible and that the coupling of the two codes is physically consistent and free from numerical artefacts.

Online material

Movie of Fig.7ad (GIF 3.50 MB) (Access here)

Movie of Fig.7eh (GIF 4.69 MB) (Access here)

All Figures