Quantifying the Magnetic Structure of a Coronal Shock Producing a Type II Radio Burst

Coronal shocks can be approximated as having a symmetric 3D bow-shock geometry (Ontiveros & Vourlidas 2009; Chen et al. 2014). In the cylindrical coordinate systems, the shape of the bow shock can be given by the following formula (Smith et al. 2003):

$\begin{eqnarray}&&Z=h-\displaystyle \frac{d}{s}\times {\left(\displaystyle \frac{\sqrt{{R}^{2}}}{d}\right)}^{s},\end{eqnarray} \tag{ 1 }$

where h describes the apex height of the shock, s determines the opening angle of the shock, and d is the semilatus rectum that controls the width of the shock. The Z-axis in the heliocentric coordinate system is assumed to be the initial symmetrical axis. Besides, we input two sets of angular parameters to fit the shock shape: the latitude and longitude (θ, ϕ) of the eruption location with the corresponding base vector e , and the latitude and longitude ($\theta ^{\prime} ,\phi ^{\prime}$) of the eruption direction with the corresponding base vector ${\boldsymbol{e}}^{\prime}$. According to (θ, ϕ) and ($\theta ^{\prime} ,\phi ^{\prime}$), the rotation matrices can be constructed to rotate the shock to the correct position. Note that there is a problem known as the gimbal lock, which is the loss of one degree of freedom in 3D space. This happens when two of the rotation axes are driven into a parallel configuration, and the system degenerates into rotating in 2D space, just like "locking."

In order to avoid the gimbal lock problem, we do not take the X-, Y-, and Z-axes in the coordinate system as the rotation axes. For example, when we rotate the fitting structure from the position (θ, ϕ) to the position ($\theta ^{\prime} ,\phi ^{\prime}$), the rotation operations are as follows: the rotation axis n is taken as ${\boldsymbol{e}}\times {\boldsymbol{e}}^{\prime}$, and the rotation angle α is taken as $\arccos ({\boldsymbol{e}}\bullet {\boldsymbol{e}}^{\prime} )$. According to n and α, the rotation matrix M_r is constructed as

$\begin{eqnarray}{M}_{r}=\left[\begin{array}{ccc}{n}_{x}{n}_{x}(1-\cos \alpha )+\cos \alpha & {n}_{x}{n}_{y}(1-\cos \alpha )+{n}_{z}\sin \alpha & {n}_{x}{n}_{z}(1-\cos \alpha )-{n}_{y}\sin \alpha \\ {n}_{x}{n}_{y}(1-\cos \alpha )-{n}_{z}\sin \alpha & {n}_{y}{n}_{y}(1-\cos \alpha )+\cos \alpha & {n}_{y}{n}_{z}(1-\cos \alpha )+{n}_{x}\sin \alpha \\ {n}_{x}{n}_{z}(1-\cos \alpha )+{n}_{y}\sin \alpha & {n}_{y}{n}_{z}(1-\cos \alpha )-{n}_{x}\sin \alpha & {n}_{z}{n}_{z}(1-\cos \alpha )+\cos \alpha \end{array}\right],\end{eqnarray} \tag{ 2 }$

where n _i is the projection of n on the X-, Y-, and Z-axes. Thus, it needs to be rotated once from (θ, ϕ) to ($\theta ^{\prime} ,\phi ^{\prime}$), and the gimbal lock problem can be avoided through the above operations. Then, we can adjust the rotation parameters (θ, ϕ) and ($\theta ^{\prime}$, $\phi ^{\prime}$) and the shock shape parameters (h, s, and d) until Equation (1) best matches the shock wave fronts in the EUV observations of SDO and STEREO. The wavelengths of 193 and 211 Å for AIA and 195 Å for the EUVI are suitable for observations of coronal shocks (Ma et al. 2011; Su et al. 2015). Therefore, the observations of these wavelengths are used to fit the shock surface. The fitted surfaces of the shock at 09:24:42 and 09:26:18 UT are represented by the blue isolines in Figure 4.

Before determining the magnetic structure around the radio source, the uncertainty of the 3D location of the radio source needs to be discussed. A conventional approach is to estimate the radial height of the radio source using the coronal density model (Zucca et al. 2014b; Su et al. 2015; Morosan et al. 2019, 2020). The propagation effect of electromagnetic (EM) waves in the plasma is considered in this work. Since the turbulent plasma can cause dispersion and scattering when EM waves propagate in plasmas (Su et al. 2021), the propagation effect may affect the apparent positions of the radio sources (Kontar et al. 2017; Chen et al. 2020). In this work, we use ray-tracing simulation results (Zhang et al. 2021) to estimate the uncertainty of the radio source location. The refraction and anisotropic scattering effect of the solar radio emission are considered in the calculation. With the longitude of the radio source θ = 65° (consistent with the eruption direction of the shock), the relative density fluctuation variance = 0.3, and the anisotropic parameter α = 0.4, we find that the location uncertainties of the radio source in the plane of the sky at the harmonic frequencies 173.2 and 150.9 MHz are about 0.024 and 0.023R_⊙, respectively.

