© ESO, 2013

1. Introduction

Coronal mass ejections (CMEs) are the main drivers of space weather, a term used to describe the effect of plasmas and magnetic fields on the near Earth environment. Although the precise mechanism that causes a CME on the Sun is still unclear, the ejection of a magnetic flux rope from the solar corona successfully describes many of the general features of CMEs (Forbes & Isenberg 1991; Amari et al. 2000; Fan & Gibson 2007). After the ejection of the magnetic flux rope, the CME propagates through the solar corona and into interplanetary space. Understanding CME propagation is a key element in space weather.

A key characteristic of a CME is the three-part structure (Hundhausen 1987). A CME is normally composed of a dense bow front, followed by a low-density region and, finally, a very dense core placed approximately at the centre of the curved front. The propagation of this structure is normally used to infer the trajectory and speed of the CME. Several studies have analysed the kinematics of CMEs, and quoted speeds span from 100 km s^-1 to 3300 km s^-1, with an average speed of about 500 km s^-1 (Gopalswamy 2004). Some CMEs undergo an impulsive acceleration phase followed by a propagation at constant speed (Zhang et al. 2001), while others are subject to relatively long acceleration (Chen & Krall 2003). While CMEs may propagate at different speeds, all CMEs eventually couple with the solar wind once they reach a height of about 4  R_⊙ (Gopalswamy et al. 2000). Similarly, CMEs expand while travelling outward (Ciaravella et al. 2001; Lee et al. 2009).

As CMEs propagate in the solar corona they interact with the pre-exsiting plasma and magnetic fields. The solar corona is a highly dynamic environment composed of magnetized plasma whose global physical evolution can be described by magnetohydrodynamics (MHD) theory. In the solar corona, plasma motions are primarily driven by the Lorentz force. However, under equilibrium conditions we can assume a zero Lorentz force acting on the plasma. In such cases, the plasma distribution in the corona can be described as being stratified by the effect of gravity. The solar corona is inherently a multi-temperature environment, although highly efficient thermal conduction tends to thermalize the plasma, especially in the quiet corona.

The ratio between the thermal and magnetic pressure, the β of the plasma determines whether the local dynamics is governed by the magnetic field (β < 1) or the plasma (β > 1). Gary (2001) describes the β distribution in the corona and shows that in the lower corona, β < 1. In contrast in the outer corona, the β value has a more diversified distribution where β < 1 regions alternate with regions where β > 1. Moreover, when a CME occurs, it carries outward plasma and magnetic flux, which compress the plasma and magnetic field ahead of it, thereby changing the β profile. Despite this interaction, it is important to note that the closed magnetic field of the ejected plasmoid allows neither the mixing nor thermalization of the ejected plasma with the surrounding corona (Ciaravella et al. 2001; Pagano et al. 2007).

Thermal pressure is generally considered to play no role in the initial stages of the ejection of the flux rope, but does become relevant during the propagation phase of the CME when the flux rope travels at large radial distances. Similarly, gravity tends to obstruct the ejection, although it can significantly affect only fragments of the eruption. Innes et al. (2012) observed fragments of the eruption falling back to the Sun after an eruption.

Here we specifically address the role of gravitational stratification on the CME propagation. In our framework, the ejection is caused by an magnetic field configuration where the Lorentz force is not zero. This magnetic field is added to the stratified solar corona, which is initially in hydrostatic equilibrium and decoupled from the magnetic field. In particular, we start from the model of Pagano et al. (2013) where an eruptive magnetic configuration is used, and the full life span of a magnetic flux rope is described by coupling the global non-linear force free model of Mackay & van Ballegooijen (2006) with MHD simulations run with the AMRVAC code (Keppens et al. 2012). In Pagano et al. (2013) a rather simple background density and thermal pressure profile allowed a detailed study of the dynamics of the ejection. Here, we extend that model to include solar gravity and density stratification in the initial conditions. We explore how the parameter space affects the eruption characteristics by tuning both the temperature of the solar corona and the intensity of the magnetic field.

Several studies have included density stratification and gravity in CME modelling and propagation. Pagano et al. (2007) study the role of the ambient magnetic field topology in the CME expansion and thermal insulation, Zuccarello et al. (2009) focus on the CME initiation mechanisms including gravity and the solar wind. Archontis et al. (2009) describe a flux rope ejection in the presence of gravity following a magnetic flux emergence event. Finally, Manchester et al. (2004) model a CME considering the solar wind and gravity up to 1 AU, while Reeves et al. (2010) simulate a solar eruption in 2.5D including solar gravity and the effects of thermal conduction, radiation, and coronal heating. Moreover, Roussev et al. (2012) study solar eruptions with a treatment of thermodynamics and an eruption initiated from below the corona, while the effect of solar wind heating has recently been studied by Pomoell & Vainio (2012).

Although many observations have focussed on studying the properties of CMEs, studies that simultaneously consider the density stratification, and the ejected plasma and magnetic flux are difficult, owing the complexity of the different diagnostics involved. Some recent work that describes the propagation of CMEs in the solar corona include Gopalswamy et al. (2012) and Bemporad et al. (2011). Other studies give diverse examples of how the CME propagation can be modified according to the environment where it develops; e.g., Temmer et al. (2012) describe the interaction of a CME with another CME that preceded it. Separately, some studies have been carried out to analyse the consequence of the gravitational stratification of the solar corona (Guhathakurta et al. 1999; Antonucci et al. 2004; Verma et al. 2013). As observational studies are difficult to carry out, our simulations will shed some light on physical problems that are still difficult to investigate from an observational view point.

