Unveiling the Initiation Route of Coronal Mass Ejections through Their Slow Rise Phase

Understanding the early evolution of coronal mass ejections (CMEs), in particular their initiation, is the key to forecasting solar eruptions and induced disastrous space weather. Although many initiation mechanisms have been proposed, a full understanding of CME initiation, which is identified as a slow rise of CME progenitors in kinematics before impulsive acceleration, remains elusive. Here, with a state-of-the-art thermal magnetohydrodynamics simulation, we determine a complete CME initiation route in which multiple mainstream mechanisms occur in sequence yet are tightly coupled. The slow rise is first triggered and driven by the developing hyperbolic flux tube (HFT) reconnection. Subsequently, the slow rise continues as driven by the coupling of the HFT reconnection and the early development of torus instability. The end of the slow rise, i.e., the onset of the impulsive acceleration, is induced by the start of the fast magnetic reconnection coupled with the torus instability. These results unveil that CME initiation is a complicated process involving multiple physical mechanisms, thus being hardly resolved by a single initiation mechanism.


INTRODUCTION
The Sun frequently produces violent plasma eruptions such as coronal mass ejections (CMEs), which are released into the interplanetary space and induce geomagnetic storms (Gosling 1993) and damage aerospace equipment, satellite communications and power grids (Elovaara 2007) when colliding with the Earth's magnetosphere.The forecast of CMEs is thus extremely important for preventing space disasters and exploring habitable exoplanets (Khodachenko et al. 2007).
The evolution characteristics of CMEs and their pre-eruptive structures (e.g.sigmoids and filaments (Cheng et al. 2017); hereafter named as CME progenitors (Chen 2011)) in kinematics at various stages lay the foundation for CME predictions.For hours to days before the eruption, CME progenitors are stable and evolve quasi-statically, rising with a small velocity (< 1 km s −1 ) and tiny acceleration (Xing et al. 2018).In contrast, during the fast eruption, CMEs show an impulsive increase in velocity (up to hundreds to thousands of km s −1 ) and acceleration (up to hundreds of m s −2 ) in a short time Xing et al. (tens to hundreds of minutes) (Zhang & Dere 2006).These two stages of kinematics are named as the quasi-static phase and the impulsive acceleration phase, respectively (Zhang et al. 2001;Zhang & Dere 2006;Xing et al. 2018).The understanding of the kinematics during these two phases is relatively thorough, where the former is believed to be a response to the photospheric motion considering their similar characteristic speeds (Yang et al. 2004;Kaithakkal & Solanki 2019) and the latter is driven by the coupling of magnetohydrodynamics (MHD) instabilities (Kliem & Török 2006;Aulanier 2014;Green et al. 2018;Aulanier 2021) and fast magnetic reconnection (Lin & Forbes 2000;Zhang & Dere 2006;Jiang et al. 2021).
The initiation of CMEs refers to the transition between the quasi-static state of CME progenitors and the impulsively eruptive state of CMEs.In theoretical works, it usually corresponds to a turning point at the start of the impulsive acceleration phase, on the premise that the evolution of CME progenitors and CMEs is composed of the two phases mentioned above.Such a CME initiation is mainly explained by three mainstream mechanisms which are able to break the stable equilibria of CME progenitors (that is, to trigger the CME eruption) (Aulanier 2014(Aulanier , 2021)), i.e., (1) the torus instability (Kliem & Török 2006;Zuccarello et al. 2015) which refers to a pre-eruptive flux rope losing its equilibrium when the background constraint field decreases rapidly enough with increasing altitude, (2) the breakout reconnection (Antiochos et al. 1999;Masson et al. 2019) in which the external reconnection at a null point above the pre-eruptive structure removes the constraint of the overlying field, and (3) the tether-cutting reconnection (Moore et al. 2001;Jiang et al. 2021) in which the CME is caused by the fast reconnection with outflow jets.
However, some observational works demonstrate that CME progenitors actually always leave their quasi-static phase and experience a slow rise phase in a few to tens of minutes prior to the impulsive acceleration phase (Kahler et al. 1988;Zhang et al. 2001;Cheng et al. 2020;Prasad et al. 2023), specifically, manifesting as slowly rising with an increasing velocity (tens of km s −1 ) and a moderate acceleration (tens of m s −2 ) in this period (Cheng et al. 2020).Meanwhile, there is always a weak soft X-ray enhancement as the precursor of the main flare (Simnett & Harrison 1985;Zhang et al. 2001) and the CME progenitors often appear as high-temperature structures (known as hot channels) in extreme ultraviolet (EUV) bands (Cheng et al. 2012;Zhang et al. 2012).These observations indicate that the CME initiation, as the transition between the quasi-static state of CME progenitors and the impulsively eruptive state of CMEs, is actually a process before the impulsive acceleration phase, in which various precursors including the slow rise occur.
Different conjectures have been proposed in observations to understand the CME initiation as a process.On the one hand, it is argued that the eruption-trigger mechanism (e.g., torus instability, breakout reconnection, and tether-cutting reconnection) occurs at the start of the initiation process and drives the slow rise during the initiation (Zhang et al. 2001;Zhang & Dere 2006), considering the necessary early development of the eruption-trigger mechanism before inducing significant acceleration of CMEs (e.g., Schrijver et al. 2008).On the other hand, it is suggested that the CME initiation is related to the weak magnetic reconnection prior to the onset of eruption-trigger mechanisms, as the magnetic reconnection builds up the CME progenitor and reduces the overlying field restraining it (Sterling et al. 2007;Savcheva et al. 2012;Cheng et al. 2020; also see the review by Green et al. 2018).Specifically, Cheng et al. (2023) shown that the slow rise and the heating of CME progenitors are caused by a moderate precursor reconnection above the top of pre-flare loops.However, it is obvious that these two types of interpretations are disjointed, indicating that the initiation of a CME event is likely more complex than those described in above-mentioned conjectures and thus any single mechanism may only reveal a part of nature of the CME initiation.To determine a complete understanding of the CME initiation and to reveal how different physics are coupled with each other, we here perform a state-of-the-art three-dimensional (3D) observationallyinspired thermal-MHD simulation.The comprehensive analyses disclose a complete CME initiation that involves multiple physical mechanisms and even a coupling between two mechanisms.

