A Modified Spheromak Model Suitable for Coronal Mass Ejection Simulations
Talwinder Singh, Mehmet S. Yalim, Nikolai V. Pogorelov, Nat Gopalswamy
AA Modified Spheromak Model Suitable for CoronalMass Ejection Simulations
Talwinder Singh , Mehmet S. Yalim , Nikolai V. Pogorelov , and NatGopalswamy Department of Space Science, The University of Alabama in Huntsville, AL 35805, USA Center for Space Plasma and Aeronomic Research, The University of Alabama in Huntsville, AL35805, USA NASA/Goddard Space Flight Center, Greenbelt, MD 20771, USA
Abstract
Coronal Mass Ejections (CMEs) are one of the primary drivers of extreme spaceweather. They are large eruptions of mass and magnetic field from the solar corona andcan travel the distance between Sun and Earth in half a day to a few days. Predictionsof CMEs at 1 Astronomical Unit (AU), in terms of both its arrival time and magneticfield configuration, are very important for predicting space weather. Magnetohydro-dynamic (MHD) modeling of CMEs, using flux-rope-based models is a promising toolfor achieving this goal. In this study, we present one such model for CME simulations,based on spheromak magnetic field configuration. We have modified the spheromaksolution to allow for independent input of poloidal and toroidal fluxes. The moti-vation for this is a possibility to estimate these fluxes from solar magnetograms andextreme ultraviolet (EUV) data from a number of different approaches. We estimatethe poloidal flux of CME using post eruption arcades (PEAs) and toroidal flux fromthe coronal dimming. In this modified spheromak, we also have an option to controlthe helicity sign of flux ropes, which can be derived from the solar disk magnetogramsusing the magnetic tongue approach. We demonstate the applicability of this modelby simulating the 12 July 2012 CME in the solar corona. a r X i v : . [ a s t r o - ph . S R ] F e b Introduction
Coronal Mass Ejections (CMEs) are one of the most violent events in our solar system. Thetotal energy released in these events can range between 10 to 10 Joules [Vourlidas et al.,2002]. A CME can be ejected with speeds ranging between 20 km/s and 3500 km/s. CMEshave an average speed of 300 km/s during solar minima and 500 km/s during solar maxima[Yashiro et al., 2004]. Thomson scattering of white light from CME electrons allows us tocalculate the total mass of CMEs, which has been found to be between 10 and 4 × g,with an average of 10 g [Gopalswamy , 2010]. One of the important CME features is theirmagnetic flux rope structure, which is primarily responsible for the CME’s geoeffectiveness.Particularly, a CME is more geoeffective if its flux rope has a large negative Bz component ofmagnetic field at Earth. This is because of the favorable conditions for magnetic reconnectionbetween the flux-rope field and Earth’s magnetic field in the day-side magnetosphere.CME predictions at 1 AU remain an area of active research. MHD simulations are clearlyof importance for achieving better accuracy as compared with the simple empirical models[e.g. Vandas et al., 1996, Brueckner et al., 1998, Gopalswamy et al., 2001, 2005, Wang etal., 2002, Manoharan et al., 2004]. Many case studies have been done using MHD modelsthat show reasonable agreement between the simulated and observed properties of CMEs[Manchester et al., 2004, Jin et al., 2017, Singh et al. , 2018, 2019, and references therein].Current MHD CME models are broadly divided into two categories: (1) over-pressuredplasmoid models, such as the blob model [e.g., see, Chane et al., 2005, Odstrcil & Pizzo,1999] and (2) flux-rope-based models, such as the Titov–Demoulin model [Titov & Demoulin,1999], the Gibson–Low (GL) [Gibson & Low, 1998] model, and their variations. Since themagnetic flux rope of a CME is primarily responsible for its geoeffectiveness, flux-rope-basedmodels are clearly more realistic and promising for space weather predictions.Accurate predictions at 1 AU are impossible without constraining a CME model withobservations. The parameters that must be constrained during a CME eruption are: 1)CME speed, 2) direction, 3) orientation (the tilt angle between the horizontal and the axisof the flux rope at its apex), 4) poloidal flux, 5) toroidal flux, 6) helicity sign, and 7) mass.In this study, we have modified the force balanced spheromak solution, so that the above-mentioned CME properties can be constrained substantially in simulations. We also discusssome of the existing methods being used to derive these parameters from the observations.In Sec. 2, we present the data-driven solar wind model and the modified spheromak model.Section 3 describes various data and methods being used to find CME parameters. In Sec.4, we present an example of the application of our simulation model. Our conclusions arepresented in Sec. 5. 1 Simulation models
In this study, we perform the CME simulations in two steps. First, we create a solar windbackground from 1.03 R (cid:12) to 30 R (cid:12) using solar synoptic magnetograms. Then the flux ropemodel is inserted into the domain and it erupts as a CME due to pressure imbalance. Weuse Multi Scale Fluid Kinetic Simulation Suite (MS-FLUKSS, Pogorelov et al. [2014, 2017]).MS-FLUKSS is a highly parallelized code that can be used for MHD treatment of plasmaand fluid or kinetic treatment of neutral hydrogen atoms. Both the solar wind model andthe flux rope model are described in the following subsections. In this study we use the global MHD solar corona model [Yalim et al., 2017] driven byradial synoptic maps from the Solar Dynamics Observatory’s Helioseismic and MagneticImager (SDO-HMI). The ideal MHD equations are solved in the solar co-rotating frame, withvolumetric heating terms providing the required acceleration to the solar wind [Nakamizo etal., 2009]. These volumetric heating terms consist of an exponential heating function thattakes the expansion factor into account and a Spitzer-type thermal conduction term thatconducts heat along the magnetic field lines. The initial distribution of the magnetic field isfound using the Potential Field Source Surface (PFSS) model [Toth et al., 2011]. The rest ofthe initial plasma parameters are computed using Parker’s 1D isothermal solar wind solution[Parker, 1958]. We describe the implementation of boundary conditions in more detail inSec. 4.
Here, we describe a flux rope model based on the spheromak solution in which magneticforces are balanced by plasma pressure gradient forces to create a spherical spheromak inequilibrium [Gibson & Low, 1998, Lites, 1995]. The ∇ · B = 0 condition is specifically takeninto account in this procedure. The magnetic field morphology thus achieved can be seenin Fig. 9 of Lites [1995], Fig. 4 in Gibson & Low [1998] and Fig. 1 in Singh et al. [2018].The analytical solution for the spheromak model, as given in Apendix B2 of Gibson & Low[1998] , is: (cid:126)b = 1 r sin θ (cid:16) r ∂A∂θ ˆ r − ∂A∂r ˆ θ + α A ˆ φ (cid:17) , (1) A = 4 πa α (cid:104) r g ( α r ) g ( α r ) − r (cid:105) sin θ, (2)2 ( α r ) = sin( α r ) α r − cos( α r ) , (3)with α and r related as α r = 5 . J / .Here, r is the spheromak radius. The magnetic field strength in the spheromak is controlledby the parameter a . The plasma pressure in a spheromak is given by P = a A . The originof the spherical coordinate system r , θ and φ is placed at the spheromak center. We canperform coordinate transformations to shift this sphere to some off-center position and alsorotate it. Gibson & Low [1998] built on this solution by including a stretching parameterthat can turn a spherical torus into a tear drop shape. Several studies have been done toshow the applicability of this model to simulate flux-rope-driven CMEs [e.g. Manchester etal., 2004, Lugaz et al., 2005, Jin et al., 2017, Singh et al. , 2018]. When the force-balancedflux-rope model is superimposed with the background solar wind, the pressure imbalanceinside and outside the flux rope results in its eruption.There are two major drawbacks of using this model for CME simulations. Firstly, thepoloidal and toroidal magnetic fluxes cannot be controlled independently in this model,since there is only one parameter a that controls the magnetic field strength of the fluxrope. For different CME sizes and magnetic strength parameters, we find that the poloidaland toroidal fluxes do not differ more than 10% from each other. This is not necessarily truein an actual CME, so using the correct magnetic fluxes in the model is a key requirementfor B z -prediction at 1 AU. Secondly, the spheromak solution is unable to control the helicitysign of the flux rope. The helicity sign defines the direction of magnetic field line winding inthe flux rope. The change in its sign can result in a completely different magnetic field at 1AU.To address these shortcomings in the original spheromak solution, we propose to modifyit by introducing two extra parameters γ and δ in Eq. 1 so that (cid:126)b = 1 r sin θ (cid:16) γ r ∂A∂θ ˆ r − γ ∂A∂r ˆ θ + δα A ˆ φ (cid:17) This new solution is no longer in the force-balance condition. This means that this modelcannot be used to simulate pre-eruption, force-balanced flux ropes and their initiation phase.This model, however, can be readily used to simulate flux ropes of erupting CMEs, which arealready in a force-imbalance condition. It should be noted that erupted flux ropes of CMEsare considerably different from pre-eruptive flux ropes. This is because magnetic reconnectionoccurring during an eruption modifies the poloidal flux of a CME considerably [Longcopeet al. , 2007, Qiu et al., 2007, Gopalswamy et al., 2018]. Using such a force imbalancedflux rope model not only makes the simulation more robust, it also facilitates the input3f observed magnetic flux into it. Force–imbalanced models have been successfully used inprevious works as well. E.g., Manchester et al. [2004] modified a force balanced Gibson-Lowflux rope by reducing its density by 20% before superimposing it on the background solarwind. They further modified density and pressure inside the flux rope such that they do notdrop below 25% of the background values. These modifications result in force imbalance inthe initial flux rope itself. Similarly, Lugaz et al. [2007] used a modified Titov-Demoulin fluxrope by removing the strapping magnetic field lines from flux rope that were keeping it inforce balanced condition. Thus, their initial flux rope was in force–imbalanced condition aswell.Our modified flux rope still satisfies the ∇ · B = 0 condition. Plasma density insidethe flux rope is assumed constant initially, and it depends on the mass of the simulatedCME. Now, in this model, the toroidal flux is proportional to a while the poloidal flux isproportional to the product γa . This gives us an independent control over the two fluxesby varying a and γ . We have explained the method of deriving the poloidal and toroidalflux of a flux rope in Appendix A. Figure 1 shows the effect of varying γ on the magneticconfiguration of the flux rope. The poloidal flux increases proportionally to γ while thetoroidal flux remains unchanged in all cases. We show only a few magnetic field lines here todemonstrate the increase in the number of turns in response to an increasing poloidal flux.The translucent slice is colored by plasma density. Since we initially assume constant densityinside the flux rope when it is added to the background solar wind, one can see the flux ropeedge as an enhanced density region. In this example, the size parameter r is 1 R (cid:12) . The fluxrope is shifted by r = 1 R (cid:12) from the Sun center. This means that we are using roughly onlyone half of the spheromak in our simulation. This geometry roughly resembles a bent–tubeflux rope. Therefore, the force responsible for its eruption will be the Lorentz hoop force,arising due to the curvature of the flux rope. This force is at least partially responsible forthe propagation of the eruption in actual CMEs [Chen , 2017, Green et al., 2018]. In thisexample, we have kept the flux rope in the direction S12W06 while the orientation angle was53 degrees with respect to the solar equator. All this can be done easily by transformingthe spheromak solution from a local coordinate system to a shifted and rotated coordinatesystem.The helicity sign of the flux rope can be controlled by δ . The only values it can takeare ±
1. The plus sign will result in a positive helicity flux rope. Using δ = − left to right : flux rope with γ = 1, 2, and 4 respectively. The yellow sphererepresents the Sun, with magnetic field lines given by arrowed black lines. The slice throughthe flux rope is colored by the plasma density.Figure 2: Flux ropes with positive ( left panel ) and negative ( right panel ) helicity.The plasma β is very low in flux ropes. Therefore, we do not include the plasma pressureinto our model. The magnetic pressure dominates and is primarily responsible for CMEeruption. However, we observe that introducing thermal pressure comparable to magneticpressure in the flux rope can significantly change its eruption speed. This is equivalent toincreasing the plasma energy density in the flux rope. This can be used to constrain simulatedCME’s speed. We can also specify the CME mass uniformly distributed throughout itsvolume. The stretching operation, similar to the one described by Gibson & Low [1998], canstill be applied to this spheromak to convert the shape from sphere to a tear–drop shape.This shape conversion can bring the two legs of the flux rope closer to each other, thusreducing the width of the flux rope. This makes it possible to match the initial shape of oursimulated CME with the observed CME, as shown in Singh et al. [2018].5 CME observations
In this section we discuss the measurable properties of CMEs, which can be used to con-strain our CME model. We will briefly discuss how each quantity can be derived fromobservations. We use data from Solar Terrestrial Relations Observatory (STEREO) andSolar and Heliospheric Observatory (SOHO) coronagraphs, SDO magnetograms and SDOextreme ultraviolet (EUV) images to get these CME properties.
Coronagraphs are one of the primary instruments used to study CME evolution in the corona.However, a single viewpoint image is not sufficient to resolve the 3-D structure of a CME. Toovercome these limitations, the twin STEREO spacecraft [Kaiser et al., 2008] were launchedto see the same CME from multiple viewpoints. The triangulation techniques can then beused to study the CME in 3D space. This helps us to remove the projection effects andget the true kinematic properties of a CME. One of the best models for utilizing the threeviewpoints of STEREO A&B and SOHO coronagraphs to find 3D CME evolution is theGraduated Cylindrical Shell (GCS) model [Thernisien , 2011]. This model fits a typicalCME shape, i.e., the curved front with conical legs, to observations, giving us an estimatefor the CME height, direction and orientation with respect to solar equator. Figure 3 showsthis model applied to a CME that erupted on July 12, 2012. Fitting this model for a timeseries can be used to estimate its speed from the height-time profile [see e.g. Hess & Zhang, 2014].
The poloidal magnetic flux of a CME can be measured from the reconnected flux either usingflare ribbons or Post Eruption Arcades (PEAs). Gopalswamy et al. [2017] showed that bothof these methods are highly correlated. The PEA method requires EUV data only at thetime when the PEA structure has fully matured, typically in the decay phase of the flare,whereas the flare ribbon method requires 1600 ˚A data throughout the flaring time. Since thePEA method is more robust and easy to implement, Gopalswamy et al. [2018] propose to useit in their Flux Rope from Eruption Data (FRED) model. When the overlaying magneticfield lines gets stretched and reconnected during a flux rope eruption, one half of the linesreconnect down in the active region forming PEAs and the other half contributes to thepoloidal flux of the flux rope. This process can be understood more clearly looking at theflux rope eruption cartoon in Fig. 1 of Klein [2017]. Gopalswamy et al. [2017] show that the6
Figure 3: (
Top panel, from left to right ) July 12, 2012 CME seen in STEREO B Cor2,SOHO C2 and STEREO A Cor2 coronagraphs respectively. (
Bottom panel ) The sameimages overlapped with GCS model.poloidal flux of the CME is half the unsigned flux in the area covered by the PEAs.The left panel of Fig. 4 shows the PEA for 12 July 2012 eruption in SDO AIA 94, 131, and193 ˚A composite data at 22:30 UT, when the PEA has fully matured, ∼ . × Mx.Including the green area in our calculations increases the poloidal flux to 1 . × Mx, a9% increase. This gives us an estimate of the subjective error possible due to the manualPEA area selection in this method. Gopalswamy et al. [2018] found the poloidal flux for thesame event to be 1 . × Mx using just the 193 ˚A data, a very similar result compared7ith the one we obtained using the composite image. Therefore, it may be sufficient to usejust the 193 ˚A data for this analysis.Figure 4: (
Left panel ) SDO AIA 94, 131, and 193 ˚A composite data on 12 July 2012 22:30UT. (
Right panel ): SDO HMI magnetogram on 12 July 2012 16:10 UT. In both the cases,the area enclosed by PEAs is shown with red contours. The poloidal flux of the eruptedCME is half the unsigned flux in this area. The area enclosed by the green contours hassome ambiguity on whether it is a part of the PEA or some pre-eruption loops. Includingthis area into the PEA calculation increases poloidal flux by 9%. The cutout size is 300 Mm ×
300 Mm.
The toroidal flux of a CME, also known as the axial flux, or the core flux in the literature,can be found if the footpoints of the erupting CME are found [see, e.g., Webb et al., 2000].Dissauer et al. [2018] describe a way to do this using the coronal dimming. When a fluxrope erupts, a clear dimming is seen in EUV images. This is due to the mass loss duringthe eruption. Dissauer et al. [2018] show that this dimming can be separated into twocategories, 1) core dimming and 2) secondary dimming using appropriate thresholds. Thesecore dimming regions show higher mass loss and are seen during the first 30 minutes of theeruption. Once the core dimming regions have been found, the toroidal flux of a CME isgiven by the unsigned average magnetic flux in the positive and negative polarity footpoints.In Fig. 5, we show the procedure we followed to find the toroidal flux in the 12 July 2012CME. We use 193 ˚A data in our analysis, which is shown to be most suitable to study coronaldimming by Dissauer et al. [2018]. A pre-event image is found by taking a pixel-by-pixelmedian of 10 images within a half hour before the eruption starts (left panel of Fig. 5).Then, during 1 hour after the eruption has started, we use 2-minute-cadence data to detect8he pixels where the logarithm of the ratio of any post-event image data and the pre-eventimage data falls below -0.19. These are the regions showing coronal dimming. We also recordthe lowest value of logarithmic ratios during the selected time interval among these pixels.This forms a minimum intensity logarithmic base ratio image which has non-zero values onlyin the coronal dimming pixels. This is shown in the middle panel of Fig. 5. Similarly, aminimum intensity base difference image is formed by taking differences at each pixel, ratherthan logarithmic ratios. Once these two minimum maps have been formed, we find a subsetof the pixels that reside in the core dimming regions using the following thresholds fromDissauer et al. [2018]. A = ¯ I BD − . σ BD B = ¯ I LBR − . σ LBR
Here, ¯ I BD and ¯ I LBR are the mean intensities over all pixels in the minimum-intensity loga-rithmic base ratio image and the minimum-intensity base difference image, respectively. σ values give the corresponding standard deviations. The pixels flagged as core dimming pixelsare shown in red color in the right panel of Fig. 5. We can see that one group (inside a violetcircle) represents the positive footpoint, whereas the other (inside a green circle) representsthe negative footpoint of the erupting flux rope. The unsigned flux in these core dimmingpixels can be used as an estimate for the toroidal flux. The toroidal flux calculated usingthese pixels is 2 . × Mx. By changing the thresholds A and B by ± . × Mx and 2 . × Mx. Therefore, the toroidal flux varies by ≈ ± ≈
15% of its poloidal flux. Such a CME cannot bereliably simulated by the original spheromak model because the poloidal and toroidal fluxesdo not differ by more than 10% from each other in this model when it is inserted near theSun as discussed in Sec. 2.2.
The helicity sign of a flux rope can be determined from the pre-eruptive active region mag-netic field configuration [Bothmer and Schwenn , 1998]. For example, Luoni et al. [2011]show how the magnetic tongues in the ARs can be used to estimate the helicity sign in theoverlaying flux ropes. The helicity sign is easy to determine if one can estimate the directionof magnetic field lines on the axis of the flux rope and the overlaying loops, which eventuallycontribute to the poloidal flux of a CME during an eruption, as discussed in Sec. 3.2. InFig. 6, we show how the helicity sign can be found if we know the flux rope footpoints andthe neutral line (NL) above which the flux rope exists in corona. This NL is typically in the9igure 5: (
Left panel ) SDO AIA 193 ˚A data showing the pre-eruption median image ofthe active region. (
Middle panel ) The minimum intensity logarithmic ratio image showingcoronal dimming after the eruption. ( right ) The core dimming pixels are plotted over theSDO HMI magnetogram. This data is for July 12, 2012 15:40 UT. The violet and greencircles enclose the areas containing the footpoints of the flux rope, predominantly in thepositive and negative flux regions, respectively. The flux in these core dimming regionsprovides an estimate for the toroidal flux inside a CME. (Cutout size is 1000 arcsec)middle of the PEA area. The method to find footpoints has already been discussed in Sec.3.3. If we follow the NL in the direction from a positive footpoint to a negative footpoint withthe thumb and curve the fingers in such a way that they follow the overlaying field, we cantell whether the helicity sign is positive or negative based on whether the right or left handfollows the lines properly. In the example shown in Fig. 6 for July 12, 2012 magnetogram,we know that the axial field is directed from the right to the left, since the positive footpointis on the right as seen in Fig. 5. This determines the direction of orange arrow in Fig. 6.The overlying field lines are shown in green color, while the direction is simply from positiveto negative polarity. We can see that this configuration is consistent with the right-handrule. Therefore, the helicity sign is positive in the erupting flux rope. The helicity sign ofthis event was found by Gopalswamy et al. [2018] using a different method.
Mass of a CME can be found using the CME brightness in the coronagraph images. Thebrightness is due to the Thomson scattering of photospheric light by the plasma electrons[Billings, 1966]. By integrating over the CME area and removing the projection effects usingmultiple coronagraph viewpoints, we can calculate the true mass of a CME [Colaninno,2009]. 10
Figure 6: Helicity sign estimates can be made by following the axial field lines (orange) withthe thumb and curving the fingers along overlying field lines (green). This example is for 12July 2012 CME and shows that the erupted flux rope will have positive helicity., since thefield lines are traced properly by right hand.
We will show now the applicability of our approach to an observed CME. We choose the July12, 2012 CME as an example. This CME erupted from AR 11520 at 15:54 UT. The locationof the AR was 15 degrees south and 1 degree west at the time of eruption. The reason forchoosing this event is that the source active region is close to solar disk center, thus enablingmore accurate estimation of the AR magnetic field strengths. At this time, the STEREOspacecraft were observing at nearly right angles to the Sun-Earth line. This increases theaccuracy of the GCS method, thus allowing us to accurately determine the speed, direction,and orientation of the CME. Moreover, this event has been studied in detail by many authors[see e.g. Gopalswamy et al. , 2013, Hess & Zhang , 2014, Shen et al., 2014, Gopalswamy etal., 2018], thus allowing us to verify the accuracy of our derived quantities.We determined the properties of this CME using the methods described in the previoussection.1. The CME speed was found to be 1210 km/s at the height of 15 R (cid:12) . This was foundby fitting a quadratic function to the height-time profile found using the GCS method.A linear fit to the height-time plot gives us a speed of 1265 km/s for this CME.11. The direction of the CME was found to be at 12 degrees south and 8 degrees west.This is a little different from the source active region location, thus the CME showeda slight deflection in the lower corona.3. The CME flux rope orientation, found using the GCS method, was 53 degrees withrespect to solar equator. These GCS results are consistent with the ones reported inGopalswamy et al. [2018].4. The poloidal flux of the CME was found to be between 1 . × Mx and 1 . × Mx using the PEA method.5. The toroidal flux of the CME was found to be between 1 . × Mx and 2 . × Mx using the coronal dimming method.6. As shown in Fig. 6, the CME flux rope has a positive helicity sign.7. The mass of the CME was found to be 1 . × g. This is the true mass of the CMEwith projection effects removed using the multiple viewpoints of STEREO.The CME simulation is carried out in two steps. First we create a solar coronal MHDbackground and then introduce the flux rope model in it. The coronal background is createdin a fully spherical simulation domain, which extends from 1.03 R (cid:12) to 30 R (cid:12) . The SDOHMI synoptic map for Carrington rotation 2125 is used at the inner boundary of the domainlocated just above the transition region at 1.03 R (cid:12) . We relax the initial PFSS magnetic fieldto obtain a steady state solution. We used the Total Variation Diminishing (TVD), finitevolume Rusanov scheme [Kulikovskii et al., 2001] to compute the numerical fluxes and theforward Euler scheme for time integration. In order to satisfy the solenoidal constraint, weuse the Powell et al. [1999] approach. All simulations are performed in the frame corotatingwith the Sun. MS-FLUKSS is built on Chombo library, which ensures a highly parallelizedimplementation of our numerical schemes [Pogorelov et al., 2017].At the inner boundary of the computational domain, at 1.03 R (cid:12) , we specify the dif-ferential rotation [Komm et al., 1993a] and meridional flow [Komm et al., 1993b] formulaefor determining the horizontal velocity components at the ghost cell centers. Density andtemperature are kept constant as n = 1 . × cm − and T = 1 . × K, respectively.The radial velocity component is imposed to be zero at the boundary surface. The radialmagnetic field component is imposed from the magnetogram data. The transverse magneticfield components are extrapolated from the domain to the physical ghost cells below theinner radial boundary. No boundary conditions are required at the outer boundary of the12omain, since it is located beyond the critical point, where the plasma flow is superfastmagnetosonic.We introduce the modified spheromak into the domain so that the sphere center rests atthe inner boundary. This ensures that only one side of the spheromak is introduced into thedomain, which resembles a flux rope with two legs (see Fig. 1). The flux rope is introducedso that it replaces the magnetic field in the background corona and the calculated mass ofthe CME is uniformly distributed over it. We also adjust the parameters a and γ so thatthe poloidal and toroidal fluxes introduced into the model match the observations. We stillhave one free parameter, namely the radius of the spheromak, r . This parameter can beconstrained by the speed of the erupting CME. We find that the speed of eruption dependsinversely on the radius of the spheromak, keeping poloidal and toroidal fluxes constant. Thiscan be easily understood considering the source of flux rope eruption is pressure imbalance.If we want to have same flux in a smaller flux rope, we need to increase the magnetic fieldstrength in it. This increases the magnetic pressure inside the flux rope. Therefore, when thesolution is evolving in time, smaller flux ropes erupt at greater speeds. Singh et al. [2019]reported that the speed of this eruption also depends inversely on the magnetic pressurein the background solar wind in the region in the direction of CME propagation. We keepthe initial spheromak in the direction and with the orientation found using GCS method.By varying the value of r between 0.7 and 1.0 R (cid:12) , we find that r = 0 . R (cid:12) results inthe speed of the CME as 1200 km/s at 15 R (cid:12) , which is similar to the one calculated fromobservations. For this value of r , the a was 0.224 Gauss/ R (cid:12) and γ was 10.3 for the poloidaland toroidal fluxes to be 1 . × Mx and 2 . × Mx, respectively, which is the averageof the observed range. The parameter δ was kept as +1, since the flux rope was observedto have positive helicity. In Fig. 7, we show synthetic white-light images to compare thesimulated CME shape with the one observed 1.5 hours after eruption. Since we started witha uniform-density flux rope, the CME core is not seen clearly. We do, however, see the brightfront and the cavity of this CME. The overall shape of the simulated CME also agrees verywell with the observations. This also implies that the initial force distribution in our model isrealistic, since a non realistic force distribution would have resulted in considerably differentCME shapes compared to observations. Figure 8 shows our simulation results for this CME.The CME evolution is shown through snapshots at different time steps. The simulatedCME erupts with the proper speed, direction, orientation and magnetic field properties, akey requirement if we want to use MHD modeling for CME predictions.13 Figure 7: (
Top ) CME as seen by STEREO A (left) and B (right) COR 2 at 12-July-201217:54 UT. (
Bottom ) Synthetic white light images of the simulated CME at same height, withobserver fixed at the location of STEREO A (left) and B (right). These images are created byusing the ratio of line-of-sight-integrated brightness of the post-event and pre-event images.The color map is based on this ratio.Figure 8: (
From left to right ) The colorplot of the radial speed of simulated CME, 17, 34and 51 minutes after its initial insertion. 14
Conclusions
In this study, we modify the spheromak flux rope model so that the poloidal and toroidalfluxes in it can be controlled independently. The motivation for this is the possibility ofdetermination of both these fluxes independently from observations. We show how thismodel can be used for MHD simulations of CMEs. The pressure imbalance between theflux rope and the surrounding solar wind background leads to its eruption, resembling thecharacteristics of a CME. In the example shown in this paper, we were able to constrain thespeed, direction, orientation, and magnetic properties of the flux rope. We believe that thisapproach can become a viable option for predicting CMEs, not only in their arrival time,but in their magnetic field properties at 1 AU as well.TS acknowledges the graduate student support from NASA Earth and Space ScienceFellowship. The authors acknowledge the support from the UAH IIDR grant 733033. Thiswork is partly supported by the PSP mission through the UAH-SAO agreement SV4-84017.We also acknowledge NSF PRAC award OAC-1811176 and related computer resources fromthe Blue Waters sustained-petascale computing project. Supercomputer allocations werealso provided on SGI Pleiades by NASA High-End Computing Program award SMD-16-7570 and on Stampede2 by NSF XSEDE project MCA07S033. NG was supported in partby NASA’s LWS TR&T program.This work utilizes data from
SOHO which is a project of international cooperation be-tween ESA and NASA. The HMI data have been used courtesy of NASA/
SDO and HMIscience teams. The
STEREO /SECCHI data used here were produced by an internationalconsortium of the Naval Research Laboratory (USA), Lockheed Martin Solar and Astro-physics Lab (USA), NASA Goddard Space Flight Center (USA), Rutherford Appleton Lab-oratory (UK), University of Birmingham (UK), Max-Planck-Institute for Solar System Re-search (Germany), Centre Spatiale de Li`ege (Belgium), Institut d’Optique Th´eorique etAppliqu´ee (France), and Institut d’Astrophysique Spatiale (France). This work uses SOHOCME catalog which is generated and maintained at the CDAW Data Center by NASA andThe Catholic University of America in cooperation with the Naval Research Laboratory.
A Calculating the poloidal and toroidal flux of a fluxrope
In Fig. 9, we show a flux rope anchored on the Sun, with its curved axis in the x - y plane.We shaded the areas used to calculate the magnetic fluxes. The blue region is in the x - y plane and represents the area where B z >
0, i.e., field lines are coming out of the plane in15his region. The red region represents the cross section of the flux rope and belongs to the x - z plane. The poloidal flux of a flux rope in this configuration isΦ p = (cid:90) Blue area B z dA. The toroidal flux is Φ t = (cid:90) Red area B y dA. Figure 9: A flux rope with its curved axis in the z = 0 plane is shown anchored to the Sun.The blue shaded region is in the z = 0 plane and represents the area where B z >