The paper is structured as follows. In Sect. 2, we describe the numerical model, in Sect. 3 we describe the results of our simulations, in Sect. 4 we discuss the results, and the conclusions are drawn in Sect. 5.

2. Model

To study the effect of gravitational stratification on the propagation of a CME, we start from the work of Pagano et al. (2013) where a magnetic configuration has been shown to be suitable for ejecting a magnetic flux rope. In the present work, we perform a number of changes to the modelling technique of Pagano et al. (2013), in order to focus on the role of gravitational stratification. We have carried out a set of MHD simulations that consider variations in the temperature of the background coronal atmosphere and the magnetic field intensity.

2.1. MHD simulations

We used the AMRVAC code developed at the KU Leuven to run the simulations (Keppens et al. 2012). The code solves the MHD equations, and the terms that account for gravity are included:

where t is the time, ρ the density, v the velocity, B the magnetic field, and p the thermal pressure. The total energy density e is given by (5)where γ = 5/3 denotes the ratio of specific heats. The expression for the solar gravitational acceleration is (6)where G is the gravitational constant, and M_⊙ denotes the mass of the Sun. In order to gain accuracy in the description of the thermal pressure, we make use of the magnetic field splitting technique (Powell et al. 1999), as explained in Sect. 2.3 of Pagano et al. (2013).

The initial magnetic field condition of the simulations is chosen such that it initially produces the ejection of a flux rope and so the subsequent CME propagation can be studied. The configuration of the magnetic field is Day 19 in the simulation of Mackay & van Ballegooijen (2006). Pagano et al. (2013) in Sect. 2.2 explain in detail how the magnetic field distribution is imported from the global non-linear force free field (GNLFFF) model of Mackay & van Ballegooijen (2006) to our MHD simulations. Since the GNLFFF model contains no plasma, it needs to be defined. We assume an initial atmosphere of a constant temperature corona stratified by solar gravity, where we also allow for low background pressure and density: where we use ρ₀ to tune the density at the lower boundary, μ = 1.31 is the average particle mass in the solar corona, m_p is the proton mass, k_b Boltzmann constant, and T the temperature of the corona. We tune ρ₀ depending on the temperature T in order to always have the same density at the bottom of our computational domain. When ρ_bg = p_bg = 0, we have a stratified atmosphere in hydrostatic equilibrium. In some simulations, we need to have non-zero values for ρ_bg and p_bg to avoid extremely low β values at the outer radial boundary of the simulation. These values are chosen such that p_bg = ρ_bg/(μm_p)k_bT. In our simulations such a departure from equilibrium implies a negligible inflow of plasma (in contrast to the real solar corona that is outflowing). More details are given in the appendix. Finally, we set v = 0 as initial condition for the velocity.

Since the results of Mackay & van Ballegooijen (2006) are specified in terms of the potential vector A, we can uniformly multiply our initial magnetic field distribution without hindering the validity of the results of Mackay & van Ballegooijen (2006). In the present work, we use the maximum value of the magnetic field intensity of the initial magnetic field configuration, B_max, as a simulation parameter. In our simulation, the magnetic field intensity is maxmium, | B| = B_max, at the centre of the right-hand side bipole on the lower boundary.

The simulation domain extends over 3 R_⊙ in the radial dimension starting from r   =   R_⊙. The colatitude, θ, spans from θ   =   30° to θ = 100° and the longitude, φ, spans over 90°. This domain extends to a larger radial distance than the domain used in Mackay & van Ballegooijen (2006) from which we import the magnetic configuration. To define the MHD quantities in the portion of the domain from 2.5 R_⊙ to 4 R_⊙, we use Eqs. (7) and (8) for density and thermal pressure and the magnetic field for r > 2.5 R_⊙ is assumed to be purely radial (B_θ = B_φ = 0) where the magnetic flux is assumed to be conserved (9)It should be noted that the initiation of the ejection is not affected in any way by the extension of the magnetic field, since the initial dynamics are a result of the flux rope that lies at about 1.2 R_⊙ out of magnetohydrostatic equilibrium.

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. The θ boundary condition is designed to not allow any plasma or magnetic flux through, while the φ boundary conditions allow the plasma and magnetic field to freely evolve across the boundaries. These boundary conditions match those used in Mackay & van Ballegooijen (2006). In our simulations, the expanding and propagating flux rope only interacts with the θ and φ boundaries near the end of the simulations, thus they do not affect our main results regarding the initiation and propagation of the CME. At the lower boundary, we impose constant boundary conditions taken from the first four θ-φ planes of cells derived from the GNLFFF model. The computational domain is composed of 256 × 128 × 128 cells distributed in a uniform grid. Full details of the grid can be found in Sect. 2.3 of Pagano et al. (2013).

2.2. Parameter space investigation

To analyse the role of the background stratified corona we ran a set of nine simulations using three different temperatures (T   =   1.5, 2, 3 MK) of the corona and three different maximum intensity values of the magnetic field (B_max = 7, 21, 42 G), as illustrated in Fig. 1.

