Shock Heating Energy of Umbral Flashes Measured with Integral Field Unit Spectroscopy

TL;DR

This study measures the shock heating energy of umbral flashes in sunspot umbrae using novel integral field spectroscopy, revealing that shock energy is insufficient to sustain the chromosphere and that shocks propagate at sound speed.

Contribution

It introduces a new integral field spectrograph for detailed shock analysis in sunspot umbrae and quantifies shock properties and energies using radiative transfer modeling.

Findings

Shock heating energy per unit mass is 2 x 10^{10} erg g^{-1}.

Shocks propagate upward at the sound speed.

Mach number does not depend on upstream temperature.

Abstract

Umbral flashes are periodic brightness increases routinely observed in the core of chromospheric lines within sunspot umbrae and are attributed to propagating shock fronts. In this work we quantify the shock heating energy of these umbral flashes using observations in the near infrared He I triplet obtained on 2014 December 7 with the SpectroPolarimetric Imager for the Energetic Sun (SPIES), which is a novel integral field unit spectrograph at the Dunn Solar Telescope. We determine the shock properties (the Mach number and the propagation speed) by fitting the measured He I spectral profiles with a theoretical radiative transfer model consisting of two constant property atmospheric slabs whose temperatures and macroscopic velocities are constrained by the Rankine-Hugoniot relations. From the Mach number, the shock heating energy per unit mass of plasma is derived to be 2 x 10^{10} erg g^{-1}, which is insufficient to maintain the umbral chromosphere. In addition, we find that the shocks propagate upward with the sound speed and the Mach number does not depend on the temperature upstream of the shocks. The latter may imply suppression of the amplification of the Mach number due to energy loss of the shocks.

Equations

$$q=\frac{1}{\gamma-1}\frac{R}{\mu}\,T_{{\rm u}}\frac{\rho_{{\rm u}}}{\rho_{{\rm d}}}\Bigl{[}\frac{p_{{\rm d}}}{p_{{\rm u}}}-\Bigl{(}\frac{\rho_{{\rm d}}}{\rho_{{\rm u}}}\Bigr{)}^{\gamma}\Bigr{]},$$

$$\frac{{\rho }_{\mathrm{d}}}{{\rho }_{\mathrm{u}}}=\frac{V_{\mathrm{u}}}{V_{\mathrm{d}}}=\frac{(\gamma +1)M^2}{(\gamma -1)M^2+2},$$

$$\frac{p_{\mathrm{d}}}{p_{\mathrm{u}}}=\frac{2\gamma M^2-\gamma +1}{\gamma +1},$$

$$\frac{T_{\mathrm{d}}}{T_{\mathrm{u}}}=\frac{[2\gamma M^2-\gamma +1][(\gamma -1)M^2+2]}{(\gamma +1{)}^2{M}^2},$$

$$q=\frac{\gamma}{\gamma-1}\frac{R}{\mu}\,T_{{\rm u}}\frac{\rho_{{\rm u}}}{\rho_{{\rm d}}}\Bigl{[}\Bigl{(}\frac{p_{{\rm d}}}{p_{{\rm u}}}\Bigr{)}^{\frac{1}{\gamma}}-\frac{\rho_{{\rm d}}}{\rho_{{\rm u}}}\Bigr{]}.$$

$$I(\lambda )=S(1-e^{-{\tau }_{\lambda }})+I_{\mathrm{C}}{e}^{-{\tau }_{\lambda }},$$

$$\tau_{\lambda}=\tau\sum_{i=1}^{3}\omega_{i}\exp{\Bigl{\{}-\Bigl{(}\frac{\lambda-\lambda_{i}-\lambda_{{\rm D}}}{\Delta\lambda}\Bigr{)}^{2}\Bigr{\}}},$$

$$I(\lambda )=S_{\mathrm{u}}(1-e^{-{\tau }_{\mathrm{u}}(\lambda )})$$

$$\tau_{{\rm u}}(\lambda)=\tau_{{\rm u}0}\sum_{i=1}^{3}\omega_{i}\exp{\Bigl{\{}-\Bigl{(}\frac{\lambda-\lambda_{i}-\lambda_{{\rm Du}}}{\Delta\lambda_{{\rm u}}}\Bigr{)}^{2}\Bigr{\}}},$$

$$\tau_{{\rm d}}(\lambda)=\tau_{{\rm d}0}\sum_{i=1}^{3}\omega_{i}\exp{\Bigl{\{}-\Bigl{(}\frac{\lambda-\lambda_{i}-\lambda_{{\rm Dd}}}{\Delta\lambda_{{\rm d}}}\Bigr{)}^{2}\Bigr{\}}},$$

