Excitation of Quasiperiodic Fast-propagating Waves in the Early Stage of the Solar Eruption

We propose a mechanism for the excitation of large-scale quasiperiodic fast-propagating magnetoacoustic (QFP) waves observed on both sides of the coronal mass ejection. Through a series of numerical experiments, we successfully simulated the quasi-static evolution of the equilibrium locations of the magnetic flux rope in response to the change of the background magnetic field, as well as the consequent loss of the equilibrium that eventually gives rise to the eruption. During the eruption, we identified QFP waves propagating radially outward of the flux rope, and tracing their origin reveals that they result from the disturbance within the flux rope. Acting as an imperfect waveguide, the flux rope allows the internal disturbance to escape to the outside successively via its surface, invoking the observed QFP waves. Furthermore, we synthesized the images of QFP waves on the basis of the data given by our simulations and found consistency with observations. This indicates that the leakage of the disturbance outside the flux rope could be a reasonable mechanism for QFP waves.


INTRODUCTION
Solar atmosphere is a highly structured medium capable of supporting various types of waves, which might belong to the families of the Magnetohydrodynamic (MHD) waves in the solar corona including fast-mode (Williams et al. 2001(Williams et al. , 2002;;Nakariakov et al. 2003) and slow-mode waves (De Moortel et al. 2000;Wang et al. 2009;Ofman & Wang 2002) , and they have been extensively investigated by many authors (e.g., see also Li et al. 2020;Nakariakov & Kolotkov 2020;Banerjee et al. 2021;Nakariakov et al. 2021;Wang et al. 2021b).They impose additional constraints on the explanation or understanding of fundamental physical processes occurring in solar activity, such as particle acceleration and magnetic energy release (McLaughlin et al. 2018).In addition, these perturbation may interact with the plasma and the magnetic filed nearby, delivering the energy to the background atmosphere, which may eventually compensate the radiative losses of thermal energy in the atmosphere (Porter et al. 1994;Van Doorsselaere et al. 2020).
Global large scale EUV waves, usually observed in the EUV wavelengths, are considered a ubiquitous phenomenon (Thompson et al. 1998;Wang et al. 2009;Lin et al. 2010) associated with the disturbance to the atmosphere by the eruption.They typically manifest as bright, ring-like structures that gradually propagate outward, sweeping across a significant portion of the solar surface and leaving behind a dimming area (Thompson & Myers 2009;Mei et al. 2012b).In their propagation, refraction and reflection may occur to them as they go through regions of different densities or through the boundaries separating magnetic fields of different topologies, like that between active regions and coronal holes (Gopalswamy et al. 2009;Shen et al. 2019;Wang et al. 2009Wang et al. , 2015;;Xie et al. 2019).Furthermore, if the wave is energetic enough, the nearby magnetic fields of opposite polarity could even be pushed to approach each other, leading to the occurrence of magnetic reconnection (Zhou et al. 2020).cylinder).Their results revealed that the straight cylinder could behave as waveguides, allowing long-period (57 s period) and short-period (5.8 s period) oscillations to occur within the waveguide.They used a fast gyrosynchrotron code to produce synthetic signals in the microwave band, explaining the generation of quasi-periodic pulsations (QPPs) in solar flares.It is worth noting that although not explicitly mentioned in their study, their results suggest that when short-period oscillation signals pass through the boundary of the cylinder, a portion of the oscillation leaks outside the cylinder and forms multiple wavefronts, which is similar to the multiple QFP waves observed in the solar eruption.Analyzing the observational data from the 1.6-meter Goode Solar Telescope, Yuan et al. (2023) discovered persistent transverse waves in the chromospheric umbral fibrils of a sunspot with strong magnetic field.They reproduced these waves using a two-fluid MHD simulation.Interestingly, their results indicate that perturbations confined within the chromospheric fibrils (acting as waveguide structures) oscillated back and forth, and leaked out at the boundaries, giving rise to multiple wavefront signals.
In order to investigate the QFP phenomena and mechanisms for their generation during solar eruptions, we will first examine the evolution in the coronal magnetic structure including a magnetic flux ropes through a set of quasistatic states to eventually losing the equilibrium.Then, we follow the consequent dynamic evolution in the system, investigate the disturbance to the background environment and the occurrence of QFP waves in this process, look into the physical property and mechanism triggering the QFP wave.We will provide a detailed introduction of the model and the relevant numerical techniques used in the numerical simulations in Section 2, and display our results in Section 3. Synthetic images deduced according to the numerical experiments and the associated comparisons with observations will be given in Section 4. Finally, we summarize this work in Section 5.

