INITIATION OF CORONAL MASS EJECTIONS IN A GLOBAL EVOLUTION MODEL

To model the evolution of the large-scale magnetic field in the solar corona, we make the following simplifying assumptions.

• 1.  
  The magnetic field in the solar corona may be separated into a large-scale, mean, component on lengthscales ∼50 Mm, and a fluctuating component on lengthscales ∼1 Mm.
• 2.  
  On the timescale of large-scale photospheric motions, the large-scale coronal magnetic field is in equilibrium, except perhaps during flares or CME eruptions.
• 3.  
  Above the photosphere, and in the low corona up to 2.5 R_☉, magnetic forces dominate over plasma forces (i.e., it is a low-β plasma). In the photosphere magnetic flux tubes are passively advected.
• 4.  
  The large-scale magnetic field in the photosphere is predominantly radial, and its horizontal lengthscale is larger than that of the supergranular convection.
• 5.  
  For purposes of the large-scale magnetic field evolution, emerging active regions may be represented by twisted magnetic bipoles, with an idealized form.

Assumption 1 allows the application of mean field theory to the corona (Seehafer 1994; van Ballegooijen et al. 2000). In that case, only the mean magnetic field is solved for, and the effect of significant small-scale structure—such as braiding and current sheets produced by interaction with convective flows in the photosphere (Parker 1972)—is included through a turbulent electromotive force (e.m.f.). The mean-field induction equation takes the form

$$\begin{equation} \frac{\partial \mathbf {A_0}}{\partial t} = \mathbf {v}_0\times \mathbf {B}_0 + \mathcal {E}_0, \end{equation} \tag{ 1 }$$

where A₀(r, t) is the vector potential for the mean magnetic field B₀ = ∇ × A₀, v₀(r, t) is the mean plasma velocity, and $\mathcal {E}_0(\mathbf {r},t)$ is the turbulent e.m.f..

Assumption 2 is valid because magnetic disturbances from the photosphere propagate up into the corona at Alfvén speeds of the order 1000 km s⁻¹ or greater (Régnier et al. 2008), well above the maximum flow speed of surface rotation (of the order 0.2 km s⁻¹). The evolution of the coronal magnetic field may then be understood as a continual distortion of existing equilibrium by footpoint motions, with subsequent fast relaxation to a neighboring equilibrium (Seehafer 1994). By assumption 3, the form of this equilibrium is well approximated by a force-free field j₀ × B₀ = 0, where j₀ = ∇ × B₀/μ₀ is the mean current density.

In contrast to the corona, the photospheric field is assumed to be passively advected by plasma flows. Assumption 4 then justifies use of the standard surface flux transport model for the evolution of the radial component B_0r on the solar surface (Sheeley 2005, and references therein). That the large-scale photospheric field should be approximately radial is supported both by vector magnetic measurements (Martinez Pillet et al. 1997), and by theoretical considerations, due to the concentration of convection zone magnetic flux into isolated kilogauss flux tubes (van Ballegooijen & Mackay 2007). In the surface flux transport model, B_0r is advected by the large-scale motions of differential rotation and meridional circulation, and the influence of supergranular convection on this large-scale field is modeled by a surface diffusion term (Leighton 1964). In the model described in this paper, surface flux transport, coupled with the emergence of new active regions, acts as the lower boundary condition to drive the coronal field evolution.

Assumption 5 is the most uncertain, due to the idealized nature of the magnetic bipoles. Previous studies have shown that such idealized bipoles are adequate for reproducing the distribution of large-scale radial field B_0r on the solar surface (Paper I). The main uncertainty is the inclusion of twist. While our implementation in this paper is simplistic, it is better than not including any twist. Such twist has been shown to be required both by vector magnetic field measurements (Démoulin 2007) and by our previous simulations (Paper II). Part of the aim of the present work is to determine whether such a simplified model can yield useful information about the global magnetic field, or whether more detailed modeling of individual regions is necessary to capture even the simplest large-scale phenomena.

As in Paper II, we solve Equation (1) in a domain extending from 0° to 360° in longitude, −80° to 80° in latitude, and R_☉ to 2.5 R_☉ in radius. We assume the turbulent e.m.f. to take the form

