Magnetohydrodynamic Non-linearities in Sunspot Atmospheres: Chromospheric Detections of Intermediate Shocks

The formation of shocks within the solar atmosphere remains one of the few observable signatures of energy dissipation arising from the plethora of magnetohydrodynamic waves generated close to the solar surface. Active region observations offer exceptional views of wave behavior and its impact on the surrounding atmosphere. The stratified plasma gradients present in the lower solar atmosphere allow for the potential formation of many theorized shock phenomena. In this study, using chromospheric Ca II 854.2nm spectropolarimetric data of a large sunspot, we examine fluctuations in the plasma parameters in the aftermath of powerful shock events that demonstrate polarimetric reversals during their evolution. Modern inversion techniques are employed to uncover perturbations in the temperatures, line-of-sight velocities, and vector magnetic fields occurring across a range of optical depths synonymous with the shock formation. Classification of these non-linear signatures is carried out by comparing the observationally-derived slow, fast, and Alfv\'en shock solutions to the theoretical Rankine-Hugoniot relations. Employing over 200,000 independent measurements, we reveal that the Alfv\'en (intermediate) shock solution provides the closest match between theory and observations at optical depths of log(tau) = -4, consistent with a geometric height at the boundary between the upper photosphere and lower chromosphere. This work uncovers first-time evidence of the manifestation of chromospheric intermediate shocks in sunspot umbrae, providing a new method for the potential thermalization of wave energy in a range of magnetic structures, including pores, magnetic flux ropes, and magnetic bright points.


