4D Modeling of CME expansion and EUV dimming observed with STEREO/EUVI

This is the first attempt to model the kinematics of a CME launch and the resulting EUV dimming quantitatively with a self-consistent model. Our 4D-model assumes self-similar expansion of a spherical CME geometry that consists of a CME front with density compression and a cavity with density rarefaction, satisfying mass conservation of the total CME and swept-up corona. The model contains 14 free parameters and is fitted to the 2008 March 25 CME event observed with STEREO/A and B. Our model is able to reproduce the observed CME expansion and related EUV dimming during the initial phase from 18:30 UT to 19:00 UT. The CME kinematics can be characterized by a constant acceleration (i.e., a constant magnetic driving force). While the observations of EUVI/A are consistent with a spherical bubble geometry, we detect significant asymmetries and density inhomogeneities with EUVI/B. This new forward-modeling method demonstrates how the observed EUV dimming can be used to model physical parameters of the CME source region, the CME geometry, and CME kinematics.


Introduction
It has become a generally agreed concept that the EUV dimming observed during the onset of a coronal mass ejection (CME) manifests the coronal mass loss of the CME, and thus we basically expect a one-to-one correlation between the detections of CMEs and EUV dimmings, unless there exist special circumstances. For instance, the CME could originate behind the limb, in which case the EUV dimming is obscured, or the CME could start in the upper corona, where there is little EUV emission because of the gravitational stratification. The latter case would imply very low masses com-Correspondence to: Aschwanden (email: aschwanden@lmsal.com) pared with a CME that originates at the base of the corona, i.e., ≈ 10% at two thermal scale heights. However, there exists a case with an average CME mass that did not leave any footprints behind in EUV (Robbrecht et al. 2009). A statistical study on the simultaneous detection of EUV dimmings and CMEs has recently been performed by Bewsher et al. (2008). This study based on SOHO/CDS and LASCO data confirms a 55% association rate of dimming events with CMEs, and vice versa a 84% association rate of CMEs with dimming events. Some of the non-associated events may be subject to occultation, CME detection sensitivity, or incomplete temperature coverage in EUV and soft X-rays. Perhaps the CME-dimming association rate will reach 100% once the STEREO spacecraft arrive at a separation of 180 • and cover all equatorial latitudes of the Sun.
A number of studies have been carried out by using the detection of coronal dimming to identify CME source regions, focusing on transient coronal holes caused by filament eruptions (Rust 1983;Watanabe et al. 1992), EUV dimming at CME onsets (Harrison 1997;Aschwanden et al. 1999), soft X-ray dimming after CMEs (Sterling & Hudson 1997), soft X-ray dimming after a prominence eruption (Gopalswamy & Hanaoka 1998), simultaneous dimming in soft X-rays and EUV during CME launches (Zarro et al. 1999;Harrison & Lyons 2000;Harrison et al. 2003), determinations of CME masses from EUV dimming from spectroscopic data (Harrison & Lyons 2000;Harrison et al. 2003) or from EUV imaging data (Zhukov and Auchere 2004;Aschwanden et al. 2009b). All these studies support the conclusion that dimmings in the corona (either detected in EUV, soft X-rays, or both) are unmistakable signatures of CME launches, and thus can be used vice versa to identify the mutual phenomena.
In this study here we attempt for the first time to model the kinematics of a CME and the resulting EUV dimming quantitatively, which provides us unique physical parameters of the CME source region and on the CME kinematics in the initial acceleration phase. Aschwanden: CME Expansion and EUV Dimming

Model Assumptions
The dynamics of a CME can often be characterized by a rapid expansion of a magnetically unstable coronal volume that expands from the lower corona upward into the heliosphere. Different shapes have been used to approximately describe the 3D geometry of a CME, such as a spherical bubble, an ice-cone, a crescent, or a helical flux rope, which expand in a self-similar fashion and approximately maintain the aspect ratio in vertical and horizontal directions during the initial phase of the expansion. Here we develop a four-dimensional (4D=3D+T) model that describes the 3D evolution of the CME geometry in time (T) in terms of 4D electron density distributions n e (x, y, z, t) that allow us also to predict and forward-fit a corresponding EUV intensity image data cube I EUV (x, y, t) in an observed wavelength.
For the sake of simplicity we start in our model here with the simplest case, assuming: (1) spherical 3D geometry for the CME front and cavity; (2) self-similar expansion in time; (3) density compression in the CME front and adiabatic volume expansion in the CME cavity; (4) mass conservation for the sum of the CME front, cavity, and external coronal volume; (5) hydrostatic (gravitational stratification) or superhydrostatic density scale heights; (6) line-tying condition for the magnetic field at the CME base; and (7) a magnetic driving force that is constant during the time interval of the initial expansion phase. This scenario is consistent with the traditional characterization of a typical CME morphology in three parts, including a CME front (leading edge), a cavity, and a filament (although we do not model the filament part). The expanding CME bubble sweeps up the coronal plasma that appears as a bright rim at the observed "CME front" or leading edge. The interior of the CME bubble exhibits a rapid decrease in electron density due to the adiabatic expansion, which renders the inside of the CME bubble darker in EUV and appears as the observed "CME cavity". The assumption of adiabatic expansion implies no mass and energy exchange across the outer CME boundary, and thus is consistent with the assumption of a low plasma β-parameter in the corona with perfect magnetic confinement, while the CME bubble will become leaking in the outer corona and heliosphere, where the plasma β-parameter exceeds unity (not included in our model here).

