Data-constrained Magnetohydrodynamic Simulation of an Intermediate Solar Filament Eruption

Solar eruptive activities could occur in weak magnetic field environments and over large spatial scales, which are especially relevant to eruptions involving intermediate or quiescent solar filaments. To handle the large scales, we implement and apply a flux rope embedding method using regularized Biot–Savart laws in the spherical coordinate system. Combined with a potential field source surface model and a magneto-frictional method, a nonlinear force-free field comprising a flux rope embedded in a potential field is constructed. Using the combined nonlinear force-free field as the initial condition, we then perform a zero-β data-constrained magnetohydrodynamic (MHD) simulation for an M8.7 flare at 03:38 UT on 2012 January 23. The MHD model reproduces the eruption process, flare ribbon evolution (represented by the quasi-separatrix layer evolution), and kinematics of the flux rope. This approach could potentially model global-scale eruptions from weak field regions.


INTRODUCTION
Solar eruptive activities in the solar atmosphere, such as coronal mass ejections (CMEs), flares and filament/prominence eruptions, are the sources of disastrous space weather.Observing and modeling these eruptions are an effective means to understand, predict, and reduce their adverse impacts.
Through long-lasting studies in the past, a type of magnetic field configuration, either a flux rope or sheared arcade, can become unstable and erupt, driving various aforementioned solar activities (Chen 2011;Green et al. 2018;Toriumi & Wang 2019).A flux rope is a bundle of magnetic field lines winding around an axis.If the twist is low, say, less than 1 turn, or the toroidal flux is much larger than the poloidal flux, the bundle of field lines becomes a sheared arcade (Patsourakos et al. 2020).
Local coronal plasma may cool down and condense on a flux rope or sheared arcade, which appears as a filament or prominence (Chen et al. 2020).A successful filament/prominence eruption is usually accompanied by a CME and flare with the erupting magnetic field and plasma, which involves an induced shock, more complex electric currents, and accelerated particles.
Many studies revealed that magnetohydrodynamic (MHD) instabilities and magnetic reconnection play important roles in triggering and driving flux rope eruptions (Antiochos et al. 1999;Chen & Shibata 2000;Moore et al. 2001;Török & Kliem 2005;Kliem & Török 2006;Keppens et al. 2014;Ishiguro & Kusano 2017).MHD instabilities and reconnection are regarded as ideal and resistive processes, respectively.The origin of the MHD instabilities is frequently the Lorentz force, and that of the reconnection is the electron and ion dynamics and electric current dissipation mechanism.A powerful tool to study the realistic triggering and driving processes is data-driven or data-constrained MHD models (Kliem et al. 2013;Jiang et al. 2016;Amari et al. 2018;Inoue et al. 2018;Guo et al. 2019a;Zhong et al. 2021).These models incorporate observations to constrain the initial and boundary conditions, and can reproduce many key observational features, such as the morphology of flare ribbons, kinematics of erupting filaments, and so on (Jiang et al. 2022).
Depending on the locations of their source regions, solar filaments are divided into three categories, namely, active region, intermediate, and quiescent ones (Engvold 1998).An intermediate filament has one footpoint rooting in an active region and the other footpoint in a quiet region.Traditional nonlinear force-free field (NLFFF) extrapolation techniques need extra flux rope embedding strategies, as they will not easily produce a proper flux rope model for an intermediate or quiescent filament as that for an active region filament.Since the magnetic field under an intermediate or quiescent filament is usually not very strong and they mostly detach from the photosphere, the transverse components of the vector magnetic field are too weak to constrain highly twisted features in the corona.
For example, statistics conducted by Ouyang et al. (2017) found that 96% of quiescent filaments are supported by twsited flux ropes, while only 4% of them are supported by sheared arcades.However, as for active-region filaments, 40% of them are supported by sheared arcades.There are only a few models that have reconstructed flux ropes in weak field regions (Su & van Ballegooijen 2012;Jiang et al. 2014;Mackay et al. 2020).The regularized Biot-Savart laws (RBSL) proposed by Titov et al. (2018) provide a recipe to address the issue.The key idea is to embed a magnetic flux rope (B FR ) into a potential field (B P ).In this case, the combined magnetic field is not force-free, since the Lorentz force component exerted by B P on the flux rope generally cannot balance the hoop force of the flux rope due to its curvature.It should be further relaxed to a force-free state by, e.g., a magneto-frictional method (Guo et al. 2016a,b) or a zero-β MHD relaxation method (Titov et al. 2021;Kucera et al. 2022).Guo et al. (2019b) implemented the RBSL method in the Message Passing Interface Adaptive Mesh Refinement Versatile Advection Code (MPI-AMRVAC; Keppens et al. 2003Keppens et al. , 2012Keppens et al. , 2021Keppens et al. , 2023;;Porth et al. 2014;Xia et al. 2018), and thus far used it in a Cartesian coordinate system.For example, Guo et al. (2021b) simulated a flux rope eruption using the RBSL model embedded in a potential field as the initial condition.Here, we extend the RBSL model to the spherical coordinate system and present a data-constrained zero-β MHD model to study the eruption of an intermediate filament.The motivating observation is described in Section 2. Our MHD simulation setup is described in Section 3.
The results are presented in Section 4. A summary and discussion are provided in Section 5.