Equations
In this work, we study the CME initiation by performing a 3D observationally-inspired thermal-MHD simulation labeled as "Simulation Driven-eruption" (abbreviated as "Simulation De") with the code MPI-AMRVAC (Xia et al. 2018), which solves the following equations in Cartesian coordinates: where ρ, p, T , e int , v, B, J represent the mass density, thermal pressure, temperature, internal energy, velocity, magnetic field and current density, respectively.g = −ge z is the gravity acceleration.µ is the dynamic viscosity coefficient and set to 10 −4 (in dimensionless unit) throughout the simulation, S is the strain tensor where , and I is the unit tensor.η is the resistivity coefficient whose value is set differently in different phases.κ || = 8 × 10 −7 T 5/2 erg cm −1 s −1 K −1 is the parallel conductivity coefficient, and b = B/B is the normalized magnetic field.We introduce a generalized Lagrange multiplier ψ to maintain the ∇ • B = 0 condition (i.e., generalized Lagrange multiplier (GLM) method; Dedner et al. 2002).The evolution of ψ follows Equation 6 and c h and c p are constants.As a thermal-MHD simulation, we consider compression heating, viscous dissipation, Ohmic dissipation, and thermal conduction in the energy equation.Especially, we solve the internal energy equation instead of total energy equation to avoid the heating from the numerical resistivity and numerical viscosity.

Numerical Methods
The equations are solved with the HLL scheme, the third-step Runge-Kutta time discretization method, and the fifth-order mp5-limited reconstruction (Suresh & Huynh 1997).We also use the magnetic field splitting method (Tanaka 1994;Xia et al. 2018) in this simulation: the magnetic field B is split into two parts, the invariable background field B 0 and the deviation B 1 .In this simulation,

Xing et al.
B 0 is set as the initial magnetic field, and thus the initial B 1 is 0. This method is conducive to eliminating numerical errors related to the divergence of magnetic field and achieving more precise results, especially in regions where the magnetic field evolves with minor change from the initial value, although it is relatively less effective in regions where the change of magnetic field is considerable (e.g., the centers of polarities during the converging phase).
The equations are solved in the dimensionless form, with the magnetic permeability equal to 1.The atmosphere is set to a fully ionized ideal gas with a hydrogen-helium abundance ratio of 10:1 (Xia et al. 2012).The dimensionless unit of the length, time, mass density, thermal pressure, temperature, velocity and magnetic field strength is 10 Mm, 67.89 s, 2.34 × 10 −15 g cm −3 , 0.51 erg cm −3 , 1.6 MK, 147.30 km s −1 and 2.53 G, respectively.
The physical domain of the simulation is in the range of −7 ≤ x ≤ 7, −7 ≤ y ≤ 7, and 0 ≤ z ≤ 14 (in dimensionless unit).We set a stretched grid (nx × ny × nz = 144 × 144 × 96), which is symmetric stretched in x and y directions with stretched ratio of 1.029, and unidirectional stretched in z direction with stretched ratio of 1.028.The spatial resolution in three directions is in the range of 0.0297 ≤ dx ≤ 0.2262, 0.0297 ≤ dy ≤ 0.2262, and 0.0298 ≤ dz ≤ 0.4103 (in dimensionless unit).In dimensional unit, the finest resolutions in three directions are about 300 km, comparable to those of Solar Dynamics Observatory/Atmospheric Imaging Assembly (AIA; Lemen et al. 2012).This spatial resolution is lower than those in other simulations of CME eruptions (Aulanier et al. 2010;Jiang et al. 2021), but it is a compromise for the sake of saving computation resources for solving the thermal-MHD equations.