Analytical Model
A spherical 3D geometry can be characterized by one single free parameter, the radius R of the sphere. The self-similar expansion maintains the spherical shape, so the boundary of the CME bubble can still be parameterized by a single timedependent radius R(t). The time-dependence of the CME expansion is controlled by magnetic forces, e.g., by a Lorentz force or hoop force. For sake of simplicity we assume a constant force during the initial phase of the CME expansion, The definition of the self-similar geometry of the CME model is depicted for four different times, consisting of a cyclidrical base and a spherical shell. The height of the centroid of the spherical CME volume is h(t), the outer radius of the CME sphere is R(t), and the inner radius of the CME front is r(t). These parameters increase quadratically with time. The circular footpoint area of the CME with radius r0 stays invariant during the self-similar expansion in order to satisfy the line-tying condition of the coronal magnetic field at the footpoints.
which corresponds to a constant acceleration a R and requires three free parameters (R 0 , v R , a R ) to characterize the radial CME expansion, where R 0 = R(t = t 0 ) is the initial radius at starting time t 0 , v R is the initial velocity and a R is the acceleration of the radial expansion. For the motion of the CME centroid at height h(t) we assume a similar quadratic parameterization, where h 0 = h(t = t 0 ) is the initial height at starting time t 0 , v h is the initial velocity and a h is the acceleration of the vertical motion. This parameterization is consistent with a fit to a theoretical MHD simulation of a breakout CME (Lynch et al. 2004) as well as with kinematic fits to observed CMEs (Byrne et al. 2009). Further we constrain the CME geometry with a cylindrical footpoint area of radius r 0 , which connects from the solar surface to the lowest part of the spherical CME bubble. In order to ensure magnetic line-tying at the footpoints, the CME bubble should always be located above the cyclidrical footpoint base, which requires that the initial height satisfies h 0 > r 0 and the acceleration constants are a h > a r . We visualize the model geometry in Fig. 1.
The time-invariant CME footprint area allows us a simple mass estimate of the CME from the cylindrical volume integrated over a vertical scale height, since the spherical CME bubble will eventually move to large heights with no additional mass gain (at time t n ≫ t 0 , see right-hand panel in Fig.1).
Assuming adiabatic expansion inside the CME cavity, the electron density in the confined plasma decreases reciprocally to the expanding volume, i.e., so it drops with the third power as a function of time from the initial value n 0 (of the average density inside the CME).
For the mass distribution inside the CME we distinguish between a compression region at the outer envelope, containing the CME front, and a rarefaction region in the inside, which is also called CME cavity. We define an average width w 0 of the CME front that is assumed to be approximately constant during the self-similar expansion of the CME. While the radius R(t) demarcates the outer radius of the CME front, we denote the inner radius of the CME front or the radius of the cavity with r(t), which has an initial value of The total volume V 0 of the CME is composed of a spherical volume with radius R(t) and the cylindrical volume beneath the CME with a vertical height of (h 0 − r 0 ), which has an initial volume value of V 0 , The volume of the CME front V F (t) is essentially the difference between the outer and inner sphere (neglecting the cylindrical base at the footpoint) while the volume V C of the cavity is, We have now to define the time-dependent densities in the CME, for both the CME front, which sweeps up plasma during its expansion, as well as for the CME cavity, which rarifies due to the adiabatic expansion. The total mass m E (t) of the plasma that is swept up from the external corona in a CME corresponds to the total CME volume V (t) minus the initial volume of the CME cavity, where m p is the mass of the hydrogen atom and < n E > is the electron density in the external corona averaged over the CME volume. The same mass has to be contained inside the volume V F of the CME front, Thus, mass conservation yields a ratio of the average electron density < n F > in the CME front and the average external density < n E > of This density ratio amounts to unity at the initial time, i.e., q n (t = t 0 ) = 1 and monotonically increases with time. The maximum value of the density jump in MHD shocks derived from the Rankine-Hugoniot relations (e.g., Priest 1982) is theoretically q n,max = 4. Fast CMEs are expected to be supersonic and will have a higher compression factor at the CME front than slower CMEs. Thus we keep the maximum compression factor q n,max as a free parameter, keeping in mind that physically meaningful solutions should be in the range of 1 < ∼ q n,max < ∼ 4. Since we prescribe both the width w 0 of the CME front as well as a maximum density compression ratio q n,max we obtain a definition of the critical value ρ(t) for the cavity radius r(t) when the prescribed maximum density compression q n,max is reached (using Eq. 11), which yields the critical radius ρ(t), Therefore, only plasma outside this critical radius ρ(t) can be compressed in the CME front, while the plasma inside this critical radius dilutes by adiabatic expansion and forms the cavity, yielding an average density ratio q n,cav inside the cavity (according to Eq. 3), Our numerical model of a spherical CME expansion has a total of 14 free parameters: 3 positional parameters (the heliographic coordinates (l, b) and height h 0 of the initial CME centroid position), 5 kinematic parameters (starting time t 0 , velocities v h , v R , accelerations a h , a R ), 2 geometric parameters (initial radius r 0 and width w 0 of the CME front), and 4 physical parameters (coronal base density n 0 , maximum density compression factor q n,max in the CME front, the mean coronal temperature T 0 at the observed wavelength filter), and a vertical density scale height factor (or superhydrostaticity factor) q λ that expresses the ratio of the effective density scale height to the hydrostatic scale height Thus, assuming an exponentially stratified atmosphere (Eq. 15), a density compression factor q n (t) ≤ q n,max in the CME front (Eq. 12), and adiabatic expansion inside the CME cavity (Eq. 14), we have the following time-dependent 3D density model: where d is the distance of an arbitrary location with 3D coordinates (x, y, z) to the instantaneous center position [x 0 (t), y 0 (t), z 0 (t)] of the CME, which is located at height h(t) vertically above the heliographic position (l, b).