Fig. 1 Grid summarizing the parameter space we investigate. We ran 9 simulations with all the combinations of T = 1.5, 2, 3 MK with Bmax = 7, 21, 42 G. The numbers in the cell represent the log10(β) at the lower boundary (lower number) and the log10(β) at the upper boundary (higher number). Red are negative values, and blue positive ones.

In our model, a higher corona temperature implies a more uniform density and pressure gradient from the photosphere to the outer corona and a higher amount of mass that constitutes the solar corona, as in Fig. 2a. At the same time, the higher temperature leads to remarkably higher β in the outer corona, as shown in Fig. 2b for the B_max = 7 G case, while in the lower corona (r < 1.2  R_⊙) β is clearly below ~10^-1 regardless of the temperature. For the B_max = 21 and 42 G cases, the β value is lower. The outer corona switches from a low to a high β regime when the temperature increases from 1.5 to 3 MK (Fig. 2b).

We use the temperature to define a part of the parameter space since it is the appropriate parameter to consistently tune the profile of density and thermal pressure in the solar corona. Since we always assume the same value for the photospheric density, an increase in temperature implies a heavier column of mass placed above the magnetic flux rope. By computing the integral of Eq. (7) from r = 1.12 R_⊙ to r = 4 R_⊙, we get a column of mass of 7.5 × 10^-6 g/cm² above the flux rope for the simulation with T = 1.5 MK and 25 times more for the simulation with T = 3 MK.

At the same time, the pressure scale height reduces when lower temperatures are considered. This implies that the ejected flux rope encounters less resistance from the compression produced ahead when propagating. The pressure length scale is around 0.06 R_⊙ when the temperature is T = 1.5 MK, which doubles to 0.12 R_⊙ when the temperature is T = 3 MK. It should be clarified, however, that in our simplified model for the gravitational stratification, we do not aim to realistically describe the entire multi-temperature structure of the solar corona, but rather the corona surrounding an active region. In all of the simulations, the density and pressure drop steeply in comparison to the full extent of the computational domain out to 4 R_⊙.

Fig. 2 a) Profiles of log10(ρ) along r above the centre of the LHS bipole for different temperatures of the solar corona, T. b) Profiles of log10(β) along r above the centre of the LHS bipole for different temperatures of the solar corona, T, and with Bmax = 7 G. c) Profiles of log10(β) along r above the centre of the LHS bipole with different values of the parameter Bmax and with T = 2 MK.

Simulations with different B_max (Fig. 2c) basically differ in the plasma β, which uniformly decreases as B_max increases. Thus, by changing the parameter B_max we generally modify the dominant forces and subsequently the capacity of the solar corona to react to the flux rope ejection.

In some of the simulations (T = 1.5 MK; B_max = 21 and B_max = 42 G), we use a non-zero p_bg and ρ_bg to avoid numerical problems due to extremely low β at the external boundary. For these simulations, the background pressure and density are four orders of magnitude slower than the values near the lower boundary where the dynamics originate. The resulting pressure and density profile departs from the analytic hydrostatic equilibrium profile, but we point out that this can slow down the ejection only slightly, because more plasma has to be displaced as the CME propagates due to the background density (ρ_bg). However, the kinetic energy carried by the ejected plasma in the flux rope is much greater than the kinetic energy produced by these motions, and these effects can in no way produce artificial ejections. The tests described in the appendix show that such departures from hydrostatic equilibrium only lead to appreciable effects on the magnetic configuration on much larger time scales than the one for the CME to escape the solar corona in the present simulations.

3. Results

The initial magnetic field configuration is identical in all of the simulations and it is chosen to produce the ejection of the magnetic flux rope due to an initial excess of the Lorentz force. Figure 3 shows the initial magnetic configuration common to all the simulations. The only difference from the one used in Pagano et al. (2013) is the larger extension of the domain in the radial direction. The flux rope that is about to erupt lies under the arcade system, and external magnetic field lines are shown above. Some of the external magnetic field lines belong to the external arcade. while others are open. A full description of the initial magnetic field configuration is given in Sect. 3.1 of Pagano et al. (2013).

Fig. 3 Magnetic field configuration used as the initial condition in all the MHD simulations. Red lines represent the flux rope, blue lines the arcades, green lines the external magnetic field. The lower boundary is coloured according to the polarity of the magnetic field from blue (negative) to red (positive) in arbitrary units.

In Pagano et al. (2013) the dynamics of the ejection and the production of the initial condition in the GNLFFF simulation is discussed in detail. Here, we only focus on the propagation of the flux rope once the ejection is triggered.

We first describe the characteristics of a typical ejection in this framework. Following this we compare some key features between different simulations in order to highlight the role of the temperature of the stratified corona (T) and the intensity of the magnetic field (B_max).

3.1. Typical simulation, (T =  2 MK, B_max = 21 G)

A typical simulation is one with T = 2 MK and B_max = 21 G, which is the central one in the parameter space grid shown in Fig. 1. In this simulation the flux rope escapes the computational domain at 4 R_⊙ and a CME occurs. In Fig. 4 we show the ejection and expansion of the flux rope in 3D to illustrate the simulation with some of the magnetic field lines drawn from the flux rope footpoints and from the axis of the flux rope at different times. In future figures we draw cuts in 2D planes to focus more clearly on specific aspects.