OBSERVATION
In this study, we focus on an M8.7 flare starting at 03:38 UT on 2012 January 23 in active region 11402, which is located in an active region cluster comprising three other regions (active regions 11401, 11405 and 11407).The coronal magnetic field of this active region cluster was extrapolated before by Guo et al. (2012) with the optimization method (Wheatland et al. 2000;Wiegelmann 2004).Guo et al. (2016a) developed another NLFFF model for this active region cluster using the magnetofrictional method (Yang et al. 1986) developed in MPI-AMRVAC.A flux rope was found in active region 11402.However, if we compare carefully the constructed flux rope with observations by the Atmospheric Imaging Assembly (AIA; Lemen et al. 2012) aboard the Solar Dynamics Observatory (SDO; Pesnell et al. 2012), as shown in Figure 9 of Guo et al. (2016a), the length of the rope is much shorter than the observed filament.To improve on this aspect, we try to make a new model in this study.Moreover, we want to simulate the dynamic process in active region 11402 that causes the M8.7 flare.
Figure 1 shows the evolution of the M8.7 flare on 2012 January 23.A prominent filament is suspended in the corona at 03:00 UT (Figures 1a, 1b, and 1c), well before the flare starts at 03:38 UT.There are some coronal loops overlying the filament as shown in the SDO/AIA 171 Å image (Figures 1a).These loops lean towards the east compared to the radial direction in the active region, which is towards the west as indicated by the red and green arrows in Figure 1a.The filament starts to lift up at 03:35 UT, and brightened flare ribbons can be seen in different EUV wavebands (Figures 1d,1e and 1f).At 03:42 UT, the erupting filament has reached and stretched the visible overlying coronal loops (Figures 1g, 1h and 1i).The eruption process has been analyzed by previous observational studies (Cheng et al. 2013;Joshi et al. 2013;Sterling et al. 2014).