STEREO/EUVI Observations
One CME event observed with STEREO that appears as a spherically expanding shell, and thus is most suitable for fitting with our analytical model, is the 2008-Mar-25, 18:30 UT, event. This CME occurred near the East limb for spacecraft STEREO/Ahead, and was observed on the frontside of the solar disk with spacecraft STEREO/Behind. Some 6 Aschwanden: CME Expansion and EUV Dimming preliminary analysis of this event regarding CME expansion and EUV dimming can be found in Aschwanden et al. (2009a), the CME mass was determined in white light with STEREO/COR-2 (Colaninno and Vourlidas 2009) and with STEREO/EUVI (Aschwanden et al. 2009b), and the 3D geometry was modeled with forward-fitting of various geometric models to the white-light observations (Thernisien et al. 2009;Temmer et al. 2009;Maloney et al. 2009;Mierla et al. 2009a,b). While most previous studies model the white-light emission of this CME, typically a few solar radii away from the Sun, our model applies directly to the CME source region in the lower corona, as observed in EUV. We follow the method outlined in Aschwanden et al. (2009a). Our model also quantifies the geometry and kinematics of the CME, as well as the EUV dimming associated with the launch of the CME. In order to determine the positional parameters of the CME as a function of time we trace the outer envelope of the CME bubble (by visual clicking of 3 points) in each image and each spacecraft and fit a circle through the 3 points in each image. The selected points for fitting the position of the CME bubble were generally chosen in the brightest features of the lateral CME flanks, but could not always been traced unambiguously in cases with multiple flank features. In those cases we traced the features that were closest to a continuously expanding solution. The radii and ypositions of the circular fits are fully constrained from the STEREO/A images, so that only the x-positions of the centroid of the spherical shell need to be fitted in the epipolar STEREO/B images. We note that the fits of the CME bubble roughly agree with the envelopes of the difference flux in the STEREO/B images initially (up to 18:48 UT), while there is a discrepancy later on. Apparently the CME has a more complex geometry than our spherical bubble model, which needs to be investigated further.

