mydate []

# The role of non-axisymmetry of magnetic flux rope in constraining solar
eruptions

Ze Zhong School of Astronomy and Space Science, Nanjing University, Nanjing
210023, People’s Republic of China Key Laboratory for Modern Astronomy and
Astrophysics (Nanjing University), Ministry of Education, Nanjing 210023,
People’s Republic of China. ${}^{\textsuperscript{\Letter}}$email:
<EMAIL_ADDRESS><EMAIL_ADDRESS>Yang Guo${}^{\textsuperscript{\Letter}}$
School of Astronomy and Space Science, Nanjing University, Nanjing 210023,
People’s Republic of China Key Laboratory for Modern Astronomy and
Astrophysics (Nanjing University), Ministry of Education, Nanjing 210023,
People’s Republic of China. ${}^{\textsuperscript{\Letter}}$email:
<EMAIL_ADDRESS><EMAIL_ADDRESS>M. D. Ding${}^{\textsuperscript{\Letter}}$
School of Astronomy and Space Science, Nanjing University, Nanjing 210023,
People’s Republic of China Key Laboratory for Modern Astronomy and
Astrophysics (Nanjing University), Ministry of Education, Nanjing 210023,
People’s Republic of China. ${}^{\textsuperscript{\Letter}}$email:
<EMAIL_ADDRESS><EMAIL_ADDRESS>
## Abstract

Whether a solar eruption is successful or failed depends on the competition
between different components of the Lorentz force exerting on the flux rope
that drives the eruption. The present models only consider the strapping force
generated by the background magnetic field perpendicular to the flux rope and
the tension force generated by the field along the flux rope. Using the
observed magnetic field on the photosphere as a time-matching bottom boundary,
we perform a data-driven magnetohydrodynamic simulation for the 30 January
2015 confined eruption and successfully reproduce the observed solar flare
without a coronal mass ejection. Here we show a Lorentz force component,
resulting from the radial magnetic field or the non-axisymmetry of the flux
rope, which can essentially constrain the eruption. Our finding contributes to
the solar eruption model and presents the necessity of considering the
topological structure of a flux rope when studying its eruption behaviour.

## Introduction

Magnetic flux ropes (MFRs) are a bundle of twisted field lines with electric
currents flowing inside. They play a critical role in explaining a variety of
phenomena in astrophysics[2015Kuijpers], solar physics[2016Liu, 2018Amari,
2018Inoue, 2019Cheung], space physics[2000Webb], and laboratory
plasmas[1986Taylor]. In the solar atmosphere, MFRs are regarded as one of the
basic magnetic structures to host prominences[1998Aulanier, 2017Fan],
sigmoidal structures[2009Green, 2010Aulanier], coronal cavities[2006Gibson,
2013Vourlidas], and hot channels[2011Cheng, 2013Patsourakos]. MFRs can erupt
under some conditions, which can in turn drive prominence eruptions and
coronal mass ejections (CMEs)[2018Green], as well as solar flares when
magnetic reconnection occurs[2011Shibata]. The propagation of CMEs into the
interplanetary space can disturb the space environment and is a potential
hazard to the high technology facilities of the modern society if they direct
at the Earth. Therefore, study on the physical mechanisms of MFR eruptions are
one of the central topics in solar and space physics.

MFR eruptions are usually thought to be driven by ideal magnetohydrodynamic
(MHD) instabilities. For a single isolated MFR, both the helical kink
instability and torus instability could drive its eruption, probably at
different stages. The helical kink instability occurs when the safety factor,
$q$, is lower than a critical value, $q_{\mathrm{c}}$, or when the twist
number measured in turns, $\mathcal{T}$, is greater than a threshold,
$\mathcal{T}_{\mathrm{c}}$, where $q=1/\mathcal{T}$ and
$q_{\mathrm{c}}=1/\mathcal{T}_{\mathrm{c}}$ (see Methods subsection Magnetic
field diagnosis). The classical Kruskal-Shafranov critical value for helical
kink instability is $q_{\mathrm{c}}=1$ or
$\mathcal{T}_{\mathrm{c}}=1$[1954Kruskal]. In fact, the critical values depend
on the details of the MFR equilibrium, such as the radial profile of the
magnetic field, plasma $\beta$, and the flux rope aspect ratio,
$\epsilon=R_{0}/a$, where $R_{0}$ and $a$ are the major and minor radii of the
MFR, respectively. The critical values vary in a range of
$q_{\mathrm{c}}\in[0.6,1.0]$ or $\mathcal{T}_{\mathrm{c}}\in[1.0,1.7]$
depending on the parameter settings of the MFR[1979Hood, 2004Torok]. On the
other hand, the torus instability occurs when the external poloidal magnetic
field decays fast enough along the eruption path of the MFR, namely, the decay
index, $n$, is larger than a critical value, $n_{\mathrm{c}}$ (see Methods
subsection Magnetic field diagnosis). Assuming that the MFR current path is
semicircular and without any internal structures, the stability analysis
yields $n_{\mathrm{c}}=1.5$[2006Kliem]. The critical value for the torus
instability also depends on the details of the MFR equilibrium, including the
shape of the current path. A range of $n_{\mathrm{c}}\in[0.8,1.5]$ is derived
with different prescriptions of the MFR structure[1978Bateman, 2010Olmedo,
2010Demoulin, 2010Fan].

The MFR eruptions can be confined in the lower corona or they fail to
propagate out of the solar corona. A well known mechanism is the failed kink
instability. Namely, the twist of the erupting flux rope is high enough to
incur helical kink instability, but the decay index of the background field
does not exceed the critical value for torus instability along the erupting
path[2005Torok]. Is there a case where the decay index exceeds the critical
value, but the eruption of a flux rope is still confined? Observations have
reported quite a number of such examples[2015Sun, 2019Zhou]. Recently, a
theoretical explanation has been presented in a series of papers[2015Myers,
2017Myers]. The authors pointed out that the previous torus instability model
misses a key component of the Lorentz force, namely, the magnetic tension
force generated by the toroidal magnetic field (also known as the guide field)
and the poloidal electric current in the MFR. The magnetic tension force could
direct downward, similarly to the strapping force from the external poloidal
field perpendicular to the MFR axis. The combination of these two downwardly
directed forces counterparts the upwardly directed hoop force of the MFR. Such
a theory was termed as the failed torus regime[2015Myers], in which a dynamic
magnetic tension force might confine an MFR eruption. A typical signature for
the failed torus instability to work is that the poloidal magnetic field is
rapidly converted to the toroidal field when the MFR erupts, therefore
decreasing the twist number and increasing the safety factor. We note that,
all the previously proposed models assume a certain degree of axisymmetry of
the MFR in which only the poloidal and toroidal components of the magnetic
field are considered while the radial component is neglected. In real
observations, the magnetic topology is very complex and the MFR is usually
non-axisymmetric, so that all the components should be considered when
evaluating the Lorentz force acting on the MFR.

In order to clarify the causes of failed torus instability, MHD simulations
are needed to be performed to reproduce observations. However, there are a
number of difficulties in this problem. First, we deal with a fully three-
dimensional (3D), dynamically evolving plasma system, which rules out the
stationary extrapolation techniques that are usually used. Second, it is very
difficult to extract the key parameters responsible for the failed torus
instability with pure observations, and thus the parameter space is too large
to explore if parameterized simulations are used to fit the observations.
Considering these, a state-of-the-art technique, namely, a data-driven or
data-constrained MHD simulation, should be adopted to explore an observed
failed eruption. Even so, there are still some technical difficulties in
computing the key parameters, such as how to identify the boundary of an MFR
(see Methods subsection Identification of the MFR) and how to compute the
decay index and twist number more accurately.

Here we apply a data-driven MHD model[2019Guo] in the framework of Message
Passing Interface Adaptive Mesh Refinement Versatile Advection Code (MPI-
AMRVAC 2.0[2018Xia]) to a confined MFR eruption[2019Zhong] and diagnose in
detail the physical causes leading to the failure of the eruption with some
elaborately designed diagnostic tools[2017Guo]. We show a Lorentz force
component, induced by the non-axisymmetry (or the radial magnetic field
component) of the MFR, plays a major role in preventing the MFR from erupting
successfully.

## Results

### Overview of the event

The eruption event under study was hosted in NOAA Active Region (AR) 12268,
which appeared on the visible side of the Sun on 21 January 2015. The AR was a
bipolar cluster on the east solar limb. With the emergence and cancellation of
the photospheric magnetic field, a parasitic positive polarity gradually
emerged near the center of the negative polarity in the AR for several days.
Subsequently, the magnetic configuration of the AR became complex. It produced
six Geostationary Operational Environmental Satellite (GOES) M- and C-class
flares[2019Zhong] within 36 hours from January 29 to 30. However, it is
puzzling that none of the six flares produced any CME although the 3 M-class
flares were very intense and erupted in a short interval. Here, we focus on
the intense M2.0 flare (Figure 1) that occurred at 00:32 UT on 30 January
2015, which was observed by the Atmospheric Imaging Assembly[2012lemen] (AIA)
on board Solar Dynamics Observatory[2012pesnell] (SDO). The rationale for our
choice is that a special magnetic structure, namely, a bifurcated MFR, was
formed before the eruption onset[2019Zhong].

