The density compression ratio of shock fronts associated with coronal mass ejections

We present a new method to extract the three-dimensional electron density profile and density compression ratio of shock fronts associated with Coronal Mass Ejections (CMEs) observed in white light coronagraph images. We demonstrate the method with two examples of fast halo CMEs ($\sim$ 2000 km s$^{-1}$) observed on 2011 March 7 and 2014 February 25. Our method uses the ellipsoid model to derive the three-dimensional (3D) geometry and kinematics of the fronts. The density profiles of the sheaths are modeled with double-Gaussian functions with four free parameters and the electrons are distributed within thin shells behind the front. The modeled densities are integrated along the lines of sight to be compared with the observed brightness in COR2-A, and a $\chi^2$ approach is used to obtain the optimal parameters for the Gaussian profiles. The upstream densities are obtained from both the inversion of the brightness in a pre-event image and an empirical model. Then the density ratio and Alfv\'{e}nic Mach number are derived. We find that the density compression peaks around the CME nose, and decreases at larger position angles. The behavior is consistent with a driven shock at the nose and a freely-propagating shock wave at the CME flanks. Interestingly, we find that the supercritical region extends over a large area of the shock and last longer (several tens of minutes) than past reports. It follows that CME shocks are capable of accelerating energetic particles in the corona over extended spatial and temporal scales and are likely responsible for the wide longitudinal distribution of these particles in the inner heliosphere. Our results also demonstrate the power of multi-viewpoint coronagraphic observations and forward modeling in remotely deriving key shock properties in an otherwise inaccessible regime.


INTRODUCTION
Shocks associated with Coronal Mass Ejections (CMEs) are one of the sources responsible for highly energetic particles, called Solar Energetic Particles (SEPs; e.g., see reviews by Reames 1999;Desai and Giacalone 2016). While particle acceleration by flares (see Reames 1999) is expected to occur in a limited volume (i.e., magnetic reconnection site), fast-mode coronal shocks are able to directly inject SEPs over a broad range in heliolongitudes (e.g., Cliver et al. 1995Cliver et al. , 2005. SEPs are an important component of space weather as the cause of radiation hazards to astronauts and satellites, orbital degradation of satellites, communication disruptions, and electrical blackouts (e.g., Balan et al. 2014, and references therein). To better understand how coronal shocks accelerate particles over a wide range of heliolongitudes and hence improve our ability to predict their intensity, duration and energy spectrum, it is important to understand the properties of coronal shocks associated with CMEs and their temporal and spatial relationship with the SEPs measured in interplanetary space (IP).
In IP, the association between CME-driven shocks and particles is known for a long time thanks to direct in-situ measurements of shocks and SEPs (e.g., Cane et al. 1988). This has been difficult in the corona because direct measurements are currently unavailable, although the situation will soon be remedied (McComas et al. 2016). Understanding the formation and evolution of shocks in the corona is crucial for discriminating between flare and CME origins, particularly for the highest energy (GeV) SEPs ejected when CMEs are at heliocentric distances of around 2-10 solar radii (R ⊙ ; Tylka et al. 2005).
Coronal shocks are observed via remote-sensing observations in Extreme Ultraviolet (EUV), white light, and radio imaging and via spectroscopy (Vourlidas and Bemporad 2012). They form large-scale spherical fronts seen as halos by white light observations, or as EUV waves propagating against the solar coronal base (e.g., Kwon and Vourlidas 2017, see also recent reviews by Patsourakos and Vourlidas (2012), Warmuth (2015), and Long et al. (2017) for debates on the nature of EUV waves). Occasionally, type II radio bursts over a large spectral range accompany these shocks, particularly in IP space (e.g., Gopalswamy et al. 2008), while metric Type II emission seems to originate from the flanks of shock waves close to the Sun (e.g., Démoulin et al. 2012).
There is a growing amount of evidence for the association of CME-driven shocks with SEPs in the corona (see recent papers, e.g., Rouillard et al. 2012Rouillard et al. , 2016Carley et al. 2013;Lario et al. , 2016Salas-Matamoros et al. 2016, and references therein) but it remains unclear whether such shock waves are capable of accelerating particles at the observed energies Mancuso 2010, 2011). An important shock parameter that could be accessible from coronagraphic observations is the ratio between the downstream and upstream electron densities or density compression ratio, X. According to the diffusive shock acceleration theory, the slope of the SEP intensity spectrum depends only on the density compression ratio (e.g. Eq. (11) in Desai and Giacalone 2016). Therefore, measurements of the compression ratio and its temporal evolution in the corona can go a long way in understanding SEP in situ observations in the inner heliosphere.
Because the coronagraph images provide only the projected density, some model or knowledge of the density distribution along the line of sight (LOS) is necessary to obtain an estimate of the true volumetric electron density across the shock. Ontiveros and Vourlidas (2009) were the first to extract estimates of the compression ratio for a number of shock fronts in the Large Angle and Spectrometric Coronagraph (LASCO; Brueckner et al. 1995) field of view using forward modeling techniques. However, those were single viewpoint observations and hence the shock reconstructions might be subject to considerable uncertainties.
Here, we return to this issue employing multi-viewpoint observations from the Solar Terrestrial Relations Observatory (STEREO ; Kaiser et al. 2008) Sun-Earth Connections Coronal and Heliospheric Investigation (SECCHI; Howard et al. 2008) COR2 coronagraph and more sophisticated forward models and analysis methods. We reconstruct the threedimensional (3D) electron density at the shock sheath using an ellipsoidal forward model and populating it with a thin shell of Gaussian distributed electron density. The upstream density is derived with two methods: using the Leblanc et al. (1998) model and via the inversion of a pre-event polarized brightness image. The method is applied to two CMEs, on 2011 March 7 and 2014 February 25, associated with intense SEP events (Park et al. 2013;Lario et al. 2016).
The paper is organized as follows. In Section 2, we present the model and method used to determine the electron density distributions at the shock sheaths. In Section 3, we present the resulting sheath (downstream) and upstream electron densities and density ratios and infer Alfvénic Mach numbers. Our summary and conclusion are given in Section 4.  (b)) on the image plane is shown with the dotted line. The diamond refers to the geometric center of the ellipsoid model projected onto the same plane. Two examples of radial cuts (relative to the geometric center), used to take the brightness profile, are shown by the two rectangles marked as Cut A and Cut B. Thick lines in the rectangles show the length of L used for the calculation of χ 2 (see Figure 4). (b) The excess brightness in panel (a) with the 3D shock front (green lines). The 3D geometry was modeled with the ellipsoid model in Kwon and Vourlidas (2017). (c) The geometric relation among the Sun, shell-like sheath, and LOS. A partial circle around the origin O is the solar disk. A shell-like sheath is represented in gray color. Arrows in black, blue and red are the LOS, the projected shock normal on the image plane, and the actual shock normal in 3D, respectively.