Initial Conditions
The initial magnetic field (similar to that in OHM; Aulanier et al. 2010) is composed of two potential bipolar fields: where (c 1 = 60, x 1 = 0.9, The first (last) two monopoles are close to (far from) each other and locate shallow (deep), producing a strong (weak) potential field that decays rapidly (slowly) with altitude within the physical domain.The first potential field mimics the background field in an active region while the second one mimics the background field of the Sun on a global scale.The second potential field is set to ensure that the Alfvén speed at the upper part of the domain is large enough and thus to avoid shocks caused by the driving motion.The maximum dimensionless field strength in the plane z = 0 (which is named as box bottom surface in the following) is about 45.Such a strength is smaller than that in active regions in observations, while this is also a compromise for saving computation resources.The initial atmosphere is set to a hydrostatic corona with a uniform dimensionless temperature of 1.The box bottom surface corresponds to 3 Mm above the solar surface, and the initial dimensionless mass density at the box bottom surface is set to 1.The initial distribution of thermal pressure and mass density with altitude follows: where R sun = 69.55 is the dimensionless solar radius and g 0 = −0.126 is the dimensionless gravity acceleration at the solar surface.With such an initial atmosphere, the plasma β along the vertical slit through the origin is less than 1 in the region z < 3.7 (in dimensionless unit) at t = 0. While, as the magnetic field evolves, the plasma β further decreases: the plasma β along the vertical slit through the origin is less than 1 in the region z < 4.8 at t = 18, in the region z < 6.0 at t = 58, and in the region z < 7.5 at t = 68, which ensures that the modelled CME (progenitor) is always located in a low plasma β region which mimics the corona.The initial velocity is set to zero and the parameter ψ is also set to zero in the whole physical domain.

Driving Motions and Line-tied Conditions
In this work, we impose line-tied motions on the bottom boundary to drive the magnetic field and set special bottom boundary conditions to ensure the line-tied condition.Here, the line-tied condition means that the footpoint of the field line can only move horizontally according to the prescribed motion in the case of ideal MHD (Aulanier et al. 2005).
Inspired by the characteristics of the photospheric magnetic field evolution in observations (Yang et al. 2004;Green & Kliem 2009;Schrijver et al. 2011), two driving motions are considered in this simulation: a shearing motion that drives the initial potential field to a highly sheared state, and a converging motion that facilitates the flux cancellation near the polarity inversion line (PIL).The "Simulation De" is composed of two phases: the shearing phase (0 ≤ t ≤ 18) with only shearing motion applied and the converging phase (18 < t ≤ 68) with only converging motion applied.Such a setup differs from that in the previous work (Zuccarello et al. 2015), as here the converging motion is never switched off.We note that imposing only one type of motion at each phase makes the simulation easier to control and analyze: the shearing phase can be understood as a preparatory stage of the simulation to create a sheared magnetic field for the converging phase, during the latter of which we analyze the formation and eruption of the flux rope; while the converging phase is designed to fulfill the flux cancellation model.
The motions are imposed at the cell-center bottom surface which is a horizontal surface at the altitude of the cell center of the first layer of the physical domain (cell layers in physical domain indexed by k = 1, 2, 3, ..., 96 from the bottom to the top) at each time step to directly control the evolution of the magnetic field in the layer k = 1.For "Simulation De", we set: and the converging motion (v c x , v c y , v c z ): where and Ψ 1 = 5.5 for the shearing motion and Ψ 1 = 27.5 for the converging motion.The maximum speed of the driving motion v max 0 is set to 0.16 (in dimensionless unit) for both the shearing and the converging motions.The function Ψ 0 (t) is used to keep the maximum of the term V 2 x + V 2 y always v max 0 .As shown in Fig. 1, the shearing motion is along the tangent direction of the contour of B z (k = 1), while the converging motion is perpendicular to that tangent direction.Near the PIL, the shearing motion is antiparallel on two sides of the PIL, while the converging motion is towards the PIL.
It should be noted that the maximum speed of the driving motions is set to roughly 20 times of that in observations to speed up the simulation.While, it is still acceptable as it is less than the Alfvén speed and the sound speed in the core-area affected by the driving motion: at t = 0, in the domain that −4.20 ≤ x ≤ 4.20, −4.20 ≤ y ≤ 4.20, and 0 ≤ z ≤ 8.36 (in dimensionless unit), the ratio of v max 0 to the average Alfvén speed is 0.101 and the ratio of v max 0 to the maximum Alfvén speed is 0.004; the ratio of v max 0 to the sound speed is constant of 0.124 in this domain.Several settings are adopted to achieve the line-tied condition at the cell-center bottom surface.First, the GLM parameter ψ is fixed to zero in the layer k = 1, and we set a symmetric boundary condition centered at the cell-center bottom surface for ψ (see Section 2.5 for more details), so all components of the term ∇ • (ψI) is zero in the layer k = 1.Second, for "Simulation De", in the layer k = 1, the dissipation term in the Equation 3 is modified into: During the shearing phase of "Simulation De" (0 ≤ t ≤ 18), the dissipation term in this layer equals zero so that it has no effect on the magnetic field.The resistivity η is fixed to zero in the layer k = 1 and set to be uniform (η = 10 −4 ; in dimensionless unit) in the whole region (including both the physical domain and the ghost cell layers) except the layer k = 1 during the shearing phase.However, during the converging phase of "Simulation De" (18 < t ≤ 68), we adopt a two-dimensional (2D) dissipation term (Aulanier et al. 2010) in the layer k = 1 to allow the flux cancellation close to the PIL.A larger, uniform resistivity (η = 4 × 10 −4 ; in dimensionless unit), which facilitates the flux cancellation, is set in the whole region including the layer k = 1.
In addition, the shearing motion in Equation 9satisfies that the z-component of the term ∇ • (vB − Bv) is zero, so analytically B z in the cell-center bottom surface should remain invariable during the shearing phase.To achieve this feature more accurately, we enforce that B z remains invariable in the layer k = 1 during this phase.

Boundary Conditions
The mp5-limited reconstruction method uses three ghost cell layers at each boundary.The mesh in the ghost cell layer is also stretched in the same way as that in the physical domain.In all ghost cell layers, the background magnetic field B 0 remains invariable; the two components of the magnetic field B 1 which are parallel to the boundary surface are derived by one-sided third-order equal-gradient extrapolation, and the component of B 1 which is normal to the boundary surface is derived under the constraint of the ∇ • B = 0 condition.Such a setting will keep the divergence of the magnetic field zero in the layer k = 1, so it is reasonable to set ψ to zero in this layer.
We assume that at each time step, the atmosphere in both the bottom ghost cell layer and the layer k = 1 is hydrostatic.Thus, the thermal pressure in the bottom ghost cell layer can be derived with centered difference.We further assume that the dimensionless temperature in the ghost cell layer is constant of 1, same as the initial temperature, and then the mass density is derived by the ideal gas law.The thermal pressure and the mass density in the top ghost cell layers are derived with the same method.In the four side ghost cell layers, the thermal pressure and the mass density are fixed to their initial values.
The velocity in the bottom ghost cell layer k = n (layer indexed by k = −1, −2, −3 from the top to the bottom of the bottom ghost cells) is derived from the velocity in the layer k = 1 and that in the layer k = 1 − n by the linear extrapolation.This boundary condition is the same as the antisymmetric boundary condition for v z and the pseudo antisymmetric boundary condition for v x and v y at the bottom boundary in the previous work (Aulanier et al. 2005).The velocity in top ghost cell layers is derived by the one-sided second-order zero-gradient extrapolation, and the z-component is forced to be no less than zero to avoid possible downward flows.The velocity in the side ghost cell layers is set by the antisymmetric boundary condition centered on the side surface of the box.
To achieve the goal of removing the contribution of ∂ z ψ in the layer k = 1, a symmetric bottom boundary condition for ψ is set centered on the cell center of this layer.The aim is fulfilled more precisely by replicating the flux of ψ from the top surface of the layer k = 1 to its bottom surface (box bottom surface).In the top and the side ghost cell layers, ψ is fixed to zero.

Overview of Modelled CME Event
The evolution of the modelled CME progenitor and CME in "Simulation De" is shown in Fig. 2.During the shearing phase (0 ≤ t ≤ 18), the initial potential field is driven to a highly sheared state under the shearing motion (Fig. 2b).After that, during the converging phase (18 < t ≤ 68), a pre-eruptive flux rope composed of twisted field lines is first formed by magnetic reconnection under the converging motion, and then it erupts as a CME (Fig. 2b).The pre-eruptive flux rope is hot and visible in synthetic EUV images mimicking AIA 335 Å observations (Fig. 2a and Fig. 10).
The initial magnetic topology of the reconnection region is a bald patch (BP; Titov et al. 1993; during 20 ≤ t ≤ 45), characterized by (B • ∇)B z > 0 where B z = 0.As shown by the quasiseparatrix layers (QSLs; Priest & Démoulin 1995), such a BP topology later gradually bifurcates into a hyperbolic flux tube (HFT; Titov et al. 2002; since t = 46; see Fig. 5 and Fig. 8f).An obvious difference between the BP reconnection and the HFT reconnection is that the former forms flux rope field lines tangent to the bottom surface at their dips (Fig. 3a), while the latter forms flux rope field lines whose dipped parts are detached from the bottom surface and low-lying loops (Fig. 3b).

KINEMATICS OF CME PROGENITOR AND CME
The kinematics of CME progenitor and CME in "Simulation De" is estimated by measuring the height of the apex of an overlying field line right above the flux rope.Such an estimation method in simulations is consistent with the estimation method in observations where the kinematics of the hot channel is measured at its leading edge.The overlying field line is selected as one traced from a same fixed point at the center of the positive polarity and on the cell-center bottom surface.The velocity at this fixed point is extremely small, very close to or equal to zero, in all phases of "Simulation De", which ensures that the footpoint of this field line hardly moves and its height evolution can well reflect that of the CME (progenitor).
The kinematic of CME progenitor and CME, as shown in Fig. 4, is composed of three phases which are comparable to the three phases in observations, even though the simulation is merely designed for CME progenitor formation and CME eruption and certainly not to fit the kinematics a priori.The CME progenitor first rises quasi-statically with a tiny and fluctuating acceleration (during 35 ≤ t < 52; Fig. 4f), consistent with the behavior in the quasi-static phase in observations (Xing et al. 2018).As a result, the velocity-time curve is in a quasi-linear form in this phase (Fig. 4e).We note that the decaying fluctuation of the acceleration may reflect the process of the CME progenitor gradually approaching the quasi-equilibrium state under the driving motion.Later (during 52 ≤ t ≤ 65), the acceleration of the flux rope continuously increases to tens of m s −2 (Fig. 4f) and thus its velocity increases quickly (Fig. 4e).This phase corresponds to the slow rise phase in observations and thus marks the initiation process of the modelled CME, as this phase is just before the impulsive acceleration of the CME and the evolution of acceleration in it is significantly different from those in phases both before and after it.Even, the acceleration in this phase is comparable to that in the slow rise phase in observations (also tens of m s −2 ; Cheng et al. 2020).The last phase (65 < t ≤ 68) is considered to correspond to the impulsive acceleration phase in observations, as the modelled CME acceleration impulsively rises to hundreds of m s −2 in this period (Fig. 4c).
The eruption of CME is most likely triggered by the torus instability (TI), which occurs in the later stage of the slow rise phase (at a moment during 59 ≤ t < 63; see the analysis of torus instability in Appendix A).Other potential eruption-trigger mechanisms are ruled out: (1) The HFT reconnection can not trigger the eruption as the flux rope fails to erupt in the presence of the HFT in the control simulation "Simulation Ne" (see Appendix A). ( 2) The long and thin current sheet, which is usually necessary for the fast magnetic reconnection (Jiang et al. 2021), does not show up before t = 65 in this simulation.Therefore, the fast magnetic reconnection, even if it sets in this simulation, occurs later than the onset of the torus instability and does not play a key role in triggering the eruption.
Lastly, it should be noted that the behavior of the acceleration of the CME progenitor is less affected by the absolute speed of the driving motion during the quasi-static phase in this simulation.The reason is that the rising of the CME progenitor is a response to the driving motion in this phase, while the maximum speed of driving motions is a constant.Therefore, even if the absolute speed of the driving motion is large, the acceleration of the CME progenitor, as the first order derivative of the velocity with respect to time, is still almost close to 0. Therefore, the quasi-static phase in our simulation can be compared to that in observations.

Triggering Mechanism of Slow Rise Phase
The transition from the quasi-static phase to the slow rise phase occurs (at t = 52) just after the first appearance of HFT (at t = 46).This indicates that the HFT reconnection is likely the mechanism triggering the slow rise phase.
This speculation is supported by detailed analyses on the relation between the reconnection and the kinematics.Fig. 5a,f reveal that the upward outflow of the reconnection is very slow at the stage of BP (e.g., t = 45) and the early stage of HFT (46 ≤ t < 52), but it is markedly accelerated as the HFT is rapidly built up (e.g., from t = 52 to t = 58).Quantitatively, the growth rate of the reconnection upward outflow velocity is temporally consistent with the acceleration of the CME progenitor, and both of them start to increase shortly after the first appearance of HFT (Fig. 5g).We consider that the velocity of the reconnection upward outflow reflects the ability of reconnection to contribute to the flux rope.Therefore, the above results indicate that the slow rise is triggered when the contribution from reconnection starts to quickly increase after the HFT forms.In other words, the slow rise is triggered by the HFT reconnection.We eliminate the possibility that the torus instability triggers the slow rise, since the failed eruption in the "Simulation Ne" implies that the torus instability, as the triggering mechanism of the eruption, does not occur before t = 58 in the "Simulation De" (see Appendix A).

Driving Mechanisms of Slow Rise Phase
We further investigate the mechanisms driving the slow rise of the flux rope.Considering that the torus instability occurs during the slow rise phase, we divide the slow rise phase into two stages, i.e., the earlier stage (52 ≤ t ≤ 58) and the later stage (59 ≤ t ≤ 65), and analyze the driving mechanisms in these two stages separately.
For the earlier stage, the motion of the flux rope is controlled by a net upward force on its lower part and a net downward force on its upper part (see Fig. 6b, for an example at t = 58).The Lorentz force and the thermal pressure gradient force are the major contributors to the net force, while the gravity and the viscous force are of little importance (Fig. 6).Especially, the maxima of the Lorentz force, the net force, and the velocity in the flux rope region are all concentrated near the HFT (Fig. 6a,b) during this stage.This implies that it is the HFT reconnection that powers the flux rope rise in the earlier stage of the slow rise phase.In addition, the synchronized evolutions of the flux rope acceleration and the growth rate of the reconnection upward outflow velocity in Fig. 5g also indicate that the slow rise is driven by the HFT reconnection during the earlier stage.
The HFT reconnection drives the slow rise in the earlier stage, on the one hand, by adding twisted concave-upward field lines to the lower part of the flux rope (Fig. S3).These field lines provide magnetic tension to the flux rope (Fig. 6a; known as the slingshot effect), enhance the force imbalance in the flux rope, and finally accelerate the flux rope.On the other hand, the reconnection upward outflow in the HFT is convected upwardly into the flux rope (Fig. 5a,e), and thus it can directly promote the flux rope rise.Specifically, the plasma in the reconnection upward outflow is accelerated by a net upward force which is dominated by the Lorentz force and the thermal pressure gradient force (Fig. 5).The upward Lorentz force consists of magnetic tension force mainly at the HFT top and magnetic pressure gradient force mainly at the HFT center (Fig. 7), the former of which is provided by the newly formed concave-upward field lines and the latter of which is associated to the accumulation of the sheared/guide component of the magnetic field.The upward thermal pressure gradient force is due to an increase in the gas pressure at the flanks of the flux rope bottom, which is a result of the density pileup pinched by the reconnection inflow (Fig. 7).
Furthermore, since the HFT appears, the two footprints of the HFT continuously separate from each other in the direction perpendicular to the PIL, and meanwhile they extend in the direction parallel to the PIL to both sides until they are completely separated (Fig. 8a-c; shown by the QSL footprint, which manifests as the flare ribbon in observations during the eruption (Zhao et al. 2016), but probably invisible before the eruption due to the weak reconnection rate).The reconnection electric field within the HFT, which represents the reconnection rate (Forbes & Priest 1984), also shows a gradually increasing trend since t = 46 and until t = 59 (Fig. 8h; also see Appendix C for more details).These results strongly suggest that the HFT and its associated reconnection are gradually built up during the earlier stage of the slow rise phase.Both the slingshot effect and the outflow effect of the HFT reconnection thus become stronger, which leads to the continuous increase of acceleration of the CME progenitor in this stage.
For the later stage of the slow rise phase, we argue that the HFT reconnection still has a considerable contribution to the rising motion, since the HFT reconnection becomes even more stronger in this stage, shown by the quickly increasing reconnection electric field (see that during 59 ≤ t ≤ 65 in Fig. 8h).Moreover, the torus instability also promotes the slow rise in the later stage since it starts.On the one hand, the hoop force could directly accelerate the flux rope (Kliem & Török 2006); on the other hand, the torus instability could indirectly contribute to the slow rise by enhancing the HFT reconnection, evidenced by that the HFT reconnection electric field immediately shows a quick increase once the torus instability sets in (Fig. 8h).This enhancement may be due to that the current structure in the HFT could become longer as the flux rope rises higher and faster with the contribution from the torus instability, which well reflects the close coupling of torus instability and HFT reconnection in the later stage.

Triggering and Driving Mechanisms of Impulsive Acceleration Phase
A fast magnetic reconnection likely occurs during the impulsive acceleration phase, as a thin and long current sheet, which is usually the necessary condition for the fast magnetic reconnection, is formed at the start of the impulsive acceleration phase (t ∼ 66) and soon built up since then (Fig. 9).The other evidence for the fast magnetic reconnection is the impulsively increasing reconnection electric field in the reconnection region during the impulsive acceleration phase (Fig. 8g).Therefore, we consider that the impulsive acceleration phase (65 < t ≤ 68) is likely triggered as the magnetic reconnection coupled to the torus instability transforms from the relatively weaker HFT reconnection to the relatively stronger fast magnetic reconnection, and that the phase is continuously driven by the close coupling of the fast magnetic reconnection and the torus instability which induces an impulsive acceleration.

HEATING OF CME PROGENITOR
Benefitting from the thermal-MHD simulation, we also study the thermal properties of the CME progenitor.The modelled CME progenitor is significantly heated and visible in EUV synthetic images at the AIA 335 Å (Fig. 2a and Fig. 10; see Appendix D for more details).It is channel-like in the face-on view, elliptical in the edge-on view, and sigmoidal in the top view.It is very similar to preeruptive hot channels discovered in observations (Zhang et al. 2012;Patsourakos et al. 2013;Cheng et al. 2014), except that our modelled one is somewhat less hot than those in observations due to the weak magnetic field adopted in our simulation.The temperature of the modelled CME progenitor is expected to increase with a stronger magnetic field, as the former is roughly proportion to the square of the later.Regardless of the specific temperature, the modelled bright structure basically corresponds to the pre-eruptive flux rope (Fig. 10), which supports the view that the hot channel represents the hot flux rope in observations (Zhang et al. 2012).
In addition, taking the snapshot at t = 58 as an example, an analysis of the heating sources reveals that the CME progenitor is mainly heated by the Ohmic dissipation (Fig. 11).Especially, one can find that at the bottom and side boundaries of the flux rope, the Ohmic dissipation is relatively stronger and its profile is highly consistent with that of the HFT and the QSLs (Fig. 11).This result indicates that the magnetic reconnection in the HFT and the QSLs surrounding the flux rope effectively heats the local plasma while it forms the pre-eruptive flux rope, which confirms the conjecture in the previous 3D CME model (Janvier et al. 2014).In addition, the bulk heating within the flux rope results from the setup of a uniform resistivity in the simulation domain.

SUMMARY: A COMPLETE ROUTE OF CME INITIATION
Based on our simulation, we present a complete route of the CME initiation as summarized in Fig. 12.In short, during the CME initiation, the slow rise motion is first triggered and driven by the HFT reconnection and later driven by an early coupling of the HFT reconnection and the torus instability.
In detail, the kinematic evolution around the CME initiation is described as the following.With the magnetic topology of the reconnection region changing from the BP to the HFT, the HFT reconnection starts to promote the flux rope rising with its slingshot effect and reconnection upward outflow.As the HFT reconnection is developing, the effect of the HFT reconnection to the flux rope rise gradually increases and eventually starts the slow rise phase/CME initiation.Afterwards, at the earlier stage of the slow rise phase/CME initiation, the pre-eruptive flux rope rises slowly with an increasing acceleration owing to the progressively enhanced HFT reconnection.With the HFTdominated slow rise going on, the flux rope later reaches a critical height where the torus instability sets in.The onset of the instability marks the onset of the CME eruption in the physical sense, i.e., the flux rope being out of equilibrium.Then, at the later stage of the slow rise phase/CME initiation, the torus instability induces a stronger HFT reconnection, and they together promote the slow rise of the flux rope.Finally, the impulsive acceleration of CME starts (the CME initiation ends) when the magnetic reconnection coupled to the torus instability transforms from the HFT reconnection to the fast magnetic reconnection, which is later than the eruption onset in the physical sense.

DISCUSSION
In previous studies, the slow rise motion during the initiation process was explained by a single mechanism such as the eruption-trigger mechanism (Zhang et al. 2001;Zhang & Dere 2006) or the weak magnetic reconnection before the onset of eruption-trigger mechanisms (Sterling et al. 2007;Savcheva et al. 2012;Cheng et al. 2020Cheng et al. , 2023)).In contrast, our work presented here reveals that the CME initiation is actually a multiple-physics coupled-process and its complete route is hardly to be explained by the previous single mechanism.Firstly, the HFT reconnection triggers the CME initiation and drives the CME progenitor to slowly rise in quasi-equilibrium.Secondly, after the torus instability triggers the eruption, the less energetic destabilization and the HFT reconnection jointly contribute to the slow rise/CME initiation for a while until the coupling of the destabilization and the fast magnetic reconnection begins and induces the impulsive acceleration of the CME.
For the HFT reconnection, our simulation especially provides insight into the question how it drives the flux rope during the initiation process.We find that both the Lorentz force and the thermal pressure gradient force play important roles in the slow rise phase, at least in its earlier stage.Quantitatively, for the flux rope field lines at the upper part of the HFT, the ratio of the maximum of the thermal pressure gradient force to that of the Lorentz force is roughly 40% at t = 52 and 12% at t = 58 (Fig. 5c,d; the plasma β is roughly 0.01 at regions of the above-mentioned maximum of the force).This result for the first time confirms the previous conjecture on the key role of Lorentz force in causing the rise of pre-eruptive structures in the "tether-cutting reconnection" model (Moore et al. 2001).However, the contribution of thermal pressure gradient force to the slow rise is also not negligible, while this result could be further influenced by the plasma β in the HFT reconnection region.This underscores the indispensability of a thermal-MHD simulation, beyond the zero-β MHD simulations (Aulanier et al. 2010;Zuccarello et al. 2015), in the study of the CME initiation.
Another important finding in our work is that the impulsive acceleration of CMEs does not immediately starts at the onset of the torus instability but commences at a later moment when the coupling of the torus instability and the fast magnetic reconnection starts.Such a time delay may explain why the critical decay indices estimated in observations (Zou et al. 2019;Cheng et al. 2020Cheng et al. , 2023) ) are sometimes larger than those derived theoretically (Kliem & Török 2006;Démoulin & Aulanier 2010), as the observationally determined eruption onset time is probably later than the real onset time of the torus instability in these works.It should be mentioned that a larger magnetic field strength in observations compared to that in this simulation may lead to a larger kinematic acceleration, and thus the above-mentioned time delay may become shorter in observations.In addition, such a time delay indicates that the torus instability can only induce a relatively small acceleration at its beginning, compared to that in the impulsive acceleration phase.Lastly, we emphasize that the finding on the time delay benefits from the setup of our simulation, i.e., continuously imposing driving motion after the eruption is triggered.Such a setup is distinguished from that in previous works (Aulanier et al. 2010;Zuccarello et al. 2015), in latter of which the driving motion is switched off to relax the system during the eruption.Our simulation avoids the extra relaxation effect in previous works, which does not occur in observations but could also contribute to the time delay in simulations and pose a significant difficulty for analysis (as discussed in Appendix A.2).
Our simulation sheds light on the heating process of CME progenitors.However, it should be noted that the thermodynamics in our simulation is a simplification of that of the Sun for the sake of clear and concise results.First, the density and temperature profile of the atmosphere was simplified in the simulation by omitting the lower atmosphere (from the photosphere to the transition region) and setting a uniform coronal temperature.Second, some processes in flares, e.g., the radiative cooling and the chromospheric evaporation, were not taken into account in the simulation.Even so, we emphasize that these simplifications do not critically affect our result, i.e., the magnetic reconnection takes place in the current sheet surrounding the flux rope and builds up and heats the CME progenitor.
In summary, we successfully reproduce the complete early kinematics of a CME event with a 3D thermal-MHD simulation that is more physically realistic in this work.For the first time, we explicitly disentangle the complex CME initiation process by systematically investigating how various mechanisms play roles at different stages of the kinematics of a CME event.
Xing et al.      (progenitor).The red, green, purple, blue curves show the four accelerations contributed by the torus instability, the HFT reconnection, the fast magnetic reconnection, and the driving motion, respectively.The onset times of the HFT, the slow rise, the torus instability, and the impulsive acceleration (abbreviated as acc.) of the CME are displayed by four dashed lines and labeled as t HFT-formation , t slow-rise-beginning , t onset-TI , and t impulsive-acc.-beginning, respectively.In our simulation, these four onset times correspond to t = 46, t = 52, a moment during the period 59 ≤ t < 63, and t = 65, respectively.The light grey and dark grey regions show the earlier stage and the later stage of the slow rise phase/CME initiation process, respectively.The red curve is dashed during a part of the slow rise phase and the impulsive acceleration phase, indicating that the dominant mechanism for the acceleration is unknown whether it is torus instability or magnetic reconnection in these periods.