Figure 2a shows the AIA 171 Å extreme-ultraviolet image at 00:00 UT, overlaid
with the full view of the reconstructed coronal magnetic fields before the
flare onset. The reconstructed fields are extrapolated using a nonlinear
force-free field model[2016Guo], based on the photospheric vector magnetic
field taken 32 minutes before the flare onset by the Helioseismic and Magnetic
Imager[2012scherrer] (HMI) on board SDO. The extrapolated fields reach a
stable state, which is as close as possible to an equilibrium state, and
satisfy the force-free and divergence-free conditions very well. The magnetic
configuration of the AR is appealing because it is complicated enough. As
shown in Figure 2b, the AR consists of three main polarities. A positive
polarity (P1) in the center is partially enclosed by a negative one (N1
extending to N1e). Another remote positive polarity (P2) also connects to N1
through a bifurcated MFR. We mark the bifurcated and unbranched parts of the
MFR with yellow and orange field lines, respectively. There also exists a
dome-like structure in the corona marked by cyan field lines. Although the
bifurcated MFR has three footpoints, the strong current channel is mainly
concentrated in the central unbranched part of the MFR. This current channel
has two footpoints rooted in the photosphere. We also make the scatter plots
of the twist number of all sample magnetic field lines against the distance
from the MFR axis (Supplementary Figure 1). The maximum value of the absolute
twist number of the MFR is larger than 1.6 turns, which is close to the
maximum threshold for helical kink instability.

Observations at SDO/AIA 94, 171 and 1600 Å passbands (Figure 3) show a clear
process of the confined eruption. In particular, we present the running-
difference between two consecutive AIA 171 Å images in order to reveal more
clearly the eruption structure. The emission at AIA 171 Å originates from the
upper transition region and the quiet corona, which unveils the plasma with a
characteristic temperature of $6.3\times 10^{5}$ K. When the M2.0 flare
occurred, a clear plasmoid-like structure ejected (Figure 3a), and it was
subsequently flowed to the polarity P2 as shown in the 171 Å images. We make
four slices to track the motion of the plasma. As shown in Figure 3b, the
time-distance diagrams for slices 1 and 2 display plasma flows along the
coronal loops towards polarity P2 near the flare peak time, while those for
slices 3 and 4 perpendicular to the coronal loop show that the plasma rises to
a certain height and then stops rising, indicating a failed eruption. Compared
to emission at 171 Å, emission at 94 Å has a higher characteristic temperature
of $6.3\times 10^{6}$ K, which reveals a large current channel in the corona.
The composite snapshot in Figure 3c clearly displays that the emission at 94 Å
is stronger than that at 171 Å in the region of the bifurcated MFR. One can
see that the two contours delineating the strongest 94 Å emission coincide
with the two footpoints of the current channel integrated along the vertical
direction deduced from the extrapolated magnetic field. It is noted that the
flare emission at 1600 Å shows a complex structure with six ribbons (Figure
3d). Ribbons R1, R2 and R5 are semicircular in shape that enclose the polarity
P1. Ribbons R2 and R4, positioned on both sides of the polarity inversion line
(PIL), exhibit a clear shear and separate from each other with flare
development, as predicted in the standard flare model. The intensity of ribbon
R3 is weak, which appears later than the aforementioned ribbons. The emission
characteristics of ribbon R6 at polarity P2 are similar to that of ribbons R2
and R4, suggesting that these ribbons are heated through the same energy
release mechanism.

### 3D MHD simulations

To clarify the mechanism for the failed eruption, we perform a data-driven MHD
simulation using the zero-$\beta$ approximation. The initial condition is
reconstructed from a nonlinear force-free field[2016Guo] (see Methods
subsection Initial and boundary conditions) based on the vector magnetic field
observed by SDO/HMI. The boundary condition on the bottom is provided by the
time series of vector magnetic fields and vector velocities (see Methods
subsection Initial and boundary conditions). Figure 4a–d exhibits the 3D
dynamic evolution of the magnetic field and the electric current density at
four typical moments covering the main period of the event. First, the
bifurcated part of the MFR disappears when it rises initially. Subsequently,
the unbranched part of the MFR continues to rise and further reconnects with
the sheared field lines that surround it. During this process, several field
lines below the MFR become a part of the MFR by magnetic reconnection. Their
footpoints at one end, originally located in the parasitic polarity P1, are
gradually transferred to the remote polarity P2. Meanwhile, the dome-like
structure above is also stretched by the rising MFR. Although the whole MFR
then expands to a very large size, the spatial distributions of four
quantities, including the electric current, quasi separatrix layer, twist
number, and magnetic flux, indicate that the main body of the MFR is still
complete and keeps its coherence. It is noted that a current sheet is
stretched out under the MFR, which has been predicted in the standard solar
flare model. Detailed dynamics of the MFR are shown in Supplementary Movie 1.
We also measure the 3D height of the MFR front and its velocity from the MHD
simulation. As shown in Figure 5a, the MFR goes up at an almost constant
velocity after its initial acceleration, and it then decelerates gradually
after the flare peak, mimicking what has actually happened. We further study
the temporal evolution of the toroidal flux of the MFR calculated by
$\phi_{\mathrm{T}}=\int\mathbf{B}_{\mathrm{T}}\cdot\mathrm{d}\mathbf{S}$,
where $\mathbf{B}_{\mathrm{T}}$ is the total toroidal magnetic field and
$\mathbf{S}$ is the cross-sectional area of the MFR. We then get the poloidal
flux by $\phi_{\mathrm{P}}=\mathcal{T}\phi_{\mathrm{T}}$ (see Methods
subsection Magnetic field diagnosis for computation of the twist number
$\mathcal{T}$). We find that during the eruption of the MFR, the toroidal flux
increases continuously, while the poloidal flux first increases, then
decreases after the flare peak time, and keeps almost unchanged in the late
phase (Figure 5b).

We then characterize the magnetic topology through comparing the simulation
results with the AIA extreme-ultraviolet observations. First, we compare the
3D magnetic field with the AIA emission at 304 Å (Figure 6a). We delineate
typical field lines with different colors, which represent the MFR and the
overlying field. The footpoints of field lines coincide well with the
brightening ribbons of the flare. Second, we calculate the electric current
density and identify some important magnetic topological structures, such as
separators[2010Parnell], hyperbolic flux tube[2002Titov] (HFT), and, more
commonly, quasi-separatrix layers[1995Priest] (QSLs). QSLs are defined as
places with strong gradients of magnetic connectivity. We usually identify the
QSLs through a topological parameter called the squashing factor $Q$. Figure
6b shows a side view of the MFR, which clearly displays the positions of high
electric current density and that of QSLs ($\mathrm{log}(Q)>3$) . It is seen
that a part of high electric current densities coincide well in space with
high $Q$ values. Figure 6c shows the extreme-ultraviolet emission at 94 Å at
the flare peak time, which indicates a hot channel structure along the PIL. We
also make a synthetic map of the AIA emissivity by integrating the electric
current density along the vertical direction at the flare peak time (Figure
6d). It is seen that the emission along the PIL is more prominent at the
position of the MFR.

It is known that QSLs (Figure 7a and 7d) are also favorable locations where
energetic particles can be accelerated and move downward along the field lines
to heat the lower solar atmosphere, resulting in the brightening of flare
ribbons. Previous studies have shown the spatial correspondence between flare
ribbons and photospheric QSLs by using extrapolation models[1997Demoulin],
magneto-frictional models[2015Savcheva, 2016Janvier] or MHD
simulations[2006Aulanier, 2014Dudik]. However, when a flare occurs, the
magnetic topology in the core of the flare region changes rapidly due to
magnetic reconnection. Thus, the magnetic topology itself is very dynamic in
nature. In practice, it is difficult to trace precisely the simulated QSLs
following the observed flare ribbons due to an unavoidable gap between
simulations and observations, such as some unpredictable resistivity in the
real corona but spurious numerical resistivity in simulations[2018Jiang]. In
spite of these difficulties, we try to find the dynamic relationship between
the QSLs and flare ribbons. To this end, we compute the squashing factor $Q$
on the bottom and identify the signed QSLs with $|\mathrm{log}(Q)|>3$, which
are then overlaid on the flare ribbons (Supplementary Figures 2 and 3). Four
typical moments, corresponding to the pre-flare, flare onset, flare peak and
end times, respectively, are shown in order to reveal the dynamic evolution.
We find that the evolution of the QSLs matches the flare ribbons very well.
This means that our simulation reaches a high degree of realism at least in
the dynamic evolution of the magnetic topology. Further quantitative analysis
shows that the separation of the QSLs undergoes two stages. The separation
velocity before the flare peak time is larger than that after the peak time, a
dynamic behaviour which is almost the same as what is revealed by the observed
flare ribbons. More details on the dynamics of the flare ribbons compared with
the QSLs are shown in Supplementary Movie 2.