INTRODUCTION
The desire to understand wave energy transportation, and subsequent dissipation in the solar atmosphere, is a major driver behind much of solar physics research. Owing to significant advancements in instrumentation, post-processing techniques and numerical simulations, our understanding of wave activity in the atmosphere has vastly improved in recent years (e.g., Roberts 2000;Nakariakov & Verwichte 2005;Bard & Carlsson 2010;Felipe 2012;Mathioudakis et al. Corresponding author: S. J. Houston shouston22@qub.ac.uk 2013; Jess et al. 2015;Felipe 2019, to name but a few). It is becoming increasingly apparent that highly magnetic solar regions, such as those associated with sunspots, pores and magnetic bright points, are able to play a prominent role in atmospheric heating (Sobotka et al. 2013;Grant et al. 2018).
In recent times there has been a drive to more closely examine the energy dissipation occurring in the solar chromosphere, with readily developing shock fronts seen as a likely mechanism for such dissipation. In a magnetohydrodynamic (MHD) framework, the propagation of slow magnetoacoustic waves, and their subsequent development into prominent shocks, has attained widespread examination since their ubiquitous detection inside highly magnetic sunspot umbrae (Beckers & Tallant 1969). This phenomenon is commonly referred to as umbral flashes (UFs), and are a consequence of the steepening of slow magnetoacoustic waves as they traverse the rapid density stratification of the umbral atmosphere Henriques et al. 2017). These events provide local temperature enhancements on the order of 1000 K in the low-to-mid chromosphere, and are able to manipulate the geometry of the embedded magnetic fields through increased adiabatic pressure (Houston et al. 2018). Indeed, slow acoustic-type shocks are so widespread that they are regularly visible in smaller scale structures, such as those associated with Ca II grains (Rutten & Uitenbroek 1991;Carlsson & Stein 1997;Cauzzi et al. 2009;Vecchio et al. 2009;Martínez-Sykora et al. 2015).
The study of more elusive forms of shock phenomena, including fast-mode and intermediate (Alfvén) shocks, has started to become more prevalent in recent years. Utilizing the theoretical work of Montgomery (1959) and Hollweg et al. (1982), observational evidence for the non-linear steepening of elliptically polarized Alfvén waves (in the form of resonantly coupled fast-mode shocks) has recently been put forward by Grant et al. (2018). On the contrary, purely incompressible Alfvén waves are much more resistant to energy dissipation, and thus the observational signatures associated with their thermalization remains unclear. Early studies focusing on the energy dissipation of Alfvén waves employed 1.5D models to study coronal energetics and the subsequent driving of the solar wind (Hollweg 1981(Hollweg , 1992. More recently, Matsumoto & Suzuki (2014) concluded that shock heating, arising from a photospheric Alfvénic driver, was likely the dominant mechanism in the chromosphere. Arber et al. (2016) utilized 1.5D numerical models to show that Pedersen resistivity is able to directly dissipate highfrequency Alfvén waves, while Snow et al. (2018) revealed theoretical evidence for how vortex motion applied to magnetic flux tubes is able to drive intermediate shocks that propagate upward with speeds of approximately 50 km s −1 , hence transporting energy and momentum into the upper layers of the solar atmosphere. Recently, Snow & Hillier (2019) demonstrated how long-lived intermediate shocks can form within the confines of a traditional slow-mode shock, with their extended lifetimes arising due to the collisional coupling between species in a partially ionized plasma like the solar chromosphere. Hence, there has been a rapid improvement in our theoretical understanding of intermediate shocks, which naturally inspires the search for these signatures in cutting-edge observational image sequences.
Here, we present the first observational detection of intermediate shocks manifesting at the chromospheric umbral/penumbral interface of a sunspot. We use Ca II 8542Å spectropolarimetric data products obtained with the Dunn Solar Telescope (DST), in conjunction with modern inversion techniques and analytical theory, to provide unique insights into the dynamic plasma fluctuations associated with the manifestation of intermediate shock fronts in the Sun's magnetic atmosphere.

OBSERVATIONS
The data presented in this study represents an observational sequence carried out during 13:39 -16:43 UT on 2016 May 20, using the National Solar Observatory's DST at Sacramento Peak, New Mexico, USA. The telescope was pointed at the active region NOAA 12546, which is one of the largest sunspots to emerge on the solar surface in the past 20 years. The sunspot was positioned very close to disk center at the time of observing, at heliocentric coordinates (33 , −83 ), corresponding to a heliocentric angle of 5.37 • (µ 0.997), or S07.0W02.0 in the conventional heliographic coordinate system. Observations were obtained with the Interferometric BIdimensional Spectrometer (IBIS ;Cavallini 2006;Reardon & Cavallini 2008).
The IBIS instrument was utilized to obtain a long time series of high-spatial and -temporal resolution spectropolarimetric imaging scans of the photospheric Fe I 6173Å and chromospheric Ca II 8542Å spectral lines. Twenty-one discrete, equidistant wavelength steps were used across each of the Fe I 6173Å and Ca II 8542Å lines, with the Fe I 6173Å line covering the range 6173.14 − 6173.54Å using a spectral sampling of 20 mÅ, while the Ca II 8542Å line covered a range 8541.50 − 8542.70Å with a sampling of 60 mÅ. Sampling the full-Stokes profiles, using an exposure time of 80 ms, of AR12546 across the Fe I 6173Å and Ca II 8542Å spectral lines leads to a total cadence of 48 s, with a spatial sampling of 0 . 098 per pixel. The same data set was employed in Stangalini et al. (2018) who analysed circular polarisation oscillations to detect propagating MHD surface modes within the sunspot. To learn more about the magnetic structure in this data set we direct the reader to the paper by Murabito et al. (2019), who performed a complete analysis of the magnetic field geometry by using spectropolarimetric inversions.
The long duration of the observing period (∼3 hours), coupled with the relatively short temporal cadence, makes the data ideal for studying oscillatory phenomena in the vicinity of the active region. The near-simultaneous observations of both the photosphere and the chromosphere makes it possible to search for signatures of wave propagation through the atmosphere, with the full-Stokes information allowing the application of inversion routines to locations of interest. A whitelight camera, synchronised with the narrow-band feed, was employed to enable processing of the narrowband image sequences. High order adaptive optics (Rimmele 2004) were employed throughout the data acquisition, with the large cen- Figure 1. Left: IBIS Ca II red wing image of active region NOAA 12546, acquired at 8542.6Å (line core +0.5Å). Middle-left: Ca II 8542Å line core image, where the over-plotted colored pixels represent the locations where reversals in either Stokes Q/Ic , U/Ic or V /Ic spectropolarimetric profiles are detected between neighboring frames. The color scale indicates the total number of polarimetric reversals occurring in that pixel throughout the duration of the time series. Middle-right: Locations and occurrences of pixels that exhibit a change in their Stokes I/Ic magnitude exceeding 3σ above their quiescent value. Right: The red lines represent plasma-β = 1 isocontours spanning optical depths of (∼450 km) −3 > log 10 τ > −4 (∼625 km). The blue contour represents an outer boundary that encompasses all locations that are deemed 'active' pixels (see Section 3.1). tral sunspot chosen as the lock-point. While data reduction of the observations followed standard calibration techniques (i.e., dark subtraction, flat fielding and polarimetric calibration), the images obtained were also subjected to Multi-Object Multi-Frame Blind Deconvolution (MOMFBD;van Noort et al. 2005) techniques in order to mitigate the effects of atmospheric aberrations. Simultaneous broad band images were restored and narrow band images were destretched using them as a reference.
A contextual full-disk continuum image was obtained from the Helioseismic and Magnetic Imager (HMI; Schou et al. 2012) on board the Solar Dynamics Observatory (SDO; Pesnell et al. 2012) at 13:36 UT, which was utilized for the purpose of co-aligning the IBIS images with the HMI reference frame. A sub-field of 400 × 400 was extracted from the full-disk image, with a central pointing close to that of the ground-based observations. The HMI continuum image was then used to define absolute solar coordinates, with the IBIS observations subsequently subjected to crosscorrelation techniques to provide sub-pixel co-alignment accuracy. The composition and pointing of fully-calibrated IBIS images is displayed in Figure 1 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, we hence 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 hashed 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 dynam- Figure 2. Left panels represent sample Ca II 8542Å quiescent (i.e., pre-shock; solid black line) spectropolarimetric profiles for Stokes I/Ic , Q/Ic , U/Ic, and V /Ic. The right panels represent the the corresponding Stokes profiles associated with shocked plasma. The red dashed lines in each panel show the synthetic profiles generated from the NICOLE inversions. The blue shaded regions represent the spatially and temporally averaged standard deviations between the input IBIS and synthesized intensities across all pixels employed in the analysis.
ically 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; and 2. There needs to be a distinct increase, ∆I, in the integrated Stokes I/I c intensity originating from within the pixel location, with ∆I > 3σ set as a threshold to distinguish quiescent pixels from their active counterparts, where σ is the standard deviation of the in- Employing these criteria ensured that all 'active' pixel detections were statistically significant and not a result of smallamplitude waves or instrumental noise. Pixels that satisfied the individual criteria mentioned above were then cospatially and co-temporally 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 (>3σ) 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.

Inversions
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 datasets, solving multi-level, 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; ∆λ = 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 con-sists 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 (log 10 τ ). 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 active pixels identified in our dataset 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, B x and 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 wavelengthdependent 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). 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 < log 10 τ < −2.0, is shown using colored shaded regions in Figure 3.

RESULTS
Following the completion of the 6964 non-LTE spectropolarimetric inversions using NICOLE, we are provided with a number of key plasma parameters as a function of optical depth. These include the vector magnetic fields (B x , B y and B z ), temperatures, LOS velocities and densities for the pixels of interest, both during and immediately prior to their 'active' stage. Below we discuss the relationships between these constituent components of the plasma parameter space.

Velocity and Temperature Changes
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 display 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 log 10 τ = −2 corresponds to the photosphere, while log 10 τ = −6 is indicative of upper-chromospheric locations. We note that the atmospheric solution contains larger uncertainties in the log 10 τ = −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 red-shifted (i.e., downflowing) plasma, and negative values correspond to blue-shifted (i.e., upflowing) plasma.
Inspection of the results shows that the plasma is predominantly red-shifted immediately prior to the formation of a shock, but this changes abruptly to blue-shifted 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) non-linearities. 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 log 10 τ = −6 and −5, it has been seen that 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 log 10 τ = −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, ∆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 K < ∆T < 2500 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 log 10 τ = −5 (purple data points), where quies- The shock LOS Doppler velocities plotted as a function of their quiescent (i.e., pre-shock) Doppler velocities for the same pixel location. The black dashed lines are located along velocities of 0 km s −1 to provide easier visual segregation of the directional characteristics of the bulk motions. The background blue-red color scheme helps visualize the Doppler velocities corresponding to each quadrant of the plot, with progressively more blue and red colors representing larger up-and down-flowing material, respectively. Right: Shock temperature changes displayed as a function of the pre-shock background temperature. The dashed black line represents a zero change in temperature (i.e., ∆T = 0 K). The background blue-red color scheme provides a visual representation of temperature, with more red colors corresponding to both hotter quiescent and shock-induced temperatures. In both panels the colored data points correspond to the optical depths at which the plasma parameters are extracted, as defined in the legends located in the upper-left corner of each panel.
cent plasma with temperatures of around 3000 K experience ∆T ∼ 1500 K, while background temperatures of 5000 K only provide ∆T ∼ 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 ∆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 ∆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 K, ∼3650 K, ∼3465 K, ∼6020 K and ∼10845 K for optical depths corresponding to log 10 τ = −2 (∼250 km; Maltby et al. 1986), log 10 τ = −3 (∼450 km), log 10 τ = −4 (∼625 km), log 10 τ = −5 (∼1150 km) and log 10 τ = −6 (∼1850 km), respectively. Interestingly, the largest average ∆T (+19.6%) occurs at an optical depth of log 10 τ = −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 < log 10 τ < −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., 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.