Fig. 4 The ejection of the magnetic flux rope. The red lines are all magnetic flux rope lines drawn from both the footpoints of the flux rope and from the centre of its axis.

In the initial magnetic configuration (Fig. 3), the flux rope lies above the centre of the left-hand side (LHS) bipole, three magnetic arcades connect adjacent and opposite polarities, and a larger arcade connects the opposite external polarities of the two bipoles. The flux rope that is in non-equilibirium lies in the asymmetric part of the global configuration. The density profile falls off with radial distance where ρ ~ 10^-15 g/cm³ at the lower boundary and ρ ~ 10^-20 g/cm³ at the top boundary (Fig. 2a). The plasma β is approximately 10^-2 near the flux rope and it increases to 10^-1 radially above it (Fig. 2c). We note that there are regions where β reaches as much as 10² in the corona (such as near null points), where the magnetic field is weak.

To follow the evolution of the simulation, in Fig. 5 we show two quantities in the r−φ plane passing through the centre of the bipoles. Figures 5a–c show the density contrast, (10)that indicates the variation in density with respect to the initial condition, which is useful for marking the ejection of plasma. Figures 5d–5f show the ratio B_θ/|B | that illustrates the evolution of the axial magnetic field of the flux rope, which initially lies in the θ-direction. As explained in Mackay & van Ballegooijen (2006), the formation of the flux rope due to magnetic diffusion leads to a strong axial component of the magnetic field along the PIL of the LHS bipole. Simultaneously, the increased tilt of the bipole leads to a quasi-antiparallel magnetic field around the magnetic flux rope. Because of the construction of the system of a E-W bipole and N-S PIL, at t = 0 the axial magnetic field of the flux rope is mostly along the θ direction and the blue region between pink and green in Figs. 5d–f, across which B_θ/|B | changes sign, marks the borders of the magnetic flux rope and the overlying arcade. As long as no major reconnection occurs, the quantity B_θ/|B | is clearly positive above the bipoles where the magnetic flux rope propagates and expands.

Fig. 5 Simulation with T = 2 MK and Bmax = 21 G. a)–c) Maps of density contrast (Eq. (10)) in the (r−φ) plane passing through the centre of the bipoles at t = 0 h, t = 0.58 h, t = 1.74 h. Superimposed are magnetic field lines plotted from the same starting points (green lines) and the contour line of β = 1 (white line). d)–f) Maps of Bθ/|B | on the same plane and at the same times. Maps show the full domain of our simulations from r = 1  R⊙ to r = 4  R⊙. In panel d) the yellow dashed line is the cut for the plots in Fig. 6. The temporal evolution is available in the on-line edition.

In this simulation, the flux rope is ejected outwards, and it leads to an increase in density at larger radii (Figs. 5a–c), and the propagation and expansion of the region where the magnetic field is mostly axial, i.e. B_θ/|B| ~ 1 (Figs. 5c–f). The shape of the density propagation roughly reproduces the typical three-component structure of a CME, which is clearly visible in Fig. 5b. A higher density bow front propagates upwards, and behind it lies a region with lower density. Finally, a dense core is located at the centre of the expanding dome. In our simulation, the region with highest density coincides with the region where B_θ/|B | is positive. The quantity B_θ/|B | is maximum just ahead of the high- density region, while the peak of density is highest roughly at the centre of the expanding structure where B_θ/|B | remains positive (Fig. 5b). This suggests that the magnetic flux rope is propagating upwards, perpendicular to its axial magnetic field lines, lifting coronal plasma to produce the high density region. The same process is seen in Pagano et al. (2013), and it is considered a standard process in filament eruptions.

At the same time, the region where a positive axial component of the magnetic field is dominant expands and propagates upwards, roughly covering the same volume of space as that of the high-density region (white and pink region in Figs. 5d–f). Since the region where B_θ/|B| > 0 roughly reproduces the high- density region that corresponds to the CME, we use this quantity to track the ejection and expansion of the CME. Figure 6 shows the profile of B_θ/|B | radially from the centre of the LHS bipole at different times. The radial position where B_θ/|B| = 0 along the radial direction vertically from the centre of the LHS bipole is defined to be the top of the magnetic flux rope in our representation.

In the plot in Fig. 6 we see the top of the flux rope at 1.4 R_⊙ at t = 0 h (black line), and it reaches 3.15 R_⊙ after 1.16 h (red line). This clearly shows the ejection of the flux rope.

Fig. 6 Profile of Bθ/|B | above the centre of the LHS bipole, at different times (different colours). Dashed lines of different colours indicate where Bθ/|B| = 0 and where we locate the top of the flux rope at a given time.

3.2. Parameter space investigation

3.2.1. Quenched ejections

Following the location of the top of the flux rope we are able to track the ejection of the flux rope and thus identify whether the flux rope ejection produces a CME (i.e. extends beyond r   =   4 R_⊙) or whether it is a quenched ejection. We find that the simulations with T = 2 MK and T = 3 MK where B_max = 7 G do not show any CME propagation. The flux rope only rises up to a given height and then stops. In Fig. 7, we show a quenched ejection for the simulation with T = 2 MK and B_max = 7 G. In this simulation the magnetic flux rope is stopped and remains positioned at r = 1.8 R_⊙ after 1.74 h.