Data
The 3D reconstruction of these CMEs and shocks has been reported in Kwon and Vourlidas (2017) using the Graduated Cylindrical Shell (GCS; Thernisien et al. 2006) model and the ellipsoid model model , respectively. The fits were based on three viewpoint observations from STEREO-A, -B and SDO (Solar Dynamics Observatory; Pesnell et al. 2012)/SOHO (SOlar and Heliospheric Observatory; Domingo et al. 1995). For the density ratio analysis, we use the COR2-A images because of their higher signal-to-noise ratio. The COR2 white light images capture Thomson-scattered photospheric light by the electrons in the corona in a field of view (FOV) of 2.5-15 R ⊙ . We do not analyze any COR1 images for these events because CME-associated shock waves are generally too faint in these heights to allow extraction of profiles.

3D Model of Density Structure in Sheaths
Our density analysis method uses the 3D geometry of shocks. Figure 1 shows the geometric relation among the 2D image plane, the spherical shock front in 3D and the LOS. Panel (a) shows an excess brightness due to the CME and shock, for example, taken at 2011 March 7 20:24 (UT) by COR2-A. The halo front on the image plane (dotted line in panel (a)) is the projection of a bubble-shaped shock front (green lines in panel (b)). The 3D geometry of the bubble-shaped shock fronts is modeled with the ellipsoid model . The sheath structure is given as a shell behind the ellipsoid model. Panel (c) shows a cross-section of the shell-like sheath. The image plane of observations (panel (a)) is defined as the Y -Z plane, and the X-axis is towards the observer. The origin O is at Sun center. L with an arrow in blue color is the distance normal to the projected shock front on the image plane. L is measured from the edge of the shock front on the image plane (dashed LOS arrow in this panel). Note that the shock normal on the image plane differs from the actual shock normal in 3D. The shock normal in 3D is shown by arrows in red, and τ denotes the distance normal to the shock front in 3D from the edge of the shock where the electron density due to the shock begins increasing. The excess brightness on the image plane results from the integration over the LOS passing through the shell-like sheath.
We would like to emphasize that the electron density jump along the 3D shock normal τ cannot be determined directly from the brightness profile along the projected shock normal L. As shown in Figure 1(c), a LOS passes through various τ , and the observed excess brightness is the integration over the LOS. A way to overcome this is to model the excess brightness I ′ (L) integrating the sheath electron density (N e (τ )) along the LOS (s) and compare it   . Because the Gaussian tail is not zero even at large distances, Ne(τ < 0) is not exactly zero, if d front = 0. By definition, we simply set Ne(τ < 0) = 0. The 3D geometry of the shock used to calculate I ′ is shown in Figure 1 with the observed excess brightness I(L). Once we obtain the best I ′ (L) for I(L), N e (τ ) along the shock normal is determined. To model N e (τ ), we use Gaussian function that is generally used for the density structure of a wave, but allow asymmetry by using two Gaussian functions, namely, where ρ e and d are constants. τ ′ is the full width at half maximum of the Gaussian function with d front . The outermost edge of the density structure is defined at L = 0, so that N e (τ < 0) = 0. ρ e is the peak electron density excess (density jump). The upper panels in Figure 2 show the resulting sheath electron density distributions N e (τ ) from various sets ρ e , d front , and d sheath . For instance, panel (a) shows N e (τ ) obtained when varying ρ e while the other two parameters are fixed. In general, N e (τ ) starts increasing at τ = 0, reaches its maximum at τ = τ ′ , and decreases gradually where τ > τ ′ . The key features of the function and the physical interpretations can be summarized as follows.
1. If d front = d sheath , it becomes a linear wave.
2. If d sheath is very large (dash-dot-dotted line in panel (b)) and d front = 0 (solid line in panel (c)), N e (τ ) becomes a step function near the shock front, followed by a gradually decreasing tail. This profile can represent steepening waves and shock waves.
3. if d front = 0 and d front < d sheath (panel (a)), it resembles the in-situ measurements and the theoretical expectations of collisionless shocks (see a review by Treumann 2009). In this case, d front serves as the density in the foot and ramp, with d sheath being the density in the sheath including overshoot and turbulence (if any).
The real density distribution may deviate from our double-Gaussian model. Note, however, that the total number of electrons will largely depend on the density excess ρ e , rather than the detailed shape of the function or fine scale structures (e.g., foot and overshoot). Since we fit the density model to observed brightnesses that result from the LOS integrations, we expect that the details of the profile is negligible when determining the density jump at shocks (see Figure 5). Note that Equation (2) is similar to the function in Thernisien et al. (2006). Thernisien et al. (2006) used the function to reproduce the enhanced brightness due to the pileup plasma of an expanding flux rope, whereas the function is used for the density structure of waves and shocks in this paper. Since observations do not resolve density distributions along the shock normal as discussed above, we calculate the Thomson-scattering brightness corresponding to the given density distribution N e (τ ) to be compared with the actual observation. We use the ellipsoid fits in Kwon and Vourlidas (2017) as the 3D geometry of the shock fronts. Given the 3D geometry as illustrated in Figure 1(c), the 3D coordinates of all points along the LOS are known, and thus τ along the LOS is also known.
Once the sheath electron distribution N e (τ ) and the coordinates of the LOS are given, Thomson-scattering brightnesses resulting from the integration can be calculated using the standard Thomson-scattering equations (e.g. Billings 1966). We assume that the sheath electron density structure varies slowly along the shock surface so that a single density distribution model can represent the density structure for the region where the LOS are passing through (see Figure 1(c)). This assumption is valid because the fits of the model to observations are done only close to the shock front (see Section 3.1).
The lower panels of Figure 2 show the calculated Thomson-scattering brightnesses I ′ (L) from electron density distributions N e (τ ) given in the upper panels, together with the 3D geometry of the shock shown in Figure 1(b). Note that the shapes of I ′ (L) differ from the sheath density distributions N e (τ ). For instance, I ′ (L = 0) at the shock front is zero, although the electron density has the maximum at the shock front (see the solid line in Figure 2(c)). These differences are due to the integration along the LOS.
As shown in Figure 2, the shape of I ′ varies with the density model N e (τ ) for the given geometries. In this sense, we find the optimal parameters of N e (τ ) by minimizing χ 2 defined as where I and I ′ are the observed and modeled excess brightnesses, respectively. m is the total number of the points compared, and w is the weight function. ∆L is the displacement on the sky plane of the true shock front from the initially estimated one, i.e., L = L ′ -∆L. Because of this, our fit has the four parameters, ρ e , d front , d sheath , and ∆L. Note that the difference between I and I ′ is normalized by I ′ . Instead of the error in I, we use a weight function w for the denominator, and the weight function is defined as,

