Numerical Investigations of Catastrophe in Coronal Magnetic Configuration Triggered by Newly Emerging Flux

We performed 2D magnetohydrodynamical numerical experiments to study the response of the coronal magnetic configuration to the newly emerging magnetic flux. The configuration includes an electric-current-carrying flux rope modeling the prominence floating in the corona and the background magnetic field produced by two separated magnetic dipoles embedded in the photosphere. Parameters for one dipole are fixed in space and time to model the quiet background, and those for another one are time dependent to model the new flux. These numerical experiments duplicate important results of the analytic solution but also reveal new results. Unlike previous works, the configuration here possesses no symmetry, and the flux rope could move in any direction. The non-force-free environment causes the deviation of the flux rope equilibrium in the experiments from that determined in the analytic solution. As the flux rope radius decreases, the equilibrium could be found, and it evolves quasi-statically until the flux rope reaches the critical location at which the catastrophe occurs. As the radius increases, no equilibrium exists at all. During the catastrophe, two current sheets form in different ways. One forms as the surrounding closed magnetic field is stretched by the catastrophe, and another one forms as the flux rope squeezes the magnetic field nearby. Although reconnection happens in both the current sheets, it erases the first one quickly and enhances the second simultaneously. These results indicate the occurrence of the catastrophe in asymmetric and non-force-free environment, and the non-radial motion of the flux rope following the catastrophe.


Introduction
The prominences floating in the solar corona contain cool and dense plasma (Mackay & van Ballegooijen 2009). Coronal mass ejections (CMEs)-the most violent eruptive process on the Sun-are considered strongly related to the eruptive prominence (Gopalswamy et al. 2003;Jing et al. 2004). Revealing the nature of the magnetic structure surrounding eruptive prominence is therefore important for studying observational behaviors of CMEs. Three well-accepted models of the magnetic configuration process related to the eruptive prominence were proposed, namely, the sheared arcade model (e.g., Kippenhahn & Schlüter 1957;Antiochos et al. 1994;Mikić & Linker 1994;DeVore & Antiochos 2000), the flux rope model (e.g., Kuperus & Raadu 1974;Lin & Forbes 2000;Gibson et al. 2004;Amari et al. 2014;Liu 2020), and the breakout model (Antiochos et al. 1999).
Since the coronal magnetic configuration that includes the prominence anchors in the photosphere with their footpoints, the unceasing motion of the photospheric plasma brings the footpoints of the coronal configuration to move together continually, changing the state of the equilibrium in the configuration successively. Eventually, the configuration loses its mechanical equilibrium and the eruption is invoked. The transition from the equilibrium to the loss of equilibrium constitutes the catastrophe. This process consists of two stages. First, the system accumulates the energy as the magnetic structure in the corona is deformed owing to the motion of the footpoints, and then the excess magnetic energy is released as the loss of equilibrium in the system takes place. According to the catastrophe model of the solar eruption, the loss of the equilibrium in the system thrusts the flux rope that is used to model the prominence rapidly outward, stretching severely the magnetic field lines passing around the flux rope with their two ends anchored in the photosphere. A current sheet develops below the flux rope, separating magnetic fields of opposite polarity; magnetic reconnection takes place in the current sheet and helps the flux rope escape from the Sun smoothly (e.g., Lin & Forbes 2000;Lin 2002;Lin et al. 2003). In the framework of this model, the magnetic configuration possesses symmetry with respect to the direction perpendicular to the solar surface (also known as the radial direction) either before or after losing equilibrium.
In reality, on the other hand, the eruptive prominence does not necessarily escape ideally in the radial direction of the Sun. Gopalswamy et al. (2003) statistically studied 186 events of the eruptive prominence and classified them into radial events and transverse ones depending on the ejecta trajectory. Here, for radial events the ejecta moves radially, and for the transverse ones the motion is roughly tangential to the solar limb. The Gopalswamy et al. (2003) study found that 82% of the events were radial and 18% were transverse. McCauley et al. (2015) studied 904 events of high spatial resolution and found that 75% of events were radial, 12% were oblique, 11% were transverse, and the remaining were ambiguous. These results suggest that the motion patterns of the ejecta in reality are apparently more complex than those described in a simple model given by Lin & Forbes (2000), Lin (2002), Lin (2004), Lin & Soon (2004), and .
When looking into the triggering of the catastrophe due to the newly emerging magnetic flux, Lin et al. (2001) suggested that the asymmetry in the magnetic configuration might be able to account for the complex motion patterns of the eruptive prominences and/or CMEs, although their model could be considered simple. For the system that they studied, on the other hand, Lin et al. (2001) could not determine whether the loss of equilibrium of system would develop to a successful eruption, or what the motion of the flux rope would look like. This is because the consequent evolution following the catastrophe in an asymmetric configuration is already beyond the scope of the analytic solution, and a numerical approach should be considered instead.
On the basis of the work of Lin et al. (2001), we performed numerical experiments in the present work to look into the detailed evolution in the coronal magnetic configurations driven by the newly emerging magnetic flux. An electriccurrent-carrying flux rope is included in the configuration. We present the quasi-static evolution in the system, and we also briefly discuss the dynamic behavior of the system after the catastrophe. In the next section, we introduce the magnetic configuration constructed by Lin et al. (2001) and the numerical method used in this work. In Section 3 we present the numerical results, in Section 4 we give a discussion about this work, and in Section 5 we summarize this work and give our conclusions.