The coordinates of the centroids of the NRH radio sources at 173.2 and 150.9 MHz are (−1.12, −0.04)R_⊙ and (−1.23, −0.17)R_⊙ in the plane of the sky, respectively. We take the beam size of the NRH observations as the location uncertainties. For the 173.2 MHz radio map, the location uncertainties are about ±0.09R_⊙ and ±0.22R_⊙ on the x-axis and y-axis in the plane of the sky (see the left panel of Figure 5), respectively; for the 150.9 MHz radio map, the location uncertainties are about ±0.11R_⊙ and ±0.25R_⊙ in the x-axis and y-axis in the plane of the sky (see the right panel of Figure 5), respectively. In the Sun–Earth direction, we cannot get the radio centroid locations from the NRH observations directly, considering the projection effect, so we assume that the longitude of the radio source is approximate to that of the corresponding active region (θ = 65°); thus, the radio centroids in the Sun–Earth direction are 0.47R_⊙ and 0.52R_⊙ at 173.2 and 150.9 MHz, respectively. We take the location uncertainties along the Sun–Earth direction to be the mean size of the beam in the plane of the sky, which are about ±0.17 and ±0.18R_⊙ at 173.2 and 150.9 MHz, respectively. In this way, the 3D uncertainty volumes of the radio source locations can be approximated as elliptic cylinders. The radio source centroids with location uncertainties in the plane of the sky are shown in Figure 6, and the shock surface in the 3D uncertainty volumes of the radio source at 173.2 and 150.9 MHz are shown as the red shadows in Figures 7(b), (c), (e), and (f). We can see that most parts of the shock surface are in the 3D uncertainty volume of the radio source.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The PFSS model is used to extrapolate the magnetic field around the shock front. Based on the extrapolated magnetic field, we can trace the magnetic field lines in the corona, and select the lines around the source region of the burst. The magnetic field lines around the shock surface are shown as the white and green lines in Figure 6, representing closed and open field lines, respectively. Thus, the 3D magnetic structure around the shock front is constructed. Since the PFSS model is a potential-field model, it is only suitable for describing the magnetic field in the region undisturbed by the shock (upstream), not the magnetic field in the region disturbed by the shock (downstream).

Next, in order to distinguish whether the shock front at the source region of the radio burst is quasi-perpendicular or quasi-parallel, we need to describe the magnetic structure around the shock surface quantitatively. We overlay the magnetic field lines around the radio source region onto the EUV images observed by SDO and STEREO. We can further obtain the intersection points between the fitted surface of the shock and the magnetic field lines. Then, we calculate the normal base vector ( e _n) of the fitted shock surface and the tangent base vector ( e _t ) of the magnetic field lines at the intersection points. The angle θ_Bn between e _n and e _t can be calculated as arccos( e _t • e _n ). Usually, if the quantity of θ_Bn is less than 45°, the shock front at the radio source is regarded as a quasi-parallel shock; if the quantity of θ_Bn is larger than 45°, the shock front is considered a quasi-perpendicular shock. We sample the points on the shock surface evenly, and get θ_Bn of these points on the shock surface.

The distributions of θ_Bn on the shock surface at 09:24:42 and 09:26:18 UT are shown in Figures 7(b) and (e), where the noses of the shock are indicated by the white asterisks. We also plot histograms of θ_Bn of the radio source at 09:24:42 and 09:26:18 UT in Figure 8. We find that the distributions of θ_Bn can be roughly separated into two parts, one for θ_Bn ≲ 45°, the other for θ_Bn ≳ 45°, meaning that θ_Bn has an obvious bimodal distribution of θ_Bn. The bimodal distribution of θ_Bn implies that both quasi-parallel and quasi-perpendicular shocks may be at work in the generation of the type II radio burst. The average values of θ_Bn of the radio sources on the shock surface are 36° ± 16° and 42° ± 20° at 09:24:42 and 09:26:18 UT, respectively.

Zoom In Zoom Out Reset image size

