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

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é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 ∼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 has 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 narrowband Fabry–Pérot-based imaging spectrograph. Using similar observations, Joshi & de la Cruz Rodríguez (2018) reported that the temperature at the chromosphere increases from 3700 K up to 6200 K. Grant et al. (2018) also measured significant temperature enhancements up to a maximum of ∼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 ±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 × 10⁶ erg cm⁻² s⁻¹ based on the net radiative cooling rate for a semiempirical 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 × 10⁷ erg cm⁻² s⁻¹ in the photosphere and 8.3 × 10⁴ erg cm⁻² s⁻¹ 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 the 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, give our conclusions (Section 6).

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¹ ) of the National Science Foundation'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 × 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 × 60 spatial pixels, corresponding to 15 × 16 arcsec² 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 mosaicked by steering the optical field with a two-axis field steering mirror. In our observations below, we perform a 4 × 4 mosaic for a total field of view of 59 × 65 arcsec² and a temporal cadence of 14 s.

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 the 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 with 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 that 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) on board 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 Figures 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.

Zoom In Zoom Out Reset image size

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 (Schatzman 1949; Weymann 1960) 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 a constant specific volume until its entropy matches that of the unshocked gas, after which it adiabatically expands to the unshocked pressure and specific volume. Following Bray & Loughhead (1974, Section 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

$$\begin{eqnarray}&&q=\displaystyle \frac{1}{\gamma -1}\displaystyle \frac{R}{\mu }\,{T}_{{\rm{u}}}\displaystyle \frac{{\rho }_{{\rm{u}}}}{{\rho }_{{\rm{d}}}}\left[\displaystyle \frac{{p}_{{\rm{d}}}}{{p}_{{\rm{u}}}}-{\left(\displaystyle \frac{{\rho }_{{\rm{d}}}}{{\rho }_{{\rm{u}}}}\right)}^{\gamma }\right],\end{eqnarray} \tag{ 1 }$$

where γ is the polytropic index (the ratio of the specific heats), R is the gas constant, μ is the mean molecular weight, T is the temperature, ρ 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

$$\begin{eqnarray}&&\displaystyle \frac{{\rho }_{{\rm{d}}}}{{\rho }_{{\rm{u}}}}=\displaystyle \frac{{V}_{{\rm{u}}}}{{V}_{{\rm{d}}}}=\displaystyle \frac{(\gamma +1){M}^{2}}{(\gamma -1){M}^{2}+2},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{p}_{{\rm{d}}}}{{p}_{{\rm{u}}}}=\displaystyle \frac{2\gamma {M}^{2}-\gamma +1}{\gamma +1},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{T}_{{\rm{d}}}}{{T}_{{\rm{u}}}}=\displaystyle \frac{\left[2\gamma {M}^{2}-\gamma +1\right]\left[\left(\gamma -1\right){M}^{2}+2\right]}{{\left(\gamma +1\right)}^{2}{M}^{2}},\end{eqnarray} \tag{ 4 }$$

where V_u and V_d are velocities measured in the reference frame comoving 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