DESCRIPTIONS OF THE NUMERICAL METHOD
We utilize the MPI-AMRVAC code (Keppens et al. 2012;Porth et al. 2014;Xia et al. 2018) to study the evolution in the coronal magnetic configuration through a set of quasi-static equilibrium states to the loss of equilibrium.Our calculations are involved in solving the following MHD equations, which include the gravity and the thermal conduction: where ρ, v, B, p, e, and T are density, velocity, magnetic field, pressure, energy, and temperature, respectively.In the description of thermal conduction in the plasma, we employ the classical Spitzer model (Spitzer 1962) with κ ∥ = 9000 W m −1 K −1 for the coefficient of the thermal conductivity parallel to the magnetic field.The gravity is g(y) = g 0 /(1 + y/r 0 ) 2 ŷ with g 0 (g 0 = 274 m s −2 ) being the value of the gravity on the solar surface, r 0 (r 0 = 6.9551 × 10 8 m) the radius of the Sun, and ŷ the unit vector in the y-direction.Due to numerical diffusion in calculations, we set η to 0 for simplicity.We note here that the thermal conduction parallel to the magnetic field is included in our calculations.The point of doing so is that the thermal conduction in the coronal environment is a dominant process of the energy transportation, and that the field-aligned conduction is powerful compared to that perpendicular to the field in the corona.Therefore including a field-aligned thermal conduction in calculations here is of consideration for the model to be as realistic as possible.
In addition, the energy e and the pressure p are given as: where γ = 5/3.To suppress spurious oscillations in our computations, we employ a third-order asymmetric TVD limiter and utilize a "three-step" scheme of the third-order predictor-corrector.Regarding the divergence-free of magnetic Hu et al.
field, we apply the "typedivbfix=linde" condition in practical calculations to maintaining the numerical conservation of the magnetic field divergence (∇ • B = 0).According to Equation ( 5), in low-beta plasma atmospheres, numerical errors in the magnetic energy may potentially exceed the internal energy, leading to negative values for the internal energy or gas pressure when solving the complete energy equation.To address this issue, we simultaneously solve both the complete energy equation and the internal energy equation.

the gravitationally stratified atmosphere
Calculations are conducted in the Cartesian coordinate system (x, y), with the origin located at the bottom of the photosphere.The x-axis lies on bottom of the photosphere, and the y-axis is the radial direction away from the Sun.The simulation domain spans over −1.6L 0 ≤ x ≤ 1.6L 0 and 0 ≤ y ≤ 3.2L 0 , where L 0 = 2 × 10 10 cm.In our simulation, we employ a two-layer atmospheric structure to model the solar atmosphere, where the region with higher density (y < h p ) represents the photosphere, and the region with lower density (y > h p ) represents the corona.The temperature of the gravitationally stratified atmosphere is set to: For the photospheric layer, we follow the setting of Mei et al. (2012b), and take h p = 2.5 × 10 7 cm.The isothermal property of the atmosphere in the low corona yields the gas pressure: where p c and p p are the pressure at the bottom of the corona and the bottom of the photosphere, respectively, with Here, λ p and λ c are the pressure scale heights for the photosphere and the corona, respectively: where H ea , k b , and m p represent the helium abundance, Boltzmann constant, and proton mass, respectively.The density distribution could be determined according to the isothermal property of the atmosphere.

