# A fast method to identify mean motion resonances

E. Forgács-Dajka1,2, Zs. Sándor1,3, and B. Érdi1
1Department of Astronomy, Eötvös Loránd University, Pázmány Péter sétány 1/A,
H-1117 Budapest, Hungary
2Wigner RCP of the Hungarian Academy of Sciences, 29-33 Konkoly-Thege Miklós
Str, H-1121 Budapest, Hungary
3Konkoly Observatory, Research Centre for Astronomy and Earth Sciences,
Hungarian Academy of Sciences,
H-1121 Budapest, Konkoly Thege Miklós út 15-17., Hungary E-mail: E.Forgacs-
<EMAIL_ADDRESS>
(Accepted 2018 March 7. Received 2018 March 7; in original form 2018 January
3)

###### Abstract

The identification of mean motion resonances in exoplanetary systems or in the
Solar System might be cumbersome when several planets and large number of
smaller bodies are to be considered. Based on the geometrical meaning of the
resonance variable, an efficient method is introduced and described here, by
which mean motion resonances can be easily find without any a priori knowledge
on them. The efficiency of this method is clearly demonstrated by using known
exoplanets engaged in mean motion resonances, and also some members of
different families of asteroids and Kuiper-belt objects being in mean motion
resonances with Jupiter and Neptune respectively.

###### keywords:

celestial mechanics – methods: numerical – (stars:) planetary systems

††pubyear: 2017††pagerange: A fast method to identify mean motion resonances–A
fast method to identify mean motion resonances

##  1 Introduction

As a result of the ongoing discoveries, the number of exoplanets is
continuously growing; up to now more than 3500 planets are known in more than
2600 planetary systems. These numbers indicate that a significant amount of
planets form multi-planet systems, in which at least two planets are revolving
around a central star. Planets form in gas rich protoplanetary disks,
therefore beside gravity, forces arising from the ambient gaseous material
play important role as altering typically their semi-major axes during a
process called orbital migration. If the orbital distance between two planets
decreases in time due to their migration, the planets can be orbitally locked
at certain ratio of their semi-major axes, where the ratio of the mean angular
velocities (mean motions) can be expressed as a ratio of small positive
integer numbers. If certain dynamical conditions are fulfilled, a resonant
capture may happen meaning that the ratio of the semi-major axes and that of
the mean motions does not change during the further migration of the planets.

Before the Kepler-era, mean motion resonances (MMRs) outside the Solar System
have been detected between giant planets. One of the most remarkable resonant
systems discovered at that time is around the M dwarf star GJ 876 in which two
planets are in 2:1 MMR (Marcy et al., 2001). Later on three planets in
Laplace-type resonance have been identified (Rivera et al., 2010), being the
first known resonant chain of exoplanets resembling to the Jupiter’s Galilean
satellites Io-Europa-Ganymede. Other examples for the 2:1 MMR are the system
HD 128311, in which two massive planets are in resonance but with only one
critical argument librating (Vogt et al., 2005), and HD 73526, also with one
librating critical argument (Tinney et al., 2006). In the system HD 60532 two
giant planets are in 3:1 MMR (Desort et al., 2008). A further example for a
planetary system that might be in the Laplace resonance is the system HR 8799
in which three massive planets are orbiting around an A-type star (Goździewski
& Migaszewski, 2014). More recently, more planets have been discovered to be
dynamically locked in chains of MMRs among the Kepler systems. In Kepler-60
three planets with masses around $\sim 4M_{\oplus}$ are in a 5:4:3 Laplace-
type MMR (Goździewski et al., 2016). The four Neptune-mass planets in
Kepler-223 have periods in ratios close to 3:4:6:8 (Mills et al., 2016) being
another clear indication of planetary migration.

Beside the exoplanetary systems, mean motion resonances also play an important
role in shaping the dynamics of the Solar System bodies. Mean motion
resonances in the Solar System usually occur between a planet and small
bodies: e.g. the members of the Hilda group of asteroids are in a 3:2, while
the Trojan asteroids are in a 1:1 MMR with Jupiter. MMRs between terrestrial
planets and particular members of the asteroid family Hungaria can also be
found. Moreover, as the existence of the Kirkwood gaps clearly indicates, the
dynamical structure of the main belt is shaped by the MMRs between the
asteroids and Jupiter (though their depletion is not only due to MMRs).
Regarding larger objects, between Neptune and the dwarf planet Pluto there is
a 3:2 MMR, which is a protective resonance keeping them on non-near-
approaching orbital positions. This is of high importance, since due to the
high eccentricity of Pluto’s orbit ($e=0.25$) the projection of its perihelion
to Neptune’s orbital plane lies inside the orbit of Neptune thus the two
orbits are crossing each other (in projection). Interestingly, there are many
other Kuiper-belt objects that are locked in the 3:2 MMR with Neptune, these
are the plutinos sharing their orbits with Pluto.

When studying the dynamics of exoplanetary systems with planets engaged in
chains of mean motion resonances, or the motion of small celestial bodies,
asteroids or trans-Neptunian objects (TNOs), one often finds the problem how
to identify the possible mean motion resonances among the bodies involved. For
example, during the formation of the dynamical structure of the trans-
Neptunian region, resonances swept through the primordial disk due to the
migration of the planets, and planetesimals were captured into mean motion
resonances. Modelling this process, a monitoring is necessary to see how the
captures take place, which planetesimals are captured to which resonances
(Hahn & Malhotra, 2005; Levison et al., 2008). Similarly, studying the
dynamical evolution of an asteroid family, one has to check possible
resonances with the planets (Galiazzo et al., 2013). This needs the
computation of the resonance variables, to check their possible libration, for
all reasonable (often very high order) resonant ratios of the mean motions of
all small bodies with all considered planets, and this has to be repeated
several times during the investigated time-span resulting in a considerable
computational workload.

In order to make easier the identification of the different mean motion
resonances, we present a method we name _FAIR_ (as _FA_ st _I_ dentification
of mean motion _R_ esonances) that is easy to use and by which the
identification of the MMRs is possible without any a priori knowledge on them.
Our paper is organized as the following: in the next section we introduce
_FAIR_ for inner and outer type MMRs, then in two further sections particular
examples are presented for Solar System objects, and for known exoplanetary
systems, also including the study of resonant chains of planets.

##  2 The method _FAIR_

Let us denote by $a$, $a^{\prime}$ the semi-major axes of two celestial bodies
locked in a MMR, and by $n$, $n^{\prime}$ the corresponding mean motions. For
the sake of simplicity let us deal with the situation of an asteroid and a
planet, in this case the primed quantities refer to the planet. This notation
enables us to distinguish between inner and outer resonances from the point of
view of the asteroid. In the case of two giant planets neither of the bodies
are in such distinguished positions, the non-primed quantities refer to the
planet being the body to which respect the MMR is studied.

Considering an asteroid and a planet, they can be either in

(i) _inner_ mean motion resonance, if $a<a^{\prime}$ and the ratio of the mean
motions is approximately

$\frac{n}{n^{\prime}}=\frac{p+q}{p},$ (1)

or in

(ii) _outer_ mean motion resonance, if $a>a^{\prime}$ and approximately

$\frac{n^{\prime}}{n}=\frac{p+q}{p},$ (2)

where $p$ and $q$ are relative prime integers. Here $q$ is the order of the
resonance. The resonant perturbations are proportional to the $q$th power of
the orbital eccentricities, thus for small eccentricities low order resonances
are the most important. However, for large eccentricities high order
resonances can also be significant. (These refer to the so-called
eccentricity-type resonance, when the longitudes of the perihelion are also
involved, as in the cases studied below.)

Resonances can be studied by using resonance variables. For a given mean
motion resonance there can be several types of resonance variables (Murray &
Dermott, 1999). From among them we consider those for which the method _FAIR_
is applicable.

###  2.1 Inner resonance

Since beside the selected asteroid and planet there are also other planets in
a planetary system, we should also take into consideration the secular rates
$\dot{\varpi}$, $\dot{\varpi}^{\prime}$ of the perihelion longitudes, caused
by the mutual gravitational perturbations.

In this case one can obtain the following resonance variables:

$\theta_{1}\equiv(p+q)\lambda^{\prime}-p\lambda-q\varpi,$ (3)
$\theta_{2}\equiv(p+q)\lambda^{\prime}-p\lambda-q\varpi^{\prime},$ (4)

where $\lambda=M+\varpi$, $\varpi=\omega+\Omega$, $\lambda$ is the mean
orbital longitude, $M$ the mean anomaly, $\omega$ the argument of perihelion,
$\Omega$ the longitude of the ascending node, $\varpi$ the longitude of the
perihelion of the asteroid. The primed variables refer to the planet with
similar meaning.

Eqs (3) and (4) refer to the case of exact resonance. Near a resonance, Eqs
(3) and (4) are satisfied approximately, and the resonance variables oscillate
(librate) around a mean value, which is usually $0^{\circ}$ or $180^{\circ}$,
with amplitudes $\Delta\theta<180^{\circ}$.

Actually, in the neighbourhood of a resonance, the behaviour of a resonance
variable can be quite complex. Depending on the actual system and initial
conditions, it can exhibit libration, circulation, alternation between
libration and circulation, or chaotic fluctuations. When searching for
resonance, our aim is to establish that the resonant ratio of the mean motions
(or their combinations with the secular frequencies) is maintained and the
corresponding resonance variable librates. These are necessary and sufficient
conditions that a resonance exist. We note that the above introduced resonance
variables appear as critical arguments in the series expansion of the
perturbed two-body potential in terms that contain the $q\mathrm{th}$ powers
of $e$ and $e^{\prime}$. There are mixed types of eccentricity-type resonances
too, corresponding to critical arguments (and resonance variables) that
contains coefficients of $ee^{\prime}$ type. In this work, however, we do not
consider the latter cases of mixed eccentricity-type resonances.

Now let us consider Eq. (3) in two cases. First, when
$\lambda=\lambda^{\prime}$, that is at the conjunction of the asteroid and the
planet, it follows from (3) that