Descriptions of Numerical Approach
The dynamical evolution of the physical system being investigated is governed by the magnetohydrodynamics (MHD) equations given below: where ρ, v, e, p, B, and, g are mass density, plasma velocity, total energy density, gas pressure, magnetic field, and gravity, respectively. The total energy e = ε + ρv 2 /2 + B 2 /2μ includes the thermal energy ε, the kinetic energy, and the magnetic energy; τ = ν d [∇v + (∇v) T − 2(∇ · v)I/3] is the stress tensor, where ν d = ρ ν k is the dynamic viscosity and ν k = 10 7 m 2 s −1 is the kinematic viscosity. Assuming that the system of interest follows the ideal gas law and that the plasma is the completely ionized hydrogen plasma, we then have p = (γ − 1)ε, with γ = 5/3 being the ratio of the specific heats. The physical magnetic diffusivity η = 5 × 10 7 m 2 s −1 is set for the whole simulation domain for simplicity. This value of η is apparently large compared to that of the true corona in order to match the grid size used for the calculations in the present work. For the large-scale activities occurring in the corona, like eruptive prominences and CMEs, the values of the Alfvén speed, v A , and the length scale, L 0 , are typically 10 3 km s −1 and 10 5 km, respectively, which leads to the magnetic Reynolds number R m = v A L 0 /η = 2 × 10 4 for the above η value. To justify this result, we estimate the effective R m after considering the impact of the grid size as below.
In the true corona, because of the low diffusivity and the large scale, both the Reynolds number, R e = UL/ν k , and R m possess huge values, say, up to 10 12 . In numerical simulations, on the other hand, the values of R e and R m cannot reach such a high value owing to the finite size of the grid. This is because any fine structure of size less than the grid size is inevitably diffused, which is equivalent to introducing an extra (artificial) diffusive process into calculations. Therefore, the effective value of R e and/or R m in simulations is somehow governed by the grid size. In MHD simulations, the impact of the finite size on the effective R m could be evaluated in the following way.
Normalizing η to v A L 0 gives R m = 1/η. According to Forbes & Priest (1983), the grid size δx should be less than the Kolmogorov dissipation scale l k = (η 3 /ò) 1/4 , with ò being the rate of the energy dissipation. All parameters here are dimensionless. For the system with given η, l k depends on ò, which reaches maximum, say, 0.1, in the fast reconnection process and yields the smallest l k . Therefore, to avoid apparent numerical diffusion, δx max = (10η 3 ) 1/4 . Conversely, if δx is selected, then η related to the numerical diffusion is consequently determined, η = (0.1δx 4 ) 1/3 , which yields the effective R m = (0.1δx 4 ) −1/3 . Thus, for the grid used in the present work, δx = 10 −3 , we have the effective R m = 2.15 × 10 4 , which is consistent with the value of R m associated with η mentioned above. The effective R e could be estimated in the same way.

The Initial Configuration
We start with the magnetic configuration that was analytically derived by Lin et al. (2001). The analytic solution provides the information of the evolution in the equilibrium location of the flux rope in response to the gradual change in the bottom boundary, which is used to model the surface of the photosphere. The flux rope is used to model the filament or the prominence floating in the corona. Interested readers may refer to Lin et al. (2001) for details of the analytical description of the related magnetic configuration that is used as the initial condition in the present work. Figure 1 displays the sketch of the initial magnetic configuration. The background magnetic field is produced by two magnetic line dipoles beneath the photosphere, accounting for the preexisting magnetic field (namely, Dipole 1) and the newly emerging magnetic flux (Dipole 2). In the present work, Dipole 1 remains unchanged and provides an unchanged background field, and Dipole 2 changes its location or strength to model the newly emerging magnetic field. To describe the initial magnetic structure easily, we duplicate some works of Lin et al. (2001) and write down the total magnetic field as follows:  where M and (0, −d), S and (x d , −y d ), J and (x h , y h ), and − J and (x h , − y h ) represent the dimensionless strengths and the positions of Dipole 1, Dipole 2, the center of the electric current inside the flux rope, and its mirror image, with d fixed at d = 4 × 10 7 m. We note here that the magnetic field governed by Equations (6) and (7) is valid for the magnetic configuration with a thin flux rope required for the validity of the analytic solution. For the case studied here, the flux rope possesses finite radius r, and a smooth distribution of the electric current should be included inside the flux rope in a more realistic fashion. Following the strategy of Forbes (1990), we take the profile of the electric current density j z inside the flux rope as where j 0 has the dimension of electric current density, J is the electric current intensity inside the flux rope in units of I 0 , I 0 = j 0 [π(r 2 + Δ 2 /4 − 2Δ 2 /π)], j 0 J is the electric current density at the center of the flux rope, Δ = r/5 is the width of the thin shell around the edge of flux rope, and To perform numerical experiments, we rewrite the magnetic field governed by Equations (6) and (7) as follows by including the electric current distribution given by Equations (8)-(10):