The structure of magnetic field and filament
Following Isenberg et al. (1993), we model the filament with a current-carrying magnetic flux rope levitating in the corona.To ensure a smooth transition of the parameters from the flux rope to the background corona, we impose a thin shell outside the flux rope with a radius of R.This shell is known as the PCTR (prominence-corona transition region, e.g., see also Mei et al. 2012a), with a thickness of ∆ = R/5.Based on the distance r − from any point in space to the center of the flux rope (r − = x 2 + (y − h) 2 , where h is the initial height of the flux rope), the entire simulation domain is divided into three regions: the interior of the flux rope (Ar1), for 0 ≤ r − < (R − ∆/2); the PCTR (Ar2), for (R − ∆/2) ≤ r − ≤ (R + ∆/2); and the exterior of the flux rope (Ar3), for r − > (R − ∆/2).According to Mei et al. (2012a), the current density of the flux rope is: where j 0 is a constant.
The initial magnetic field results from three sources: the internal current of the flux rope centered at y = h, the mirror current centered at y = −h, and the buried magnetic quadrupoles beneath the photosphere at y = −d.So the initial magnetic field is described by: where with r + and r d being the distances from any point in space (x, y) to the mirror current and the quadrupole, respectively, and σ denoting the relative strength of the background magnetic field (interested readers refer to Forbes 1990).
In the magnetohydrostatic equilibrium, the initial gas pressure within the flux rope is given by: The second term on the right-hand side of equation ( 18) is the pressure generated by the azimuthal magnetic field, and the third term on the right-hand side is the pressure generated by the axial magnetic field.The external equilibrium condition of the flux rope determines the intensity of the electric current within the flux rope, while the pressure due to the azimuthal magnetic field inside the flux rope governs its internal equilibrium.Prior to numerical experiments, we adjust the axial magnetic pressure according to specific requirements to regulate the internal pressure of the flux rope so that its temperature and density could be as close as possible to those in the true corona.
To smoothly connect the temperature inside the flux rope to that in the background, we set the temperature distribution as follows: where T f = 5 × 10 4 K and T c = 10 6 K, the density distribution inside the flux rope can be obtained through the ideal gas law.
Since we need to study the gradual evolution in the background magnetic field driving the coronal magnetic field to evolve through a set of stable equilibrium to a critical state in a quasi-static fashion, to ensure consistency across all sets of numerical experiments, the simulation grid for each experiment is uniformly set to 3840 × 3840.Figure 1(a) illustrates the initial magnetic configuration and the density distribution in the gravitationally stratified atmosphere, Figure 1(b) displays the distribution of fast speed.We notice that the speed around the flux rope is almost uniform, which allows the wave around the flux rope to propagate almost isotropically.This issue will be further discussed later.3. RESULTS et al. (1993) studied analytically the change in the background magnetic field driving the coronal magnetic field structure to evolve through a set of stable equilibrium states to a critical state in a quasi-static approach.They described the equilibrium height h of the flux rope as a function of the background magnetic field σ (see the black solid line in Figure 2(a)).Before investigating the quasi-periodic fast-propagating (QFP) wave in the early stage of the solar eruption, it is essential to verify whether the process investigated by Isenberg et al. (1993) can indeed occur.In other words, we need to confirm that the coronal magnetic configuration would indeed lose the equilibrium after evolving quasi-statically through a set of equilibrium states to the critical state in response to the change in the background field.This is important for two reasons: first, an unstable state should evolve from the stable state, so the magnetic configuration in reality needs to be in equilibrium for a while prior to the eruption.Second, numerical experiments starts with a configuration in non-equilibrium might bring unphysical issues into calculations, and yielding the unrealistic results and conclusions.