$$\begin{equation} \mathcal {E}_0 = -\eta \mathbf {j}_0, \end{equation} \tag{ 2 }$$

with the turbulent diffusivity

$$\begin{equation} \eta = \eta _0\left(1 + 0.2\frac{|\mathbf {j}_0|}{|\mathbf {B}_0|}\right) \end{equation} \tag{ 3 }$$

as in Mackay & van Ballegooijen (2006a). The first term is a uniform background value and the second term is an enhancement in regions of strong current density. This enhancement was introduced because earlier simulations (Mackay & van Ballegooijen 2001) produced highly twisted structures which are typically not observed (Bobra et al. 2008; Su et al. 2009).

To follow the evolution through a sequence of nonlinear force-free states in response to flux emergence and surface shearing (Section 2.1), the MHD momentum equation is approximated by the magnetofrictional method (Yang et al. 1986), setting

$$\begin{equation} \mathbf {v}_0 = \frac{1}{\nu }\frac{\mathbf {j}_0\times \mathbf {B}_0}{B^2} + v_{\rm out}(r)\hat{\mathbf {r}}. \end{equation} \tag{ 4 }$$

The second term in Equation (4) is a radial outflow imposed near the upper boundary, where it simulates the effect of the solar wind in opening up field lines in the radial direction (Mackay & van Ballegooijen 2006a). It has a peak value v₀ = 100 km s⁻¹ and falls off exponentially away from the upper boundary, so that only structures near this boundary are affected.

We assume closed boundaries in the latitudinal direction and an open boundary at the upper source surface (r = 2.5 R_☉), where the radial outflow ensures that the magnetic field is radial. At the lower boundary (r = R_☉) the radial magnetic field B_0r is evolved according to the surface flux transport model with flux emergence, as described in Paper I.

In this paper, we simulate the solar corona over 177 days in 1999, the same period as Paper II. The initial condition is a potential field extrapolation for 1999 April 16 (day of year 106) taken from a Kitt Peak synoptic magnetogram corrected for differential rotation (Paper I). The coronal magnetic field is then evolved continuously using Equations (1) and (4), with surface flux transport on the lower boundary and flux emergence, until 1999 October 10 (day of year 283). During this period, 119 new active regions are inserted in the form of 3D twisted magnetic bipoles, to be described in Section 2.3. The model then follows the quasi-static buildup of magnetic energy in the solar corona.