$\theta_{c}=q(\lambda-\varpi)=qM,$ (5)

where $\theta_{c}$ is the value of $\theta_{1}$ at conjunction. Writing
$\bar{\theta}=\theta_{10}$ and denoting the amplitude of $\theta_{1}$ by
$\Delta\theta$, since

$\bar{\theta}-\Delta\theta\leq\theta_{c}\leq\bar{\theta}+\Delta\theta,$

it follows that

$\bar{\theta}-\Delta\theta\leq qM\leq\bar{\theta}+\Delta\theta.$

Here $M$ is the mean anomaly of the asteroid at the moment of conjunctions,
and $qM$ should be taken mod $2\pi$. Thus

$\frac{\bar{\theta}+k2\pi-\Delta\theta}{q}\leq
M\leq\frac{\bar{\theta}+k2\pi+\Delta\theta}{q},$

where $k$ is an integer. This gives $q$ centres

$\frac{\bar{\theta}}{q}+k\frac{2\pi}{q},\quad k=0,1,\ldots,q-1$

along the orbit of the asteroid around which the conjunctions of the asteroid
with the planet can take place within regions of half-size $\Delta\theta/q$.
For example, for $q=3$ and $\bar{\theta}=0^{\circ}$ the centres are at
$0^{\circ}$, $120^{\circ}$, and $240^{\circ}$, while for
$\bar{\theta}=180^{\circ}$ these are at $60^{\circ}$, $180^{\circ}$, and
$300^{\circ}$.

In the second case, when $M=0^{\circ}$, that is at the perihelion of the
asteroid, $\varpi=\lambda$, and it follows from (3) that

$\theta_{p}=(p+q)(\lambda^{\prime}-\lambda),$

where $\theta_{p}$ is the value of $\theta$ at the perihelion. Since

$\bar{\theta}-\Delta\theta\leq\theta_{p}\leq\bar{\theta}+\Delta\theta,$

it follows that

$\bar{\theta}-\Delta\theta\leq(p+q)(\lambda^{\prime}-\lambda)\leq\bar{\theta}+\Delta\theta.$

As before, it also follows that

$\frac{\bar{\theta}+k2\pi-\Delta\theta}{p+q}\leq\lambda^{\prime}-\lambda\leq\frac{\bar{\theta}+k2\pi+\Delta\theta}{p+q}$

giving $p+q$ centres

$\frac{\bar{\theta}}{p+q}+k\frac{2\pi}{p+q},\quad k=0,1,\ldots,p+q-1$

around which in regions with half-size $\Delta\theta/(p+q)$ the asteroid can
be found with respect to the planet when the former is at perihelion.

The previous considerations can be summed up in stating that plotting
$\lambda^{\prime}-\lambda$ against $M$ in a rectangular coordinate system with
$M$ as horizontal and $\lambda^{\prime}-\lambda$ as vertical axis, there will
be $q$ centres on the horizontal, and $p+q$ centres on the vertical axis.

The same reasoning can be repeated with Eq. (4), concluding to that plotting
$\lambda^{\prime}-\lambda$ against $M^{\prime}$ there will be $q$ centres on
the horizontal, and $p$ centres on the vertical axis.

###  2.2 Outer resonance

In the case of an outer resonance, the resonance variables for the
eccentricity-type resonances are

$\theta_{3}=(p+q)\lambda-p\lambda^{\prime}-q\varpi,$ (6)

and

$\theta_{4}=(p+q)\lambda-p\lambda^{\prime}-q\varpi^{\prime}.$ (7)

By similar considerations as before, it follows from Eq. (6) that plotting
$\lambda-\lambda^{\prime}$ against $M$, there will be $q$ centres on the
horizontal, and $p$ centres on the vertical axis. In the case of Eq. (7),
plotting $\lambda-\lambda^{\prime}$ against $M^{\prime}$, there will be $q$
centres on the horizontal, and $p+q$ centres on the vertical axis. Table 1
gives a summary of the different cases in the inner and outer resonances.

###  2.3 Inclination-type resonances of inner and outer type

We have investigated the eccentricity-type resonances so far, however, it can
be easily shown that with minor modifications, the method _FAIR_ also works
for inclination-type resonanances. Let us consider the non-mixed inclination-
type resonances of order $q$, where the critical arguments appear in terms
containing the $q$th power of either $I$ or $I^{\prime}$ (being the
inclinations of the body and of the perturber, respectively) as coefficients
in the series expansion of the perturbing potential:

$\theta_{I,1}=(p+q)\lambda^{\prime}-p\lambda-q\Omega,$ (8)
$\theta_{I,2}=(p+q)\lambda^{\prime}-p\lambda-q\Omega^{\prime},$ (9)

where $\Omega$ and $\Omega^{\prime}$ are the longitude of nodes of the body
and the perturber, respectively. Equating the mean longitudes,
$\lambda=\lambda^{\prime}$, similarly to Equation (5), conjunctions happen
when

$\theta_{I,c}=q(\lambda-\Omega)=q(\omega+M).$ (10)

This suggests that in the inclination-type resonances one should replace $M$
with $M+\omega$, and by a similar argumentation to the eccentricity-type
resonances, there will be $q$ centres on the horizontal, and $(p+q)$ centres
on the vertical axis when plotting $\lambda^{\prime}-\lambda$ versus
$M+\omega$ in a rectangular coordinate system. The above considerations can
also be applied to the outer inclination-type resonances, when
$\lambda^{\prime}-\lambda$ should be plotted either versus $M+\omega$ or
$M^{\prime}+\omega^{\prime}$. Thus by replacing in Table 1 $M$ and
$M^{\prime}$ with $M+\omega$ and $M^{\prime}+\omega^{\prime}$, moreover
$\varpi$ and $\varpi^{\prime}$ with $\Omega$ and $\Omega^{\prime}$
respectively, one can easily obtain the resonance variable of the
corresponding inclination-type resonance.

Table 1: Resonance variables. The third column shows the variables to be plotted versus each other, the fourth and fifth the number of centres on the horizontal and vertical axes. type | resonance variable | plot | hor | vert
---|---|---|---|---
inner | $(p+q)\lambda^{\prime}-p\lambda-q\varpi$ |  $\lambda^{\prime}-\lambda$ vers $M$ | $q$ | $p+q$
inner | $(p+q)\lambda^{\prime}-p\lambda-q\varpi^{\prime}$ |  $\lambda^{\prime}-\lambda$ vers $M^{\prime}$ | $q$ | $p$
outer | $(p+q)\lambda-p\lambda^{\prime}-q\varpi$ |  $\lambda-\lambda^{\prime}$ vers $M$ | $q$ | $p$
outer | $(p+q)\lambda-p\lambda^{\prime}-q\varpi^{\prime}$ |  $\lambda-\lambda^{\prime}$ vers $M^{\prime}$ | $q$ | $p+q$

##  3 Identification of resonances

The properties of the resonance variables, described in Section 2, can be used
to find easily whether a resonance exists between a planet and an asteroid.
Deciding from the semi-major axes that the looked for resonance is inner or
outer, then searching for inner resonances, one plots
$\lambda^{\prime}-\lambda$ versus $M$ and $M^{\prime}$, (or in the cases of
inclination-type MMRs versus $M+\omega$ and $M^{\prime}+\omega^{\prime}$)
computed by numerical integration of the equations of motion for sufficiently
long time. If there is a resonance, some stripes will appear on the plots (see
the examples below), and counting the number of their intersections with the
horizontal and vertical axis will provide the $p$ and $q$ values of the
resonance. Searching for outer resonances, the procedure is the same with
plotting $\lambda-\lambda^{\prime}$ versus $M$ and $M^{\prime}$ and counting
the numbers of intersections on the horizontal and vertical axis. Since in the
case of a resonance the appearance of the proper number of intersections is a
necessary condition, once $p$ and $q$ are found, the ratio of the mean motions
and the libration of the resonance variable should be checked with them.

When processing the above plots, there can be raised two technical questions:
(i) how long the numerical integrations should be carried out, and (ii) how
frequently should the points be plotted in order to obtain good quality plots
suitable for determination of a MMR. Without giving a definite answer, we have
empirically found that regarding the plotting frequency, the one hundredth of
the period of the inner body certainly will give good result when using an
integration time that corresponds to one hundred periods of the outer body.

###  3.1 Solar System examples

In this section we provide a few examples for the application of the method
_FAIR_ to Solar System bodies involved in various mean motion resonances. The
orbits have been numerically calculated by using an own-developed adaptive
step size Runge-Kutta-Nystrom 6/7 N-body integrator (Dormand & Prince, 1978).
In our integrations all planets except Mercury have been taken into account,
and as initial conditions we used heliocentric coordinates and velocities
taken from JPL Horizons at the epoch JD 2457754.50.

####  3.1.1 153 Hilda - Jupiter 3:2 MMR

Figure 1: Asteroid 153 Hilda in a 3:2 MMR with Jupiter.

The main dynamical characteristic of the Hilda group of asteroids (named after
153 Hilda) is that its members are locked in a 3:2 MMR with Jupiter (Brož &
Vokrouhlický, 2008). As a first example on how the method _FAIR_ works is
given by plotting $\lambda^{\prime}-\lambda$ versus $M$ (see Figure 1, left
panel), where $\lambda^{\prime}$ is the mean longitude of the perturbing body
(Jupiter, in this case), and $\lambda$, $M$ are the mean longitude, and mean
anomaly of the asteriod 153 Hilda, respectively. Since this is a MMR of inner
type, and $\lambda^{\prime}-\lambda$ is plotted against $M$, according to the
first row of Table 1 the counting of the number of intersecting stripes with
the horizontal and vertical axis gives $q=1$ and $p+q=3$, see the left panel
of Figure 1. Thus, the ratio of the mean motions is $(p+q)/p=3/2$, and the
corresponding resonance variable (or resonant angle) is
$\theta_{1}=3\lambda^{\prime}-2\lambda-\varpi$, librating in this case around
$0^{\circ}$ as shown in the right panel of Figure 1. We note that in this
case, the ratio of the mean motions and the resonance variable are determined
solely by counting the crossings of the stripes with the horizontal and
vertical axes.

