Magnetohydrodynamic Nonlinearities in Sunspot Atmospheres: Chromospheric Detections of Intermediate Shocks

The region of interest for the identification of dynamic sunspot phenomena encompasses both umbral and penumbral locations. The observed sunspot is very large, and as a result, weak photon flux is present in the photospheric Fe i λ6173 spectral line at the central core of the umbra; hence, we exclude this location from subsequent analyses with the Ca ii λ8542 data to be on the safe side and avoid any possible effects of low signal-to-noise ratio (see the hatched region in Figure 1 of Stangalini et al. 2018). Following common convention, all Stokes profiles are normalized to the average Stokes I continuum intensity, I_c, providing values of $I/{I}_{c}$, $Q/{I}_{c}$, $U/{I}_{c}$, and $V/{I}_{c}$ for subsequent study.

In previous chromospheric sunspot umbral investigations, shock locations are normally identified through the application of running mean subtraction methodologies in combination with observing intensity variations in the blue wing of the spectral profiles above a given threshold (Rouppe van der Voort et al. 2003; Madsen et al. 2015). However, to exploit the imaging spectropolarimetry provided by IBIS, we adopt a distinctly different approach in the identification of dynamically evolving sunspot features. Here we set the following criteria for selecting "active" pixels:

• 1.  
  A spectropolarimetric reversal in any of the Stokes $Q/{I}_{c}$, $U/{I}_{c}$, or $V/{I}_{c}$ profiles needs to be identified from one frame to the next. Examples of such spectropolarimetric flips are shown in Figure 2, with a two-dimensional map of their locations shown in the middle left panel of Figure 1.
• 2.  
  There needs to be a distinct increase, ${\rm{\Delta }}I$, in the integrated Stokes $I/{I}_{c}$ intensity originating from within the pixel location, with ${\rm{\Delta }}I\gt 3\sigma$ set as a threshold to distinguish quiescent pixels from their active counterparts, where σ is the standard deviation of the integrated Stokes $I/{I}_{c}$ fluctuations for that pixel across all time. Integrated Stokes $I/{I}_{c}$ intensities, rather than Stokes $I/{I}_{c}$ values at a particular wavelength, are used to mitigate against variable Doppler shifts producing intensity fluctuations at different wavelengths. Similar changes in the Stokes $I/{I}_{c}$ profiles have been observed in MHD shocks in previous Ca ii λ8542 investigations (de la Cruz Rodríguez et al. 2013; Grant et al. 2018). A two-dimensional representation of the locations of large ${\rm{\Delta }}I$ fluctuations is shown in the middle right panel of Figure 1.

Employing these criteria ensured that all "active" pixel detections were statistically significant and not a result of small-amplitude waves or instrumental noise. Pixels that satisfied the individual criteria mentioned above were then cospatially and cotemporally cross-correlated to identify the locations and times when spectropolarimetric reversals and large intensity fluctuations were observed simultaneously. The blue contour in the right panel of Figure 1 displays the time-integrated boundary that encompasses all established "active" pixels, with 3482 individual pixels identified over the ∼180-minute observational period. Due to the heightened ($\gt 3\sigma$) emission found in the blue wing of the Ca ii λ8542 spectral line, the isolated "active" pixels are likely to correspond to shocked plasma, similar to that identified by Houston et al. (2018) and Grant et al. (2018), only now with prominent and simultaneous spectropolarimetric reversals.

Zoom In Zoom Out Reset image size

The Non-LTE Inversion Code using the Lorien Engine (NICOLE; Socas-Navarro et al. 2015) was used to examine the changes in atmospheric parameters when transitioning from a pre-active to active state. NICOLE is a parallelized inversion code that can invert large data sets, solving multilevel, non-LTE radiative transfer problems using the pre-conditioning methods outlined in Socas-Navarro & Trujillo Bueno (1997). The inversion process requires an initial input model atmosphere, containing physical parameters such as temperature, line-of-sight (LOS) velocity, magnetic field, gas pressure, density, and microturbulence. The initial model used was the "M" sunspot model of Maltby et al. (1986), which is then perturbed to minimize the difference between the observed and synthetic Stokes profiles.

To prepare the observational data for use with NICOLE, the normalized Stokes $I/{I}_{c}$, $Q/{I}_{c}$, $U/{I}_{c}$, and $V/{I}_{c}$ profiles were interpolated onto a more dense wavelength grid (41 points; ${\rm{\Delta }}\lambda =30$ mÅ). This allows NICOLE to better fit the synthetic spectra and enables the use of the cubic DELO-Bezier formal solver outlined in de la Cruz Rodríguez & Piskunov (2013). To ensure that interpolated points, i.e., those not corresponding to a physically observed wavelength, do not contribute to the synthetic outputs, such points were assigned a negligible weighting. The Ca ii atom used consists of five bound levels plus a continuum as detailed in the works of Shine & Linsky (1974) and Socas-Navarro et al. (2000), with inversions carried out following methods outlined in previous studies (e.g., Henriques et al. 2017; Kuridze et al. 2018). The effect of Ca ii isotropic splitting was also included in the inversions (Leenaarts et al. 2014).