Xing et al.
index for the torus instability onset should be in a range quite similar to 1.5227 ≤ n ≤ 1.6239.In the following, we take the range 1.5227 ≤ n ≤ 1.6239 for the critical decay index in "Simulation De".In Fig. S1d, it is clear that the height of the flux rope axis apex in "Simulation De" is already higher than the height of the upper limit of the critical decay index at t = 63, indicating that the torus instability could set "Simulation De" and that its onset time should be before t = 63.The earliest onset time of the torus instability in "Simulation De" is still t = 59, as the apex of the flux rope axis reaches the height of the lower limit of critical decay index at this moment.Therefore, we summarize that the torus instability sets in "Simulation De" during 59 ≤ t < 63.
We note that it is not an accident that the time range of the torus instability onset in "Simulation De" is earlier than that in "Simulation Ue".In the period 60 < t ≤ 64, the driving motion becomes stronger and stronger from "Simulation Ue" to "Simulation He" to "Simulation De", with the maximum speed of the driving motion varying from 0 to 0.08 to 0.16 (in dimensionless unit).In Fig. S1e,f, we show that, on the one hand, for these three simulations, the flux rope in the simulation with a stronger driving motion could rise to a higher altitude at the same moment (taking examples at t = 63 (see panel e) and at t = 64 (see panel f)).On the other hand, in comparison among the three simulations, the one with a stronger driving motion has a higher height corresponding to the same critical decay index at the same moment (see panels e and f).However, it is clear that such a difference in the height of the flux rope axis is much larger than the difference in the height of a certain decay index, when comparing each two of these three simulations.This suggests that in comparison among different simulations, the stronger driving motion imposed, the earlier the flux rope axis reaches the height of the critical decay index and thus the earlier the torus instability starts, if regardless of the change in the critical decay index itself.