####  3.1.2 Neptune - Pluto 3:2 MMR

Figure 2: The 3:2 MMR between Neptune and Pluto

As a further well known example, we consider the case of Neptune and Pluto,
being in a 3:2 MMR (Cohen & Hubbard, 1965). In this case, the more massive
Neptune (the perturber) orbits closer to the Sun, thus the MMR under study is
of outer type. In the left panel of Figure 2 we displayed
$\lambda-\lambda^{\prime}$ versus $M$. Since this is a MMR of outer type,
according to the third row of Table 1 the counting of the crossings of the
stripes with the horizontal and vertical axis gives $q=1$ and $p=2$. Thus the
ratio of the mean motions is $(p+q)/p=3/2$, and the corresponding resonance
variable, $\theta_{3}=3\lambda-2\lambda^{\prime}-\varpi$ librates around
$180^{\circ}$.

We note that although the above two well known mean motion resonances can be
identified easily without the application of the method _FAIR_ , we were able
to identify them without knowing their character. In what follows, we show
some more sophisticated cases, in which the identification of the resonances
and the corresponding resonance variables might be more complicated.

####  3.1.3 Neptune - 2005 TN74 5:3 MMR

Figure 3: The 2005TN74 trans-Neptunian object in a 5:3 MMR with Neptune.

2005 TN74 is a trans-Neptunian object (TNO) discovered by Sheppard et al.
(2005). In order to demonstrate the effectiveness of the method _FAIR_ , let
us forget for now that this body is in a 5:3 MMR with Neptune. The only
information we suppose to have is that this TNO orbits outside Neptune,
therefore when seeking for a resonant behaviour we should consider a MMR of
outer type. To do so we display the plot $\lambda-\lambda^{\prime}$ versus
$M$, where $\lambda^{\prime}$ is the mean longitude of Neptune, while
$\lambda$ and $M$ are the mean longitude and mean anomaly of 2005 TN74,
respectively. Studying the left panel of Figure 3, we can count the numbers of
crossings of the stripes with the axes, yielding $q=2$ and $p=3$ (see the
third row of Table 1). Thus the ratio of the mean motions is $(p+q)/p=5/3$,
and the resonance variable $\theta_{3}=5\lambda-3\lambda^{\prime}-2\varpi$
librates around $180^{\circ}$, see the right panel of Figure 3.

####  3.1.4 Neptune - 136108 Haumea 12:7 MMR

Figure 4: Dwarf planet 136108 Haumea in a 12:7 MMR with Neptune

136108 Haumea is a dwarf planet orbiting beyond Neptune, being discovered
independently by Brown et al. (2005) and Ortiz et al. (2005). It is known that
there is a 12:7 MMR of outer type between Haumea and Neptune. This is a 5th
order MMR, and its identification might not be too straightforward being the
ratio of the mean motions close to the 2:1 ratio. Seeking a resonant behaviour
with Neptune of outer type, we plot $\lambda-\lambda^{\prime}$ against $M$,
see Figure 4. As usual, $\lambda^{\prime}$ refers to Neptune’s mean orbital
longitude, while the non-primed quantities to Haumea’s orbital elements.
Investigating Figure 4, one can identify $q=5$ crossings with the horizontal,
and $p=7$ crossings with the vertical axis. This yields to a mean motion ratio
of $(p+q)/p=12/7$, and a resonance variable
$\theta_{3}=12\lambda-7\lambda^{\prime}-5\varpi$ librating around
approximately $120^{\circ}$.

####  3.1.5 Neptune - 2001 QR322 1:1 MMR

Figure 5: Neptune and 2001QR322 in 1:1 MMR

The method _FAIR_ can also be applied to Trojan-type motions. The main
characteristic of the Trojan-type, or in other words co-orbital motion is the
1:1 MMR meaning that the bodies involved are sharing similar orbits. The best
known examples for Trojan-type motion are the Trojan asteroids populating the
neighbourhood of the stable triangular Lagrangian points $L_{4}$ and $L_{5}$
of the Sun-Jupiter system. The resonance variable in these cases is
$\lambda-\lambda^{\prime}$ which librates either around $60^{\circ}$, or
$300^{\circ}$ depending on whether the asteroid is in the vicinity of the
$L_{4}$, or $L_{5}$ point. In the $\lambda^{\prime}-\lambda$ versus $M$ plots
of the Neptune’s Trojan 2001 QR (see Figure 5, left panel) we find one strip
parallel to the horizontal axis of the coordinate system yielding $q=0$ (see
Table 1, first row), e.g. a zeroth order resonance. The only crossing with the
vertical axis means $p+q=1$, that is $p=1$. Thus the resonance is
$(p+q)/p=1/1$, and the resonance variable is
$\theta=\lambda^{\prime}-\lambda$, librating around $300^{\circ}$ (Figure 5,
right panel).

###  3.2 Mean motion resonances in exoplanetary systems

Finally, with two additional examples we demonstrate the applicability of the
method _FAIR_ to identify resonances between pairs of exoplanets, or chains of
resonances in which many of them are involved. The numerical integration of
orbits have been performed by using the already mentioned RKN 6/7 integrator.

####  3.2.1 HD 60532

Figure 6: The HD 60532 system in 3:1 inner MMR

Until the discovery of the system around HD 60532, an F-type star with mass
$M_{*}=1.44~{}M_{\odot}$ by Desort et al. (2008), there were not known any
giant planets involved in the 3:1 MMR. The final confirmation that giant
planets can be in the 3:1 MMR was given by Laskar & Correia (2009), while a
formation study favouring a migration based scenario was performed by Sándor &
Kley (2010). In the latter study the authors demonstrated that the convergent
migration, thus capture in the 3:1 MMR was only possible for larger planetary
masses corresponding to Fit II of Laskar & Correia (2009). In this work we
checked both Fit I with low and Fit II with high planetary masses to show the
efficiency of the method _FAIR_. Our results are displayed in Figure 6, where
in the upper left panel the $(\lambda^{\prime}-\lambda)$ versus $M$ plot is
displayed for Fit I, while in the bottom left panel for Fit II. These
resonances are of inner types (as the primed orbital elements are refering to
the outer giant planet), and by calculating the number of the crossings of the
strips with the horizontal and vertical axis we have (see Table 1, first row)
$q=2$ and $p+q=3$. Thus $p=1$, and the ratio of the mean motions is
$(p+q)/p=3/1$. This means that without any a priori knowledge, just by a
careful analysis of the $(\lambda^{\prime}-\lambda)$ versus $M$ plot we could
identify the mean motion resonance, and also write down the corresponding
resonant variable, which in these cases is
$\theta_{1}=3\lambda^{\prime}-\lambda-2\varpi$, librating around $180^{\circ}$
(Figure 6, right panels). It is noteworthy that due to the relatively large
libration amplitude, the strips in the $(\lambda^{\prime}-\lambda)$ versus $M$
plot are quite broad.

####  3.2.2 Kepler 60

Figure 7: The resonant pairs of the system Kepler 60, the three planets being
in 5:4:3 Laplace-type MMR

As we have already mentioned, there are planetary systems, typically
discovered by the Kepler mission, in which the planets are captured in chains
of mean motion resonances. A prominent example of these systems is Kepler 60,
where three planets with masses around $4~{}M_{\oplus}$ are in the 5:4:3
Laplace-type MMR (Goździewski et al., 2016). Here we apply the method _FAIR_
to the Fit II of the cited paper, in which the pairs of the planets in MMRs
are investigated. First, without any a priori knowlegde on the MMRs, we can
study planets b and c. We display $\lambda_{1}-\lambda_{2}$ as the function of
$M_{1}$, see the upper left panel of Figure 7. We note that the index 1 stays
for the inner, while index 2 for the outer planet, planets b and c in this
configuration. Considering that $\lambda$, $\varpi$, $M$ correspond to
$\lambda_{1}$, $\varpi_{1}$, $M_{1}$ and $\lambda^{\prime}$,
$\varpi^{\prime}$, $M^{\prime}$ correspond to $\lambda_{2}$, $\varpi_{2}$,
$M_{2}$, and that this is an outer type MMR, according to Table 1 (third row)
the numbers of crossings with the horizontal and vertical axis give $q=1$ and
$p=4$. The mean motions ratio in this case is thus $(p+q)/p=5/4$ corresponding
to a resonant variable $\theta_{3,12}=5\lambda_{1}-4\lambda_{2}-\varpi_{1}$,
librating around $180^{\circ}$ (see the upper right panel of Figure 7).

As a next step, we can also check, whether the planets c and d are in a MMR
too. To do so, we plot $\lambda_{2}-\lambda_{3}$ as the function of $M_{2}$,
shown in the middle left panel of Figure 7. Now the non-primed elements
correspond to those with index 2, and the primed ones to those with index 3.
According to Table 1 (third row), we have $q=1$ and $p=3$, that is there is a
$(p+q)/p=4/3$ MMR with a resonant variable
$\theta_{3,23}=4\lambda_{2}-3\lambda_{3}-\varpi_{2}$, librating around
$180^{\circ}$ (middle right panel of Figure 7).

Finally, we also plot $\lambda_{1}-\lambda_{3}$ as the function of $M_{1}$
(bottom left panel of Figure 7). Counting the numbers of crossings of the
stripes with the axes, we find $q=2$ and $p=3$ corresponding to a mean motion
ratio $(p+q)/p=5/3$, and indeed, the resonant variable
$\theta_{3,13}=5\lambda_{1}-3\lambda_{3}-2\varpi_{1}$ librates around
$~{}225^{\circ}$ (bottom right panel of Figure 7). Thus using the system
Kepler 60, we gave an evidence that the method _FAIR_ is also applicable for a
quick identification of MMRs between pairs of exoplanets being in resonant
chains. The above analysis can also be performed for _inner_ type resonances.