MHD SIMULATION
Observations show that the pre-eruptive phenomenon of the source region is dominated by an intermediate filament.We first construct a flux rope model for the filament using the RBSL method.
Then, we embed the RBSL flux rope into a potential field derived by the PFSS model.The combined RBSL and PFSS magnetic field is relaxed by the magneto-frictional method to derive an NLFFF model, which further serves as the initial condition for the following zero-β MHD simulation.
3.1.Regularized Biot-Savart Laws in Spherical Coordinates Titov et al. (2018) proposed the RBSL method for constructing a flux rope B FR with an axis path C, minor radius a, magnetic flux F , and electric current I: where x is the position vector of the magnetic field, l is the arc length of the axis path, R(l) is the position vector of the axis path, r ≡ r(l) = (x − R(l))/a(l) is the separating vector from the source point R(l) to the field point x and normalized by a(l), R ′ = dR/dl is the tangential unit vector, and C * is the mirror image of C referred to the photosphere.For an axial current density distributed parabolically along the flux rope minor radius, the integration kernels are obtained by Titov et al. (2018): (5) Once the free parameters including the axis path C, minor radius a, magnetic flux F , and electric current I are determined, we could derive the flux rope magnetic field B FR using Equations ( 1)-(3).
The RBSL method was proposed in the Cartesian coordinate system (Titov et al. 2018).To implement it in a spherical coordinate system, we require that the distance between the two footpoints of the modeled flux rope is as small as possible, though a distance less than or close to the solar  and F2) of the filament on the photosphere, the X-axis is perpendicular to both OO ′ and the line joining the two footpoints, the Y -axis is parallel to the line joining the two footpoints, and the Z-axis constitutes a right-hand system together with the X-and Y -axes and pointing radially outward.
radius is tolerable.Therefore, the two footpoints are close enough to keep the radial magnetic field not changing too much on the photosphere when we make the coordinate transformation.We define two coordinate systems, one is the heliocentric O-xyz and the other is the local O ′ -XY Z as shown in Figure 2. Since the two coordinate systems are determined, we can convert any useful values between them, such as the coordinates themselves and the magnetic field vecctor components.The idea is to provide the intended spherical coordinates (r, θ, ϕ), which are converted to the coordinates in O-xyz and further to O ′ -XY Z, then to compute the RBSL magnetic field at these locations in O ′ -XY Z, and finally to convert the magnetic field back to O-xyz and the spherical coordinates.
Computing the RBSL magnetic field in O ′ -XY Z is exactly the same as done in a Cartesian coordinate system (Titov et al. 2018), which has been implemented in MPI-AMRVAC (Guo et al. 2019b).We have to determine the three dimensional (3D) path of the filament C, its minor radius a, magnetic flux F , and electric current I. Figure 3 shows how the path is determined.Two footpoints are pinpointed first, then the path is depicted along the solar surface by visually following the projection positions of the filament.Note that since the path is a projection of the 3D filament on the solar surface, it does not follow the filament body because of this projection effect (Figure 3a).The final 3D path is determined by assigning a height distribution to the projected path, which is a circular arc profile determined by the two footpoints and the apex height, 30 Mm, of the filament.We have to repeat the measurement by trial-and-error, which is 5 times in this case, and compare the constructed 3D flux rope with the observed filament to determine the validity of the selected path.To where µ is the permeability of the plasma.The plus and negative signs are determined by the helicity sign of the active region.If we stand on the positive polarity as shown in Figure 3b, the axial field of the flux rope points to the right hand side.Therefore, the chirality of the axial field is dextral, which implies a negative helicity.Then, we select the negative sign in Equation (6).

Potential Field Source Surface Model
The

Nonlinear Force-Free Field Model
With the previous two steps, a RBSL flux rope and a PFSS model are derived.We add them together, or equivalently, embed the RBSL flux rope into the PFSS model.This combined magnetic field recovers the radial magnetic field on the photosphere at r = 1R ⊙ , but it is not in an exact numerical force-balanced equilibrium state.Indeed, we find that the force-free and divergence free metrics (Wheatland et al. 2000) are σ J = 0.38 and ⟨|f i |⟩ = 2.4 × 10 −5 , respectively.Using the magneto-frictional method built in MPI-AMRVAC (Guo et al. 2016a,b), we relax the magnetic field.
The bottom boundary is the preprocessed vector magnetic field B r , B θ , and B ϕ on the photosphere, which are prescribed in the inner ghost layers as the boundary conditions.After iteration of 1.8 × 10 5 steps, the force-free and divergence-free metrics decrease to σ J = 0.32 and ⟨|f i |⟩ = 1.2 × 10 −5 .We note that σ J = 0.51 and ⟨|f i |⟩ = 9.9 × 10 −4 in Guo et al. (2016a).Both metrics are improved in this case.This relaxed NLFFF model serves as the initial magnetic field for the following zero-β MHD simulations.