The magnetic field strength at the shock upstream (B₁) is here obtained by the PFSS model. The distributions of B₁ on the shock surface at 09:24:42 and 09:26:18 UT are shown in Figures 7(c) and (f), respectively. The mean values of B₁ of the radio source on the shock surface are 3.05 ± 0.74 and 1.70 ± 0.39 Gauss at 09:24:42 and 09:26:18 UT, respectively. The histograms of B₁ of the radio source at 09:24:42 and 09:26:18 UT are shown in the second row of Figure 8. It is seen that, different from the bimodal distributions of θ_Bn, the value of B₁ shows a roughly unimodal distribution, but with asymmetry.

The Alfvén Mach number, M_A = v/v_A , is defined as the ratio between the speed v of the shock and the Alfvén speed v_A . v_A is expressed as

$\begin{eqnarray}&&{v}_{A}=\displaystyle \frac{{B}_{1}}{\sqrt{{\mu }_{0}\rho }},\end{eqnarray} \tag{ 3 }$

where μ₀ is the magnetic permeability of the vacuum, and B₁ and ρ are the magnetic field strength and the mass density of the shock upstream (Priest 2014). Here, we expect to get the distributions of v_A and M_A on the shock front. As mentioned in Section 3.2, the distribution of B₁ on the shock surface can be derived from the PFSS model, and the results are shown in Figures 7(c) and (f).

Besides the distribution of B₁, we also need to obtain the distribution of ρ in order to estimate the distribution of v_A on the shock surface. In theory, the fundamental frequency of the type II radio bursts can be deemed as the Langmuir frequency of the shock upstream (Cairns & Melrose 1985; Cairns 1986, 1988; Knock & Cairns 2005). With the relation between the Langmuir frequency f and the number density n_e, $f=8980\sqrt{{n}_{{\rm{e}}}}$, we can then derive the value of n_e. Note that the frequencies corresponding to the red markers in Figure 1 are the first harmonic frequencies of the type II radio burst, thus the corresponding fundamental frequencies are approximately 173.2/2 and 150.9/2 MHz, and the corresponding n_e are 9.3 × 10⁷ and 7.0 ×10⁷ cm⁻³, respectively.

As most of the coronal density models are only in the radial direction (e.g., Newkirk 1961; Mann et al. 1999), the density variations in the longitude and latitude directions are not available. Zucca et al. (2014a) used the DEM method to calculate 2D distributions of n_e around the shock background in the corona. Rouillard et al. (2016) further updated this approach to estimate the 3D distributions of n_e around the shock background. This powerful approach has also been applied in other works (Zucca et al. 2018; Frassati et al. 2019). Based on the observations by SDO/AIA at six EUV wavelengths—93, 131, 171, 193, 211, and 335 Å—we can use the DEM method to derive the distributions of n_e in the corona (Aschwanden et al. 2001). A widely used DEM method proposed by Weber et al. (2004) is applied (Cheng et al. 2012; Su et al. 2016, 2018a) in this work. According to the approach of Rouillard et al. (2016), we obtain the distributions of n_e in the corona every 6 hr using the DEM method, which is equivalent to rotating the corona in the plane of the sky by 33 each time. From 2014 March 3 to March 5, we repeat the process mentioned above eight times, covering the entire volume of the shock surface shown in Figure 7, which spans nearly 30° in longitude. Then, based on the distributions of n_e on the eight meridional planes, we get the distribution of n_e on the shock surface by interpolation, and the results are shown in Figures 9(a) and (d). The mean values of n_e of the radio source on the shock surface are 10.0 ± 0.6 × 10⁷ and 7.9 ± 0.2 × 10⁷ cm⁻³ at 09:24:42 and 09:26:18 UT, respectively.

Zoom In Zoom Out Reset image size

Based on the distributions of B₁ and ρ, we can estimate the distribution of v_A on the shock surface using Equation (3). The results are shown in Figures 9(b) and (e). Similar to the distributions of B₁, the distribution of v_A is smooth on the shock surface. The mean values of v_A of the radio source on the shock surface are 585 ± 156 and 366 ± 88 km s⁻¹ at 09:24:42 and 09:26:18 UT, respectively. The histograms of v_A on the shock surface are shown in the fourth row of Figure 8, and roughly show a skewed unimodal distribution.

In order to derive the distribution of M_A on the shock surface, we first need to get the shock speed v_sh. We select nine slices, starting from the eruption source site (marked as a white asterisk in Figure 3), across the shock front, and the neighboring slices are separated by 6°, as shown in Figure 3(a). In the clockwise sequence, we denote these nine slices as S1, S2, S3, S4, S5, S6, S7, S8, and S9. The time–distance diagrams along the slices S1–S9 are displayed in the bottom panels of Figure 3, from which we measure the shock speed v. The shock speeds along the nine directions are 878 ± 86, 900 ± 76, 955 ± 76, 985 ± 52, 977 ± 109, 1000 ±91, 977 ± 95, 909 ± 100, and 863 ± 102 km s⁻¹, respectively.