The emerging magnetic bipoles, which constitute the primary observational input to our simulations, are based on active regions observed in synoptic normal-component magnetograms from the National Solar Observatory (NSO) at Kitt Peak. For each of 119 bipoles emerged during this 177 day period, we recorded the day of central meridian passage, latitude, size, tilt angle, and magnetic flux from the Kitt Peak data (Paper I). Rather than using the actual magnetic field distribution of each region, we insert idealized magnetic bipoles—of the mathematical form given in Paper II—into the simulation, with parameters chosen to match the above observed properties. The advantage of this approach is that the bipoles may readily be inserted in 3D with photospheric and coronal components, and nonzero magnetic helicity. To model each observed active region more accurately would require the availability of vector magnetograms, and an extrapolation to model the unknown coronal magnetic field. Methods are under development to extrapolate nonlinear force-free fields from photospheric vector magnetograms, but there are still serious problems with this procedure (see DeRosa et al. 2009), apart from the major computational cost for a large number of regions. Our idealized bipoles cannot reproduce the detailed structure and dynamics within individual active regions. However, the model has been shown to reproduce both the long-term global dispersal of active regions across the solar surface (Paper I) and the statistical hemispheric pattern of filament chirality (Paper II). The reproduction of the observed surface magnetic field is illustrated in Figure 1, which shows the observed and simulated B_0r(R_☉, θ, ϕ) on day 147 (after 41 days' evolution) and on day 283 (at the end of the simulation). The accuracy of the simulated field may be seen from Figures 1(c) and (d), where the zero contours from the observed magnetograms are overlaid on the simulated flux distribution. We now consider two aspects of the idealized bipoles which are poorly determined from observations.

Zoom In Zoom Out Reset image size

2.3.1. Bipole Insertion

2.3.2. Bipole Twist

The 3D magnetic bipoles are given a nonzero magnetic helicity through a parameter β (defined in Equations (6) and (8) of Paper II), which describes the initial twist of the field lines in the corona, and does not affect the bipole's radial magnetic footprint on the solar surface. Figure 2 shows the structure of a bipole with two different magnitudes and signs of β, taken from two simulation runs in this paper. The effect of β on the flux ropes forming in a simple two-bipole configuration has previously been considered by Mackay & van Ballegooijen (2006a). Here coronal flux ropes formed not only within the two bipoles at sites of flux emergence, but also between the bipoles where flux was canceling. The magnitude of β affected the formation rate of magnetic flux ropes, and hence the time until loss of equilibrium, but did not affect whether or not a flux rope formed. Mackay & van Ballegooijen (2005) showed that the sign of β also affects the chirality type (direction of shear) in this simple configuration. Subsequent studies of the chirality generated in global simulations have found that both emerging helicity and shearing by surface motions can be important at different locations on the Sun (Paper III). For the 109 filament locations examined in that study, only 32% changed chirality when the sign of β in nearby regions was reversed.

Zoom In Zoom Out Reset image size

Ideally, the magnitude of helicity in each bipole would be based on vector magnetograph data for the real active region in 3D space, but such measurements are unavailable in the corona. Although a number of authors have attempted to estimate helicity in active regions using vector magnetic field measurements in the photosphere (see Démoulin 2007), the horizontal magnetic field components have high uncertainties—including a 180° ambiguity (Metcalf et al. 2006)—which are exacerbated by taking derivatives to compute the vertical current density. Thus we are unable to determine observationally the best value of β to model each active region. We therefore assume that all bipoles in each hemisphere have the same value of β, and run simulations with different values to determine the effect this has on the behavior of the model. Note that taking the same value of β for all regions in each hemisphere does not lead to uniform magnetic helicity in space as the bipoles have different sizes and field strengths. This is well illustrated by Yeates et al. (2008b) in plots of the "current helicity" α = ∇ × B₀ · B₀/B²₀, a measure of the local twist of field lines in a force-free field. One seemingly robust feature of magnetic helicity in the solar corona is the hemispheric pattern, whereby the majority of active regions in the northern hemisphere have negative helicity, and in the southern hemisphere, positive helicity. This is suggested not only by the photospheric vector magnetic measurements (Pevtsov et al. 1995; Abramenko et al. 1996; Hagino & Sakurai 2004), but also by proxy observations such as Hα images of active region structure (Hale 1927; Balasubramaniam et al. 2004), X-ray sigmoids (Canfield et al. 2007), and in situ heliospheric measurements (Smith & Bieber 1993). In accordance with this pattern we choose opposite signs of β in each hemisphere.

To determine the effect of the β parameter on the magnetic flux ropes forming in the global simulation, we consider in this paper a number of simulation runs, as summarized in Table 1. These include runs where the emerging bipole twist takes the observed majority sign in each hemisphere (A2, A4, and A6), runs where it takes the opposite sign (Am2, Am4, and Am6), and a run where the bipoles emerge untwisted (A0). For comparison, we also include a run with no emerging bipoles (AN), a run with lower coronal diffusivity (D4), and a run with a lower outflow speed at the upper boundary (V4).

Table 1. Summary of Parameters in Different Simulation Runs

Run	β in N. hemisphere	β in S. hemisphere	η0 (km2 s−1)	v0 (km s−1)
AN	No emerging regions		45	100
Am6	0.6	−0.6	45	100
Am4	0.4	−0.4	45	100
Am2	0.2	−0.2	45	100
A0	0	0	45	100
A2	−0.2	0.2	45	100
A4	−0.4	0.4	45	100
A6	−0.6	0.6	45	100
D4	−0.4	0.4	22.5	100
V4	−0.4	0.4	45	50

Download table as:  ASCIITypeset image

The basic idea is to search for the distinctive Lorentz force structure of a flux rope. This is illustrated in Figures 3(e)–(h) for the example given in Section 3. The Lorentz force may be decomposed as

$$\begin{equation} \mathbf {j}_0\times \mathbf {B}_0 = \frac{1}{\mu _0}\left(\mathbf {B}_0\cdot \nabla \right)\mathbf {B}_0 - \nabla \left(\frac{B_0^2}{2\mu _0}\right), \end{equation} \tag{ 6 }$$

where the first term on the right-hand side is the magnetic tension T and the second is the magnetic pressure gradient P. Figures 3(e)–(h) show the vertical components of each of these terms T_z (dashed line) and P_z (solid line) as a function of height above the center of the PIL. Note that the two terms are approximately equal and opposite, consistent with a force-free equilibrium with vanishing Lorentz force. The distinctive flux rope structure consists of a sign reversal in both P_z and T_z at the height of the flux rope axis, denoted by the vertical dashed line in Figures 3(f) and (g). In a flux rope, the helical structure is such that the tension force acts inward toward the rope axis, but the magnetic pressure force acts outward (Mackay & van Ballegooijen 2006a, see their Figure 13).

The automated procedure for detecting flux ropes then consists of two stages.

• 1.  
  Point testing. Test individual points on the numerical grid for the characteristic flux rope structure of sign changes in the vertical magnetic pressure and tension forces.
• 2.  
  Clustering. Use a hierarchical clustering algorithm (Johnson 1967) to associate neighboring points which form part of the same flux rope structure.

First consider stage 1. At each point (r_i, θ_i, ϕ_i) in the computational grid the vertical components P_z and T_z of the magnetic pressure and tension forces are computed from the magnetic field B₀. We then require that

$$\begin{eqnarray} &P_z(r_{i-1},\theta _i,\phi _i) < 0, \end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray} &P_z(r_{i+1},\theta _i,\phi _i) > 0, \end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray} &T_z(r_{i-1},\theta _i,\phi _i) > 0, \end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray} &T_z(r_{i+1},\theta _i,\phi _i) < 0, \end{eqnarray} \tag{ 10 }$$