Zero-β MHD Model
The zero-β MHD model considers three equations, namely, the mass conservation, momentum conservation, and magnetic induction equations (Guo et al. 2019a(Guo et al. , 2021b)): where ρ, v, and B are the density, velocity, and magnetic field, respectively, η is the resistivity and the electric current J = ∇ × B in the dimensionless form.In the momentum conservation equation, only the Lorentz force is included to change the momentum, while no gas pressure gradient and other external forces (e.g., gravity) are considered.The energy equation is neglected.This model could approximate the evolution of magnetic field in low-β circumstances, such as those in the corona, although it neglects some important MHD processes.For example, no slow-mode MHD waves are incorporated, and dynamical aspects like Petschek-type reconnection involving slow-mode shocks can never be resolved.The computation domain is selected as pole and is the complementary angle of the latitude, and ϕ is the longitude measured from the central meridian on the backside of the Sun.The computation domain is resolved by a uniform mesh of 360 × 420 × 400 cells.
The initial condition for the magnetic field is the NLFFF model relaxed from the RBSL and PFSS combined magnetic field at 03:00 UT.The initial condition for the density is provided by a stratified atmosphere model, which is very similar to that in Guo et al. (2019a).Namely, we first prescribe a step-like function for the temperature T (r): where Then, we compute the density distribution ρ strat (r) by: with p = ρT in the dimensionless form.The density at 1R ⊙ is assumed to be ρ(R ⊙ ) = 1.0 × 10 8 ρ 0 , where ρ 0 = 2.34×10 −15 g cm −3 .Similar to the simulation in Guo et al. (2021b), in order to mimic the slow evolution stage due to the heavy mass of the filament, the initial density at 03:00 UT is increased to ρ fil = 1.0 × 10 5 ρ 0 = 2.34 × 10 −10 g cm −3 wherever the density is smaller than this value, namely, ρ(r) = max(ρ strat (r), ρ fil ).The density ρ fil is used to mimic the inertia of the filament material.
Although the mass density of filaments still has large range of values due to the uncertainties in measuring hydrogen ionization degree (Labrosse et al. 2010), a typical range is 10 −14 -10 −12 g cm −3 (Hirayama 1986;Heinzel et al. 1996).The mass density ρ fil is about two orders of magnitude larger than the upper limit of the density range.This is necessary in the zero-β model because there is no gravity in the momentum equation.In the second stage between 03:34 UT to 03:44 UT, the initial density at 03:34 UT is reset back to the original density profile ρ strat (r) with the minimum density at the coronal top down to 3.4ρ 0 , such a manipulation is aimed to simulate the drainage of filament material.Without the density reset, the slow phase and fast eruption phase of the filament eruption cannot be reproduced with the zero-β model.Finally, we note that the initial velocity is always zero everywhere in the computation box.
The boundary condition is a data-constrained case as in Guo et al. (2019a).The density on all six boundaries is provided by the initial values.The velocity on all six boundaries is zero.The magnetic field in the inner ghost layer on the bottom boundary is the observed vector magnetic field B r , B θ , and B ϕ at 03:00 UT, and it is provided by a zero-gradient extrapolation for all the other boundaries.
To control the magnetic reconnection and the evolution process of the flux rope, we adopt an artificial resistivity in the first two minutes of the simulation (03:00-03:02 UT), where η = η 0 [(J − J c )/J c ] 2 wherever J ⩾ J c .We set η 0 = 8.1 × 10 12 cm 2 s −1 and J c = 1.1 × 10 −8 A cm −2 .At other times, there is no explicit resistivity, so only numerical dissipation is present.The resistivity is needed at the beginning to decrease the electric current, so the flux rope does not writhe too much during the eruption.It is turned off at later times to keep the flux rope coherent and not to be dissipated away.Finally, we adopt a three-step time integration, HLL Riemamn solver and Koren limiter in MPI-AMRVAC to solve the zero-β MHD equations.