The nodes used for the calculation of each parameter are equally spaced along the optical depth scale at 500 nm (${\mathrm{log}}_{10}\tau$). Perturbations to the background model are applied at these locations, with the correction between them performed using cubic Bezier interpolation. NICOLE inversions are computationally intensive, although fortunately the number of active pixels identified in our data set was relatively small. In total, 6964 pixels were inverted (3482 pre-active and 3482 active). To improve the fit of the synthetic profiles to the observations, the inversions were carried out in three cycles, with subsequent cycles having increased numbers of nodes to improve the quality of convergence, as suggested by Ruiz Cobo & del Toro Iniesta (1992). The node points used for each cycle are summarized in Table 1. In between the first and second cycle the atmosphere was smoothed, both horizontally and vertically, with the smoothing only performed on perturbations between the input and generated atmospheres. The atmospheres from the previous cycle were subsequently used as inputs for the next inversion cycle. Throughout all cycles a weighting of 1 is applied to Stokes $I/{I}_{c}$ and $V/{I}_{c}$ profiles, as this resulted in the best synthetic profile fits. In the initial cycle we include one node for the transverse components of the magnetic field, ${{\boldsymbol{B}}}_{x}$ and ${{\boldsymbol{B}}}_{y}$, due to the need to generate accurate Stokes $I/{I}_{c}$ and $V/{I}_{c}$ fits from which to base the following cycles off. The weights across cycles relative to Stokes $I/{I}_{c}$ in $Q/{I}_{c}$ and $U/{I}_{c}$ were 0.2, 1, and 5, respectively, while in $V/{I}_{c}$ they were weighted the same as $I/{I}_{c}$ in each cycle. We apply such weights to Stokes $Q/{I}_{c}$ and $U/{I}_{c}$ to better constrain the transverse component of the magnetic field, since this is an important parameter in the classification of umbral shocks and dynamics (Houston et al. 2018). The weighting of $I/{I}_{c}$ and $V/{I}_{c}$ was kept the same throughout, as this resulted in the best overall fits across all Stokes profiles. Figure 2 displays the Stokes $I/{I}_{c}$, $Q/{I}_{c}$, $U/{I}_{c}$, and $V/{I}_{c}$ profiles corresponding to a pixel in a quiescent phase (left panels) and that same pixel during a shock (right panels). The black lines represent the observed spectra obtained from the IBIS instrument, with the red dashed lines showing the best-fit synthetic profiles generated from the NICOLE inversion process. The shaded regions indicate the spatially and temporally averaged standard deviations corresponding to offsets between the input and synthetic profiles. The confinement of the wavelength-dependent standard deviations (blue shaded regions in Figure 2) shows statistically a high degree of precision throughout the inversion process, something that is also highlighted by Houston et al. (2018).

Table 1.  Number of Nodes Used for Each Cycle of the NICOLE Inversions

Physical Parameter	Cycle 1	Cycle 2	Cycle 3
Temperature	3	5	7
LOS Velocity	1	3	3
${{\boldsymbol{B}}}_{x}$	1	2	3
${{\boldsymbol{B}}}_{y}$	1	2	3
${{\boldsymbol{B}}}_{z}$	1	2	3
Microturbulence	1	1	1

Download table as:  ASCIITypeset image

Figure 3 displays the fractional uncertainties for the derived NICOLE parameters. The uncertainties were determined by performing numerous inversions on 100 randomly extracted pixels from our list of active locations using different initial conditions. The mean standard deviation across all pixels, at all optical depth points, for the different initial models was then determined, with the values normalized by their respective parameter mean to generate a fractional uncertainty. The calculated fractional uncertainties, for each parameter across the optical depth range spanning $-6.0\lt {\mathrm{log}}_{10}\tau \lt -2.0$, is shown using colored shaded regions in Figure 3.

Zoom In Zoom Out Reset image size

The left panel of Figure 4 details the changes experienced by the Ca ii λ8542 line core Doppler velocity when transitioning from a quiescent atmosphere to a shock environment. Each of the data points displays the relationship between the quiescent (i.e., pre-shock; x-axis) and active (i.e., shock; y-axis) states, where the color represents the optical depth at which the measurement was made. Here an optical depth of ${\mathrm{log}}_{10}\tau =-2$ corresponds to the photosphere, while ${\mathrm{log}}_{10}\tau =-6$ is indicative of upper chromospheric locations. We note that the atmospheric solution contains larger uncertainties in the ${\mathrm{log}}_{10}\tau =-6$ regime (Quintero Noda et al. 2016), where the layers could be affected by extrapolation effects from gradients deeper down in the atmosphere. Note that the sign of the Doppler velocities is linked to the induced wavelength shift, whereby positive values represent redshifted (i.e., downflowing) plasma and negative values correspond to blueshifted (i.e., upflowing) plasma.

Zoom In Zoom Out Reset image size