$$\begin{eqnarray}&&q=\displaystyle \frac{\gamma }{\gamma -1}\displaystyle \frac{R}{\mu }\,{T}_{{\rm{u}}}\displaystyle \frac{{\rho }_{{\rm{u}}}}{{\rho }_{{\rm{d}}}}\left[{\left(\displaystyle \frac{{p}_{{\rm{d}}}}{{p}_{{\rm{u}}}}\right)}^{\tfrac{1}{\gamma }}-\displaystyle \frac{{\rho }_{{\rm{d}}}}{{\rho }_{{\rm{u}}}}\right].\end{eqnarray} \tag{ 5 }$$

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. (2000a, 2000b, 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.

Zoom In Zoom Out Reset image size

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

$$\begin{eqnarray}&&I(\lambda )=S(1-{e}^{-{\tau }_{\lambda }})+{I}_{{\rm{C}}}\,{e}^{-{\tau }_{\lambda }},\end{eqnarray} \tag{ 6 }$$

where λ is the wavelength, S is the source function of the slab, and I_C is the intensity in the continuum. Optical thickness of the He i triplet as a function of the wavelength, τ_λ, can be written as

$$\begin{eqnarray}&&{\tau }_{\lambda }=\tau \displaystyle \sum _{i=1}^{3}{\omega }_{i}\exp \left\{-{\left(\displaystyle \frac{\lambda -{\lambda }_{i}-{\lambda }_{{\rm{D}}}}{{\rm{\Delta }}\lambda }\right)}^{2}\right\},\end{eqnarray} \tag{ 7 }$$

where τ is the optical thickness of the slab at the spectral line center, ω_i are the statistical weights derived from the number of the magnetic sublevels of the excited state under a constraint of ${\sum }_{i=1}^{3}{\omega }_{i}=1$, λ_i are the center wavelengths of the triplet lines, λ_D is the Doppler shift, and Δλ 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_C to 1 as we normalize each intensity profile by I_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 τ to be 0.5 (see Figure 6 of Schad et al. 2015). The quality of the single-component fits is 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 ∼3.7 × 10⁵ 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 (∼99.8%) of the Doppler shifts and the line widths vary within the ranges −0.2 Å < λ_D < 0.35 Å and 0.2 Å < Δλ < 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⁻¹ 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.

Zoom In Zoom Out Reset image size

Figure 4 shows temporal variations of the inferred Doppler velocity, ${V}_{\mathrm{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 λ₀ = 10830 Å is the wavelength at the line center. When V_Dop is larger than −1.85 km s⁻¹, 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).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In order to examine phase relations between the Doppler velocity and the line width statistically, we construct two-dimensional histograms of their temporal variations versus the phase as displayed in Figure 6. We used values of the umbral data enclosed within the contour in 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.

Zoom In Zoom Out Reset image size

The temporal variation in the V_Dop reveals a sawtooth pattern, that is, the phase of the maximum V_Dop is smaller than the half-period. In addition, the line width has its maximum value during the time V_Dop is increasing. The sawtooth pattern and the enhancement of the line width during the V_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 θ ∼ 30°, 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 θ ∼ 140°, the Doppler velocity reaches its maximum, because the falling upstream component exits the line-forming region and stops to contribute the line formation. For θ ∼ 320°, the line starts to broaden as the next shock reaches the He i formation layer.

Zoom In Zoom Out Reset image size

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 magnetohydrodynamic 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 nonthermal 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

$$\begin{eqnarray}\begin{array}{rcl}I(\lambda ) & = & {S}_{{\rm{u}}}(1-{e}^{-{\tau }_{{\rm{u}}}(\lambda )})+{S}_{{\rm{d}}}{e}^{-{\tau }_{{\rm{u}}}(\lambda )}(1-{e}^{-{\tau }_{{\rm{d}}}(\lambda )})\\ & & +\,{I}_{{\rm{C}}}\,{e}^{-{\tau }_{{\rm{u}}}(\lambda )-{\tau }_{{\rm{d}}}(\lambda )},\end{array}\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{\tau }_{{\rm{u}}}(\lambda )={\tau }_{{\rm{u}}0}\displaystyle \sum _{i=1}^{3}{\omega }_{i}\exp \left\{-{\left(\displaystyle \frac{\lambda -{\lambda }_{i}-{\lambda }_{\mathrm{Du}}}{{\rm{\Delta }}{\lambda }_{{\rm{u}}}}\right)}^{2}\right\},\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{\tau }_{{\rm{d}}}(\lambda )={\tau }_{{\rm{d}}0}\displaystyle \sum _{i=1}^{3}{\omega }_{i}\exp \left\{-{\left(\displaystyle \frac{\lambda -{\lambda }_{i}-{\lambda }_{\mathrm{Dd}}}{{\rm{\Delta }}{\lambda }_{{\rm{d}}}}\right)}^{2}\right\},\end{eqnarray} \tag{ 10 }$$

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, Δλ_th, are related by ${\rm{\Delta }}{\lambda }_{s}^{2}={\rm{\Delta }}{\lambda }_{{\rm{I}}}^{2}+{\rm{\Delta }}{\lambda }_{\mathrm{th},{\rm{s}}}^{2}$, where s indicates the upstream or downstream plasma. Therefore, using Equation (4)

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Delta }}{\lambda }_{{\rm{d}}} & = & \left[{\rm{\Delta }}{\lambda }_{\mathrm{th},{\rm{u}}}^{2}\displaystyle \frac{\left(2\gamma {M}^{2}-\gamma +1\right)\left[\left(\gamma -1\right){M}^{2}+2\right]}{{\left(\gamma +1\right)}^{2}{M}^{2}}\right.\\ & & {\left.+{\rm{\Delta }}{\lambda }_{{\rm{I}}}^{2}\right]}^{0.5},\end{array}\end{eqnarray} \tag{ 11 }$$