Isenberg
Figure 2(a) shows the continuous curve taken from Isenberg et al. (1993), depicting the equilibrium height h of the magnetic flux rope as a function of the strength of the background field σ.Due to discrepancies between analytical and numerical solutions, when we tentatively take the parameters from the continuous curve (black dots) for the initial magnetic configuration in our numerical experiments, it is not surprising that the equilibrium locations (as indicated by the row of blue inverted triangles located at the bottom) of the flux rope in the numerical experiment are not located on the continuous curve, instead the equilibrium locations determined in the numerical experiment always deviate from those deduced from the analytical solution.Moreover, the closer the tentative positions are to the critical point on the analytic curve (near the black cross), the larger the deviation is.The four black dots on the curve correspond to σ = 0.88, 0.86, 0.84, and 0.82 respectively.The results indicate that when the magnetic flux rope is initially located at heights corresponding to these four points on the equilibrium curve, it undergoes oscillations for a while and eventually finds a stable equilibrium position at the height marked by the triangles closest to the black equilibrium curve after self-adjustment.This demonstrates that the stable equilibrium states and the associated quasi-static evolution in the configuration described by Isenberg et al. (1993) can indeed be realized.Driven by the change in the photospheric magnetic field and plasma, the system evolves quasi-statically from one stable equilibrium state to another.
However, a change occurs near the critical point on the equilibrium curve (marked by the cross).The critical point is the joint where the stable section on the equilibrium curve connects to the unstable one.In reality and numerical experiment, no magnetic configuration could exist if its parameters are located on the unstable section of the equilibrium curve.So the loss of equilibrium in the magnetic configuration inevitably occurs as the system evolves to the critical point along the stable section of the equilibrium.In this study, we found that when the initial position is located at point σ = 0.80 (the location marked by symbol "+" in Figure 2(a)), the flux rope could no longer find equilibrium location through self-adjustment during the consequent evolution, and the motion of the flux rope becomes highly dynamic.Evolution behaviors of the system at σ = 0.82 and σ = 0.80 suggest the existence of a critical value between these two σ values, to which the corresponding magnetic configuration reaches to the critical state where stable and unstable equilibria join.To locate the critical point, we conducted four sets of test experiments between σ = 0.82 and σ = 0.80, using σ= 0.815, 0.8125, 0.810, and 0.8075.The change in the height of the flux rope versus time for each group of experiments is shown in Figure 2(b).From the figure, we can see that when σ = 0.820, 0.815, and 0.8125, the flux rope quickly finds a new equilibrium several rounds of oscillation.When σ = 0.8075 and 0.805, the flux rope fails to find an equilibrium and rapidly moves upward.When σ = 0.810, the flux rope oscillates and remains at a certain height for a while but eventually loses its equilibrium.This indicates that the catastrophe occurs at this point and it is the critical point found in our numerical experiments, denoted by the red pentagram on the analytical curve in Figure 2(a).The height at which the flux rope could stay for a while in the numerical experiment is marked by the red triangle in the inset of Figure 2(a).Thus, connecting those triangles closest to the equilibrium curve in Figure 2(a) yields the equilibrium curve for our numerical experiments.Figure 2(a) also shows that the equilibrium curve obtained from numerical experiments in the present work differs slightly from the analytical results given by Isenberg et al. (1993).This difference is primarily due to the existence of the gas pressure, gravity, and numerical dissipation in simulations (e.g., see also Chen et al. 2022).However, as expected, these issues do not introduce significant deviations between the numerical and the analytical results because the lower corona is very close to a force-free environment.Moreover, in order to find the regions near the equilibrium curve where the flux rope can stably stay, we performed a set of test calculations, and roughly outlined a range of values of the relevant parameters that govern the equilibrium property of the system.We find that the system could eventually reach the equilibrium if the initial value of the relevant parameters were taken in the region below the red triangles and above the solid curve shown in Figure 2(a).If the initial values were taken somewhere below the solid curve, the flux rope will move upward and find the equilibrium position in that region as well.The system will never find the equilibrium as the initial value is taken at or above the red triangle, or in the region at the left to the red asterisk on the curve.
So far, we have confirmed that the magnetic configuration shown in Figure 1 (a) can evolve, in response to the change in the background field, quasi-statically through a set of stable equilibrium configurations to the critical configuration, and then loses its equilibrium, triggering the eruption.The main purpose of this work is to investigate the wave phenomenon in the eruption, so to save computational resources, we directly start our numerical experiments at σ = 0.81 and focus on the issues related to the waves occurring in the eruption.

The global evolution
Figure 3 illustrates evolutionary features in the system at different times.In the early stages of evolution at t = 16 s, we see an arc-shaped wavefront (labeled as "w 1 ") that includes the flux rope propagating outward while expanding and extending in various directions simultaneously.As it extends to the bottom boundary, an echo occurs as indicated by black arrows in Figure 3(g), a phenomenon commonly reported in both numerical experiments and observation (e.g., see also Xie et al. 2019;Wang et al. 2009Wang et al. , 2015;;Mei et al. 2012b;Liu et al. 2012;Asai et al. 2012).In addition to the density distribution, we also examine the distributions of gas pressure p, horizontal velocity component v x , and vertical velocity component v y (see rows 2 to 4 in Figure 3).We notice that, comparing with the distribution of the density, those of the velocity and the pressure reveal more fine structures.These features may not be easily recognized in images that basically display the brightness of targets.If the information about the Doppler shift could be obtained, say via spectroscopic approach, they may become apparent.
At t = 33 s, as shown in the second column of Figure 3, two semicircular wave fronts ("w1" and "w2") are propagating forward of the flux rope.In contour plots for velocity (Figure 3(e) and 3(h)), the wave patterns are quite rich and impressive.Previous numerical experiments (e.g., see also Xie et al. 2019;Wang et al. 2009Wang et al. , 2015;;Mei et al. 2012b) usually displayed only one wave front, while the eruption process investigated in the present work produces three sets of wave fronts, which can even be clearly distinguished in the density distribution plot (Figure 3(c)).
In general, the generated waves seem to propagate isotropically, i.e., the local propagation speed does not seem to depend on the direction.This is governed by the fact that the fast speed is almost uniform in the region near the flux as shown in Figure 1(b).These wave fronts appear strikingly similar to the quasi-periodic fast magnetoacoustic wave trains displayed by Liu et al. (2012).We shall discuss these results later in more detail and compare them with observation.Furthermore, density contours in Figure 3(a) through Figure 3(c) showed that the wavefront at the low altitude is more significant than that at the high altitude .This difference may be attributed to the gravitationally stratified atmosphere, in which the plasma is dense at low altitudes and tenuous at high altitudes, therefore, the response of low atmosphere to the wave is more apparent than that of high atmosphere.