A key difference between the simulations where a CME develops (Fig. 5a) and where the ejection is quenched (Fig. 7a) is that a high-β region lies above the flux rope in the simulation where the ejection is quenched. The fast CME can only develop in a low-β environment. In particular, the β = 1 contour overlies the flux rope in Fig. 7a, while it only surrounds the null-point region in Fig. 5a. Although this is not a necessary or sufficient condition to determine the onset of a full or quenched ejection, whether a high-β region lies over the flux rope is one of the factors that can determine the CME evolution.

Fig. 7 Simulation with T = 2 MK and Bmax = 7 G. Maps of Bθ/|B | in the (r−φ) plane passing through the centre of the bipoles at t = 0 h and t = 1.74 h. Superimposed is the contour line of β = 1 (white line). Black crosses indicate where the top of the flux rope is positioned.

3.2.2. Parameter space of CMEs

All the remaining simulations show an ejection of the flux rope in a similar manner to the one in Sect. 3.1. We do not present a detailed analysis here for each simulation, but plot the position of the top of the flux rope as a function of time. In Fig. 8a we show the position of the top of the flux rope for the three simulations where B_max = 21 G (central column in Fig. 1). We stop tracking the top of the flux rope as soon as it is close to the outer boundary and its propagation is affected by boundary effects. The speed quoted in Fig. 8 is the average speed of propagation. The average speed of the CME spans a wide range from 69 km s^-1 to 498 km s^-1. The lower the temperature, the faster the CME. Also, the CMEs propagating in a 2 MK or 1.5 MK corona are at near constant speed, while the speed varies slightly with radial distance for the slowest CME (T = 3 MK).

Figure 8b shows the height of the top of the flux rope as a function of time for the three simulations with T = 2 MK (central row in Fig. 1). The higher the magnetic field (B_max), the stronger the initial force due to the unbalanced Lorentz force under the magnetic flux rope, and the CME travels faster. With B_max = 7 G, we have a quenched ejection, while we have a 718 km s^-1 fast CME with B_max = 42 G. By changing the parameter B_max, we change the Alfvén speed in the system and thus the evolution time scale. It should be noted that the result of the simulations with different B_max cannot be inferred solely through timescale arguments owing to the presence of the plasma. For instance, this is demonstrated by the qualitatively different evolution of the simulation with B_max = 7 G from the other simulations.

Fig. 8 a) Position of the top of the flux rope as a function of time in the three simulations with Bmax = 21 G. b) Position of the top of the flux rope as function of time in the three simulations with T = 2 MK.

We note that the speed of the top of the flux rope is not necessarily the speed of the CME. The top of the flux rope represents the leading edge of the CME, which is subjected to the combined effect of propagation and expansion. However, the results approximately indicate the speed of the CME and its dependence on the parameters T and B_max. In Fig. 9 we summarize the speed of the CMEs for all of the simulations. The fastest CME in our simulations can reproduce speeds of 838 km s^-1, and only two simulations show a quenched ejection.

Fig. 9 Grid summarizing the average speed of the CME in the parameter space we investigate. Crosses indicate a quenched ejection, and in the other boxes we write the value of the CME speed computed from tracking the top of the flux rope as in Fig. 8.

We also use the quantity B_θ/|B | to follow the φ extension of the flux rope and thus its expansion. We draw the contour of the flux rope by tracking where the cut of B_θ/|B | along the φ direction on the r-φ plane passing through the centre of the bipoles changes sign at different radial distances. Then, we consider the maximum φ extension of the resulting flux rope contour. Figure 10a shows the extension of the flux rope as a function of the height of the top of the flux rope for the simulations where unquenched eruptions occur. In this plot we group the simulations by temperature with different colours. The intensity of the magnetic field B_max does not significantly influence the expansion of the CME, i.e. its shape. This can be seen because lines with the same colour lie close to one another in Fig. 10a. In contrast, the background temperature, T, plays a significant role in shaping the ejection.

Below the height of 2 R_⊙ the expansion proceeds in a similar manner for all the simulations. However, around this radial distance the simulations with T = 3 MK show a larger expansion for higher radii. In the simulations with T = 1.5 MK and T = 2 MK, the expansion of the flux rope proceeds nearly uniformly until the flux rope touches the upper boundary where the extension is approximately 50° and 70°, respectively. In contrast, in the simulations with T = 3 MK, the expansion rate increases at about 2 R_⊙, and the flux ropes covers an angle that is close to 80° beyond 3 R_⊙. At the height of 2.5 R_⊙, the simulations with T = 3 MK show an expansion that is about twice as large as the expansion found in the simulations with a cooler stratified solar corona.

In summary, Fig. 10a shows that higher temperatures favour higher expansion of the CME, which in turn results in a flattened shape. The higher pressure and density above the ejecting flux rope slow down the radial propagation and favour the expansion of the CME towards the sides. An example of the different shape of the CME is given in Fig. 10b where we show a snapshot of the simulation with T = 3 MK and B_max = 21 G (to be compared with Fig. 5e).

