© ESO, 2014

1. Introduction

Coronal mass ejections (CMEs) are a common occurrence in the Sun’s atmosphere that are known to release giga-tons of plasma into interplanetary space. Some of the ejected plasma can reach the space environment of the Earth and have a strong and complex influence on space activity by inducing geospace disruptions that can severely impact spacecraft, power grids, and communication (Baker et al. 2013). While CMEs are quite commonly observed (Evans et al. 2013), especially during the peak of the solar cycle, they are still poorly understood. Some of the biggest CME mysteries pertain to their origin, propagation, and relation to flares.

The initiation of CMEs has been widely studied and yet remains largely unexplained (see reviews by Forbes et al. 2006; Chen 2011). However, many observational studies of associated features have led to clues about how they occur and what factors contribute to their destabilization (see review by Gopalswamy et al. 2006). Prior to an eruption, large-scale shear motions are often observed in photospheric images, especially about the magnetic neutral line (Krall et al. 1982) and in the form of sunspot rotations (Tian & Alexander 2006). In addition, patches of magnetic flux are found to emerge, expand, move, fragment, coalesce, and cancel over a wide range of length and time scales (Sheeley 1969; Zwaan 1985; Centeno et al. 2007; Parnell et al. 2009). It is believed that shear motions, sunspot rotation, and the emergence of new flux are all related to the injection of magnetic helicity into coronal magnetic structures that could be directly involved in the eruption (Chae 2001; Kusano et al. 2002; Démoulin et al. 2002; Pariat et al. 2006; Magara & Tsuneta 2008).

In addition to the growing body of observational studies that have improved our understanding of CMEs, many new insights have also emerged from theoretical and numerical efforts. CMEs have been modeled in two and three dimensions using both simple analytical methods and sophisticated magnetohydrodynamic (MHD) simulations (see Jacobs & Poedts 2011, and references therein). These models differ widely in physical and numerical details, each making its own choice of how to address the trade-off between complexity and computational feasibility.

Early theoretical models explained CMEs as a loss of equilibrium, due to magnetic buoyant instabilities (e.g., van Tend & Kuperus 1978; Low 1981; Demoulin & Priest 1988), as well as MHD flows (Low 1984) and reconnection (Forbes & Isenberg 1991). Forbes & Priest (1995) proposed a CME model based on the movement of magnetic footpoints (sources) below a flux rope and the subsequent development of a singular current sheet, through which a large magnetic energy release should take place as the flux rope moves continually outwards. Lin & Forbes (2000) refined their model and computed exact solutions for the energy release, flux rope height, current sheet length, and reconnection rate. The Lin & Forbes (hereafter, “LF”) model, while simplistic, provides an important step forward in CME modeling because it offers exact solutions to the time-dependent nonlinear problem of a flux rope eruption and includes more than a heuristic treatment of magnetic reconnection. Furthermore, it predicts many features (e.g., morphology, current sheet, post-flare loops, flows, energetics) confirmed by observations (Ciaravella et al. 2002; Ko et al. 2003; Lin et al. 2005).

A similar two-dimensional flux rope model was proposed by Chen & Shibata (2000). Like LF, the Chen & Shibata (“CS”) model consists of a two-dimensional configuration in which a flux rope sits above the photosphere, surrounded by a line-tied coronal arcade. In both models, the magnetic equilibrium is destabilized by photospheric driving, causing a current sheet to form in the flux rope’s wake as it moves outwards. However, whereas the LF model calls for a somewhat manufactured mechanism for destabilization via large-scale convergence of the sources, the CS model improves upon the LF model by incorporating flux emergence as the driver.

While it does not lend itself to a purely analytical treatment, the CS model is suitable for numerical simulation. The authors report four very different outcomes based on the position and direction of the driving, showing that the location of the emergence per se is not a critical factor for destabilizing the coronal flux rope but rather that the relative orientation of the emerging flux determines whether the flux rope moves outwards/upwards (CME-like) or inwards/downwards (failed eruption).

Several subsequent studies have built upon the CS model. For example, Chen et al. (2004), Shiota et al. (2003), and Shiota et al. (2004) produced synthetic emission images from CS simulations to compare morphological features, reconnection in-flows, and coronal dimmings found in actual CME observations. Moreover, Shiota et al. (2003, 2005) were able to identify the formation, structure, and location of slow and fast shocks in the CMEs produced in these simulations. Gravitational density stratification in an isothermal atmosphere was considered by Chen et al. (2004); Shiota et al. (2004) and also in a later study by Dubey et al. (2006) in spherical coordinates with axisymmetry.

In this study, we re-examine the CS model using a more sophisticated numerical tool, a more realistic atmosphere, and higher spatial resolution than previous studies. Simulations are performed using a high-order spectral element method with numerically accurate, self-consistent treatments of diffusive transport (i.e., resistivity, viscosity, and thermal conduction). In addition, we reformulate the initial conditions to have magnetic fields that are everywhere continuous and differentiable, and to include a solar-like temperature profile with a sharp transition region and density stratification.

Through an exploration of physical parameters, we find that the CS magnetic equilibrium can be unstable even without a flux emergence driver. Linear theory has shown that sufficient perturbation of the field lines near an X-point by waves or motion can disrupt the balance between magnetic pressure and magnetic tension, causing the X-point to collapse and form a reconnecting current sheet (Priest & Forbes 2000, Chap. 2). Our simulations demonstrate that under a wide range of conditions the CS equilibrium is susceptible to such a collapse via nonlinear accumulation of fast magnetosonic waves at the X-point (McLaughlin et al. 2009). However, we also show that in a sufficiently incompressible plasma due, for example, to the presence of a background “guide” magnetic field co-aligned with the axis of the flux rope, the X-point collapse does not take place and the CS magnetic equilibrium can be stabilized. For both stable and unstable configurations, we investigate the impact of the resistivity model enabling magnetic reconnection below the flux rope, as well as the plasma parameters in the low solar atmosphere, on the flux rope’s response to the flux emergence driver. We show that flux emergence can produce a rising flux rope both in a stratified and an unstratified atmosphere, though the resulting ejection speed, as well as the plasma dynamics around the X-point, can be strongly effected by the magnitude of the guide field and the atmospheric stratification.

2. Model

The CS model has a two-dimensional domain with motion and magnetic field allowed perpendicular (as well as parallel) to the plane of the domain (v,B ∈ R³). Therefore, we can write the magnetic field, normalized to some value B₀, in terms of a scalar potential ψ representing the in-plane flux, and an out-of-plane scalar field: ${B}\mathrm{\left(}\mathit{x,y}\mathrm{;}\mathit{t}\mathrm{\right)}\mathrm{=}\mathrm{\nabla }\mathrm{(}\mathrm{-}\mathit{\psi }\mathrm{)}\mathrm{\times }\hat{{e}}\mathit{z}\mathrm{+}{\mathit{b}}_{\mathit{z}}\hat{{e}}z$(1)All quantities are normalized in terms of the first three constants found in Table 1: L₀ = 5 Mm, which is the unit of length; B₀ = 10 G, the unit of magnetic field strength; and N₀ = 10⁹ cm^-3, the unit of number density. Given the Alfén velocity ${\mathit{v}}_{\mathrm{A}}\mathrm{\equiv }{\mathit{B}}_{\mathrm{0}}\mathit{/}\sqrt{{\mathit{\mu }}_{\mathrm{0}}{\mathit{m}}_{\mathrm{p}}{\mathit{N}}_{\mathrm{0}}}$, where m_p is the proton mass, we define the unit time as τ ≡ L₀/v_A, unit temperature as ${\mathit{T}}_{\mathrm{0}}\mathrm{\equiv }{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}\mathit{/}\mathrm{\left(}{\mathit{\mu }}_{\mathrm{0}}{\mathit{k}}_{\mathrm{B}}{\mathit{N}}_{\mathrm{0}}\mathrm{\right)}$, and unit pressure as ${\mathit{P}}_{\mathrm{0}}\mathrm{\equiv }{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}\mathit{/}{\mathit{\mu }}_{\mathrm{0}}$. The solar surface gravity g_S = 274 m/s² is similarly normalized as g ≡ g_S(τ/v_A) = 2.88 × 10^-3.

Due to symmetries intrinsic to the model, only half of the domain in the horizontal direction has to be resolved (x> 0). Thus, simulations are performed in a computational domain (x,y) ∈ [0,L_x] × [0,L_y], with the solar convection zone assumed to be located below the domain (y< 0).