where Δλ_I = 29 mÅ is the instrumental line width. Assuming a monoatomic ideal gas, we fix γ to 5/3. The line width and the Mach number of the upstream component also determine velocities in the shock frame, V_u and V_d, of the upstream and downstream components as

$$\begin{eqnarray}&&{V}_{{\rm{u}}}=M{\left\{\displaystyle \frac{\gamma {{Rmc}}^{2}}{2\mu {k}_{{\rm{B}}}{\lambda }_{0}^{2}}({\rm{\Delta }}{\lambda }_{{\rm{u}}}^{2}-{\rm{\Delta }}{\lambda }_{{\rm{I}}}^{2})\right\}}^{0.5},\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{V}_{{\rm{d}}}={V}_{{\rm{u}}}\displaystyle \frac{(\gamma -1){M}^{2}+2}{(\gamma +1){M}^{2}},\end{eqnarray} \tag{ 13 }$$

where m is the mass of helium and k_B is the Boltzmann constant. We derived the mean molecular weight, μ, 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

$$\begin{eqnarray}&&{\lambda }_{\mathrm{Du}}={\lambda }_{0}({V}_{{\rm{u}}}+U)/c,\end{eqnarray} \tag{ 14 }$$

and

$$\begin{eqnarray}&&{\lambda }_{\mathrm{Dd}}={\lambda }_{0}({V}_{{\rm{d}}}+U)/c.\end{eqnarray} \tag{ 15 }$$

In summary, the parameters of the radiative transfer equation are M, U, Δλ_u, S_u, S_d, τ_u0, τ_d0, I_C, and Δλ_I.