###  3.3 Application of the method _FAIR_ to inclination-type MMR

Similarly to the eccentricity-type MMRs, the method _FAIR_ can also be applied
to detect inclination-type resonances. The best known example for inclination-
type MMR in the Solar System is between the Saturnian satellites Mimas and
Thetis corresponding to a 2:1 commensurability. Inclination-type MMRs can also
be developed between migrating giant planets if the eccentricity damping of
the outer migrating planet is modest, allowing the rapid growth of the inner
planet’s eccentricity. This can additionally lead to the fast growth of the
inner planet’s inclination (see Libert & Tsiganis, 2009). Here we consider the
following example that leads to a 3:1 MMR: the inner planet with mass
$m_{1}=1M_{J}$ (where $M_{J}$ is the mass of Jupiter) is started from 5 au
from nearly circular and planar, initially non-migrating orbit. An outer
planet with mass $m_{2}=2M_{J}$ is started from 16.5 au from a circular orbit
with negligible inclination with respect to the plane of reference of the
coordinate system, and is forced to migrate in timescale of $\tau_{a}=7\times
10^{5}$ years, with the same eccentricity damping timescale. According to
Libert & Tsiganis (2009), the two planets enter first into a 3:1 eccentricity-
type, and later on into an inclination-type MMR, as the libration of the
following eccentricity-type and inclination-type variables clearly indicates:

$\theta_{1}=3\lambda_{2}-\lambda_{1}-2\varpi_{1},$ (11)
$\theta_{I,1}=3\lambda_{2}-\lambda_{1}-2\Omega_{1}.$ (12)

The resonance variable $\theta_{1}$ begins its libration roughly after $t\sim
3.5\times 10^{5}$ years being the time when the resonant capture happens (see
Figure 8). Interestingly, the center of libration is shifted from
$180^{\circ}$ to lower values between $100^{\circ}-150^{\circ}$. The
inclination of the inner planet gets excited around $t\sim 9.5\times 10^{5}$
years, that coincides to the libration of the resonance variable
$\theta_{I,1}$ (see Figure 9). After this epoch the system is both in 3:1
eccentricity-type and inclination-type MMR, see the right panels of Figures 8
and 9. Studying the left panels of these figures that display the plots
$(\lambda_{2}-\lambda_{1})$ versus $M_{1}$ for the eccentricity-type, and
$(\lambda_{2}-\lambda_{1})$ versus $(M_{1}+\omega_{1})$ for inclination-type
MMRs, one can see that the method _FAIR_ is able to identify both
eccentricity-type and inclination-type MMRs providing the correct resonance
variables. We note that the plots have been made for the whole timespan of the
numerical integration, also before the planets are got captured in resonances,
thus there are scattered points in the figures being not yet settled in the
stripes indicating the MMRs.

Figure 8: Two migrating giant planets captured into a 3:1 eccentricity-type
MMR. The method _FAIR_ is applied to the whole length of numerical
integration. Figure 9: Two migrating giant planets captured into a 3:1
inclination-type MMR. The method _FAIR_ is applied to the whole lenght of
numerical integration.

###  3.4 Temporary capture into a MMR

Figure 10: Temporary capture of 2007 RW10 into a 1:1 MMR, being a co-orbital
companion of Neptune for a while.

Finally, we also show that the method _FAIR_ is able to detect temporary
capture into a MMR. It can happen, for instance, that during their migration
planets enter into a MMR, see our results in the previous part. The left
panels of Figures 8 and 9 contain scattered points corresponding to that
epochs when the planets are not engaged into the 3:1 MMR. On the other hand it
is important that the scattered points should not repress those settled in the
stripes.

Beside the cases of migrating planets, there are examples of temporary capture
of bodies in co-orbital orbits, too. This happens to the asteroid 2007 RW10,
which was assumed as Neptune Trojan (de la Fuente Marcos & de la Fuente
Marcos, 2012). The right panel of Figure 10 shows the temporary libration of
the synodic longitude $\lambda^{\prime}-\lambda$, around $300^{\circ}$. The
left panel of this figure shows the $(\lambda^{\prime}-\lambda)$ versus $M$
plot, in which one parallel stripe to the horizontal axis becomes visible
around $\lambda^{\prime}-\lambda=300^{\circ}$ among the scattered points
indicating the case of the 1:1 MMR, see additionally Figure 5.

##  4 Conclusions

In this paper we presented a novel method which is suitable to easily decide
whether two planets are involved in some mean motion resonance, being either
of eccentricity or inclination type, without any a priori knowledge of its
character. Based on the geometrical meaning of the resonance variables, first
we described in detail how our method works. Next we demonstrated through a
few examples of real bodies of the Solar System and of exoplanetary systems,
involved in various MMRs, the straightforward applicability and efficiency of
the method _FAIR_. We have found that this method is able to detect MMRs also
in those cases when the involved bodies are only temporarily captured, such as
migrating pairs of planets or co-orbital companions to giant planets. Our
method can also serve as a practical help for future analysis of celestial
bodies which may be envolved in mean motion resonances.

## Acknowledgements

We thank the anonymous referee for her/his suggestions and comments helping us
to improve the work. This research has been supported by the Hungarian
National Research, Development and Innovation Office, NKFIH grant K-119993 and
the HAS Wigner RCP – GPU-Lab. Zs. Sándor thanks the support of the János
Bolyai Research Scholarship of the Hungarian Academy of Sciences.

## References

  * Brož & Vokrouhlický (2008) Brož M., Vokrouhlický D., 2008, MNRAS, 390, 715
  * Brown et al. (2005) Brown M. E., et al., 2005, ApJ, 632, L45
  * Cohen & Hubbard (1965) Cohen C. J., Hubbard E. C., 1965, AJ, 70, 10
  * Desort et al. (2008) Desort M., Lagrange A.-M., Galland F., Beust H., Udry S., Mayor M., Lo Curto G., 2008, A&A, 491, 883
  * Dormand & Prince (1978) Dormand J. R., Prince P. J., 1978, Celestial Mechanics, 18, 223
  * Galiazzo et al. (2013) Galiazzo M. A., Bazsó Á., Dvorak R., 2013, Planet. Space Sci., 84, 5
  * Goździewski & Migaszewski (2014) Goździewski K., Migaszewski C., 2014, MNRAS, 440, 3140
  * Goździewski et al. (2016) Goździewski K., Migaszewski C., Panichi F., Szuszkiewicz E., 2016, MNRAS, 455, L104
  * Hahn & Malhotra (2005) Hahn J. M., Malhotra R., 2005, AJ, 130, 2392
  * Laskar & Correia (2009) Laskar J., Correia A. C. M., 2009, A&A, 496, L5
  * Levison et al. (2008) Levison H. F., Morbidelli A., Van Laerhoven C., Gomes R., Tsiganis K., 2008, Icarus, 196, 258
  * Libert & Tsiganis (2009) Libert A.-S., Tsiganis K., 2009, MNRAS, 400, 1373
  * Marcy et al. (2001) Marcy G. W., Butler R. P., Fischer D., Vogt S. S., Lissauer J. J., Rivera E. J., 2001, ApJ, 556, 296
  * Mills et al. (2016) Mills S. M., Fabrycky D. C., Migaszewski C., Ford E. B., Petigura E., Isaacson H., 2016, Nature, 533, 509
  * Murray & Dermott (1999) Murray C. D., Dermott S. F., 1999, Solar system dynamics
  * Ortiz et al. (2005) Ortiz J. L., et al., 2005, Minor Planet Electronic Circulars, 2005-O36
  * Rivera et al. (2010) Rivera E. J., Laughlin G., Butler R. P., Vogt S. S., Haghighipour N., Meschiari S., 2010, ApJ, 719, 890
  * Sándor & Kley (2010) Sándor Z., Kley W., 2010, A&A, 517, A31
  * Sheppard et al. (2005) Sheppard S. S., Trujillo C. A., Marsden B. G., 2005, Minor Planet Electronic Circulars, 2005-U97
  * Tinney et al. (2006) Tinney C. G., Butler R. P., Marcy G. W., Jones H. R. A., Laughlin G., Carter B. D., Bailey J. A., O’Toole S., 2006, ApJ, 647, 594
  * Vogt et al. (2005) Vogt S. S., Butler R. P., Marcy G. W., Fischer D. A., Henry G. W., Laughlin G., Wright J. T., Johnson J. A., 2005, ApJ, 632, 638
  * de la Fuente Marcos & de la Fuente Marcos (2012) de la Fuente Marcos C., de la Fuente Marcos R., 2012, A&A, 545, L9


# MHD simulation of Solar Eruption from Active Region 11429 Driven by
Photospheric Velocity Field

Xinyi Wang  SIGMA Weather Group, State Key Laboratory for Space Weather,
National Space Science Center, Chinese Academy of Sciences,Beijing 100190, PR
China  College of Earth and Planetary Sciences, University of Chinese Academy
of Sciences, Beijing 100049, PR China  Chaowei Jiang Institute of Space
Science and Applied Technology, Harbin Institute of Technology, Shenzhen
518055, PR China<EMAIL_ADDRESS>Xueshang Feng SIGMA Weather Group, State
Key Laboratory for Space Weather, National Space Science Center, Chinese
Academy of Sciences,Beijing 100190, PR China  Institute of Space Science and
Applied Technology, Harbin Institute of Technology, Shenzhen 518055, PR China,
<EMAIL_ADDRESS>
###### Abstract