for the point (r_i, θ_i, ϕ_i) to be selected. In practice, we do not need to test all points, but rather test points only up to height r = 1.44 R_☉, and use a coarser "testing grid" of (21, 93, 120) points in the (r, θ, ϕ) directions. In the ϕ direction at the equator this corresponds to one third of the computational grid resolution. The testing grid is chosen such that each point represents an equal 3D volume, by taking equal steps in cos θ, ϕ, and r³. Using a coarser grid reduces computational effort and does not affect the results since we are only interested in well resolved flux rope structures. In order to concentrate the sample on twisted flux ropes, we have implemented the fifth condition of a sufficiently strong parallel current at each point tested, requiring that

$$\begin{equation} \frac{\left|\mathbf {j}_0\cdot \mathbf {B}_0\right|}{B_0^2} > \alpha ^*, \end{equation} \tag{ 11 }$$

with the threshold α* = 0.7 × 10⁻⁸ m⁻¹ determined by experiment. By way of example, Figure 4 (a) shows the flux rope points identified by this first stage on day 202 of run A2, mid-way through the simulation. The points are projected on a plot of B_0r on the photosphere, showing the PIL dividing regions of positive (white) and negative (gray) magnetic polarity. The identified points tend to lie above PILs, as expected from our basic theory of flux rope formation (Section 3). Furthermore, the points are mostly grouped into larger structures. Automated detection of these groupings forms stage 2 of the procedure.

Zoom In Zoom Out Reset image size

The basic clustering idea is simple. Starting with each point as an individual entity, the two closest points are grouped together, and the process repeated until the shortest inter-group distance (in 3D space) is above some threshold. For this threshold we choose 5ΔR_☉, where Δ is the heliographic angle in radians of a computational grid cell at the equator. After running this clustering algorithm, any groups with fewer than 8 points on the testing grid are removed. Again this value was determined by visual inspection of the selected structures. The results after clustering for day 202 of run A2 are shown in Figure 4(b). Each individual group of points, or "flux rope," is identified by color and a number. The actual magnetic field structures corresponding to these flux ropes are illustrated in Figure 4(c), which shows 3D magnetic field lines traced from the selected flux rope points in each structure.