Fig. 10 a) Angular extension of the flux rope as a function of the radial distance of the top of the flux rope for all the simulations that show a CME. Simulations are grouped according to the simulation parameter T given with different colours. b) Map of Bθ/|B | in the (r−φ) plane passing through the centre of the bipoles at t = 1.55 h for the simulation with T = 3 MK and Bmax = 21 G. The map shows the full domain of our simulations from r = 1 R⊙ to r = 4 R⊙.

4. Discussion

In this paper we have considered the role of gravitational stratification on the propagation of a CME. In order to do so, we investigated the parameter space by tuning the temperature of the corona (i.e. gravity stratification) and the maximum value of the magnetic field (i.e. the entire β of the corona). The techniques and the initial magnetic configurations are identical to those of Pagano et al. (2013) where a CME initiation was successfully modelled.

4.1. Comparison with the flux rope ejection simulated in Pagano et al. (2013)

The purpose of the present work is not to analyse the mechanism that initiates the flux rope ejection in detail, since it does not show any significant difference from what is described in Pagano et al. (2013). However, once the onset of the ejection occurs, some differences appear in the CME propagation depending on the characteristics of the background corona, and in the present paper we try to get some physical insight into the underlying physics involved. In contrast to the background coronal model in Pagano et al. (2013), the plasma pressure and density profile in the present simulations are a consequence of gravitational stratification.

A drop in thermal pressure implies that the ejection encounters less resistance from the solar corona, which should aid the production of a CME. While this is true, it should be noted that the ejected plasma is now subject to the gravity force that was not considered previously. The two effects nearly balance out since the speed of the ejection is generally comparable to the ones measured in Pagano et al. (2013). It should be noted, however, that slightly higher speeds are obtained in the present work. Another aspect is that because of the lower thermal pressure in the outer solar corona, the interaction between the ejected flux rope and the surrounding plasma is reduced, and this leads to a clearer three-component structure of the CME. Our present simulations clearly show the light-bulb structure that was not seen in Pagano et al. (2013). For the same reason, significantly less magnetic reconnection due to numerical diffusion occurs between the flux rope magnetic field and the external arcade magnetic field, resulting in higher coherence of the expelled flux rope and less mixing between the flux-rope magnetic field lines and the external arcade magnetic field lines.

4.2. Role of plasma β on CME speed

In this study, the plasma β plays a major role in determining the evolution of the system. Firstly, the plasma β distribution determines whether an ejection is quenched or if it will escape the solar corona. In simulations where only low β plasma covers the initial flux rope, it quickly reaches a height of 4 R_⊙. In contrast to the two simulations (Fig. 9) where the ejection stops, the flux rope is surrounded by high β plasma (see Fig. 7a). In our study we have used coronal values for magnetic field intensity, density of plasma, and temperature, and our simulations seem to confirm that the onset of solar eruptions is caused by a purely magnetic process, but the evolution of a flux rope ejection does depend on the plasma parameters. This is an important factor to be considered when we wish to address the issue of CME generation, especially in the space weather context.

Secondly, a lower plasma β seems to favour higher CME speed. In our investigation of the parameter space this was done by increasing the intensity of the flux rope magnetic field by a factor of 3 or 6. Such an increase leads to an increase of a factor of 5 or 8 in the CME speed (simulations with T = 1.5 MK in Fig. 9). Similar results are found when comparing the CME speed in the ejecting simulations at T = 2 MK and T = 3 MK,

More generally, Fig. 8b shows that we simulate two main evolution branches: some simulations show a CME and others show a quenched ejection. The plasma β in the outer corona seems to be the parameter that determines into which branch the system evolves (see Fig. 2c). Of course, with a continuous variation of the physical properties of the system, we expect intermediate cases that show a mixture of the two regimes. This occurs when the plasma β ~ 1 in the outer corona. An example of such a regime is the simulation with T = 3 MK and B_max = 21 G where an extremely slow CME (69 km s^-1) does not travel at a constant speed (Fig. 8a). In this simulation, high plasma β lies above the flux rope, but at greater radial distance than in the other simulations where the ejection stops. Similar behaviour was found in Pagano et al. (2013).

4.3. Role of gravitational stratification on CME speed

In the framework of the present study, another way of tuning the plasma β in the outer corona is to change the coronal temperature. In particular, the colder the atmosphere, the lower the plasma β, and the more likely a full ejection occurs. However, tuning the temperature only affects the capacity of the outer corona to react to the ejection, but not the initial force to which the flux rope is subject. For this reason, in our framework the temperature of the corona has a weaker effect compared to the magnetic field intensity in producing CMEs at higher speeds. At the same time, the role of the temperature stratification deserves more detailed analysis in future studies. In fact, it is remarkable that simulations where the coronal temperature only differs by 0.5 MK show such a different behaviour. This could be the reason for the enhancement of some CMEs in the solar corona. Similar to the work of Lin (2004), we find that the gravitational force can be effective on the flux rope ejection when the magnetic field is relatively weak. In particular, it can only slow down the ejected flux rope when the magnetic field is relatively weak. In future, we need to investigate the cause of the flux rope ejections more. In our study, the ejection is caused by the build up of an excess of the Lorentz force; however, the question remains open whether a MHD instability, like the Torus instability (Kliem & Török 2006), could explain the ejection. In our case, a treatment like the one of Kliem & Török (2006) should include the role of thermal pressure and gravity. In doing so, the critical index for the enhancement of the ejection could be determined as could whether it differs from previous studies.