Characteristics of Wavefronts
According to Figures 3(a)-3(c), these fronts appear sequentially in roughly symmetric shapes ahead of the flux rope.To look into their physical properties, we further analyze their apparent behaviors.Figure 4(a) displays the density distribution at t = 40 s, with black solid lines for magnetic field lines.From the figure, we see that these two wavefronts propagate successively on the solar surface to a distance of a few 10 4 km, and an obvious deflection of the magnetic field occurs on the wavefront with an apparent enhancement in the tangential component of the magnetic field.This suggests that these waves are the fast-mode shock.
To quantitatively reveal the property of these fronts, we analyze the changes in the relevant parameters before and after sweeping of the front along two white lines, s1 and s2.Here line s1 crosses the front horizontally, with the upstream (undisturbed medium) on the right side of the front and the downstream (disturbed medium) on the left side, s2 crosses the shock front vertically, with the upstream on the top side and the downstream on the bottom side of the front.Figures 4(b)-4d show the distribution of density, pressure, and temperature along s1, while Figures 4(e)-4(g) show the distributions along s2.
We found that the density and the pressure along s1 exhibit clear discontinuities, while the temperature change is relatively smooth.After the passage of the first wavefront, the plasma density, pressure, and temperature increase by 1.25, 1.23, and 1.1 times, respectively; after the passage of the second wave, the increases become 1.17, 1.18, and 1.08 times, respectively.This indicates that the first wavefront is stronger and impacts the environment more significantly than the second one.The distributions of the above parameters along s2 are similar to those along s1, and also exhibit clear discontinuity layers in density and pressure.After the passage of the first wavefront, the plasma density, pressure, and temperature increase by 1.48, 2.25, and 1.7 times, respectively; and on the second wavefront, the corresponding changes become 1.48, 2.14, and 1.41 times, respectively.In addition, we also noticed that the properties of the plasma and magnetic field along s1 and s2 show slight differences because of the effect of gravitational stratification.Furthermore, apparent similarities could be recognized when comparing the magnetoacoustic shock structure shown in Figure 4 with the formation of the shock as a result of the leakage of fast-mode waves from a plasma slab studied by Pascoe et al. (2017).In summary, the observed discontinuity layers along s1 and s2 suggest again that these disturbances are fast-mode shock.
To further confirm the fast-mode properties of these perturbations, we compare their velocities with the local fast magnetosonic wave speeds.Figures 5(a) and 5(b) show the variations of the velocity versus time for the first and second wavefronts, respectively.The blue symbols "x" are for the velocities of the wavefronts, while the red dots for the local fast speed that could be deduced from Figure 1(b) accordingly.In general, as the propagation distance increases, the velocities of the wavefronts decrease.From Figure 5(a), we can see that in the early stages, the velocity of the first wavefront is greater than the local fast magnetosonic wave speed, with an average Mach number of 1.12 (relative to the fast magnetosonic wave speed); on the other hand, as the wave further propagates, it gradually degrades into the fast magnetosonic wave.As for the second wavefront, its initial velocity (approximately 1700 km s −1 ) is lower than that of the first wavefront (approximately 1900 km s −1 ), and its velocity is approximately equal to the fast magnetosonic wave speed.Moreover, the wave patterns shown in Figure 3 are suggestive of the periodical generation of these waves.To investigate this process, we conducted a wavelet analysis of the energy variation versus time at a specific location outside the flux rope (labeled as "A" in Figure 1).The results are presented in Figure 6(a), revealing a disturbance signal with a period of approximately 20 seconds.Simultaneously, we also performed a similar analysis of the energy disturbance at the center of the flux rope, and the results are displayed in Figure 6  internal disturbance also exhibits clear periodical feature, and the period is roughly the same as that of the external disturbance.This suggests the connection of the origins of the two disturbances to each other, which we will further discuss below.