As for the single-slab-model fitting (Section 3.1), the optical depths (τ_u0 and τ_d0) and the intensity in the continuum, I_C, are fixed to be 0.5 and 1, respectively. In addition, we assume the upstream line width, Δλ_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 Δλ_u. The horizontal dashed line in 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_u, and S_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 a double slab (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 preceding 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

$$\begin{eqnarray}&&U=\displaystyle \frac{{\lambda }_{\mathrm{Du}}}{{\lambda }_{0}}c-M\sqrt{\displaystyle \frac{\gamma {{RT}}_{{\rm{u}}}}{\mu }}.\end{eqnarray} \tag{ 16 }$$

Because T_u is equal to $\tfrac{1}{2}{{mc}}^{2}{k}_{{\rm{B}}}^{-1}{\lambda }_{0}^{-2}({\rm{\Delta }}{\lambda }_{{\rm{u}}}^{2}-{\rm{\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.

Zoom In Zoom Out Reset image size

Our initial values used in the fitting are M = 1.001, $U=c{\lambda }_{\mathrm{Du}}/{\lambda }_{0}-M\sqrt{\gamma {{RT}}_{{\rm{u}}}/\mu }$, S_u = S, and S_d = 0.98, where S is derived from the single-slab-model fitting.

We performed multiple nonlinear least-squares fits using a gradient-expansion algorithm by randomly selecting the λ_Du, which determines the initial estimate of the shock speed, within a range λ_D < λ_Du < λ_D + Δλ, until the number of attempts to fit a profile reaches 1000 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 successful 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⁻¹ of the shock speed.

Zoom In Zoom Out Reset image size

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 1(d)) 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_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.

Zoom In Zoom Out Reset image size

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)–(4).

Zoom In Zoom Out Reset image size

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⁻¹ at the peak of the PDFs. The PDF of the Mach number is fitted with a logarithmic normal distribution $2.60\exp \{-0.5[(\mathrm{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 Equations (1) and (5) (Figure 12(d)). Their PDFs are almost identical, and they have peaks at 2.3 × 10¹⁰ erg g⁻¹ and 1.8 × 10¹⁰ erg g⁻¹. Their expected values ($\int q\,\mathrm{PDF}(q)\,{dq}$) are 1.8 × 10¹⁰ erg g⁻¹ and 1.4 × 10¹⁰ erg g⁻¹. We conclude typical shock heating energy per unit mass of plasma is approximately equal to 2 × 10¹⁰ erg g⁻¹.

Figure 13 shows two-dimensional histograms of the temperature upstream of the shocks versus the Mach number (13(a)) and the shock speed (13(b)). In Figure 13(a), the crosses denote the most likely Mach number for each upstream temperature histogram bin. We find no relation between T_u and M (linear Pearson correlation coefficient of −0.02). Conversely, we do find a weak linear relation in 13(b), 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⁻¹ to take into account 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 ± 0.6 km s⁻¹, we conclude that the sawtooth shocks propagate upward at the sound speed: they are in the weak shock regime.

Zoom In Zoom Out Reset image size

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 ∼2 × 10¹⁰ erg g⁻¹ 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 cm⁻³ s⁻¹ (Avrett 1981). Assuming the density, ρ, as 1 × 10⁻¹¹ g cm⁻³ at the formation layer according to an umbral flash model of Bose et al. (2019), and a typical shock period of τ = 180 s from Figure 5, we determine the shock heating rate as ρq/τ = 1 × 10⁻³ erg cm⁻³ s⁻¹, 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 (τ_u0 = τ_d0 = 0.2 and τ_u0 = τ_d0 = 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 nonthermal motions to the spectral line broadening. However, Bose et al. (2019) derived the nonthermal velocity to be ∼4 km s⁻¹ 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 nonthermal velocity, V_nth. First, we calculated the Mach number, M', solving Equation (4) from ratios of redefined temperatures as ${T}_{{\rm{u}}}^{{\prime} }\equiv {T}_{{\rm{u}}}-{{mV}}_{{nth}}^{2}/2{k}_{B}$ and ${T}_{{\rm{d}}}^{{\prime} }\equiv {T}_{{\rm{d}}}-{{mV}}_{{nth}}^{2}/2{k}_{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' and ${T}_{{\rm{u}}}^{{\prime} }$, 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 nonthermal velocity, the rates at the nonthermal velocity of 4 km s⁻¹ are insufficient to maintain the umbral chromosphere.

Zoom In Zoom Out Reset image size

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 ranges −0.2 Å < λ_D < 0.35 Å and 0.2 Å < Δλ < 0.65 Å (Figure 3). The upstream line width, Δλ_u, is fixed to be a randomly distributed parameter within a range 0.2 Å < Δλ_u < Δλ. 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.

Zoom In Zoom Out Reset image size

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 Å < Δλ < 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 λ_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

$$\begin{eqnarray}&&U\sim \displaystyle \frac{{\lambda }_{{\rm{D}}}}{{\lambda }_{0}}c-M\sqrt{\displaystyle \frac{\gamma {{RT}}_{{\rm{u}}}}{\mu }}.\end{eqnarray} \tag{ 17 }$$

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 ρ^−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⁻¹ (Figure 13(b)).

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

$$\begin{eqnarray}&&M\approx 1+\displaystyle \frac{{v}_{0}}{{C}_{s}}\exp \left(\displaystyle \frac{2h}{5{\rm{\Lambda }}}\right),\end{eqnarray} \tag{ 18 }$$

where h is the height above the photosphere, v₀ is the amplitude of the velocity fluctuation at the photosphere, Λ is the density scale height (=R T/μ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₀/C_s = 0.05, then the Mach number at h = 1000 km is approximately equal to 1.2 for T = 1.25 × 10⁴ K and 1.1 for T = 2.50 × 10⁴ K. Therefore, the Mach number would decrease with increasing temperature. However, we observed the Mach number is independent of the upstream temperature. Because v₀ 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) show 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én speed, as produced by the shocks, and Shelyag et al. (2016) shows efficient ambipolar dissipation of Alfvé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).

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 DST. 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 ∼2 × 10¹⁰ erg g⁻¹ per shock cycle, which is insufficient to maintain the umbral chromosphere.

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.

Facility: DST(SPIES). -

Software: Solar Soft.