Gravitationally Stratified Atmosphere
The 2D simulation in this work is performed in the region (− 5L 0 x 7L 0 , 0 y 12L 0 ) with L 0 = 10 8 m. The bottom of the photosphere lies on y = 0. The gravitationally stratified atmosphere is divided into three layers: the photosphere (0 y h p ), the chromosphere (h p y h c ), and the corona (h c y 12L 0 ), where h p = h c = 10 6 m. We set the temperature of the photosphere as T p = 4300 K and that at the bottom of the corona as T c = 1.06 × 10 6 K at y = h p + h c . The temperature distribution in the chromosphere via a simple interpolation follows the practice of Mei et al. (2012a): Combining with the hydrostatic equilibrium equation, and the equation of state for the ideal gas, we find the pressure and the density distribution in the photosphere and the chromosphere, where k = 1.38 × 10 −23 J K −1 is the Boltzmann constant, m p = 1.67 × 10 −27 kg is the proton mass, andˆ( ) = -+ g g y y R 1 0 2  , with R e = 6.961 × 10 8 m and the gravity near the solar surface g 0 = 274 m s −2 . Hereŷ is the unit vector in the y-direction. For the extended corona, we describe the plasma density through the S&G model (Sittler & Guhathakurta 1999) in the form of ρ(y) = ρ 0 f (y), where ρ 0 = 1.67 × 10 −11 kg m −3 is the plasma density at the bottom of the corona and f (y) is given as where z(y) = 1/(1 + y/R e ), a 1 = 0.001292, a 2 = 4.8039, a 3 = 0.29696, a 4 = −7.1743, and a 5 = 12.321. Then, we have the initial temperature and density distribution in the whole domain by using Equations (15)-(17). In our simulations, the temperature in the flux rope is T f = 5 × 10 4 K, and that in the ambient plasma is T amb = 2 × 10 6 K. We assume that the cold material inside the flux rope comes from the chromosphere (Mackay & van Ballegooijen 2009) and that the flux rope includes a thin and smooth layer for temperature: The pressure inside the flux rope p f is calculated through

Boundary Conditions
The boundary conditions in our simulation are as follows. On the three boundaries (x = −5L 0 , x = 7L 0 , and y = 12L 0 ), the simple extrapolation method is utilized to realize the outflow boundary condition. Specifically, the zero-gradient boundary condition is used to extrapolate the density, magnetic field, and pressure, and we force the plasma at the aforementioned three boundaries to flow outside. We note here that the extrapolation of the magnetic field in this paper is realized via the magnetic vector potential, A, instead of magnetic field, B. According to Ye et al. (2019), if the B-field is extrapolated, the motion of the flux rope, especially in the eruptive stage, would be eventually ceased at the upper boundary as a result of the boundary effect. If the A-field is used instead, the flux rope could pass through the top boundary smoothly. This indicates that utilizing extrapolation of the A-field would help elude the boundary effect.
For the bottom boundary, on the other hand, we follow the practice of Mei et al. (2012a) to realize the line-tied condition for the photospheric magnetic field. We fix all the parameters at the bottom and utilize two dense layers to approximate the linetied condition described by Forbes (1990) and Shen et al. (2011). The boundary conditions are thus given as follows: where p c = 2ρ 0 kT c /m p = 0.29 Pa is the pressure at y = h p + h c The velocity at the bottom follows the no-slip and the wall conditions. As discussed in Mei et al. (2012a), it is a challenge to realize the line-tied condition with small plasma β at the bottom boundary. Here two dense layers yield the plasma β = 900 at y = 0, and the motions of the magnetic field lines are mainly dominated by the bulk flow of the no-slip plasma. Consequently, the magnetic field lines are rooted at the bottom and the line-tied condition can be handled.