Inspection of the results shows that the plasma is predominantly redshifted immediately prior to the formation of a shock, but this changes abruptly to blueshifted material during the development of the shock itself. Such a change is consistent with previous MHD shock studies, including those linked to magnetoacoustic (Joshi & de la Cruz Rodríguez 2018; Anan et al. 2019) and resonantly amplified fast-mode (Grant et al. 2018) nonlinearities. This trend is consistent across all optical depths, although the magnitudes of the velocities increase with atmospheric height, as expected, due to the reduced plasma pressure in these locations. At optical depths of ${\mathrm{log}}_{10}\tau =-6$ and −5, it has been seen that the vector shock velocities are averaged (de la Cruz Rodríguez et al. 2012), resulting in an underestimation of the true shock speed. In this study, we focus on optical depths ${\mathrm{log}}_{10}\tau =-4$ and −3, which we observe to closely match the simulations; however, we note that the derived speed is likely to be a lower limit of the true shock velocity.

The right panel of Figure 4 displays the changes in temperature, ${\rm{\Delta }}T$, associated with the transition from a pre-shock phase to a shock environment as a function of the quiescent plasma temperature. The derived temperature changes are in the range of $-250\,{\rm{K}}\lt {\rm{\Delta }}T\lt 2500\,{\rm{K}}$. Generally, a greater temperature perturbation is provided to plasma with a cooler quiescent temperature at each optical depth. In the right panel of Figure 4 this is particularly visible at an optical depth of ${\mathrm{log}}_{10}\tau =-5$ (purple data points), where quiescent plasma with temperatures of around 3000 K experiences ${\rm{\Delta }}T\sim 1500$ K, while background temperatures of 5000 K only provide ${\rm{\Delta }}T\sim 250$ K.

As a result, we suggest that the plasma shocks identified have the ability to contribute to local atmospheric heating when a sufficient temperature gradient is present, with a ${\rm{\Delta }}T=0$ K value suggesting equilibrium between the developing shock and the quiescent background. To formalize the background temperature at each optical depth that promotes a ${\rm{\Delta }}T=0$ K equilibrium, we fit a linear trend line through the data points spanning each optical depth and calculate the intersection of the best-fit line with the x-axis. This provides shock/background equilibrium temperatures of ∼3690, ∼3650, ∼3465, ∼6020, and ∼10,845 K for optical depths corresponding to ${\mathrm{log}}_{10}\tau =-2$ (∼250 km; Maltby et al. 1986), ${\mathrm{log}}_{10}\tau =-3$ (∼450 km), ${\mathrm{log}}_{10}\tau =-4$ (∼625 km), ${\mathrm{log}}_{10}\tau =-5$ (∼1150 km), and ${\mathrm{log}}_{10}\tau =-6$ (∼1850 km), respectively. Interestingly, the largest average ΔT (+19.6%) occurs at an optical depth of ${\mathrm{log}}_{10}\tau =-5$, corresponding to an approximate geometric height of 1150 km (Maltby et al. 1986), which is consistent with examinations of resonantly amplified fast-mode shocks, where Grant et al. (2018) found the largest temperature perturbations to be within the range of $-5.3\,\lt {\mathrm{log}}_{10}\tau \lt -4.6$. With this in mind, the pixels we identify as "active" have clear similarities to previously detected MHD shock phenomena, with characteristics related to the temperature and LOS velocity perturbations closely resembling the signatures synonymous with magnetoacoustic (e.g., de la Cruz Rodríguez et al. 2013; Houston et al. 2018) and fast-mode (Grant et al. 2018) shocks. However, we can now employ the high-precision vector magnetic fields to further categorize the underlying shock behavior.

In a similar manner to how the temperature fluctuations are depicted in the right panel of Figures 4 and 5 displays the shock vector magnetic fields as a function of their pre-shock values, with the transverse (${{\boldsymbol{B}}}_{x}$ and ${{\boldsymbol{B}}}_{y}$) and vertical (${{\boldsymbol{B}}}_{z}$) components displayed in the top, middle, and bottom panels, respectively. For completeness, the x- and y-directions represent the two orthogonal directions within the plane coincident with the solar surface, with x representing the east–west direction and y the north–south direction with respect to the heliographic coordinate system. Examining the transverse magnetic field fluctuations (top and middle panels of Figure 5), it is clear that a reversal occurs between the quiescent (i.e., pre-shock) environment and the shocked plasma state. This reversal is indicated by the gradient of the best-fit lines (dotted red lines in the top and middle panels of Figure 5) being close to −1, with gradients of −0.99 and −1.00 found for the ${{\boldsymbol{B}}}_{x}$ and ${{\boldsymbol{B}}}_{y}$ fits, respectively. To make use of the large number statistics at our disposal, the relationships illustrated in Figure 5 for ${{\boldsymbol{B}}}_{x}$ and ${{\boldsymbol{B}}}_{y}$ are further examined through the calculation of their corresponding Spearman's rank correlation coefficients (Spearman 1904), where the probabilistic p-values are calculated under the following hypotheses: (1) there is zero linear correlation between the pre- and post-shock variables, and (2) the correlation coefficient is not equal to zero. The Spearman's rank coefficients shown in Table 2 highlight a strong negative linear association between pre- and post-shock values that is statistically significant at the 5% level, hence reiterating the strong anticorrelation found between pre- and post-shock transverse magnetic field fluctuations.

