Recurrent CME-like Eruptions in Emerging Flux Regions. I. On the Mechanism of Eruptions

The formation of active regions (ARs) is often associated with the emergence of magnetic flux from the solar interior (e.g., Parker 1955). Many explosive phenomena observed on the Sun, such as flaring events and coronal mass ejections (CMEs), are associated with ARs. In fact, it has been observed that a single AR can produce several CMEs in a recurrent manner (e.g., Nitta & Hudson 2001; Zhang et al. 2008; Wang et al. 2013).

Solar eruptions have been studied extensively in the past. Observational studies have reported on the pre-eruptive phase of the eruption (e.g., Canou & Amari 2010; Vourlidas et al. 2012; Syntelis et al. 2016), the triggering of the eruptions (e.g., Zuccarello et al. 2014; Chintzoglou et al. 2015; Reeves et al. 2015), and the propagation of the erupting structures in the interplanetary medium (e.g., Colaninno et al. 2013) and toward the Earth (e.g., Patsourakos et al. 2016).

Often, eruptions are associated with the formation of a twisted magnetic field structure, which is commonly referred to as a magnetic flux rope (FR; e.g., Cheng et al. 2011; Green et al. 2011; Zhang et al. 2012; Patsourakos et al. 2013). Still, various aspects regarding the process of formation, destabilization, and eruption of FRs are up for debate.

Numerical models studying the formation of magnetic FRs in the solar atmosphere have extensively demonstrated the role of shearing, rotation, and reconnection of field lines in the buildup of magnetic twist. As an example, magnetic flux emergence experiments (e.g., Magara & Longcope 2001; Archontis & Török 2008; Fan 2009) have shown that shearing motions along a polarity inversion line (PIL) can lead to the reconnection of sheared field lines and the gradual formation of FRs, which may erupt in a confined or ejective manner (e.g., Archontis & Hood 2012). Furthermore, experiments where rotational motions are imposed at the photospheric boundary (symmetric and asymmetric driving of polarities; DeVore & Antiochos 2008; Aulanier et al. 2010) have shown that the shearing motions can form a pre-eruptive FR and destabilize the system.

Once an FR is formed, it may erupt in an ejective manner toward outer space (e.g., Leake et al. 2014) or remained confined, for instance, due to a strong overlying field (e.g., Leake et al. 2013). There are two main proposed mechanisms that might be responsible for the triggering and/or driving of the eruption of magnetic FRs. One is the non-ideal process of magnetic reconnection and the other is the action of an ideal magnetohydrodynamic (MHD) instability.

One example of reconnection that leads to the eruption of a magnetic FR is the well-known tether-cutting mechanism (Moore & Labonte 1980; Moore & Roumeliotis 1992). During this process, the footpoints of sheared field lines reconnect along a PIL, forming an FR. The FR slowly rises, dragging in the magnetic field from the sides, and a current sheet is formed underneath the FR. Eventually, fast reconnection of the field lines that envelop the FR occurs at the current sheet. Then, the upward reconnection outflow assists the further rise of the FR. In this way, an imbalance is achieved between (a) the upward magnetic pressure and tension force and (b) the downward tension force of the envelope field lines. This leads to an ejective eruption of the FR. Another example is the so-called breakout reconnection between the envelope field and a pre-existing magnetic field. If the relative orientation of the two fields is antiparallel, (external) reconnection between them becomes very effective when they come into contact (e.g., Antiochos et al. 1999; Archontis & Hood 2012; Karpen et al. 2012; Leake et al. 2014). This reconnection releases the downward magnetic tension of the envelope field, and the FR can "break out," experiencing an ejective eruption. We should highlight that the relative orientation and field strengths of the interacting magnetic systems are important parameters that affect the eruption of the FR. In previous studies, it has been shown that depending on the value of these parameters, the rising FR could experience an ejective eruption or be confined by the envelope field or even become annihilated through the interaction with the pre-existing magnetic field (e.g., Galsgaard et al. 2007; Archontis & Hood 2012; Leake et al. 2014).

Solar eruptions can also be triggered by ideal processes, for example, the helical kink instability (Anzer 1968; Török et al. 2004), which occurs when the twist of the FR exceeds a critical value that depends on the configuration of the FR (e.g., cylindrical, toroidal) and the line-tying effect (e.g., Hood & Priest 1981; Török et al. 2004). During the instability, the axis of the rising FR develops a helical shape. The eruption of the helical magnetic field could be ejective or confined, depending, e.g., on how strong the overlying magnetic field is (Török & Kliem 2005).

Another crucial parameter that affects the eruption of an FR is how the external constraining magnetic field drops along the direction of height. This is related to the so-called torus instability (Bateman 1978; Kliem & Török 2006). In this model, a toroidal current ring with major radius R is placed inside an external magnetic field. This external magnetic field drops along the direction of the major radius as ${R}^{-n}$. Due to the current ring's curvature, a hoop force acts on the current ring. This force is directed away from the center of the torus. An inwards Lorentz force acts on the current ring due to the external magnetic field. Previous studies (Bateman 1978; Kliem & Török 2006) showed that, if the decrease rate of the external field (i.e., $n=-\partial {B}_{\mathrm{external}}/\partial \mathrm{ln}R$) exceeds a critical value (${n}_{\mathrm{crit}}=1.5$), the current ring becomes unstable. The decrease rate of the external field is commonly referred to as the torus or decay index.

The range of values of the critical torus index is still under debate. For instance, studies of emerging flux tubes with an initial arch-like configuration have reported higher values of the torus index ($n=1.7\mbox{--}2$; Fan & Gibson 2007; Fan 2010). An & Magara (2013), in a flux emergence simulation of a straight, horizontal flux tube, reported values of the torus index well above 2. Démoulin & Aulanier (2010) found that the torus index can vary depending on a range of parameters, such as the thickness of the current channel (the axial current of a twisted FR is a current channel). In cases of thin current channels, the index was found to be 1 (1.5) for straight (circular) channels. Also, the expansion of the FR during its eruption affects the critical value of the torus instability. For thick channels, the critical index for circular and straight channels does not vary much—it takes values ranging from 1.1 to 1.3 (with expansion of the FR) and 1.2–1.5 (without expansion). Zuccarello et al. (2015) investigated the role of line-tying effects on the eruption. They performed a series of simulations with a setup similar to that of Aulanier et al. (2010), but with different velocity drivers at the photosphere. They found that the critical index did not depend greatly on the pre-eruptive photospheric motions, and it was found to take values within the range of 1.1–1.3.