B. FLUX ROPE BOUNDARY
The flux rope boundary and the reconnection region are identified with the squashing degree Q which measures the mapping of the field lines.The squashing degree Q is derived by (Titov et al. 2002): with the code FastQSL (Zhang et al. 2022), where (x, y) and (X, Y ) are coordinates of two footpoints of a field line.
The QSLs refer to the locations where logQ is much larger than 2 (Titov et al. 2002), and here we define the QSLs as the region where logQ ≥ 3.In the following we take the flux rope at t = 58 in "Simulation De" as an example to explain how we determine the flux rope boundary.In the plane y = 0, the bottom and side boundaries of the flux rope can be exactly determined by the QSLs, while the top boundary is blurred as the QSLs are not closed there (Fig. S3).We determine the top boundary qualitatively at the top edge of the strong downward Lorentz force region, i.e., the apex of the flux rope is roughly at z = 3.3 (in dimensionless unit; see Fig. 6a).This determination is reasonable for the following reasons.First, as shown in Fig. S3, the weakly twisted field lines traced from the downward Lorentz force region (yellow field lines) are mostly anchored in the regions partially surrounded by the QSL footprints.They compose a flux rope configuration together with the highly twisted field lines traced from the upward Lorentz force region (green field lines).Second, the tops of the high mass density region (represented by the gravity in Fig. 6c) and of the high current density region (represented by the Ohmic heating rate in Fig. 11) are also around z = 3.3 in the plane y = 0, also indicating that the flux rope apex is around z = 3.3.