$$\displaystyle\Delta\lambda_{{\rm d}}=\Biggl{[}\Delta$$

$$V_{{\rm u}}=M\biggl{\{}\frac{\gamma Rmc^{2}}{2\mu k_{{\rm B}}\lambda_{0}^{2}}(\Delta\lambda_{{\rm u}}^{2}-\Delta\lambda_{{\rm I}}^{2})\biggr{\}}^{0.5},$$

$$V_{\mathrm{d}}=V_{\mathrm{u}}\frac{(\gamma -1)M^2+2}{(\gamma +1)M^2},$$

$${\lambda }_{\mathrm{Du}}={\lambda }_0(V_{\mathrm{u}}+U)/c,$$

$${\lambda }_{\mathrm{Dd}}={\lambda }_0(V_{\mathrm{d}}+U)/c.$$

$$U=\frac{{\lambda }_{\mathrm{Du}}}{{\lambda }_0}c-M\sqrt{\frac{\gamma R{T}_{\mathrm{u}}}{\mu }}.$$

$$U\sim \frac{{\lambda }_{\mathrm{D}}}{{\lambda }_0}c-M\sqrt{\frac{\gamma R{T}_{\mathrm{u}}}{\mu }}.$$

$$M\approx 1+\frac{v_{0}}{C_{s}}{\rm exp}\Bigl{(}\frac{2h}{5\Lambda}\Bigr{)},$$

Full text

line: profiles — shock waves — Sun: chromosphere — Sun: oscillations — (Sun:) sunspots — techniques: imaging spectroscopy

††facilities: DST(SPIES)††software: Solar Soft

1 Introduction

In the solar chromosphere, the temperature increases with height and the radiative emission is larger than that expected from radiative equilibrium (Vernazza et al., 1981; Maltby et al., 1986). Many kinds of heating mechanisms have been suggested to balance the large excess of radiative cooling energy, for example, shock heating (Schwarzschild, 1948; Carlsson & Stein, 1997; Beck et al., 2008), magnetic energy dissipation in current sheets (Parker, 1983; Solanki et al., 2003; Socas-Navarro, 2005; Tritschler et al., 2008), viscous and Ohmic dissipations in vortical magnetic structures (Moll et al., 2012), Alfv${\rm\acute{e}}$n wave turbulence (van Ballegooijen et al., 2011), and dissipation resulting from collisions between magnetized ions and unmagnetized neutral atoms (Osterbrock, 1961; Khomenko & Collados, 2012). Among the challenges to distinguish which heating mechanisms dominate is establishing remote sensing techniques to estimate each mechanism’s heating energy.

Evidence that shocks play a role in chromospheric heating comes from umbral flashes. The umbral flashes are a ubiquitous feature of the dynamic chromosphere above sunspot umbrae. They manifest as periodic brightness increases in the core of chromospheric spectral lines with a period of $\sim 3$ minutes (Beckers & Tallant, 1969; Wittmann, 1969; Khomenko & Collados, 2015). Temporal Doppler velocity fluctuations are correlated with the umbral flash (Kneer et al., 1981). Using spectroscopic observations, Giovanelli et al. (1978) interpreted the velocity fluctuations as upward propagating waves, and Lites (1984) found that the propagating waves develop into shock waves. The explanation of the umbral flash as upward propagating shocks have been confirmed by some observations (e.g. Kentischer & Mattig, 1995; Yoon et al., 1995; Brynildsen et al., 1999, 2003). Lites (1984) also suggested the possibility that umbral shocks contribute to the heating of the umbral chromosphere.

Some observations have provided evidence for temperature enhancements associated with the umbral flashes (e.g. Shibasaki, 2001; Iwai et al., 2017). de la Cruz Rodríguez et al. (2013) derived an enhancement of 1000 K from chromospheric observations in Ca II 8542 Å obtained using a narrow-band Fabry-P${\rm\acute{e}}$rot based imaging spectrograph. Using similar observations, Joshi & de la Cruz Rodríguez (2018) reported that the temperature at the chromosphere increases from 3,700 K up to 6,200 K. Grant et al. (2018) also measured significant temperature enhancements up to a maximum of $\sim 20$ % through Ca II 8542 Å observations obtained by another filter-based spectropolarimeter. Using the He I 10830 Å triplet, Houston et al. (2018) obtained temperature fluctuations of $\pm 10$ % from observations with a grating-based spectrograph.

The heating energy flux required to maintain the umbral chromosphere has been estimated to be $2.6\times 10^{6}$ erg cm*-2* s*-1* based on the net radiative cooling rate for a semi-empirical sunspot model (Avrett, 1981). To evaluate the role of shock heating, some authors have investigated energy transported by acoustic waves. Upward propagating acoustic waves will eventually develop shocks due to the decrease of the density with height, and consequently, their energy will dissipate via shock heating (Schwarzschild, 1948). Felipe et al. (2011) reproduced observed chromospheric wave signatures through a three-dimensional magnetohydrodynamic numerical simulation perturbed at the photosphere by the observed photospheric velocity fluctuations, and they obtained the acoustic energy flux as a function of the height. However, the average value of the derived acoustic energy was insufficient to sustain the umbral atmosphere. Even in the upper photosphere, low acoustic energy fluxes were estimated from Doppler velocity measurements (Giovanelli et al., 1978; Kneer et al., 1981; Chae et al., 2017).

In contrast to the above measurements, some recent observations have suggested sufficient enough acoustic energy damping in the chromosphere. Kanoh et al. (2016) obtained an upward acoustic energy flux of $2.0\times 10^{7}$ erg cm*-2* s*-1* in the photosphere and $8.3\times 10^{4}$ erg cm*-2* s*-1* within the lower transition region from Doppler velocity measurements and density estimations obtained using the Hinode (Kosugi et al., 2007) and the Interface Region Imaging Spectrograph (IRIS, De Pontieu et al., 2014). In agreement with Kanoh et al. (2016), Krishna Prasad et al. (2017) reported a decrease with height of the energy flux at multiple atmospheric heights from intensity oscillations in multiple spectral lines. Furthermore, Grant et al. (2015) interpreted observed intensity and area fluctuations as upward propagating sausage-mode oscillations and estimated their energy flux, which decreases with height similar to Krishna Prasad et al. (2017).

It is currently not clear if the observed reduction of acoustic flux with height implies wave dissipation and heating via shocks, or if it represents loss of wave energy to some other form, such as distortions of the magnetic field, energy transfer to some other wave modes that have less compression, or cascade to unresolved spatial or temporal scales (e.g. Cally & Goossens, 2008; Reardon et al., 2008). In order to assess the ability of shocks to heat the umbral chromosphere, we propose to model the thermodynamic and radiative properties of the observed shock plasma. Our approach to calculate shock heating rates uses spectroscopic observation, while Lee & Yun (1985) uses empirical atmospheric models.

Here, we determine measurements of the shock heating energy rate using spectroscopic observations of the He I 10830 Å  triplet obtained with a fiber-optic-based Integral Field Unit (IFU). Our analysis allows us to calculate the contribution of the shocks to the umbral chromospheric heating. To account for the fact that dynamic flashes rapidly change the shape of the spectral line profiles, we take advantage of IFU-fed diffraction-grating-based spectrograph that enables us to measure spectral profiles in one exposure, i.e. without scanning in wavelength as done by Fabry-Perot based instruments, or scanning in space as done by traditional slit-spectrographs. The IFU also allows us to observe the entire sunspot efficiently with a short enough cadence to resolve the shock dynamics, and it enables us to compare the spatial distribution of shock properties with that of heating signatures. This IFU technique for such high cadence imaging spectroscopy represents a new class of solar instruments based on multiplexed diffraction-grating-based spectrographs (Jaeggli et al., 2010; Lin, 2014; Schad & Lin, 2017; Jurčák et al., 2019). Below we describe details regarding the observations (Section 2), the determination of the Mach number and the temperature upstream of the shock (Section 3), and the results of the shock heating energy rate derived from the shock characteristics (Section 4). We then discuss the contributions of the shocks to the umbral chromospheric heating (Section 5) and finally our conclusions (Section 6).

2 Observations and data processing

The leading sunspot in NOAA active region 12227 was observed in a near infrared spectral window containing the He I 10830 Å triplet with the SpectroPolarimetric Imager for the Energetic Sun (SPIES, Lin & Versteegh, 2006; Lin, 2012; Schad et al., 2014) on the NSO’s 76 cm Dunn Solar Telescope (DST, Dunn, 1969) between 22:37 and 23:12 UT on 2014 December 7. The DST high-order adaptive optics system achieved real-time seeing correction and image stabilization (Rimmele et al., 2004). The He I triplet, which is formed in the upper chromosphere, provides a clear signature of shock waves via the temporal change of the Doppler shift (Lites, 1986; Centeno et al., 2006).

SPIES is a prototype instrument for the Diffraction Limited Near Infrared SpectroPolarimeter (DL-NIRSP, https://www.nso.edu/

telescopes/dkist/instruments/dl-nirsp/) of the National Science Fundation’s Daniel K. Inouye Solar Telescope (DKIST, Rimmele & ATST Team, 2008; Rimmele et al., 2015; Tritschler et al., 2016). Both SPIES and DL-NIRSP employ a fiber-optic-based integral field unit to obtain dispersed spectra within a contiguous two dimensional spatial field of view in a single exposure. SPIES’s IFU (Schad et al., 2014) contains 15,360 birefringent rectangular fibers, and reconfigures a near-square array at a focus of the DST into four slit arrays at the entrance of the near Littrow spectrograph, which are dispersed and imaged simultaneously on a $2048\times 2048$ infrared detector. Each individual fiber core has an aspect ratio of 4:1, and thus four adjacent fiber cores are used to form one square spatial sampling pixel. SPIES therefore allows us to observe a field of 64 $\times$ 60 spatial pixels, corresponding to $15\times 16$ arcsec2 simultaneously with a high spatial resolution of 0.48 arcsec and a high spectral resolution of 48 mÅ. The instrumental line width is derived as $29$ mÅ  from comparisons of observed telluric line profiles and those of a very high resolution solar flux atlas (Kurucz et al., 1984). In addition, a larger field of view can be mosaiced by steering the optical field with a two-axis field steering mirror. In our observations below, we perform a $4\times 4$ mosaic for a total field of view of $59\times 65$ arcsec2 and a temporal cadence of 14 seconds.

After the subtraction by a dark frame, the SPIES dispersed spectral images were reduced via division by a flat field made as follows. Frames for the flat field were acquired by observing quiet regions around disk center while the telescope conducted a random small amplitude scanning pattern to smear out structures. The flat field data were subtracted by a dark frame and averaged. Next, the averaged frame was divided by another flat field obtained with a lamp, which we inserted in the optical light path to calibrate the transmission profile of a blocking filter. Moreover, the transmission profile was calibrated by comparing spectral profile in continuum with that of an atlas profile (Kurucz et al., 1984). Finally, we applied two-dimensional principal component analysis to the flat frame to remove the solar or telluric spectral features, in analogy to removal of fringes (Casini et al., 2012).

To facilitate the reconstruction of the spatial image of the IFU, the relative position of each pixel in the detector plane to the location in the IFU entrance plane is determined by placing an undersized square field stop in an upstream location conjugate to the IFU entrance plane. Orthogonal scans of this field stop using the field steering mirror can then be used to pinpoint the location of each pixel in the field-of-view (see, e.g. Lin & Versteegh, 2006). Figure 1 shows images of the sunspot reconstructed from the SPIES observations as well as that acquired by the Helioseismic and Magnetic Imager (HMI, Scherrer et al., 2012; Schou et al., 2012) onboard the Solar Dynamics Observatory (SDO). The sunspot has a light bridge (Bray & Loughhead, 1964; Lites et al., 2004; Katsukawa et al., 2007). Chromospheric features, for example superpenumbral fibrils (Schad et al., 2013, 2015), appear in the SPIES He I 10830 Å image. The IFU has some dead fibers and gaps between fibers where the instrument is not sensitive (Figure 1 d). In Figure 1 (b) and (c), the values at the places of the dead fibers and spaces are interpolated by the Delaunay triangulation. We do not use the interpolated values to derive any quantitative results in the sections below.

3 Diagnosis of shock properties

To obtain an observationally based estimate for the shock heating energy, one in principle needs to remotely sense the properties of the shocks themselves. The umbral flashes are clearly slow shocks (Centeno et al., 2006) and the magnetic pressure is larger than the gas pressure in the umbral atmosphere (Mathew et al., 2004). In that case, MHD jump conditions (ratios of the state variables) reduce to the hydrodynamic jump for parallel propagating shocks (Goedbloed et al., 2010).

We use two different models (Weymann, 1960; Schatzman, 1949) for the thermodynamic cycle of a parcel of chromospheric plasma that experiences a shock. In the model of Weymann (1960), the shocked plasma radiatively cools at constant specific volume until its entropy matches that of the un-shocked gas, after which it adiabatically expands to the un-shocked pressure and specific volume. Following Bray & Loughhead (1974) §6.5.3, the heat energy, $q$, dissipated by a hydrodynamic shock per unit mass of gas in this situation can be expressed in terms of the upstream and downstream properties of the plasma as

[TABLE]

where $\gamma$ is the polytropic index (the ratio of the specific heats), $R$ is the gas constant, $\mu$ is the mean molecular weight, $T$ is the temperature, $\rho$ is the density, $p$ is the pressure, and the u and d subscripts denote upstream and downstream components of the shock. Introducing the Mach number of the upstream component, $M\equiv V_{{\rm u}}/\sqrt{\gamma RT_{{\rm u}}/\mu}$, the hydrodynamic jump conditions are derived from the Rankine-Hugoniot relations as

[TABLE]

[TABLE]

[TABLE]

where $V_{{\rm u}}$ and $V_{{\rm d}}$ are velocities measured in the reference frame co-moving with the shock (Landau & Lifshitz, 1959). In the other model (Schatzman, 1949) where the gas is heated and compressed due to the shock, expands adiabatically, and radiates back to the original state, Bray & Loughhead (1974) give the $q$ as

[TABLE]

Substituting equations (2) and (3) into equations (1) and (5) respectively, we are able to derive the shock heating energy per unit mass of gas for either model based only on the Mach number and the temperature upstream of the shock.

In this section, we describe a method to determine the Mach number and the temperature from the observed spectral profiles in the He I 10830 Å triplet. First, the Doppler shift and line width were inferred by fitting the spectra with theoretical profiles computed with a radiative transfer equation using an atmospheric model based on a single constant property slab, i.e. a homogeneous plasma (Section 3.1). We demonstrate using this single slab model that there are phase relations between the measured Doppler shift and the line width. As such relations can be interpreted as contributions of the upstream and downstream components of a shock to a line profile, we then fit the profiles with a model consisting of two constant-property slabs for which the temperatures and the macroscopic velocities are constrained by the Mach number as in Equations (2) and (4) (Section 3.2).

Socas-Navarro et al. (2000b, a, 2001) and Centeno et al. (2005) also analyzed umbral spectra in chromospheric spectral lines using two-component models. Our method differs from these earlier approaches by taking into account restrictions of the Rankine-Hugoniot relations for state variables of the two components. In addition to the Mach number, we derive the shock heating energy to assess the ability of shocks to heat the umbral chromosphere.

3.1 Single Slab model

In order to infer temporal variations of the Doppler shift and the line width of the He I triplet, we first fit the spectral profiles obtained with SPIES considering a single component model; an example of the fitting results is shown in Figure 2. As the linear dispersion and wavelength reference point vary from fiber to fiber, we determine a wavelength calibration for each spectrum by fitting the telluric lines (10832.11 Å  and 10833.98 Å) with two Gaussian profiles. We fit the Si I 10827 Å profile with a Voigt function to reduce the effect of its line wing on the fitting for the He I 10830 Å triplet.

Intensity, $I$, of the theoretical profile of the He I triplet is obtained by solving a radiative transfer equation using an atmospheric model based on a constant property slab as,

[TABLE]

where $\lambda$ is the wavelength, $S$ is the source function of the slab, and $I_{{\rm C}}$ is the intensity in the continuum. Optical thickness of the He I triplet as a function of the wavelength, $\tau_{\lambda}$, can be written as

[TABLE]

where $\tau$ is the optical thickness of the slab at the spectral line center, $\omega_{i}$ are the statistical weights derived from the number of the magnetic sub-levels of the excited state under a constraint of $\sum_{i=1}^{3}\omega_{i}=1$, $\lambda_{i}$ are the center wavelengths of the triplet lines, $\lambda_{{\rm D}}$ is the Doppler shift, and $\Delta\lambda$ is the line width. Even though the constant property slab model is simple, it is a suitable model for reproducing spectral profiles of the He I 10830 Å line (Asensio Ramos et al., 2008).

Here we fix $I_{{\rm C}}$ to 1 as we normalize each intensity profile by $I_{{\rm C}}$ prior to the fitting. The fitting cannot determine the optical depth simultaneously with the source function in the optically thin case. Because the He I 10830 Å triplet is generally optically thin in the solar spectrum (Fontenla et al., 1993), we therefore fixed $\tau$ to be 0.5 (see figure 6 of Schad et al., 2015). The quality of the single component fits are satisfactory, although there are systematic residuals. The residuals at 10831.5 Å and around the telluric line at 10833.98 Å  are likely a result of the limited quality of the flat-field calibration.

The umbra surrounded by the black solid contours in Figure 1 encompasses roughly 2450 pixels. Each pixel within this area was fitted at each of the 150 time steps resulting in a total of $\sim 3.7\times 10^{5}$ spectral fits. The results are summarized in Figure 3 which shows the distribution of the fitted He I (a) source function, (b) Doppler shift, and (c) line width. We do not take into account the values at the dead fibers or spaces to make the histograms. Since the optical depth is fixed, the source function at the peak of the histogram is determined, and it is approximately equal to 0.5. Most ($\sim 99.8\%$) of the Doppler shifts and the line widths vary within the range $-0.2$ Å $<\lambda_{{\rm D}}<0.35$ Å and $0.2$ Å $<\Delta\lambda<0.65$ Å, respectively. Since the mode of the Doppler shifts is not equal to zero, the sunspot may be receding from the observer with a speed of $1.85$ km s*-1* due to the revolution of the Earth and the spins of the Earth and the Sun, although these contributions are not calibrated in our analysis.

Figure 4 shows temporal variations of the inferred Doppler velocity, $V_{{\rm Dop}}\equiv-c\lambda_{{\rm D}}/\lambda_{0}$, and the line width of the He I triplet for a part of the entire time series at a position with heliocentric coordinates of (S66*∘, E126∘) at 22:55UT, where $c$ is the light speed and $\lambda_{0}=10830$ Å is the wavelength at the line center. When $V_{{\rm Dop}}$ is larger than $-1.85$ km s-1*, the plasma may be traveling locally upward. Although the sunspot observed here has a penumbra, the amplitude of the Doppler shift is similar to that of the pore reported in Centeno et al. (2009). Figure 5 displays the spatial variations of the dominant oscillation periods in the Doppler shift and the line width. Three minute oscillations dominate the umbral chromosphere as many works have confirmed (e.g. Tziotziou et al., 2007; Socas-Navarro et al., 2009; Reznikova et al., 2012).

In order to examine phase relations between the Doppler velocity and the line width statistically, we construct two dimensional histograms of their temporal variations vs. the phase as displayed in Figure 6. We used values of the umbral data enclosed within the contour in the Figure 1 excluding those at the dead fibers or spaces. First, we extracted their temporal variations between the time of a minimum Doppler velocity and the next following minimum. We define the oscillation period as the time interval of the minimum Doppler velocities, and the phase as the time starting from the extracted temporal variation normalized by the period to 360*∘*. To show only phase relations, the measured amplitude is normalized by the peak-to-peak value of the extracted time series and is stacked to make the two dimensional histograms. Each histogram at each phase space bin is normalized by the total number.

The temporal variation in the $V_{{\rm Dop}}$ reveals a sawtooth pattern, that is, the phase of the maximum $V_{{\rm Dop}}$ is smaller than the half period. In addition, the line width has its maximum value during the time $V_{{\rm Dop}}$ is increasing. The sawtooth pattern and the enhancement of the line width during the $V_{{\rm Dop}}$ increases indicate that a shock wave was upwardly propagating through the formation layer of the He I triplet during this phase. We interpret the increase of the line width as the superposition of newly shocked rising plasma and falling plasma from the passage of the previous shock. Both components contribute to the formation of the line profile during the phase, as Tian et al. (2014) demonstrated for an IRIS observation in the Si IV line formed at the transition region.

Figure 7 shows schematic plots to explain this interpretation. Because the sunspot was located near disk center (Figure 1), the Doppler velocity corresponds to a velocity component normal to the solar surface. The sawtooth shocks propagate upward in the umbral chromosphere. About $\theta\sim 30^{\circ}$, the shock should be about the middle of the formation layer of the He I triplet, because the shocked rising plasma and the falling plasma equally contribute to form the broadest line in this period. For $\theta\sim 140^{\circ}$, the Doppler velocity reaches its maximum, because the falling upstream component exits the line-forming region and stops to contribute the line formation. For $\theta\sim 320^{\circ}$, the line starts to broaden as the next shock reaches the He I formation layer.

3.2 Double Slab model with the Rankine-Hugoniot relations

Since the temporal variations in the Doppler velocity and the line width indicate that both the upstream and downstream components of the shock should contribute to the line formation, we now repeat the fits to the He I spectral profiles of the umbra with theoretical profiles using an atmospheric model based on two constant property slabs, of which temperatures and macroscopic velocities are constrained by equations (2) and (4). We limit our two-slab fitting analysis to the line profiles from the time during each oscillation when the width reaches its maximum, that is, the contribution to the line intensities from the plasma on each side of the shock is maximized, in the umbra surrounded by the black contour in the Figure 1.

Here, we assume that (1) the shock propagates parallel to the magnetic field. It is reasonable because the umbral flash is suggested to be a slow magnetic hydrodynamic shock (Centeno et al., 2006) and reduces the magneto-hydrodynamic shock to a hydrodynamic shock. Since the sunspot was located near disk center (Figure 1), (2) the shocks are also assumed to propagate along the line-of-sight. These assumptions greatly simplify the analysis. In addition, we assume that (3) the broadening of the triplet absorption profiles are due only to the thermal and the instrumental effects. Effects of the assumption of the zero non-thermal line broadening on the results are discussed in Section 5.

Modeled profiles for the Si I line and the telluric lines remain as before, whereas the two-slab radiative transfer equation for the upstream and downstream components of the shock, which propagates in the formation layer of the He I triplet, is now written as:

[TABLE]

[TABLE]

[TABLE]

where u and d indicate upstream and downstream components, respectively. Thermal broadening is proportional to the square root of the temperature. Here we do not consider unresolved motions and assume the observed line widths and thermal line width, $\Delta\lambda_{{\rm th}}$, are related by $\Delta\lambda_{s}^{2}=\Delta\lambda_{{\rm I}}^{2}+\Delta\lambda_{{\rm th,s}}^{2}$, where $s$ indicates the upstream or downstream plasma. Therefore, using equation (4),

[TABLE]

where $\Delta\lambda_{{\rm I}}=29$ mÅ is the instrumental line width. Assuming a monoatomic ideal gas, we fix $\gamma$ to $5/3$. The line width and the Mach number of the upstream component also determine velocities in the shock frame, $V_{{\rm u}}$ and $V_{{\rm d}}$, of the upstream and downstream components as

[TABLE]

[TABLE]

where $m$ is the mass of helium and $k_{{\rm B}}$ is the Boltzmann constant. We derived the mean molecular weight, $\mu$, as $1.29$ from the solar atomic abundance (Asplund et al., 2009), and umbral atmospheric models (Maltby et al., 1986) at the He I triplet formation height (Felipe et al., 2010). Denoting propagation speed of the shock, $U$, the Doppler shift of both components are calculated from velocities in the observer reference frame as

[TABLE]

and

[TABLE]

In summary, the parameters of the radiative transfer equation are $M$, $U$, $\Delta\lambda_{{\rm u}}$, $S_{{\rm u}}$, $S_{{\rm d}}$, $\tau_{{\rm u}0}$, $\tau_{{\rm d}0}$, $I_{{\rm C}}$, and $\Delta\lambda_{{\rm I}}$.

As for the single slab model fitting (Sec. 3.1), the optical depths ($\tau_{{\rm u}0}$ and $\tau_{{\rm d}0}$) and the intensity in the continuum, $I_{{\rm C}}$, are fixed to be 0.5 and 1, respectively. In addition, we assume the upstream line width, $\Delta\lambda_{{\rm u}}$, to be the minimum value of the line width found from the single slab model fitting just from the previous cycle at each spatial position, that is, each fitting uses a different fixed value for $\Delta\lambda_{{\rm u}}$. The horizontal dashed line in the Figure 4 marks the fixed value of the upstream line width to fit a profile obtained at the time indicated by the vertical dotted line. This approach can be justified by the fact the shock does not affect the temperature upstream in the chromosphere. Therefore, the variable parameters in the fit are only $M$, $U$, $S_{{\rm u}}$, and $S_{{\rm d}}$. Before doing the fitting, we subtracted time averaged residual profiles of the single slab model fitting from the observed spectra to avoid fitting the residuals at 10831.5 Å  with one of double slabs (Figure 2).

An example of the fittings with the double-slab model is shown in Figure 8. The separate slab contributions, as shown by the black dashed and dotted lines, are calculated by solving Equation (6) using parameters of each slab. The line width of the redder component is fixed by using the estimate described above, for the case indicated by the horizontal dashed line in the Figure 4(b), which is the value of the line width at time $t=100$ s, the minimum immediately preceeding the shock peak at time $t=170$ s. The fitting of profiles with double slabs gives us the ratio of the line widths between the upstream and downstream components and their Doppler velocities. From the line-width ratio, we are able to determine the Mach number in the upstream component using equation (11). Substituting the definition of the Mach number into equation (14), the shock propagation speed is given as

[TABLE]

Because $T_{\rm u}$ is equal to $\frac{1}{2}mc^{2}k_{{\rm B}}^{-1}\lambda_{0}^{-2}(\Delta\lambda_{{\rm u}}^{2}-\Delta\lambda_{{\rm I}}^{2})$, the shock speed is derived from the inferred Mach number, the Doppler shift of the upstream component, and the fixed upstream line width.

Our initial values used in the fitting are $M=1.001$, $U=c\lambda_{{\rm Du}}/\lambda_{0}-M\sqrt{\gamma RT_{\rm u}/\mu}$, $S_{{\rm u}}=S$, and $S_{{\rm d}}=0.98$, where $S$ is derived from the single slab model fitting.

We performed multiple non-linear least squares fits using a gradient-expansion algorithm by randomly selecting the $\lambda_{{\rm Du}}$, which determines the initial estimate of the shock speed, within a range $\lambda_{{\rm D}}<\lambda_{{\rm Du}}<\lambda_{{\rm D}}+\Delta\lambda$, until the number of attempts to fit a profile reaches to 1,000 or the number of success fits reaches 50 where success is defined both by a converged least-squares solution and that the two source functions are both less than 0.9, because the source function of the He I triplet is generally smaller than 0.9 when the optical depth is fixed to be 0.5 (Figure 3). We obtained solutions from 19,717 profiles, of which 19,125 profiles provide us 50 successfull fittings with changing the initial estimate of the shock speed, $U$. Here, the fitting solution for a line profile is defined as a parameter set of a fitting that gives us a median value in the Mach number.

The error is defined as the larger of the two values of standard deviation given by the individual solution and the outputs of the successfully randomized fittings. Figure 9 shows histograms of errors of the 19,717 profiles. Peaks of the histograms are at 0.05 in the error of the Mach number and at 1 km s*-1* of the shock speed.

4 Results

We use the SPIES instrument to determine the properties of shocks associated with umbral flashes. Spectral observations taken with each fiber allow us to directly determine the properties of each shock, while the spatial coverage of the IFU allows us to study variations in the spatial and temporal characteristics of the shocks. Figure 10 shows the spatial and temporal variations of the Mach number over the umbra. For this figure, the Mach number at the place of the dead fibers or spaces in the umbra (Figure 1d) is determined from a profile expected by a triangulation. The two dimensional plus temporal data show expanding phase fronts that originate near center of the right umbra, and expanded quasi-circularly to the edge of the umbra similar to what has been reported for umbral flashes (e.g. Rouppe van der Voort et al., 2003; Nagashima et al., 2007). In the figure 10 animation, the evolution of the shocks show spiral patterns (López Ariste et al., 2013, 2016; Su et al., 2016; Jess et al., 2017; Felipe et al., 2019; Kang et al., 2019). The radial propagation of the pattern may be apparent propagation rather than physical wave propagation (Bogdan & Judge, 2006).

Spatial distributions of the average Mach number and average shock speed at each spatial point are displayed in Figure 11. The average values are derived without the values at the dead fibers or spaces. Because the sunspot image moves across the IFU field during the observation, we are able to derive the average values across almost the entire umbra. For the shock speed, negative value means upflow, since the definition of the sign is the same as that of the Doppler shift (opposite to that of the Doppler velocity $V_{{\rm Dop}}$). Roughly speaking, higher Mach numbers are associated with faster shocks. However, the spatial distribution of the Mach number is slightly different from that of the shock speed. The light bridge is the near gap between two lobes of the sunspot. The Mach number seems to be high about the center of each lobe of the sunspot. On the other hand, the shock speed seems to be largest near the light bridge and gradually decreases away from the light bridge.

Figure 12 plots the probability density functions (PDFs) of the (a) Mach number, (b) temperature enhancement, pressure enhancement, density enhancement at the shock front, the (c) shock speed, and the (d) shock heating energy per unit mass of plasma. The PDFs are derived without the values at the dead fibers and spaces. The temperature, pressure, and density enhancements at the shock front are calculated from the Mach number using equations (2), (3) and (4).

The temperature enhancement at the peak of the PDF is approximately equal to 0.24, which is consistent with those reported by Grant et al. (2018) and Houston et al. (2018). The Mach number is 1.24, and the shock speed is $-12$ km s*-1* at the peak of the PDFs. The PDF of the Mach number is fitted with a logarithmic normal distribution $2.60\exp{\{-0.5[(\ln{|M|}-0.226)/0.116]^{2}\}}$, and the PDF of the shock speed is fitted with a normal distribution $0.168\exp{\{-0.5[(U+12.62)/2.29]^{2}\}}$. We derive the shock heating energy per unit mass of plasma with the two models of thermodynamic cycles discussed in Section 3 using the equations (1) and (5) (Figure 12 d). Their PDFs are almost identical, and they have peaks at $2.3\times 10^{10}$ erg ${\rm g^{-1}}$ and $1.8\times 10^{10}$ erg ${\rm g^{-1}}$. Their expected values ($\int q\,{\rm PDF}(q)\,dq$) are $1.8\times 10^{10}$ erg ${\rm g^{-1}}$ and $1.4\times 10^{10}$ erg ${\rm g^{-1}}$. We conclude typical shock heating energy per unit mass of plasma is approximately equal to $2\times 10^{10}$ erg ${\rm g^{-1}}$.

Figure 13 shows two-dimensional histograms of the temperature upstream of the shocks vs. the Mach number (13a) and the shock speed (13a). In Figure 13a, the crosses denote the most likely Mach number for each upstream temperature histogram bin. We find no relation between $T_{{\rm u}}$ and $M$ (linear Pearson correlation coefficient of -0.02). Conversely, we do find a weak linear relation in 13b, with greater upstream temperatures having greater shock speeds (Pearson coefficient -0.32). The red dashed line marks the temperature dependent sound speed ($C_{s}=\sqrt{\gamma RT_{{\rm u}}/\mu}$), with an offset of $1.85$ km s*-1* to take into account the the relative velocity between the sunspot and the observers (Section 3.1). Because the offset of the shock speed for the maximum density at each temperature bin (red crosses) is equal to $-0.5\pm 0.6$ km s*-1*, we conclude that the saw-tooth shocks propagate upward at the sound speed: they are in the weak shock regime.

5 Discussion

We use spectra of the He I 10830 Å triplet from the SPIES instrument to determine the properties of shocks that pass through the umbral chromosphere. From the measured Mach number of the shock and the temperature of the upstream plasma we found shock heating, $q$, of $\sim 2\times 10^{10}$ erg ${\rm g^{-1}}$ per shock cycle, independent of the assumed thermodynamic shock cycle. Moreover, we find that the shocks propagate with the sound speed, and the Mach number does not depend on the temperature.

Does the shock heating energy balance the radiative energy losses in the umbral chromosphere? Since the formation height of the He I 10830 Å triplet is estimated by Felipe et al. (2010) as 1038 - 1208 km from a height where the continuum optical depth at 500 nm is unity, the net radiative cooling rate at the formation layer can be $0.02-0.1$ erg ${\rm cm^{-3}}$ ${\rm s^{-1}}$ (Avrett, 1981). Assuming the density, $\rho$, as $1\times 10^{-11}$ g cm*-3* at the formation layer according to an umbral flash model of Bose et al. (2019), and a typical shock period of $\tau=180$ s from Figure 5, we determine the shock heating rate as $\rho q/\tau=1\times 10^{-3}$ erg ${\rm cm^{-3}}$ ${\rm s^{-1}}$, which is 1 - 5 % of the required amount of energy to compensate the radiative energy losses in the umbral chromosphere. This conclusion is consistent with that of a calculation of the shock heating rates for two umbral atmospheric models (Lee & Yun, 1985). We performed the same fittings as that described in section 3.2 but using different values of the optical depths from 0.5 ($\tau_{{\rm u}0}=\tau_{{\rm d}0}=0.2$ and $\tau_{{\rm u}0}=\tau_{{\rm d}0}=0.8$). Their results do not change the conclusion that the shock heating rate is insufficient to maintain the umbral chromosphere.

We neglected the contribution of non-thermal motions to the spectral line broadening. However, Bose et al. (2019) derived the non-thermal velocity to be $\sim 4$ km s*-1* at the formation layer of the He I triplet for the umbral flash. Figure 14 shows expected values of the shock heating energy rate as a function of the non-thermal velocity, $V_{nth}$. First, we calculated the Mach number, $M^{\prime}$, solving equation 4 from ratios of redefined temperatures as $T^{\prime}_{{\rm u}}\equiv T_{{\rm u}}-mV_{nth}^{2}/2k_{B}$ and $T^{\prime}_{{\rm d}}\equiv T_{{\rm d}}-mV_{nth}^{2}/2k_{B}$ with an assumption that $V_{nth}$ upstream and downstream of the shocks are the same. Next, the shock heating energy per unit mass of plasma for the two models of thermodynamic cycles is derived with equations (1) and (5) from $M^{\prime}$ and $T^{\prime}_{{\rm u}}$, and its expected values are obtained as described in section 4. Finally, the plotted shock heating energy rates are estimated as discussed above using the expected values. Although the heating energy rates increase with the non-thermal velocity, the rates at the non-thermal velocity of 4 km s*-1* are insufficient to maintain the umbral chromosphere.

The derived upstream temperature exhibits a significant linear relation with the shock speed but not with the Mach number. To demonstrate the relations are not systematic results from the fitting, we performed a Monte Carlo simulation by synthesizing and fitting 10,000 spectral profiles calculated with randomly distributed parameters within the measured typical range $-0.2$ Å $<\lambda_{{\rm D}}<0.35$ Å and $0.2$ Å $<\Delta\lambda<0.65$ Å (Figure 3). The upstream line width, $\Delta\lambda_{{\rm u}}$, is fixed to be a randomly distributed parameter within a range $0.2\,{\rm\AA}<\Delta\lambda_{{\rm u}}<\Delta\lambda$.Approximately 10% of the synthetic spectra could not be fit successfully according to the criteria described in section 3.2. Figure 15 displays the fit parameters of the successfully fitted profiles.

To understand how the model varies with each parameter we can consider isoparametric lines by holding one parameter fixed at a time. For instance, the upper right boundary in Figure 15 (a) corresponds to solutions with fixed line width of the single slab model, which we use for the synthesis: red pluses mark solutions for which $0.64$ Å $<\Delta\lambda<$ $0.65$ Å. For Figure 15 (b), the lower left edge is formed by solutions for the synthetic spectral profiles calculated with Doppler shifts between $-0.20$ Å and $-0.19$ Å. We explain this behavior by noting that the Doppler shift of the upstream component $\lambda_{{\rm Du}}$ is nearly the Doppler shift of the spectral line itself (see the dashed and solid lines in Figure 8). Making that approximation in Equation (16), we find that

[TABLE]

Since the Mach number does not have significant dependence on the upstream temperature (Figure 13 a), the absolute value of the shock speed increases with the upstream temperature for a fixed Doppler shift.

There is no systematic relation within the distribution of the test solutions without boundaries, and their boundaries are determined from the observed parameter ranges. In addition, the shock speed is roughly determined from the independent parameters, Doppler shift and the upstream temperature (Equation 17). Therefore, we conclude that the observed significant relations among the shock speed, the Mach number, and the upstream temperature do not systematically result from the fitting method.

Bogdan et al. (2003) performed two dimensional MHD simulations of wave propagation in a variety of solar-like magnetized atmosphere. They introduced waves using sinusoidal photospheric driving. The wave amplitude of velocity fluctuations normalized by the sound speed increased as $\rho^{-2/5}$ due to the decrease in density with height in the stratified atmosphere. Eventually, the waves gave rise to shocks. They found that slow MHD shocks propagate slightly faster than the sound speed in strong magnetic field concentrations, such as sunspots. Although it is not significant, our measured shock propagation speed tends to be faster than the sound speed by $0.5$ km s*-1* (Figure 13b).

In the Bogdan et al. (2003) simulations, after the formation of the shocks, the velocity saturates due to the numerical viscosity. If the viscosity was negligibly small as expected in the real solar atmosphere (Vranjes & Krstic, 2013), the Mach number could be written as

[TABLE]

where $h$ is the height above the photosphere, $v_{0}$ is the amplitude of the velocity fluctuation at the photosphere, $\Lambda$ is the density scale height ($=RT/\mu g$), and $g$ is the gravity acceleration. If the formation height of the He I triplet is $1000$ km (Felipe et al., 2010) and $v_{0}/C_{s}=0.05$, then the Mach number at $h=1000$ km is approximately equal to 1.2 for $T=1.25\times 10^{4}$ K and 1.1 for $T=2.50\times 10^{4}$ K. Therefore, the Mach number would decrease with increasing temperature. However, we observed the Mach number is independent of the upstream temperature. Because $v_{0}$ should not depend on $H$ and $T$, we propose that the Mach number may saturate in the umbral atmosphere due to energy loss of the shocks along a flux tube. As examples, using magnetohydrodynamic simulations of magnetic flux tubes, Takasao et al. (2013) shows the slow shock energy should be carried across magnetic field lines by the fast-mode magnetohydrodynamic waves generated at a place where the sound speed is equal to the Alfv${\rm\acute{e}}$n speed, as produced by the shocks, and Shelyag et al. (2016) shows efficient ambipolar dissipation of Alfv${\rm\acute{e}}$n waves transformed from acoustic shocks. In addition, if the shock does not propagate parallel to the magnetic field lines and it has not yet reached the steady state, intermediate shock substructures within the shock may dissipate the acoustic energy (Snow & Hillier, 2019).

6 Summary

We derived the shock heating energy per unit mass of plasma from the spectra in the near infrared He I triplet measured with the IFU spectrometer, SPIES, on the Dunn Solar Telescope. The SPIES, which is a prototype instrument of the DL-NIRSP of the forthcoming DKIST, allows us to compare shock properties with that of heating signatures over the entire umbra, even though they change rapidly in time.

In order to determine the shock parameters, we fit the measured spectral profiles in the He I 10830 Å triplet with theoretical profiles computed with the radiative transfer equation using an atmospheric model based on two constant property slabs with temperatures and macroscopic velocities constrained by the Rankine-Hugoniot relations. As a result, the typical shock heating energy per unit mass of plasma is $\sim 2\times 10^{10}$ erg ${\rm g^{-1}}$ per shock cycle, which is insufficient to maintain the umbral chromosphere.

A much stronger conclusions are that shock propagation at the sound speed is consistent with them being weak shocks , and the Mach number does not depend on the temperature upstream of the shocks. If the viscosity is negligibly small as expected in the actual solar atmosphere, the Mach number should decrease with the temperature increase. Therefore, we propose that energy loss of the shocks may suppress the amplification of the Mach number in the umbral atmosphere.

We extend thanks to Prof. J. R. Kuhn, Dr. G. I. Dima, Dr. X. Sun, Dr. M. Kramar, and Dr. A. Fehlmann who discussed with us every week. The authors are also grateful to Dr. H. Lin for providing the SPIES, and Dr. A. Hiller for discussing the interpretation. The National Solar Observatory (NSO) is operated by the Association of Universities for Research in Astronomy, Inc. (AURA), under cooperative agreement with the National Science Foundation. The SDO data are provided courtesy of NASA/SDO and the AIA and HMI science teams.