### Mechanism of the confined eruption

Our model successfully reproduces the process of the failed eruption and can
thus reveal which physical mechanisms finally halt the eruption. In a low
$\beta$ plasma, where the gravity and thermal pressure can be omitted, whether
the MFR eruption is successful or failed depends on the competition between
various Lorentz forces exerting on it, including the hoop force from the
gradient of the magnetic pressure, the tension force from the toroidal
magnetic field, the strapping force from the external poloidal magnetic field
and the force caused by the radial magnetic field of the MFR. If the net
Lorentz force is downward, the MFR cannot break through the overlying magnetic
field to form a CME. To explore this problem quantitatively, we calculate the
Lorentz force along the vertical direction,
$L_{z}=\mathbf{e}_{z}\cdot(\mathbf{J}\times\mathbf{B})$. Figure 7b and 7e show
the spatial distribution of $L_{z}$ on a vertical cut of the simulation domain
at 00:16 and 02:00 UT, respectively. We find that the direction of $L_{z}$
within the radius of the MFR is initially upward, facilitating the eruption,
and it later becomes downward, thus preventing the MFR from further moving up.
The overall distribution of $L_{z}$ at different positions in the MFR is
downward after 00:52 UT. We also notice that the Lorentz force is still upward
but the MFR has already a tiny deceleration around 00:40 UT. A possible reason
is that the Lorentz force acts directly on the local plasma, while the
apparent velocity of the MFR is related not only to the plasma velocity but
also to the change of magnetic structure of the MFR, the latter of which is
changing significantly during this period. The detailed temporal evolution of
$L_{z}$ is shown in Supplementary Movie 3.

We then decompose the total vertical Lorentz force $L_{z}$ into different
components contributed by various combinations of the magnetic field and
electric current components, as listed in Supplementary Table 1 and
Supplementary Figure 4a and 4b. Figure 7c and 7f shows the distribution of the
different force components along a vertical line perpendicular to the MFR axis
at 00:16 and 02:00 UT, respectively. As in previous models, the hoop force
($F_{\mathrm{H}}$) always directs upward below the MFR axis, acting as a
driving force of the eruption. More interestingly, we find that at the time of
02:00 UT, three force components direct downward below the MFR axis, including
the magnetic tension force ($F_{\mathrm{T}}$), and the forces induced by the
radial magnetic field component of the MFR ($F_{\mathrm{N1}}$ and
$F_{\mathrm{N2}}$). The former force does not change its direction as shown in
Supplementary Figure 4c. The latter forces are not included in previous
models. Note that the radial field component deforms the cross section and
thus leads to a non-axisymmetry of the MFR. We thus call these two forces the
non-axisymmetry induced forces. In particular, $F_{\mathrm{N1}}$ changes its
direction during the eruption process. It is initially upward, pushing the MFR
rising, but it becomes downward in the end, thus constraining the MFR
eruption, as shown in Figure 7c, 7f and Supplementary Figure 4d. It is seen
that for this event, such a non-axisymmetry induced force plays the major role
in causing the failed eruption.

Another interesting finding is that the direction of the strapping force is
reversed from downward at 00:16 UT to upward at 02:00 UT as shown in
Supplementary Figure 5. This means that in the later phase, the strapping
force, from the external poloidal field, does not act as the main factor
constraining the eruption, contrary to what is predicted in the traditional
torus instability model. In fact, this active region has a quadrupole field in
large scale, namely, the external poloidal field is mainly contributed by
polarities P2 and N1 at 02:00 UT (long blue arrow in Supplementary Figure 5),
compared to that by P1 and N1e at 00:16 UT (long white arrow in Supplementary
Figure 5). Therefore, the external poloidal field changes its relative
direction with respect to the axis of the MFR when the MFR moves upward.
Consequently, the strapping force,
$F_{\mathrm{S}}=\mathbf{e}_{z}\cdot(\mathbf{J}_{\mathrm{TR}}\times\mathbf{B}_{\mathrm{P,ex}})$,
changes its direction from downward to upward due to the change of the
external poloidal field.

Previous models have mostly invoked helical kink instability and/or torus
instability as the mechanisms triggering and driving the MFR eruption. The key
parameters characterizing these two MHD instabilities are the twist number of
the MFR and the decay index of the external poloidal field. Usually, either
one of the two parameters is not sufficient for judging whether an MFR
successfully erupts or not. An example is that the traditional torus
instability fails to explain the final failure of some eruptions, in which the
twist number decreases because of a conversion of the poloidal magnetic field
to the toroidal field, but the decay index still exceeds the critical value.
In such cases, the increasing tension force makes the eruptions finally
failed[2015Myers]. Therefore, a phase diagram has been proposed by combining
the twist number and decay index[2015Myers] that can significantly improve our
understanding of the eruption behaviour. Here, we also calculate the temporal
evolution of the mean twist of the MFR and the decay index of the external
poloidal magnetic field, which are plotted in the phase diagram of the event
(Figure 8). Based on the critical values for helical kink instability and
torus instability derived from theories[1954Kruskal, 2006Kliem, 1978Bateman],
experiments[2015Myers], observations[2019Zhou, 2018Jing] and MHD
simulations[2010Demoulin, 2005Torok], we divide the parameter space into four
distinct regimes. We can readily find that the MFR undergoes different phases,
starting from the stable regime followed by the failed kink regime, then
moving to the eruptive regime and eventually stopping in the failed torus
regime. In fact, it is the failed torus regime that can explain quite some
failed eruptions that the traditional torus instability cannot[2018Jing]. In
conclusion, our data-driven MHD simulation provides the phase diagram for a
real solar eruption and reveals the failed torus regime for confined
eruptions, which has only been proposed through laboratory
experiments[2015Myers]. More importantly, we propose a scenario for the failed
torus regime, in which a Lorentz force component, resulting from the non-
axisymmetry of the MFR, can play a major role in halting the eruption.

## Discussion

We have performed a data-driven MHD simulation using the measured photospheric
magnetic field to study a realistic confined eruption. The event under
investigation has a complex magnetic configuration including a bifurcated MFR
before the eruption onset. A good match is reached between the simulation
results and observations: the magnetic field and the electric current density
coincide with the extreme-ultraviolet emission both in space and time, and the
dynamics of the magnetic topology is consistent with the evolving flare
ribbons. Our results indicate that the MFR eruption is confined in the lower
corona even if the decay index of the external poloidal field exceeds the
criterion for torus instability.

Therefore, this work reports a confined eruption in the failed torus regime
that is really observed and reproduced with a data-driven MHD model. Such a
failed torus regime has been proposed through laboratory
experiments[2015Myers], but has rarely been found in real solar eruption
events, including statistical studies[2018Jing] and case studies[2015Sun,
2019Zhou]. The reasons are as follows. In the ideal torus instability
regime[1978Bateman], the external field is assumed perpendicular to the axis
of the MFR and the eruption path of the MFR is straightly upward. Thus, the
horizontal component of the external magnetic field, usually from the
potential field extrapolation, is often used to calculate the decay index.
However, in a realistic eruption, the external field is not always
perpendicular to the axis of the MFR, and the eruption path of the MFR is not
always in the radial direction. Therefore, we must consider the relative
orientation between the external magnetic field, the axis of the MFR and the
eruption path rather than simply adopt a horizontal field
approximation[2018Amari, 2018Inoue, 2015Sun, 2019Zhou, 2018Jing, 2017WangD].
On the other hand, previous studies are based heavily on the extrapolation
model to reconstruct the coronal magnetic field for one or limited time
instants, mostly before the eruption onset, which basically requires a quasi-
static approximation. However, in a real situation, the magnetic field
structure may undergo a significant change during the eruption process as
shown in this data-driven MHD simulation. Therefore, at least for some cases,
it seems invalid to use the characteristic parameters (like the twist number
and decay index) at only one or few moments to judge whether an eruption can
occur or not and whether it is successful or failed. Instead, a time-dependent
twist number and decay index may be more appropriate in studying the dynamics
of the MFR evolution. In addition, even in theory, the critical value,
$n_{\mathrm{c}}$, of the decay index depends on the aspect ratio, as well as
the ratio between the height and the footpoint distance of the MFR, which vary
from case to case. Therefore, there is no universal critical value,
$n_{\mathrm{c}}$, that can be applied to observations.