An immediate problem is how to define a flux rope ejection within our simulations. We adopt the practical definition of a large enough radial velocity in the magnetofrictional code. Thus we include both losses of equilibrium following gradual build-up of axial field, and sudden rises caused by nearby bipole emergence. We also include partial lift-offs where only one end of a flux rope opens up. Typically the other end is held down by overlying field from nearby regions. The majority of ejections remove the main, twisted, part of the flux rope through the top boundary of the computational box, usually within several days of the onset of rapid acceleration.

Our automated procedure to find ejections is straight- forward.

• 1.  
  Record the radial velocity v_0r in the magnetofrictional code at each of the selected flux rope points over the simulation, and select those points where v_0r > 0.5 km s⁻¹.
• 2.  
  Cluster these points into separate ejection events.

Although a similar algorithm is used, this clustering procedure differs from that used to define flux ropes on any particular day. Firstly, the clustering is carried out both in space and time, as ejections can take up to a few days to leave the computational domain. Secondly, for the distance measure we consider only the heliocentric angle between points, neglecting any radial distance. This allows for radial movement between daily snapshots. We require a minimum intergroup angle of 4Δ between clusters, and a minimum separation in time of 4 days. After clustering, all groups of points with fewer than 4 points are discounted. These optimal parameter values have been determined by testing the detection procedure over a 13 day test period from day 202 to day 214 in run A2. Locations and days of detected ejections in this test period were compared to those found in a careful manual study of the magnetic field structure. The performance of the automated detection algorithm in this period (with the optimal choice of parameters found) is demonstrated by Table 2, and also illustrated in Figure 5(a), which shows all flux rope points during the period superimposed on the normal photospheric magnetic field for day 208. Panels (b), (c), and (d) correspond to other simulation runs and will be discussed in Section 5.2. The first column of Table 2 lists the days of significant vertical motion of the flux rope in each actual ejection, as detected manually. The approximate locations in latitude and longitude are listed in the second column, and plotted as black triangles in Figure 5(a). The third and fourth columns of Table 2 give the results of the automated detection algorithm (with the optimal parameters). The clusters of points selected in this procedure are colored red in Figure 5(a). Note that only points falling within the 13 day period are shown in the figure, even if the cluster overlaps with the period. In Table 2, the number of clustered points refers to the total number in that cluster, whether in the 13 day period or not. The code detects 16 of the 19 ejections (84%), and the number of clustered points (fourth column) gives an estimate of the size of each event. Of the three missed ejections, two are small, and the third (located at −60, 5) is quite large but rises slowly and gradually, with no sudden loss of equilibrium. We note that there are also four "false positives," where spurious ejections are detected. However, two of these are small "flux rope" structures associated with temporary coronal reconfiguration following new bipole emergence, and the other two were actually considered on manual inspection to be part of other ejections. We conclude that this procedure performs well and the number of ejections detected is accurate to within ±15%.

Zoom In Zoom Out Reset image size

Figure 7 shows how the properties of flux ropes evolve over time in each simulation run. These properties include the total number of flux rope points (Figure 7(a)), the mean latitude of flux rope points (Figure 7(b)), the number of flux ropes found by the clustering stage (Figure 7(c)), and the mean size of each of these clusters (Figure 7(d)).

Zoom In Zoom Out Reset image size

A general feature in all of these properties is an initial increase, followed in most cases by a saturation after 50 or more days. This illustrates the gradual build-up of helicity in the corona as the simulation evolves away from the initial potential field. Because run AN (with no bipole emergence) also shows an initial increase in the number of flux rope points (Figure 7(a)), we see that surface shearing alone creates flux ropes over this timescale. However, the faster increase for the other runs shows that emerging helicity speeds up the formation of flux ropes. Indeed, Figure 7(b) shows that more flux ropes form at lower latitudes in the runs with bipole emergence, supporting this point. The emergence of new flux also acts to sustain the coronal magnetic field against the diffusive decay which would otherwise occur over hundreds of days. This decay is evident in run AN as reformation of ejected flux ropes and production of new flux ropes declines after day 210.