Data-driven simulation is becoming an important approach for realistically
characterizing the configuration and evolution of solar active regions,
revealing the onset mechanism of solar eruption events and hopefully achieving
the goal of accurate space weather forecast, which is beyond the scope of any
existing theoretical modelling. Here we performed a full 3D MHD simulation
using the data-driven approach and followed the whole evolution process from
quasi-static phase to eruption successfully for solar active region NOAA
11429. The MHD system was driven at the bottom boundary by photospheric
velocity field, which is derived by the DAVE4VM method from the observed
vector magnetograms. The simulation shows that a magnetic flux rope was
generated by persistent photospheric flow before the flare onset and then
triggered to erupt by torus instability. Our simulation demonstrates a high
degree of consistency with observations in the pre-eruption magnetic
structure, the time scale of quasi-static stage, the pattern of flare ribbons
as well as the time evolution of magnetic energy injection and total unsigned
magnetic flux. We further found that an eruption can also be initiated in the
simulation as driven by only the horizontal components of photospheric flow,
but a comparison of the different simulations indicates that the vertical flow
at the bottom boundary is necessary in reproducing more realistically these
observed features, emphasizing the importance of flux emergence during the
development of this AR.

Magnetohydrodynamic (MHD) — Sun: corona — Methods: numerical — Sun: magnetic
fields

##  1 Introduction

As driven by solar eruptions, the solar-terrestrial environment often
experiences variations, which are known as space weather, and forecasting the
space weather precisely is not only an important scientific topic but can also
avoid damage of the sensitive on-ground and space-based critical
infrastructures. Though many theoretical models have been proposed and
significant process has been made in understanding the triggering mechanism of
solar eruptions (Forbes, 2000; Chen, 2011; Schmieder et al., 2013; Priest,
2014), reproducing the whole life-span from quasi-static stage to eruption
using numerical models constrained and driven by observed vector magnetogram
possess unprecedented capabilities in revealing the onset mechanism of the
real eruption events, and can potentially be used for accurate space weather
forecast (Jiang et al., 2022; Jiang, 2022).

Previous study reproduced the whole process of energy accumulation and release
successfully (Jiang et al., 2016), showing an MHD system can be driven to
erupt by inputting time series of vector magnetograms at the bottom boundary
(B-driven). There are also other data-driven models, in which the evolution of
MHD system is driven by the electric field (E-driven, e.g., Cheung & DeRosa
2012; Hayashi et al. 2018; Pomoell et al. 2019; Price et al. 2019) or the
velocity field (V-driven, e.g., Hayashi et al. 2019; Guo et al. 2019; Liu et
al. 2019; He et al. 2020; Zhong et al. 2021) on the photosphere (bottom
boundary). Though the B-driven method can fully match the magnetogram, it will
introduce considerable errors of magnetic divergence from the bottom boundary.
This shortage will vanish in E-driven model, however, deriving both the
induction and potential components of the electric field on the photosphere is
not an easy task (Fisher et al., 2010, 2012, 2015, 2020), and the photospheric
flow also needs to be properly set to follow the Ohm’s law. In most of the
theoretical models of solar eruption, the key structure in favor of eruption
is assumed to be formed through the movement of the footpoints of magnetic
field lines (Moore & Labonte, 1980; Moore & Roumeliotis, 1992; Antiochos et
al., 1999; Lin & Forbes, 2000; Jiang et al., 2021b), which is driven by the
horizontal flow, and the vertical component of the photospheric flow will be
responsible for the flux emergence process. Therefore, with the photospheric
velocity field determined, the photospheric magnetic field can be generated
self-consistently. Furthermore, in the V-driven approach, there is no need to
solve the complex momentum equation at the bottom boundary (which is the most
time-consuming part in solving MHD equations). Due to these advantages, the
V-driven method has attracted many previous studies to focus on this topic.
For example, with the velocity field derived by DAVE4VM method (Schuck, 2008),
Hayashi et al. (2019) used the projected normal characteristics method to
update the physical variables other than velocity at the bottom boundary.
However, the total magnetic energy kept almost the same level of that of the
initial state without obvious magnetic energy injection. Jiang et al. (2021a)
updated the magnetic field by solving directly the magnetic induction equation
at the bottom boundary and their model can inject the magnetic energy from
bottom boundary successfully. He et al. (2020) drove the magnetic evolution by
inputting the DAVE4VM velocity field and vector magnetogram simultaneously.
The formation process of an magnetic flux rope (MFR) was obtained in their
model. Unfortunately, these simulations didn’t drive the system to erupt. The
only work we know that obtained an eruption by the velocity field derived from
observation was shown in Kaneko et al. (2021). The velocity field was derived
from the electric field on the photosphere using Ohm’s law and then was input
at the bottom boundary in their zero-beta model to drive the system to erupt.
Two eruptions they produced were identified from the evolution curves of the
magnetic and kinetic energies, however, the kinetic energy showed an overall
increase without obvious quasi-static evolution during the pre-eruption stage.
Before the first eruption, the kinetic energy was comparable with its peak
during the first eruption and the amount of the release of the magnetic energy
was also too small, which didn’t show the feature of a typical eruption event,
i.e., a large amount of magnetic energy was converted into kinetic energy
impulsively. As described above, the whole energy accumulation process from
quasi-static evolution (of typically tens of hours) to impulsive eruption has
not been realized in a self-consistent way using the V-driven MHD model, and
this is one of the motivations of this work.

In this Letter, we applied a V-driven model to investigate the evolution and
eruption of a well-studied active region (AR) NOAA 11429. Previous studies
found persistent shearing flow and flux cancellation near the main polarity
inversion line (PIL) of this AR (Shimizu et al., 2014; Zheng et al., 2017),
which were suggested to be responsible for the eruptions on 2012 March 7, 9
and 10 (Dhakal et al., 2018, 2020). An analysis of the MFR reconstructed from
vector magnetograms using the nonlinear force-free field (NLFFF) model
suggests that the homologous eruptions are triggered by torus instability (TI)
of the MFR (Chintzoglou et al., 2015). Zhang et al. (2021) suggested the
helical kink instability may also take effect. Nevertheless, whether these
mechanisms were at work requires to be further studied using dynamic modeling
of this AR evolution and eruption, which is absent in all the previous studies
and is the other motivation of this work. In our simulation, the dynamic
process from the beginning of 2012 March 4 to the eruption on March 5 in this
AR was reproduced self-consistently as driven by the photospheric velocity
derived from vector magnetogram using DAVE4VM method. Our simulation shows
that an MFR was generated near the main PIL before the flare onset and the
initiation of this eruption event depended mainly on TI of the preformed MFR.

##  2 Data and Model

AR 11429 showed a complex $\beta\gamma\delta$ configuration and is very flare-
productive, which has produced 3 X-class flares from 2012 March 5 to 7. It
first appeared on the eastern solar limb on March 4 and was located on the
eastern part of the solar disk before March 8. During this period, the AR kept
developing as characterized by the increasing total unsigned magnetic flux,
indicating the obvious flux emergence by vertical flow on the photosphere.
Since it was the first X-class flare of this AR, here we focus on the
initiation process of the X1.1 flare (which is accompanied with a halo CME
moving at a speed of 1531 $\rm km~{}s^{-1}$) around 04:00 UT on March 5 as
shown in the white box in Figure 1A. Before the flare onset, there was
persistent shearing flow near the main PIL (Figure 1B) and as a result, the
horizontal magnetic field there was highly sheared (Figure 1C), which
indicates that a large amount of free magnetic energy is stored ready for an
eruption. An hot loop first erupted away as shown in AIA 94 $\rm\AA$ (Lemen et
al., 2012) at around 03:31 UT (Figure 1D and E) and after that a pair of hook
shape flare ribbons appeared near the main PIL in AIA 1600 $\rm\AA$ at 03:36
UT (Figure 1F), i.e., the flare event started.

To understand the formation of the pre-eruptive coronal magnetic field and the
triggering mechanism of this eruption, we used the DARE-MHD model (Jiang et
al., 2016) to study the dynamic evolution of this AR. For saving the
computational time, the strength of magnetic field from the magnetograms were
reduced by a factor of 25 before being input into our code. The initial plasma
density was set as a hydrostatic isothermal model with a fixed temperature as
$T=1\times 10^{6}$ K and a modified solar gravity (Jiang et al., 2021b) to get
a plasma background that mimics the real environment in the solar corona
basing on two key parameters, the plasma $\beta$ and the $\rm Alfv\acute{e}n$
speed (in particular, the minimum $\beta$ is $6.3\times 10^{-4}$ and maximum
$\rm Alfv\acute{e}n$ speed $V_{\rm A}\sim 4800$ km s-1 in the final
equilibrium we obtained below). We chose to use the magnetograms of HMI SHARP
data set (Schou et al., 2012; Bobra et al., 2014) at 00:00UT on March 4 to
reconstruct the initial magnetic field since there were no obvious MFR at that
time. This magnetogram was smoothed first using Gaussian smoothing with FWHM
of 6 pixels and a NLFFF model from this smoothed magnetogram was extrapolated
by our CESE-NLFFF-MHD code (Jiang & Feng, 2012, 2013). Since the code (like
many other NLFFF codes) does not give a perfect force-free solution but with
residual Lorentz forces, we input the extrapolated field into the DARE-MHD
model, along with the initial background plasma, to let the MHD system relax
until the kinetic and magnetic energies were almost unchanged, i.e., an MHD
equilibrium was obtained, and the initial state was ready.

At the bottom boundary, we solve the magnetic induction equation to update the
magnetic field with the velocity field derived by DAVE4VM method to update the
magnetic field on the photosphere. The DAVE4VM velocity was strengthened by a
factor of 13.7 (determined by the ratio of the time series magnetograms’
original time cadence as 720 seconds to the time cadence in our simulation as
$0.5\times 105$ seconds, i.e., $\frac{720}{0.5\times 105}$) to speed up our
simulation and thus the time scale of quasi-static evolution prior to the
eruption onset is shorten by the same times. At the side and top boundaries,
all the variables are extrapolated from the neighboring inner points with zero
gradient along the normal direction of the boundary surface and the normal
component of magnetic field is further modified by the divergence-free
condition. The Powell-source terms and diffusion control terms was used to
deal with the divergence error of the magnetic field as described in Jiang et
al. (2010). We set the computational domain sufficiently large as $[-368,368]$
Mm in both $x$ and $y$ direction and $[0,736]$ Mm in $z$ direction with grid
resolution varies from $1^{\prime\prime}$ to $8^{\prime\prime}$ using adaptive
mesh refinement. The highest resolution is used mainly for the regions with
strong magnetic gradients and current density, in particular, the current
sheets (CSs). Explicit value of magnetic resistivity was not used in our
simulation and the magnetic reconnection was controlled by the resistivity of
the numerical method only, which mimicked the low-resistivity plasma better.
As the total unsigned flux of this AR kept increasing before the flare onset,
we energized the MHD system by full 3D DAVE4VM velocity field
($v^{D}_{x},v^{D}_{y},v^{D}_{z}$) (will be referred to as V3D simulation) and
horizontal component ($v^{D}_{x},v^{D}_{y},0$) (V2D) respectively, and a
comparison of the results for these two simulations will show the importance
of flux emergence through the vertical flow on the photosphere.