4.4. Role of gravity stratification on CME expansion

Another result of our study is that the gravitational stratification can be responsible for shaping the CME, since the thermal pressure ahead of the CME front can stimulate an expansion at its sides. We identified this effect in particular in the simulation with T = 3 MK and B_max = 21 G where the expanding CME acquires a shape that is significantly different from the ones found in other simulations. However, the result is more general, as shown in Fig. 10b. From the analysis of our simulations, we find that the width of the CME bulb can be significantly different from one simulation to the next if different coronal temperatures are considered.

Further studies are required, but in principle different stratifications may be responsible for the variety of CME shapes observed. In particular, the CME we describe in the simulation with T = 3 MK and B_max = 21 G seems to be very similar to the ejection that occurred on June 13, 2010 as described by Gopalswamy et al. (2012). It is remarkable how the CME shape in our simulation matches the snapshot of the ejection at 5:40:54, where not only is the front profile reproduced, but the northern side of the ejection also seems to behave as if there is another bipole besides the one involved in the ejection as in our simulation (Fig. 10b). It should, however, be noted that in Gopalswamy et al. (2012) the CME is observed at 1.4 R_⊙, significantly lower than our snapshot in Fig. 10b. Also Patsourakos et al. (2010) studied the expansion of a CME and found the expansion to be initially very rapid and characterized by a self-similar evolution at later times. Our results differ from their findings, since in our simulations the CME bubble expansion agrees with a constant decrease in the aspect ratio (radial position/radius), and it never reaches the self-similar regime. Such a difference in our results could be explained by the lack of an outflowing solar wind in our simulations that would contribute to narrowing the propagation of the CME. This feature will be investigated in future studies.

5. Conclusions

The present study aims at understanding the role of gravitational stratification in the propagation of CMEs. To identify the role of gravitational stratification, we ran several simulations that vary two quantities: the stratification temperature (T) and plasma β. In all simulations an identical magnetic field configuration is used that has been proven to be suitable for an ejection.

Our work appeared to produce a reasonable model for the early evolution of actual CMEs, since observations show that the CME acceleration tends to vanish above 4 R_⊙ (Vršnak 2001), when the CME kinematics couple to the solar wind (Gopalswamy et al. 2000). Therefore our model covers the domain of initiation of a CME, and we can reasonably assume that all the ejections that leave our domain would not need any further acceleration mechanisms to travel through interplanetary space dragged by the solar wind.

This study showed that gravitational stratification has an important effect on the propagation of CMEs in the solar corona through the way it specifies how large the plasma β becomes. We also find that the plasma β distribution is a crucial parameter that determines whether a flux rope ejection escapes the solar corona, turning into a CME, or if it just makes the flux rope find a new equilibrium at a greater height. Similarly, we find that a cooler solar corona (T ~ 1.5 MK) can help the escape of the CME and make it travel faster.

Both of these results – first the importance of a low β region above the magnetic flux rope to allow the ejection and secondly the role of the coronal temperature in the CME speed – can be tested with observations where 3D magnetic field reconstructions are carried along with simultaneous density and temperature diagnostics to infer the plasma stratification. However, magnetic field reconstructions are more reliable from disk observations, which makes it very difficult to infer coronal stratification because of line-of-sight integration effects. However, future missions of the STEREO type with this capacity could be required.

Acknowledgments

The authors would like to thank Dr. Bernhard Kliem for useful discussions and suggestions. The authors would also like to thank the referee for the constructive and positive feedback on the manuscript. 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, http://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 No. 269299, http://www.solspanet.eu), SPACECAST (project No. 262468, fp7-spacecast.eu), and eHeroes (project No. 284461, http://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: Equilibrium tests

To check the validity of our approach, we performed some tests with a magnetic field configuration in force balance and a stratified atmosphere. In Pagano et al. (2013) a similar test was performed to verify that the Lorentz force is properly transported from the GNLFFF to the MHD model without generating artificial forces. We did not repeat a similar test, and we consider those results as valid. Instead, we want to check that the stratified atmosphere and the extension of the domain to 4 R_⊙ does not introduce significant spurious forces able to alter the stability of the system.

This test is performed using an initial force-free magnetic condition similar to the case we wish to investigate, consisting of two magnetic bipoles that lie close to one another but do not overlap. We refer to Pagano et al. (2013) for further details on the magnetic configuration.

We consider several stratified density and pressure profiles that cover the parameter space described in Sect. 2.2. Namely, we show here in Fig. A.1 simulations with (T = 3 MK, B_max = 7 G), (T = 2 MK, B_max = 21 G), and (T = 1.5 MK, B_max = 42 G).

Fig. A.1 Maps of log10(ρ [g/cm3]) in the (r−φ) plane passing through the centre of the bipoles at different times for different test simulations. Superimposed are magnetic field lines plotted from the same starting points (green lines) and the contour line of β = 1 (white line). a) and d)T = 3 MK, Bmax = 7 G; b) and e)T = 2 MK, Bmax = 21 G; c) and f)T = 1.5 MK, Bmax = 42 G.