2.1. Initial conditions

2.1.1. Magnetic configuration

The initial magnetic configuration prescribed in the CS model consists of a coronal flux rope of radius r₀ surrounded by an arcade of “loops” that are line-tied in the photosphere.

The flux rope contains a current channel that is mirrored by an image current far below the photosphere (outside the computational domain), and four line currents produce a potential quadrupolar field just below the photosphere. Since the bottom boundary of the numerical domain coincides with the photosphere, the only visible current initially is that within the flux rope.

In the original CS study, the coronal flux rope is given by a flux function ψ_l that results in a discontinuous current density at the edge of the flux rope. Therefore, we propose the following alternative: $\begin{array}{ccc}{\mathit{\psi }}_{\mathit{l}}\mathrm{=}\{\begin{array}{c}\\ & \\ & \end{array}& & \end{array}$where r² = x² + y², and the center of the flux rope lies at (x = 0, y = h).

Our formulation for ψ_l lends itself to a continuous current density: ${\mathit{j}}_{\mathit{l}}\mathrm{=}\mathrm{-}{\mathrm{\nabla }}^{\mathrm{2}}{\mathit{\psi }}_{\mathit{l}}\mathrm{=}\{\begin{array}{c}\\ \frac{\mathrm{4}{\mathit{r}}^{\mathrm{2}}}{{\mathit{r}}_{\mathrm{0}}^{\mathrm{3}}}\mathrm{-}\frac{\mathrm{4}}{{\mathit{r}}_{\mathrm{0}}}\mathit{,}& \mathit{r}\mathrm{\le }{\mathit{r}}_{\mathrm{0}}\\ \mathrm{0}\mathit{,}& \mathit{r}\mathit{>}{\mathit{r}}_{\mathrm{0}}\mathit{.}\end{array}$(3)The other flux components of the initial configuration, representing the image current and line currents, respectively, are kept as originally defined $\begin{array}{ccc}& & \\ & & \end{array}$with r₀ = 0.5, h = 2 and c = 0.25628. All three flux functions are summed to produce the initial magnetic equilibrium, shown in Fig. 1: $\mathit{\psi }\mathrm{=}{\mathit{\psi }}_{\mathit{l}}\mathrm{+}{\mathit{\psi }}_{\mathit{i}}\mathrm{+}{\mathit{\psi }}_{\mathit{b}}\mathit{.}$(6)In addition to the line currents, which produce the in-plane magnetic field, we also allow for a uniform background magnetic field out of the plane ${\mathit{b}}_{\mathit{z}\mathrm{0}}\hspace{0.17em}\hat{{e}}{}_{\mathit{z}}$. This “guide” field contributes magnetic pressure but no current.

Fig. 1 Contours of ψ in the initial conditions.

Fig. 2 Contours of ψ (black) and bz (magenta) at t = 0 for bz as a function of r (left) and as a function of ψ (right).

In the original CS study, a density spike is applied to support the flux rope against radial compression: an outward-acting pressure gradient force offsets the inward-acting Lorentz force due to the flux rope’s poloidal field. Equivalently, the flux rope can be supported against the radial Lorentz force by magnetic pressure, as in Shiota et al. (2005): in addition to the background guide field b_z0, we apply an additional axial field in the flux rope which is highest in the center and diminishes over a radius of r₀.

We note that if the axial field b_z is specified as a function of the flux rope radius alone the flux rope will not be force-free, as the contours of ψ are not perfectly circular due to the small but finite contributions of ψ_b and, to a lesser extent, of ψ_i to the total flux in the coronal flux rope. It can be seen from the left panel of Fig. 2 that, given such a function b_z(r), the contours of ψ and b_z would not be well-aligned.

To avoid the misalignment and minimize unbalanced Lorentz forces in the initial condition, we instead choose to specify b_z as a function of ψ, as follows: $\begin{array}{ccc}& & \\ & & \end{array}$(The derivation of the above equations can be found in Appendix A.) Note that ζ = 0 when ψ = ψ₀ − r₀/ 4, and ζ = r₀ when ψ = ψ₀ + r₀/ 2, where ψ₀ ≡ [ψ_i + ψ_b] | _{(x,y) = (0,h)}.

2.1.2. Unstratified atmosphere

In the case of an unstratified atmosphere, the number density field is initialized to a uniform value of n = n₀ = 1(N₀). The pressure field of an electron-proton plasma can be determined by the following equation of state: $\mathit{p}\mathrm{=}\mathrm{2}\mathit{nT}\mathit{.}$(9)We choose a uniform initial temperature, so the initial pressure p = p₀ is also uniform. The free parameter p₀ is chosen variably in the simulations to yield temperatures close to coronal values, as well as low plasma β ≡ 2p/B².

An unstratified atmosphere has the advantage of isolating the flux rope dynamics from the thermodynamics. By controlling p₀, one essentially explores different regimes of the plasma β.

2.1.3. Stratified atmosphere

We also attempt to simulate a solar-like atmosphere by modeling the average vertical temperature profile as a hyperbolic tangent function, as in Leake & Arber (2006); Leake & Linton (2013): $\mathit{T}\mathrm{\left(}\mathit{y}\mathrm{\right)}\mathrm{=}\frac{{\mathit{T}}_{\mathrm{p}}}{{\mathit{T}}_{\mathrm{0}}}\mathrm{+}\left(\frac{{\mathit{T}}_{\mathrm{c}}\mathrm{-}{\mathit{T}}_{\mathrm{p}}}{\mathrm{2}{\mathit{T}}_{\mathrm{0}}}\right)\left[\mathrm{1}\mathrm{+}\tanh \left(\frac{\mathit{y}\mathrm{-}{\mathit{y}}_{\mathrm{TR}}\mathit{/}{\mathit{L}}_{\mathrm{0}}}{\mathrm{\Delta }\mathit{y}\mathit{/}{\mathit{L}}_{\mathrm{0}}}\right)\right]$(10)with photospheric temperature T_p = 5000 K, coronal temperature T_c = 10⁶ K, transition region height y_TR = 2.5 Mm, and transition region width Δy = 0.5 Mm.

Given this temperature profile, we seek compatible density and pressure profiles such that the plasma is in hydrostatic equilibrium: $\frac{\mathrm{d}\mathit{p}}{\mathrm{d}\mathit{y}}\mathrm{+}\mathit{ng}\mathrm{=}\mathrm{0}$(11)with the constant gravitational acceleration g pointed in the $\mathrm{-}\hat{{e}}{}_{\mathit{y}}$ direction.

We solve Eq. (11)using Eqs. (9)and (10)(see the derivation in Appendix B). The resulting pressure profile is $\begin{array}{ccc}\mathit{p}\mathrm{\left(}\mathit{y}\mathrm{\right)}& \mathrm{=}& {\mathit{p}}_{\mathrm{0}}\exp \{\frac{\mathit{g}\mathrm{\Delta }\mathit{y}\mathit{/}{\mathit{L}}_{\mathrm{0}}}{\mathrm{2}{\mathit{T}}_{\mathrm{c}}\mathit{/}{\mathit{T}}_{\mathrm{0}}}[\mathrm{-}\frac{\mathit{y}\mathrm{-}{\mathit{y}}_{\mathrm{TR}}\mathit{/}{\mathit{L}}_{\mathrm{0}}}{\mathrm{\Delta }\mathit{y}\mathit{/}{\mathit{L}}_{\mathrm{0}}}\\ & & \mathrm{+}\frac{{\mathit{T}}_{\mathrm{c}}\mathrm{-}{\mathit{T}}_{\mathrm{p}}}{\mathrm{2}{\mathit{T}}_{\mathrm{p}}}\ln \left(\frac{{\mathit{T}}_{\mathrm{p}}}{{\mathit{T}}_{\mathrm{0}}}\exp \left[\mathrm{-}\frac{\mathrm{2}\mathrm{(}\mathit{y}\mathrm{-}{\mathit{y}}_{\mathrm{TR}}\mathit{/}{\mathit{L}}_{\mathrm{0}}\mathrm{)}}{\mathrm{\Delta }\mathit{y}\mathit{/}{\mathit{L}}_{\mathrm{0}}}\right]\mathrm{+}\frac{{\mathit{T}}_{\mathrm{c}}}{{\mathit{T}}_{\mathrm{0}}}\right)]\}\mathrm{·}\end{array}$(12)Here the parameter p₀ corresponds to a constant of integration that shifts the entire pressure profile of the atmosphere. As in the unstratified case, this affects the values of β, which should be high (~10) in the photosphere and low (~10^-2) in the corona. Therefore, we do not vary p₀ for the stratified simulations.

Fig. 3 Logarithm (base 10) of normalized number density n (green), normalized pressure p (blue), and temperature T (red) as a function of height in the initial conditions for a stratified atmosphere.

Profiles of the initial density, pressure, and temperature in a stratified atmosphere are plotted in Fig. 3. It is evident that all three quantities vary smoothly and by multiple orders of magnitude. Furthermore, this transition occurs well below the height of the flux rope core (y = 2).

2.2. Numerical method

We implement this initial configuration in the high-fidelity numerical simulation framework, HiFi, which makes use of high-order spectral elements and implicit time-stepping (Lukin 2008; Lukin & Linton 2011). As a strong condition, HiFi requires all variables to be represented by continuous functions in the initial and boundary conditions. Therefore, the magnetic field must be everywhere differentiable, implying . In addition, the boundary driving of flux needs to be differentiable in both space and time, in order for the electric and magnetic fields to be smooth. These conditions are well satisfied by the initial conditions described above.

In this work, HiFi is used to integrate in time the following equations of visco-resistive MHD: $\begin{array}{ccc}& & \frac{\mathit{\partial n}}{\mathit{\partial t}}\mathrm{+}\mathrm{\nabla }\mathrm{·}\mathrm{(}\mathit{n}{{v}}^{\mathrm{)}}\mathrm{=}\mathrm{0}\\ & & \frac{\mathit{\partial }\mathrm{(}\mathrm{-}\mathit{\psi }\mathrm{)}}{\mathit{\partial t}}\mathrm{=}\mathrm{-}{v}\mathrm{\times }{B}\mathrm{+}\mathit{\eta }{\mathit{j}}_{\mathit{z}}\\ & & \frac{\mathit{\partial }{\mathit{b}}_{\mathit{z}}}{\mathit{\partial t}}\mathrm{+}\mathrm{\nabla }\mathrm{·}\left({\mathit{b}}_{\mathit{z}}{v}\mathrm{-}{\mathit{v}}_{\mathit{z}}{B}\right)\mathrm{=}\mathrm{\nabla }\mathrm{·}\left(\mathit{\eta }\mathrm{\nabla }{\mathit{b}}_{\mathit{z}}\right)\\ & & \frac{\mathit{\partial n}{v}}{\mathit{\partial t}}\mathrm{+}\mathrm{\nabla }\mathrm{·}\left\{\mathit{n}{v}{v}\mathrm{+}\mathit{p}{I}\mathrm{-}\mathit{\mu n}\mathrm{[}\mathrm{\nabla }{v}\mathrm{+}{\left(\mathrm{\nabla }{v}\right)}^{\mathit{T}}\mathrm{]}\right\}\mathrm{=}{j}\mathrm{\times }{B}\\ & & \frac{\mathrm{3}}{\mathrm{2}}\frac{\mathit{\partial p}}{\mathit{\partial t}}\mathrm{+}\mathrm{\nabla }\mathrm{·}\left(\frac{\mathrm{5}}{\mathrm{2}}\mathit{p}{v}\mathrm{-}\mathit{\kappa }\mathrm{\nabla }\mathit{T}\right)\mathrm{=}{v}\mathrm{·}\mathrm{\nabla }\mathit{p}\mathrm{+}\mathit{\eta }{\mathit{j}}^{\mathrm{2}}\mathrm{+}\mathit{\mu n}\mathrm{[}\mathrm{\nabla }{v}\mathrm{+}{\left(\mathrm{\nabla }{v}\right)}^{\mathit{T}}\mathrm{]}\mathrm{:\nabla }{v}\end{array}$with an auxiliary equation (Ampère’s law): $\mathrm{\nabla }\mathrm{\times }{B}\mathrm{=}{j}$(13f)The normalized transport coefficients found in Eqs. (13)–namely κ, μ, η – control the level of dissipation of the MHD fluid quantities through molecular diffusion: temperature, velocity, and current, respectively. Each of these three transport parameters is chosen to be compatible with the resolution and objective of each simulation. Further, for some of the simulations (see below), we allow the resistivity η to be a function of local current density, η = η_bg + η_anom(j), where $\begin{array}{ccc}{\mathit{\eta }}_{\mathrm{anom}}\mathrm{\left(}{j}\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ & \\ & \\ & \end{array}& & \end{array}$η_bg is the uniform and time-independent background resistivity, and η_anom is some “anomalously enhanced” effective resistivity, , occurring due to microphysics not captured by the MHD model whenever the current density rises above the critical current density j_c.

2.2.1. Boundary conditions

The bottom boundary of the simulation domain, representing the photosphere, is perhaps the most important boundary condition affecting the outcome of a simulation. Flux emergence is achieved by varying the flux function at this boundary in time, which is equivalent to applying an electric field. This electric field determines the evolution of the magnetic field, which can be advected in or out of the domain or resistively dissipated, as described by Ohm’s Law: $\frac{\mathit{\partial \psi }}{\mathit{\partial t}}\mathrm{=}{\mathit{E}}_{\mathit{z}}\mathrm{=}\mathrm{-}\hat{{e}}\mathit{z}\mathrm{·}{v}\mathrm{\times }{B}\mathrm{+}\mathit{\eta }{\mathit{j}}_{\mathit{z}}$(14)The resistive component, the second term on the right-hand side of Eq. (14), is determined by the geometry of the magnetic field at any given time. Therefore, by varying the flux at the boundary (∂ψ/∂t) in a prescribed way, we also induce cross-field plasma motions (v × B).

Chen & Shibata (2000) prescribe two cases of localized boundary driving, namely, over a region | x − x₀ | ≤ 0.3 centered at x₀ = 0 (case A) and at x₀ = 3.9 (case B). We apply the same method only for the case of x₀ = 0, but use a formulation that is smoother in time and in space: $\begin{array}{ccc}& & \mathit{\psi }\mathrm{\left(}\mathit{x,}\mathrm{0}\mathrm{;}\mathit{t}\mathrm{\right)}\mathrm{=}\mathit{\psi }\mathrm{\left(}\mathit{x,}\mathrm{0}\mathrm{;}\mathrm{0}\mathrm{\right)}\mathrm{+}\frac{{\mathit{\psi }}_{\mathrm{e}}\mathrm{\left(}\mathit{x}\mathrm{\right)}}{\mathrm{2}}\left[\frac{\mathit{t}}{{\mathit{t}}_{\mathrm{e}}}\mathrm{-}\frac{\sin \mathrm{(}\mathrm{2}\mathit{\pi t}\mathit{/}{\mathit{t}}_{\mathrm{e}}\mathrm{)}}{\mathrm{2}\mathit{\pi }}\right]\mathit{,t}\mathrm{\le }{\mathit{t}}_{\mathrm{e}}\\ & & \end{array}$(15)where t_e is the duration over which the electric field drive is applied at the boundary. For t>t_e, the photospheric boundary is treated as a perfect conductor.

We do not allow any in-plane flow on the bottom photospheric boundary (v_x = 0, v_y = 0) and force the normal gradients of b_z,j_z, and temperature to be zero: i.e., . The unstratified atmosphere also has , while the stratified case imposes a fixed value of the density scale height .

The left boundary is a symmetry boundary, with odd symmetry required for the horizontal and out-of-plane components of flow, v_x and v_z, and even symmetry imposed on all other dependent variables. At the outer boundaries (top and right), the gradients of density, b_z, j_z, and pressure are zero, and flow is only allowed in the vertical direction. Table 2 provides a simple reference for the various boundary conditions applied in the simulations.

2.2.2. Dissipative boundary layers

Chromosphere.

The flux emergence represented by Eq. (15)changes the flux function just at the boundary but has no direct effect on ψ anywhere else, including just above it. While Ohm’s Law Eq. (14)does relate flux evolution to fluid transport, it does not guarantee that the flux function and other quantities will be well-behaved for all time, particularly in the y-direction. Therefore, to allow flux to slip more easily through the region just above the photospheric driving (effectively, the “chromosphere”), we apply a resistive boundary layer by enhancing η locally according to the following function: $\mathit{\eta }\mathrm{=}{\mathit{\eta }}_{\mathrm{bg}}\mathrm{+}{\mathit{\eta }}_{\mathrm{anom}}\mathrm{+}{\mathit{\eta }}_{\mathrm{ph}}\exp \left(\mathrm{-}\frac{{\mathit{y}}^{\mathrm{2}}}{{\mathit{y}}_{\mathrm{0}}^{\mathrm{2}}}\right)$(16)where η_ph is the photospheric value of resistivity, and y₀ = 0.2 corresponds to the height of 1 Mm above the photosphere. Conceptually this boundary layer emulates the enhanced resistivity of the chromosphere due to collisional impedance by neutrals (Leake & Arber 2006).

Outer corona.

Separately, to mimic an open domain boundary with no wave reflections in the corona, a viscous boundary layer is prescribed close to the outer coronal (top and right) boundaries. Starting at a distance of d = 0.5 inward from the boundary of the computational domain, μ is increased gradually towards the domain boundary (according to a cosine profile) from a background value of μ_bg up to the outer boundary value of μ_out = 1. In the same fashion, near the outer boundaries with d = 0.5, resistivity is ramped up from η_bg to a boundary value of η_out = 10^-2.

3. Results

In this section, we discuss simulations of the undriven equilibrium, as well as driven simulations of flux emergence in both the stratified and unstratified cases. While the CS initial conditions describe an approximate equilibrium, we find that this equilibrium can be unstable to small perturbations. We discuss the role of MHD waves in destabilizing the flux rope via X-point collapse and the role played by plasma compressibility in stabilizing the X-point and the flux rope. Finally, we discuss the driven simulations of stable and unstable equilibria with different resistivity models, as well as details of the resulting eruption process.

Fig. 4 No driving. Four snapshots of current density jz showing outward propagating MHD waves. Notice the trapping and interference of waves at the X-point; compounding of these waves here precipitates an X-point collapse, leading to the formation of a current sheet and initiation of magnetic reconnection.

3.1. Flux rope stability

In Chen & Shibata (2000), the coefficient c (see Eq. (5)above) was determined by trial and error to yield a magnetic equilibrium such that the center of the flux rope did not move “for long enough time”. (It can also be derived by requiring that the vertical component of the Lorentz force is zero at the center of the flux rope, ŷ· [j × B] | _{(x,y) = (0,h)} = 0.) Our numerical experiments confirm that this value for c is indeed appropriate for equilibrium, but we find that the equilibrium itself is tentative and unstable to perturbations.

3.1.1. X-point collapse

At the beginning of each simulation, the flux rope – which is approximately force-free – goes through a small adjustment to settle into an actual force-free equilibrium. In Sect. 2, we explained how formulating b_z as a function of ψ forces the contours of b_z and ψ to be aligned. Although this is an improvement on the original setup, the configuration is still not perfectly force-free because the contours of j_z = −∇²ψ = −∇²ψ_l are still circular and therefore not aligned with the contours of ψ, producing a small but finite Lorentz force. The adjustment to correct the misalignment, however small it may be, is sufficient to generate waves that may destabilize the flux rope.

Figure 4 illustrates the oscillation of current density induced by the fast magnetosonic waves emanating from the flux rope as it adjusts to the initial conditions. The distribution of the waves is not uniform because the initial adjustment is one where the flux rope is squeezed in one direction (horizontally) while expanding in the other direction (vertically). Therefore, the waves propagating horizontally are out of phase with those propagating vertically. It is interesting to note that the propagation of both types of waves is initially radial but eventually becomes oblique at the flanks of the active region due to the inhomogeneity of the magnetic pressure in the corona.

The fast waves themselves do not directly destabilize the flux rope. However, the entire equilibrium can be destabilized when the waves reach the X-point below the flux rope and cause it to collapse. The role of fast waves in X-point collapse for a zero-β plasma has been studied by McLaughlin et al. (2009) and their behavior in the neighborhood of X-points has been investigated by similar earlier studies (McLaughlin & Hood 2004, 2005, 2006). As described in these studies, we find that the fast waves approach the X-point but tend to get trapped there if their phase speed becomes too low at the X-point. The trapping occurs because the waves are refracted towards the null and then wrap around it if they cannot pass through it. As a result, the waves push current towards the null where it accumulates exponentially in the linear regime. The buildup of waves at the null, however, quickly leads to nonlinear behavior, forming shocks and jets, which deform the X-point into a cusp-like geometry, which flattens and forms a current sheet (McLaughlin et al. 2009). This fast-wave accumulation and resulting collapse of the X-point are evident in the sequence of figures in Fig. 4.

The consequence of X-point collapse occurring as a result of fast wave accumulation at the null is that it forms a current sheet separating anti-parallel fields. The formation of the current sheet, in turn, kicks off magnetic reconnection that drives itself for as long as there is free magnetic energy available in the system. When the collapse forms a horizontal current sheet, the reconnection process draws in the flux rope from above, and it pulls itself down towards the photosphere to draw in flux from below, destroying the original configuration. Formation of a vertical current sheet similarly leads to a CME eruption.

Other factors that may contribute to a collapse of the X-point in the absence of driving include boundary flows, likely related to the reflection of waves, and the asymmetric resistivity model, which intentionally biases η in the y-direction in order to allow magnetic flux to slip through the photosphere (see previous section). However, we found these effects to be sub-dominant to the fast wave accumulation at the X-point in destabilizing the flux rope.

Fig. 5 Dependence of flux rope stability on the free parameters p0 and bz0 (left), as well as on the plasma β and compressibility measure Γ (right). All simulations are performed without driving. The different symbols signify that the flux rope was stable (black dots); was wobbly but on average did not rise or sink (black stars); moved upwards (blue triangles); or moved downwards (red triangles).

3.1.2. Sensitivity to compressibility

We find that sufficient magnetic pressure and/or gas pressure at the X-point can suppress the X-point collapse. In the absence of a guide field (b_z0 = 0), the magnetic pressure drops to zero at the X-point. Then, with low gas pressure, the fast wave speed is reduced drastically at the X-point causing wave refraction and accumulation. However, a parametric exploration of p₀ and b_z0 in an unstratified atmosphere revealed that increasing either of these parameters helps to stabilize the flux rope. The left panel of Fig. 5 is a graphical chart of the many simulations that were performed scanning the parameter space of p₀ and b_z0. Black dots represent simulations in which the flux rope was stable over many hundreds of Alfvén times (no X-point collapse); blue triangles represent those in which the flux rope experienced a slow rise (vertical collapse); red triangles represent those in which the flux rope descended (horizontal collapse); and black stars represent those in which the flux rope moved up and down but on average maintained the same height in the atmosphere (oscillatory X-point collapse).

One could argue heuristically that increasing either p or b_z effectively decreases the compressibility of the plasma (or increases the stiffness of the medium), so any motions at the X-point need to do more work against the gas pressure or magnetic pressure to force a collapse of the magnetic topology. Therefore, we propose a generalized measure of two-dimensional plasma compressibility: $\mathrm{\Gamma }\mathrm{=}\frac{{\mathit{b}}_{\mathrm{\perp }}^{\mathrm{2}}}{\mathrm{2}\mathit{p}\mathrm{+}{\mathit{b}}_{\mathit{z}}^{\mathrm{2}}}$(17)and relate the free parameters of the simulations, p₀ and b_z0, to the magnitude of the in-plane field b_⊥ ≈ 1 (B₀), as well as to the initial background plasma β: $\mathit{\beta }\mathrm{=}\frac{\mathrm{2}{\mathit{p}}_{\mathrm{0}}}{{\mathit{b}}_{\mathrm{\perp }}^{\mathrm{2}}\mathrm{+}{\mathit{b}}_{\mathit{z}\mathrm{0}}^{\mathrm{2}}}\mathrm{·}$(18)The right panel of Fig. 5 provides an alternative way to assess the effect of the initial parameters on the stability of the flux rope and the X-point in terms of dimensionless quantities. Note that β and Γ are related to b_z such that some combinations of the two are impossible (denoted by gray shading in the figure): ${\mathit{b}}_{\mathit{z}}^{\mathrm{2}}\mathrm{=}\frac{\mathrm{(}\mathrm{1}\mathrm{-}\mathrm{\Gamma }\mathit{\beta }\mathrm{)}{\mathit{b}}_{\mathrm{\perp }}^{\mathrm{2}}}{\mathrm{\Gamma }\mathrm{(}\mathrm{1}\mathrm{+}\mathit{\beta }\mathrm{)}}\mathit{,}$(19)which implies Γβ ≤ 1.

Within the accessible parameter space we observe that above a certain level of compressibility (approximately 8, determined empirically), the X-point tends to collapse horizontally and causes the flux rope to descend. Within the range 4.5 < Γ < 8, the X-point collapses vertically, causing the flux rope to move upwards out of equilibrium (though much more slowly than in a driven eruption). However, if the plasma is “stiffened” beyond a threshold, Γ ≲ 3, the fast waves are able to pass through the X-point as their phase speed is no longer close to zero. Since the waves no longer accumulate at the X-point, they do not cause it to collapse and the equilibrium is preserved.

While it is possible that different perturbations might produce different empirical thresholds of stability, it has not been the goal of the present study to determine particular values but rather to show that the X-point stability can be fundamentally related to the accumulation of fast waves at the X-point, which can be moderated by changing the background compressibility of the plasma. Similarly, while the magnitude of the dissipative transport coefficients within the visco-resistive MHD simulations can have some impact on the specific stability thresholds via damping of the fast waves emanating from the flux rope, such damping does not qualitatively change the conclusion of this parameter study.

Fig. 6 Out-of-plane current density jz (color, saturated high and low values), with magnetic flux contours (solid black), in a simulation of flux emergence into an unstratified atmosphere. The four panels show snapshots of the simulation at 100τ = 12 min, 240τ = 29 min, 480τ = 58 min, and 1800τ = 217.5 min.

3.2. CME eruptions driven by flux emergence

The premise of the CS model is that a stable pre-existing flux rope can be driven to eruption by magnetic flux emergence. Flux emergence is achieved through photospheric boundary driving (see Eqs. (15)): a small amount of flux is effectively emerged through the photospheric boundary by applying a time-dependent electric field. Emerging flux can cause the flux rope to move in either direction by forcing a destabilization of the X-point, similarly to the fast waves but more predictably. Within the underlying arcade, when the emergent flux is oppositely oriented to the local flux, it causes a vertical collapse of the X-point, leading to a rising flux rope. Oriented in the same sense as the local flux rope, it causes a horizontal collapse of the X-point, which forces the flux rope downwards. For emergence outside the arcade, likewise, it is possible to choose values for the coefficient c_e in Eqs. (15)such that the simulation results in a vertical X-point collapse, and when the sign of c_e is reversed, the X-point collapses horizontally. However, the sign of c_e must be carefully chosen based on topological and geometric considerations, including the sign of the local overlying flux. In this study, we restrict ourselves to discussing emergence at x₀ = 0 alone, with c_e = 1.1 as in the original CS model¹.

To evaluate the impact of reconnection microphysics, stability of the initial condition and atmospheric stratification on the system’s response to flux emergence in the CS model, the simulation study described below has been performed by changing one model parameter at a time with otherwise identical numerical and dissipation parameters. In the reference simulation with an unstratified atmosphere, the background magnetic “guide” field is set to b_z0 = 1, equivalent to 10 G and of the same order as b_⊥, such that the plasma compressibility measure Γ is less than unity and the initial configuration is stable for any plasma pressure profile. To minimize the impact of the size of the computational domain or the dissipative boundary layers on the results of the simulations, the domain boundaries are placed at L_x = 4 and L_y = 10. The computational grid spanning the (x,y) ∈ [0,L_x] × [0,L_y] domain has 864 and 1536 spatial degrees of freedom in the x and y-directions, respectively, distributed non-uniformly in such a way that the vertically elongated X-point reconnection current sheet is well-resolved in the x-direction, while both magnetic and thermodynamic structures associated with flux emergence through the chromosphere can be well resolved in the y-direction.

The background resistivity throughout the domain is set to η_bg = 10^-5, the photospheric resistivity is set to η_ph = 10^-2, there is no anomalous resistivity , the background kinematic viscosity coefficient is set to μ_bg = 10^-4, and the heat conduction is set to κ = 10^-5. The duration of the flux emergence is taken to be t_e = 300, equivalent to 36.25 min.

Fig. 7 Unstratified atmosphere. Left: height and speed of flux rope center. Smoothing is performed using a Hanning window of 12 points. The speed is computed by finite-differencing the smoothed (blue) curve. Right: temperature at the X-point (center of current sheet) below the flux rope during the same period.

3.2.1. Flux emergence in an unstratified atmosphere

To approximate the coronal conditions in the unstratified simulation, the initial pressure is set to p₀ = 10^-2, such that the initial β is ~1% throughout the domain. To produce an eruption, the photospheric electric field drive is applied within the arcade below the X-point to generate B_x opposite to the magnetic field of the arcade. As a result, the magnetic pressure above the photospheric boundary is reduced causing a local downflow towards the photosphere. This in turn reduces the plasma pressure below the X-point, which forces an in-flow at the sides of the X-point, bringing about its collapse and formation of a reconnection current sheet (e.g., see Fig. 6).

As shown in Fig. 6, the X-point collapse in this simulation is observed to occur at t ≈ 100, forming a current sheet that reaches its maximum length and strength near t = 200. Current density then also increases along the separatrices and the field lines connected to the current sheet. When the new flux stops emerging (t>t_e), the current sheet persists at approximately half to a third of its peak magnitude, slowly diminishing over time for the duration of the simulation.

As reconnection ensues, the flux rope is nudged out of equilibrium (in the $\mathrm{+}\hat{{e}}{}_{\mathit{y}}$ direction) by the reconnection outflow and continues to move outwards as reconnection proceeds. The left panel of Fig. 7 tracks the height of the flux rope center during the eruption by measuring the position of the magnetic O-point (black dots). The height measurements are smoothed (blue curve) using a Hanning window convolution over 12 point windows, and the speed (red curve) is computed by finite-differencing the smoothed height. The maximum speed of the flux rope is observed to be only about 0.7 km s^-1, quickly slowing down further as the reconnection loses steam. In the right panel of Fig. 7, the temperature at the X-point, or the current sheet center, is plotted in mega-Kelvin showing rapid heating early in the eruption due to Joule heating at the current sheet.

We note that this reference simulation results in a very slowly rising flux rope which is inconsistent with the original Chen & Shibata (2000) simulation where the flux rope rise speed of approximately 70 km s^-1 was observed. To study the sensitivity of this result to the magnitude of the background magnetic guide field and the microphysics of reconnection at the X-point, represented here by the anomalous resistivity model similar to that of Chen & Shibata (2000), a series of further simulations has been performed. Figure 8 shows traces of the height of the flux rope center for a set of five simulations with three different values of the guide magnetic field b_z0 = { 1.0,0.5,0.25 } and two resistivity models, one with and another with and j_c = 10, both using the constant background resistivity value η_bg = 10^-5.

Fig. 8 Unstratified atmosphere. Height of flux rope center for a set of five simulations with varying magnitude of initial background magnetic guide field and resistivity models. Three values of the guide magnetic field bz0 = { 1.0,0.5,0.25 } and two resistivity models, one with labeled as “eta const”, and one with and jc = 10 labeled as “eta anom”, both using the constant background resistivity value ηbg = 10-5, are considered.

The comparison of the five traces clearly demonstrates that the outcome of the simulations is much more sensitive to the magnitude of the background guide field, i.e. the global structure and stability of the magnetic configuration, than to the resistivity model. The two traces with b_z0 = 1.0, the initially stable magnetic configuration, are virtually indistinguishable from each other despite very different resistivity models. The two traces with b_z0 = 0.5 initialized from a marginally stable configuration (see Fig. 5) do show small differences during the acceleration phase. Here the simulation with anomalous resistivity allows for slightly faster rise, but both rise much faster than the b_z0 = 1.0 cases. And the initially unstable b_z0 = 0.25 case demonstrates yet faster rise of the flux rope that is comparable to the rise speed observed in the Chen & Shibata (2000) simulation. (Only the anomalous resistivity b_z0 = 0.25 simulation trace is shown in Fig. 8 because the corresponding uniform resistivity simulation produces a very intense X-point current sheet that breaks up due to secondary instabilities (Loureiro et al. 2007), leading to formation of further spatial sub-structure which we have chosen not to attempt to resolve. Detailed investigation of such multi-scale reconnection cases is left for future work.) We note that the choice of critical current density j_c = 10 for onset of anomalous resistivity is such that all five simulations achieve | j | >j_c at the X-point during the acceleration phase of the flux rope, yet that does not result in significant acceleration of the flux rope for the b_z0 = 1.0 and b_z0 = 0.5 cases.

It is also of interest that the rapid rise of the flux rope in the b_z0 = 0.25 case is followed by stagnation at the height of approximately 19Mm. Such stagnation is indicative of the system finding a new stable magnetic equilibrium where the upward force on the flux rope is balanced by the magnetic tension distributed throughout the overlying magnetic arcade.

3.2.2. Flux emergence in a stratified atmosphere

Introduction of atmospheric stratification, as described in Sect. 2.1.3, leads to a more realistic equilibrium plasma configuration that is much denser at the photosphere than in the unstratified corona-like case.

The impact of the flux emergence at the bottom boundary, with and without the atmospheric stratification, is reflected in the traces of height and speed of the respective flux rope eruptions. For the stratified atmosphere, the height and speed of the flux rope as functions of time are shown in the left panel of Fig. 9 and can be compared to the equivalent traces for the reference simulation in the left panel of Fig. 7. (Note the different ranges of the time axes of the two panels.) The two time histories are qualitatively similar, both showing rapid acceleration of the flux-rope during flux emergence, with a reduction of the ejection speed by approximately a factor of two once the driving is turned off. However, both the peak and the post-driving ejection speed of the CME in the stratified atmosphere are less than half of that obtained in the unstratified case.

Fig. 9 Stratified atmosphere. Left: height and speed of flux rope center. Smoothing is performed using a Hanning window of 12 points. The speed is computed by finite-differencing the smoothed (blue) curve. Right: temperature at the X-point (center of current sheet) below the flux rope during the same period.

Another significant difference between the two cases of flux emergence is observed by comparing the time traces of the X-point plasma temperature, shown in the right panels of Figs. 7 and 9. While in the unstratified atmosphere there is a notable temperature increase at the X-point at the time of eruption, in the stratified simulation the temperature decreases instead. Furthermore, as the flux rope begins to rise between 35 and 110 min into the simulation (300 ≲ t ≲ 900) the stratified simulation shows an oscillatory X-point temperature as long as the flux rope is within ≈1 Mm of its original position. The root cause of the overall X-point cooling can easily be explained as upflows of cold chromospheric plasma being advected into the coronal reconnection region. Nevertheless, the observed self-induced quasi-periodic oscillatory behavior of the X-point temperature is somewhat unexpected.

Fig. 10 Evolution in time of the reconnection site behind the CME flux rope in a stratified atmosphere. Each panel shows a snapshot of the temperature structure on the left, the density structure on the right, select contours of the magnetic flux ψ (the same contour values, denoting the same magnetic field lines, have been chosen for each panel), and arrows denoting the in-plane plasma flow. The snapshots are made 348τ = 42 min, 528τ = 64 min, 708τ = 86 min, and 978τ = 118 min into the simulation. Note that for illustration purposes both plasma temperature and number density are plotted using logarithmic color scales with saturated high and low values. Arrows showing the plasma flow have been scaled by a factor of 25 with respect to the linear dimensions of the domain so that an arrow of unit length corresponds to flow of 1.7 × 104 km s-1.

Figure 10 shows the evolution in time of plasma temperature, density, and flows around the X-point during the period of quasi-periodic temperature oscillations. Continuous upflows of dense cool plasma convected along the magnetic field lines and into the reconnection region around the X-point are apparent throughout the evolution. The lower-right panel of this figure makes clear that this continuous chromospheric upflow results in quasi-periodic striations of cool dense material alternating with hotter, lower-density plasma on the recently reconnected field lines rising towards (and with) the flux rope located above. These striations are the signatures of the same oscillatory behavior observed on the X-point temperature trace in Fig. 9. While the origins and parametric robustness of the observed quasi-periodic phenomenon require further in-depth study that is outside of the scope of this article, a heuristic explanation of the basic physical mechanism is straightforward. It results from the competition between the upward directed tension force in newly reconnected magnetic field lines and the downward directed gravity acting on the dense, cold plasma deposited onto these same field-lines by the chromospheric upflows. As in the formation of water droplets, whenever sufficient amount of plasma accumulates in a small enough volume in the V-shaped dip of a set of recently reconnected field lines, the gravitational pull on that plasma overcomes the field’s tension force and a droplet of plasma forms and falls vertically through the reconnection site itself. As a result, those flux-rope destined field lines that produce the droplets end up with lower density hotter plasma, while the field lines that pass through in between the droplets contain colder and heavier plasma. The temperature at the X-point, where the reconnection is regularly disrupted by the droplets, is similarly modulated when the plasma that has been heated by the reconnection process is periodically replaced by the cold plasma of the droplets.

Below the reconnection site, the pattern of chromospheric upflows along the magnetic separatrices and vertical downflows through the X-point creates a circulation of plasma between reconnection’s outflow and inflow. How, and whether or not, this circulation pattern contributes to the formation of the quasi-periodic temperature and density structure described above is left as a topic for future study.

4. Discussion and conclusions

Coronal mass ejections are eruptive solar events of enormous proportions that shed plasma and magnetic flux into interplanetary space. The Chen & Shibata model is a good starting point for understanding how such an eruption can originate from the destabilization of a global magnetic configuration by local flux emergence. It helps us to see a connection between flux emergence, a phenomenon at the solar surface, and flux rope ejection, a phenomenon in the corona. Many observational studies have shown spatio-temporal correlations between flux emergence and eruptive events, but few theoretical models to date have identified a precise single mechanism or sequence of processes whereby producing magnetic flux at the photosphere dynamically triggers an eruption. The CS model may assume an oversimplified solar atmosphere and a somewhat manufactured magnetic topology, but it does proffer a complete story. To determine the effects of a more realistic solar atmosphere, we have undertaken an effort to repeat the study using a different numerical suite and allowing for a stratified atmosphere with the density variation of over four orders of magnitude, as well as a sharp temperature transition between the chromosphere and the corona.

We have found that even in the absence of stratification the initial equilibrium can be unstable to small perturbations. The initial adjustment of the magnetic equilibrium to slight force imbalances can generate fast waves that may not be able to propagate through the X-point below the flux rope. In these cases, the fast waves accumulate in such a way as to collapse the X-point and initiate reconnection. Thus, the equilibrium can be destabilized before any photospheric driving is applied. However, we also found that the stability of the CS equilibrium can be controlled by varying the compressibility of the plasma, which in a two-dimensional system is determined by the combination of thermal pressure and the magnitude of the out-of-plane component of the magnetic field. To quantify this effect, we defined a generalized measure of compressibility Γ and have empirically determined the equilibrium’s stability boundaries in terms of Γ.

When emulating flux emergence by applying an electric field at the photospheric boundary, in the unstratified atmosphere, the results of our simulations are qualitatively similar to those of the original study. However, there are also important differences and new findings. As opposed to the original study, when initialized in a stable configuration, our simulations show little evidence of significant flux rope acceleration or Joule heating associated with the reconnection current sheet. Notably, this result appears to be insensitive to the microphysics of the reconnection region. By varying the magnitude of the background out-of-plane magnetic field component and thus changing the stability of the global magnetic configuration, we also show that flux rope rise speeds comparable to the original result are possible but require an unstable magnetic configuration as the initial condition.

We further show that the microphysics of reconnection is more likely to slow down than to accelerate the flux rope by comparing simulations with and without anomalous resistivity. It is well known that current-dependent anomalous resistivity allows for “fast” magnetic reconnection with only weak dependence on the magnitude of resistivity itself (Malyshkin et al. 2005). Yet, for both initially stable and quasi-stable magnetic configurations, allowing for anomalous resistivity did not result in a substantial increase of the flux rope rise speed. That is, merely allowing for faster reconnection did not lead to faster reconnection and faster flux rope ejection. On the other hand, in magnetic configurations where fast flux-rope ejection is possible, the simulations with low guide field indicate that the inability of the magnetic reconnection process to occur sufficiently fast could limit the rise speed of the flux rope.

In the flux emergence simulations with stable magnetic configuration and realistic atmospheric stratification, the weakness of the X-point heating and the slowness of the ejected flux rope are reproduced, and amplified. In these simulations, changes in the magnetic field structure due to flux emergence generate persistent chromospheric upflows of cold, dense material that is convected into and dramatically cools the reconnection current sheet.

In addition to the steady state upflows and cooling, the stratified simulations also produce another type of behavior: self-induced quasi-periodic oscillations in the X-point temperature, density, and other fluid quantities. The quasi-periodic oscillations observed in the stratified simulation are of transient nature, appearing after the flux emergence drive has been completed and lasting for just over an hour while the flux rope is within ≈1 Mm of its initial location. The robustness of this phenomenon will be a subject of future research, but our initial investigation indicates that a critical balance between the upward tension force of the reconnected magnetic field and the downward gravitational pull on the dense chromospheric plasma convected into the reconnection region has to be achieved in order for the quasi-periodic oscillations to appear in a simulation. While that may seem to be a prohibitive constraint, we speculate that in the three-dimensional parameter space spanned by (1) the height of the X-point; (2) the strength of the magnetic fields and (3) the horizontal location of the emerging flux relative to the separatrices of the pre-existing magnetic configuration, all quantities that can vary greatly throughout the lower solar atmosphere, there is likely embedded a two-dimensional parameter space where such balance can, indeed, be achieved.

We note that there is also extensive observational evidence for what has been called quasi-periodic pulsations (QPP) in solar and stellar flares (e.g., see Nakariakov & Melnikov 2009; Mitra-Kraev et al. 2005, and references therein) with the QPP periodicity time scale varying from fractions of a second to several minutes, comparable to the period of the oscillations produced in our simulation. In fact, Nakariakov & Melnikov (2009) have previously resorted to the water drop formation analogy in describing what they refer to as a class of “load/unload” models of long multi-minute period QPPs. The plasma droplet mechanism described in Sect. 3.2.2 above is a much more direct, and novel, analogy to the same physical process with the potential to provide a new alternative explanation for the long-duration QPPs.

Finally, we point out that the limitations of the two-dimensional MHD model used here for modeling a region of potential flaring activity embedded into a stratified solar atmosphere are many. It is well known that laminar resistive reconnection cannot account for the observed rates of magnetic energy release, particle acceleration, or radiation from solar flares, while three-dimensional effects can substantially alter both the flux-rope stability properties and the microphysics of reconnection. Nevertheless, we believe that the careful and systematic study described in this article is a prerequisite for performing more complete, and also substantially more challenging and complicated, studies of CME initiation by flux emergence in the future.

Appendix A

We derive the uniform pressure magnetic flux rope equilibrium with axial field from a familiar form of the Grad–Shafranov equation: $\frac{\mathrm{d}}{\mathrm{d}\mathit{\psi }}\left(\frac{{\mathit{b}}_{\mathit{z}}^{\mathrm{2}}}{\mathrm{2}}\right)\mathrm{=}\mathrm{-}{\mathrm{\nabla }}^{\mathrm{2}}\mathit{\psi }\mathrm{=}{\mathit{j}}_{\mathit{z}}\mathit{.}$(A.20)In particular, assuming ψ(r) = ψ_l(r) given by Eq. (2a): $\begin{array}{ccc}\frac{{\mathit{b}}_{\mathit{z}}^{\mathrm{2}}}{\mathrm{2}}& \mathrm{=}& \int {\mathit{j}}_{\mathit{z}}\mathit{d\psi }\mathrm{=}\mathrm{\int }{\mathit{j}}_{\mathit{z}}\frac{\mathrm{d}\mathit{\psi }}{\mathrm{d}\mathit{r}}\mathrm{d}\mathit{r}\\ & \mathrm{=}& \int \left[\frac{\mathrm{4}{\mathit{r}}^{\mathrm{2}}}{{\mathit{r}}_{\mathrm{0}}^{\mathrm{3}}}\mathrm{-}\frac{\mathrm{4}}{{\mathit{r}}_{\mathrm{0}}}\right]\left[\frac{\mathit{r}}{{\mathit{r}}_{\mathrm{0}}}\mathrm{-}\frac{\mathit{r}\mathrm{(}{\mathit{r}}^{\mathrm{2}}\mathrm{-}{\mathit{r}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{)}}{{\mathit{r}}_{\mathrm{0}}^{\mathrm{3}}}\right]\mathrm{d}\mathit{r}\\ & \mathrm{=}& \mathrm{-}\frac{\mathrm{2}}{\mathrm{3}}{\left(\frac{\mathit{r}}{{\mathit{r}}_{\mathrm{0}}}\right)}^{\mathrm{6}}\mathrm{+}\mathrm{3}{\left(\frac{\mathit{r}}{{\mathit{r}}_{\mathrm{0}}}\right)}^{\mathrm{4}}\mathrm{-}\mathrm{4}{\left(\frac{\mathit{r}}{{\mathit{r}}_{\mathrm{0}}}\right)}^{\mathrm{2}}\mathrm{+}\mathrm{constant,}\hspace{0.17em}\mathrm{for}\mathit{r}\mathrm{\le }{\mathit{r}}_{\mathrm{0}}\mathit{.}\end{array}$(A.21)Requiring that b_z = b_z0 for r ≥ r₀, we can determine the constant of integration such that b_z is continuous: ${\mathit{b}}_{\mathit{z}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ \sqrt{{\mathit{b}}_{\mathit{z}\mathrm{0}}^{\mathrm{2}}\mathrm{+}\frac{\mathrm{10}}{\mathrm{3}}\mathrm{-}\mathrm{8}{\left(\frac{\mathit{r}}{{\mathit{r}}_{\mathrm{0}}}\right)}^{\mathrm{2}}\mathrm{+}\mathrm{6}{\left(\frac{\mathit{r}}{{\mathit{r}}_{\mathrm{0}}}\right)}^{\mathrm{4}}\mathrm{-}\frac{\mathrm{4}}{\mathrm{3}}{\left(\frac{\mathit{r}}{{\mathit{r}}_{\mathrm{0}}}\right)}^{\mathrm{6}}}\hspace{0.17em}& \mathit{r}\mathrm{\le }{\mathit{r}}_{\mathrm{0}}\\ {\mathit{b}}_{\mathit{z}\mathrm{0}}\mathit{.}& \mathit{r}\mathit{>}{\mathit{r}}_{\mathrm{0}}\mathit{.}\end{array}$(A.22)Suppose, however, we wish to find b_z as a function of ψ, rather than of r. We then solve for the inverse function $\mathit{\zeta }\mathrm{\equiv }{\mathit{\psi }}_{\mathit{l}}^{-1}\mathrm{\left(}\mathit{r}\mathrm{\right)}$ by replacing r with ζ in Eq. (2a) and rearranging terms: ${\mathit{\zeta }}^{\mathrm{4}}\mathrm{-}\mathrm{4}{\mathit{r}}_{\mathrm{0}}^{\mathrm{2}}{\mathit{\zeta }}^{\mathrm{2}}\mathrm{+}{\mathit{r}}_{\mathrm{0}}^{\mathrm{4}}\mathrm{+}\mathrm{4}{\mathit{r}}_{\mathrm{0}}^{\mathrm{3}}{\mathit{\psi }}_{\mathit{l}}\mathrm{=}\mathrm{0.}$(A.23)Solving this quadratic equation for ζ², we find ${\mathit{\zeta }}^{\mathrm{2}}\mathrm{=}\mathrm{2}{\mathit{r}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{\pm }\sqrt{\mathrm{3}{\mathit{r}}_{\mathrm{0}}^{\mathrm{4}}\mathrm{-}\mathrm{4}{\mathit{r}}_{\mathrm{0}}^{\mathrm{3}}{\mathit{\psi }}_{\mathit{l}}}\mathit{.}$(A.24)We recover the form of Eq. (8)by rejecting the positive root (to permit small values of ζ), replacing ψ_l by the full functional form of ψ = ψ_l + ψ_i + ψ_b to approximate a force-free initial condition with well-aligned contours of constant ψ and b_z, and allowing for gauge freedom.

Appendix B

Here we present a derivation of the pressure profile used in simulations of a stratified solar atmosphere, given a particular temperature profile Eq. (10). To be physically relevant, we use here dimensional quantities, rather the normalized code variables.

We begin with the first-order differential equation governing hydrostatic equilibrium: $\frac{\mathrm{d}\mathit{p}}{\mathrm{d}\mathit{y}}\mathrm{+}{\mathit{m}}_{\mathrm{p}}\mathit{n}{\mathit{g}}_{\mathit{S}}\mathrm{=}\mathrm{0}\mathit{,}$(B.25)which we divide by p = 2nk_BT$\frac{\mathrm{d}\ln \mathit{p}}{\mathrm{d}\mathit{y}}\mathrm{+}\frac{{\mathit{m}}_{\mathrm{p}}{\mathit{g}}_{\mathit{S}}}{\mathrm{2}{\mathit{k}}_{\mathrm{B}}\mathit{T}}\mathrm{=}\mathrm{0.}$(B.26)Then $\ln \frac{\mathit{p}}{{\mathit{p}}_{\mathrm{0}}}\mathrm{=}\mathrm{-}\frac{{\mathit{m}}_{\mathrm{p}}{\mathit{g}}_{\mathit{S}}}{\mathrm{2}{\mathit{k}}_{\mathrm{B}}}\mathrm{\int }\frac{\mathrm{d}\mathit{y}}{\mathit{T}}\mathrm{·}$(B.27)We use the profile for temperature T given by Eq. (10), but with the following variable substitution: $\mathit{u}\mathrm{\equiv }\frac{\mathit{y}\mathrm{-}{\mathit{y}}_{\mathrm{TR}}}{\mathrm{\Delta }\mathit{y}}\mathit{,}$(B.28)leading to $\begin{array}{ccc}\mathit{T}\mathrm{\left(}\mathit{u}\mathrm{\right)}& \mathrm{=}& {\mathit{T}}_{\mathrm{p}}\mathrm{+}\frac{{\mathit{T}}_{\mathrm{c}}\mathrm{-}{\mathit{T}}_{\mathrm{p}}}{\mathrm{2}}\left(\mathrm{1}\mathrm{+}\tanh \mathit{u}\right)\\ & \mathrm{=}& {\mathit{T}}_{\mathrm{p}}\mathrm{+}\frac{{\mathit{T}}_{\mathrm{c}}\mathrm{-}{\mathit{T}}_{\mathrm{p}}}{\mathrm{2}}\left(\mathrm{1}\mathrm{+}\frac{{\mathrm{e}}^{\mathit{u}}\mathrm{-}{\mathrm{e}}^{\mathrm{-}\mathit{u}}}{{\mathrm{e}}^{\mathit{u}}\mathrm{+}{\mathrm{e}}^{\mathrm{-}\mathit{u}}}\right)\\ & \mathrm{=}& \end{array}$(B.29)With algebraic manipulations, we can rewrite Eq. (B.29) as $\frac{\mathrm{1}}{\mathit{T}}\mathrm{=}\frac{\mathrm{1}}{{\mathit{T}}_{\mathrm{c}}}\mathrm{+}\frac{{\mathrm{e}}^{\mathrm{-}\mathrm{2}\mathit{u}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{T}}_{\mathrm{p}}\mathit{/}{\mathit{T}}_{\mathrm{c}}\mathrm{)}}{{\mathit{T}}_{\mathrm{p}}\hspace{0.17em}{\mathrm{e}}^{\mathrm{-}\mathrm{2}\mathit{u}}\mathrm{+}{\mathit{T}}_{\mathrm{c}}}\mathrm{·}$(B.30)Then the integral in Eq. (B.27)can be evaluated $\begin{array}{ccc}\mathrm{\int }\frac{\mathrm{d}\mathit{y}}{\mathit{T}}& \mathrm{=}& \mathrm{\Delta }\mathit{y}\mathrm{\int }\frac{\mathrm{d}\mathit{u}}{\mathit{T}}\mathrm{=}\mathrm{\Delta }\mathit{y}\left[\mathrm{\int }\frac{\mathrm{d}\mathit{u}}{{\mathit{T}}_{\mathrm{c}}}\mathrm{+}\left(\mathrm{1}\mathrm{-}\frac{{\mathit{T}}_{\mathrm{p}}}{{\mathit{T}}_{\mathrm{c}}}\right)\mathrm{\int }\frac{{\mathrm{e}}^{\mathrm{-}\mathrm{2}\mathit{u}}\hspace{0.17em}\mathrm{d}\mathit{u}}{{\mathit{T}}_{\mathrm{p}}\hspace{0.17em}{\mathrm{e}}^{\mathrm{-}\mathrm{2}\mathit{u}}\mathrm{+}{\mathit{T}}_{\mathrm{c}}}\right]\\ & \mathrm{=}& \end{array}$(B.31)Finally, substituting Eq. (B.31) into Eq. (27)yields an expression for p, in terms of u: $\mathit{p}\mathrm{\left(}\mathit{u}\mathrm{\right)}\mathrm{=}{\mathit{p}}_{\mathrm{0}}\exp \left\{\frac{{\mathit{m}}_{\mathrm{p}}{\mathit{g}}_{\mathit{S}}\mathrm{\Delta }\mathit{y}}{\mathrm{2}{\mathit{k}}_{\mathrm{B}}{\mathit{T}}_{\mathrm{c}}}\left[\frac{{\mathit{T}}_{\mathrm{c}}\mathrm{-}{\mathit{T}}_{\mathrm{p}}}{\mathrm{2}{\mathit{T}}_{\mathrm{p}}}\hspace{0.17em}\ln \mathrm{(}{\mathit{T}}_{\mathrm{p}}{\mathrm{e}}^{\mathrm{-}\mathrm{2}\mathit{u}}\mathrm{+}{\mathit{T}}_{\mathrm{c}}\mathrm{)}\mathrm{-}\mathit{u}\right]\right\}\mathit{.}$(B.32)

Acknowledgments

E.L. thanks Neil Sheeley for his valuable insights into solving the hyperbolic function integrals. This work was supported by the NASA SR&T and LWS programs, as well as ONR 6.1 program. Simulations were performed under grants of computer time from the US DOD HPC program and the National Energy Research Scientific Computing Center, which is supported by the US DOE Office of Science.

References

All Tables

All Figures

	Fig. 1 Contours of ψ in the initial conditions.
In the text

	Fig. 2 Contours of ψ (black) and bz (magenta) at t = 0 for bz as a function of r (left) and as a function of ψ (right).
In the text

	Fig. 3 Logarithm (base 10) of normalized number density n (green), normalized pressure p (blue), and temperature T (red) as a function of height in the initial conditions for a stratified atmosphere.
In the text

	Fig. 4 No driving. Four snapshots of current density jz showing outward propagating MHD waves. Notice the trapping and interference of waves at the X-point; compounding of these waves here precipitates an X-point collapse, leading to the formation of a current sheet and initiation of magnetic reconnection.
In the text

	Fig. 5 Dependence of flux rope stability on the free parameters p0 and bz0 (left), as well as on the plasma β and compressibility measure Γ (right). All simulations are performed without driving. The different symbols signify that the flux rope was stable (black dots); was wobbly but on average did not rise or sink (black stars); moved upwards (blue triangles); or moved downwards (red triangles).
In the text

	Fig. 6 Out-of-plane current density jz (color, saturated high and low values), with magnetic flux contours (solid black), in a simulation of flux emergence into an unstratified atmosphere. The four panels show snapshots of the simulation at 100τ = 12 min, 240τ = 29 min, 480τ = 58 min, and 1800τ = 217.5 min.
In the text

	Fig. 7 Unstratified atmosphere. Left: height and speed of flux rope center. Smoothing is performed using a Hanning window of 12 points. The speed is computed by finite-differencing the smoothed (blue) curve. Right: temperature at the X-point (center of current sheet) below the flux rope during the same period.
In the text

Fig. 8 Unstratified atmosphere. Height of flux rope center for a set of five simulations with varying magnitude of initial background magnetic guide field and resistivity models. Three values of the guide magnetic field bz0 = { 1.0,0.5,0.25 } and two resistivity models, one with labeled as “eta const”, and one with and jc = 10 labeled as “eta anom”, both using the constant background resistivity value ηbg = 10-5, are considered.

In the text

	Fig. 9 Stratified atmosphere. Left: height and speed of flux rope center. Smoothing is performed using a Hanning window of 12 points. The speed is computed by finite-differencing the smoothed (blue) curve. Right: temperature at the X-point (center of current sheet) below the flux rope during the same period.
In the text

Fig. 10 Evolution in time of the reconnection site behind the CME flux rope in a stratified atmosphere. Each panel shows a snapshot of the temperature structure on the left, the density structure on the right, select contours of the magnetic flux ψ (the same contour values, denoting the same magnetic field lines, have been chosen for each panel), and arrows denoting the in-plane plasma flow. The snapshots are made 348τ = 42 min, 528τ = 64 min, 708τ = 86 min, and 978τ = 118 min into the simulation. Note that for illustration purposes both plasma temperature and number density are plotted using logarithmic color scales with saturated high and low values. Arrows showing the plasma flow have been scaled by a factor of 25 with respect to the linear dimensions of the domain so that an arrow of unit length corresponds to flow of 1.7 × 104 km s-1.

In the text