##  3 Results

###  3.1 Overall Process

The evolution curves of the total magnetic and kinetic energies in the
computational domain as well as the total unsigned magnetic flux at the bottom
boundary are shown in Figure 2. For the ‘V3D’ simulation, as driven by the
time-series velocity field (V3D) for a time duration of 150 minutes, the total
magnetic energy in our simulation model experienced an overall increase
firstly and then a rapid decrease, which is associated with an eruption event.
The eruption can be identified from the energy evolution with onset time
$t_{\rm E,V3D}=120$ minutes.

At the very beginning from $t=0$ to $t=22$ minutes, the magnetic energy
injection curve (black dashed line, which is computed by the time integration
of the total Poynting flux of the ‘V3D’ simulation at the bottom surface)
matches well with the solid ‘V3D’ curve (the magnetic energy increase of the
‘V3D’ simulation) and ‘OB’ curve (the magnetic energy injection computed by
DAVE4VM velocity and magnetograms) in Figure 2A. However, in the time duration
of $t\in[97,122]$ minutes, the total magnetic energy (blue solid line) is
higher than the magnetic energy injection. Such an unphysical mismatch of the
inputted energy from the boundary and the cumulative energy in the volume is
likely owing to the insufficient resolution for the bottom boundary, since the
magnetic field at the bottom boundary has accumulated a very large gradient in
this phase (see Discussions). The kinetic energy keeps a very low value of
around $10^{-3}E_{\rm p_{0}}$ (which is the potential field energy
corresponding to the magnetic field at t=0) before the major eruption begins
at $t_{\rm E,V3D}=120$ minutes when the magnetic energy reaches about 1.85
$E_{\rm p_{0}}$. The magnetic energy decreases to about 1.7 $E_{\rm p_{0}}$ at
the peak of the total kinetic energy (i.e., $E_{\rm k}=5.6\times 10^{-2}E_{\rm
p_{0}}$), and keeps decreasing to 1.45 $E_{\rm p_{0}}$ in total through this
eruption (0.4 $E_{\rm p_{0}}$ free energy loss). That is, about one third of
the magnetic energy loss has been converted to kinetic energy in 10 minutes.
If multiplied by a factor of 13.7 determined by the rate of speeding up in our
velocity-driven simulation, the quasi-static evolution time of ‘V3D’ run is
27.4 hours. This time scale is very close to the observation one which is 27.6
hours.

We have also driven the simulation by the horizontal velocity (V2D) and found
that it can also produce an eruption with rather similar onset time. However,
comparing the different curves of magnetic energy evolution in Figure 2A and
the curves of total unsigned magnetic flux in Figure 2B respectively, the blue
solid lines labeled by ‘V2D’ are obvious lower than those from the ’V3D’
simulation. This clearly shows that though horizontal velocity can also drive
the field to erupt (with a delayed onset time for about 6 minutes compared
with ‘V3D’), the vertical velocity $v^{D}_{z}$ is necessary in accounting for
the larger increase of the total magnetic energy and total unsigned flux as
shown in observations, therefore leading to a stronger eruption. The evolution
of the simulated magnetic energy, total unsigned flux and the time scale of
the ‘V3D’ run before the major eruption are more consistent with observations,
showing the importance of the vertical photospheric plasma flow in the
numerical modeling of solar eruptions. Our simulated magnetic structure (the
first panel in Figure 3E) has reasonable consistency with observations in the
pre-eruption image of AIA 171 $\rm\AA$ (the second panel in Figure 3E), and
the synthetic image of coronal emission from current density of our simulation
(the last panel in Figure 3E) is reasonable consistent with the image of AIA
131 $\rm\AA$ (the third panel in Figure 3E). Also the quasi separatrix layers
(QSLs), where the magnetic reconnection is most likely to take place and thus
represent the position of the flare ribbons (Titov et al., 2002; Liu et al.,
2016), at the bottom boundary (Figure 4F) has approximately the same patterns
as the flare ribbons in AIA 1600 $\rm\AA$ (Figure 4E). These results confirm
the validity of our V-driven DARE-MHD model as well as the DAVE4VM method.

###  3.2 Eruption Initiation Mechanism

Since the actual velocity field must contain $v_{z}$, here we analyzed the
‘V3D’ run to study the onset mechanism of this eruption. As we can see, a
group of twisted field lines (represented by the blue solid lines in Figure 3A
and B) formed and was embedded in the surrounding shear arcades, which is
similar to an MFR in morphology. In addition, the ejection of the hot loop
(also called hot channel) as observed in AIA 94 $\rm\AA$ (shown by the white
arrows in Figure 1D and E) also suggested the existence of an MFR (Cheng et
al., 2017) before the flare onset.

To identify the formation of MFR before the eruption, we calculated the QSLs
and the twist number (Berger & Prior, 2006) of our simulated magnetic fields,
which are given in Figure 3D, Figure 5A and B, respectively. As shown in the
third panel of Figure 3D, a strong QSL appeared near the core field and grew
up to be a QSL ring in the last panel, which separated the MFR and the
background magnetic field, thus representing the existence of an MFR in
topology. The QSL ring intersects itself below the MFR, forming a typical
hyperbolic flux tube (HFT), where the CS developed and magnetic reconnection
took place subsequently to further drive the eruption (Jiang et al., 2018).
The isosurface of $T_{w}=-1$, which represents the position and shape of the
MFR, are shown in Figure 5A and B respectively. It became larger and higher,
illustrating more magnetic flux were twisted as driven by persistent
photospheric flow. However, the isosurface of $T_{w}=-2$, which is the lower
threshold of kink instability (KI) according to the statistics of Duan et al.
(2019), is barely visible and KI may take little effect.

Among different triggering mechanisms, TI is considered as an efficient way of
MFR eruption. In Figure 5D, we plot the key controlling parameter that
determines the onset of TI, i.e., the decay index of the strapping field. The
variation of the decay index $n$ was calculated along the $z$ direction of the
overlying field at $t=119$ minutes before the eruption onset. Since the
potential field is not always the good approximation of the strapping field
(especially when the strapping field is highly sheared), we calculated the
decay index of our simulated field instead and plotted the decay index of the
corresponding potential field at the same time for comparison. The critical
height (above which $n>1.5$) of the simulated field is located at 60 Mm as
labeled by the vertical black dashed line in Figure 5D. When the MFR was just
formed, it is small and it’s apex was at a low height (about 40 Mm in Figure
5A). About 14 minutes later, it had grown up to be a huge twist structure and
entered the unstable zone (above 60 Mm as shown in Figure 5B), after which it
erupted violently.

The formation and the evolution of the MFR, i.e., there was an preformed MFR
and the MFR erupted after it entered the unstable zone, strongly suggest that
TI is the triggering mechanism of this eruption. To further test this
assumption, we used the V3D-driven data at $t=100$ minutes as the initial
condition and ran our code without bottom velocity driven (by setting all the
three velocity components as zero at the bottom boundary, thus referred to as
‘V0D’) to see whether the magnetic field will erupt. The evolution of energies
of this ‘V0D’ run are shown as the red solid lines in Figure 2A. During the
time duration of $t\in[100,124]$ minutes, the magnetic energy kept almost
unchanged, while the toroidal flux of the MFR (defined as $\int_{s}B_{z}ds$,
where $s$ denotes the region of $T_{w}<-1$ and $B_{z}<0$ at the bottom
boundary) kept increasing as shown in Figure 5C. In ‘V3D’ simulation, the MFR
(the isosurface of $T_{w}=-1$) became larger and higher than the critical
height before the eruption onset (Figure 5A and B) as the toroidal flux
increasing (Figure 5C), after which TI took effect and finally led to a strong
eruption. Similarly, the toroidal flux of ‘V0D’ simulation increased before
and decreased after the eruption (Figure 5C), showing the slow reconnection
can took place spontaneously without velocity driving and made the MFR larger
until the eruption started at $t_{\rm E,V0D}=124$ minutes as identified from
the energy evolution curve. Since the current density was weaker and the
current layer was thicker (as shown in the third panel in Figure 3C) than a
true CS (the first panel in Figure 4C), the magnetic gradient was lower in the
current layer than in the CS, and thus the diffusivity was relatively uniform,
which only allowed slow reconnection (i.e., a slow dissipation of the current)
to take place without impulsive energy release (i.e., not a Pescheck-type
reconnection, Yokoyama & Shibata 1994). The eruption of this ‘V0D’ run has the
similar onset time and strength with the ‘V3D’ run, which further indicates
that after $t=100$ minutes in ‘V3D’ run, the velocity field on the bottom
boundary is not necessary and the instability is sufficient to trigger the
eruption. It follows the basic developing stage of TI, i.e., the rising MFR
stretched the overlying field and consequently the flare CS formed below the
MFR. Since we didn’t use any explicit value of magnetic resistivity, and since
the CS formed in a dynamic way, the width of CS could become very thin to
trigger the fast reconnection easily as pointed out by Jiang et al. (2021b),
which further drives the eruption. Based on the analysis of our simulation and
the supplementary numerical experiment (‘V0D’ simulation), along with the
consistency between our simulation results and observations that have been
shown above, we conclude that TI is the initiation mechanism of the X1.1 flare
in AR 11429 on 2012 March 5.