Magnetic Field Perturbations
In a similar manner to how the temperature fluctuations are depicted in the right panel of Figure 4, Figure 5 displays the shock vector magnetic fields as a function of their preshock values, with the transverse (B x and B y ) and vertical (B z ) components displayed in the upper, middle and lower 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 eastwest 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 lines of best fit (dotted red lines in the upper and middle panels of Figure 5) being close to −1, with gradients of −0.99 and −1.00 found for the B x and B y fits, respectively. To make Figure 5. Two dimensional density scatter diagrams showing the vector components (Bx, upper; By, middle; Bz, lower) of the shock magnetic fields as a function of their quiescent (i.e., preshock) values. In each panel the vertical and horizontal dashed black lines highlight magnetic field components equal to 0 Gauss. The shade of each hexagon represents the density of points within that region. For the Bx (top) and By (middle) panels, the dotted red line highlights the linear line of best fit, with the shaded red region (bounded by small dotted red lines) indicating the 1σ errors associated with each fit. In the Bz (lower) panel, the dotted red line shows a 1:1 slope, where data points lying on this line have identical Bz magnitudes in both the quiescent and shock stages. use of the large number statistics at our disposal, the relationships illustrated in Figure 5 for B x and 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 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 anti-correlation found between pre-and post-shock transverse magnetic field fluctuations.
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 B x plot in the top panel of Figure 5, when compared with the 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 B x magnitudes will be strongest due 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 B x (−0.99; CI [−1.00, −0.98]) and 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 B z component of the magnetic field. As displayed in the lower panel of Figure 5, the signs and magnitudes of the 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 lower panel of Figure 5. It might be natural to assume that a spectropolarimetric reversal in Stokes V /I c (see, e.g., the lower-right panel in Figure 2) would indicate a polarity reversal in that particular pixel location. However, spectropolarimetric reversals have been witnessed previously by de la Cruz , with such signatures not necessarily implying a physical reversal of the magnetic field polarities, as is also implied in the lower 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 Figure 6. The fluctuations in Bz (i.e., ∆Bz) arising from the development of a shock, plotted as a function of the change in the total magnetic field strength (∆Btot) also produced from the commencement of the shock. The shade of each hexagon represents the density of points within that region. The vertical and horizontal dashed black lines highlight magnetic field fluctuations equal to 0 Gauss. The dotted red line displays the linear line of best fit, with the shaded red region (bounded by small dotted red lines) indicating the 1σ errors associated with the fit. 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 B x and B y values flip, while the B z component of the magnetic field remains approximately constant with the same sign. Due to the very pronounced reversals in the B x and B y components (i.e., line of best fit gradients very close to −1 in the upper and middle panels of Figure 5), fluctuations in the B x and B y terms should only contribute to very minor changes in the total magnetic field strength, B tot . This implies that any changes in B z should produce a near-equivalent fluctuation in B tot -i.e., ∆B z = ∆B tot . Figure 6 displays a scatterplot detailing the relationship between ∆B z as a function of ∆B tot , where the dotted red line highlights a linear line of best fit between the two variables. The gradient associated with the line of best fit is 0.92 ± 0.04, indicating a very close correlation between fluctuations in the vertical component of the magnetic field (B z ) and the total overall magnetic field strength (B tot ).
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, +B x /+B y or −B x /−B y . Hence, care needs to be taken when examining the outputs of NICOLE spectropolarimetric inversions as a flip (e.g., +B x → −B x and +B y → −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., +B x remains +B x and +B y remains +B y ). This may be the cause of some small positive correlations seen in Figure 5, particularly at relatively weak magnetic field strengths (|B x/y | 500 G).
While the selection criteria (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.   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  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 dataset, 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 Figure 7. Histograms documenting the percentage changes in the plasma density that are caused by shock formation for optical depths corresponding to log 10 τ = −3 (∼450 km; lower right), log 10 τ = −4 (∼625 km; upper right), log 10 τ = −5 (∼1150 km; lower left), and log 10 τ = −6 (∼1850 km; upper left). Positive values (i.e., >0%) are representative of shock-induced density enhancements, while negative values indicate a reduction in the local plasma density following the formation of a shock. 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 criteria 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.

Density Ratios
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 log 10 τ = −3 and log 10 τ = −4 (high photosphere and low chromosphere, respectively) there are shock formation signatures, with increases 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 log 10 τ = −5, while at the extreme upper-boundary of the chromosphere (log 10 τ = −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 log 10 τ ∼ −6. However, this area of research is in its infancy, and requires dedicated shock-capturing numerical simulations, (e.g. using the Lagrangian-Eulerian Remap (LareXd) code; Arber et al. 2001 or MPI-AMRVAC code;Porth et al. 2014) to further test the effects of such phenomena.

SHOCK CLASSIFICATION
For all isolated pixels of interest, we have inversion outputs that provide us with the specific plasma conditions both before and after the manifestation of a shock. The deduced trends for active pixels (see, e.g., Figures 5 & 6) indicate a reversal of their transverse magnetic fields from quiescent to shocked states. From theory, this can be interpreted as either a rotational discontinuity or an MHD shock (Goedbloed et al. 2010). Rotational discontinuities require conservation of the total magnetic field (i.e., ∆B tot = 0 Gauss). However, from examination of Figure 6 it is clear to see that changes in the total magnetic field are commonly experienced, where ∆B tot ∼ ∆B z , suggesting that the active pixels are not best described by rotational discontinuities.
For an MHD shock to be a viable interpretation for the captured plasma dynamics, there must be evidence for shockinduced compression of the local plasma. Indeed, examination of Figure 7 clearly shows that active pixels demonstrate clear density, ρ, increases at their point of formation, i.e., ρa /ρ b > 1, where the subscripts a and b represent the shocked ('after') and quiescent ('before') stages of the shock evolution. Furthermore, the LOS Doppler velocities, v los , displayed in Figure 4 also demonstrate a clear discontinuity, whereby v los,b − v los,a = 0. As the thermodynamic properties do not depend on the frame of reference (Goedbloed et al. 2019), we first check whether entropy has increased across the shock domain. The entropy change, ∆S, from the quiescent to shocked states is evaluated at four discrete optical depths (log 10 τ = −3, −4, −5, −6) following, where p is the plasma pressure (p = ρk B T /m p µ), k B is the Boltzmann constant, T is the temperature, m p is the mass of a proton, µ is the mean molecular weight (µ = 0.5), and γ is the adiabatic index. Here, an adiabatic index of γ = 1.12 ± 0.01 is utilized, which is consistent with the spectropolarimetric investigation of another chromospheric sunspot by Houston et al. (2018). The mean entropy values for all 6964 pixels are displayed in the lower-right panel of Figure 9 as a function of optical depth. At optical depths of log 10 τ = −3 and log 10 τ = −4, which are consistent with the shock formation heights, we see a clear increase in entropy (i.e., ∆S 0). At optical depths corresponding to higher geometric heights (i.e., log 10 τ < −5), the entropy change is less severe, with the 1σ error bars associated with the upper extremity of the chromosphere (log 10 τ = −6) encompassing ∆S = 0. This implies that the shock formation represents a localized change in the MHD plasma quantities, with the biggest fluctuations experienced close to the formation height of the shock itself (i.e., −3 > log 10 τ > −5). To verify this interpretation and further classify the type of MHD shock present, it is necessary to employ the Rankine-Hugoniot (RH) relations.

Rankine-Hugoniot Classification
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 is 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, where µ 0 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;Cally & Goossens 2008). Recently, their importance has also been demonstrated for the generation of resonantly driven fastmode shocks towards 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 > log 10 τ > −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, To estimate the shock normal direction,n s , we quantify the angles which 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 ofn s on the horizontal plane, Then using the first RH condition, which states that the normal magnetic field component, B n , remains constant, i.e., B n,a = B n,b , we calculate the angle ϑ s betweenê z (the vertical) andn s , Once we haven s , we can decompose the magnetic field, B, into its normal, B n , and tangential, B t = B − B nns , 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, B, and the shock normal,n s , in both quiescent and shocked states through the relation, 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 A , θ), 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 solutions output for this shock strength. The narrow band 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.
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 This relationship states that ρv n is identical in pre-and postshock states. As solutions to the shock adiabatic relation, the corresponding Alfvén mach number (M A = v n √ ρ/B 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, 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, 3. We should expect that the observed difference in the LOS velocity, ∆v los , correlates with the observed jump in the normal velocity, v n . Hence, we check the value of the quantity, ∆v los = |v n,b − v n,a − v los,b + v los,a | . 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 ≤ log 10 τ ≤ −3). From the three RH conditions defined above, we obtain 5 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 (log 10 τ = −3, −4; where the entropy change, ∆S, is largest; lower-right panel of Figure 9). The offsets associated with the fast solution, at the optical depths of peak shock formation, are an orderof-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, suggests that the detected shock behavior is best classified by Alfvén (or intermediate) MHD shocks. At higher geometric heights (log 10 τ = −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 upwards 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., log 10 τ = −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 (log 10 τ = −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 impor-tant, 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).

CONCLUSIONS
In this paper, we have presented high temporal resolution spectropolarimetric Ca II 8542Å observations, captured by the IBIS instrument at the Dunn Solar Telescope. Through comprehensive analysis of a large sunspot near solar disk center, combined with advanced inversion techniques, the plasma evolution during the formation of shocks demonstrating spectropolarimetric reversals are investigated. We find significant changes in the temperatures, LOS velocities, densities, and vector magnetic fields between the quiescent locations and those demonstrating shock activity. The largest fluctuations occur at optical depths of log 10 τ = −4, which is consistent with a geometric height of approximately 625 km, close to the boundary between the upper photosphere and the lower chromosphere. This layer has also played host to recent observations of resonantly-amplified fast-mode shocks (Grant et al. 2018), in addition to a plethora of slow magnetoacoustic umbral shock phenomena (de la Cruz Henriques et al. 2017;Joshi & de la Cruz Rodríguez 2018).
Through examination of the associated entropy and total magnetic field strength changes, we are able to exclude rotational discontinuities as a possible explanation of the observed dynamics, instead suggesting the presence of a developing magnetohydrodynamic shock. Rankine-Hugoniot relations are then used to classify the shock type, specifically by decomposing the shock characteristics into the reference frame of the shock itself, allowing the induced normal velocities to be compared to the plasma parameters derived from modern spectropolarimetric inversion routines. Through minimization of the differences between the observationally derived characteristics and those associated with theoretical Rankine-Hugoniot relationships, we find firsttime evidence highlighting the presence of Alfvén (intermediate) shocks manifesting close to the perimeter of a sunspot umbra. The importance of finding such shock activity cannot be underestimated. Now that the manifestation of Alfvén shocks has become apparent in magnetic sunspot structures, their existence may also be important for supplying thermal energy to the atmosphere of other magnetic features, including pores, magnetic flux ropes, plumes, and magnetic bright points.
Future work will require the use of all Rankine-Hugoniot relations to make more definitive shock classifications. To do so requires the determination of the shock-induced tangential velocity fields, which to observe requires greater temporal and spatial resolutions, in conjunction with high-precision spectropolarimetry. IBIS requires a scan time on the order of 23 seconds to attain sufficient spectropolarimetric signals to allow for accurate inversion processes. Felipe et al. (2018) showed that inversions of profiles scanned from from blue to red with similar cadences to IBIS, may not accurately reproduce the physical magnetic field during the development of a shock event. Future instruments, including fibre-fed spectropolarimeters, will enable simultaneous spatial and spectral information with reduced cadences, resulting in more accurately constrained parameters from inversion processes. Furthermore, with modern numerical simulations investigating the role of partial ionization during shock development (Snow & Hillier 2019), attention will naturally turn to an observational study incorporating neutral and ionized species of the same element (e.g., Ca I & II) that, when combined with cutting-edge inversion techniques, will allow the ionization degree to be derived as a function of height in the solar atmosphere. Such studies involving partially ionized plasma will still uphold MHD Rankine-Hugoniot conditions across the shock features, but may differ greatly in the details of the actual shock variations, which are treated as discontinuous in an ideal MHD setting. With the imminent arrival of a swathe of new observing facilities, such as Daniel K. Inouye Solar Telescope (DKIST; Keil et al. 2004), the Indian National Large Solar Telescope (NLST; Hasan et al. 2010), the European Solar Telescope (EST; Collados et al. 2013), and Solar-C (Watanabe 2014), we believe that future studies will be able to unequivocally document the role of MHD shocks in supplying thermal energy to the solar atmosphere.