Origin of Wave Trains
We calculate the velocity divergence at different times (Figure 7 and the supplementary movie) to analyze the origin and property of disturbances.According to Equation (1), the velocity divergence describes the compression (∇ • v < 0) and expansion (∇ • v > 0) of the plasma, making it an effective indicator of the shock.
As mentioned earlier, in the numerical experiments, the initial position of the flux rope is close to the critical point, so it stays in this position for a short period before losing the equilibrium and then rapidly moving upward after losing the equilibrium.Simultaneously, the internal equilibrium of the flux rope is also disrupted.Because of the compressibility of the plasma, once a disturbance occurs, it propagates outward in the form of the wave.With losing the internal equilibrium, two waves are invoked inside the flux rope, and propagate inward and outward, respectively.Some of the outward-propagating waves eventually leave the flux rope and enter the surrounding corona, while other part is reflected at the edge of the flux rope, forming inward-propagating waves.The inward waves eventually converge at the center of the flux rope and then rebound outward.As a result, we can see an image of the flux rope constantly oscillating and waves continuously propagating outward from its edges.The oscillations of the flux rope eventually dies down as a result of the dissipation and the wave leakage into the external medium.
Figure 7 qualitatively describes such a process.At t = 1 s (see Figure 7(a)), an outward-propagating wave (labeled as "w1") is observed outside the flux rope, corresponding to the first wave shown in Figure 3.It is caused by the loss of equilibrium of the flux rope and its rapid upward motion.Meanwhile, a circular inward-propagating disturbance (labeled as "p") appears inside the flux rope.Between w1 and p, a circular boundary interface exists and does not change significantly with time, which is actually the boundary of the flux rope (labeled as "b").From Figure 7(b) and 7(c), we realize that w 1 continues to propagate outward, while p contracts further toward the center of the flux rope until it shrinks to a point.
Later, at t = 10 s (see Figure 7(d)), the downward portion of w1 has propagated into the lower solar atmosphere, while p starts expanding outward from the center of the flux rope.At t = 15 s (see Figure 7(e)), the downward portion of w1 produces aforementioned echoes.Simultaneously, disturbance p encounters boundary b, and a portion of it transmits through b and propagates outward, forming the second wavefront (labeled as "w2" in Figure 7(e)), which is consistent with the results of Pascoe et al. (2017).Another portion reflects and propagates back toward the center of the flux rope again.The propagation of p as described above forms a cycle that continuously triggers oscillations inside the flux rope, leaking waves on its surface until the extra energy inside the flux rope totally dissipates.This scenario well displays the creation and propagation of various waves/shocks around the flux rope during the eruption, which could find observational counterparts easily.

OBSERVATIONAL CONSEQUENCES
The quasi-periodic fast-mode magnetoacoustic waves propagating on the solar surface triggered by CME have been reported by a few authors.Kumar et al. (2017) reported several wavefronts around the flux rope observed by the white-light coronagraph on board STERO-A (Figure 8(a)).Liu et al. (2012) 2019) for more information regarding the identification and analysis of the echo on the solar surface.
To compare with observations, we create synthetic images in white-light and EUV wavelength on the base of the simulation results.A synthetic images is basically constructed according to instrument response functions and atomic databases.The instrument response function provides the response of the detector in a given wavelength to the plasma radiation at a given temperature.The detector responds to the radiation in the way of:  and 8(f)), we find that the simulation results exhibit multiple wavefronts in both wavelengths, but observations display little bit different scenario such that these waves are almost invisible in front of CME, but apparent on the side of CME.This could be due to the fact that the density on the CME side is higher than that in front of CME, and the disturbance on the CME side may produce more apparent observational consequences; in the numerical experiments, on the other hand, manifestations of the synthetic image could be manipulated by adjusting the range of color bar.
To further justify that the multiple wavefront scenario is the oscillation inside the flux rope in origin, we deduce the time-sequences of the brightness distributions in directions parallel and perpendicular to the solar surface as indicated by Figure 8(f).The results are displayed in Figures 8(g) and 8(h).Comparing Figure 8(c) with Figure 8(g), we see several wave fronts in both observational and synthetic images, and these wavefront structures originate from the core region of the eruption.Furthermore, our simulation, on which the synthetic images are created, definitely confirms that these wavefronts result from leaking of the oscillation inside the flux rope.
In addition, a set of X-like patterns inside the flux rope are observed.In Figure 8(i), we plot velocity divergence in the region inside box1 in Figure 8(g).The oscillating features in the radial direction of the flux rope could be recognized, and the leakage of the oscillation through the boundary of the flux rope is also seen clearly.In this process, the waves associated the oscillation propagate both inward and outward inside the flux rope.Concentrating and bouncing of the the inward wave at the center of the flux rope constitute the X-shape feature as shown by shown in Figure 8(g) and Pascoe et al. (2013).
Comparing Figure 8(d) with Figure 8(h), we find multiple wavefronts propagating downward in both observations and synthetic images.To study the interaction of the downward wavefronts with the solar surface, we calculate the corresponding velocity divergence in the region indicated by box2 in Figure 8(h), and the result is plotted in Figure 8(j), which shows clearly the pattern of the echo.The similar result was also reported by Xie et al. (2019), with the only difference being that Xie et al. (2019) showed a single echo, whereas several echos could be identified in Figure 8(j).