Fitting of Positional Parameters
This procedure yields the CME centroid positions [x A (t i ), y A (t i )] and [x B (t i ), y B (t i )] for the time sequence t i , i = 1, ..., 16. The images in Fig. 2 and 3 are displayed as a baseline difference (by subtracting a pre-CME image at 18:36 UT) to enhance the contrast of the CME edge. The circular fits to the CME outer boundaries are overlayed in Fig. 2 and 3. Both images have been coaligned and rotated into epipolar coordinates (Inhester et al. 2006), so that the y-coordinates of a corresponding feature are identical in the spacecraft A and B images, while the x-coordinates differ according to the spacecraft separation angle α sep , which amounts to α sep = 47.17 • at the time of the CME event. The epipolar coordinates measured from both spacecraft are then related to the heliographic longitude l, latitude b, and distance r c from Sun center as follows, which can directly be solved to obtain the spherical (epipolar) coordinates (l A , b, r c ), Therefore, using stereoscopic triangulation, we can directly determine the spherical coordinates (l i , b i , r c,i ), i = 1, ..., 16 for all 16 time frames, as well as obtain the CME curvature radii R(t i ) from the circular fits to the CMEs. We plot the so obtained observables l A (t), b(t), R(t), and h(t) = r c (t) − R ⊙ in Fig. 3 and determine our model parameters l and b from the averages. We obtain a longitude of l A = −102.4 • ± 0.9 • (for spacecraft STEREO/A), l B = l A + α sep = −54.9 • ± 0.9 • (for spacecraft STEREO/B), and a latitude b = −8.8 • ± 0.6 • . Thus, the CME source region is clearly occulted for STEREO/A. These epipolar coordinates can be rotated into an ecliptic coordinate system by the tilt angle β AB = 2.66 • of the spacecraft A/B plane. Viewed from Earth, the longitude is approximately l E ≈ −102.4 + α sep /2 ≈ −78.8 • . Thus, the CME source region is 12 • behind the East limb when seen from Earth. This explains why the EUV dimming is seen uncontaminted from post-flare loops, which are seen by STEREO/B but hidden for STEREO/A.

Fitting of Kinematic Parameters
We plot also the observables h(t) and R(t) in Fig. 4 and determine the model parameters h 0 , R 0 , a h , a R by fitting the quadratic functions R(t) (Eq. 1) and h(t) (Eq. 2), for which we obtain the starting time t 0 = 18.64 hrs (18:38 UT), the initial CME radius R 0 = 45 Mm, the initial height h 0 = 45 Mm, and the accelerations a R = 0.54 km s −2 for the CME radius expansion, and a h = 0.52 km s −2 for the height motion. The initial velocity is found to be negligible (v R ≈ 0 and v h ≈ 0). We estimate the accuracy of the acceleration values to be of order ≈ 10%, based on the uncertainty of defining the leading edge of the CME. Thus, we determined 9 out of the 14 free parameters of our model sofar.
Note that the acceleration measured here refers to the very origin of the CME in low coronal heights of < ∼ 0.6 solar radii observed in EUVI data. The acceleration is expected to be initially high and to rapidly decline further out, when the driving magnetic forces decrease at large altitudes. This explains why our values for the acceleration in low coronal heights are significantly higher than measured further out in the heliosphere, typically in the order of tens of m s −1 in height ranges of 5-22 solar radii,  as measured with SOHO/LASCO. SOHO/LASCO reported even a slightly negative acceleration at altitutes of 5-22 solar radii. The driving magnetic forces that accelerate a CME are clearly confined to much lower altitudes.

Fitting of Geometric Parameters
We model the 3D geometry of the CME bubble with the time-dependent radius R(t) and the width w 0 of the CME compression region. In Fig. 5 we show cross-sectional EUV brightness profiles across the CME in horizontal direction (parallel to the solar surface) and in vertical direction for the EUVI/A 171Å observations (indicated with dotted lines in Fig. 2). These baseline-subtracted profiles clearly show a progressive dimming with a propagating bright rim at the CME boundary, which corresponds to the density compression region at the lateral expansion fronts of the CME. The bright rims are clearly visible in the images during 18:46−18:56 UT shown in Fig. 2. The average width of the observed bright rims is w 0 ≈ 10 Mm, a value we adopt in our model.