Upstream Density Profile
We estimate the upstream electron density ρ u via inversion of a pre-event polarized brightness image. The method assumes an axisymmetric density distribution that can be described with an nth-order polynomial function (van de Hulst 1950;Hayes et al. 2001). We use the polynomial form in Leblanc et al. (1998), where A, B, and C are constants. The constants are found by the fit to the observed profiles at each position angle. Equation (4) is essentially equivalent to the standard polynomial expansion used traditionally in the field. Note that it reduces the number of constants while it includes the higher-order terms. Because of the presence of F-corona polarized component, especially above 5 Rs (Hayes et al. 2001), the determined ρ u in our height range would be overestimated. In contrast, the measured ρ u in coronal holes could be suppressed because of the calibration (e.g., Hayes et al. 2001;Thompson et al. 2010). We have checked other empirical coronal electron density models that have been used widely, i.e., Newkirk (1961), Saito et al. (1977), and Leblanc et al. (1998), to infer ρ u . Since the Leblanc et al. (1998) model provides the lowest electron density among them, we use as a lower limit. In this way, the upstream electron density ρ u is obtained as a range between our measurement and the Leblanc et al. (1998) model.

Error in Electron Density Excess
The main sources of error in ρ e are the uncertainties of brightness given by white light observations, the LOS integration relying on the ellipsoid model, the Gaussian description of the density distribution, and the relation between the white light structure and the true shock.
First of all, we use the excess brightness to determine ρ e . An excess brightness image is obtained by subtracting a previous image, which is temporally closest but at least 30 minutes apart from the image. In this manner, it is expected that the noise level, including the calibration error, the F-corona polarized component, the stray light component, and the brightness due to the background electron density, is similar in the two images, since they are temporally close. In this sense, the error in the excess brightness would be less than the error (20%; Frazin et al. 2012) of the direct brightness images. In addition, the CME and shock structure in the previous image cannot affect the excess brightness close to the shock front, because the two images are at least 30 minutes apart. Note that the shocks/waves in white light observations are generally faint and diffuse (e.g., Sheeley et al. 2000), so that the signal-to-nose level is still a issue. We use the COR2-A images because of their higher signal-to-noise ratio than the COR2-B images.
The 3D geometry is the key to deriving electron densities from white light observations. We rely on a highly idealized ellipsoid description for the 3D geometry of the shock front, and thus the LOS integration. Obviously the real shock front may deviate from that. Multi-spacecraft in-situ measurements at 1 au have shown that there exist the shock surface ripples so that the shock surface cannot be represented by a planar or spherical structure (Neugebauer and Giacalone 2005, and references therein). We believe that this effect is negligible in the corona. Instead, we determine the error in the ellipsoid model as shown in the Appendix in Kwon et al. (2014). This error gives us the difference of the ellipsoid model from the observed front. Because of the error (± δ) in the 3D geometry, we have the three ellipsoid models (see dashed lines in Figure 3), and the lower and upper estimates of the 3D geometry are used to estimate the error in ρ e .
The sheath electron density distribution model in Equation (2) is also highly idealized, and the real density distribution may deviate from our model. However, the brightness is obtained by taking the integration of the density over the LOS. It is obvious that the small scale structure will not affect much the total number of electrons and, therefore, the observed brightness. As we will see in Section in 3.1, the most significant source causing the brightness increase is the peak electron density excess ρ e (density jump).
The imperfection of the background subtractions will result in error in excess brightness. If we assume that the true ρ e varies slowly along the shock surface, the error can be estimated by comparing ρ e with those in the neighboring position angles. We repeat the same analysis for every position angles and take the average and standard deviation over 7 • . The standard deviation serves as the error. Note that we discard the cases where the brightness profiles are contaminated by the following pileup plasma or deflected streamers.
As discussed above, the two errors in ρ e are given from the uncertainties of the ellipsoid model and the excess brightness. We use the larger one.

