© R.-Y. Kwon and A. Vourlidas, Published by EDP Sciences 2018

1 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., 1995, 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 et al., 2013a, 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., 2012, 2016; Carley et al., 2013; Lario et al., 2014, 2016; Salas-Matamoros et al., 2016, and references therein) but it remains unclear whether such shock waves are capable of accelerating particles at the observed energies (Bemporad and 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 three-dimensional (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.

2 Method

2.1 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 (Kwon et al., 2014), 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.

2.2 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 (Kwon et al., 2014). 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 1c, 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 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, (1) 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.

• If d_front = d_sheath, it becomes a linear wave;

• 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;

• if  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 brightness that results from the LOS integrations, we expect that the details of the profile is negligible when determining the density jump at shocks (Fig. 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 1c, 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 brightness 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 (Fig. 1c). This assumption is valid because the fits of the model to observations are done only close to the shock front (Section 3.1).

The lower panels of Figure 2 show the calculated Thomson-scattering brightness I ′ (L) from electron density distributions N_e(τ) given in the upper panels, together with the 3D geometry of the shock shown in Figure 1b. 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 Fig. 2c). 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 χ² defined as (2) where I and I′ are the observed and modeled excess brightness, 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, (3)

Fig. 1 (a) An excess brightness COR2-A image taken at 2011 March 7 20:24 (UT). The projection of the 3D shock front (green lines in panel (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 (Fig. 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.

Fig. 2 The upper three panels show the dependency of the sheath electron density profiles (Ne(τ)) on the three parameters, ρe (panel (a)), dsheath (b), and dfront (c), while the other two parameters are fixed. The lower panels are the modeled Thomson-scattering brightness profiles I ′ (L) in mean solar brightness (B⊙), resulting from the given Ne(τ) in the same line-styles in the upper panels. Lines in solid, dotted, dashed, dash-dotted, dash-dot-dotted refers to the cases when ρe = (1.6, 3.1, 4.7, 6.2, 7.8) × 104, dsheath = 1 R⊙, and dfront = 0.25 R⊙ in panel (a); ρe = 6.2 × 104, dsheath = (0.6, 1.0, 1.4, 1.8, 8.0) R⊙, and dfront = 0.25 R⊙ in panel (b); and ρe = 6.2 × 104, dsheath = 1 R⊙, and (0.0, 0.1, 0.3, 0.4, 0.5) R⊙ in panel (c). Because the Gaussian tail is not zero even at large distances, Ne(τ < 0) is not exactly zero, if dfront ≠ 0. By definition, we simply set Ne(τ < 0) = 0. The 3D geometry of the shock used to calculate I′ is shown in Figure 1b.

2.3 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), (4) 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 R_⊙ (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.

2.4 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 min 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 min 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 an 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 Fig. 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 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 the 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.

Fig. 3 Excess brightness images used to determine the density excess ρe. The upper and lower panels show two different events on 2011 March 7 and 2014 February 25, respectively. Dots in red are the points where the ρe measurements have been made. Solid lines refer to the projected shock front determined with the ellipsoid model. Dashed lines represent the ±4% error in the geometry of the shock fronts.

3 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⁻¹ and ∼2050 km s⁻¹ 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 downstream–upstream density ratios X, using the 3D geometry of the shocks. The 3D geometries 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 min 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 Fig. 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 (Fig. 1a 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-noise ratio, we average the brightness profiles across the width of the rectangle. The width ω is chosen as, ω = 2e_rsinθ, where e_r is the average distance of the projected shock front from its geometric center, and . λ is the distance corresponding to 0.5 pixels. In this way, the difference in distance L of the curved shock front across the width is at most 0.5 pixels, 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 brightness 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 Fig. 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 minimize 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 Fig. 3). The red dots in Figure 3 refer to the position angles where we have been able to determine ρ_e using χ² minimizations. In the position angles where the following CME is too close, we simply discard them.

Figure 5 describes our χ² 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 , obtain the brightness I ′ (L), and then calculate χ² comparing I ′ (L) with the observed brightness I(L), as discussed in Section 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. Note that Δρ_e shown is normalized to the Leblanc model at the corresponding height r because the electron density in the corona varies in several orders of magnitude. Once the minimum χ² is obtained, we calculate χ² 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 χ². 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 χ² values, we check the resultant profiles visually as shown in Figure 4.

The dependency of χ² on the four parameters ρ_e, d_sheath, d_front, and ΔL is shown in Figure 5. The top and bottom panels are χ² values for Cut A and Cut B, respectively. Since χ² is normalized by the model brightness I′, the minimum χ² value is less than unity (see color bars). As shown in Figure 2a, the value and the shape of I′ are sensitive to ρ_e. Because of this, as shown in Figure 5a–c and e–g, χ² 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 χ² values (black colors) are aligned along the best ρ_e value (see the vertical dashed lines in these panels). It demonstrates that our χ² approach can provide the reasonable measure of ρ_e.

We repeat the χ² fitting at all position angles. Since it is also applied in the time series, the temporal evolution of the sheath structures can be investigated. Figure 6 shows the time evolution of the sheath propagating at the position angle ζ = 80° (Cut A). The position angle in three different images is shown in the upper panels of Figure 3. The 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 Fig. 3), we do not attempt to derive ρ_e beyond ζ∼ 50 °. As expected, ρ_e decreases with time and tends to have the 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 the 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 χ² 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.

Fig. 4 Excess brightness profiles I(L = L' − ΔL) taken from Cut A (upper panel) and Cut B (lower) shown in Figure 3. The dotted line in each panel indicates the zero brightness. The red line in each panel is the modeled excess brightness profile I ′ (L) obtained from χ2 fits shown in Figure 5. The blue lines indicate the slope of the observed brightness profiles. The two dashed lines in each panel demarcate the region where the χ2 values are calculated.

Fig. 5 Fits of the sheath density distribution model to the observed brightness profiles of Cut A (upper) and Cut B (lower). From the leftmost to the rightmost, the panels show the χ2 values varying with the parameters, dsheath − ρe, dfront − ρe, ΔL − ρe and dfront − dsheath, respectively, and the χ2 values are indicated by the color bars. Dashed lines in each panel indicate the parameters at the minimum χ2.

Fig. 6 Evolution of the sheath taken at Cut A (ζ = 80 °) of the 2011 March 7 event (rectangles in the top panels of Figure 3). From the top to bottom, each panel shows the modeled sheath electron density distribution plus the background electron density at 20:24 UT, 20:39 UT, and 20:54 UT, respectively. Note that the ranges of the abscissas are different.

Fig. 7 Determined electron density excess ρe for the 2011 March 7 and 2014 February 25 events. The results were obtained at three different times for each event, and they are shown in different colors, black, red, and cyan in chronological order (20:24 UT, 20:39 UT, and 20:54 UT for the 2011 March 7 event and 01:24 UT, 01:39 UT, and 01:54 UT for the 2014 February 25 event). The upper row shows the ρe distributions over position angle ζ. The position angle that ζ = 0° is at the projected directional axis of the ellipsoid model. The lower row shows the ρe distributions over height. Lines in solid and dashed represent two empirical electron density models in Leblanc et al. (1998) and Saito et al. (1977), respectively.

3.2 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 Fig. 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 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 (Fig. 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 than 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; Bemporad and 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⁻¹ observed by LASCO C2 coronagraph. Bemporad and 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 (Tab. 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 Bemporad and Mancuso (2010, 2011). Bemporad and Mancuso (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⁻¹. 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 at higher altitudes.

Fig. 8 Derived upstream electron densities ρu for the 2011 March 7 and 2014 February 25 events. The densities determined at three different frons are shown in black, red, and cyan, in accordance with Figure 7. Lines in solid and dashed represent two empirical electron density models in Leblanc et al. (1998) and Saito et al. (1977), respectively.

Fig. 9 Density ratios of downstream over upstream X = ρd/ρu. The left and right columns are the results for the 2011 March 7 and 2014 February 25 events, respectively. The three rows from the top show X over position angle ζ. The bottom panels show the determined X over height. Solid lines in black, red and cyan are the speeds of the shock fronts derived from the 3D ellipsoid model, in accordance with the ordinate on the right-handed side.

3.3 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) (5) where M_A⊥ = [0.5X(X + 5)/(4 − X)]^1/2 and M_A∥ = X^1/2, assuming the adiabatic γ = 5/3.

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 (Fig. 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 , where 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), (6) where A = 3.6 m s⁻¹ and r₁ = 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⁻¹. 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)², 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⁻¹. It is consistent with the low β assumption for the corona widely used in this field.

Fig. 10 Mach numbers MA (shaded regions) derived from the full ranges of X in Figure 9. Red lines indicate the critical Mach number for the determined θBn. θBn is also plotted by dashed lines (see the y-axes on the right-handed side).

Fig. 11 Upper panel: Local Alfvén speed, vA (plus symbols), inferred from the Mach numbers in 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.

3.4 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., 2012, 2016; Lario et al., 2014, 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 a 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.

4 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 three-dimensional 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⁻¹) 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 χ² 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; Bemporad and 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,

The density excess, ρ_e peaks at and around the CME leading edge and the peak is maintained during the time interval we analyzed (Fig. 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;

• The density compression ratio, X (Fig. 9), and Alfvénic Mach number, M_A (Fig. 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 (Tab. 1). It is not unreasonable because the local Alfvén speed v_A also decreases as the shocks propagate further out (Fig. 11);

• 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);

• the averages of X and M_A are 2.1–2.6 and 1.1–1.8, respectively (Tab. 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 χ² 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.

Acknowledgements

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.

References

All Tables

All Figures

Fig. 1 (a) An excess brightness COR2-A image taken at 2011 March 7 20:24 (UT). The projection of the 3D shock front (green lines in panel (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 (Fig. 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.

In the text

Fig. 2 The upper three panels show the dependency of the sheath electron density profiles (Ne(τ)) on the three parameters, ρe (panel (a)), dsheath (b), and dfront (c), while the other two parameters are fixed. The lower panels are the modeled Thomson-scattering brightness profiles I ′ (L) in mean solar brightness (B⊙), resulting from the given Ne(τ) in the same line-styles in the upper panels. Lines in solid, dotted, dashed, dash-dotted, dash-dot-dotted refers to the cases when ρe = (1.6, 3.1, 4.7, 6.2, 7.8) × 104, dsheath = 1 R⊙, and dfront = 0.25 R⊙ in panel (a); ρe = 6.2 × 104, dsheath = (0.6, 1.0, 1.4, 1.8, 8.0) R⊙, and dfront = 0.25 R⊙ in panel (b); and ρe = 6.2 × 104, dsheath = 1 R⊙, and (0.0, 0.1, 0.3, 0.4, 0.5) R⊙ in panel (c). Because the Gaussian tail is not zero even at large distances, Ne(τ < 0) is not exactly zero, if dfront ≠ 0. By definition, we simply set Ne(τ < 0) = 0. The 3D geometry of the shock used to calculate I′ is shown in Figure 1b.

In the text

Fig. 3 Excess brightness images used to determine the density excess ρe. The upper and lower panels show two different events on 2011 March 7 and 2014 February 25, respectively. Dots in red are the points where the ρe measurements have been made. Solid lines refer to the projected shock front determined with the ellipsoid model. Dashed lines represent the ±4% error in the geometry of the shock fronts.

In the text

Fig. 4 Excess brightness profiles I(L = L' − ΔL) taken from Cut A (upper panel) and Cut B (lower) shown in Figure 3. The dotted line in each panel indicates the zero brightness. The red line in each panel is the modeled excess brightness profile I ′ (L) obtained from χ2 fits shown in Figure 5. The blue lines indicate the slope of the observed brightness profiles. The two dashed lines in each panel demarcate the region where the χ2 values are calculated.

In the text

Fig. 5 Fits of the sheath density distribution model to the observed brightness profiles of Cut A (upper) and Cut B (lower). From the leftmost to the rightmost, the panels show the χ2 values varying with the parameters, dsheath − ρe, dfront − ρe, ΔL − ρe and dfront − dsheath, respectively, and the χ2 values are indicated by the color bars. Dashed lines in each panel indicate the parameters at the minimum χ2.

In the text

	Fig. 6 Evolution of the sheath taken at Cut A (ζ = 80 °) of the 2011 March 7 event (rectangles in the top panels of Figure 3). From the top to bottom, each panel shows the modeled sheath electron density distribution plus the background electron density at 20:24 UT, 20:39 UT, and 20:54 UT, respectively. Note that the ranges of the abscissas are different.
In the text

Fig. 7 Determined electron density excess ρe for the 2011 March 7 and 2014 February 25 events. The results were obtained at three different times for each event, and they are shown in different colors, black, red, and cyan in chronological order (20:24 UT, 20:39 UT, and 20:54 UT for the 2011 March 7 event and 01:24 UT, 01:39 UT, and 01:54 UT for the 2014 February 25 event). The upper row shows the ρe distributions over position angle ζ. The position angle that ζ = 0° is at the projected directional axis of the ellipsoid model. The lower row shows the ρe distributions over height. Lines in solid and dashed represent two empirical electron density models in Leblanc et al. (1998) and Saito et al. (1977), respectively.

In the text

	Fig. 8 Derived upstream electron densities ρu for the 2011 March 7 and 2014 February 25 events. The densities determined at three different frons are shown in black, red, and cyan, in accordance with Figure 7. Lines in solid and dashed represent two empirical electron density models in Leblanc et al. (1998) and Saito et al. (1977), respectively.
In the text

Fig. 9 Density ratios of downstream over upstream X = ρd/ρu. The left and right columns are the results for the 2011 March 7 and 2014 February 25 events, respectively. The three rows from the top show X over position angle ζ. The bottom panels show the determined X over height. Solid lines in black, red and cyan are the speeds of the shock fronts derived from the 3D ellipsoid model, in accordance with the ordinate on the right-handed side.

In the text

	Fig. 10 Mach numbers MA (shaded regions) derived from the full ranges of X in Figure 9. Red lines indicate the critical Mach number for the determined θBn. θBn is also plotted by dashed lines (see the y-axes on the right-handed side).
In the text

Fig. 11 Upper panel: Local Alfvén speed, vA (plus symbols), inferred from the Mach numbers in 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.

In the text