C. ESTIMATING THE RECONNECTION ELECTRIC FIELD
We estimate the reconnection electric field with a method commonly used in observations (Qiu et al. 2002).In brief, the reconnection electric field (E) is derived by E = v Q B z , where v Q is the separation speed of the QSL footprint in the direction perpendicular to the PIL and B z is the z-component of the magnetic field at the QSL footprint.This estimation method is based on a 2D flare model but is considered still applicable to a 3D situation (Forbes & Lin 2000).We refer readers to the previous paper (Qiu et al. 2002) for more details of this method.
In Fig. 8a-e, we show the separation motion of the QSL footprint on the cell-center bottom surface.To estimate the separation motion speed, we set a slit (represented by the orange dashed lines in Fig. 8a-c) on this surface, starting from the origin and along the direction perpendicular to the PIL.The time-slice plot of the normalized logQ (Fig. 8f), where the logQ on the slit is normalized by its maximum at each moment, exhibits how the QSL footprint separates from the PIL.The center of the QSL footprint, where the normalized logQ reaches its local maximum, is marked by the red symbol at each moment.The phenomenon that the QSL footprint stays on the PIL before and at t = 45 confirms that the magnetic topology of the reconnection region is a BP in this period.The separation of QSL footprint starts at t = 46, marking the first appearance of HFT at this moment.The separation speed (v Q ) is derived to be the derivative of the distance from the origin to the QSL footprint with time.Finally, the electric field of the HFT reconnection is derived by multiplying the separation speed with the local B z .The electric field during 46 ≤ t ≤ 51 (see Fig. 8g,h) is represented by that derived from the orange slit, considering the reconnection still occurs in BP rather than HFT in many other places in this period.The electric field during 52 ≤ t ≤ 68 (also see Fig. 8g,h) is represented by the average of the electric fields derived from the orange slit and three other slits (represented by the pink dashed lines in Fig. 8a-c; all of them crossing through HFT footprints rather than BP footprints since t = 52).The aim of the multiple measurements and averaging is to better show the evolution of the electric field in the entire HFT.