Results of Experiments
We are now ready to perform numerical experiments to look into the evolution of the location of the flux rope in response to the change in the background field that is governed by the behavior of Dipole 2. A root grid of dimension 480 × 480 spans the simulation domain, with the largest grid size being δx = δy = 0.025L 0 . Static mesh refinement is used in two regions: the bottom layer that possesses a sharp gradient in density, and the region above that layer that covers the rest of the simulation domain. Six levels of the refinement are performed in the bottom layer, such that refinement of level 6 is performed for the region y 0.01L 0 , that of level 5 for 0.01L 0 < y 0.02L 0 , that of level 4 for 0.02L 0 < y 0.03L 0 , that of level 3 for 0.03L 0 < y 0.04L 0 , that of level 2 for 0.04L 0 < y 0.06L 0 , and that of level 1 for 0.06L 0 < y 0.1L 0 . In addition, three rectangle regions including the flux rope are selected as follows: the level 4 refinement is conducted for region 3.1. Impact of the Initial Radius of the Flux Rope on the System Evolution The NIRVANA code 3.8 (Ziegler 2004(Ziegler , 2005(Ziegler , 2008(Ziegler , 2011) is used to carry out the calculation for solving MHD Equations (1)−(5) associated with Equations (6)−(27). Accordingly, the Godunov scheme, HLL Riemann solver, and mesh refinement technique are turned on in calculations to capture fine structures in the magnetic reconnection region.
According to Lin et al. (2001), the equilibrium curve of the flux rope is apparently affected by the initial radius r 00 . To investigate the dependence of the system evolution on the initial radius, we perform several simulations for r 00 = 0.01d, r 00 = 0.04d, and r 00 = 0.08d, respectively, while keeping the other parameters fixed. The initial location of the flux rope and the related parameters are taken from the equilibrium curve given by analytic solutions (see, e.g., Lin et al. 2001). For each value of r 00 , we check the change in the flux rope location from its initial location for S taking values of 1, 2, 3, 4, 5, and 6, respectively. As expected, behaviors of the flux rope vary from case to case. For the case of r 00 = 0.01d, the flux rope quickly leaves the initial locations in the simulations and soon finds its new equilibrium near the initial location. For the case of r 00 = 0.04d, the flux rope duplicates its behavior in the case of r 00 = 0.01d when S takes values from 1 to 5, but the flux rope cannot find an equilibrium position, as it leaves the initial location when S = 6. For any case of r 00 = 0.08d, the flux rope can never find equilibrium after leaving the initial location. Therefore, in the rest of our work, we shall not look into the case of r 00 = 0.08d since no equilibrium configuration could be reached at least for S taking the values of our interest. Here we just show one simulation with r 00 = 0.08d and S = 3 for illustration (see Figure 2). The relevant results with r 00 = 0.01d and r 00 = 0.04d will be discussed in the following sections. Figure 3 displays the location (x h , y h ) of the flux rope in equilibrium. Solid curves in Figures 3(a) and 3(b) are for the analytic results, and the separated circles are for the numerical results when r 00 = 0.01d and 0.04d, respectively. As we mentioned before, the flux rope could always find equilibrium after leaving the initial location if r 00 = 0.01d, and no equilibrium could be found if r 00 = 0.04d and S = 6. So Figures 3(a) and 3(b) include six circles for r 00 = 0.01d and five circles for r 00 = 0.04d. This indicates that the magnetic configuration of interest could remain equilibrium, as the strength of the newly emerging flux, S, varies in a large range in the case of thin flux rope, and that when the configuration includes a flux rope of medium size, it still can find equilibrium as S varies within a certain range, but no equilibrium could be reached as S goes beyond this range. Considering the fact that the coronal magnetic configuration evolves slowly prior to the eruption in response to the gradual variation in the photosphere, we conclude that the behavior of the flux rope in the case of r 00 = 0.01d displayed in Figures 3(a) and 3(b) could describe the evolutional feature of a true prominence. Furthermore, the behavior of the flux rope in the case of r 00 = 0.04d suggests that the newly emerging flux could indeed drive the coronal magnetic structure to the catastrophe point through a set of equilibrium states in a quasi-static fashion.
We note here that the deviation of the numerical results from the analytic ones is apparent. This mainly results from the existence of the gravity in the present work. We point out here that the gravity impacts all aspects of the system, including stratification of the atmosphere, shape of the flux rope, and the way the Alfvén speed varies with the altitude. Previous works either did not consider the gravity at all (see, e.g., Lin & Forbes 2000), or just considered the gravity acting on the flux rope only (see, e.g., Reeves & Forbes 2005), or considered the atmosphere stratified by the gravity but did not include the gravity in the force acting on the flux rope (see, e.g., Lin 2002;, or considered the impact of the gravity on both the atmosphere and the flux rope independently (Lin 2004). Therefore, the impact of the gravity was just considered by these works in an approximate fashion. In the present work, obviously, the impact of the gravity on the whole system is considered comprehensively, and the deviation of the results obtained here from previous ones is thus expected.
In order to quantitatively describe the gravity effect on the evolution, we follow the methodology of Reeves & Forbes (2005) by imposing the gravity on the flux rope, and we rewrite the equation for the global equilibrium in the y-direction given by Lin et al. (2001): being the ratio of the gravity f g = mg 0 to the Lorenz force f b = b 0 I 0 /c, and b 0 = I 0 /(cd) = 2.2 G. As an example, we select one point with r 00 = 0.04d and S = 5 to perform the test because this point appears to be far from the analytic equilibrium curve. Combining Equation (28) and the flux-frozen condition at the surface of the flux rope given by Lin et al. (2001), we get the new equilibrium position with the effect of gravity: (x hg = 0.50, y hg = 1.87). Comparing with the result of Lin et al. (2001), in which case no gravity is included, (x h0 = 0.31, y h0 = 1.46), one can easily see the difference. Interesting enough, for the same setting of the other parameters, the numerical experiments give (x h = 0.58, y h = 2.08), which is consistent with the analytic result determined by Equation (28).
The red triangles in Figure 3 indicate the equilibrium positions deduced analytically as the gravity is included. Apparently, analytic and numerical approaches yield the same result at least for the parameters that we used here. The deviation is due to ignoring gas pressure in the analytic solution. Therefore, the equilibrium positions of the flux rope obtained here differ from those deduced by Lin et al. (2001) apparently because of the existence of the gravity.
In addition to the impact of the non-force-free environment occurring in the numerical experiment, the asymmetric feature of the configuration plays an important role in yielding the deviation. Xie et al. (2017) performed a similar experiment for the evolution in a symmetric system. Their results did show deviation from the analytic ones, but the symmetry allows the impact of various physical parameters to compensate one another easily, so that the deviation is not very large. In the case of asymmetry, on the other hand, the system possesses more dimensions in which changes may happen, relations among these parameters to one another become much more complicated, and the compensations of the impact of them on the system behavior get more difficult simultaneously. This directly causes the numerical results to deviate from the analytic ones apparently. Fortunately, the present numerical experiments show that the coronal magnetic configurations could find stable equilibrium states under the influence of the newly emerging flux and that the catastrophe could occur in some specific situations as well. This allows us to investigate the dynamic evolution in the system after the catastrophe in a self-consistent fashion.
Before we turn to the dynamic episode, we realize that Figure 3 reveals more information of interest. First of all, positions of the flux rope in equilibrium manifest a similar evolution trend in both the x-and y-direction for both cases of r 00 = 0.01d and r 00 = 0.04d (Figures 3(a) and 3(b)), and the corresponding radius of the flux rope displays similar behavior for the two cases (Figure 3(c)). However, the electric currents inside the flux rope in these two cases show different trends such that J decreases much faster in the case of small radius (Figure 3(d)). Looking carefully into the detailed behavior of the flux rope interior, we realized that it is the numerical diffusion that should account for the evolutionary feature of J in the case of r 00 = 0.01d. This is because the region where the flux rope is located includes fewer grid lines as the radius of the flux rope gets smaller for the given grid size. In our case, about 100 grid lines are included as r 00 = 0.04d, and only about 10 grid lines are included as r 00 = 0.01d. The impact of the numerical errors on the case of small radius is apparent and is presented by the unexpected decrease in J versus S, which indicates the extra diffusion of the electric current inside the flux rope by the numerical resistivity.
We also noticed that, on the other hand, the global evolutionary behavior of the flux rope is not affected apparently by its internal behavior. Such a behavior of the system has been pointed out and discussed by Lin et al. (2001). The internal distribution and change of the electric current may alter the magnetic field around the flux rope, but the interaction between the magnetic field produced by the flux rope and that produced by the other sources allows the flux rope to find the new equilibrium via adjusting the flux rope position and the internal current in a certain range. This confirms that the global behavior is the main issue that dominates the overall evolutionary property of the system.

Internal Equilibrium of the Flux Rope
The equilibrium in the system is reached as the net force acting on the flux rope vanishes. The total force includes both the internal and the external forces. In the analytic framework, the radius of the flux rope is small compared to the global scale of the system, so that the external magnetic field could be regarded invariant on the cross section of the flux rope in the evolution; the equilibrium can be divided into two parts, the external (global) equilibrium and the internal (local) equilibrium; and the two equilibria could be realized separately (see, e.g., Forbes & Isenberg 1991;Isenberg et al. 1993; Forbes et al. 1994). Consequently, the global equilibrium yields the dependence of the flux rope location on the background field, and the local equilibrium relates the radius of the flux rope to the total electric current intensity inside the flux rope. The two equilibria connect to each other via J since the magnetic field associated with the flux rope and that produced by the image of the flux rope are governed by J. For the force-free case, Isenberg et al. (1993) improved the solution given by Parker (1974) for describing the internal equilibrium of the flux rope and presented the accurate relation of the flux rope radius to J. But the expression is implicit and is not easy to use. Lin et al. (1998) obtained a simplified relation as follows: ( ) = r r J, 2 9 0 0 0 which turns out to be a good approximation to the exact solution of Isenberg et al. (1993). In the numerical experiment, the requirement for the small radius of the flux rope to decouple the two equilibria is lifted, and it is not necessary to decouple them any longer. The dependence of r 0 on J could be obtained directly from the simulation results. Figure 4 plots variations of r 0 versus J as S varies from 1 to 6 for r 00 = 0.01d and from 1 to 5 for r 00 = 0.04d. In addition, the results given by Equation (29) are also plotted (solid curves) for comparison. We notice that both analytic and numerical results show that r 0 decreases with J, but values of r 0 obtained numerically are always larger than those obtained analytically. This is because the requirement of small radius of the flux rope in the analytic solution is lifted in the numerical solution, and r 0 could be any value determined by both internal and external conditions in the simulation as a result of the self-adjustment in both the global and local magnetic fields. However, the selfadjustment does not always work, and the catastrophe in the system occurs as long as the self-adjustment fails. For example, in the case of r 00 = 0.01d, it works well in the range of S that we selected; in the case of r 00 = 0.04d, on the other hand, the  self-adjustment works well at S 5 and fails at S 6 (more accurate calculations indicate that the turning (critical) point at which the self-adjustment fails occurs at S = 5.75); and no equilibrium exists in the system at all for r 00 = 0.08d, where the function of the self-adjustment to help keep the system in equilibrium vanishes completely. This seems to suggest that the coronal magnetic configuration including a thick flux rope (prominence) will lose the equilibrium more easily triggered by the newly emerging flux than that including a thin flux rope. We need to find the observational evidence to confirm this result in the future.

Evolution in the System after the Catastrophe
In this part of work, we look into the detailed evolution in the system for the flux rope with the initial radius r 00 = 0.04d. Correspondingly, the initial value of the electric current density inside the flux rope is J = 1, and six cases for different strengths of Dipole 2 are used to model the process of the flux emergence. The other parameters for the initial configuration of these cases are listed in Table 1. For all the cases, we set j 0 = 1.07 × 10 −2 A m −2 , M = 1, x d /d = 8, and y d /d = 4.7.
As mentioned earlier, the system manifests similar evolutionary behavior for S taking values of 1−5. Figure 5 displays the distribution of the plasma density and the magnetic field for S = 3 and the initial position (x h /d = 0.09, y h /d = 1.06). We notice that when the simulation starts, the flux rope leaves the initial position (marked by a black plus sign) and adjusts itself to a new equilibrium location within 4000 s (Figures 5(a) −5(d)), and the displacement between these two positions is quite small as shown in Figure 5(d). Overall, we find that the global configuration does not show apparent variation. To illustrate the adjustment in the flux rope position more clearly, Figures 5(e)−5(h) display a time sequence of an area around the flux rope marked by the black box in Figure 5(a). We can see easily that the flux rope initially oscillates around the equilibrium location and that the equilibrium location is soon reached. This occurs for all the cases with S between 1 and 5.
So far we have demonstrated that the system could evolve through a set of quasi-static states to the critical point in response to the newly emerging flux. Figure 6 displays the aforementioned evolution of the plasma density and the magnetic field. Figures 6(a)−6(e) show the magnetic configurations as the flux rope reaches the position of the new equilibrium. These panels display a process of the quasi-static evolution and the deformation of the magnetic configuration during the magnetic emergence. Figure 6(f) presents magnetic configuration with S = 6; no equilibrium in the configuration could be found after the flux rope leaves the initial location in this case, and the consequent evolution of the system is dynamic.
Evolutionary behaviors in the system for S = 5 and S = 6 implies that a critical point or state exists between S = 5 and S = 6. After several runs of calculations, we locate the critical point at S = S * = 5.75. When the simulation commences at this point, the flux rope also duplicates an oscillating feature with the amplitude getting larger and larger and eventually leaves the initial location far away, which suggests the occurrence of the catastrophe. Unlike the flux rope in the symmetrical system that moves in the y-direction only, the flux rope studied here moves in both x-and y-directions because of the asymmetry. In this case, the catastrophe means the loss of equilibrium of the flux rope in both the x-and y-directions. Figure 7 describes the equilibrium state of the flux rope from the quasi-static evolution process (S changes from 1 to 5) to the catastrophe (S = 5.75), where the solid curves are from analytic solutions, the circles are simulation results in the quasi-static evolution process, the black asterisks are the initial position of the flux rope as the catastrophe commences, and the red asterisks indicate the location of the flux rope in the new equilibrium. The existence of the new equilibrium suggests that the loss of equilibrium does not evolve to a plausible eruption, though magnetic reconnection takes place following the catastrophe.
We note here that the failure of the eruption is not a geometric artifact. Instead, it is indeed due to the true stabilizing factor provided by the complex magnetic structure studied in this work. Lin & Forbes (2000) and Lin (2002) pointed out that the magnetic configuration including the detached flux rope still suffers from the Aly-Sturrock paradox such that the catastrophe cannot develop to a plausible eruption if reconnection is prohibited or even if the rate of reconnection is not fast enough. This implies that the flux rope cannot escape from the confinement of the background field if the background is not erased fast enough by reconnection. Here the meaning of "the background is not erased fast enough by reconnection" is twofold: First, it means that the reconnection process is indeed slow and the magnetic tension that prevents the flux rope from escaping does not decrease fast enough. Second, the background field itself is very strong, and reconnection cannot erase it and weaken the magnetic tension at a reasonably fast rate. Hence, we argue that the failure of the eruption suggested by Figure 7 corresponds to the second case. Unlike the magnetic configuration studied by Lin & Forbes (2000) and Lin (2002), in which the background magnetic field is due to one source only, the background field in the configuration studied here is of the combination of two sources. So the confinement of the background field to the flux rope is apparently stronger than that of Lin & Forbes (2000) and Lin (2002). This is why successful development of the catastrophe to the eruption is not guaranteed here. From another perspective, the analytic solutions shown in Figure 7 more or less imply such a result since a couple of equilibrium curves exist in the same region in which the loss of equilibrium occurs in the numerical experiments. Figure 8 shows the dynamic evolution in the magnetic configuration after the flux rope loses its equilibrium. The colored shading represents the density distribution, and the black solid lines are the magnetic field lines. Unlike previous cases in which the flux rope continues to move after leaving the equilibrium (see, e.g., Mei et al. 2012a;Ye et al. 2019), the flux rope in the present case eventually stops at a new position between Dipole 1 and Dipole 2 with different (x h , y h ) from the initial one. This is a so-called failed eruption, in which the overall background magnetic field configuration is not destroyed by magnetic reconnection to a certain level, so that the magnetic tension that prevents the flux rope from escaping is still strong enough. During the fast evolution, on the other hand, the flux rope continues to expand and the radius reaches maximum of r 0 = 0.46 when the system realizes the global equilibrium at (x h /d = 3.15, y h /d = 4.06). At the earlier stage of the dynamic process, the flux rope is launched, generating a fast disturbance propagating at a speed of 220 km s −1 in front of the flux rope, which is consistent with the results of Xie et al. (2019). The difference between our results and those of Xie et al. (2019) is that the disturbance in the case of Xie et al. (2019) is symmetric about the y-axis and the disturbance in the present case is asymmetric. In addition, the disturbance to the environment caused by the catastrophe in the present case is also weaker than previous cases (see, e.g., Mei et al. 2012b;Xie et al. 2019) because the background field is produced by two dipoles in the present case, and the magnetic tension that prevents the flux rope from moving quickly is apparently stronger than that in previous cases. Figure 9 displays variations of the flux rope displacements in both the x-and y-directions (Figure 9(a)), as well as the corresponding velocities (Figure 9(b)). We notice that, at the early stage following the loss of equilibrium, the flux rope undergoes fast acceleration and quickly reaches the maximum speeds of 15.4 and 22 km s −1 in the x-and y-directions, respectively, and then the flux rope decelerates and experiences oscillation. Eventually, it reaches a new equilibrium. The average speeds before oscillation are 8.7 and 11.2 km s −1 in two directions, which is low compared to those of CMEs detected by LASCO during solar cycle 24 (≈300-400 km s −1 ) according to Lamy et al. (2019), but quite close to those of the transverse event (9.7 km s −1 ) reported by Gopalswamy et al. (2003). Since the catastrophe taking place in this case does not develop to a plausible eruption, it is necessary to invest more efforts in the future in looking into the validity to explain the transverse event according to the present work.
In addition, we also notice that two current sheets form during the catastrophe process, as shown in Figure 10. It is interesting to note that the first current sheet (CS1) is produced as the flux rope moves upward and stretches the magnetic field lines around the flux rope and that the second current sheet (CS2) is formed as the flux rope moves to the right and squeezes the magnetic field nearby. Apparently, CS1 results from the motion of the flux rope in the y-direction, and CS2 results from that in the x-direction. Magnetic reconnection obviously occurs in both CS1 and CS2, but the loss of equilibrium in the flux rope does not eventually develop to a successful eruption, at least for the case of r 00 = 0.04d.

Discussion
The catastrophe studied here is due to the ideal loss of mechanical equilibrium of the flux rope (see, e.g., discussions in Aulanier 2014) although magnetic reconnection does take place in the stage of the quasi-static evolution of the system in response to the slow (compared to the local Alfvénic speed and sound speed) change in the photosphere. Lin et al. (2001) pointed out that the evolution in the photosphere is so slow that any coronal current sheet would have plenty of time to dissipate by means of slow reconnection. In this case, reconnection occurring at any X-point plays a role in helping rearrange the magnetic structure in the global configuration only, not in triggering the eruption (comparing with the case of Mikić & Linker 1994). Lin & van Ballegooijen (2002) further noted that if the rate of reconnection in the system is faster than the rate of evolution in the photosphere, the catastrophe is considered nonideal; otherwise, the catastrophe is ideal.
Generally, in our numerical experiments, on the other hand, it is difficult to achieve such slow evolution without dissipating the flux rope and the surface currents associated with the line tying. Therefore, experiments that evolve the boundary conditions in time may also create current sheet near the X-points, and these currents might alter the subsequent evolution of the system. To elude this dilemma, we evolve the system by changing the boundary condition step by step: we allow the system to relax for a while each time the boundary changes. We observed that the system could eventually reach the new equilibrium if the variation of the boundary is within a given range, although a couple of X-points do exist in the configuration. Hence, the catastrophe occurring in the system studied here is due to the ideal loss of mechanical equilibrium according to Lin et al. (2001) and Lin & van Ballegooijen (2002).  Feynman & Martin (1995) reported that the magnetic flux that emerges outside the filament channel and possesses the reconnection-favorable orientation (namely, the orientation of the magnetic field supports the formation of the X-line between the preexisting and newly emerging fluxes) would usually result in a successful eruption (with some exceptions). Numerical experiments of Chen & Shibata (2000) seem to support this conclusion. A 3D and self-consistent simulation performed by Roussev et al. (2012) found that the magnetic reconnection between the emerging flux and the preexisting field would not only decrease the magnetic tension but also lead to the formation of the flux rope. Although their results suggested the simple rule that the emerging flux with the reconnection-favorable orientation would trigger the solar eruption, they indeed revealed a more complex physical process during the emerging and implied the complexity of the effect of the newly emerging field. The results that we obtained here clearly indicate that the polarity of the emerging flux is not the only issue that governs the eventual evolutionary behavior of the system; the physics behind the rule that determines whether a successful eruption can be triggered by the emerging flux could be very complex, which is consistent with theoretical (e.g., Lin et al. 2001) and observational (e.g., Zhang et al. 2008) results.
Furthermore, the existence of CS1 and CS2 allows the process of magnetic reconnection to take place, as indicated by Figure 10, which directly suggests the occurrence of the flare in the related coronal magnetic configuration. However, the fact that the flux rope eventually fails to escape implies that in this case no CME can be expected. This may help explain why the solar active region AR 12192 was a flare-rich one, although a CME-poor one. This active region possessed very complex magnetic structure, hosted the largest sunspot group, and produced the most violent solar flare in solar cycle 24 among six X-class flares and 29 M-class flares within 12 days (Chen et al. 2015), but only six jet-driven CMEs with velocities ranging from 200 to 300 km s −1 (Panesar et al. 2016). Sun et al. (2015) listed three issues that may account for the unexpected behavior of AR 12192 as follows: weaker nonpotentiality, strong overlying field, and small flare-related field changes. Looking into the characteristics of the magnetic configuration shown in Figure 8 or Figure 10, we pointed out that although the magnetic configuration in the present work might not be very complex, its features somehow fit these three issues well.
As a follow-up of Lin et al. (2001), the present work aims at confirming via numerical experiments that the catastrophe in reality does occur in the configuration with a complex background field. Our experiments indicate that within a given scope of parameters for the photosphere, stable equilibria do realize in the magnetic configuration of our interest, and the catastrophe indeed occurs as the system evolves to the critical point, although the difference between analytical results and numeric ones exists. As discussed earlier, Figure 5 confirms that the stable equilibrium in the magnetic configuration could be realized in the coronal environment, and Figure 6 indicates that the catastrophe could occur in such a configuration through a set of equilibria in response to the change in the background magnetic field. As for whether or not the catastrophe happening here would smoothly develop to a plausible eruption eventually is beyond the scope of the present work and will be studied in the future.

Conclusions
On the basis of the analytic model developed by Lin et al. (2001), we performed a set of numerical experiments in order to investigate the evolution in the coronal magnetic configuration in response to the newly emerging flux. Unlike previous cases in which the background magnetic field was produced by a single source and the system possessed symmetry, the background field in the present case results from two sources and no symmetry exists in the system. As expected, on the other hand, the flux rope in the present system can still find equilibrium positions, as in the analytical case, and the configuration may evolve gradually through a set of equilibrium states when the background field changes in the fashion used to describe the newly emerging flux as done by Lin et al. (2001). Deviations of the numerical results from the analytical ones exist obviously as a result of the nonforce-free environment in the simulation, but they are apparently larger than those occurring in previous cases with symmetry (e.g., Xie et al. 2017). For the flux ropes of the given sizes, however, it may not find equilibrium locations, as the parameters for the background field take values beyond specific ranges. This suggests that the catastrophe may indeed happen in the system we studied here.
The main results of this work can be summarized as follows: (1) The size of the flux rope affects the evolutionary behavior of the system in an important way; the larger the radius of  the flux rope is, the more easily the system loses the equilibrium. In the analytic solution, the radius of the flux rope is assumed to be small compared to the global scale of the system so that the external magnetic field over the cross section of the flux rope could be treated uniformly, and vanishing of the force due to the external field at the center of the flux rope gives rise to the equilibrium of the whole flux rope. In reality, on the other hand, the radius of the flux rope is not necessarily small compared to its length, and the uniformity of the external field over the cross section of the flux rope is not easy to realize. Hence, vanishing of the force at a given location over the cross section of the flux rope cannot guarantee the whole flux rope force-free. This indicates that it is hard for the configuration including a thick flux rope to find an equilibrium state.
(2) In reality, or for the magnetic configuration including a flux rope of finite radius, the internal and external equilibria of the flux rope couple to one another and cannot be dealt with separately. However, the basic equilibrium behavior of the system is mainly governed by the external equilibrium of the flux rope, and the disturbance to the internal equilibrium could be absorbed by the flux rope itself via adjusting its radius.
(3) The asymmetry in the magnetic configuration studied in this work provides the flux rope one more dimension to move. The motion of the flux rope could be in both the xand y-directions no matter whether the evolution of the system is quasi-static or dynamic. This allows us to gather much richer information about the evolutionary pattern of the system and implies more possible styles of the catastrophe, as suggested by the analytic solution (see, e.g., Lin et al. 2001), which suggests that the catastrophe is easier to occur than that in the symmetric configuration. (4) In the case of r 00 = 0.04d, the equilibrium in the system is lost and the catastrophe takes place at S = S * = 5.75 after evolving through a set of equilibrium states in the quasi-static fashion as the strength of dipole increases. The motion of the flux rope displays an apparent transverse component, as shown by many events (see, e.g., McCauley et al. 2015). However, the catastrophe does not develop to a plausible eruption in the present work and may account for the failed eruptive event (see, e.g., Ji et al. 2003 for the observational event). (5) Two current sheets form during the catastrophe; one (CS1) results from stretching of the closed magnetic field due to the loss of equilibrium of the flux rope in the ydirection, and another one (CS2) from squeezing of the magnetic field due to the loss of equilibrium in the xdirection. Magnetic reconnection obviously takes place inside both current sheets, which tends to help the loss of equilibrium develop to a successful eruption. But, at least, the eruption fails in the case of r 00 = 0.04d. This implies that no simple and universal rule exists in this case that may relate the likelihood for an eruption to the orientation of the emerging flux (see discussions in theory by Lin et al. 2001 and those in observations by Zhang et al. 2008 for more details). Figure 10. The evolution of the current and magnetic field for the catastrophe case. "CS1" and "CS2" represent two different current sheets that have opposite polarity.