##  4 Conclusions and Discussions

We have carried out a full 3D MHD simulation of an X-class flare eruption
event on 2012 March 5 in AR 11429 using the V-driven DARE-MHD model. An MHD
equilibrium was obtained by relaxing the NLFFF reconstructed by the CESE-MHD-
NLFFF code and was set as the initial condition. Then the initial state was
driven to evolve by the DAVE4VM velocity field on the bottom boundary of our
simulation box. The analysis of the quasi-static evolution stage before the
eruption onset shows the gradual formation of an MFR above the main PIL. When
the MFR first appeared, it was relatively low and then it grew up into the
torus unstable region, after which TI was triggered. Then the energy
conversion process was accomplished by reconnection in the CS that is formed
due to the stretching effect of the erupting MFR. The images of SDO
observation and the important physical quantities computed from the
magnetograms are reasonably consistent with our result in the pre-flare
magnetic structure, the morphology of flare ribbons, evolution curves of the
magnetic and kinetic energies as well as the total unsigned magnetic flux. The
time duration of the quasi-static evolution process before the simulated
eruption is also very close to the actual time scale before the flare onset.

Nevertheless, our simulation does not reproduce accurately the evolution of
magnetic field as shown in the observed magnetograms. As can be seen in the
last panel of Figure 3A, part of the magnetic flux was transported to and
concentrated at the edge of the AR, while this pileup was not observed in HMI
magnetograms. One likely reason for this flux pileup is that, the DAVE4VM
method only solves the normal component of the magnetic induction equation in
the least-square sense (Lumme et al., 2019), which may not reproduce the
velocity precisely at every point and can only be used to approximate the
overall distribution and evolution of magnetic flux in the AR. In addition,
the magnetic flux on the photosphere should be dispersed and dissipated by
granular and supergranular convection (Wang et al., 1989) as well as small
scale turbulent diffusion in practical situation. Therefore, without dealing
properly with these effects, the flux pileup will be more obvious than
observations and as a result the magnetic energy injection rate of ‘V3D’ run
will be overall higher than the ‘OB’ one as shown in Figure 2A (comparing the
black dashed line and the black solid line). This unrealistic flux pileup may
also contribute to the overshoot of the magnetic energy increase (Figure 2A,
the mismatch of the black dashed line and the green line before the eruption),
because it results in a very large magnetic gradient at the bottom boundary,
thus even the finest spatial resolution in our simulation was inadequate to
capture the too high gradient during the time duration close to the eruption
onset (i.e., $t\in[97,122]$ minutes). The insufficient grid resolution led to
the considerable numerical error there, which thus makes the magnetic energy
increase higher than the magnetic energy injection in ‘V3D’ simulation. The
proper settings of numerical diffusion and grid resolution, along with the
improved methods for deriving the photospheric flow will need to be considered
in future works for a more accurate data-driven simulation of solar eruptions.

To summarize, our simulation shows that, besides inputting the time series of
vector magnetogram, the numerical model of solar corona can be driven to erupt
by inputting the time series velocity field at the bottom boundary, in which
the driver, i.e., the bottom velocity field is derived from time series of
vector magnetograms. The numerical model we established here shows the
possibility of driving the evolution of solar corona using different physical
variables, which offers a straightforward way for revealing the eruption
mechanisms in real events. Thus, such model has a great potentiality in
forecasting the onset time as well as the strength of solar eruptions and
evaluating the quantitative impact on space weather.

This work is jointly supported by the National Natural Science Foundation of
China (NSFC 41731067, 42030204, 42174200), the Fundamental Research Funds for
the Central Universities (grant No. HIT.OCEF.2021033), and Shenzhen Science
and Technology Program (grant Nos. RCJC20210609104422048 and
JCYJ20190806142609035). The computational work was carried out on TianHe-1(A),
National Supercomputer Center in Tianjin, China, and we thank Jun Chen for his
informative and helpful discussions.

## References

  * Antiochos et al. (1999) Antiochos, S. K., DeVore, C. R., & Klimchuk, J. A. 1999, ApJ, 510, 485, doi: 10.1086/306563
  * Berger & Prior (2006) Berger, M. A., & Prior, C. 2006, Journal of Physics A Mathematical General, 39, 8321, doi: 10.1088/0305-4470/39/26/005
  * Bobra et al. (2014) Bobra, M. G., Sun, X., Hoeksema, J. T., et al. 2014, Sol. Phys., 289, 3549, doi: 10.1007/s11207-014-0529-3
  * Chen (2011) Chen, P. F. 2011, Living Reviews in Solar Physics, 8, 1, doi: 10.12942/lrsp-2011-1
  * Cheng et al. (2017) Cheng, X., Guo, Y., & Ding, M. 2017, Science China Earth Sciences, 60, 1383, doi: 10.1007/s11430-017-9074-6
  * Cheung & DeRosa (2012) Cheung, M. C. M., & DeRosa, M. L. 2012, ApJ, 757, 147, doi: 10.1088/0004-637X/757/2/147
  * Chintzoglou et al. (2015) Chintzoglou, G., Patsourakos, S., & Vourlidas, A. 2015, ApJ, 809, 34, doi: 10.1088/0004-637X/809/1/34
  * Dhakal et al. (2018) Dhakal, S. K., Chintzoglou, G., & Zhang, J. 2018, ApJ, 860, 35, doi: 10.3847/1538-4357/aac028
  * Dhakal et al. (2020) Dhakal, S. K., Zhang, J., Vemareddy, P., & Karna, N. 2020, ApJ, 901, 40, doi: 10.3847/1538-4357/abacbc
  * Duan et al. (2019) Duan, A., Jiang, C., He, W., et al. 2019, ApJ, 884, 73, doi: 10.3847/1538-4357/ab3e33
  * Fisher et al. (2020) Fisher, G. H., Kazachenko, M. D., Welsch, B. T., et al. 2020, ApJS, 248, 2, doi: 10.3847/1538-4365/ab8303
  * Fisher et al. (2012) Fisher, G. H., Welsch, B. T., & Abbett, W. P. 2012, Sol. Phys., 277, 153, doi: 10.1007/s11207-011-9816-4
  * Fisher et al. (2010) Fisher, G. H., Welsch, B. T., Abbett, W. P., & Bercik, D. J. 2010, ApJ, 715, 242, doi: 10.1088/0004-637X/715/1/242
  * Fisher et al. (2015) Fisher, G. H., Abbett, W. P., Bercik, D. J., et al. 2015, Space Weather, 13, 369, doi: 10.1002/2015SW001191
  * Forbes (2000) Forbes, T. G. 2000, J. Geophys. Res., 105, 23153, doi: 10.1029/2000JA000005
  * Guo et al. (2019) Guo, Y., Xia, C., Keppens, R., Ding, M. D., & Chen, P. F. 2019, ApJ, 870, L21, doi: 10.3847/2041-8213/aafabf
  * Hayashi et al. (2018) Hayashi, K., Feng, X., Xiong, M., & Jiang, C. 2018, The Astrophysical Journal, 855, 11, doi: 10.3847/1538-4357/aaacd8
  * Hayashi et al. (2019) Hayashi, K., Feng, X., Xiong, M., & Jiang, C. 2019, ApJ, 871, L28, doi: 10.3847/2041-8213/aaffcf
  * He et al. (2020) He, W., Jiang, C., Zou, P., et al. 2020, ApJ, 892, 9, doi: 10.3847/1538-4357/ab75ab
  * Jiang (2022) Jiang, C. 2022, Reviews of Geophysics and Planetary Physics, 53, 497, doi: 10.19975/j.dqyxx.2022-022
  * Jiang et al. (2021a) Jiang, C., Bian, X., Sun, T., & Feng, X. 2021a, Frontiers in Physics, 9, 224, doi: 10.3389/fphy.2021.646750
  * Jiang & Feng (2012) Jiang, C., & Feng, X. 2012, ApJ, 749, 135, doi: 10.1088/0004-637X/749/2/135
  * Jiang & Feng (2013) —. 2013, ApJ, 769, 144, doi: 10.1088/0004-637X/769/2/144
  * Jiang et al. (2022) Jiang, C., Feng, X., Guo, Y., & Hu, Q. 2022, The Innovation, 3, 100236, doi: https://doi.org/10.1016/j.xinn.2022.100236
  * Jiang et al. (2010) Jiang, C., Feng, X., Zhang, J., & Zhong, D. 2010, Sol. Phys., 267, 463, doi: 10.1007/s11207-010-9649-6
  * Jiang et al. (2016) Jiang, C., Wu, S. T., Feng, X., & Hu, Q. 2016, Nature Communications, 7, 11522, doi: 10.1038/ncomms11522
  * Jiang et al. (2018) Jiang, C., Zou, P., Feng, X., et al. 2018, The Astrophysical Journal, 869, 13, doi: 10.3847/1538-4357/aaeacc
  * Jiang et al. (2021b) Jiang, C., Feng, X., Liu, R., et al. 2021b, Nature Astronomy, 5, 1126, doi: 10.1038/s41550-021-01414-z
  * Kaneko et al. (2021) Kaneko, T., Park, S.-H., & Kusano, K. 2021, ApJ, 909, 155, doi: 10.3847/1538-4357/abe414
  * Lemen et al. (2012) Lemen, J. R., Title, A. M., Akin, D. J., et al. 2012, Sol. Phys., 275, 17, doi: 10.1007/s11207-011-9776-8
  * Lin & Forbes (2000) Lin, J., & Forbes, T. G. 2000, J. Geophys. Res., 105, 2375, doi: 10.1029/1999JA900477
  * Liu et al. (2019) Liu, C., Chen, T., & Zhao, X. 2019, A&A, 626, A91, doi: 10.1051/0004-6361/201935225
  * Liu et al. (2016) Liu, R., Kliem, B., Titov, V. S., et al. 2016, ApJ, 818, 148, doi: 10.3847/0004-637X/818/2/148
  * Lumme et al. (2019) Lumme, E., Kazachenko, M. D., Fisher, G. H., et al. 2019, Sol. Phys., 294, 84, doi: 10.1007/s11207-019-1475-x
  * Moore & Labonte (1980) Moore, R. L., & Labonte, B. J. 1980, in Solar and Interplanetary Dynamics, ed. M. Dryer & E. Tandberg-Hanssen, Vol. 91, 207–210
  * Moore & Roumeliotis (1992) Moore, R. L., & Roumeliotis, G. 1992, Triggering of Eruptive Flares - Destabilization of the Preflare Magnetic Field Configuration, ed. Z. Svestka, B. V. Jackson, & M. E. Machado, Vol. 399, 69, doi: 10.1007/3-540-55246-4_79
  * Pomoell et al. (2019) Pomoell, J., Lumme, E., & Kilpua, E. 2019, Sol. Phys., 294, 41, doi: 10.1007/s11207-019-1430-x
  * Price et al. (2019) Price, D. J., Pomoell, J., Lumme, E., & Kilpua, E. K. J. 2019, A&A, 628, A114, doi: 10.1051/0004-6361/201935535
  * Priest (2014) Priest, E. 2014, Magnetohydrodynamics of the Sun, doi: 10.1017/CBO9781139020732
  * Schmieder et al. (2013) Schmieder, B., Démoulin, P., & Aulanier, G. 2013, Advances in Space Research, 51, 1967, doi: 10.1016/j.asr.2012.12.026
  * Schou et al. (2012) Schou, J., Scherrer, P. H., Bush, R. I., et al. 2012, Sol. Phys., 275, 229, doi: 10.1007/s11207-011-9842-2
  * Schuck (2008) Schuck, P. W. 2008, ApJ, 683, 1134, doi: 10.1086/589434
  * Shimizu et al. (2014) Shimizu, T., Lites, B. W., & Bamba, Y. 2014, Publications of the Astronomical Society of Japan, 66, S14, doi: 10.1093/pasj/psu089
  * Titov et al. (2002) Titov, V. S., Hornig, G., & Démoulin, P. 2002, Journal of Geophysical Research (Space Physics), 107, 1164, doi: 10.1029/2001JA000278
  * Wang et al. (1989) Wang, Y. M., Nash, A. G., & Sheeley, N. R., J. 1989, Science, 245, 712, doi: 10.1126/science.245.4919.712
  * Yokoyama & Shibata (1994) Yokoyama, T., & Shibata, K. 1994, ApJ, 436, L197, doi: 10.1086/187666
  * Zhang et al. (2021) Zhang, Y., Bastian, T. S., Liu, J. H., et al. 2021, ApJ, 910, 40, doi: 10.3847/1538-4357/abded6
  * Zheng et al. (2017) Zheng, J.-C., Yang, Z.-L., Guo, J.-P., et al. 2017, Research in Astronomy and Astrophysics, 17, 081, doi: 10.1088/1674-4527/17/8/81
  * Zhong et al. (2021) Zhong, Z., Guo, Y., & Ding, M. D. 2021, Nature Communications, 12, 2734, doi: 10.1038/s41467-021-23037-8

