A Time-dependent Self-similar Reconstruction of Solar Coronal Mass Ejections Based on the Gibson–Low Model

The analytic Gibson–Low (GL) model, a time-dependent self-similar solution of the magnetohydrodynamics, is first used to directly reconstruct a coronal mass ejection (CME) through the method of forward modeling in this study. A systematic description of the GL model is presented at the beginning, and a set of parameters is introduced to define the model. Then a CME on 2011 March 7 is reconstructed by fitting of GL (FGL) of the multi-viewpoint and time series observations. The first step of FGL is the initialization of the location and orientation of the GL using the information of the CME source region. The second step is to fit the parameters of size, shape, velocity, and strength of the magnetic field of the GL to the observations of coronagraphs at 20:24 and 20:39. The GL at 20:54 and 3 R ⊙ is generated through the theory of self-similar expansion respectively. Comparisons between the synthetic images of the GL and the real observations of the CME prove the performance of FGL that the reconstructions well match the observations.


Introduction
Coronal mass ejections (CMEs) are major sources of space weather Temmer 2021), which plays an important role in the high-tech modern society. Since the first space-borne coronagraph on board the Orbiting Solar Observatory (OSO-7; Koomen et al. 1975) discovered CMEs in the 1970s, a series of coronagraphs from space and the ground contributed to the studies of CMEs. The coronagraphs record the line-of-sight (LOS) integrals of Thomson scattering from the free electrons of optically thin CMEs (Billings 1966). These white-light observations of CMEs greatly help researchers understand abundant properties of CMEs, such as morphology, mass, velocity, energy, and frequency of occurrence (Webb & Howard 2012).
An obvious limit of the coronagraph observations is the LOS integrals because the results of these observations are twodimensional (2D) images of the three-dimensional (3D) CMEs. It is necessary to reconstruct the 2D images to obtain the more realistic 3D structure of CMEs. Various methods (Mierla et al. 2010;Thernisien et al. 2011) have been developed to achieve 3D reconstructions, and one category of these methods is called forward modeling . The forward modeling assumes that CMEs can be parameterized to a specific model as a cone (Fisher & Munro 1984;Zhao et al. 2002;Xie et al. 2004;Xue et al. 2005;Zhao 2005Zhao , 2008Cremades & Bothmer 2005;Michalek 2006;Na et al. 2013;Zhang 2021), flux rope (Chen et al. 2000;Wood & Howard 2009;Wood et al. 2011), or graduated cylindrical shell (GCS; Thernisien et al. 2006Thernisien et al. , 2009Thernisien 2011;Volpes & Bothmer 2015;Kay & Gopalswamy 2017;Chi et al. 2020), etc. The outline or whitelight brightness of one such model can be projected into the plane of sky (POS) of a coronagraph to make a synthetic CME image. Some distribution of the electron density should be located on the boundary of the model CME to produce the white-light brightness. Then, some parameters of the model, like propagation direction, angular width, and height of the leading edge, are modified to match the observation of the coronagraph. The new set of parameters defines a new distribution of the CME outline or white-light brightness that is used to match the observation for the next iteration. The above fitting process is repeated until the difference between the model and observation is small enough to complete the reconstruction.
So far, the models used in the forward modeling have mainly focused on the outer shape of the CMEs and not on the internal structure. In order to obtain the internal structure, more realistic models of CMEs should be used in the forward modeling. One such model may be the one developed by Gibson & Low (1998); hereafter, this model is named GL for convenience. The GL is an analytic model of magnetohydrodynamics (MHD) that is able to construct the 3D structure of magnetized plasma, especially for three-part CMEs (Illing & Hundhausen 1985;Vourlidas et al. 2013). The GL model has been applied to numerous studies of CMEs using numerical simulations (Manchester et al. 2004a(Manchester et al. , 2004b(Manchester et al. , 2008(Manchester et al. , 2014a(Manchester et al. , 2014bLugaz et al. 2005;Jin et al. 2016Jin et al. , 2017aJin et al. , 2017bJin et al. , 2018Borovikov et al. 2017). It is convenient to employ the GL in the simulation of CMEs using the Space Weather Modeling Framework (SWMF; Tóth et al. 2005Tóth et al. , 2012Gombosi et al. 2021). Singh et al. (2018Singh et al. ( , 2019 employed the GCS to derive the geometry parameters of the GL. Recently, Singh et al. (2020aSingh et al. ( , 2020b further modified the original GL to make it possible to independently control the poloidal and toroidal fluxes. The GL was also used to verify the accuracy of other methods of CME reconstruction (Dai et al. 2020).
However, the GL has not been applied to directly reconstruct the CMEs by forward modeling. In this study, the GL replaces the models previously used in the forward modeling, such as cone or GCS, to reconstruct the CMEs including both the outer shape and the internal magnetized plasma. The details of the GL are presented in Section 2. Section 3 shows an example of CME reconstruction employing the GL in forward modeling, and a summary is given in Section 4.

The GL Model
The CMEs can be described by time-dependent 3D MHD equations: which are usually not tractable. Low (1982aLow ( , 1982bLow ( , 1984aLow ( , 1984b) illustrated a mathematical process to generate the analytical solutions of these equations under the following assumptions: 1. The magnetic field is assumed to be symmetric about an axis and does not depend on the azimuth angle in spherical coordinates (Chandrasekhar 1956;Chandrasekhar & Prendergast 1956). 2. The magnetized plasma only flows in the radial direction. 3. The radial flow propagates in a self-similar way. 4. The radial flow is isentropic and adiabatic with the polytropic index γ = 4/3, which makes the entropy-like quantity s = pρ −4/3 conservative in Equation (5) under the self-similar assumption. Gibson & Low (1998) pointed out that the conservative condition s = pρ −4/3 physically closes Equations (1) to (5) under the above assumptions. Because s is an arbitrarily prescribed value, the system of MHD equations is not yet mathematically closed, and Gibson & Low (1998) found a stretching transformation in the radial direction to mathematically close the system. The other basic transformation for the GL, also in the radial direction, is presented in Low (1982aLow ( , 1982bLow ( , 1984aLow ( , 1984b) based on the self-similar theory.
Before the introduction of these two core transformations of GL, it is necessary to show the basic magnetohydrostatic (MHS) system: where the magnetic field b and pressure Π can be explicitly solved inside a sphere with r 0 radius under the axial symmetry assumption (Prendergast 1956;Lites et al. 1995), and the analytic solutions are sin , 10 The ball-like system of magnetized plasma described by Equations (6) and (7) is represented by a circle with radius r 0 in the top-left panel of Figure 1. The center of the circle is r 1 away from the center of the Sun, which is also the origin of the coordinates. For this system, the Lorentz force and the gradient of thermal pressure are in equilibrium as shown in Equation (6). α 0 is restricted by the Bessel function J 5/2 (α 0 r 0 ) = 0, and a 1 is a free parameter that modulates both b and Π. The gravitational force is not considered in Equation (6), and in order to obtain the plasma density that is included in the gravity term, Gibson & Low (1998) constructed the following transformation in the radial direction: to solve the following MHS equations including gravity: The exact solutions of the magnetic field H, pressure P, and density D can be derived from b and Π, which are the solutions of Equations (6) and (7): where the radial component Λ is replaced by ζ. The arbitrary constant k is set to 1 in Equation (12) because it just scales the entire system, as pointed out by Gibson & Low (1998). After the stretching transformation ζ = Λ − a, the ball-like system turns to a flux rope (FR) with the shape of a teardrop as shown in the top-left panel of Figure  is represented by the red teardrop whose vertex coincides with the Sun center. The second self-similar radial transformation (Low 1982a(Low , 1982b(Low , 1984a(Low , 1984b is ss 2 ss ss which turns the MHS system of Equations (13) and (14) into a time-dependent MHD system defined by Equations (1) and (2).
Equations (22) and (23) combine the radial component and time and give the self-similar radial velocity: ss ss Then, the time-dependent solutions of the MHD equations, B, p, and ρ, can be derived from the MHS solutions, H, P, and D: The arbitrary constant α in Equations (13), (20), and (23) determines an additional force αζ originating from the selfsimilar transformation. In the similarity space (ζ, θ, f), αζ may balance the gravity, gradient of pressure, and Lorentz force for the MHS system. For the MHD system in real time and space (t, r, θ, f), the FR may expand away from the Sun with a motion of acceleration (α > 0), deceleration (α < 0), or inertia (α = 0). In the top-right panel of Figure 1, FR 1 expands to FR 2 in a self-similar way from the time t 1 to t 2 . A series of radial expansions of the FR is shown in the bottom-right panel of Figure 1 rendered with different colors. In this work, the values of α and η are set to 0 and t 1 scl 2 , respectively, which are the same as done in Gibson & Low (1998). Then the temporal function Φ ss defined in Equation (23) is equal to t/t scl , and we have ss ss ss scl where the units of t and ζ are seconds and kilometers, and t scl = 3600 s. Speeds of different parts of the FR are of the order of a few hundred or thousand kilometers per second because of the timescale normalization by setting h = t 1 scl 2 , as pointed out by Gibson & Low (1998). Equation (29) means that each particle in the self-similar radial flow moves at a constant speed because α = 0, and different particles have different speeds according to their own value of ζ.
The parameters of GL are summarized as follows: 1. Longitude and latitude of the FR. 2. Orientation of the symmetry axis of the ball-like system before the two radial transformations. 3. r 0 in Equations (10) and (21), r 1 in Equation (21), and a in Equations (12) and (19). (24). 5. a 1 in Equations (9) and (10).

v in Equation
In the next section, the values of the above parameters will be specified one by one through the fitting of GL (FGL) to real observations of a CME. For simplicity, the helicity of the GL is fixed to be positive in this paper, the same as done in Jin et al. (2017b). In future works using the FGL, the sign of the helicity can be determined by various methods (Luoni et al. 2011;Jin et al. 2017b;Gopalswamy et al. 2018;Singh et al. 2020b). The source code of the GL employed in this study can be found in the open-source version of the SWMF 3 . The exact location of the CME source region and the orientation of the symmetry axis will be extracted from the HMI magnetogram. First, it is necessary to remove the projection effect for AR 11164 located at the northwest corner of the Sun. An IDL routine WCS_2D_SIMULATE (Thompson 2006;Thompson & Wei 2010) in SolarSoft (Freeland & Handy 1998) is employed to simulate an image where AR 11164 is moved to be located at 0°longitude and 0°latitude as shown in the bottom-left panel of Figure 2. The region surrounded by a white square is enlarged and shown in the bottom-right panel. The positive and negative areas of AR 11164 are selected by two white rectangles. Point C 1 is the weighted average location of the positive magnetic field inside the white rectangle, and C 2 is the weighted average location of the negative magnetic field. The center of AR 11164 is then defined as the midpoint of the line C 1 C 2 , which is marked with C 3 located at 45°. 73 longitude and 28°. 39 latitude in the Stonyhurst-Heliographic coordinates. This is the initial location of the FR and will be adjusted according to the fitting to CME observations. Finally, the symmetry axis of the ball-like system is parallel to C 4 C 5 , which is the perpendicular bisector of C 1 C 2 . The symmetry axis is 39°. 65 with respect to the solar equator.

The FGL and Its Application
AR 11164 was marked by AR0 in Figure 5 of Jin et al. (2017b) and by a red ellipse in Figure 4 of Singh et al. (2018). The same extent of the source region is selected to calculate the location and orientation of the FR in this work. As pointed out by Jin et al. (2017b), the FR orientation is set to be parallel to the neutral line, and the FR location is determined by the intersection between the neutral line and the line connecting the weighted average centers of the negative and positive magnetic fields, C 1 C 2 . In this paper, the midpoint C 3 is directly set to be the location of the FR, and the orientation of the symmetry axis of the ball-like system is parallel to C 4 C 5 instead of the neutral line. The location of the FR is 155°Carrington longitude (46°S tonyhurst longitude) and 27°latitude in Jin et al. (2017a), which are close to 45°. 73 Stonyhurst longitude and 28°.39 latitude in this work.
The second step of the FGL focuses on the CME observations from the coronagraph. As illustrated in the topright panel of Figure 1, L 3 and L 4 on the central axis are the true apexes of FR 1 and FR 2 respectively. However, we usually cannot see L 3 and L 4 in the CME image unless the CME is a perfect limb event whose central axis lies in the POS. Actually we usually can see L 1 and L 2 , whose projections in the POS are P 1 and P 2 as shown in Figures 1 and 3. In other words, P 1 and P 2 are the projections of the visual apexes of the CME. For the same reason, the central axis exhibited with a black line in Figure 3 is the projection of the visual central axis.
We initially set r 0 = 1.0 R e , r 1 = 1.5 × r 0 = 1.5 R e , and ( ) =´=´-= a a r r 0.5 0.5 0.25 max 1 0 R e . An additional parameter, τ, is used to modify the size of the whole FR: where τ is set to be 2.8 initially. Now the 3D distribution of the FR in the self-similarity space (ζ, θ, f) is completely initialized because the exact location of the CME source region and the orientation of the symmetry axis have been specified from Figure 2 and the values of r 0 , r 1 , and a are also initialized. At this point, if the time variable, t, in Equation (28) is specified, an expanding FR in real time and space (t, r, θ, f) will be ready. For the time at 20:39 in the top-right panel of Figure 3, t 2 is initially set to be 9000 s (2.5 hours), and the corresponding time, t 1 , at 20:24 in the top-left panel should be 9000 − 60 × 15 = 8100 s. With this time setting, the corresponding outlines of FR 1 at t 1 and FR 2 at t 2 are projected into the top-left and top-right panels of Figure 3, respectively, and we can make comparisons to the real CME observations. The size of the FR is obviously larger than the real CME, and the parameters of the FR, especially τ, should be modified to fit the real observations better. The apex of the FR marked with a black plus (P 1 or P 2 ) should coincide with the apex of the real CME marked with a red plus if we have the proper value of τ as the case shown in the bottom panel of Figure 3. Because the speed of the self-similar radial flow just depends on the value of ζ as defined in Equation (29), we have where v fr is the speed of the FR at L 1 and L 2 , which have the same radial distance, ζ fr , in the self-similarity space, and P 1 and P 2 are the projections of L 1 and L 2 , respectively, in the POS as shown in Figure 1. P 1 and P 2 are also displayed in the top panel of Figure 3 marked with two black pluses. We assume that the apexes of the real CME marked with the two red pluses are the projections of ¢ L 1 and ¢ L 2 , which are located along the same radial direction as L 1 and L 2 . Thus the FR and real CME have the same initial radial distance, ζ ini , in this direction for L 1 , L 2 , ¢ L 1 and ¢ L 2 . v cme is the speed of the real CME at ¢ L 1 and ¢ L 2 who have the same radial distance, ζ cme , in the self-similarity space. If τ fr = τ cme , ¢ L 1 will coincide with L 1 , and ¢ L 2 will coincide with L 2 . Now the proper value of τ cme can be derived from Equation (31): where τ fr is already set to be 2.8 as an initial value of τ for the FR. v fr = 1065.93 km s −1 , which can be calculated through Equation (29). v cme should be calculated from the distance Δd cme between the two red pluses, the time difference Δt cme between 8:24 and 8:39, and the angle λ between the radial direction and the POS as because the CME is assumed to be inertial (α = 0), the same as the FR. For this event, Δd cme = 1.45 × 10 6 km, Δt cme = 60 × 15 = 900 s, and λ = 29°.97. Then, v cme = 1859.76 km s −1 , and so τ cme = 4.89. Thus r 0 = 4.89 × 1.0 = 4.89, r 1 =4.89 × 1.5 = 7.335, and a = 4.89 × 0.25 = 1.2225. The latitude of the FR is also modified from 28°.39 to 38°.39 to fit the CME observations better. Finally, we obtain the projections of FR outlines in the POS as shown in the bottom panel of Figure 3. As previously said, t 2 is initially set to be 9000 s and now it should be modified to match the value of τ cme . If τ fr = τ cme = 4.89, then v fr = v cme = 1859.76 km s −1 . The height of L 2 is h 2 and is equal to 5607152 km, then t 2 = h 2 /v cme = 3015 s, and t 1 = t 2 − 60 × 15 = 2115 s. The speed of the real apex of the CME, L 3 or L 4 , is equal to 2160 km s −1 .
The last undetermined parameter of the GL is the a 1 in Equations (9) and (10). Jin et al. (2017b) derived an equation to calculate a 1 based on the empirical relationships among the speed of the CME, the magnetic field of the source region, and the size of the FR. Singh et al. (2018) made a similar parametric study to calculate the a 1 . Unlike Equation (1) in Jin et al. (2017b), Equation (18) in Singh et al. (2018) includes a term of the simulated solar wind pressure instead of the observed magnetic field of the source region to calculate the a 1 .
In this paper, a direct way to obtain the value of a 1 is presented. The right-hand side of Equation (19) tells us that the density, D(ζ, θ, f), of the FR is composed of the square of the On the other hand, the density determines the mass (Vourlidas et al. 2010;Bein et al. 2013;Dai et al. 2014Dai et al. , 2015Pluta et al. 2019) of the FR: where M fr and M cme can be calculated by the SolarSoft routine SCC_CALC_CME_MASS based on the synthetic and observed white-light brightness of the FR and CME, respectively. a 1fr is an arbitrary value of a 1 for the FR, and a 1cme is the final value that fits the CME observation. The initial value of a 1fr is set to 0.001, and the corresponding value of M fr is equal to 1.27 × 10 16 grams at 20:39 while M cme = 6.64 × 10 15 grams. So we have a 1cme = 7.23 × 10 −4 according to Equation (36). The synthetic brightness of the FR when a 1 = a 1cme is shown in the bottom-right panel of Figure 4.  Singh et al. (2018) and the FGL in this work, the value of a 1 will be changed to 6.81 and 1.01 (at 20:39), respectively, according to the inversely proportional relationship with r 0 4 . Aside from the CME mass, the parameter α can also impact the calculation of a 1 . According to the relationship between α and D(ζ, θ, f) as shown in Equations (19) and (20), D(α > 0) <D(α = 0) and D(α < 0) > D(α = 0). In other words, M fr (α > 0)< M fr (α = 0) and M fr (α < 0) > M fr (α = 0). Finally, we have a 1cme (α > 0) > a 1cme (α = 0) and a 1cme (α < 0) < a 1cme (α = 0) based on Equation (36). That is to say, the value of a 1cme will be underestimated with the assumption of inertia (α = 0) if the CME is  actually in acceleration (α > 0), and the value of a 1cme will be overestimated with the assumption of inertia (α = 0) if the CME is actually in deceleration (α < 0). Thus the assumption of inertia (α = 0) may lead to the underestimation or overestimation of the magnetic field, pressure, and density, because they are proportional to a 1 , a 1 2 , and a 1 2 , respectively, as shown in Equations (8) to (20). In future studies, the more accurate value of α may be derived from the time sequence of coronagraph observations if the CMEs are accelerating or decelerating, and we may obtain a more accurate FGL reconstruction.
The results of FGL from the views of the Solar and Heliospheric Observatory (SOHO; Domingo et al. 1995) and :39 are shown in Figure 5. The C3, a coronagraph from the Large Angle Spectroscopic Coronagraph (LASCO; Brueckner et al. 1995) package on board SOHO, captures the CME at 20:30 while the SECCHI/ COR2 captures the same CME at 20:24 as shown in the first row of Figure 5. Because of this slightly different observation timing between C3 and COR2, we decide to directly fit the GL with the observations from COR2. It means that the outer shape of the FR may be slightly smaller than the observed CME in the field of view (FOV) of C3. The corresponding brightness images of the FR are shown in the second row. The outer shapes and brightness of the FR and CME at 20:39 are shown in the last two rows. It is able to adjust the parameters of GL based on the multi-viewpoint and time series observations of the CME until the fitting is good enough.
It is feasible to constrain the FR parameters r 1 and a with a fixed value of r 0 using the angular width of the CME seen in the coronagraph images. The width of the FR becomes narrower from the left to right panels in Figure 6 due to the larger stretching transformation, which is controlled by the values of r 1 and a as described in Equations (12) and (21).

Because
[ ] Î a a 0, max and =a r r max 1 0 , we should choose a combination of r 1 and a like the case in the first panel of Figure 6 if the width of the observed CME is large. On the other hand, we should choose a combination of r 1 and a like the case in the third panel if the width of the observed CME is small.
The results of the FGL and GCS reconstruction of this CME are shown together in Figure 7. The GCS with the values of longitude (45°.73), latitude (38°.39), and tilt angle (39°.65), the same as the FGL, is superposed on the observed CME in the first row. We can see that the outer shape and position of the FGL and GCS are comparable and can be mutually verified by comparison to the observed CME. In order to obtain a more consistent reconstruction between the FGL and GCS, the latitude of the GCS is reduced to 28°.39 in the middle row, and the tilt angle of the GCS is further increased to 90°in the bottom row. It is not feasible to modify the tilt angle in FGL according to the comparison between the outer shape of the FR and the observed CME because the aspect ratio of the current GL model is always equal to 1.0. So, the modification of the tilt angle should be carried out with the assistance of other techniques, like GCS, if the CME rotates about the radial direction during its outward propagation. The GCS can be achieved through the user-friendly IDL routine rtsccguicloud in the Solarsoft at stereo/secchi/idl/scraytrace.
The agreements of the outer shapes and brightness between the FR and the CME as shown in Figures 4 and 5 tell us that the self-similar expansion of the FR successfully captures the outward motion of the CME in the FOVs of coronagraphs. The self-similar expansions of the velocity, magnetic field, pressure, and density are described by Equations (24) to (27), respectively, and the two-dimensional distributions of these variables at 20:39 are shown in Figure 8. We can see a typical  Figure 9, respectively. As a time-dependent self-similar model, the GL can expand along the radial direction as shown in the bottom-right panel of Figure 1. In other words, we are able to obtain the reconstructed CME at any time during the self-similar expansion once the complete set of parameters is determined. For example, the FR at 20:54 is shown in the third column of Figure 9, which is derived from Equations (24) to (29) by setting t = t 3 = t 2 + 60 × 15 = 3015 + 900 = 3915 s. It means that the FR at 20:54 is the result of the self-similar expansion of the FR at 20:39 after 15 minutes.
An obvious self-similar expansion with a constant speed can be seen from the distributions of speed at 20:24, 20:39, and 20:54 as shown in the first row of Figure 9, which is constructed from Equation (24). The evolutions of the magnetic field, density, and pressure are shown in the second to fourth rows, which are constructed from Equations (25) to (27), respectively. In order to verify the time-dependent self-similar result of reconstruction at 20:54, the observed brightness of the CME and the synthetic brightness of the FR in the COR2 FOVs of STEREO-A and B are shown in Figure 10. We find that the agreements of the outer shapes and brightness between the CME and the FR at 20:54 are good enough, the same as the results at 20:24 and 20:39 as shown in Figures 4 and 5. Low (1982a) pointed out that the motions of CMEs are nonself-similar in the low corona below 2 R e . Cremades et al. (2020) analyzed multi-viewpoint data of CMEs and demonstrated that their initial expansions are considerably asymmetric and non-self-similar below 3 R e . It is not recommended to casually apply the FGL below this critical height because the CMEs may no longer propagate in a self-similar way and the FGL is not applicable. The initial FR, whose apex is located at 3 R e , is demonstrated in Figure 11 to show the typical structure of the twisted magnetic field. In order to obtain this FR, we need to calculate the time t ini , at which the apex of the FR reaches 3 R e . Because the speed of the CME apex v = 2610 km s −1 , then t ini = 3 × R e /v = 3 × 6.955 × 10 5 /2160 = 799 s.