RESULTS
Figure 4 and the attached animation show the evolution of magnetic field lines in the MHD simulation.A magnetic flux rope relaxed from the RBSL and PFSS combined magnetic field suspends in the corona as shown in Figure 4a.The south footpoint is rooted in the major negative polarity, and its north footpoint in a quiescent region with a weak positive polarity (see also Figure 3b).If we look from a side view as shown in Figure 4b To compare the MHD model with observations, we overlay the magnetic field lines on SDO/AIA 193 Å images as shown in Figure 5. Since the MHD model is computed in the spherical coordinate system, the alignment between it and SDO/AIA observations is straightforward.We put both the MHD and observational data in the heliocentric coordinate system, O-xyz.They are aligned in Figure 5.At 03:00 UT, the magnetic field is our initial condition and the shape of the magnetic flux rope coincides with the filament well, since the path is determined by the filament itself.The good alignment also demonstrates that our method to determine the flux rope path as described  in Section 3.1 is reliable.At 03:35 UT, the filament starts to erupt, so it shows a bright thread as indicated by the arrow in Figure 5d.The magnetic flux rope in the MHD model also starts to accelerate (Figure 5c).Filament matter drains along both the northern and southern legs to their footpoints.At 03:42 UT, the magnetic flux rope has risen outside the field of view (Figure 5e), and only the filament material in the northern leg could be seen (Figure 5f).
Quasi-separatrix layers (QSLs) are 3D thin volumes indicating a drastic change of magnetic field line linkage (Priest & Démoulin 1995;Démoulin et al. 1996).QSLs are potential regions for electric current accumulation and magnetic reconnection.Titov et al. (2002) proposed the squashing degree Q to locate QSLs, which are regions where Q ≫ 2. Tassev & Savcheva (2017) and Scott et al. (2017) proposed new methods using transported displacement vectors or tangents along field lines to calculate the squashing degree Q. Kai E. Yang implemented a QSL locator code (K-QSL) using the formulation of Scott et al. (2017), which is released on GitHub1 .One advantage of this code is that it can also be applied to a spherical coordinate system.Therefore, we apply this QSL locator code to our MHD model.7b).The slices are aligned with each other.Then, we measure the positions of the filament and flux rope at different times for both the SDO/AIA 171 Å data series and MHD simulation snapshots.We  phases in the evolution, including a slow rising phase and a fast erupting phase.The velocity reaches 200.6 km s −1 in the fast erupting phase of the MHD simulation.The evolution in the fast erupting phase of the observation shows a slightly different transition, where the velocity is lower than the simulation at ∼03:36 UT but higher than it at ∼03:44 UT.The average twist of the flux rope is −3.4 turns at 03:02 UT, which is computed by the formula in Berger & Prior (2006).We think this magnetic structure is feasible for the studied case, since multiwavelength emission is usually different from magnetic field structure as noticed by other researchers (e.g., Fan 2017).High twist flux rope is also prone to appear in intermediate and quiescent prominence regions (e.g., Su et al. 2015;Mackay et al. 2020;Guo et al. 2021a).
There are other MHD simulations in a spherical wedge, where they either use more advanced unstructured adapted mesh (Amari et al. 2018), or sophisticated MHD equations including full thermodynamics (Fan 2017).Our model only uses uniform grid and zero-β model.Taking the full advantage of MPI-AMRVAC, we could extend the simulation using adaptive mesh refinement or stretched grids, and including more physics with full thermodynamic MHD equations.Besides, Amari et al. (K-QSL), based on the method proposed by Scott et al. (2017).The advantage of this code is that it can be applied to static magnetic field models (e.g., PFSS) and MHD simulations in both the Cartesian and spherical coordinate systems, as demonstrated here.
This paper illustrates that the RBSL model together with the PFSS, magneto-frictional method and MHD models in spherical coordinates provides a solid approach to simulate intermediate and quiescent filaments in a large field of view.With these techniques, we could extend these simulations to larger scales, all the way to the interplanetary space.

Figure 1 .
Figure 1.Evolution of the M8.7 flare on 2012 January 23 in different SDO/AIA EUV wavebands.Panels (a), (d) and (g) show the 171 Å images at 03:00, 03:35 and 03:42 UT, respectively.The red arrow in panel (a) points to the apex of the coronal loops.The green arrow points to the local radial direction projected on the plane of sky.Panels (b), (e) and (h) show the 304 Å images.Panels (c), (f) and (i) show composite images constructed by 211, 193 and 171 Å observations representing the red, green and blue channels, respectively.An animation is attached to this figure, which shows the high cadence evolution of the M8.7 flare in the same wavebands covering a larger time range from 03:00 to 05:00 UT on 2012 January 23.

Figure 2 .
Figure 2. Heliocentric coordinate system and the local coordinate system of the filament.Point O is located at the center of the Sun, the x-axis points to the central meridian on the Sun's backside and on the equator, the y-axis points to the east solar limb and on the equator, and the z-axis constitutes a right-hand system together with the x-and y-axes.O ′ is located at the middle point of the two footpoints (indicated by the two blue dots F1 and F2) of the filament on the photosphere, the X-axis is perpendicular to both