Fitting of Physical Parameters
Finally we are left with the four physical parameters T 0 , q λ , n 0 , and q n,max . Since we show here only data obtained with the 171Å filter, the mean temperature is constrained by the peak temperature of the instrumental EUVI response function, which is at T 0 = 0.96 MK. This constrains the thermal scale height to λ T = 47, 000 × 0.96 = 45, 000 km. The remaining free parameters q λ , n 0 , and q n,max need to be determined from best-fit solutions by forward-fitting of simulated EUV brightness images (or profiles, as shown in Fig. 5) to observed EUV brightness images (or profiles). The EUV emission measure in each pixel position (x, y) can be computed by line-of-sight integration along the z-axis in our 3D density cube n e (x, y, z) per pixel area dA for each time from which the intensity I 171 (x, y) in the model image in units of DN s −1 can be directly obtained by multiplying with the instrumental response function R 171 (T ) of the 171Å filter, where R 171 = 3632 × 10 −44 DN s −1 cm 3 MK −1 and the FWHM of the 171 filter is ∆T 171 = 1.25 − 0.51 = 0.74 MK.
In Fig. 6 we show best-fit solutions of horizontal and vertical brightness profiles. The absolute flux level is proportional to the coronal base density squared, which we obtain by minimizing the mean flux difference between simulated and observed flux profiles. We obtain a best-fit value of n 0 = 6.5 × 10 8 cm −3 . The super-hydrostaticity factor is most sensitive to the vertical flux profile (Fig. 5, right-hand side panels), for which we find a best-fit value of q λ = 1.45. Thus, the average density scale height in the CME region is slightly super-hydrostatic, as expected for dynamic processes. These values are typical for quiet-Sun corona and active region conditions (see Figs. 6 and 10 in Aschwanden and Acton 2001).
The last free parameter, the maximum density compression factor q n,max , affects mostly the brightness of the CME rims. Fitting the brightness excess at the CME rims at those times where bright rims are visible are consistent with a value of q n,max ≈ 2.

Comparison of Numerical Simulations with Observations
After we constrained all 14 free parameters (listed in Table 1) of our analytical 4D model by fitting some observables, such as measured coordinates (Fig. 4) and cross-sectional horizontal and vertical brightness profiles (Fig. 5), we are now in the position to inter-compare the numerically simulated images with the observed images, as shown in Fig. 6 for 5 selected times, for both the STEREO/A and B spacecraft. The comparison exhibits a good match for the extent of the dimming region the the bright lateral rims, both extending over about 1.5 thermal scale heights above the solar surface. The basedifference images of EUVI/A reveal a fairly symmetric CME (as the model is by design), surrounded by spherical bright rims at the northern and southern CME boundaries (as the model is able to reproduce it). The model, however, is less in agreement with the observed EUVI/B images. The extent of the EUV dimming region matches relatively exactly, although the observed dimming region is somewhat cluttered with bright postflare loops that appear in the aftermath of the CME, which are mostly hidden in the EUVI/A observations. The biggest discrepancy between the model and the EUVI/B observations is the location of the brightest rim of the CME boundary. The combination of projection effects and gravitational stratification predict a brighter rim on the west side, where we look through a longer and denser column depth tangentially to the CME bubble, which is not apparent in the observations of EUVI/B. Instead, there is more bright emission just above the eastern limb that cannot be reproduced by the model. Apparently there exists stronger density compression on the eastern side of the CME bubble than the model predicts. Another inconsistency is the bright loop seen in EUVI/B at 18:51 UT, which does not match the surface of the modeled CME sphere as constrained by EUVI/A. Apparently, there are substantial deviations from a spherically symmetric CME bubble model that are visible in EUVI/B but not in EUVI/A. Perhaps a flux rope model could better fit the observations than a spherical shell model. These discrepancies between the observations and our simple first-cut model provide specific constraints for a more complex model (with more free parameters) that includes inhomogeneities in the density distribution of the CME.

Estimate of the CME Mass
Our model allows us, in principle, to estimate the CME mass by integrating the density n e (x, y, z, t) over the entire CME sphere, which is of course growing with time, but expected to converge to a maximum value once the CME expands far out into the heliosphere. A simple lower limit can analytically be obtained by integrating the density in the cylindrical volume above the footpoint area, M CME = m p n e (x, y, z)dV CME ≈ m p πR 2 0 n 0 λ T q λ . (22) From our best-fit values R 0 = 45 Mm, q λ = 1.45, n0 = 6.5 × 10 8 cm −3 and the thermal scale height of λ T = 47 Mm, we obtain a lower limit of M CME ≥ 0.47 × 10 15 g. However, this CME appears to expand in a cone-like fashion in the lowest density scale height, so the total volume and mass is likely to be about a factor of ≈ 2 higher. Moreover, the mass detected in 171Å amounts only to about a third of the total CME mass (Aschwanden et al. 2009b), so a more realistic estimate of the total CME mass is about a factor 6 higher than our lower limit, i.e., M CME ≈ 3 × 10 15 g, which brings it into the ballpark of previous CME mass determinations of this particular event, i.e., m CME = 2.9 × 10 15 g from STEREO/COR-2 white-light observations (Colaninno and Vourlidas 2009), or m CME = (4.3 ± 1.4) × 10 15 g from STEREO/EUVI observations (Aschwanden et al. 2009b).