D. SYNTHETIC EUV IMAGES OF THE CME PROGENITOR
We synthesize the EUV images at the AIA 335 Å as observed from three side views through integrating the emissivity along each direction under the optically thin emission assumption.The passband is selected at 335 Å to best demonstrate the structure of the hot CME progenitor.The emissivity is derived from the temperature and the number density of electron (which are dimensionalized with the dimensionless units; the number density is derived under the assumption of a fully ionized ideal gas with a hydrogen-helium abundance ratio of 10:1) with the AIA response function (performed by SSW package code "aia get response.pro"),and the integration is performed in a box where −3.2 ≤ x ≤ 3.2, −4.9 ≤ y ≤ 4.9, and 0.015 ≤ z ≤ 3.86 (in dimensionless unit).To better show the flux rope structure, we show the percentage difference intensity rather than the original synthetic intensity, which is derived by: where D i and I i are the percentage difference intensity and the original synthetic intensity at t = t i , respectively, and I 0 is the original synthetic intensity at t = 0. We note that the percentage difference Xing et al.
intensity in the flux rope, although sometimes negative, is still larger than that of its surrounding, clearly showing the presence of the hot CME progenitor.

Figure 1 .Figure 2 .Figure 3 .Figure 4 .Figure 5 .
Figure 1.Distribution of B z and imposed horizontal driving motions on the cell-center bottom surface.The black (white) contours in the positive (negative) polarity are the contours of B z , and the purple ones show the PILs.The red (green) arrows in the positive (negative) polarity in the left panel show the shearing flow at t = 15, and those in the right panel show the converging flow at t = 21.The arrow orientation denotes the direction of the motion, and its length represents the magnitude of the speed (v h ) .The scale of v h is shown at the bottom right corner of each panel.All parameters are in dimensionless units. 1