The speeds we measured from the EUV images are the speeds of the outermost edge of the 3D shock surface along the nine slices projected in the plane of the sky. Obviously, in 3D space, the outermost edge is where the line of sight is tangential to the shock surface. For each slice, we find the heliocentric coordinates (x, y, z) of each tangent point of the line of sight and the 3D shock surface (tangent point). Taking into account the projection effect, and combining the shock speed (v_sh) measured from the time–distance diagrams (the bottom panels of Figure 3), we get the revised shock speeds, v, along each direction, based on the heliocentric coordinates (x, y, z) of each point. According to the shape of the shock surface (Equation (1)) and the coordinate transformation relationship (Equation (2)), we can get the value of h (defined in Equation (1)) for each intersection point through the heliocentric coordinates (x, y, z) of each tangent point. We assume that the shock speeds of the points on the shock surface with the same h are the same. Then, through linear interpolation, we can get the distribution of the shock speed on the shock surface.

With the distributions of v and v_A , we can estimate the distribution of M_A = v/v_A on the shock surface, and the results are shown in Figures 9(c) and (f). Similar to the distributions of B₁ and v_A , the distribution of M_A is smooth on the shock surface as well. The mean values of M_A of the radio source on the shock surface are 1.64 ± 0.32 and 2.61 ± 0.47 at 09:24:42 and 09:26:18 UT, respectively. The histograms of M_A of the radio source on the shock surface at 09:24:42 and 09:26:18 UT are shown in Figure 8. Similarly, the distributions of M_A show a roughly unimodal distribution, with asymmetry.

In theory, an MHD shock can be described by Rankine–Hugoniot (R–H) jump conditions, which result from the conservation of mass, momentum, energy, and magnetic flux (Priest 2014). The plasma parameters of the shock downstream can be determined in terms of the plasma parameters of the shock upstream using R–H relations (Ruan et al. 2018). The relation between the Alfvén Mach number M_A and the compression ratio (X) can be derived from the R–H relations, which are, for a perpendicular shock, ${M}_{{\rm{A}}\perp }=\sqrt{X(5+X)/[2(4-X)]}$, and, for a parallel shock, ${M}_{{\rm{A}}\parallel }=\sqrt{X}$. For a general oblique shock, M_A∠ can be given under a first-order approximation as follows (Bemporad & Mancuso 2011, 2013):

$\begin{eqnarray}&&{M}_{{\rm{A}}\angle }=\sqrt{{\left({M}_{{\rm{A}}\perp }\sin {\theta }_{\mathrm{Bn}}\right)}^{2}+{\left({M}_{{\rm{A}}\parallel }\cos {\theta }_{\mathrm{Bn}}\right)}^{2}},\end{eqnarray} \tag{ 4 }$

where θ_Bn is the angle between the shock normal and the upstream magnetic field vector, which has been determined in Section 3.2. Thus, combined with θ_Bn, X can be derived from M_A∠. The distributions of X on the shock surface at 09:24:42 and 09:26:18 UT are shown in Figures 10(a) and (c). The mean values of X of the radio source on the shock surface are 2.10 ± 0.53 and 2.99 ± 0.45 at 09:24:42 and 09:26:18 UT, respectively. The histograms of X of the radio source on the shock surface at 09:24:42 and 09:26:18 UT are shown in Figures 8(a) and (b).

Zoom In Zoom Out Reset image size

From the R–H jump conditions (Priest 2014), the ratio between the shock downstream and the upstream magnetic field can be given as:

$\begin{eqnarray}&&\left|\displaystyle \frac{{B}_{2}}{{B}_{1}}\right|={\left[\cos {\theta }_{\mathrm{Bn}}^{2}+{\left(\displaystyle \frac{{M}_{A}^{2}-1}{{M}_{A}^{2}-X}\right)}^{2}{X}^{2}\sin {\theta }_{\mathrm{Bn}}^{2}\right]}^{1/2},\end{eqnarray} \tag{ 5 }$

where the subscripts 1 and 2 denote the values in the up- and downstream regions, respectively. With the known θ_Bn, M_A , X, and B₁ from above, the magnetic field strength at shock downstream, B₂, can be calculated from Equation (5). The distributions of B₂ on the shock surface at 09:24:42 and 09:26:18 UT are shown in Figures 10(b) and (d). The mean values of B₂ of the radio source on the shock surface are 6.25 ± 0.89 and 4.09 ± 0.98 at 09:24:42 and 09:26:18 UT, respectively. The histograms of B₂ on the shock front at 09:24:42 and 09:26:18 UT are shown in Figures 8(a) and (b).