Recently, some studies suggest that the degree of electric current
neutralization may be related to the ability of active regions to produce
CMEs[2017Liu, 2020Avallone]. It provides us another view to explaining the
failed eruption by observational vector magnetograms. We calculate the ratio
of direct current (DC) to return current (RC) in a series of vector
magnetograms for a long time containing 15 moments from 23:48 UT on 29 January
2015 to 02:36 UT on 30 January 2015. As shown in Supplementary Figure 6a, the
active region shows a coherent, middle-scale pattern of DC around its PIL.
Quantitatively, the ratio $|$DC$/$RC$|$ in both positive and negative
polarities only slightly deviates from unity (Supplementary Figure 6b). The
recent statistical study[2020Avallone] implies that a value of $|$DC$/$RC$|$
that is marginally above unity as in our case cannot determine whether the
active region is more likely CME-productive or CME-quiet. Therefore,
prediction of flare/CME production based on the parameter $|$DC$/$RC$|$ alone
is difficult. We further measure the total axial current in the MFR as shown
in Supplementary Figure 7. One can see that the absolute values of both the DC
and RC decrease rapidly while the net current increases significantly before
the flare onset. After the flare peak, the RC is close to zero without obvious
change while both the DC and the net current decrease gradually with time. In
the end, the value of DC is reduced to only half of the original value but
dominates the net current of the MFR. Our event thus displays electric current
non-neutralization much more clearly in the MFR than in the whole active
region.

Our simulation highlights the importance of the ideal MHD instability in
driving the eruption. However, we do not preclude the role of magnetic
reconnection in shaping the behaviour of the eruption. There are mainly three
kinds of reconnection, namely internal reconnection between twisted field
lines within the MFR body, tether-cutting reconnection between sheared field
lines under the MFR, and breakout-type reconnection between the MFR and the
overlying field lines. From Figure 5b, we can judge which kinds of magnetic
reconnection play the main role in different stages. First, before the flare
peak, both toroidal and poloidal fluxes increase rapidly, implying a
consequence of the tether-cutting reconnection. Then, the poloidal flux
decreases while the toroidal flux still increases between 00:40 UT and 01:20
UT. It indicates that in this stage, although the tether-cutting reconnection
still works, the internal reconnection plays a more important role, which
transfers some of the poloidal flux to the toroidal flux. Finally, in the
later stage, the poloidal flux keeps almost unchanged but the toroidal flux
continues increasing. This implies a rough balance between the internal
reconnection and the tether-cutting reconnection in changing the poloidal
flux, since the former decreases but the latter increases the poloidal flux,
respectively. However, both reconnections tend to increase the toroidal flux.
In all the stages, the breakout-type reconnection seems to play a minor role,
since it would decrease the toroidal flux of the MFR, which is always
increasing as shown in Figure 5b. Besides, the absolute value of the mean
twist number, $\mathcal{T}$, of the MFR is shown to decrease in the later two
stages, which is contrary to the scenario of breakout-type reconnection that
tends to increase the value of $\mathcal{T}$. In addition, the Lorentz force
densities are small on the upper boundary of the MFR. This implies that even
if the breakout-type reconnection would be present, it has a tiny influence on
the overall distribution of the Lorentz force.

In our simulation, the field lines writhe and rotate when the MFR rises. The
MFR changes its connectivity in the multipolar magnetic configuration. It
initially connects polarities N1 and P1, but N1 and P2 in the end, with one
footpoint changing from P1 to P2. The MFR writhing could be caused by the
helical kink instability or the Lorentz force exerted on the MFR legs from the
external sheared field[2012Kliem]. Several plausible pieces of evidence
indicate that the internal magnetic reconnection might occur in the twisted
field lines if the MFR writhes and rotates. For instance, we find a 3D bald-
patch separatrix surface configuration[2006bGibson] in Supplementary Figure 8,
showing that two typical field lines meet each other in a region depicted by a
yellow dashed box. This region has a high $Q$ value ($\mathrm{log}(Q)>3$) and
a high electric current density, and thus it might be a possible site for
internal reconnection.

We also find a strong shear motion with a velocity up to 4 km s-1 at the PIL.
When inputting the time series of vector velocities to the bottom boundary to
drive the MHD simulation, such a shear motion may lead to magnetic
reconnection in the lower atmosphere, which produces a mildly imbalanced force
and makes the MFR rise up slowly at the early stage. On the other hand,
although the iteration process in the magneto-frictional method ensures that
the initial condition is as close as possible to an equilibrium state, the
residual force is not strictly zero, which may also incur a very slow motion
of the MFR with the simulation going on.

In conclusion, our data-driven MHD model reproduces well the observations of a
failed solar eruption in both the morphology of the flare and the kinematics
of the MFR. The failure of the eruption can be well interpreted in terms of
the failed torus regime, which cannot be explained by the traditional torus
instability model. In particular, we renovate the physical framework of the
failed torus regime by adding two important ingredients. We find that a force
component, resulting from the radial field of the MFR, can halt the eruption,
in addition to the magnetic tension force previously found[2015Myers]. The
radial field is reflected by the non-axisymmetry of the MFR, which more or
less exists in real events. We also find that the strapping force is not
always constraining but sometimes facilitating the rise of the MFR, in
particular in a quadrupole field. These findings are crucial in more properly
predicting the space weather effect of CMEs that originate in complex and fast
evolving magnetic structures.

## Methods

MHD model. Our MHD model is based on the zero-$\beta$ MHD approximation and
omits the gravity and gradient of gas pressure[2019Guo], which is described as
follows:

$\dfrac{\partial\rho}{\partial t}+\nabla\cdot(\rho\mathbf{v})=0,$ (1)
$\dfrac{\partial(\rho\mathbf{v})}{\partial
t}+\nabla\cdot\bigg{(}\rho\mathbf{v}\mathbf{v}+\frac{1}{2}B^{2}\mathbf{I}-\mathbf{B}\mathbf{B}\bigg{)}=0,$
(2) $\dfrac{\partial\mathbf{B}}{\partial
t}+\nabla\cdot(\mathbf{v}\mathbf{B}-\mathbf{B}\mathbf{v})=0,$ (3)
$\mathbf{J}=\nabla\times\mathbf{B},$ (4)

where $\mathbf{\rho}$ represents the density, $\mathbf{v}$ the velocity,
$\mathbf{B}$ the magnetic field, and $\mathbf{I}$ the unit tensor. Our model
also omits density diffusion in the mass conservation equation, dynamic
viscosity in the momentum conservation equation and resistivity in the
magnetic induction equation. Although there is no explicit resistivity in the
magnetic induction equation, magnetic reconnection still occurs due to
existence of numerical resistivity. Equations (1)–(4) are written in the
dimensionless form. The length, density, magnetic field, velocity and time are
normalized by the characteristic quantities $L_{0}=1.0\times 10^{9}$ cm,
$\rho_{0}=2.3\times 10^{-15}\ \mathrm{g\ cm^{-3}}$, $B_{0}=2.0$ G, $v_{0}$ =
$B_{0}/(\mu_{0}\rho_{0})^{1/2}$ and $t_{0}=L_{0}/v_{0}$, respectively. The
vacuum permeability, $\mu_{0}$, is assumed to be unity. We solve the MHD
equations numerically using the open source MPI-AMRVAC 2.0[2018Xia] and employ
Powell’s method[1999powell] to clean the divergence of the magnetic field. Our
computation domain has a size of $412\times 286\times 286$ $\mathrm{Mm^{3}}$,
which is divided uniformly by $288\times 200\times 200$ cells. The full
computation box, which corresponds to $X$ = [$-$3′′, 573′′] and $Y$ =
[$-$295′′, 105′′] in the $X$–$Y$ plane, is larger than the field of view shown
in Figure 2b and Supplementary Movie 1. Note that the fast mode wave produced
by the eruption has not arrived at the top and side boundaries when the
eruption is being confined, so that the boundary conditions can hardly affect
the dynamics of the MFR.

The temperature distribution of the initial atmosphere is simplified as a
stepwise function:

$T=\begin{cases}T_{0}&0.0\leq h<h_{0}\\\ k_{T}(h-h_{0})+T_{0}&h_{0}\leq
h<h_{1}\\\ T_{1}&h_{1}\leq h<28.6\end{cases},$ (5)

where $T_{0}=1.4\times 10^{4}$ K, $T_{1}=1.0\times 10^{6}$ K, $h_{0}=3.5\times
10^{8}$ cm, $h_{1}=1.0\times 10^{9}$ cm and
$k_{T}=(T_{1}-T_{0})/(h_{1}-h_{0})$ in cgs units. Such a distribution
approximately mimics the thermal structure of the actual atmosphere. The
initial density profile is obtained by assuming hydrostatic equilibrium:

$\frac{\mathrm{d}p}{\mathrm{d}h}=-\rho g,$ (6)

where the gas pressure is expressed as $p=\rho T$ dimensionlessly. The
density, $\rho$, is chosen to be $1.0\times 10^{8}$ (namely, $2.3\times
10^{-7}$ g cm-3 before normalization) on the bottom boundary, and it decreases
with height until it is fixed to be $5.0\times 10^{3}$ above $h=7.5\times
10^{9}$ cm.