FLorentzFigure 6 .Figure 7 .
Figure6.Component of forces and velocity in the z-direction in the plane y = 0 at t = 58.They are in order as follows: Lorentz force (F Lorentz ), magnetic tension force (F mag−tens ), magnetic pressure gradient force (F mag−pres ), thermal pressure gradient force (F ther−pres ), net force (F net ), velocity (v z ), gravity (F grav ), viscous force (F visc ), and force related to the convection (F convect ).F net is the sum of F Lorentz , F ther−pres , F grav , and F visc .Here the F mag−tens (F mag−pres ) is derived by first decomposing B • ∇B (−∇(B 2 /2)) to the normal direction of magnetic field and then further decomposing to the z-direction, given that the tangential component does not contribute to the acceleration.All of forces and velocity are in dimensionless units.The green contours represent the profiles of QSLs in the half panel of x ≤ 0. See Appendix B for more details for locating the flux rope boundary.

Figure 8 .Figure 10 .
Figure 8. Evolutions of HFT footprint and HFT reconnection electric field.(a)-(c) Squashing degree logQ on the cell-center bottom surface at t = 52, t = 55, and t = 58.The field of view is marked by the red dashed box in panel d.(d)-(e) Squashing degree logQ on the cell-center bottom surface at t = 58 and t = 63.The two hook-shape QSL footprints around (X, Y ) = (1, 2) and (X, Y ) = (−1, −2) mark the boundaries of two footpoints of the flux rope, respectively.At the ends of straight parts of QSL footprints (around (X, Y ) = (−2, 2) and (X, Y ) = (2, −2) in panel e), the curved QSL footprints demonstrate the rapid change in the connectivity of the field lines anchored there as a result of the driving motion.(f) Time-slice plot of the normalized logQ.The slit is denoted by the orange dashed line in panels a-c.At each moment, logQ is normalized by the maximum logQ along the slit, and the red symbol marks the location of the local maximum of normalized logQ which corresponds to the HFT footprint.(g) Reconnection electric field E during 46 ≤ t ≤ 68.The reconnection electric field during 46 ≤ t ≤ 51 is derived from the orange slit in panels a-c; the reconnection electric field during 52 ≤ t ≤ 68 is averaged from the electric fields derived from four slits (represented by the orange and the pink dashed lines in panels a-c).The vertical dashed lines, light/dark grey regions, and yellow oblique-line-region have the same meanings as those in Fig. 4. (h) Same as panel g but for the reconnection electric field E during 46 ≤ t ≤ 65.

Figure 11 .Figure 12 .
Figure11.Heating rates on the plane y = 0 at t = 58.They are in order as follows: rates of Ohmic heating (Q Ohmic ), viscous heating (Q visc ), compression heating (Q compress ), and net heating (Q net ) which is the sum of the first three.All heating rates are in dimensionless unit.The field of view is the same as that of the synthetic image along the y-direction in Fig.2a.The green contours represent those of QSLs at the half panel of x ≤ 0.

Figure S2 .
Figure S2.Top views of magnetic field lines in Fig. 2b.The left panel is at t = 57 and the right panel is at t = 64.The vertical direction is the y-direction.