Zoom In Zoom Out Reset image size

Table 2.  Spearman Rank Correlation Coefficients for Pre- and Post-shock Locations Corresponding to the Transverse Magnetic Field Components, ${{\boldsymbol{B}}}_{x}$ and ${{\boldsymbol{B}}}_{y}$

	Correlation	p-value	95% Confidence Interval
${{\boldsymbol{B}}}_{x}$	−0.8402	<0.0001	(−0.845, −0.835)
${{\boldsymbol{B}}}_{y}$	−0.8143	<0.0001	(−0.820, −0.809)

Download table as:  ASCIITypeset image

Of course, the exceptionally clear trends depicted here may not come as a complete surprise, since our pixel identification methodology required a spectropolarimetric reversal in the observed Stokes profiles (see, e.g., Figure 2). It must be noted that the grouping of the data points is slightly more extended (i.e., more measurements reaching larger magnetic field strengths) for the ${{\boldsymbol{B}}}_{x}$ plot in the top panel of Figure 5, when compared with the ${{\boldsymbol{B}}}_{y}$ scatterplot depicted in the middle panel of Figure 5. This is a consequence of the identified pixels predominantly residing on the eastern side of the sunspot (Figure 1), where ${{\boldsymbol{B}}}_{x}$ magnitudes will be strongest owing to the natural orientation of the magnetic fields along that direction. However, to account for the spread in the data, robust linear regression (Lange et al. 1989), assuming t-distributed residuals, was utilized to provide a better fit to the heavy tails of the data distribution. Similar gradients and confidence intervals (CI) were found for ${{\boldsymbol{B}}}_{x}$ (−0.99; CI [−1.00, −0.98]) and ${{\boldsymbol{B}}}_{y}$ (−1.00; CI [−1.01, −0.99]) comparisons, reiterating the strong negative association between pre- and post-shock values that are highly significant.

Interestingly, no such dominant polarity reversal is identified in the ${{\boldsymbol{B}}}_{z}$ component of the magnetic field. As displayed in the bottom panel of Figure 5, the signs and magnitudes of the ${{\boldsymbol{B}}}_{z}$ terms are consistent between quiescent and shocked states. This can be seen through the relatively tight grouping of data points along the 1:1 linear trend plotted as a red dotted line in the bottom panel of Figure 5. It might be natural to assume that a spectropolarimetric reversal in Stokes $V/{I}_{c}$ (see, e.g., the bottom right panel of Figure 2) would indicate a polarity reversal in that particular pixel location. However, spectropolarimetric reversals have been witnessed previously by de la Cruz Rodríguez et al. (2013), with such signatures not necessarily implying a physical reversal of the magnetic field polarities, as is also implied in the bottom panel of Figure 5. Joshi & de la Cruz Rodríguez (2018) found that magnetic field perturbations are not the result of opacity changes, and therefore the observed reversals are not consistent with opacity effects. Instead, the spectropolarimetric reversals found in Stokes $V/{I}_{c}$ may be the consequence of a developing shock creating a two-component atmosphere, where independent bulk motions of the peak opacity-forming regions give rise to shifts in the polarimetric profiles, similar to that observed by Socas-Navarro et al. (2000).

From examination of Figure 5, it is clear that under the development of a shock the ${{\boldsymbol{B}}}_{x}$ and ${{\boldsymbol{B}}}_{y}$ values flip, while the ${{\boldsymbol{B}}}_{z}$ component of the magnetic field remains approximately constant with the same sign. Due to the very pronounced reversals in the ${{\boldsymbol{B}}}_{x}$ and ${{\boldsymbol{B}}}_{y}$ components (i.e., best-fit line gradients very close to −1 in the top and middle panels of Figure 5), fluctuations in the ${{\boldsymbol{B}}}_{x}$ and ${{\boldsymbol{B}}}_{y}$ terms should only contribute to very minor changes in the total magnetic field strength, ${{\boldsymbol{B}}}_{\mathrm{tot}}$. This implies that any changes in ${{\boldsymbol{B}}}_{z}$ should produce a near-equivalent fluctuation in ${{\boldsymbol{B}}}_{\mathrm{tot}}$—i.e., Δ${{\boldsymbol{B}}}_{z}$ = Δ${{\boldsymbol{B}}}_{\mathrm{tot}}$. Figure 6 displays a scatterplot detailing the relationship between Δ${{\boldsymbol{B}}}_{z}$ and Δ${{\boldsymbol{B}}}_{\mathrm{tot}}$, where the dotted red line highlights a linear best-fit line between the two variables. The gradient associated with the best-fit line is 0.92 ± 0.04, indicating a very close correlation between fluctuations in the vertical component of the magnetic field (${{\boldsymbol{B}}}_{z}$) and the total overall magnetic field strength (${{\boldsymbol{B}}}_{\mathrm{tot}}$).

Zoom In Zoom Out Reset image size