Initial and boundary conditions. The initial condition for the MHD model is
provided by the nonlinear force-free field (NLFFF) derived from the magneto-
frictional method[2016Guo]. On the bottom boundary, we use the vector magnetic
field data at 00:00 UT (about 32 minutes before the onset time of the M2.0
flare) on 30 January 2015, which were observed by SDO/HMI and provided by the
Joint Science Operations Center. The series name for the magnetic field data
is hmi.ME$\\_$720s$\\_$fd10. Some preprocessing procedures need to be done
before the magnetic field data can be used in the model, which include removal
of the 180∘ ambiguity[2006Metcalf] in the transverse field, correction of the
projection effect, and removal of the residual Lorentz force and torque in the
original data[2006Wiegelmann]. We project the computation box to the
heliocentric coordinate system in Figure 2a and plot the reconstructed
magnetic field for the initial atmosphere in Figure 2b. The figure provides a
comparison between the reconstructed field and SDO/AIA observations.

After obtaining a self-consistent initial model, we then run the MHD
simulation with the bottom boundary driven by the time series of the vector
magnetic fields and the velocity fields, which contains 11 moments from 00:00
UT to 02:00 UT with an interval of 12 minutes. The magnetic fields should be
preprocessed as above and the velocity fields can be calculated by using the
Differential Affine Velocity Estimator for Vector Magnetograms[2008Schuck]. We
directly set the magnetic field data to the inner ghost layer and make a zero-
gradient extrapolation to the outer ghost layer of the bottom boundary. Since
the vector magnetograms have been acquired with a cadence of 12 minutes, which
is not sufficient for simulations. We then generate a new series of
magnetograms with a higher cadence by linear interpolation in time. The
magnetic field on the top and all side boundaries is given by the zero-
gradient extrapolation. Similarly, we set the velocity data to the inner ghost
layer and make a zero-gradient extrapolation to the outer ghost layer of the
bottom boundary. In addition, the velocity is set to be zero on the top and
four side boundaries. The density is assumed not to change on all the
boundaries.

Identification of the MFR. An MFR consists of a bundle of helical field lines
with electric currents flowing inside. First, we calculate the twist number of
all sample field lines as a function of the distance, $r$, from a pre-assumed
MFR axis. An MFR is confirmed, if there is a group of field lines with the
twist number larger than 1 turn. Then, we identify the boundary of the MFR by
calculating the 3D QSLs. For example, the yellow semi-transparent isosurfaces
in Supplementary Figure 9a delineate the 3D boundary of the MFR at 00:00 UT.
We choose a high $Q$ value with $\mathrm{log}(Q)>3$ as the boundary of the
MFR. In a 2D view shown in Supplementary Figure 9b, the MFR boundary is more
clearly delineated by QSLs, together with some selected field lines that
represent the main body of the MFR.

Magnetic field diagnosis. The magnetic twist[2017Guo] of an MFR is defined as

$\mathcal{T}=\frac{1}{2\pi}\int_{s}\mathbf{T}(s)\cdot\mathbf{V}(s)\times\frac{\mathrm{d}\mathbf{V}(s)}{\mathrm{d}s}\mathrm{d}s,$
(7)

where $\mathbf{T}(s)$ denotes the unit vector tangent to the axis of the MFR
and $s$ is the arc length from a reference starting point on the axis curve.
$\mathbf{V}(s)$ is a unit vector normal to $\mathbf{T}(s)$ and pointing from
the axis curve to the secondary curve. Obviously, the twist number
$\mathcal{T}$ can be computed as a function of the distance to the axis curve,
which in practice can be chosen to be an arbitrary field line within the MFR.
Here, we define the axis curve as the field line to which the absolute value
of the average twist number is the largest. We use this largest twist number
as a proxy for judging if it meets the critical condition for helical kink
instability.

The decay index of the external poloidal magnetic field is defined as

$n=-\frac{\mathrm{d}(\log\left(B_{\mathrm{P,ex}}\right))}{\mathrm{d}(\log(r))},$
(8)

where $B_{\mathrm{P,ex}}$ is the magnitude of the poloidal component of the
external field perpendicular to the axis and the eruption path of the MFR, and
$r$ is the distance measured along the eruption path of the MFR. Here, the
poloidal component of the potential field at the apex of the MFR is adopted to
represent the value of $B_{\mathrm{P,ex}}$ for calculating the decay index. We
use this decay index as a proxy in the phase diagram to see if it reaches the
threshold for torus instability.

Data availability. The SDO/AIA multiwavelength data and SDO/HMI vector
magnetograms can be downloaded from
http://jsoc.stanford.edu/ajax/exportdata.html. The GOES X-ray flux data can be
downloaded from https://www.ngdc.noaa.gov/stp/satellite/goes/dataaccess.html.
All snapshots of MHD simulation data used in the present study are available
on request to the corresponding author. The data for all figures are avaliable
at the website https://doi.org/10.6084/m9.figshare.12949349.v2.

Code availability. The MHD simulation code we used is the open source code of
MPI-AMRVAC 2.0, which can be download from its website
https://github.com/amrvac/amrvac/releases/tag/v2.0/. The diagnostic tools of
solar magnetic fields in the above analyses are available at the website
https://github.com/njuguoyang/magnetic_modeling_codes.