Figure 3 .
Figure 3. Path of the filament.(a) Filament in the SDO/AIA 304 Å image at 03:00 UT. Green diamonds represent the projected path of the filament body on the solar surface.(b) Projected path overlaid on SDO/HMI magnetic field at 03:00 UT. Green and yellow filled circles represent the two footpoints of the filament.Red and blue lines depict the projected path of the filament on the solar surface.
keep the current circuit closed, we set the mirror image of the observed flux rope path in O ′ -XY Z as the complementary circuit beneath the photosphere.Due to this setting, we use the RBSL kernels derived byTitov et al. (2018) to configure the magnetic flux rope model.Next, the minor radius is estimated by the width of the filament, which is 6 Mm.Then, the axial magnetic flux F 0 is estimated by the average of the absolute values of the two footpoints, |F + | and |F − |, as shown by the green and yellow filled circles in Figure 3b.It is found that F 0 = (|F + | + |F − |)/2 = 1.75 × 10 20 Mx.This flux can be considered as the axial flux of the rope only if the footpoint distance is short.In practice, it is often multiplied by a factor since the method does not require F = F 0 .In this case, the difference of the longitudinal flux at the footpoints of the flux rope can be compensated by adding the RBSL field and PFSS field together.We finally use F = 4F 0 to simulate the flux rope eruption since it can better match the evolution of the filament eruption compared to the numerical experiments with other values such as F = F 0 and F = 2F 0 .The electric current is computed by Equation (12) in Titov et al. (2018): RBSL flux rope constructed with the previous parameters (axis path C, minor radius a, magnetic flux F , and electric current I) is in an internal quasi-equilibrium state.It represents the magnetic field generated by local electric currents in the corona.But it does not include the background coronal magnetic field generated by magnetic polarities on the photosphere other than the flux rope itself.To simulate the magnetic field in the whole active region, we have to embed it in a large-scale magnetic field produced by photospheric magnetogram excluding the flux rope, which could be well approximated by a potential field.It is specifically modelled by the PFSS model in the spherical coordinate system.The PFSS model solves the Laplace equation, which governs the potential field, with the spherical harmonic expansion(Altschuler & Newkirk 1969;Schrijver & De Rosa 2003).The solution is represented by the sum of a series of polynomials with harmonic coefficients.The harmonic coefficients are computed by the radial magnetic field on the photosphere.It is usually derived from a synoptic map constructed with magnetic field observations in a full solar rotation period.Sometimes, a synoptic frame is modified by inserting a snapshot of magnetic field observations into the corresponding synoptic map to increase the temporal resolution on the visible solar disk.In this paper, we adopt the same synoptic frame as used inGuo et al. (2012), which is constructed by the synoptic map for Carrington rotation 2119 from 2012 January 9 to February 6 and the vector magnetic field observed by the Helioseismic and Magnetic Imager (HMI; Scherrer et al. 2012) aboard SDO.However, before using it as the boundary condition for the PFSS model, one more step is needed, i.e., to subtract the radial magnetic field produced by the RBSL flux rope on the two footpoints.According to the design of the RBSL method, the vertical magnetic field B Z on the local bottom plane O ′ -XY only has values at the two footpoints within minor radius a.Note that if the two footpoints are not separated too far, the RBSL radial magnetic field B r,RBSL (R ⊙ , θ, ϕ) at r = 1R ⊙ is also concentrated at the two footpoints.We have to subtract B r,RBSL (R ⊙ , θ, ϕ) from the synoptic frame to prepare the final boundary condition for the PFSS model, which is computed by the PFSS module in MPI-AMRVAC(Porth et al. 2014) throughout this paper.
, the flux rope displays a cavity morphology as often observed beyond the solar limb.As we reset the density distribution by decreasing the high density to a lower value at 03:34 UT, the flux rope starts to accelerate.The snapshots at 03:35 UT in Figures 4c and 4d are typical ones at the beginning of the acceleration.The body of the flux rope rises and expands radially outward.At 03:42 UT, the flux rope rises up to a certain height, close to the top boundary of the computation box.The flux rope axis writhes a little due to its large twist.

Figure 4 .
Figure 4. Evolution of magnetic field lines in the MHD simulation.Panels (a), (c) and (e) show the SDO view at 03:00, 03:35 and 03:42 UT, respectively.Panels (b), (d) and (f) show a side view from the south end of the flux rope to the north end.An animation is attached to this figure, showing the evolution of the magnetic field in the time range of 03:00-03:44 UT.The time cadence is 2 minutes between 03:00-03:34 UT, and it is 12 seconds between 03:34-03:44 UT.