In our paper, we show the results of a simulation of magnetic flux emergence, which occurs dynamically from the solar interior to the outer solar atmosphere. We focus on the formation of magnetic FRs in the emerging flux region and their possible eruption. In particular, we show how reconnection leads to the formation of FRs and how/why these FRs erupt. We find that the emergence of a single subphotospheric magnetic flux tube can drive recurrent eruptions, which are produced due to the combined action of the torus instability and reconnection of the envelope field lines in a tether-cutting manner. We find that, at least in the first eruption, the fast ejection phase of the torus-unstable FR is triggered by tether-cutting reconnection. A geometrical extrapolation of the size of the eruptions showed that they can develop into large-scale structures, with a size comparable to small CMEs. The plasma density and temperature distributions reveal that the structure of the erupting fields consist of three main parts: a "core," a "cavity," and a "front edge," which is reminiscent of the "three-part" structure of CMEs. We find that the plasma, at the close vicinity of the "core," is hotter and denser when the envelope field lines reconnect with themselves in a tether-cutting manner during the eruption. The same area appears to be cooler and less dense when the envelope field lines reconnect with some other neighboring (e.g., sheared J-like) field lines.

In Section 2, we describe the initial conditions of our simulations. Section 3.1 is an overview of the dynamics occurring in our simulation leading to four recurrent eruptions. In Section 3.2, we show the morphology of the magnetic field (before, during, and after the eruptions) and the triggering mechanism of these eruptions. In Section 3.3, we show the distribution of the various properties of the erupting fields, such as the density, temperature, velocity, and current profiles. In Section 3.4, we perform an extrapolation of the size of the erupting structures. In Section 4, we summarize the results.

To perform the simulations, we numerically solve the three-dimensional (3D) time-dependent, resistive, compressible MHD equations in Cartesian geometry using the Lare3D code of Arber et al. (2001). The equations in dimensionless form are

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}+{\rm{\nabla }}\cdot (\rho {\boldsymbol{v}})=0,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}\displaystyle \frac{\partial (\rho {\boldsymbol{v}})}{\partial t} & = & -{\rm{\nabla }}\cdot (\rho {\boldsymbol{vv}})+({\rm{\nabla }}\times {\boldsymbol{B}})\times {\boldsymbol{B}}-{\rm{\nabla }}P+\rho {\boldsymbol{g}}\\ & & +{\rm{\nabla }}\cdot {\boldsymbol{S}},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial (\rho \epsilon )}{\partial t}=-{\rm{\nabla }}\cdot (\rho \epsilon {\boldsymbol{v}})-P{\rm{\nabla }}\cdot {\boldsymbol{v}}+{Q}_{\mathrm{joule}}+{Q}_{\mathrm{visc}},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t}={\rm{\nabla }}\times ({\boldsymbol{v}}\times {\boldsymbol{B}})+\eta {{\rm{\nabla }}}^{2}{\boldsymbol{B}},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\epsilon =\displaystyle \frac{P}{(\gamma -1)\rho },\end{eqnarray} \tag{ 5 }$$

where ρ, ${\boldsymbol{v}}$, ${\boldsymbol{B}}$, and P are the density, velocity vector, magnetic field vector, and gas pressure. Gravity is included. We assume a perfect gas with a specific heat of $\gamma =5/3$. Viscous heating Q_visc and Joule dissipation Q_joule are also included. We use an explicit anomalous resistivity that increases linearly when the current density exceeds a critical value J_c:

$$\begin{eqnarray}\eta =\left\{\begin{array}{ll}{\eta }_{b}, & \mathrm{if}\ | J| \lt {J}_{c}\\ {\eta }_{b}+{\eta }_{0}\left(\tfrac{| J| }{{J}_{c}}-1\right), & \mathrm{if}\ | J| \gt {J}_{c}\end{array}\right.,\end{eqnarray} \tag{ 6 }$$

where ${\eta }_{b}=0.01$ is the background resistivity, J_c = 0.005 is the critical current, and ${\eta }_{0}=0.01$.

We use a normalization based on the photospheric values of density ${\rho }_{{\rm{c}}}=1.67\times {10}^{-7}\ {\rm{g}}\ {\mathrm{cm}}^{-3}$, length ${H}_{{\rm{c}}}=180\ \mathrm{km}$, and magnetic field strength ${B}_{{\rm{c}}}=300\ {\rm{G}}$. From these, we get the pressure ${P}_{{\rm{c}}}=7.16\times {10}^{3}\ \mathrm{erg}\ {\mathrm{cm}}^{-3}$, temperature ${T}_{{\rm{c}}}=5100\,{\rm{K}}$, velocity ${v}_{0}=2.1\ \mathrm{km}\ {{\rm{s}}}^{-1}$, and time ${t}_{0}=85.7\ {\rm{s}}$.

The computational box has a size of $64.8\times 64.8\times 64.8\ \mathrm{Mm}$ in the x-, y-, z-directions, in a $417\times 417\times 417$ grid. We assume periodic boundary conditions in the y-direction. Open boundary conditions are at the two yz plane boundaries and at the top of the numerical box. The domain consists of an adiabatically stratified subphotosheric layer at $-7.2\ \mathrm{Mm}\leqslant z\lt 0\ \mathrm{Mm}$, an isothermal photospheric–chromospheric layer at $0\ \mathrm{Mm}\leqslant z\lt 1.8\ \mathrm{Mm}$, a transition region at $1.8\ \mathrm{Mm}\leqslant z\lt 3.2\ \mathrm{Mm}$, and an isothermal coronal at $3.2\ \mathrm{Mm}\leqslant z\lt 57.6\ \mathrm{Mm}$. We assume a field-free atmosphere in hydrostatic equilibrium. The initial distribution of the temperature (T), density (ρ), and gas (P_g) pressure is shown in Figure 1.

Zoom In Zoom Out Reset image size

We place a straight, horizontal FR at $z=-2.1\ \mathrm{Mm}$. The axis of the FR is oriented along the y-direction, so the transverse direction is along x and the height is in the z-direction. The magnetic field of the FR is

$$\begin{eqnarray}&&{B}_{y}={B}_{0}\exp (-{r}^{2}/{R}^{2}),\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{B}_{\phi }=\alpha {{rB}}_{y},\end{eqnarray} \tag{ 8 }$$

where R = 450 km is a measure of the FR's radius, r is the radial distance from the FR's axis, and $\alpha =0.4$ (0.0023 km⁻¹) is a measure of twist per unit of length. The magnetic field's strength is ${B}_{0}=3150\,{\rm{G}}$. Its magnetic pressure (P_m) is overplotted in Figure 1. Initially, the FR is in pressure equilibrium. The FR is destabilized by imposing a density deficit along its axis, similar to the work by Archontis et al. (2004):

$$\begin{eqnarray}&&{\rm{\Delta }}\rho =\displaystyle \frac{{p}_{{\rm{t}}}(r)}{p(z)}\rho (z)\exp (-{y}^{2}/{\lambda }^{2}),\end{eqnarray} \tag{ 9 }$$

where p is the external pressure and p_t is the total pressure within the FR. The parameter λ is the length scale of the buoyant part of the FR. We use $\lambda =5$ (0.9 Mm).

3.1. Overall Evolution: A Brief Overview

In the following, we briefly describe the overall evolution of the emerging flux region during the running time of the simulation. At t = 25 minutes, the crest of the subphotospheric FR reaches the photosphere. It takes 10 minutes for the magnetic buoyancy instability criterion (see Acheson 1979; Archontis et al. 2004) to be satisfied and thus, for the first magnetic flux elements to emerge at and above the solar surface. Eventually, the emerging magnetized plasma expands as it rises, due to the magnetic pressure inside the tube and the decreasing gas pressure of the background stratified atmosphere. Because of the expansion, the outermost expanding field lines adopt a fan-like configuration, forming an envelope field that surrounds all of the upcoming magnetized plasma. As we discuss later in this paper, the characteristics and dynamical evolution of this envelope field play an important role in understanding the eruptions coming from the emerging flux region.

At the photosphere, the emergence of the field forms a bipolar region with a strong PIL. Similarly to previous studies (e.g., Manchester 2001; Archontis & Török 2008; Leake et al. 2013), we find that the combined action of shearing, driven by the Lorentz force along the PIL, and the reconnection of the sheared field lines leads to the formation of a new magnetic FR, which eventually erupts toward outer space. In fact, this is an ongoing process, which leads to the formation and eruption of several FRs during the evolution of the system. Since these FRs are formed after the initial flux emergence at the photosphere, we will refer to them as post-emergence FRs.

Figure 2 shows the temperature distribution (vertical xz-midplane) and selected field lines at the times of four successive eruptions in our simulation (panels (a)–(d)). The temperature distribution delineates the (bubble-shaped) volume of the erupting field, which is filled by cool and hot plasma. In Section 3.3, we discuss the physical properties (e.g., temperature, density) of the erupting plasma in more detail. The field lines are drawn in order to show a first view of the shape of the envelope field (green) and the core of the erupting FRs (yellow). Notice the strongly azimuthal nature of the envelope field and the S-shaped configuration of the FR's field lines in the first eruption (Figure 2(a), top view). In the following eruptions, the orientation of the envelope field changes (in a counterclockwise manner; Figures 2(b)–(d), top view). The morphology of the field lines during the four eruptions is discussed in detail in Section 3.2. We find that all of the eruptions are fully ejective (i.e., they exit the numerical domain from the top boundary).

Zoom In Zoom Out Reset image size

To further describe the overall dynamical evolution of the eruptions, we calculate the total magnetic and kinetic energies (black and red lines, respectively; Figure 3) above the mid-photosphere (z = 1.37 Mm). The first maximum of the kinetic energy at $t=45.7\,\mathrm{minutes}$ corresponds to the initial emergence of the field. Then, we find four local maxima of the magnetic and kinetic energies, which correspond to the four eruptions (e.g., kinetic energy peaks at $t=74.3,85.7,117.1,194.3\,\mathrm{minutes}$, marked by the vertical lines in the figure). As expected, the magnetic (kinetic) energy decreases (increases) after each eruption. Notice that this is less pronounced for the magnetic energy in the first eruption because of the continuous emergence of magnetic flux, which increases the total amount of magnetic energy above the mid-photosphere. Also, the local maximum of the kinetic energy at $t=205.7\,\mathrm{minutes}$ corresponds mainly to the fast reconnection upflow underneath the erupting FR, which is about to exit the numerical domain.

Zoom In Zoom Out Reset image size

In a similar way, we compute the self-helicity (Figure 4). For a single twisted flux tube, the self-helicity is assumed to correspond to the twist within the flux tube. For the calculation, we used the method described in Moraitis et al. (2014). Overall, we find that the temporal evolution of the self-helicity is similar to that of the kinetic energy (e.g., they reach local maxima at the same time), which indicates that the erupted field is twisted. We also find that between eruptions, self-helicity increases because of the gradual build up of the twist of the post-emergence FRs.

Zoom In Zoom Out Reset image size

3.2. Flux Rope Formation and Eruption Mechanisms

3.2.1. First Eruption

The formation of the post-emergence FR occurs in the low atmosphere due to the combination of (a) shearing and converging motions along the PIL, (b) rotation of the polarities of the emerging flux region, and (c) reconnection of the sheared and rotated field lines along the PIL.

First, we would like to focus on the role of shearing along the PIL and the rotation of the polarities during the pre-eruptive phase. For this reason, we present a side view (Figures 5(a), (b)) and a top view (Figures 5(c), (d)) of a close-up of the emerging flux region. We plot the sheared arcade field lines (blue), the $| J/B|$ isosurface, and the photospheric B_z component of the magnetic field (black/white plane). On the photospheric plane, we also plot the planar component of the velocity field vector (yellow arrows) and the ${\omega }_{z}$ component of the vorticity (red contours). The visualization of the velocity field reveals (a) the shearing motion along the PIL (the yellow arrows are almost antiparallel on the two sides of the PIL) and (b) the converging motions toward the PIL and close to the two main polarities, due to their rotation. These motions (shear and rotation) are also apparent by looking at the vertical component of the vorticity (red contours). Notice that ${\omega }_{z}$ is strong close to the two polarities, where the rotation is fast. Along and sideways of the PIL, there is only apparent "vorticity" due to the velocity, which is developed by the shearing.

Zoom In Zoom Out Reset image size

The footpoints of the sheared arcade field lines are rooted in both sides of the PIL (e.g., blue lines in Figures 5(a) and (c)). Due to the shearing, their footpoints move toward the two polarities where they undergo rotation (e.g., see the footpoints of the blue field lines, which go through the red contours close to the two opposite polarities; Figures 5(b) and (d)). Due to the rotation, the sheared field lines adopt the characteristic hook-shaped edge, forming J-like loops. The isosurface of high values of $| J/B|$ shows the formation of a strong current between the J-like loops. When the J-like field lines reconnect at the current sheet, new twisted field lines, with an overall sigmoidal shape, are formed.

Figure 6 is a visualization of a series of selected field lines during the slow-rise (panels (a) and (b)) and the fast-rise (panels (c)–(f)) phase of the first eruption. In a similar manner to Figure 5, Figure 6(a) shows the sheared field lines (blue) and the $| J/B|$ isosurface (purple). Reconnection between the sheared field lines forms a new set of longer field lines (yellow), which connect the distant footpoints of the sheared field lines. Thus, the longer field lines produce a magnetic loop above the PIL. As time goes on (panel (b)), further reconnection between the J-like sheared field lines (blue) forms another set of field lines, which wrap around the magnetic loop, producing the first (post-emergence) magnetic FR. The red and green field lines are not reconnected field lines. They were traced from arbitrary heights above the yellow field lines. They belong to the emerging field, which has expanded into the corona. In that respect, they create an envelope field for the new magnetic FR.

Zoom In Zoom Out Reset image size

Eventually, the envelope field lines just above the FR (e.g., red lines; Figure 6(b)) are stretched vertically, and their lower segments come into contact and reconnect at the flare current sheet underneath the FR in a tether-cutting manner. Hereafter, for simplicity, we call the reconnection between envelope field lines the EE-TC reconnection (i.e., Envelope–Envelope–Tether-Cutting reconnection). This reconnection occurs in a fast manner, triggering an explosive acceleration of the FR. During this process, the plasma temperature at the flare current sheet reaches values up to 6 MK. The rapid eruption is followed by a similar type of reconnection of the outermost field lines of the envelope field (green lines, Figure 6(c)). Figures 6(c), (d), and (e) show the side, front, and top views of the field-line morphology at t = 74 minutes. Figure 6(f) is a close-up of the reconnection site underneath the erupting FR.

Notice that, due to EE-TC reconnection, the red field lines are wrapped around the central region of the erupting field (yellow field lines). They make at least two turns around the axis, becoming part of the erupting FR. During the eruption, these field lines may reconnect more than once, and thus, have more than two full turns around the axis. The close-up in Figure 6(f) shows that a post-reconnection arcade (light blue field lines) is formed below the flare current sheet. At the top of the arcade, the plasma is compressed and the temperature increases up to 10 MK.

The time evolution of the post-emergence FR can be followed by locating its axis at different times. To find the axis, we use a vertical 2D cut (at the middle of the FR, along its length), which is perpendicular to the field lines of the FR. Then, we locate the maximum of the normal component of the magnetic field (B_n) on this 2D plane for every snapshot. We also found that the location of the axis of the FR is almost identical to the location of the maximum plasma density within the central region of the FR. The latter can be used as an alternative tracking method for the location of the FR's axis. Using the above method(s), we are able to plot the height–time profile of the erupting FR (see Figure 7(a) black line) and its derivative (blue line). The height–time profile shows a phase of gradual upward motion (slow-rise phase), which is followed by an exponential period (fast-rise phase). The terminal velocity before the FR exits the numerical box is 170 $\mathrm{km}\,{{\rm{s}}}^{-1}$. During the eruptive phase, the FR is not very highly twisted and it also does not have the characteristic deformation of its axis that results from the kink instability. As a result, kink instability does not seem to play a role in this case. To study whether torus instability is at work, we follow the torus index calculation method of Fan & Gibson (2007) and Aulanier et al. (2010). We first estimate the external (envelope) field by calculating the potential magnetic field (B_p). This is done based on the calculations made to derive the helicity (details in Moraitis et al. 2014). To solve the Laplace equation for the calculation of the potential field, it is assumed that both the magnetic field and the potential field have the same normal component at the boundaries (Neumann conditions). The lower xy-plane boundary is the photosphere at z = 0.51 Mm, and the rest of the boundaries are the sides of the numerical domain. Having calculated B_p, we then compute the torus index as $n=-z\partial \mathrm{ln}{B}_{p}/\partial z$. Then, we find the value of the torus index at the position of the FR center by measuring the value of n along the height–time profile. We plot the results in Figure 7(b).

Zoom In Zoom Out Reset image size

From the height–time profile (black line and inset in Figure 7(a)), we find that the FR enters an exponential rise phase just after $t=61.4\,\mathrm{minutes}$ (first vertical line). The torus index at this time is n = 1.81, which lies within the estimated range of values for the occurrence of the torus instability (see the Introduction and the highlighted region in Figure 7(e)). Therefore, we anticipate that the FR in our simulation becomes torus unstable at $t\geqslant 61.4\,\mathrm{minutes}$.

We should highlight that the envelope field lines above the FR start to reconnect in a TC manner at $t\geqslant 67.9\,\mathrm{minutes}$ (second vertical line, Figure 7). As a result, the mean tension of the envelope field lines (Figure 7(c)) decreases while the FR height and velocity increase dramatically (Figure 7(a)). We also find that the fast reconnection jet (V_z up to 550 $\mathrm{km}\,{{\rm{s}}}^{-1}$), which is ejected upward from the flare current sheet, is transferring momentum to the FR and contributes to its acceleration (Figure 7(d)).

3.2.2. Second Eruption

In the following, we focus on the dynamics of the second eruption. Figures 8(a) and (b) are close-ups of the area underneath the first erupting FR at t = 74 minutes and t = 79 minutes, respectively. In a similar manner to the formation of the first FR, the second FR (yellow field lines, Figure 8(b)) is formed due to reconnection between J-loops (blue field lines). The post-reconnection arcade (green and red field lines in Figures 8(a) and (b)), which was formed after the first eruption (cyan lines, Figure 6(f)), overlies the yellow field lines and thus, it acts as an envelope field for the second FR. Above and around this envelope field, there are field lines (gray) that belong to the first eruptive flux system but have not exited the numerical domain yet. Hereafter, we refer to this field as the external, pre-existing field.

Zoom In Zoom Out Reset image size

As the second post-emergence FR moves upwards, the envelope field lines are stretched vertically and their footpoints move toward the current sheet (pink isosurface). However, they do not reconnect in an EE-TC manner. Instead, the lower segments of the envelope field lines reconnect with the J-like loops. Hereafter, for simplicity, we refer to this as EJ-TC reconnection (i.e., Envelope–J–Tether-Cutting reconnection). This difference is due to the different orientation of the envelope field lines. As we previously showed (green lines, Figure 2(b), top view), the envelope field lines in the second eruption do not have a strongly azimuthal nature. They are mainly oriented along the y-direction. Therefore, their lower segments come closer to the J-like loops and reconnect with them (e.g., bottom-right red lines, Figure 8(c)).

To better illustrate the EJ-TC reconnection, in Figure 8(d) we show a close-up of this region. Here, the envelope field lines (green) reconnect with the J-like loops (blue) to form the hook-shaped field lines (red). Eventually, this process occurs in both footpoints of the envelope field lines, forming new field lines such as the red ones in Figure 8(e). Notice that these new reconnected field lines are winding around the footpoints of the rising FR and therefore, they become part of the erupting field. In general, the EJ-TC reconnection removes flux from the envelope field and adds flux to the FR. Also, the downward tension of the envelope field decreases during EJ-TC reconnection. Before the FR exits the box (Figure 8(f)), most of the envelope field has been subject to EJ-TC reconnection. We should emphasize that we do not find evidence of EE-TC reconnection during the second eruption.

The EJ-TC and EE-TC reconnections produce field lines with a different shape. In the first eruption, the EE-TC reconnected field lines (red, Figures 6(c)–(e)) are ejected toward the FR center, adopting a "W-shaped" configuration. The concave-upward segments of the W-like field lines (arrows, Figure 6(d)) bring dense plasma from the low atmosphere into the central region of the FR. In the second eruption, the EJ-TC reconnected field lines have hook-like segments in their footpoints (arrows, Figure 8(f)). In this case, the tension in the reconnected field lines ejects hot and dense plasma sideways (mainly along the y-direction) and not toward the center of the FR. Thus, due to the different way that the envelope field lines reconnect, the temperature and density distributions within the erupting field show profound differences between the first and subsequent eruptions. This is discussed in more detail in Section 3.3.

We plot now the height–time profile and its derivative for the second FR (black and blue lines, Figure 9(a)). To calculate the torus index, we consider the potential magnetic field B_p. As discussed earlier, the calculation of the potential field takes into account all of the boundaries of the numerical domain. This means that the potential field solution will not approximate the envelope field everywhere. It will approximate the envelope field up to a height where the solution of the Laplace equation will be strongly influenced by the lower boundary (photosphere). Above that height, the potential solution will be influenced by the upper boundary and will describe the external field. So, we examine the values of the potential field along the height. We expect them not to change drastically in the region of the envelope field. We do find that the potential field solution does not describe the envelope field accurately above certain heights (different heights in the different snapshots). Below those heights, the potential field describes the envelope field well. This transition happens around $z\approx 15\mbox{--}20$ Mm. Therefore, when we calculate the torus index, we do not take into account the values of the torus index when the FR is located above z = 15 Mm.

Zoom In Zoom Out Reset image size

According to the height–time profile, we find that the FR enters the exponential rise phase at $t=79.3\,\mathrm{minutes}$ (first vertical line, Figure 9(b)). The torus index at this time is n = 1.22 and lies in the estimated range of values for the occurrence of a torus instability. During this phase, the maximum $| J/B|$ does not increase dramatically (Figure 9(c)). The current sheet becomes more elongated, and the reconnection outflow becomes more enhanced after t = 81 minutes (Figure 9(d)).

When the EJ-TC reconnection starts, we find that the tension above the FR starts to decrease drastically (second vertical line, Figure 9(e)). Also, after the initiation of the EJ-TC reconnection, the current density of the current sheet becomes more enhanced.

Due to the above, one possible scenario is that the torus instability is responsible for the onset of the exponential phase of the height–time profile, and the EJ-TC reconnection occurs during the rapid rise of the FR. Another possible scenario is that both processes are at work during the eruptive phase, and it is the interplay between them that leads to the fast eruption of the FR.

In terms of the energy, we found that the kinetic energy of the second eruption is larger than that of the first eruption (red line, Figure 3). This difference is not necessarily associated with the different TC reconnection processes. For instance, the downward magnetic tension of the envelope field above the second FR is less. As a result, the upward motion of the FR is faster. Also, the photospheric unsigned magnetic flux increases between the two eruptions due to the continuous emergence. Thus, there is more available flux at the photosphere for the second eruption. Similarly, the magnetic energy in the corona (black line, Figure 3) increases between the two eruptions, indicating that more energy is available for the second eruption.

3.2.3. Third and Fourth Eruptions

After the second FR exits the numerical domain, the overall field-line morphology is similar to the first post-eruption phase. There is an external field, a post-reconnection arcade, that acts as an envelope field and also the J-like loops. At the photosphere–chromosphere, we also find sea-serpent field lines (dark green lines, Figure 10(a)), similar to those in the previous work by Fan (2009) and Archontis et al. (2013). Most of these field lines originate from the partial emergence of the subphotospheric field at different locations along the PIL. These field lines reconnect at many sites along the PIL during the early FR formation. Still, the major role in the FR formation is played by the reconnection of J-like loops (blue and yellow lines, Figure 10(b)).

Zoom In Zoom Out Reset image size

In comparison to the second eruption, we find that the morphology of the external field is different. The second eruption (with a kinetic energy peak at t = 87 minutes) happened right after the first eruption (with a kinetic energy peak at t = 72 minutes). Thus, the external field that the second eruption had to push through was more horizontal (Figure 8(c), gray lines are almost parallel to the photosphere). The third eruption, during which the kinetic energy has its maximum value at t = 119 minutes, happens after the second FR exits the numerical box. As a result, the external field is more vertical to the photosphere and, consequently, it has a very small downward tension (gray lines, Figure 10(b)).

EJ-TC reconnection also occurs during the third eruption (Figure 10(c)). However, we find that only some of the envelope field lines reconnect in both their footpoints (Figures 10(c) and (d)) before they exit the numerical domain. The implication of this difference will be discussed in Section 3.3. We do not find evidence of EE-TC reconnection during the third eruption.

Regarding the torus instability, we should mention that at $t\approx 100\mbox{--}104\,\mathrm{minutes}$, the FR is located very close to the photosphere, at heights $z\approx 1.5\mbox{--}3$ Mm. We find that the value of B_p (and hence n) at these heights depends on the choice of the lower boundary (i.e., the exact height of the photospheric layer, which is used to calculate the potential field). Thus, the values of the torus index for heights up to $z\approx 3$ Mm are different. Above that height, all of the solutions converge. We conjecture that the main reason for the change in the values of B_p and n is the build up of a complex external field after each eruption.

However, from the height–time profile (Figure 11(a)), we find that for $t\simeq 104\,\mathrm{minutes}$, the FR is located just above $z\approx 3$ Mm, where the value of the torus index is well-defined. Also, we find that $n\geqslant 1$ for $t\gt 104\,\mathrm{minutes}$ (first vertical line, Figure 11(b)). This is an indication (although not conclusive) that the torus instability might be associated with the onset of the eruption.

Zoom In Zoom Out Reset image size

Notice that during the time period $t\approx 104\mbox{--}110\,\mathrm{minutes}$, there is no direct evidence that effective reconnection (e.g., EJ-TC reconnection) is responsible for driving the eruption. Figures 11(c), (d), and (e) show that the reconnection upflow underneath the FR undergoes only a small increase (due to the reconnection between J-like field lines), and J/B experiences a limited drop. The tension in the envelope field lines decreases mainly because of the 3D expansion and not because of vigorous EJ-TC reconnection. Therefore, due to the above limitations, we cannot reach a definite conclusion about the exact contribution of reconnection at the onset of the eruption in this initial phase.

In contrast, for $t\gt 110\,\mathrm{minutes}$, there is a clear correlation between the increase of the reconnection outflow and J/B and the decrease of the tension. This is due to effective EJ-TC reconnection, which releases the tension of the envelope field and boosts the acceleration of the erupting field. A preliminary comparison between the second and third eruptions shows that the the maximum values of the current and reconnection outflows are similar, while the length of the CS and the extent of the jet are much smaller. The fourth eruption is very similar to the third eruption.

3.3. Temperature, Density, Velocity, and Current

There are some remarkable similarities and differences among the four eruptions, as illustrated in Figure 12. All panels in this figure are 2D cuts at the vertical xz-midplane, at times just before the erupting structures exit the numerical domain.

Zoom In Zoom Out Reset image size

The density distribution (first row) shows that all eruptions adopt an overall bubble-like configuration, due to the expansion of the magnetic field as it rises to largeratmospheric heights. We notice that the erupting field consists of three main features, which are common in all eruptions. For simplicity, we mark them only in the first eruption (panel a1). These features are (a) the innermost part of the bubble, which is located at and around the center of the erupting field (marked by an asterisk), filled with dense plasma—we refer to this part as the "core" of the eruption, (b) the low-density area that immediately surrounds the "core"—we refer to this as the "cavity" and is the result of the cool adiabatic expansion of the rising magnetic field, and (c) the "front" of the erupting structure, which is a thin layer of dense material that envelops the "cavity" and demarcates the outskirts of the erupting field. To some extent, the shape of the eruptions in our simulations is reminiscent of the "three-part" structure of the observed small-scale prominence eruptions (mini or micro CMEs; e.g., Innes et al. 2010; Raouafi et al. 2010; Hong et al. 2011) and/or CMEs (e.g., Reeves et al. 2015). Because of this, we refer hereafter to the simulated eruptions as CME-like eruptions.

Now, by looking at the temperature distribution (second row), we notice that there is a mixture of cold and hot plasma within the erupting field (in all cases, b1–b4). In fact, in the first eruption, there is a noticeable column of hot plasma, which extends vertically from x = 0 Mm, z = 10 Mm up to z = 40 Mm. Thus, in this case, the "core" of the erupting field appears to be hot, with a temperature of about 8 MK. On the contrary, the "core" of the following eruptions is cool (5000–20,000 K) and dense, but is surrounded by hot (0.5–2 MK) plasma. In all cases, the origin of the hot plasma is the reconnection process occurring at the flare current sheet underneath the erupting field. The distribution of $\sqrt{| J/B| }$ is shown in the fourth row of Figures 12(d1)–(d4). The flare current sheet is the vertical structure with large values of $\sqrt{| J/B| }$ and is located at around x = 0 Mm and between z = 12 Mm and z = 25 Mm. The velocity distribution (panels (c1)–(c4)) shows that a bidirectional flow is emitted from the flare current sheet. This flow is a fast reconnection jet, which transfers the hot plasma upwards (into the erupting field) and downwards (to the flare arcade located below z = 10 Mm).

Thus, a marked difference between the first and the following eruptions is that, in the first eruption, the upward reconnection jet shoots the hot plasma vertically into the "core" of the erupting field, while in the following eruptions, the upward jet only reaches lower heights, arriving below the "core." In the latter cases, the jet is diverted sideways at heights below the center of the erupting FR, adopting a Y-shaped configuration (e.g., see (c2)–(c4)). In the first eruption, the EE-TC reconnection creates field lines that have a highly bent concave-upward shape (i.e., toward the central region of the erupting bubble; see the red lines in Figure 6(d)). It is the strong (upward) tension of these field lines that leads to the hot plasma being ejected at large heights and into the "core" of the field. In the following eruptions, the tension force that accelerates the hot jet upflow is weaker. This is because the reconnected field lines of the jet is the result of reconnection between Js (e.g., see the blue lines in Figure 8(e)), which are not as vertically stretched as the envelope field lines during the EE-TC reconnection. Thus, the upward tension of the reconnected field lines at the flare current sheet is weaker. Therefore, the hot reconnection jet is not strong enough to reach large atmospheric heights and to heat the central region of the erupting field. When it reaches close to the heavy core of the erupting FR, it is diverted sideways (where the pressure is lower), and the embedded hot plasma runs along the reconnected field lines.

In general, the temperature distribution within the overall volume of the erupting field correlates well with the distribution of $\sqrt{| J/B| }$, which implies that heating occurs mainly at sites with strong currents. As we mentioned above, one such area is the flare current sheet underneath the erupting FR. Another example is the heating that occurs at a thin current layer formed between the "core" and the "cavity" of the erupting bubble (e.g., see panel b1). This current layer is formed after the onset of the fast-rise phase of the erupting FR. The uppermost field lines of the erupting core are moving upwards with a higher speed than the field lines within the "cavity," which rise due to the expansion of the emerging field. Thus, at the interface between the two sets of field lines, the plasma is compressed and heated locally (up to 1 MK). This process occurs in all cases (panels (d1)–(d4)), although it is more clearly visible in the first eruption (panel (b1)). The reconfiguration of the field after the first eruption leads to a more complex field-line morphology, distribution of $\sqrt{| J/B| }$, and heating within the rising magnetized volume (panels (b1)–(b4)).

3.4. Geometrical Extrapolation

Coronographic observations of CMEs show that they usually exhibit a constant angular width (i.e., the flanks of the erupting structure move upward, along two approximately straight lines; e.g., Moore et al. 2007). Based on that, we perform a geometrical extrapolation of the size of the first eruption. For this, we find the location of the flanks of the structure at consecutive times and fit a straight line. First, we mark the location of the flank of the erupting structure at a time t_i, when the flank is very distinguishable (diamond on the left flank, Figure 13(a)). Next, we select the flank location prior to t_i (marked with ${t}_{i-6},{t}_{i-5}$, etc.) and after (marked with ${t}_{i+1}$), and fit a straight line through these points (blue line). We then do the same for the other flank. The point where they intersect is approximately the height of the initiation. These extrapolated lines are also plotted in the 3D volume of our numerical box for better visualization (Figure 13(b)).

Zoom In Zoom Out Reset image size

After we find these lines, we extrapolate them to $0.6\,{R}_{\odot }$. For size comparison, we plot them on the solar limb (blue lines, Figure 13(b)). The box at the bottom of the extrapolations shows the size of our numerical box. It is clear that although the eruptions originate from a small-scale region, they grow in size, and it is not unlikely that they will evolve into considerably larger-scale events. We should note that the above method is a first-order approximation regarding the spatial evolution of the first eruption, assuming that the erupting field will continue to rise and expand even after it leaves the numerical domain.

The maximum value of the magnetic energy in the simulated eruptions is 1 × 10²⁸ erg and the kinetic energy varies in the range $3\times {10}^{26}\mbox{--}1.5\times {10}^{27}$ erg. Based on the size of our numerical box and the aforementioned values of energies, the eruptions in this simulation could describe the formation and ejection of small-scale CME-like events. Most CMEs have typical values of kinetic energies around ${10}^{28}\mbox{--}{10}^{30}$ erg (Vourlidas et al. 2010).

In this work, we studied the formation and triggering of recurrent eruptions in an emerging flux region using numerical simulations. The initial emergence of the subphotospheric flux tube formed a bipolar region at the photosphere. The combination of shearing motions and the rotation of the two opposite polarities formed J-like field lines, which reconnected to create an FR that eventually erupted ejectively toward the outer solar atmosphere. In total, four successive eruptions occurred in the simulation. We found that the strength of the magnetic envelope field above the eruptive FRs dropped fast enough that the FRs became torus unstable.

The initial slow-rise phase of the first FR started due to the torus instability. The rising FR pushed the envelope field upwards. The field lines of the envelope field reconnected in a tether-cutting manner and as a result, the tension of the overlying field dropped exponentially. At that time, the FR entered the fast-rise phase. The field lines formed due to the reconnection of Js and turned about once around the axis of the FR, while the field lines resulting from the tether-cutting of the envelope field turned about at least twice around the axis of the FR. The reconnected field lines that were released downwards formed a post-reconnection arcade.

After the eruption of the first FR, reconnection of J-like field lines continued to occur and another FR, which eventually erupted, was formed. This process of FR formation occurred two more times in a similar manner. In all cases, the post-reconnection arcade acted as a new "envelope" field for the next FR. We found that the envelope field was decaying fast enough to favor torus instability. The envelope fields between the second, third, and fourth eruptions differed mostly at the height where the FRs became torus unstable ($n\approx 1\mbox{--}2$). However, we should highlight that our calculation of the torus index is approximate because the envelope field evolves dynamically (e.g., it undergoes expansion). The derivation of the torus instability criteria based on previous analytical studies took into account perturbations of a static configuration. Thus, a more accurate estimate of the torus index in our simulations would be to let the envelope field relax at each time step and then calculate n. This can only be done if the driver of the system could be stopped, letting the overall magnetic flux system reach an equilibrium (e.g., Zuccarello et al. 2015). However, in our dynamical simulations, there is a certain amount of available magnetic flux, which can emerge into the photosphere and above. The driver of the evolution of the system (i.e., magnetic flux emergence) cannot be stopped before the available magnetic flux is exhausted. Therefore, on this basis, we study the continuous evolution of the system. Still, in our experiments, the magnitude of the current inside the envelope field is at least 10 times lower than the one in the FR core, so we expect that the envelope field is not far away from the potential state.

The removal of the downward tension of the envelope field is important for the erupting FRs. In the first eruption, the removal of the envelope tension occurred through the reconnection of the envelope field with other envelope field lines (EE-TC reconnection). In the other three eruptions, the envelope field reconnected with J-like field lines (EJ-TC reconnection). The differences between the EE-TC reconnection and the EJ-TC reconnection were found to be significant for the density and temperature distribution within the erupting structure. After the EE-TC reconnection, the reconnected field lines underneath the erupting FR adopted a W-like shape, with two concave-upward regions (red lines in Figure 6(d); see arrows). After the EJ-TC reconnection, the lower segments of the reconnected field lines adopted a hook-like shape (red lines in Figure 8(f); see arrows).

In the case of EE-TC reconnection, the upward tension of the reconnected field lines (as illustrated by the upward-stretched segments in the middle of the W-shaped field lines) pushed hot plasma from the flare current sheet into the erupting field via a hot and fast collimated jet. Due to this process, the temperature of the central region of the erupting FR changed during the eruption, from low to high values (Figure 12(b1)).

In the case of EJ-TC reconnection, the plasma transfer from the flare current sheet to the erupting field was mainly driven by the reconnection of Js, and therefore the resulting reconnection jet was not as collimated as the one in the EE-TC reconnection. In the second eruption, this post-reconnection hot jet collided with the FR and became diverted into two side jets (Figures 12(a2) and (c2)). In the third and fourth eruptions, the jets were not fast enough to enter the region of the erupting core of the field (Figures 12(a3) and (c3), and (a4) and (d4)).

Thus, the study of the temperature distribution revealed that due to EE-TC reconnection, the erupting field develops a "three-part" structure consisting of a hot front "edge," a cold "cavity," and a hot and dense "core." In the following eruptions, the temperature of the plasma within the central region of the FRs remained low. Therefore, we suggest that observations of erupting FRs, which are heated, e.g., from 10³ to 10⁶ K, during their eruptive phase, might indicate that EE-TC reconnection is at work. We should mention that heat conduction is not included in our simulation. Therefore, the exact value of the temperature within the erupting field may change if heat conduction were to be included in the numerical experiment.

Overall, we report that the physical mechanism behind the formation of recurrent ejective eruptions in our flux emergence simulation is a combination of torus-unstable FRs and the onset of tether cutting of the overlying field through a flare current sheet. Both the EE-TC reconnection and the EJ-TC reconnection were found to remove the downward tension of the overlying field and thus assist the eruptions. In the first eruption, it is likely that torus instability occurs first, and the rapid exponential rise phase of the erupting FR comes after the EE-TC reconnection. For the other eruptions, where the structure of the magnetic field above the FR has a more intricate morphology, it is difficult to conclude which process is responsible for the onset of the various phases of the eruptions.

Comparing our results with previous studies, the formation of all FRs in our simulation is due to the reconnection of sheared J-like field lines, in a similar manner to earlier simulations (e.g., Aulanier et al. 2010; Archontis & Hood 2012; Leake et al. 2013, 2014).

It is also interesting to note that the velocity and current profiles of our first eruption (Figures 12(c1), (d1)) are very similar morphologically to the ones produced from the flare reconnection in the breakout simulation of Karpen et al. (2012), who used a (different) 2.5D adaptive grid code. Such similarities indicate that the resulting morphologies might be generic and indicative of the EE-TC reconnection.

Moreno-Insertis & Galsgaard (2013) performed a flux emergence simulation of a highly twisted flux tube into a magnetized atmosphere and found recurrent eruptions. In comparison to our simulation, the subphotospheric flux tube in the work of Moreno-Insertis & Galsgaard (2013) had a higher magnetic field strength (${B}_{0}=3.8$ kG) and greater length of the buoyant part of the flux tube (λ = 20 in comparison to our λ = 5), and was located closer to the photosphere ($z=-1.7$ Mm). In their work, their first FR is formed, similarly to our simulation, by the reconnection of sheared arcade field lines. The higher λ leads to the formation of a more elongated emerging FR and a longer sigmoid. The eruption mechanism, though, is very different. It involves reconnection between the sheared arcade field lines and the open field lines of the ambient field. Also, it involves the reconnection of the sheared arcade with a magnetic system produced from the reconnection of the ambient field with the initial emerging envelope field. Their second and third eruptions are off-centered eruptions of segments of the initial flux tube, which eventually become confined by the overlying field. In our case, the flux tube axis emerges only up to 2–3 pressure scale heights above the photosphere (z = 0), and the erupting FRs are all formed due to the reconnection of J-loops.

Murphy et al. (2011) discussed possible heating mechanisms for the dynamic heating of CMEs, one of which is heating from the CME flare current sheet. Taking into account the results of previous studies (e.g., Lin et al. 2004), they reported that the reconnection of hot upward jets from the flare current sheet could reach the cool central region of the erupting FR and heat it. In fact, this leads to some mixing of hot and cool plasma within the central erupting volume.

From the two different tether-cutting reconnections found in our simulation, only the EE-TC reconnection allows the effective transfer of hot plasma from the flare current sheet into the FR central region through the reconnection outflow. This might account for a process similar to the aforementioned mixing of hot and cold plasma, as suggested by Lin et al. (2004). On the other hand, during EJ-TC reconnection, hot plasma is mainly found at the periphery of the central region of the FR.

The physical size of our simulated emerging flux region was 23.4 Mm, and the size of the FRs was up to 64.8 Mm (the length of the y-axis). The height of our numerical box was 57.6 Mm. The kinetic energies of the eruptions were $3\times {10}^{26}\,\mbox{--}\,1.5\times {10}^{27}$ erg and the magnetic energies around 1 × 10²⁸ erg. These values suggest that our numerical experiment describes an emerging flux region, which hosts relatively low-energy eruptions in comparison to CMEs. Based on the sizes and the energetics, these eruptions can describe the formation and eruption of small-scale eruptive events. For instance, such an eruption in terms of physical size, and not magnetic configuration, was reported by Raouafi et al. (2010) and Reeves et al. (2015). Still, the results on the plasma transfer for the different flare reconnections (EE-TC reconnection and EJ-TC reconnection) should be scale invariant.

Having reproduced a CME-like configuration (Figures 6(a1) and (b1)), we extrapolated the expansion of the flanks of the erupting "bubble" and estimated its size in $0.6\,{R}_{\odot }$. We found that these eruptions have the potential to become comparable to small-sized CMEs (Figure 13(c)), but with one order of magnitude lower kinetic energy. We aim to study the parameters that would increase the energies of the produced eruptions. For this, in our next paper, we will present the results of a parametric study of the magnetic field strength of the subphotospheric flux tube. Our aim is to study the differences in the energetics, physical size, and recurrence of the eruptions.

The authors would like to thank the referee for the constructive comments. This project has received funding from the Science and Technology Facilities Council (UK) through the consolidated grant ST/N000609/1. This research has been co-financed by the European Union (European Social Fund—ESF) and Greek national funds through the Operational Program "Education and Lifelong Learning" of the National Strategic Reference Framework (NSRF)—Research Funding Program: Thales Investing in Knowledge Society through the European Social Fund. The authors acknowledge support by the Royal Society. This work was supported by computational time granted from the Greek Research & Technology Network (GRNET) in the National HPC facility—ARIS. This work used the DIRAC 1, UKMHD Consortium machine at the University of St Andrews and the DiRAC Data Centric system at Durham University, operated by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (http://www.dirac.ac.uk). This equipment was funded by BIS National E-infrastructure capital grant ST/K00042X/1, STFC capital grant ST/H008519/1, and STFC DiRAC Operations grant ST/K003267/1, and Durham University. DiRAC is part of the National E-Infrastructure.