It must be noted that the NICOLE spectropolarimetric inversions performed in this study harness the Zeeman effect to estimate the various magnetic field parameters. As a result, Zeeman-induced Stokes inversions produce a 180° azimuthal ambiguity in the direction of the transverse magnetic field. Therefore, identical observational Stokes I/Q/U/V profiles can produce inversion outputs of either, for example, +${{\boldsymbol{B}}}_{x}$/+${{\boldsymbol{B}}}_{y}$ or −${{\boldsymbol{B}}}_{x}$/−${{\boldsymbol{B}}}_{y}$. Hence, care needs to be taken when examining the outputs of NICOLE spectropolarimetric inversions, as a flip (e.g., +${{\boldsymbol{B}}}_{x}$ $\to -$ ${{\boldsymbol{B}}}_{x}$ and +${{\boldsymbol{B}}}_{y}$ $\to -$ ${{\boldsymbol{B}}}_{y}$) in the transverse magnetic field could be a consequence of this Zeeman-based ambiguity and not a result of a physical change in the Sun's vector magnetic field. However, in the present work we strive to alleviate this concern by implementing stringent selection criteria for our active pixels (see Section 3.1), which requires a distinct reversal in the observed Stokes profiles, hence highlighting pixel locations where shock-induced morphological changes are indeed present. Of course, even with clear reversals in the observed Stokes profiles, the Zeeman-induced ambiguity associated with the subsequent NICOLE inversions may provide incorrect transverse magnetic field information (e.g., +${{\boldsymbol{B}}}_{x}$ remains +${{\boldsymbol{B}}}_{x}$ and +${{\boldsymbol{B}}}_{y}$ remains +${{\boldsymbol{B}}}_{y}$). This may be the cause of some small positive correlations seen in Figure 5, particularly at relatively weak magnetic field strengths ($| {{\boldsymbol{B}}}_{x/y}| \lesssim 500$ G).

While the selection criterion (see Section 3.1) for Stokes profiles helps to identify regions of the solar atmosphere that are undergoing shock-induced morphological changes, it does not provide spectral constraints related to the dynamic source functions in these rapidly evolving locations. In particular, a major challenge facing both observers and theoreticians is to understand how changing gradients of the source function (e.g., when a shock causes an absorption Ca ii λ8542 spectral line to transition into emission) also effects the subtle variations seen in optically thick chromospheric Stokes Q/U/V spectra (e.g., López Ariste et al. 2001; de la Cruz Rodríguez et al. 2013; Joshi & de la Cruz Rodríguez 2018). Numerical modeling by Felipe et al. (2014) demonstrated how synthetically generated Ca ii λ8542 spectra, following the creation of magnetoacoustic shocks, often displayed Stokes Q/U/V reversals—a consequence of the source function no longer monotonically decreasing throughout the chromosphere (see also the Stokes V spectral discussions by de la Cruz Rodríguez et al. 2013). However, the spectra generated by Felipe et al. (2014) are further complicated by the presence of additional turning points within the Stokes Q/U/V profiles, something that is not observed in our identified active pixel locations. Such modeling endeavors are at the forefront of current solar physics research, and as a result, more detailed radiative transfer calculations need to be undertaken in order to isolate the specific mechanism(s) responsible for Stokes Q/U/V reversals witnessed during shock formation in the solar chromosphere. As such, with the present data set, we cannot exclude the possibility that some of the Stokes Q/U/V reversals captured may be emphasized to a degree by variations in the associated gradients of the contributing source function.

Future examinations of the magnetic field perturbations caused by shocks may wish to make use of both Zeeman and Hanle diagnostics to minimize ambiguities caused by the inversion process (Asensio Ramos et al. 2008; Centeno 2018). Furthermore, to disambiguate the cause of the reversals within our Stokes Q/U/V profiles requires higher-sensitivity polarimetric measurements of the corresponding spectra. In conjunction, improved modeling of radiative transfer effects within the optically thick lower atmosphere will be required, since Grant et al. (2018) have shown that the origin of developing shocks within sunspot umbrae can span more than 1000 km in geometric height, which likely has implications for the subtle variations captured in the associated spectral profiles. However, we must emphasize that our selection criterion for active pixels requires an observed and measurable reversal of the spectropolarimetric Stokes profiles captured by IBIS and, as a consequence, does not rely solely on the transverse magnetic field outputs from the NICOLE inversions.