Figure 5 .
Figure 5. MHD model overlaid on SDO/AIA 193 Å observations.(a, c and e) Background images show the 193 Å observations at 03:00, 03:35 and 03:42 UT from top to bottom panels.Colored lines represent the magnetic field lines in the MHD model.(b, d and f) SDO/AIA 193 Å images highlighting the erupting filament indicated by the cyan arrows at 03:00, 03:35 and 03:42 UT from top to bottom panels.
Figure 6 and the attached animation show the distribution of the logarithm of the squashing degree (log Q).We find that the erupting flux rope stretches overlying field lines, and new QSLs on the solar surface form during this process as shown in the closed orange curve in Figure 6g.These new QSLs coincide with the flare ribbons observed by the SDO/AIA 1600 Å waveband.The flux rope is wrapped by QSLs, which is clearly revealed by the log Q distribution on a θ = 50.28• plane in Figures 6c, 6f, and 6i.The evolution of log Q on this vertical plane displays the eruption process and expansion of the flux rope.We compare the MHD simulation quantitatively with observations by tracking the kinematics of the erupting flux rope and filament.Figure7displays the time-distance profiles for them.We first select a slice vertically in the SDO/AIA 171 Å images and in our MHD simulation (Figures7a and

Figure 6 .
Figure 6.QSL overlaid on SDO/AIA 1600 Å images and SDO/HMI magnetic field.(a, d and g) Distribution of the logarithm of the squashing degree (log Q) on the solar surface overlaid on the SDO/AIA 1600 Å images at 03:00, 03:35 and 03:42 UT from top to bottom panels.The closed orange curve in panel (g) indicates the region where QSLs coincide with flare ribbons.(b, e and h) Distribution of log Q on the solar surface overlaid on the radial magnetic field observed by SDO/HMI.(c, f and i) Distribution of log Q on the θ = 50.28• plane.An animation is attached to this figure, showing the evolution of the squashing degree log Q in the time range of 03:00-03:44 UT.The time cadence is 2 minutes between 03:00-03:34 UT, and it is 12 seconds between 03:34-03:44 UT.

Figure 7 .
Figure 7. Kinematics of the erupting filament in observation and flux rope in MHD simulation.(a) The purple line over-plotted on the SDO/AIA 171 Å image at 03:00 UT indicates the slice that we choose to measure the time-distance profile of the erupting filament.Red dots are one example of measurements of the filament positions from 03:00-03:44 UT.(b) The purple line over-plotted on the MHD simulation at 03:24 UT indicates the slice that we choose to measure the velocity of the erupting magnetic flux rope.Yellow dots are one example of measurements of the flux rope positions from 03:00-03:44 UT.(c) Red dots and vertical line segments represent the positions and errors measured from the observation of SDO/AIA 171 Å waveband.Yellow dots and vertical line segments represent the positions and errors measured from the MHD simulation.The blue and cyan lines are linear fittings to the time-distance profile of the filament and flux rope at the fast erupting stage, which lead to speeds of 521.8 km s −1 and 200.6 km s −1 , respectively.

5.
SUMMARY AND DISCUSSION We implemented and applied a RBSL model in the spherical coordinate system, embeded it in the PFSS model, and relaxed the combined model by the magneto-frictional method in MPI-AMRVAC.The relaxed model resembles the SDO/AIA multi-wavelength observations very well, including the intermediate filament and coronal loops.Using the relaxed model as the initial condition, we performed a zero-β MHD simulation for the M8.7 flare on 2012 January 23.The simulation also reproduced the intermediate filament eruption well.The general morphology and eruption direction of the simulation are similar to those in the SDO/AIA observations.The QSLs on the bottom surface display separating flare ribbons, and the QSLs on a vertical surface with θ = 50.28• show the rising and expanding flux rope.The height-time profile measured from the MHD simulation is close to that from the SDO/AIA 171 Å observation, which indicates that the MHD model reproduces the kinematics of the observed filament assuming the magnetic configuration is a flux rope to a good extent.

(
2018) studied an active region with strong magnetic field, which can be constructed directly by an NLFFF extrapolation.Fan (2017) used a fully parameterized magnetic flux rope model, which is not constrained by observations.The model developed in this work can be applied to intermediate or quiescent filaments with weak magnetic field on the photosphere.Meanwhile, the magnetic field model is constrained by observations.When comparing simulations with observations, several features, such as the flux rope and flare ribbons, are related to QSLs, which are quantitatively parameterized by the squashing factor Q.Various QSL calculation codes have been developed.Recently,Zhang et al. (2022) compared different QSL codes, all in the Cartesian coordinate system.Locating QSLs in spherical settings has only been done in very few cases, such asTitov et al. (2011), who used a generalized definition of the squashing factor Q proposed byTitov (2007).In order to calculate the Q distribution in the spherical coordinates, one of the coauthors of this paper developed an open-source, spherical QSL locator code