SUMMARY AND DISCUSSIONS
On the basis of the analytical model of Isenberg et al. (1993), we investigated the evolution in a magnetic configuration that includes an electric current-carrying flux rope and possesses a relatively simple background magnetic field.Our results confirmed that the coronal magnetic structure described by the model can quasi-statically evolve through a set of configurations in the stable equilibrium to an unstable configuration, eventually leading to an eruption in a catastrophic fashion.Comparing with the works by Lin et al. (2001), Kolotkov et al. (2016Kolotkov et al. ( , 2018)), as well as Chen et al. (2022), we notice that the evolutionary behaviors of the magnetic structure studied here is simple and straightforward such that the displacement of the flux rope only moves in the y-direction during both quasi-static and dynamic stages.Lin et al. (2001), Kolotkov et al. (2016Kolotkov et al. ( , 2018)), and Chen et al. (2022) found that the motion of the flux rope could be in both x-and y-directions(here the y-direction is the vertical direction) as the background field is complex, and even the loss of equilibrium in the flux rope may not be able to develop to a plausible eruption eventually.This might be more realistic.
However, the purpose of this work is to look into the origin of the QFP wave during the eruption no matter how the eruption is triggered.So a simpler magnetic configuration was chosen and the straightforward evolutionary behavior in the configuration was studied.After confirming the occurrence of the loss of equilibrium of the flux rope, we subsequently investigated the occurrence of QFP waves in the solar eruption.Our results showed that arc-like wavefronts successively appear on the surface of the flux rope and propagate outward at speed of a few 10 3 km s −1 as the flux rope moves out.
We analyzed the mechanism and propagating features of these wavefronts , and found out they were moving at speeds faster or marginally faster than the local fast magneto-acoustic wave speed, which implies their nature of the fast-mode shock.Looking into the velocity divergence of the wavefront, we discovered that these QFP waves originate from disturbances within the flux rope as a result of the loss of the internal equilibrium, and part of the disturbance leaks through the flux rope boundary out to the surrounding corona, forming the QFP waves observed.When QFP waves reach the lower solar atmosphere, they generate multiple echos.In summary, the role of flux rope in the origin of the wave trains observed in the solar eruption is twofold.First, the loss of the internal equilibrium of the flux rope ignites oscillations inside the flux rope; and second, the flux rope behaves like a leaking wave guide from which the oscillation successively leaks out to the surrounding corona, invoking the wave train with period of several tens of second.The numerical simulation of Wang et al. (2021a) showed that the period of the wave train is around 30 s.
In addition to the above point regarding the origin of wave trains, Wang et al. (2021a) suggested that the wave train is a component of the large-scale EUV waves.But this mechanism is hard to account for the fact that various wave fronts come from roughly the same region, and that velocities of different wave fronts are approximately the same.
Numerical experiments performed here duplicated the large-scale QFP waves observed in the lower solar corona by Liu et al. (2012) as well.In addition, Kumar et al. (2017) also identified QFP waves in the high corona on the basis of analyzing the data from the white-light coronagraph COR1 on STEREO-A.To our knowledge, no author reports simultaneous identification of QFP waves in both the low and the high solar corona regions in a single eruptive event.Furthermore, because of the limit to the computational resources, our numerical experiments could not investigate QFP in both regions simultaneously, either; so it is very hard to reveal the relations of the wave train in the low to that in the high corona for the time being.But the possibility exists that the QFP wave in the low corona propagates upward, reaches to the large altitudes, and evolves into that in the high corona.
Overall, to definitively determine the mechanisms for their origins, it is necessary to identify and compare the processes of QFP wave generation and propagation in both the low and higher solar corona regions.This requires observational instruments with large FOV that cover both the low and the high corona regions with sufficiently high time cadence.The existing observational data may not fully meet these requirements, we hope that advanced technology and observational instruments in the future would be able to help solve this puzzle.

Figure 1 .
Figure 1.Initial distributions of (a) the plasma density and (b) fast speed.The black solid curves are the magnetic field lines and point A in (a) is selected for analyzing the period of the QFP as shown later in Fig 6.