The final plasma parameter to examine is the density, with histograms of density fluctuations, related to quiescent and shock environments for different optical depths, shown in Figure 7. The histograms reveal the percentage changes in NICOLE-derived densities resulting from shock formation. We note that NICOLE computes the gas stratification assuming hydrostatic equilibrium, with the exclusion of the Lorentz force and advection term. With the mainly vertical magnetic fields present within the sunspot, the assumption of hydrostatic equilibrium is likely to be a robust approximation. Future work may wish to consider non-steady-state model atmospheres, which would provide more accurate flow field information across a broader range of optical depths, which will be important to unequivocally constrain the magnitudes of the Lorentz and advection terms. It can be clearly seen that at ${\mathrm{log}}_{10}\tau =-3$ and ${\mathrm{log}}_{10}\tau =-4$ (high photosphere and low chromosphere, respectively) there are shock formation signatures, with increase in local densities of approximately 10%–20%. This identifies the layers of the solar atmosphere where the local plasma has been substantially compressed by the shock development. The induced density fluctuations begin to reduce at optical depths around ${\mathrm{log}}_{10}\tau =-5$, while at the extreme upper boundary of the chromosphere (${\mathrm{log}}_{10}\tau =-6$) a decrease in shock-induced density is found. This is a possible consequence of the shock developing in the upper photosphere/lower chromosphere, with the signatures becoming diluted as they traverse multiple density scale heights, where the density scale height in the chromosphere is ∼300 km (Peter & Marsch 1998). It may also be a consequence of an increase in adiabatic pressure resulting from the shock, which causes the magnetic fields to expand at higher atmospheric heights where there is less plasma pressure, hence reducing the local density (Houston et al. 2018). Finally, the generation of Prandtl–Meyer expansion fans (Chen & Feldman 2015; Cao et al. 2017) as the supersonic plasma associated with the shock interacts with the geometry of the magnetic field may initiate Mach waves, which subsequently produce low-density wakes as they traverse through the upper layers of the chromosphere, hence manifesting as decreased density perturbations at optical depths of ${\mathrm{log}}_{10}\tau \sim -6$. However, this area of research is in its infancy and requires dedicated shock-capturing numerical simulations (e.g., using the Lagrangian–Eulerian Remap code; Arber et al. 2001 or MPI-AMRVAC code; Porth et al. 2014) to further test the effects of such phenomena.

Zoom In Zoom Out Reset image size

The RH relations are typically expressed in the rest frame of the shock and are in their most condensed form when exploiting the de Hoffman–Teller frame, where the magnetic and velocity components are coplanar and aligned in both the quiescent and shocked states. However, we are unable to deduce the true vector velocity field before and during shock activation, since while the spatial sampling and temporal cadence of the IBIS data products are relatively high, the rapid evolution and creation of two-component atmospheres are prohibitive without having to rely on additional unconstrained assumptions. Hence, we employ the spectropolarimetric inversions, which provide us with vector magnetic fields and thermodynamic information, to further classify the captured shock activity.

The vector magnetic fields, together with the plasma temperatures and densities, allow us to calculate the local plasma-β values, where β is the ratio between the plasma gas pressure and the pressure of the magnetic field, defined as

$$\begin{eqnarray*}&&\beta =\displaystyle \frac{2{\mu }_{0}{n}_{{\rm{H}}}{{Tk}}_{{\rm{B}}}}{{{\boldsymbol{B}}}_{\mathrm{tot}}^{2}},\end{eqnarray*}$$

where μ₀ is the magnetic permeability and n_H is the hydrogen number density. Locations where the plasma-β are close to unity are important in the studies of wave propagation, since they offer more efficient regions for mode coupling and resonant amplification of the embedded wave amplitudes (Zaqarashvili & Roberts 2006; Zaqarashvili et al. 2006; Cally & Goossens 2008). Recently, their importance has also been demonstrated for the generation of resonantly driven fast-mode shocks toward the edges of sunspot umbrae (Grant et al. 2018). In the right panel of Figure 1 we display the plasma-β = 1 isocontours spanning the optical depths (inner contour; ∼450 km) $-3\gt {\mathrm{log}}_{10}\tau \gt -4$ (outer contour; ∼625 km), where the detected shock activity is believed to first manifest. From Figure 1 it can be seen that the active pixels are predominantly contained within the plasma-β = 1 isocontours, suggesting that these locations may play a crucial role in the development of shock phenomena displaying spectropolarimetric reversals.

From quiescent to shocked states, the calculated plasma-β values systematically become larger, remaining consistent with a shock-induced increase in the local plasma pressure. We are able to determine the sound, v_S, and Alfvén, v_A, speeds in both quiescent and shocked states using the relationships

$$\begin{eqnarray*}&&{v}_{S}^{2}=\gamma p/\rho ,\end{eqnarray*}$$

$$\begin{eqnarray*}&&{v}_{{\rm{A}}}^{2}={{\boldsymbol{B}}}_{\mathrm{tot}}^{2}/{\mu }_{0}\rho .\end{eqnarray*}$$

To estimate the shock normal direction, ${\hat{{\boldsymbol{n}}}}_{{\rm{s}}}$, we quantify the angles that give its orientation. The position of the pixel (x_p, y_p) with respect to sunspot center is used to calculate the angle φ_s between the x-axis and the projection of ${\hat{{\boldsymbol{n}}}}_{{\rm{s}}}$ on the horizontal plane,

$$\begin{eqnarray*}&&{\varphi }_{s}=\arctan \displaystyle \frac{{x}_{p}}{{y}_{p}}.\end{eqnarray*}$$

Then, using the first RH condition, which states that the normal magnetic field component, ${{\boldsymbol{B}}}_{{\boldsymbol{n}}}$, remains constant, i.e., ${{\boldsymbol{B}}}_{{\boldsymbol{n}},a}={{\boldsymbol{B}}}_{{\boldsymbol{n}},b}$, we calculate the angle ${\vartheta }_{s}$ between ${\hat{e}}_{z}$ (the vertical) and ${\hat{{\boldsymbol{n}}}}_{{\rm{s}}}$,