To test the stability, we let the system evolve for a large number of Alfvén times. In the simulation with B_max = 7 G, we estimate the Alfvén time τ_Alf = 120 s as the time it takes an Alfvén wave to travel across a single bipole with a typical speed of 800 km s^-1. Consequently, the Alfvén time is τ_Alf = 40 s and τ_Alf = 20 s for the simulations with B_max = 21 G and B_max = 42 G. Figure A.1 shows the density contrast (defined in Eq. (10)). The plasma density contrast does not show any significant evolution for any of the test simulations, even though small-scale changes occur.

The only test that shows an observable change in the magnetic structure is the one with T = 1.5 MK and B_max = 42 G. This is due to the non-zero values of p_bg and ρ_bg that lead to a downward bulk motion of plasma in the higher part of the corona, as explained in Sect. 2.2. After this occurs, numerical reconnection takes place, and some of the initially open magnetic field lines reconnect near the outer boundary. However,

only a small change in the magnetic configuration can be seen after ~3.3 h, equivalent in this particular simulation to 600 τ_alf. In the corresponding simulation in Sect. 3.2 where a CME occurs, the magnetic flux rope reaches the top of the domain in less than 0.8 h. Therefore the motions due to the non-zero values of p_bg and ρ_bg are negligible in comparison to the dynamics produced by the flux rope ejection.

These tests show that the transporting between the GNLFFF and MHD models is possible when including the effects of gravity, a stratified atmosphere, and magnetic field intensity tuning and that these changes do not generate numerical artefacts in the MHD simulation.

Online material

Movie 1 (Access here)

Movie 2 (Access here)

All Figures

	Fig. 1 Grid summarizing the parameter space we investigate. We ran 9 simulations with all the combinations of T = 1.5, 2, 3 MK with Bmax = 7, 21, 42 G. The numbers in the cell represent the log10(β) at the lower boundary (lower number) and the log10(β) at the upper boundary (higher number). Red are negative values, and blue positive ones.
In the text

	Fig. 2 a) Profiles of log10(ρ) along r above the centre of the LHS bipole for different temperatures of the solar corona, T. b) Profiles of log10(β) along r above the centre of the LHS bipole for different temperatures of the solar corona, T, and with Bmax = 7 G. c) Profiles of log10(β) along r above the centre of the LHS bipole with different values of the parameter Bmax and with T = 2 MK.
In the text

	Fig. 3 Magnetic field configuration used as the initial condition in all the MHD simulations. Red lines represent the flux rope, blue lines the arcades, green lines the external magnetic field. The lower boundary is coloured according to the polarity of the magnetic field from blue (negative) to red (positive) in arbitrary units.
In the text

	Fig. 4 The ejection of the magnetic flux rope. The red lines are all magnetic flux rope lines drawn from both the footpoints of the flux rope and from the centre of its axis.
In the text

Fig. 5 Simulation with T = 2 MK and Bmax = 21 G. a)–c) Maps of density contrast (Eq. (10)) in the (r−φ) plane passing through the centre of the bipoles at t = 0 h, t = 0.58 h, t = 1.74 h. Superimposed are magnetic field lines plotted from the same starting points (green lines) and the contour line of β = 1 (white line). d)–f) Maps of Bθ/|B | on the same plane and at the same times. Maps show the full domain of our simulations from r = 1  R⊙ to r = 4  R⊙. In panel d) the yellow dashed line is the cut for the plots in Fig. 6. The temporal evolution is available in the on-line edition.

In the text

	Fig. 6 Profile of Bθ/|B | above the centre of the LHS bipole, at different times (different colours). Dashed lines of different colours indicate where Bθ/|B| = 0 and where we locate the top of the flux rope at a given time.
In the text

	Fig. 7 Simulation with T = 2 MK and Bmax = 7 G. Maps of Bθ/|B | in the (r−φ) plane passing through the centre of the bipoles at t = 0 h and t = 1.74 h. Superimposed is the contour line of β = 1 (white line). Black crosses indicate where the top of the flux rope is positioned.
In the text

	Fig. 8 a) Position of the top of the flux rope as a function of time in the three simulations with Bmax = 21 G. b) Position of the top of the flux rope as function of time in the three simulations with T = 2 MK.
In the text

	Fig. 9 Grid summarizing the average speed of the CME in the parameter space we investigate. Crosses indicate a quenched ejection, and in the other boxes we write the value of the CME speed computed from tracking the top of the flux rope as in Fig. 8.
In the text

Fig. 10 a) Angular extension of the flux rope as a function of the radial distance of the top of the flux rope for all the simulations that show a CME. Simulations are grouped according to the simulation parameter T given with different colours. b) Map of Bθ/|B | in the (r−φ) plane passing through the centre of the bipoles at t = 1.55 h for the simulation with T = 3 MK and Bmax = 21 G. The map shows the full domain of our simulations from r = 1 R⊙ to r = 4 R⊙.

In the text

	Fig. A.1 Maps of log10(ρ [g/cm3]) in the (r−φ) plane passing through the centre of the bipoles at different times for different test simulations. Superimposed are magnetic field lines plotted from the same starting points (green lines) and the contour line of β = 1 (white line). a) and d)T = 3 MK, Bmax = 7 G; b) and e)T = 2 MK, Bmax = 21 G; c) and f)T = 1.5 MK, Bmax = 42 G.
In the text