## References

  * 1. Kuijpers, J., Frey, H. U. & Fletcher, L.  Electric current circuits in astrophysics.  _Space Sci. Rev._ 188, 3–57 (2015).
  * 2. Liu, R. _et al._ Structure, stability, and evolution of magnetic flux ropes from the perspective of magnetic twist.  _Astrophys. J._ 818, 148 (2016).
  * 3. Amari, T., Canou, A., Aly, J.-J., Delyon, F. & Alauzet, F.  Magnetic cage and rope as the key for solar eruptions.  _Nature_ 554, 211–215 (2018).
  * 4. Inoue, S., Kusano, K., Büchner, J. & Skála, J.  Formation and dynamics of a solar eruptive flux tube.  _Nat. Commun._ 9, 174 (2018).
  * 5. Cheung, M. C. M. _et al._ A comprehensive three-dimensional radiative magnetohydrodynamic simulation of a solar flare.  _Nat. Astron._ 3, 160–166 (2019).
  * 6. Webb, D. F., Cliver, E. W., Crooker, N. U., Cry, O. C. S. & Thompson, B. J.  Relationship of halo coronal mass ejections, magnetic clouds, and magnetic storms.  _J. Geophys. Res._ 105, 7491–7508 (2000).
  * 7. Taylor, J. B.  Relaxation and magnetic reconnection in plasmas.  _Rev. Mod. Phys._ 58, 741–763 (1986).
  * 8. Aulanier, G. & Démoulin, P.  3-D magnetic configurations supporting prominences. I. The natural presence of lateral feet.  _Astron. Astrophys._ 329, 1125–1137 (1998).
  * 9. Fan, Y.  MHD simulations of the eruption of coronal flux ropes under coronal streamers.  _Astrophys. J._ 844, 26 (2017).
  * 10. Green, L. M. & Kliem, B.  Flux rope formation preceding coronal mass ejection onset.  _Astrophys. J._ 700, L83–L87 (2009).
  * 11. Aulanier, G., Török, T., Démoulin, P. & DeLuca, E. E.  Formation of torus-unstable flux ropes and electric currents in erupting sigmoids.  _Astrophys. J._ 708, 314–333 (2010).
  * 12. Gibson, S. E., Foster, D., Burkepile, J., de Toma, G. & Stanger, A.  The calm before the storm: The link between quiescent cavities and coronal mass ejections.  _Astrophys. J._ 641, 590–605 (2006).
  * 13. Vourlidas, A., Lynch, B. J., Howard, R. A. & Li, Y.  How many CMEs have flux ropes? Deciphering the signatures of shocks, flux ropes, and prominences in coronagraph observations of CMEs.  _Sol. Phys._ 284, 179–201 (2013).
  * 14. Cheng, X., Zhang, J., Liu, Y. & Ding, M. D.  Observing flux rope formation during the impulsive phase of a solar eruption.  _Astrophys. J. Lett._ 732, L25 (2011).
  * 15. Patsourakos, S., Vourlidas, A. & Stenborg, G.  Direct evidence for a fast coronal mass ejection driven by the prior formation and subsequent destabilization of a magnetic flux rope.  _Astrophys. J._ 764, 125 (2013).
  * 16. Green, L. M., Török, T., Vršnak, B., Manchester, W. & Veronig, A.  The origin, early evolution and predictability of solar eruptions.  _Space Sci. Rev._ 214, 46 (2018).
  * 17. Shibata, K. & Magara, T.  Solar flares: Magnetohydrodynamic processes.  _Living Rev. Sol. Phys._ 8, 6 (2011).
  * 18. Kruskal, M. & Schwarzschild, M.  Some instabilities of a completely ionized plasma.  _Proc. R. Soc. Lond. A_ 223, 348–360 (1954).
  * 19. Hood, A. W. & Priest, E. R.  Kink instability of solar coronal loops as the cause of solar flares.  _Sol. Phys._ 64, 303–321 (1979).
  * 20. Török, T., Kliem, B. & Titov, V. S.  Ideal kink instability of a magnetic loop equilibrium.  _Astron. Astrophys._ 413, L27–L30 (2004).
  * 21. Kliem, B. & Török, T.  Torus instability.  _Phys. Rev. Lett._ 96, 255002 (2006).
  * 22. Bateman, G.  _MHD instabilities_ (Cambridge, MA: MIT Press, 1978).
  * 23. Olmedo, O. & Zhang, J.  Partial torus instability.  _Astrophys. J._ 718, 433–440 (2010).
  * 24. Démoulin, P. & Aulanier, G.  Criteria for flux rope eruption: Non-equilibrium versus torus instability.  _Astrophys. J._ 718, 1388–1399 (2010).
  * 25. Fan, Y.  On the eruption of coronal flux ropes.  _Astrophys. J._ 719, 728–736 (2010).
  * 26. Török, T. & Kliem, B.  Confined and ejective eruptions of kink-unstable flux ropes.  _Astrophys. J. Lett._ 630, L97–L100 (2005).
  * 27. Sun, X. _et al._ Why is the great solar active region 12192 flare-rich but CME-poor?  _Astrophys. J. Lett._ 804, L28 (2015).
  * 28. Zhou, Z. _et al._ Why do torus-unstable solar filaments experience failed eruptions?  _Astrophys. J. Lett._ 877, L28 (2019).
  * 29. Myers, C. E. _et al._ A dynamic magnetic tension force as the cause of failed solar eruptions.  _Nature_ 528, 526–529 (2015).
  * 30. Myers, C. E. _et al._ Quasi-static and dynamic magnetic tension forces in arched, line-tied magnetic flux ropes.  _Plasma Phys. Control. Fusion_ 59, 014048 (2017).
  * 31. Guo, Y., Xia, C., Keppens, R., Ding, M. D. & Chen, P. F.  Solar magnetic flux rope eruption simulated by a data-driven magnetohydrodynamic model.  _Astrophys. J. Lett._ 870, L21 (2019).
  * 32. Xia, C., Teunissen, J., Mellah, I. E., Chané, E. & Keppens, R.  MPI-AMRVAC 2.0 for solar and astrophysical applications.  _Astrophys. J. Suppl._ 234, 30 (2018).
  * 33. Zhong, Z., Guo, Y., Ding, M. D., Fang, C. & Hao, Q.  Transition from circular-ribbon to parallel-ribbon flares associated with a bifurcated magnetic flux rope.  _Astrophys. J._ 871, 105 (2019).
  * 34. Guo, Y. _et al._ Magnetic helicity estimations in models and observations of the solar magnetic field. III. Twist number method.  _Astrophys. J._ 840, 40 (2017).
  * 35. Lemen, J. R. _et al._ The atmospheric imaging assembly (AIA) on the solar dynamics observatory (SDO).  _Sol. Phys._ 275, 17–40 (2012).
  * 36. Pesnell, W. D., Thompson, B. J. & Chamberlin, P. C.  The solar dynamics observatory (SDO).  _Sol. Phys._ 275, 3–15 (2012).
  * 37. Guo, Y., Xia, C., Keppens, R. & Valori, G.  Magneto-frictional modeling of coronal nonlinear force-free fields. I. Testing with analytic solutions.  _Astrophys. J._ 828, 82 (2016).
  * 38. Scherrer, P. H. _et al._ The helioseismic and magnetic imager (HMI) investigation for the solar dynamics observatory (SDO).  _Sol. Phys._ 275, 207–227 (2012).
  * 39. Parnell, C. E., Maclean, R. C. & Haynes, A. L.  The detection of numerous magnetic separators in a three-dimensional magnetohydrodynamic model of solar emerging flux.  _Astrophys. J. Lett._ 725, L214–L218 (2010).
  * 40. Titov, V. S., Hornig, G. & Démoulin, P.  Theory of magnetic connectivity in the solar corona.  _J. Geophys. Res._ 107, 1164 (2002).
  * 41. Priest, E. R. & Démoulin, P.  Three-dimensional magnetic reconnection without null points. 1. Basic theory of magnetic flipping.  _J. Geophys. Res._ 100, 23443–23463 (1995).
  * 42. Démoulin, P., Bagala, L. G., Mandrini, C. H., Henoux, J. C. & Rovira, M. G.  Quasi-separatrix layers in solar flares. II. Observed magnetic configurations.  _Astron. Astrophys._ 325, 305–317 (1997).
  * 43. Savcheva, A. _et al._ The relation between solar eruption topologies and observed flare and observed flare features. I. Flare ribbons.  _Astrophys. J._ 810, 96 (2015).
  * 44. Janvier, M. _et al._ Evolution of flare ribbons, electric currents, and quasi-separatrix layers during an X-class flare.  _Astron. Astrophys._ 591, A141 (2016).
  * 45. Aulanier, G., Pariat, E., Démoulin, P. & Devore, C. R.  Slip-running reconnection in quasi-separatrix layers.  _Sol. Phys._ 238, 347–376 (2006).
  * 46. Dudík, J. _et al._ Slipping magnetic reconnection during an X-class solar flare observed by SDO/AIA.  _Astrophys. J._ 784, 144 (2014).
  * 47. Jiang, C. _et al._ Magnetohydrodynamic simulation of the X9.3 flare on 2017 September 6: Evolving magnetic topology.  _Astrophys. J._ 869, 13 (2018).
  * 48. Jing, J. _et al._ Statistical analysis of torus and kink instabilities in solar eruptions.  _Astrophys. J._ 864, 138 (2018).
  * 49. Wang, D. _et al._ Critical height of the torus instability in two-ribbon solar flares.  _Astrophys. J._ 843, L9 (2017).
  * 50. Liu, Y., Sun, X., Török, T., Titov, V. S. & Leake, J. E.  Electric-current neutralization, magnetic shear, and eruptive activity in solar active regions.  _Astrophys. J. Lett._ 846, L6 (2017).
  * 51. Avallone, E. A. & Sun, X.  Electric current neutralization in solar active regions and its relation to eruptive activity.  _Astrophys. J._ 893, 123 (2020).
  * 52. Kliem, B., Török, T. & Thompson, W. T.  A parametric study of erupting flux rope rotation. Modeling the “cartwheel CME” on 9 April 2008.  _Sol. Phys._ 281, 137–166 (2012).
  * 53. Gibson, S. E. & Fan, Y.  The partial expulsion of a magnetic flux rope.  _Astrophys. J. Lett._ 637, L65–L68 (2006).
  * 54. Powell, K. G., Roe, P. L., Linde, T. J., Gombosi, T. I. & De Zeeuw, D. L.  A solution-adaptive upwind scheme for ideal magnetohydrodynamics.  _J. Comput. Phys._ 154, 284 – 309 (1999).
  * 55. Metcalf, T. R. _et al._ An overview of existing algorithms for resolving the 180∘ ambiguity in vector magnetic fields: Quantitative tests with synthetic data.  _Sol. Phys._ 237, 267–296 (2006).
  * 56. Wiegelmann, T., Inhester, B. & Sakurai, T.  Preprocessing of vector magnetograph data for a nonlinear force-free magnetic field reconstruction.  _Sol. Phys._ 233, 215–232 (2006).
  * 57. Schuck, P. W.  Tracking vector magnetograms with the magnetic induction equation.  _Astrophys. J._ 683, 1134–1152 (2008).

## Acknowledgements

We thank the Solar Dynamics Observatory (SDO) team for providing the data for
study. SDO is a mission of NASA’s Living With a Star Program. Z.Z.
acknowledges Kai Yang for his help. Z.Z., Y.G., and M.D.D. are supported by
NSFC under grants 11733003, 11773016, 11533005, and 11961131002, and the China
Scholarship Council 201906190107. The MHD simulations were conducted on the
cluster system of the High Performance Computing Center (HPCC) in Nanjing
University.

## Author contributions

Z.Z. analysed the observational data and performed model calculations. Y.G.
contributed to the theoretical formulation of the model and supervised the
project. M.D.D. conceived the study and supervised the project. All authors
discussed the results and wrote the manuscript.

## Competing interests

The authors declare no competing interests.