$$\begin{eqnarray*}&&{\vartheta }_{s}=\arctan \displaystyle \frac{-{B}_{z,b}+{B}_{z,a}}{\left({B}_{x,b}-{B}_{x,a}\right)\cos {\varphi }_{s}+\left({B}_{y,b}-{B}_{y,a}\right)\sin {\varphi }_{s}}.\end{eqnarray*}$$

Once we have ${\hat{{\boldsymbol{n}}}}_{{\rm{s}}}$, we can decompose the magnetic field, ${\boldsymbol{B}}$, into its normal, ${{\boldsymbol{B}}}_{{\boldsymbol{n}}}$, and tangential, ${{\boldsymbol{B}}}_{{\rm{t}}}={\boldsymbol{B}}-{{\boldsymbol{B}}}_{{\boldsymbol{n}}}$ ${\hat{{\boldsymbol{n}}}}_{{\rm{s}}}$, components. This also provides us with the total magnetic field jump entirely in the tangential direction. This allows us to quantify the angle, θ, between the vector magnetic field, ${\boldsymbol{B}}$, and the shock normal, ${\hat{{\boldsymbol{n}}}}_{{\rm{s}}}$, in both quiescent and shocked states through the relation

$$\begin{eqnarray*}&&\theta =\arctan \left(\displaystyle \frac{{{\boldsymbol{B}}}_{{\rm{t}}}}{{{\boldsymbol{B}}}_{{\boldsymbol{n}}}}\right).\end{eqnarray*}$$

The shock adiabatic relation (Anderson 1963) provides the propagation speed of the shock as a function of its strength and the composition of the upstream parameters. From the outputs of the NICOLE inversions and subsequent calculation of the upstream parameters (${v}_{S},{v}_{{\rm{A}}},\theta$), we are able to compute the shock solutions corresponding to the three possible pre-shock normal speeds (slow, fast, and Alfvén).

Figure 8 displays the shock normal speeds for a typical pixel capturing a shock. The shaded regions represent the averaged standard deviations corresponding to offsets between the shock normal solutions when the most extreme input parameters are propagated through the calculations. We see in the top panels of Figure 8 that when the plasma density ratio is altered by within the statistical uncertainty range from the inversions, there is no significant change in the solutions, and for the shock strength of the given pixel (dashed black line) there is no overlap between the three solutions. The bottom panels of Figure 8 display the uncertainties arising from a change in the magnitude of the magnetic field. Again, it is clear that there is no overlap in the solution output for this shock strength. The narrowband regions show the statistical significance of the solutions. The parameters attained from the inversion process are sufficiently accurate to not affect the final shock classification, with no overlap present in the shock solutions for the slow, Alfvén, or fast cases at the shock strength of the event.

Zoom In Zoom Out Reset image size

We can then verify which of the three solutions comes closest to obeying the underlying RH conditions. First, we compute the post-shock normal velocity from the mass flux continuity requirement (Gosling et al. 1968),

$$\begin{eqnarray*}&&{v}_{{\boldsymbol{n}},a}=\displaystyle \frac{{v}_{{\boldsymbol{n}},b}\,{\rho }_{b}}{{\rho }_{a}}.\end{eqnarray*}$$

This relationship states that $\rho {v}_{{\boldsymbol{n}}}$ is identical in pre- and post-shock states. As solutions to the shock adiabatic relation, the corresponding Alfvén mach number (${M}_{{\rm{A}}}={v}_{{\boldsymbol{n}}}\sqrt{\rho }/{B}_{{\boldsymbol{n}}}$) stays below unity in both pre- and post-shock states for the "slow" solution, remains above unity for the "fast" solution, and jumps across the M_A = 1 boundary (i.e., super-Alfvénic to sub-Alfvénic) for the "Alfvén" solution.

With the normal velocity fields calculated, we are able to subsequently test the following three RH relations:

• 1.  
  The tangential magnetic field in pre- and post-shock states must be parallel according to

  $$\begin{eqnarray*}&&{{\boldsymbol{B}}}_{{\rm{t}},b}\left(\displaystyle \frac{{(\rho {v}_{{\boldsymbol{n}}})}^{2}}{{\rho }_{b}}-{{\boldsymbol{B}}}_{{\boldsymbol{n}}}^{2}\right)={{\boldsymbol{B}}}_{{\rm{t}},a}\left(\displaystyle \frac{{(\rho {v}_{{\boldsymbol{n}}})}^{2}}{{\rho }_{a}}-{{\boldsymbol{B}}}_{{\boldsymbol{n}}}^{2}\right).\end{eqnarray*}$$

  As a vector relation, each of the x, y, z directional components provides us with an independent check.
• 2.  
  The momentum flux must balance according to

  $$\begin{eqnarray*}&&\left\{\rho {v}_{{\boldsymbol{n}}}^{2}+p+\displaystyle \frac{{{\boldsymbol{B}}}_{{\rm{t}}}^{2}}{2}\right\}=0.\end{eqnarray*}$$