Figure 2 .
Figure2.Variations of the equilibrium height of the flux rope versus the strength of the background field.The solid line in panel (a) is the equilibrium curve that was duplicated from the analytical solution ofIsenberg et al. (1993), and the black cross marks the critical point.All the solid dots on the equilibrium curve are used as the initial location of the flux rope and the corresponding strength of the background field for numerical experiments.Blue triangles mark positions where the flux rope finds the equilibrium in numerical experiments, while red triangles outline the boundary above which no equilibrium location exists for the flux rope in numerical experiments.The inset in (a) displays more details of equilibrium locations obtained in both theory (continuous curve) and numerical experiments (blue triangles).The flux rope at locations marked by crosses on the curve are still in equilibrium in theory, but not in numerical experiments.Instead the red asterisk and the red triangle of the same value of σ mark the critical point for the numerical experiment.Panel (b) shows the temporal evolution of the flux rope location in the corresponding experiments with the initial locations marked by solid dots, asterisk, and crosses on the theoretical curve displayed in the inset of (a).

Figure 3 .
Figure 3. Evolutionary features of quasi-periodic fast propagating (QFP) waves at different times.The rows represent density, pressure, horizontal velocity, and vertical velocity, respectively.The arrows indicate the features described in section 3.1, where "w1" and "w2" refer to the first and second wavefronts, respectively.

Figure 4 .
Figure 4. Physical quantities variations in slices crossing wavefronts.(a) displays a snapshot of density and magnetic filed at t = 40 s.Magnetic field lines are depicted by solid gray lines.(b-d) and (e-g) illustrate changes in density, pressure and temperature along slices s1 and s2, respectively.

Figure 5 .
Figure 5. Propagation speeds for the first wavefront "w1" in (a) and the second wavefront "w2" in (b) versus time.The blue symbols "x" are for the velocities of the wavefronts, while the red dots for the local fast speed that could be deduced from Figure 1(b) correspondingly.

Figure 6 .
Figure 6.Periodicity analysis of quasi-periodic wave trains.(a) denotes the periodicity analysis of the fixed point outside the flux rope (marked as "A" in Fig 1), and (b) denotes that of the central point of the flux rope.

Figure 7 .
Figure 7.The velocity divergence at different times.The red arrow "p" is the inner perturbation, and the black arrow "b" the boundary of the flux rope."w1" and "w2" are the first and the second wavefronts, respectively.The animated velocity divergence images display the propagation of internal disturbances and the formation of QFP waves in the time interval from 0 s to 50.67 s with a duration of 2 s.Static image only presents information at a few key moments, while the animation shows the evolution.An animated version of this figure is available.
performed detailed studies of the similar phenomenon observed in EUV by SDO/AIA (Figure 8(b)).Figures 8(c) and 8(d) show the time-sequences of the brightness distributions along curve s3 and line s4 specified in Figure 8(b), revealing information about the propagation of the wavefronts in various directions.Figure 8(c) demonstrated multiple wavefronts (indicated by white arrows) are successively ejected from the boundary of CME core.Box1 in the figure highlights how the upper part of the wavefront leaks out. Figure 8(d) describes the situation in the vertical direction along s4, where multiple downward wavefronts (indicated by white arrows) reach the solar surface, and the echo could be recognized.Box2 in the figure outlines the region where the waves interact with the lower solar atmosphere.Interested readers refer to Xie et al. ( T e , n e )dl.(DN/s) (20) Here, dl is the length element along the line of sight (LOS), DN the number of electrons produced by the photons received by the instrument, n e the electron density, and f i (T e , n e ) the response function of the instrument for density n e and temperature T e .A synthetic image then could be obtained provided the form of f i (T e , n e ), and distributions of T e as well as n e in LOS are known.Figure 8(e) shows the synthetic image in white-light, and Figures 8(f) through 8(h) provide the information in 193 Å. Comparing observations (Figures 8(a) and 8(b)) with the synthetic images (Figures 8(e )

Figure 8 .
Figure 8.Comparison between observations and the simulation.The top row displays the observations while the middle row the simulation.Specifically, (a) depicts an event that occurred on May 7, 2012, where multiple wavefronts were detected on both sides of a flux rope through the white-light coronagraph COR1 (Kumar et al. 2017).(b) shows an event that occurred on September 8-9, 2010, where QFP waves were observed in the AIA image at 193 Å (Liu et al. 2012).(c) and (d) are the time-sequences of the brightness distribution along curve s3 and line s4 in (b).(i) and (j) are magnified images of the velocity divergence within the boxed regions of (g) and (h).
This work was supported by the National Key R&D Program of China No.2022YFF0503804, the NSFC grants 11933009, 12273107, U2031141, and 12073073, grants associated with the Yunling Scholar Project of the Yunnan Province, the Yunnan Province Scientist Workshop of Solar Physics, the Applied Basic Research of Yunnan Province (2019FB005, 202101AT070018), Yunnan Key Laboratory of Solar Physics and Space Exploration of Code