Figure 1:  GOES soft X-ray fluxes of the M2.0 class solar flare on 30 January
2015 at 1.0–8.0 Å (red curve) and 0.5–4.0 Å (blue curve). The purple vertical
line shows the peak time of the flare at 00:44 UT. The three letters of M, C
and B represent the flare classes, according to the X-ray peak flux at 1.0–8.0
Å.  Figure 2:  Overview of the M2.0 class solar flare observed on 30 January
2015. a Full-disk image of AIA 171 Å at 00:00 UT on January 30. Some selected
field lines (in green) at a relatively high altitude are overlaid in the
active region. b The initial reconstructed magnetic field with a zoomed-in
view of the region corresponding to the white dotted box in panel a. Yellow
and orange lines represent the bifurcated and unbranched parts of the MFR,
respectively. Cyan lines represent the dome-like structure. The white semi-
transparent isosurface represents the electric current density that is larger
than $32.9\%$ of the maximum value in the whole domain. The transparent slice
perpendicular to the MFR axis shows a cross section of the MFR delineated by
the QSLs. The background image shows the magnetogram with polarities labeled
as N1, N1e, P1 and P2, overlaid by the white and black contours with contour
levels of $B_{z}$ ($B_{z}\equiv\mathbf{e}_{z}\cdot\mathbf{B}$) being $-$50 G
and 50 G, respectively.
Figure 3:  Flare emission in AIA multiwavelength channels. a Running
difference image of AIA 171 Å showing the eruption structure at the peak time
of the M2.0 flare, 00:44 UT, on 30 January 2015. The red and blue colors
represent positive and negative differences in intensity, respectively. The
black contour in the background shows the positive polarity with $B_{z}$ being
+50 G. The cyan and orange lines show two different slices parallel to the
coronal loop, while the green and black lines refer to the slices
perpendicular to the coronal loop. b Time-distance diagrams showing the motion
of the plasma along the four slices either parallel or perpendicular to the
coronal loop. The black dotted line shows the flare peak time at 00:44 UT. c
Composite image of AIA 94 Å (red) and 171 Å (cyan) at the onset time of the
M2.0 flare, 00:32 UT, on 30 January 2015. The green and white contours
represent $60\%$ of the maximum emission at 94 Å and $50\%$ of the maximum
electric current integrated along the line of sight, respectively. d Time-
series of AIA 1600 Å emission starting from 00:32 UT. The black contour is the
same as that in panel a. The arrows point to the flare ribbons, which are
labeled as R1–R6.
Figure 4:  Snapshots showing the MFR and the electric current density. a–d
Four snapshots at 00:00, 00:32, 00:44 and 01:20 UT on 30 January 2015
corresponding to the pre-flare, flare onset, flare peak and end times,
respectively. The yellow and orange lines have the same meaning as that in
Figure 2b. The pink and olive lines within the MFR represent the field lines
connecting the negative polarity N1 and the positive polarities P1 and P2,
respectively. The vertical transparent slice displays the electric current
density. The white semi-transparent isosurface represents the electric current
density larger than $32.9\%$ of the maximum value in the whole domain. The
background shows the distribution of the vertical magnetic field component,
$B_{z}$.
Figure 5:  Kinematics and magnetic fluxes of the MFR. a Time-height
measurement of the MFR in the MHD simulation. The orange triangles and purple
circles show the height and velocity, respectively. b Temporal evolutions of
the toroidal flux (blue squares) and poloidal flux (red circles) of the MFR.
The green dashed line represents the peak time of the M2.0 flare at 00:44 UT.
Figure 6:  Comparison between MHD simulations and AIA extreme-ultraviolet
observations at the peak time of the M2.0 flare. a Selected magnetic field
lines overlaid on an AIA 304 Å image from the vertical perspective. The field
lines in green are located at a relatively high altitude. The pink and olive
lines have the same meaning as that in Figure 4. The cyan, yellow and orange
lines have the same meaning as that in Figure 2b. b The MFR from a side view.
The field lines in pink are the same as that in panel a. The vertical slice
perpendicular to the MFR axis displays the electric current density overlaid
by the contour of the QSLs with $\mathrm{log}(Q)>3$. The two red arrows
indicate the inflow at the reconnection site. The white semi-transparent
isosurface represents the electric current density larger than $32.9\%$ of the
maximum value in the whole domain. The bottom surface shows the distribution
of the vertical magnetic field component, $B_{z}$, overlaid by the polarity
inversion line in cyan. c AIA 94 Å image at 00:44 UT on 30 January 2015. d
Synthetic emission at 94 Å calculated by the integration of the electric
current density along the vertical direction at the peak time.
Figure 7:  Comparison between the QSLs and the Lorentz force, $L_{z}$. The
upper panels are for the time of 00:16 UT and the lower panels for the time of
02:00 UT, which represent the early phase and the later phase of the eruption,
respectively. a Distribution of the QSLs in the $Y$–$Z$ plane at the time of
00:16 UT. The QSLs show the boundary of the MFR. The white plus symbol refers
to the axis of the flux rope. The background at the bottom of the figure shows
the distribution of the vertical magnetic field component, $B_{z}$. The yellow
and orange semi-transparent isosurfaces represent the electric current density
larger than $21.9\%$ and $32.9\%$ of the maximum value in the whole domain,
respectively. b Distribution of the Lorentz force in the $Y$–$Z$ plane at the
time of 00:16 UT. The white plus symbol has the same meaning as that in a. The
white vertical line has a length of 31.5 Mm (covering a height range from 16.8
to 48.3 Mm), which extends from the bottom to the top of the MFR. The
background is the same as that in a. c Distribution of the vertical component
of the Lorentz force, $L_{z}$, along the white line in panel b. The gray line
is for the net force. The orange, purple, blue, red and cyan lines represent
the different components including the hoop force ($F_{\mathrm{H}}$), non-
axisymmetry induced forces ($F_{\mathrm{N1}}$ and $F_{\mathrm{N2}}$),
strapping force ($F_{\mathrm{S}}$) and tension force ($F_{\mathrm{T}}$), which
are listed in Supplementary Table 1. The black vertical line shows the
location of the axis of the MFR. d Similar to panel a but at the time of 02:00
UT. e Similar to panel b but at the time of 02:00 UT. The white vertical line
has a length of 63.0 Mm (covering a height range from 65.8 to 128.8 Mm). f
Similar to panel c but at the time of 02:00 UT.
Figure 8:  Phase diagram showing the temporal evolution of the absolute value
of the mean twist of the MFR versus the decay index of the external poloidal
magnetic field. The $X$-axis represents the twist number $\mathcal{T}$
measured in turns, while the $Y$-axis represents the decay index, which are
considered as basic parameters for helical kink instability and torus
instability, respectively. The whole region is divided into four distinct
parameter quadrants by the critical values for two instabilities. The shaded
region refers to the range (the lower and upper limits) of the critical value
for torus instability, as derived from theories[2006Kliem, 1978Bateman,
2010Olmedo, 2010Demoulin], observations[2015Sun, 2019Zhou, 2018Jing] or
experiments[2015Myers, 2017Myers]. The error bar shown in the bottom left part
of the panel refers to the uncertainty of the mean twist for a typical moment.

Supplementary Figures