Note that all runs show a temporary decrease in the number of flux rope points and in their mean latitude around day 210. This is caused by the ejection of a pair of high-latitude flux ropes, illustrated in a time sequence in Figure 8 (for run A2). The two ejections are labeled E1 and E2, and take place from days 199 and 204, respectively. Ejection E2 is also seen in Figure 5(a), but E1 occurred before the period shown there. A similar high-latitude ejection causes a dip in the number of flux rope points around day 260 for run AN.

Zoom In Zoom Out Reset image size

5.1.1. Effect of Emerging Bipole Twist

5.1.2. Effect of Radial Outflow

5.1.3. Effect of Coronal Diffusion

Using the automated procedure in Section 4.2, we find the number of ejections per day to have a ramp-up phase of about 40 days, before fluctuating considerably over the rest of the simulation, as well as varying between runs. Figure 9 shows the total number of ejections in each run, with error bars showing the estimated errors in the automated detection procedure (Section 4.2). We immediately see the importance of emerging bipoles for flux rope ejections. In run AN there are only 23 ± 4 ejections, as compared to 108 ± 16 in run A0, and more for the runs with emerging bipole twist—up to 202 ± 30 for run A6. Partly, this reflects the increased number of flux ropes in runs with emerging bipoles, and partly it reflects the influence of strong emerging magnetic fields in reconfiguring the coronal field. The increasing number of ejections with greater emerging bipole twist reflects the quicker formation of flux ropes with strong axial fields, which then lose equilibrium. For example, compare panels (a) and (c) in Figure 5, which show flux rope points and radial velocities for runs A2 and A4. The flux rope labeled M in Figure 5(c), internal to a newly emerged bipole, both forms and is ejected in run A4. However, in run A2 there is not sufficient axial magnetic field at this location for a flux rope to be detected at this time.

Zoom In Zoom Out Reset image size

There is also a tendency for the runs with the majority hemispheric sign of bipole twist (A2, A4, and A6) to have more ejections than the runs with the same magnitude but opposite sign of twist. This may be related to the larger size of mid-latitude flux ropes in the first set of runs, as described in Section 5.1. This is illustrated by comparing run Am2 with run A2 in Figure 5 (panels (b) and (a), respectively). The black boxes in Figure 5(b) show decaying bipolar regions where no flux ropes are detected in run Am2. In run A2, both of these locations have already formed strong flux ropes which are ejected during this particular period.

From run V4—shown by a square in Figure 9—we see that halving the outflow speed has no significant influence on the number of ejections (giving 178 ± 27). So, just as the outflow has no major influence on the formation of flux ropes, it has no major influence on their loss of equilibrium. As before, this is because it acts only near the upper boundary, which a flux rope reaches only after equilibrium is lost. In contrast, halving the coronal diffusivity η₀ increases the number of ejections to 233 ± 35 in run D4 (triangle in Figure 9). An example of an additional ejection not occurring in run A4 is labeled N in Figure 5(d). The higher number of ejections is explained by the larger, more twisted flux ropes which are able to form in this case.

Finally, how does the number of ejections produced in our simulations compare with the number of CMEs observed on the real Sun during this period? The CDAW catalog (Yashiro et al. 2004) is the standard manually compiled list of CMEs observed by the SOHO/LASCO coronagraphs (Brueckner et al. 1995). To avoid the initial ramp-up phase in the simulation, we compare the rates of ejections between day 183 (1999 July 2) and day 283 (1999 October 10). Table 3 shows the number of flux rope ejections per day for each simulation run over this period, in addition to the observed rate from the CDAW catalog. Two observed rates are given: the first includes all events in the catalog, and the second omits events labeled "poor" by the LASCO operator. We believe the second rate to be more appropriate for this comparison, as our global simulations relate to large-scale flux rope events, which are unlikely to be labeled "poor." This rate will be a lower estimate for two reasons: firstly not all far-side events are observed, and, secondly, there are several gaps in the instrument coverage over the period. Using this rough estimate, our simulations produce flux rope ejections at about 50% of the observed CME rate (after the initial ramp-up phase).