RESULTS AND DISCUSSION
We apply our method to two fast CMEs on 2011 March 7 and 2014 February 25. The details of the 3D reconstructions with the ellipsoid model are given in Kwon and Vourlidas (2017) but we summarize the modeling results here for completeness. The shock speeds in 3D are ∼ 2200 km s −1 and ∼ 2050 km s −1 for the March and February events, respectively. The minimum angular widths of the shock envelopes are 192 • and 252 • , whereas the CME angular widths are 58 • and 90 • , respectively. The shock speeds vary with position angle. The maximum speeds are seen near the CME noses, in which the speeds are well correlated with the CME nose speeds. The speeds in the far-flanks tend towards the local fast magnetosonic speed.
This 3D modeling considered only the outline of the CME and its outer shock envelope. Here, we take the next step and attempt to fit the observed brightness of the shock by modeling the shock envelope as a thin shell with some electron density distribution. Thus we can determine the density excess ρ e at the sheaths and the downstreamupstream density ratios X, using the 3D geometry of the shocks. The 3D geometries and kinematics and the estimated electron densities enable us to also estimate the Mach numbers M A and upstream Alfvén speeds v A .
3.1. Density Excess ρ e at Shock Front Figure 3 shows the excess brightness COR2-A images for two events on 2011 March 7 (upper panels) and 2014 February 25 (lower panels). The previous images (at least 30 minutes apart) are subtracted from these images, so the brightness is largely due to the electrons in the CME and sheath. To derive the density structure of the sheath we take radial cuts (relative to the geometric center of the shock front; see the diamond symbols in the first panel of Figure 3) along all position angles at 1 • intervals. The position angle, ζ, is measured counterclockwise from the shock leading edge (directional axis of the ellipsoid model projected onto the image plane). The two rectangles in the first panel of Figure 3 (see also Figure 1(a) zoomed in on the shock front) are two example radial cuts, far from and close to the CME nose (Cut A and Cut B, respectively). To increase the signal-to-nose ratio, we average the brightness profiles across the width of the rectangle. The width ω is chosen as, ω = 2e r sin θ, where e r is the average distance of the projected shock front from its geometric center, and θ = cos −1 (1 − λ/e r ). λ is the distance corresponding to 0.5 pixel. In this way, the difference in distance L of the curved shock front across the width is at most 0.5 pixel, and thus the curvature is negligible when taking the average over the width.
The upper and lower panels of Figure 4 show the excess brightness profiles of Cut A and Cut B, respectively. Since the background brightness has been subtracted, the profiles are flattened, and the brightness in the region where L(= L ′ − ∆L) > 0 is 0, as indicated by the horizontal dotted line in each panel. While our density model in Equation (2 is only for the excess brightness due to the shock, the brightness profiles also contain the brightnesses owing to the plasma pileup ahead of the following CME and deflected streamers. In order to discriminate the shock part from the following pileup plasma in the excess brightness profiles, we use the slope of the profiles. The brightness at the outermost part (L ∼ 0) is expected to be mostly due to the density jump, but the following pileup plasma enhances the brightness. The top panel in Figure 4 shows the case that the brightness profile is taken at a position angle far away from the CME. While the brightness due to the pure sheath (red line; see Figure 5 and the description) is expected to decrease at large distances (|L|), the observed brightness increases with distance. Since the CME is far away from the shock front at this position angle, the effect of the pileup plasma on the brightness profile would be gentle. As indicated by the three blue lines in this panel, it results in the slower increase in brightness in the intermediate part (second blue line; -2 ≤ L ≤ -0.6) than the outermost part (third blue line; 0 ≤ L ≤ -0.6). Since the innermost part (first blue line) is mostly due to the pileup plasma, the gentle slope drops off into a steep slant as the distance increases. When the CME is close to the shock front, as seen in the lower panel of Figure  4, the slope keeps increasing with distance. The three-part slope enables us to discriminate the part due to the shock and minimizes the effect of the following pileup plasma on the fits. We only use the first increase part as demarcated by the two dashed lines in each panel (see also the thick lines in the rectangles in Figure 3). The red dots in Figure  3 refer to the position angles where we have been able to determine ρ e using χ 2 minimizations. In the position angles where the following CME is too close, we simply discard them. Figure 5 describes our χ 2 approach to obtain the optimal parameters of the double-Gaussian density distribution model from the observed brightness profiles I in Figure 4. The best model brightness profiles I ′ are shown as the red lines in Figure 4. We perturb the parameters (ρ e , d front , d sheath , ∆L), obtain the brightness I ′ (L), and then calculate χ 2 comparing I ′ (L) with the observed brightness I(L), as discussed in Sec. 2.2. The various density models are obtained by the various sets of parameters which are given with the intervals of ∆ρ e = 0.5N Leblanc (r), ∆d front = 0.1 R ⊙ , ∆d sheath = 0.5 R ⊙ and ∆∆L = 0.02 R ⊙ . The full ranges of the input parameters are shown in Figure 5 Leblanc et al. (1998) and Saito et al. (1977), respectively.
in the corona varies in several orders of magnitude. Once the minimum χ 2 is obtained, we calculate χ 2 again with the different intervals of the parameters. The full ranges of the parameters used for the second calculation are given as ± 2∆ ′ , where ∆ ′ is the initial intervals, centered on the parameters at the minimum χ 2 . The intervals of the parameters are ∆ρ e = 0.1N Leblanc (r), ∆d front = 0.02 R ⊙ , ∆d sheath = 0.1 R ⊙ , and ∆∆L = 0.02 R ⊙ . Not only the minimum χ 2 values, we check the resultant profiles visually as shown in Figure 4. The dependency of χ 2 on the four parameters ρ e , d sheath , d front , and ∆L is shown in Figure 5. The top and bottom panels are χ 2 values for Cut A and Cut B, respectively. Since χ 2 is normalized by the model brightness I ′ , the minimum χ 2 value is less than unity (see color bars). As shown in Figure 2(a), the value and the shape of I ′ are sensitive to ρ e . Because of this, as shown in Figures 5(a)-(c) and 5(e)-(g), χ 2 converges at a ρ e . It is also shown that d front = 0 (panels (b), (d), (f) and (h)), being consistent with the density profile of shock waves. It may imply that the shock thickness is very small and thus cannot be resolved in our observation and method. It is evident that d sheath (panels (a) and (e)) and ∆L (panels (c) and (g)) do not affect the measure of ρ e . The low χ 2 values (black colors) are aligned along the best ρ e value (see the vertical dashed lines in these panels). It demonstrates that our χ 2 approach can provide the reasonable measure of ρ e .
We repeat the χ 2 fitting at all position angles. Since it is applied also in the time series, the temporal evolution of the sheath structures can be investigated. Figure 6 Leblanc et al. (1998) and Saito et al. (1977), respectively.
background density in this figure has been determined independently in Section 3.2. As the shock propagates outward, ρ e falls off. Figure 7 shows the derived ρ e for the 2011 March 7 and 2014 February 25 CMEs. The upper panels show ρ e versus ζ for three consecutive times (marked by the different colors). Because of the shock-streamer interaction at 20:24 UT on 2011 March 7 (see the arrow in Figure 3), we do not attempt to derive ρ e beyond ζ ∼ 50 • . As expected, ρ e decreases with time, and tends to maximum around the shock leading edge (ζ ∼ 0 • ). In the lower panels, we plot ρ e versus heliocentric distance for all position angles. ρ e decreases with height but remains high at the leading edge (marked by the downward arrows). This is a clear indication, at least to us, that the super-magnetosonic CME keeps driving shocks and it results in the higher density jump than the lateral parts that are likely decoupled from the driver and turning into a freely-propagating shock wave (Kwon and Vourlidas 2017). Kwon and Vourlidas (2017) showed that the speed of the shocks ahead of the CME noses is well correlated with the CME nose speeds while the shock speed in the far-flanks is not. The far-flank speeds tend rather towards the local fast magnetosonic speed. For comparison, the empirical density models of Saito et al. (1977) and Leblanc et al. (1998) are shown by the dashed and solid lines, respectively.
The errors in ρ e are shown by error bars in Figure 7. To determine the errors, we take into account the two sources of the error in ρ e as described in Section 2.4. First, the error can arise from the uncertainty of the excess brightness. Assuming that ρ e varies slowly along the shock surface in the real corona, the fluctuations in ρ e determined within a small range of ζ could be due to the uncertainty. We have repeated our fits for every 1 • intervals of ζ independently, and we take the average and standard deviation over ζ spanning 7 • . The other source of error is the 3D geometry of the shock modeled with the ellipsoid. The error of the ellipsoid model has been determined as shown in the Appendix of Kwon et al. (2014). The error is ∼ 4% for all images we have analyzed. The two additional ellipsoid models considering ± 4% error are obtained, and their projections on the images are shown as dashed lines in Figure 3. We repeat the same χ 2 calculation with these ellipsoid models, and the results are considered as the error in ρ e due to the 3D geometry of the shock. Once the two errors are calculated, we use the larger one for the error in ρ e .

Density Compression Ratio
The density compression ratio, X = ρ d /ρ u , where ρ d = ρ e + ρ u , can only be determined if ρ u is known. We determine ρ u using the polarized brightness images prior to each event as described in Section 2.3. In practice, we have used the polarized brightness images observed from 00:08 UT to 14:08 UT on 2011 March 7 and from 19:08 UT on 2014 February 24 to 00:08 UT on 25. The standard polarized brightness inversion method is applied to the averaged polarized brightness images. Figure 8 shows the determined ρ u along the shock fronts where ρ e has been determined. Given the inhomogeneity of the background corona, the derived ρ u varies considerably along the shock fronts resulting in the 'U-shaped' curves. The errors in ρ u are determined from the errors in r due to the uncertainty of the 3D geometry of the shock. The errors in r and ρ u are shown by the horizontal and vertical error bars, respectively. Because of the uncertainty in the background electron density (e.g., Hayes et al. 2001), we also employ an empirical coronal electron density model of Leblanc et al. (1998) that provides the lowest values among the other empirical models (solid lines in Figure 8). In general, measured ρ u is greater than the Leblanc model probably due to F-corona polarized component, but a portion of the 2011 March 7 shock propagates in the northern polar coronal hole in the sky planes, and the measured ρ u in this region is lower than that of the Leblanc model due to the calibration (e.g., Thompson et al. 2010). The measure of electron density in coronal holes is even more uncertain (e.g., Hayes et al. 2001). We also consider the error in r when taking ρ u from the empirical model. We use the full ranges covering the observations and empirical model as well as their errors, as the upstream electron densities ρ u . Figure 9 shows the ratios between the downstream and upstream electron densities, i.e., compression ratios X. The speeds are also shown in colored lines for the fronts where X is determined. The error bars are the propagation errors given by the errors in ρ e and ρ u . The compression ratio peaks around the shock leading edge, as also seen in ρ e in Figure 7. Note that there is no dramatic increase in X around the shock leading edge of the 2014 February 25 shock as the 2011 March 7 event. It may be because of the upstream electron density ρ u and the proximity of the shock front to the CME. As seen in Figure 8, ρ u at the leading edge of the February event is lower than that of the March event. The 3D angular distances between the CME noses and the the projected shock leading edges are ∼ 40 • (2011 March 7) and ∼ 50 • (2014 February 25). In addition, X seems to be correlated with speed. However, it is interesting to note that X is nearly constant in time (cf. Bemporad and Mancuso 2011), while the speed decreases with time, as clearly seen in the bottom panels. Figure 6 indicates that this is due to the decrease in ρ u with height (see also Figure 8). As the shocks propagate outward, ρ e is falling off together with ρ u . Table 1 shows the full ranges of the determined X with the averages.
The ratios X in Figure 9 are found to be larger than 4 around the shock leading edge. The upper limit of X is given as (γ + 1)/(γ − 1), where γ is the adiabatic index, and the upper limit is 4 if γ = 5/3. It might imply that the adiabatic index γ could be close to 1 in the low corona (Van Doorsselaere et al. 2011). If γ = 4/3, for example, the upper limit is 7. Alternatively, the regions where X > 4 are only the small parts of the shock fronts as shown in Figure 9. The averages and standard deviations shown in Table 1 are lower than 4. It may imply that it is due to the errors in ρ e and/or ρ u , and the actual X is less then 4. As discussed in Section 3.1, the brightness profiles taken close to the shock leading edge could be slightly contaminated by the pileup plasma of the following CME.
Several attempts have been made in the past literature to determine the density ratios X from white light (Ontiveros and Vourlidas 2009;Mancuso 2010, 2011) and EUV (Kozarev et al. 2011;Kouloumvakos et al. 2014;Long et al. 2015), and radio observations (Ma et al. 2011). Ontiveros and Vourlidas (2009) determined X, ranging between 1.2 and 2.8 for 11 CMEs that were faster than 1500 km s −1 observed by LASCO C2 coronagraph. Mancuso (2010, 2011) also used the LASCO C2 images and found that X ≈ 1.2-3.0. Although these values were based on projected brightness profiles or single-view reconstructions, they are within the range of our values determined via more sophisticated 3D reconstructions. Similarly, for the variation along the shock front, Bemporad and Mancuso (2011) concluded that the compression ratio varied from 1.2 at the lateral flank to 3.0 at the shock leading edge. Our results validate past efforts and show that the compression ratio, in the middle corona at least, generally ranges from 2 to 4 (see Table 1). In addition, compression ratios were also determined in the early stage of events from EUV (Kozarev et al. 2011;Kouloumvakos et al. 2014;Long et al. 2015) and radio (Ma et al. 2011;Kouloumvakos et al. 2014) observations. These compression ratios are slightly larger than unity and tend to be smaller than those derived from white light observations. Note that these measurements have been done close to the Sun, indicating the early stage of the events. It may imply that the shock is still developing in those heights.
While the derived X values are consistent with past results, their temporal evolution differs in that the values are more or less constant in contrast to the results in Mancuso (2010, 2011). Bemporad  (2011) showed that the peak X value declines from ∼ 3.0 to ∼ 1.5 as the shock leading edge travels from ∼ 3 R ⊙ to ∼ 6 R ⊙ . Note that our CMEs are very fast, and the speeds at the times of the latest images are still over 1500 km s −1 . Kwon and Vourlidas (2017) showed that the speeds of these driver CMEs are much faster than the local fast-mode wave speeds and are able to generate shocks during our measurements. The fast speeds may be responsible for the high X even in higher altitudes.

Alfvénic Mach Number, M A
To relate our results to modeling effort and past literature, we proceed to compute the Alfvénic Mach number, M A , under the assumption of γ = 5/3. M A can be estimated from the compression ratio, X (e.g., Bemporad and Mancuso 2011, and references therein). M A is a function of plasma β and θ Bn , and we could simplify to the case of β → 0 for the coronal medium (Kwon et al. 2013b). θ Bn is the angle between the magnetic fields and the shock. We assume that the coronal magnetic fields are purely radial in our height range and determine θ Bn using the ray-like trajectories of the shocks shown in Kwon and Vourlidas (2017). The Mach number of an oblique shock is given by (e.g., Bemporad and Mancuso 2011, and references therein) where M A⊥ = [0.5X(X + 5)/(4 − X)] 1/2 and M A = X 1/2 , assuming the adiabatic γ = 5/3.  Figure 10. Dashed and dash-dotted lines refer to vA that are modeled with the empirical electron density models of Leblanc et al. (1998) and Saito et al. (1977), respectively, together with the empirical magnetic field function of B = 2.2(r/R⊙) −1 . Lower panel: Magnetic field strength B inferred from vA in the upper panel together with the determined ρu. Lines in dashed and dash-dotted are radial magnetic field profiles given in Mann et al. (1999) above the quiet Sun and Dulk and McLean (1978) above an active region, respectively. Figure 10 shows the derived M A from the full ranges of X in Figure 9. Dashed lines are the estimated θ Bn . Because of the assumption of purely radial magnetic fields, θ Bn decreases monotonically with time. The full ranges of M A and the averages are given in Table 1. Similarly to X, M A seems to be more or less constant in time, which is contrary to the prevailing wisdom (as shown in Bemporad and Mancuso 2011). This does not seem unreasonable since the local Alfvén speed should decrease with height. As a shock propagates outward the shock speed should also fall, and it results in a (relatively) constant M A since M A = V sh /v A (Figure 11).
Since M A = V sh /v A , we can use the Mach number to estimate the Alfvénic speed across the shock front. To that end, we need to move to the shock frame by estimating the ambient solar wind speed, V SW . Then, the true shock speed V sh = V ′ sh -V SW ·cos θ Bn , where V ′ sh is the measured one. Since V SW measurements does not exist in these coronal heights, we turn to a formula given in Sheeley et al. (1997), where A = 3.6 m s −1 and r 1 = 2.1 R ⊙ . In this manner, v A can be estimated and is shown in the upper panel of Figure  11. The upstream v A decreases with height and time but our estimates are higher than those predicted by empirical models of density and magnetic field (see the dashed and dash-dotted lines). This may not be unreasonable since our two case studies are very fast CMEs with speeds higher than 2000 km s −1 . Note that the magnetic field strengths inferred from the determined v A and ρ u in the lower panel of Figure 11 are consistent with the empirical models given in Mann et al. (1999) and Dulk and McLean (1978). We need further investigations including slower CMEs to see the general characteristics of shocks associated with CMEs. In addition, plasma β can also be determined from the estimated upstream v A by a relation that β = (2/γ)(c S /v A ) 2 , if sound speed c S is known (e.g., Kwon et al. 2013b). We obtain that β = 0.06 ± 0.02 in our height range, assuming γ = 5/3 and c S = 200 km s −1 . It is consistent with the low β assumption for the corona widely used in this field.

Physical Implications for SEP Acceleration
One of the challenges in SEP research is to explain the wide longitudinal distribution of SEPs originating in a single flare and CME event (Desai and Giacalone 2016, and references therein). The large widths of SEP-associated CMEs are the obvious candidate. While the spatial and temporal relationships between the CME-driven shock wave and SEP injections seem to support this assertion (e.g., Rouillard et al. 2012Rouillard et al. , 2016Lario et al. , 2016, it remains unclear whether the shocks are capable of accelerating particles at these locations since the properties of the local plasma environment (i.e., seed particle populations, turbulence levels) in the coronal heights are largely unknown. However, we have presented a method that allows us to derive the electron density distribution (in 3D) in the shock and hence obtain the density compression ratio and consequently the Mach number, Alfvén speed and other parameters (under further assumptions). So we can, at least, deduce whether, and most importantly where, our shocks have the potential to accelerate particles by comparing our Mach numbers to the critical Mach number, M c , (e.g., Bemporad and Mancuso 2011).
Red lines in Figure 10 show the critical Mach number corresponding to the estimated θ Bn . The critical Mach number M c for a collisionless shock is taken from Figure 2 in Treumann (2009), for the case of β ≈ 0 (we have obtained β < 0.1). It seems that our shocks are supercritical (M A > M c ) over significant portion of their extent. We also see that the condition extends for longer periods than those reported by Bemporad and Mancuso (2011). Note that our case events are associated with longitudinally-wide SEP events (Park et al. 2013;Richardson et al. 2014;Lario et al. 2016). Our work here provides significant additional support to these previous works that the SEPs are produced at the CME shocks and that the wide extent of these shocks is the reason for the distribution of SEPs over a very wide range of heliospheric longitudes.

SUMMARY AND CONCLUSION
We present a new method that, using multi-viewpoint observations and forward modeling techniques, enables us to model the observed brightness of shock fronts in white light coronagraphic observations and extract the threedimensional electron density distribution across the fronts. We apply the method to two case studies; the CME events on 2011 March 7 and 2014 February 25. Both CMEs were fast (> 2000 km s −1 ) and wide (≥ 200 • ) halo events. The 3D reconstructions of the CMEs and their shock envelopes were reported in a separate study (Kwon and Vourlidas 2017). The shock envelope is based on an ellipsoidal fitting and its white light emission is assumed to arise from electrons distributed in a thin shell over the surface of the shock. The electron distribution is assumed to be a double-Gaussian described by three free parameters. By varying the parameters and integrating the resulting density profiles along the known (from the ellipsoid model) LOS, we can derive modeled brightness profiles to be compared against the observed ones. We locate the best set of free parameters via a χ 2 minimization approach. We use two methods to estimate the background density; the empirical Leblanc et al. (1998) model and the standard polarized brightness inversion technique using the pre-event polarized brightness images. In this way, we obtain an upper and lower estimate for the background density and thus we can bound the resulting density compression ratio. The final results of this exercise are the excess electron density ρ e , density ratio X and Alfvén Mach number M A across the full shock front. This is the first time that the properties of white light shocks are quantified in this way. Our results corroborate past estimates, mostly based on single viewpoint observations (Ontiveros and Vourlidas 2009;Mancuso 2010, 2011) and are fully consistent with the presence of radio emissions and SEPs in the two cases.
Our findings can be summarized as follows, 1. The density excess, ρ e peaks at and around the CME leading edge and the peak is maintained during the time interval we analyzed (Figure 7). This finding indicates that the two CMEs drive bow-type shock waves around their noses while the shock waves in the lateral flank propagate nearly freely, (i.e. as blast waves). This is the same conclusion we reached in the previous paper (Kwon and Vourlidas 2017) using the 3D speeds and a simple model. 2. The density compression ratio, X (Figure 9), and Alfvénic Mach number, M A (Figure 10), also vary with the position angle ζ, i.e., their maximum occurs around the CME nose and correlates with the shock speed. Both X and M A remain relatively constant in time and distance despite the decrease of the shock speed (see also Table  1). It is not unreasonable because the local Alfvén speed v A also decreases as the shocks propagate further out ( Figure 11).
3. The shocks are supercritical over a wider spatial range, and they last longer than those of what has been reported by Bemporad and Mancuso (2011).
4. The averages of X and M A are 2.1-2.6 and 1.1-1.8, respectively (see Table 1).
Once again, we show (this time via density analysis) that the diffuse white light emission ahead of fast CMEs outlines the shock sheath. Our 3D analysis suggests that the density compression is sufficiently high for the production of high energy particles and that the conditions last for, at least, tens of minutes over a wide range of position angles around the CME nose. It suggests that the CME-associated shocks can account both for the SEP production and their longitudinal distribution.
Our method builds upon and greatly extends past work (Thernisien et al. 2006;Ontiveros and Vourlidas 2009) by using the 3D reconstruction information of the shock envelope to deconvolve the density distribution from projection effects and to derive the 3D speed, upstream and downstream densities, Mach number and other parameters across the full shock front. The method holds promise for improving the extraction of quantitative information from shocks in the corona using remote sensing observations. It can readily use observations from different vantage points, including from the imagers (Howard et al. 2013;Vourlidas et al. 2016) aboard the upcoming Solar Orbiter (Müeller et al. 2013) and Solar Probe Plus (Fox et al. 2016) missions to be compared with direct in situ measurements from these missions. The technique will also greatly improve the results from any future joint coronagraphic imaging and off-limb spectroscopy of shocks (Bemporad and Mancuso 2011;Vourlidas and Bemporad 2012). We plan to test the double-Gaussian density profile along the shock normal, 3D geometry and the χ 2 approach with numerical simulations. Further analyses, including slower CME-shock events, will lead to a better understanding of the shock properties and its relation to SEPs and other phenomena such as EUV waves and radio bursts.
We would like to thank Rob Decker and David Lario for insightful discussions. This work of R.Y. K. and A.V. is supported by NASA Grant NNX16AG86G issued under the HSR Program. The SECCHI data are produced by an international consortium of the NRL, LMSAL and NASA GSFC (USA), RAL and Univ. Bham (UK), MPS (Germany), CSL (Belgium), IOTA and IAS (France). The editor thanks two anonymous referees for their assistance in evaluating this paper.