Figure 1: ($\mathbf{A}$): A full-disk composite image of SDO AIA 171 $\rm\AA$
and HMI magnetogram. The white box shows the AR 11429. ($\mathbf{B}$): The
vectors denote the time average horizontal velocity from 2012 March 4 00:00 UT
to March 5 03:36 UT on the photosphere. The background shows the magnetic flux
distribution, in which black and white represents the negative and positive
polarity respectively. ($\mathbf{C}$): The vectors denote the transverse
magnetic field on the photosphere. ($\mathbf{D}$): Pre-flare image of AIA
94$\rm\AA$. The white arrow shows the erupting structure. ($\mathbf{E}$): Same
as $\mathbf{D}$, but at a time closer to eruption before the flare onset.
($\mathbf{F}$): The first image of flare ribbons observed in AIA 1600$\rm\AA$.
Letter P labels the positive ribbon and N labels the negative ribbon. Figure
2: ($\mathbf{A}$): Magnetic and kinetic energy evolution during the whole
process. The unit of $x$-axis is in minutes. The solid curves in different
colors denote the corresponding evolution of different energies and the
vertical dashed lines in different colors denote the eruption onset time of
different runs. The vertical black solid line shows the flare onset time in
observations. The black curve ‘OB’ denotes the magnetic energy injection
computed by DAVE4VM velocity and magnetograms, thus representing the actual
magnetic energy evolution. In our simulation, magnetic energy of the potential
field of the initial magnetogram is $E_{\rm p_{0}}=2.52\times 10^{30}$ erg
(which is a fixed value and is used for normalization here). This value should
be multiplied by a factor of 625=$25^{2}$ if scaling to the realistic value,
thus being $1.57\times 10^{33}$ erg. ($\mathbf{B}$): Time evolution of the
total unsigned magnetic flux in our simulations. The ‘OB’ curve is computed by
the magnetograms and the other labels ‘V3D’, ‘V2D’ and the vertical black
solid line have the same meanings as in $\mathbf{A}$. All curves are
normalized by their initial value at $t=0$. The corresponding animation is
attached for the V3D run, and it starts at $t=0$ and ends at $t=147$ minutes
in simulation time, which corresponds to 33.6 hours in real-time duration.
Panel $\mathbf{a}$ in the animation corresponds to rows $\mathbf{A}$ in Figure
3 and 4. Panel $\mathbf{b}$ corresponds to the rows $\mathbf{C}$ in Figure 3
and 4 in time series. Panel $\mathbf{c}$ shows the corresponding evolution as
in the row $\mathbf{D}$ in Figure 4, and panel $\mathbf{d}$ represents the
evolution of QSLs as in Figure 4F. We only plotted the kinetic energy in panel
$\mathbf{e}$, with the vertical solid line showing the current simulation
time; the magnetic flux evolution is not included in the animation. The
cadence between each snapshot used in the animation is 210 s in the quasi-
static period ($t\in[0,119]$ minutes) and 21 s in the eruption period
($t\in(119,147]$ minutes) in simulation time, respectively. Figure 3: The
magnetic evolution during the pre-eruption stage. ($\mathbf{A}$): Top view of
3 bunches of field lines with the footpoints approximately following the
movement of the magnetic flux. The white solid line denotes the position of
the slices in $\mathbf{C}$ and in Figure 4C and D. The red solid line segment
shows the position of the slices in $\mathbf{D}$, which is almost
perpendicular to the PIL and at the middle of MFR. ($\mathbf{B}$): Side view
of the same 3D magnetic field lines as in $\mathbf{A}$. ($\mathbf{C}$):
Vertical cross section of the current layer near the main PIL. The position of
these slices are labeled by the white solid line in $\mathbf{A}$.
($\mathbf{D}$): Magnetic squashing factor $Q$ on the same cross section of
$\mathbf{C}$ and the region with high $Q$ denotes the QSLs. The position of
these slices are labeled by the red solid line segment in $\mathbf{A}$.
($\mathbf{E}$): The first 2 panels show the comparison between our simulation
and observations in AIA 171 $\rm\AA$. The last 2 panels show comparison
between the synthetic image of coronal emission from current density of our
simulation and observations in AIA 131 $\rm\AA$. Row $\mathbf{A}$ is shown at
panel $\mathbf{a}$, and rows $\mathbf{C}$ is shown at panel $\mathbf{c}$ in
the animation of Figure 2, respectively. Figure 4: The magnetic evolution of
the eruption stage. ($\mathbf{A}$), ($\mathbf{B}$) and ($\mathbf{C}$): All
settings are the same as in Figure 3A, B and C, except the footpoints of the
magnetic field lines are fixed and the current layer in Figure 3C turns into a
current sheet here. ($\mathbf{D}$): Outflows at the position of current sheet.
($\mathbf{E}$): Projection corrected images of the flare ribbons observed in
AIA 1600 $\rm\AA$. The letter P denotes the positive ribbon and N denotes the
negative ribbon. ($\mathbf{F}$): Evolution of the bottom QSL of our
simulation, where sgn$(B_{z})$ denotes the sign of $B_{z}$. The letter P
denotes the positive QSL and N denotes the negative QSL. The slices in
$\mathbf{C}$ and $\mathbf{D}$ are located at the position labeled by the white
solid line in the first panel of Figure 3A. Row $\mathbf{A}$ is shown at panel
$\mathbf{a}$, and rows $\mathbf{C}$ and $\mathbf{D}$ are shown at panel
$\mathbf{c}$ and panel $\mathbf{d}$ in the animation of Figure 2,
respectively. The QSL evolution in row $\mathbf{F}$ is shown at panel
$\mathbf{b}$. Figure 5: ($\mathbf{A}$): The isosurface of $T_{w}=-1$ and
$T_{w}$ is the twist number of the magnetic field in ‘V3D’ simulation. The red
‘X’ labels the position where we calculate the decay index. ($\mathbf{B}$):
Same as $\mathbf{A}$ but at a different time. ($\mathbf{C}$): The green solid
curve denotes the relative toroidal flux of the MFR in ‘V3D’ simulation and
the red solid curve in ‘V0D’ simulation with both curves normalized by their
initial values respectively. Two vertical dashed lines in different colors
denotes the onset time of ‘V0D’ and ‘V3D’ simulations, respectively. The
toroidal flux is defined as $\int_{s}B_{z}ds$, where ‘s’ is the region of
$T_{w}<-1$ and $B_{z}<0$ at the bottom boundary. ($\mathbf{D}$): The decay
index $n$ of the magnetic field in ‘V3D’ simulation (solid curve) and the
corresponding potential field (dashed curve) at the same time of $t=119$
minutes. The horizontal dashed line denotes the critical value of $n=1.5$ and
the vertical dashed line denotes the critical height, above which $n>1.5$.