Supplementary Figure 1:  Magnetic properties of the initial flux rope. a The
distribution of the twist number as a function of the distance to the MFR axis
calculated from the initial reconstructed magnetic field. The strength of
$B_{z}$ measured in the positive polarity at the footpoints of the field lines
is marked by different colors. b The position of the MFR axis labeled as a
white plus symbol. Other field lines are labeled as red plus symbols, which
are surrounded by the boundary of the MFR delineated by the QSLs shown on the
transparent slice in Figure 2b.
Supplementary Figure 2:  Comparison between the QSLs on the bottom and the
extreme-ultraviolet emission in AIA 304 Å. a–d Time series of AIA 304 Å images
at 00:00, 00:32, 00:44 and 01:20 UT on 30 January 2015, corresponding to the
pre-flare, flare onset, flare peak and end times, respectively. The line of
white plus symbols in panel a shows a slice crossing the flare ribbons. The
QSLs with $|\mathrm{log}(Q)|>3$, signed with positive and negative polarities,
are overlaid on the 304 Å images. e Quantitative comparison between the signed
QSLs with $|\mathrm{log}(Q)|>3$ and the flare ribbons. The black and purple
dotted lines show the separation motion of the flare ribbon in the positive
polarity, with velocities being 1.9 km s-1 and 0.5 km s-1, respectively. The
white and yellow dotted lines show the motion of the flare ribbon in the
negative polarity, with velocities being 2.99 km s-1 and 0.84 km s-1,
respectively. The cyan dotted line indicates the flare peak time at 00:44 UT.
Supplementary Figure 3:  Comparison between the QSLs on the bottom and the
extreme-ultraviolet emission in AIA 1600 Å. a–d Time series of AIA 1600 Å
images at 00:00, 00:32, 00:44 and 01:20 UT on 30 January 2015, corresponding
to the pre-flare, flare onset, flare peak and end times, respectively. The
line of white plus symbols in panel a shows a slice crossing the flare
ribbons. The QSLs with $|\mathrm{log}(Q)|>3$, signed with positive and
negative polarities, are overlaid on the 1600 Å images. e Quantitative
comparison between the signed QSLs with $|\mathrm{log}(Q)|>3$ and the flare
ribbons. The black and purple dotted lines show the separation motion of the
flare ribbon in the positive polarity, with velocities being 1.9 km s-1 and
0.5 km s-1, respectively. The white and yellow dotted lines show the motion of
the flare ribbon in the negative polarity, with velocities being 2.99 km s-1
and 0.84 km s-1, respectively. The cyan dotted line indicates the flare peak
time at 00:44 UT.
Supplementary Figure 4:  Schematic picture of magnetic fields, electric
currents and Lorentz forces in the MFR. a The MFR depicted by the red semi-
transparent surface. The magnetic field components ($\mathbf{B}_{\mathrm{R}}$,
$\mathbf{B}_{\mathrm{P}}$, $\mathbf{B}_{\mathrm{T}}$) are shown in cyan while
the electric current components ($\mathbf{J}_{\mathrm{PT}}$,
$\mathbf{J}_{\mathrm{TR}}$, $\mathbf{J}_{\mathrm{RP}}$) are shown in yellow.
The background image in black and white displays the SDO/HMI $B_{z}$
component. b The local coordinate system used to define the three components.
c Sketch map of the hoop and tension forces. d Sketch map of the strapping
force and non-axisymmetry induced forces 1 and 2 at two moments before and
after the flare onset, respectively.
Supplementary Figure 5:  Sketch map of the strapping force. The yellow, orange
and cyan lines represent the main body of the MFR. The white semi-transparent
contour represents the electric current density that is larger than $17\%$ of
the maximum value in the whole domain. The direction of the electric current
in the MFR points inwards. The long white and blue arrows represent the
direction of the initial and final external poloidal fields, respectively.
Also drawn at the bottom are the Cartesian coordinate system and two sketch
maps of the strapping force corresponding to the initial and final fields. The
fields associated with a downward force are shown in white, while the fields
associated with an upward force are shown in blue. The background image shows
the magnetogram with polarities labeled as N1, N1e, P1 and P2, overlaid by the
white and black contours with contour levels of $B_{z}$ being $-$50 G and 50
G, respectively.
Supplementary Figure 6:  Temporal evolution of the vertical electric current
density $J_{z}$. a The vertical electric current density $J_{z}$ at 23:48 UT
on 29 January 2015, which is derived from SDO/HMI vector magnetograms,
overlaid by the black and green contours with contour levels of $B_{z}$ being
$-$300 G and 300 G, respectively. b The evolution of the unsigned magnetic
flux (gray line) and the ratio of direct current to return current,
$|$DC$/$RC$|$, in both positive (red line) and negative (purple line)
polarities. The green dashed line represents the peak time of the M2.0 flare
at 00:44 UT.
Supplementary Figure 7:  Temporal evolution of the axial electric currents of
the MFR. The red, orange and purple lines represent the positive, negative and
net electric currents, respectively. The green vertical line shows the peak
time of the M2.0 flare at 00:44 UT.
Supplementary Figure 8:  Selected magnetic field lines and electric current
density at 01:20 UT showing a possible internal reconnection. a A side view
along the $Y$-axis. The vertical transparent slice displays the electric
current density depicted by a high $Q$ value of $\mathrm{log}(Q)>3$ (in pink).
The red and green lines represent a configuration similar to the bald-patch
separatrix surface. The orange lines represent a group of lines close to the
red line, while the olive lines are a group close to the green line. The white
isosurface represents the electric current density larger than $32.9\%$ of the
maximum value in the whole domain. The background shows the distribution of
the vertical magnetic field component, $B_{z}$. b A side view along the green
field line at its apex. The yellow dashed box depicts a region where the green
and red field lines meet each other.
Supplementary Figure 9:  Overview of the MFR structure shown with magnetic
field lines and the $Q$ value. The 3D semi-transparent yellow boundary
represents a high $Q$ value of $\mathrm{log}(Q)>3$. a Selected field lines of
the MFR. The yellow and orange lines represent the bifurcated and unbifurcated
parts of the MFR, respectively. The green lines show the main body of the MFR.
The background image in black and white displays the SDO/HMI $B_{z}$
component. b An 2D slice displaying the boundary of the MFR depicted by a high
$Q$ value (in red). The slice and the field lines are the same as that in
panel a but with a different viewing angle.

Supplementary Table

Supplementary Table 1: Different components of the magnetic field, electric
current and Lorentz force

Quantity |  |  |  |  |  | Symbola,b |  |  |  |  | Expression
---|---|---|---|---|---|---|---|---|---|---|---
External radial magnetic field |  |  |  |  |  | $\mathbf{B}_{\mathrm{R,ex}}$ |  |  |  |  | …
Rope radial magnetic field |  |  |  |  |  | $\mathbf{B}_{\mathrm{R,in}}$ |  |  |  |  | …
External poloidal magnetic field |  |  |  |  |  | $\mathbf{B}_{\mathrm{P,ex}}$ |  |  |  |  | …
Rope poloidal magnetic field |  |  |  |  |  | $\mathbf{B}_{\mathrm{P,in}}$ |  |  |  |  | …
External toroidal magnetic field |  |  |  |  |  | $\mathbf{B}_{\mathrm{T,ex}}$ |  |  |  |  | …
Rope toroidal magnetic field |  |  |  |  |  | $\mathbf{B}_{\mathrm{T,in}}$ |  |  |  |  | …
Total radial magnetic field |  |  |  |  |  | $\mathbf{B}_{\mathrm{R}}$ |  |  |  |  | $\mathbf{B}_{\mathrm{R,ex}}+\mathbf{B}_{\mathrm{R,in}}$
Total poloidal magnetic field |  |  |  |  |  | $\mathbf{B}_{\mathrm{P}}$ |  |  |  |  | $\mathbf{B}_{\mathrm{P,ex}}+\mathbf{B}_{\mathrm{P,in}}$
Total toroidal magnetic field |  |  |  |  |  | $\mathbf{B}_{\mathrm{T}}$ |  |  |  |  | $\mathbf{B}_{\mathrm{T,ex}}+\mathbf{B}_{\mathrm{T,in}}$
Radial field induced current |  |  |  |  |  | $\mathbf{J}_{\mathrm{PT}}$ |  |  |  |  | $\nabla\times\mathbf{B}_{\mathrm{R}}/\mu_{0}$
Poloidal field induced current |  |  |  |  |  | $\mathbf{J}_{\mathrm{TR}}$ |  |  |  |  | $\nabla\times\mathbf{B}_{\mathrm{P}}/\mu_{0}$
Toroidal field induced current |  |  |  |  |  | $\mathbf{J}_{\mathrm{RP}}$ |  |  |  |  | $\nabla\times\mathbf{B}_{\mathrm{T}}/\mu_{0}$
Hoop force density |  |  |  |  |  | $F_{\mathrm{H}}$ |  |  |  |  | $\mathbf{e}_{\mathrm{z}}\cdot(\mathbf{J}_{\mathrm{TR}}\times\mathbf{B}_{\mathrm{P,in}})$
Strapping force density |  |  |  |  |  | $F_{\mathrm{S}}$ |  |  |  |  | $\mathbf{e}_{\mathrm{z}}\cdot(\mathbf{J}_{\mathrm{TR}}\times\mathbf{B}_{\mathrm{P,ex}})$
Tension force density |  |  |  |  |  | $F_{\mathrm{T}}$ |  |  |  |  | $\mathbf{e}_{\mathrm{z}}\cdot(\mathbf{J}_{\mathrm{RP}}\times\mathbf{B}_{\mathrm{T}})$
Non-axisymmetry induced force 1 |  |  |  |  |  | $F_{\mathrm{N1}}$ |  |  |  |  | $\mathbf{e}_{\mathrm{z}}\cdot(\mathbf{J}_{\mathrm{PT}}\times\mathbf{B}_{\mathrm{P}})$
Non-axisymmetry induced force 2 |  |  |  |  |  | $F_{\mathrm{N2}}$ |  |  |  |  | $\mathbf{e}_{\mathrm{z}}\cdot(\mathbf{J}_{\mathrm{PT}}\times\mathbf{B}_{\mathrm{T}})$

  * a

The subscript “ex” refers to the field components of the potential field,
while the subscript “in” refers to the field components that are obtained
after subtracting the “ex” components from the total field.

  * b

Notations for the total poloidal magnetic field, $\mathbf{B}_{\mathrm{P}}$,
the total toroidal magnetic field, $\mathbf{B}_{\mathrm{T}}$, and their
decompositions are the same as that used by Myers et al. (2015).

Supplementary Movies

Supplementary Movie 1:  Temporal evolution of the magnetic field lines and the
electric current density. The movie displays a central part of the full
computation domain, namely, $X$ = [171′′, 382′′] and $Y$ = [$-$210′′, 105′′]
in the $X$–$Y$ plane. The vertical transparent slice displays the electric
current density. The cyan, yellow and orange lines have the same meaning as
that in Figure 2b, only that more field lines are drawn here. The olive lines
refer to field lines connecting the negative polarity N1 and the positive
polarity P2 at the first moment. The white isosurface represents the electric
current density larger than $32.9\%$ of the maximum value in the whole domain.
The background shows the distribution of the vertical magnetic field
component, $B_{z}$.
Supplementary Movie 2:  Temporal evolution of the SDO/AIA 304 Å images and the
signed QSLs with $|\mathrm{log}(Q)|>3$.
Supplementary Movie 3:  Temporal evolution of the QSLs (left) and the Lorentz
force, $L_{z}$ (right). The background at the bottom shows the distribution of
the vertical magnetic field component, $B_{z}$. The yellow and orange semi-
transparent isosurfaces represent the electric current density larger than
$21.9\%$ and $32.9\%$ of the maximum value in the whole domain, respectively.