• 3.  
  We should expect that the observed difference in the LOS velocity, ${\rm{\Delta }}{v}_{\mathrm{los}}$, correlates with the observed jump in the normal velocity, ${v}_{{\boldsymbol{n}}}$. Hence, we check the value of the quantity

  $$\begin{eqnarray*}&&{\rm{\Delta }}{v}_{\mathrm{los}}=\left|{v}_{{\boldsymbol{n}},b}-{v}_{{\boldsymbol{n}},a}-{v}_{\mathrm{los},b}+{v}_{\mathrm{los},a}\right|.\end{eqnarray*}$$

Deviations from any of these RH relations would indicate an incompatibility with that particular shock classification (i.e., slow, fast, and Alfvén). On the contrary, a set of measured parameters that provide minimal offsets (i.e., numerical values tending to zero) between the generalized RH relationships would indicate a more robust classification of the detected shock environment. As such, we investigate the level of agreement between each of the three RH relations defined above and our isolated shock pixels, with the best agreements (i.e., minimizing any offsets between what is expected and what is measured) providing important information to classify the specific type of shock present.

We evaluate the validity of each RH condition across a range of optical depths spanning the mid-photosphere through to the upper chromosphere ($-6\leqslant {\mathrm{log}}_{10}\tau \leqslant -3$). From the three RH conditions defined above, we obtain five measurements per pixel for each of the slow, Alfvén, and fast shock solutions across four discrete optical depths. This leads to 60 individual values for each pixel, where the quantities represent the offsets between the idealized RH relationships and those extracted from the observations. With 3482 shock pixels identified over the ∼180-minute observing period, this equates to more than 200,000 individual measurements that can be used to robustly identify the type of MHD shock manifesting in our observational time series.

The mean of the values for each optical depth and MHD shock type (slow, fast, and Alfvén) is calculated, with the results displayed in the form of box-and-whisker plots in Figure 9. Examination of Figure 9 shows how the Alfvén solution systematically provides the lowest numerical offsets between the observations and the three generalized RH conditions. In particular, the median value (represented by the red horizontal line in each box) is three times lower for the Alfvén solution than the corresponding slow solution at depths consistent with shock formation (${\mathrm{log}}_{10}\tau =-3,-4;$ where the entropy change, ΔS, is largest; bottom right panel of Figure 9). The offsets associated with the fast solution, at the optical depths of peak shock formation, are an order of magnitude greater than that of the Alfvén solution, providing strong evidence that the shocks observed are not fast MHD shocks. Ideally, the numerical offset values between the observations and the three RH conditions should be zero for complete certainty when characterizing the embedded shocks. However, due to instrumental, atmospheric seeing, and inversion constraints, this level of precision is not possible. Nevertheless, the much-reduced numerical offsets for the Alfvén solutions, at optical depths consistent with the formation of the shock phenomena, suggest that the detected shock behavior is best classified by Alfvén (or intermediate) MHD shocks.

Zoom In Zoom Out Reset image size

At higher geometric heights (${\mathrm{log}}_{10}\tau =-5,-6$), we observe a clear increase in the Alfvén solution offsets between the RH criteria values and those computed from the observations. At these optical depths, which are consistent with atmospheric heights pushing the upper chromosphere, the plasma still experiences perturbations in its temperature, LOS velocity, and vector magnetic field. However, these perturbations are the result of the shock forming much deeper in the solar atmosphere, with the ensuing dynamics propagating upward across multiple density scale heights and subsequently impacting the plasma conditions in the upper chromosphere. As a result, the RH conditions are not expected to be satisfied at these optical depths since the shock boundaries, where the RH relations should be evaluated, form at much lower geometric heights (i.e., ${\mathrm{log}}_{10}\tau =-3,-4$). Therefore, comparisons between the idealized RH conditions and our observational parameter space reveal that MHD shocks, demonstrating spectropolarimetric polarity inversions, form at optical depths consistent with the upper photosphere/lower chromosphere (${\mathrm{log}}_{10}\tau =-3,-4$), and that the subsequent shock dynamics are best characterized by the formation of Alfvén (or intermediate) shock types.

We have provided observational evidence of Alfvén shocks manifesting within a sunspot umbra, where the shocks have the ability to perturb the local plasma (e.g., density, temperature, magnetic field, etc.) parameters. Finding evidence of such phenomena has implications for the supply of thermal energy to chromospheric umbral regions. Recently, Anan et al. (2019) suggested that traditional magnetoacoustic shock behavior in sunspot umbrae is unable to supply sufficient thermal energy to maintain the umbral chromosphere. However, here we demonstrate that in addition to traditional magnetoacoustic shocks, there exists an abundance of Alfvén shock developments also able to provide substantial thermalization in the solar chromosphere. Snow & Hillier (2019) have suggested that such shock behavior has the potential to occur in a wide range of phenomena in the solar atmosphere in which partial ionization effects are important, including magnetic reconnection (e.g., Ellerman bombs, spicules) and wave steepening (e.g., umbral flashes). Snow & Hillier (2019) highlight that the shock effects are likely to be most significant in the lower chromospheric regions, which is consistent with what we demonstrate in the present study. The plethora of possible Alfvén shock environments implies that these events may provide a significant contribution to the overall heating requirements of the chromosphere, consistent with the ideas put forward by Matsumoto & Suzuki (2014).