Properties of shock waves in the quiet Sun chromosphere

Short-lived (100s or less), sub-arcsec to a couple of arcsec size features of enhanced brightenings in the narrowband images at the $\mathrm{H_{2V}}$ and $\mathrm{K_{2V}}$ positions of the Ca II H and K lines in the quiet Sun are known as bright grains. With simultaneous observations of a quiet Sun internetwork region in the Fe I 6173 {\AA}, Ca II 8542 {AA}, and Ca II K lines acquired by the CRisp Imaging Spectro-Polarimeter and the CHROMospheric Imaging Spectrometer instruments on the Swedish 1-m Solar Telescope, we performed multi-line non-local thermodynamic equilibrium inversions using the STockholm inversion Code to infer the time-varying stratified atmosphere's physical properties such as the temperature, line-of-sight (LOS) velocity, and microturbulence. The Ca II K profiles of bright grains show enhancement in the $\mathrm{K_{2V}}$ peak intensities with absence of the $\mathrm{K_{2R}}$ features. At the time of maximum enhancement in the $\mathrm{K_{2V}}$ peak intensities, we found average enhancements in temperature at lower chromospheric layers (at $\log\tau_{500}$ $\simeq$ $-$4.2) of about 1.1 kK with maximum enhancement of about 4.5 kK. These temperature enhancements are colocated with upflows, as strong as $-$6 $\mathrm{km\;s^{-1}}$, in the direction of the LOS. The LOS velocities at upper chromospheric layers at $\log\tau_{500}$<$-$4.2 show consistent downflows greater than $+$8 $\mathrm{km\;s^{-1}}$. The retrieved value of microturbulence in the atmosphere of bright grains is negligible at chromospheric layers. The study provides observational evidence to support the interpretation that the bright grains observed in narrowband images at the $\mathrm{H_{2V}}$ and $\mathrm{K_{2V}}$ positions of the Ca II H and K lines are manifestations of upward propagating acoustic shocks against a background of downflowing atmospheres.


Introduction
Waves are omnipresent in the solar atmosphere.Understanding the source, nature, and energy dissipation mechanism of waves in the solar atmosphere is essential to know possible contributions to chromospheric and coronal heating.
Narrowband images at the H 2V and K 2V positions of the Ca ii H&K lines in the internetwork region of chromosphere show short-lived (100 s or less), sub-arcsec to a couple of arcsec sized brightenings called grains (Beckers 1964).These bright grains have also been observed in the Transition Region and Coronal Explorer near-UV passbands (Rutten et al. 1999;Handy et al. 1999;Tian et al. 2010).The typical size of these grains in the Ca ii K filtergrams observed with the German Vacuum Tower Telescope (von der Lühe 1998) and narrowband imaging spectroscopy in the Ca ii 8542 Å line from the Dunn Solar Telescope (DST; Dunn 1969) are on the order of 1 .95 (Wöger et al. 2006) and 0 .8 (Vecchio et al. 2009), respectively.Decades of observations have suggested that these are manifestations of acoustic-shock-like disturbances that propagate upward in the solar atmosphere (Cram 1974(Cram , 1976;;Zirin 1974;Cram et al. 1977;Carlsson & Stein 1997;Beck et al. 2008).These grains are not known to be caused by magnetism-related processes (Kalkofen 1990;Rutten & Uitenbroek 1991;Lites et al. 1999).We refer the reader to Rutten & Uitenbroek (1991; and the references therein) for an extensive review on early efforts of the observation, simulation, and interpretation of grains.Carlsson & Stein (1992, 1997) reproduced the bright grains in the Ca ii K line using one-dimensional hydrodynamic models with radiation in non-local thermodynamic equilibrium (non-LTE) assuming complete redistribution (CRD).Quiet-Sun acoustic shock waves were studied by Wedemeyer et al. (2004) using three-dimensional simulation assuming LTE conditions.These simulations showed that upward-propagating acoustic waves from the lower atmosphere turn into shocks in the chromosphere due to a drop in the gas density by several orders of magnitude (Carlsson & Stein 1992;Wedemeyer et al. 2004;Danilovic et al. 2014).
The density decreases while the shock continues to travel upward from the chromosphere to the transition region, leading to the steepening of the shock and an increase in the shock amplitude.As a result, one would expect to see a transition region counterpart of the bright grains.However, observational studies in the 1990's (Judge et al. 1997;Steffens et al. 1997) from the Solar Ultraviolet Measurements of Emitted Radiation (SUMER, Wilhelm et al. 1995) instrument on board the Solar and Heliospheric Observatory (Domingo et al. 1995) suggested that the bright grains do not have a transition-region counterpart.Thus, it led to the conclusion that upward-propagating acoustic shock waves do not play a significant role in the heating of the transition region (Judge et al. 1997).More recently, however, using observations from the Interface Region Imaging Spectrograph (IRIS, De Pontieu et al. 2014) and the Swedish 1-m Solar Telescope (SST, Scharmer et al. 2003a), Martínez-Sykora et al. (2015) showed that some bright grains have signatures in the Si iv 1403 Å line, suggesting that some acoustic shock waves can in fact reach the transition region.Using imaging spectroscopic observations of the quiet internetwork chromosphere in the Ca ii 8542 Å line, Vecchio et al. (2009) suggested that some acoustic waves may be unable to propagate to the transition region because horizontal magnetic field lines (the canopy) prevent acoustic shock waves from propagating to the outer solar atmosphere.
Semi-empirical models and magneto-hydrodynamics (MHD) simulations have also suggested that upward-propagating acoustic shocks cause brightenings in the millimeter continuum (Loukitcheva et al. 2006;Wedemeyer-Böhm et al. 2007;da Silva Santos et al. 2020;Eklund et al. 2021a).Through analysis of data from the Atacama Large Millimeter Array (ALMA; Wootten & Thompson 2009), Eklund et al. (2020) found a typical chromospheric temperature enhancement of 400-750 K during the shock with a maximum enhancement of 1200 K.Moreover, Eklund et al. (2021b) suggested that the ALMA observations studied by Eklund et al. (2020) had limited spatial resolution, which could lead to reduced brightness temperatures and an underestimation of the temperature enhancement.
The spectra of grains in the Ca ii H line recorded by the POlarimetric LIttrow Spectrograph (POLIS, Beck et al. 2005) was inverted by Beck et al. (2013) using the LTE approximation.They found an increase in temperature of 200-300 K in the quiet-Sun chromosphere.
The limitations in these earlier observational studies, such as the limited spatial resolution of ALMA and the LTE approximation used in inversions of the Ca ii H spectra, merit further observational work with advanced methods.In this study, we use high-resolution imaging spectropolarimetric observations of the bright grains in the Ca ii K line from the SST to infer the time-varying stratified atmospheric properties, such as the temperature, Doppler shift, and microturbulence using a non-LTE inversion code.

Observations
The observations were made with SST on 6 June 2019 between 10:26:20-10:40:06 UT.The observed field of view (FOV) consists of a quiet-Sun scene centered at solar (x, y) = (631 , −23 ), corresponding to a viewing angle cos θ = µ = 0.7.Here, θ represents the angle with respect to the solar normal.The observations were obtained with the CRisp Imaging Spectro-Polarimeter (CRISP; Scharmer et al. 2008) and CHROMospheric Imaging Spectrometer (CHROMIS; Scharmer 2017) instruments simultaneously.The observed FOV of CHROMIS covers an area of 50 × 33 , which overlaps with the CRISP FOV of 41 × 37 , as shown in Fig. 1.
The CRISP instrument sampled the Ca ii 8542 Å line with 20 non-equidistant wavelength points spanning from −1.69 Å to +2.36 Å, and the Fe i 6173 Å line was observed in 14 nonequidistant wavelength points between −0.35 Å and +0.64 Å.
We acquired observations in these lines with full-Stokes polarimetry.The pixel scales of the CRISP data are 0 .059 and the observation sequences have a cadence of 37 s.
With CHROMIS, we observed the Ca ii K line in 29 wavelength positions around the line center (3933.682Å), with 65 mÅ sampling.The CHROMIS instrument also observed one additional wavelength position in the continuum at 4000 Å.The data from the CHROMIS instrument have a cadence of 8.26 s and a pixel scale of 0 .038.
Wide-band (WB) images using an auxiliary channel with a filter centered at 3950 Å (FWHM = 6.5 Å) and a filter centered at the Hα line (FWHM = 4.9 Å) are also observed using the CHROMIS and CRISP instruments, respectively.The details of the optical setup, along with the passbands of the different filters, are given in Löfdahl et al. (2021).
High-quality spatial resolution was achieved by the combination of good seeing conditions, the adaptive optics system, and the excellent CRISP and CHROMIS re-imaging systems (Scharmer et al. 2003b(Scharmer et al. , 2019)).The raw data were then processed using the CRISPRED (de la Cruz Rodríguez et al. 2015b) and CHROMISRED (Löfdahl et al. 2021) pipelines for CRISP and CHROMIS data, respectively.The reduction process includes image restoration by using Multi-Object Multi-Frame Blind Deconvolution (MOMFBD;van Noort et al. 2005).The images from both the instruments were de-rotated to account for diurnal field rotation, and aligned and de-stretched to remove warping due to seeing effects, before making them ready for scientific analyses.The CRISP data were co-aligned with CHROMIS data using cross-correlation of photospheric WB channels.Before the alignment, the CRISP data were resampled to the CHROMIS pixel scale, and their cadences were adjusted to CHROMIS cadences by using nearest-neighbor interpolation.The last step in the post-process procedure was to perform absolute intensity calibration of the data by comparing average line profiles observed at the disk center with the Hamburg atlas spectrum (Neckel & Labs 1984).We also performed Milne-Eddington (ME) inversions of the Fe i 6173 Å data to infer the LOS magnetic field (B LOS ) utilizing the pyMilne code, a parallel C++/Python implementation 1 (de la Cruz Rodríguez 2019).
Figure 1  Ca ii 8542 Å line (at −15.8 mÅ, see panel g) but not in the red wing images (at +64.2 mÅ, see panel h).These compact brightenings are prevalent in the internetwork region and their typical lifetime is found to be shorter than 60 s.These small-scale and short-lived brightenings in the Ca ii K narrowband images are referred to as chromospheric bright grains (CBGs), which are the focus of this study.Prior to the onset of CBG, the K 1 and the K 2 features are absent, and only the K 3 absorption is present.The nomenclature of the K 1 , K 2 , and K 3 spectral features in the Ca ii K line are described in Fig. 1 of Rutten & Uitenbroek (1991).At the onset of the CBG, we see an enhancement in the intensity of the K 2V peak, which increases as the CBG reaches its maximum with the K 3 feature redshifted (∆λ = +153 mÅ).At the end of the CBG, the Ca ii K profiles appear similar to the

Identification of the CBG-like spectral profiles
High-resolution and high-cadence spectropolarimetric solar observations typically generate large data volumes.In the recent past, many authors have used algorithmic ways to reduce the dimensionality of the data, which facilitates efficient, qualitative, and statistical analysis of millions of data points.The k-means clustering is one such technique that has been extensively used to analyze solar spectra (e.g., Panos et al. 2018;Sainz Dalda et al. 2019;Nóbrega-Siverio et al. 2021;Joshi & Rouppe van der Voort 2022).We employed the k-means technique to partition the Ca ii K spectral profiles into k = 100 clusters as described in Appendix A. Each k-means cluster is represented by the mean of all the profiles belonging to that cluster, which we refer to as the RP.The CRISP data have lower cadence and lower spatial resolution compared to the CHROMIS data, and therefore the evolution of the several CBGs could have been missed in the CRISP data, as illustrated in Fig. 2. Hence, we performed the k-means clustering solely on the Ca ii K data from CHROMIS and did not include the Ca ii 8542 Å and Fe i 6173 Å lines from CRISP.We projected the clusters over the co-spatial Ca ii 8542 Å and the Fe i 6173 Å spectra to calculate the A153, page 5 of 21 shown by blue shaded areas in Fig. 5, which is close to the mean profiles in each cluster (shown in black), implying a good partitioning through the k-means algorithm.

Inversions
We applied the MPI-parallel STockholm inversion Code (STiC, de la Cruz Rodríguez et al. 2016, 2019) to retrieve the atmospheric parameters and the evolution of the selected CBGs.Based on a modified version of the RH radiative transfer code (Uitenbroek 2001), STiC uses the cubic Bezier solvers to solve the polarized radiative transfer equation (de la Cruz Rodríguez & Piskunov 2013).It can fit multiple spectral lines simultaneously in non-LTE assuming statistical equilibrium.It takes into account partial redistribution (PRD) effects using a fast approximation method (for more details see Leenaarts et al. 2012).It also fits the intensity in each pixel, assuming a plane-parallel atmosphere (also called the 1.5D approximation).The equation-of-state utilized in STiC is obtained from the library functions in Spectroscopy Made Easy (SME) package code (Piskunov & Valenti 2017).
We inverted the Ca ii K, Ca ii 8542 Å, and Fe i 6173 Å lines simultaneously.We used a six-level Ca ii atom.The Ca ii H&K profiles were synthesized in PRD (Milkey & Mihalas 1974;Uitenbroek 1989), and the Ca ii IR profiles were modeled using the CRD approximation.The atomic parameters for the Fe i 6173 Å line were retrieved from Kurucz's line lists (Kurucz 2011) and synthesized in the LTE approximation.The line broadening cross-sections of radiative transitions of the Ca ii atom for collisions with neutral Hydrogen were retrieved from the table published by Barklem et al. (2000).The Hydrogen broadening Cross-Section Calculator (abo-cross2 ) code was used to calculate line broadening cross-sections of the Fe i 6173 Å line (Anstee & O'Mara 1995;Barklem & O'Mara 1997;Barklem et al. 1998a,b).
The STockholm inversion Code iteratively perturbs the atmospheric parameters such as temperature (T ), LOS velocity (V LOS ), and microturbulence (V turb ) to minimize the χ 2 distance between the synthesized and the observed profile.The stratification of the atmospheric parameters is given on the optical depth scale at 5000 Å (500 nm), hereafter log τ 500 .The gas pressure and density stratification was computed assuming hydrostatic equilibrium.The atmosphere was perturbed at predefined specific locations of the log τ 500 scale, called nodes, then interpolated to the full depth grid.Providing an initial guess atmospheric model close to the actual solution helps reduce the time it takes for inversion codes to converge.Therefore, to find a good guess value of atmospheric parameters and finalize the node positions, we first inverted the RPs.
We only inverted the Stokes I of the Ca ii 8542 Å and Fe i 6173 Å, along with the Ca ii K spectrum, to infer the stratification of T , V LOS , and V turb .We did not include the other Stokes parameters due to a negligible signal and the selected ROI s having little magnetic activity.Due to the lower cadence of the CRISP data with respect to the CHROMIS data, we assigned a one-fourth weightage to the Ca ii 8542 Å and a one-half weightage to the Fe i 6173 Å compared to the Ca ii K data from CHROMIS while inverting the RPs and the ROI s.The quality of inversion fits are discussed in Appendix B.
To find suitable inversion node positions, we started with the FAL-C (Avrett et al. 1985;Fontenla et al. 1993) model atmosphere as an initial guess for the inversions of the RPs.We started with three equally spaced nodes in T and V turb , and one in V LOS , and proceeded with inversions.
Step by step, we added nodes for those RPs for which we did not find a good fit with the initial setup.After achieving a reasonable fit, we calculated response functions for all the RPs.Response functions contain information on how a spectral profile is sensitive to perturbations in atmospheric parameters at different log τ 500 .Depending on the response functions, we repositioned nodes for different atmospheric parameters, and in cases where a satisfying fit was not achieved, we added more nodes.With this practice, we concluded that different RPs require different number of nodes along with their different positioning along the optical depth scale in order to achieve the best fit.Thus, we classified the RPs into three groups based on the number of nodes required and their position, one group having "quiescent" profiles (30 RPs), one having "emission" profiles (68 RPs), and the third group consisting of two RPs, numbered 18 and 78, which are the strongest CBG RPs.Quiescent RPs are those RPs that do not show any clear K 2 features.Profiles that show either one or both the K 2 features are referred to as emission RPs.Table 1 shows the node positions (in log τ 500 scale) used for different atmospheric parameters for the three categories of RPs.We describe the guess atmospheric parameters inferred from the inversions of the RPs for one time step of ROI A in Fig. 6.There is a temperature rise of about 1 kK (at log τ 500 = −4.2) compared to neighboring pixels and a blueshift of about −2 km s −1 at log τ 500 = −3 at the pixels identified as CBGs.We used the guess atmosphere with the same node positions as the corresponding RPs to invert all the time steps of all the ROI s.The inversions of the selected time sequences for all eight ROI s (78 × 50 × 50 = 195 000 pixels), starting from the guess atmospheres derived above, took 48750 CPU hours using a cluster with each node having two Intel Xeon E5-2683 v4 processors with 16 cores each and a 2.10 GHz base frequency.
For the absolute velocity calibration, we refer to the convective blueshift of the Fe i 6173 Å line at µ = 0.7 measured by Löhner-Böttcher et al. (2019).We calculated a weighted average of the velocities inferred from the RPs in the log τ 500 range of [−1, 0] with weights proportional to the number of occurrences of the RPs in the entire FOV and subtracted it from the convective blueshift to derive the calibration velocity.Finally, we added the calibration velocity (−1.12 km s −1 ) to the V LOS stratification inferred from the inversion of the ROI s.

Evolution of atmospheric parameters
We present the evolution of the atmospheric parameters inferred with the inversions of the eight ROI s that are marked in Fig. 1.We note here that the data presented above were recorded at a viewing angle of µ = 0.7 with the solar normal.The stratification of the atmosphere inferred from inversions may be affected by this geometrical effect because the spatial resolution of the data (0 .038 ∼ 25 km) is extreme compared to the typical formation height of the shocks (∼1000 km; Carlsson & Stein 1997).If one assumes that the shocks propagate vertically upward (Carlsson & Stein 1997;Vecchio et al. 2009), non-shocking atmosphere may have some contribution to the inferred stratification of atmospheric parameters of the CBGs when viewed from a slanted geometry.
The inversion results for ROI A are displayed in Fig. 7. Before the onset of the CBG, the retrieved T map at log τ 500 = −4.2looks homogeneous with minimal variations in T across the ROI and with an average value of ∼3.8 kK.At the onset of the CBG, at t = 41.3 s, a weak compact enhancement in T appears at the center of the ROI .The CBG reaches the maximum area at t = 49.6 s where clear structures of enhancements in T can be seen.At the peak of the CBG, these enhancements  in T are of about 1-2 kK, and these enhancements are found to be maximum at the core of the CBG (inner contour in navy blue color).After t = 66.1 s, the area and T enhancements start decreasing.In contrast, the T map at log τ 500 = −3 shows a reduction in T at the location of the CBG.At the peak of the CBG (t = 49.6 s), this decrease in T is maximum in the range of about 0.5-1 kK.At log τ 500 = −4.2,before the onset of the CBG, the ROI in general shows weak downflow.Between t = 49.6-66.1 s, we see upflows in the CBG at the locations corresponding to higher enhancements in T .The retrieved V LOS exhibits stronger upflows in the CBG at log τ 500 = −3 than at log τ 500 = −4.2.At log τ 500 = −4.2 and −3, the maximum upflow is about −4 km s −1 and −6 km s −1 , respectively.In the photospheric layer, at log τ 500 = −0.1, a typical granulation pattern is visible in the T maps.
At log τ 500 = −4.2, the pixels corresponding to the CBG, generally show less microturbulence compared to the surrounding pixels.Moreover, at log τ 500 = −3 and −0.1, there is almost no microturbulence at the CBG.In contrast, we see constant and consistent microturbulence of 3-5 km s −1 outside of the CBG.One could argue that the T increase inferred in the CBGs atmosphere is an artifact caused by such low values of microturbulence.Inversions are known to have degeneracy between T and V turb , which are competing mechanisms (T and V turb ) that increase the width of the spectral line.Díaz Baso et al. ( 2022) demonstrated that when the line is in emission, the T and the V turb are anticorrelated, which means that increasing the T broadens the spectral line, and thus to maintain the same width, the V turb must be decreased.We note here that the enhancement in T at log τ 500 −4 is necessary to produce intensity enhancement at the K 2V wavelength position; only uncertainty in the amplitude of the enhancement in the T due to such low values of V turb is not very well determined.To better understand the degeneracy between T and V turb , we performed an experiment in which we inverted the brightest pixel of ROI B, but with different values of the V turb , as described in Appendix C. We concluded that as the value of the V turb increases, the intensity at the K 2V wavelength position decreases, causing the fit of the Ca ii K profile to worsen.For higher values of V turb , the inverted profiles show the K 2R peak that is absent in the observed profile.Thus, minimal values of V turb are required to achieve a satisfactory fit of the Ca ii K profile, and the uncertainty in T enhancement due to the degeneracy between T and V turb is not very significant.
The vertical cut before the onset and at the peak of the CBG activity is shown in panel d.The T , compared to the pixels outside the CBG, starts to rise during the CBGs at log τ 500 −3.8, which is co-located with the maximum upflows in V LOS .These upflows are of about −4 km s −1 with a temperature increase of about 1 kK with respect to before the onset of the CBG.The atmosphere before the occurrence of the CBG displays a weak downflow at log τ 500 = −4 and in the atmosphere above.During the CBGs, the magnitude of these downflows show enhancements simultaneous with the enhancements in T and the appearance of upflows between log τ 500 = −3 and −4.At the core of the CBG, no microturbulence is present at any continuum optical depth.
The stratification and evolution of atmospheric parameters of the CBG in ROI B, shown in Fig. 8, are similar to those in ROI A. Before the onset of the CBG, the atmospheres at and above log τ 500 = −4.2show a weak downflow.The start of the enhancement in T relative to surrounding pixels is co-located with upflows at log τ 500 −3.8 (see panel (d) of Fig. 8).There is relatively less microturbulence in the CBG than the surrounding pixels at all optical depths.
The structure and evolution of the atmosphere obtained for the CBG in ROI C are shown in In Fig. 10, we present the variation in the T and V LOS with time in the core of the CBG (also the brightest pixel) from each ROI , A-H.For example the brightest pixel is marked with a "+" in panel a of Fig. 2.During the evolution of a CBG, the T at log τ 500 = −4.2increases and attains a maximum enhancement with respect to the atmosphere before the onset of the CBG, before falling to the values before the onset.This enhancement in T in general is about 1.5-2 kK.The inferred V LOS at log τ 500 = −4.2 at the time of large temperature enhancement in T in the core of the CBGs show the signatures of upflows.For example, the brightest pixel in ROI E and ROI F shows maximum enhancement in T of about 3.5 and 4.5 kK, with upflows of about −5.5 and −3 km s −1 , respectively.However, sometimes these signatures of upflows are missing.For example, in ROI G and ROI H, the maximum enhancement in T is of 2.5 and 1 kK, with downflows of about +0.5 and +2 km s −1 , respectively.The evolution of T at log τ 500 = −3 is opposite to that at log τ 500 = −4.2, and as the CBG evolves, the T decreases, attains minimum, and then recovers to the approximate value before the onset.The retrieved V LOS at log τ 500 = −3 majorly shows signatures of upflows, but sometimes also of weak downflows (ROI C).The upflows at log τ 500 = −3 at the time of peak enhancement in T are stronger than those at log τ 500 = −4.2.

Ca II K response functions
To study the optical depths of formation of the CBGs, we present the response functions of the Ca ii K line with V LOS , V turb , and T in Fig. 11, in the same format as shown in Kianfar et al. (2020).Beckers & Milkey (1975) define the response function (RF) for a physical parameter X as RF X (τ, λ) = δI(λ)/δX(τ).
In the following part of this section, we discuss the response of the Ca ii K Stokes I profile to perturbations in T , V LOS , and, V turb in detail.We have presented RFs for three time steps (rowwise) for the brightest pixel of the CBG in ROI A: before the onset of the CBG, at the onset, and when the CBG is at its maximum.Each plot in the 3×3 grid shows the response of the Stokes I profile to the perturbation in an atmospheric parameter (V LOS , V turb , and T column-wise) as a function of log τ 500 and wavelength.The inferred atmospheric parameters from the inversions of the observed profile (red) and median profile (gray) are also shown for comparison.The observed Stokes I profile (dotted red) and the median profile (dotted gray) are shown in the right column.The synthesized intensities derived from inversions for the observed and the median profile are shown as red and gray dashed curves, respectively.RF T : At t = 33.0s, before the onset of the CBG, the stratification of T retrieved from the multi-line inversions is similar to the T stratification obtained from the median profile.The Ca ii K line within the observed spectral range is sensitive to perturbations in T at log τ 500 between −6 and −2, with line wings forming between log τ 500 = [−4, −2] and line core forming between log τ 500 = [−6, −4].As the CBG evolves, at t = 41.3 s, the K 2V peak has a strong response around log τ 500 = −5, where we see an enhancement of T of ∼0.8 kK compared to that in the atmosphere before the onset of the CBG.The optical depth where the response to T is maximum has decreased from log τ 500 = −5 to at the time (at t = 49.6 s) of strongest enhancement of the K 2V peak intensity.The inferred V LOS near the optical depth with the strongest response (log τ 500 −3.9) shows an upflow of ∼−3 km s −1 .RF V turb : Before the onset of the CBG (at t = 33.0s), the line is sensitive to the perturbations in V turb between log τ 500 = [−6.5,−4.5] with the inferred V turb of ∼3.8 km s −1 , which is higher than that from median profile (∼1.3 km s −1 ).During the onset of the CBG (at t = 41.3 s), the response to perturbations in V turb gets limited between log τ 500 = [−5.5,−4] and the inferred V turb remains nearly unchanged.When the enhancement in intensity of the K 2V peak is the strongest (at t = 49.6 s), the maximum response to perturbations in V turb is restricted about log τ 500 = −4, and the inferred V turb reduces to ∼1 km s −1 above log τ 500 = −4 and vanishes below.
In general, the Ca ii K line is sensitive to perturbations in T between log τ 500 = [−6, −2].However, during the strongest enhancement in the K 2V peak intensity, the sensitivity to perturbations gets restricted to atmospheres below log τ 500 = −4.2,which is also the optical depth of the start of the rise in T .The sensitivity to perturbations in V LOS gets restricted between log τ 500 = [−4.5,−3.5] during the strongest enhancement of the K 2V peak intensity, while before the onset of the CBG with no K 2 features, the V LOS is sensitive to a wide range of optical depths between log τ 500 = [−7, −2].The maximum response to perturbations in V turb during the maximum enhancement of the K 2V peak intensity is limited about log τ 500 = −4.2, with weak sensitivity up to log τ 500 = −5.5, while before the onset of the CBG, the sensitivity to perturbations in V turb is limited between the narrow range of log τ 500 = [−6.5,−4.5].Although we have shown here the response of the Ca ii K line to the perturbations of T , V LOS , and V turb , we included the Fe i 6173 Å and Ca ii 8542 Å lines in inversions, which gives us confidence in the stratification of the inferred atmospheric parameters from the photospheric layers (log τ 500 0) to the chromosphere (log τ 500 −5).

Relationship of velocity and temperature enhancement
We present the scatter and probability distribution plots for maximum T in the CBGs with respect to the T before the onset of the CBGs in panel a of Fig. 12. Panel b of Fig. 12 presents the scatter and probability distribution plots for maximum T enhancement (∆T ) in the CBGs with respect to the inferred V LOS at the time of the maximum enhancement of T .We used the pixels belonging to RPs 18 and 78 in ROI A-H, for example the pixels shown in the middle contours (azure blue) in panel a of Fig. 2, to make the above plots.
At log τ 500 = −4.2,most pixels classified as CBGs have a near-constant maximum T of ∼5.5 kK, irrespective of the T before the onset of the CBG.Few pixels show a large T enhancement with maximum T values reaching up to 8.5 kK.The average quiescent T and enhancement of T is about 4.75 kK and 0.75 kK, respectively.The T at log τ 500 = −3 for pixels in CBG is typically less than the T before the onset of the CBG with an average decrease of −0.75 kK and a maximum reduction in T of −1.5 kK.The T at log τ 500 = −0.1 nearly remains identical to that before the onset of CBGs.
The V LOS retrieved at log τ 500 = −4.2generally show upflows for ∆T values greater than 1 kK and downflows when ∆T is between 0.5-1 kK.The stronger the upflows, the higher the value of ∆T .The inferred V LOS at log τ 500 = −3 typically show upflows, with stronger upflows corresponding to more negative values of ∆T .At the photospheric layers, at log τ 500 = −0.1, the ∆T nearly remains zero, irrespective of the value of V LOS .Panel a shows the scatter plot between T when the CBG is at its maximum brightness (maximum enhancement in the K 2V peak intensity) with respect to T before the onset of the CBG.Panel b shows the scatter plot between the T enhancement with respect to the V LOS when the CBG is at its maximum brightness.

Discussions and conclusions
We have presented the evolution of stratified atmospheric parameters, namely temperature, LOS velocity, and microturbulence from imaging spectropolarimetric observations in CBGs observed in the Ca ii K line with the highest known spatial, spectral, and temporal resolution data.The Ca ii K profiles in the CBGs show an enhanced K 2V peak and a hint of redshift of the K 3 feature.In the center of the CBG structure, when it is at its maximum, we find an average enhancement in temperature at log τ 500 = −4.2 of 1.1 kK and a maximum enhancement of up to ∼4.5 kK.These enhancements in temperature are co-located with upflows in the LOS velocity between log τ 500 = [−4.2,−3].The average strength of these upflows is −2.5 km s −1 and can be as large as −6 km s −1 .However, the extreme values of temperature enhancements and upflows appear only in a few pixels in an acoustic shock region (see Fig. 12).Above log τ 500 = −4.2,we found strong downflows greater than +8 km s −1 , which are stronger than the downflows observed in the quiescent atmospheres prior to the onset of the CBGs.The retrieved value of microturbulence in the atmosphere of CBGs is negligible at chromospheric layers.As explained in Sect.4.1, there could be some degeneracy between the temperature and the microturbulence, which may lead to an overestimation of the values of the temperature enhancement.However, following our experiment described in Appendix C, we conclude that the uncertainty in temperature enhancement is not very significant.An analysis of temperature response functions suggests that during the peak phase of the CBGs, the Ca ii K line is most sensitive to a perturbation in the temperature at log τ 500 = −4.2, in contrast to the quiescent atmosphere, which is sensitive at about log τ 500 = −5.5.We corroborate the findings by da Silva Santos et al. (2020), who found that at the time of shock propagation, using simultaneous inversions of IRIS and ALMA data, the k 2V peak of the Mg ii k line has a predominant contribution from the atmospheric layers that are sensitive around log τ 500 = −4.2.During the CBG activity, we found temperature enhancements with the upflows in the lower chromosphere and downflows in the upper chromosphere in the direction of the LOS.As the CBG progresses, the upper chromosphere (log τ 500 < −4.2) becomes more and more downflowing and the lower chromosphere (log τ 500 = [−4.2,−3]) becomes more and more upflowing.The maximum response to perturbations in temperature and LOS velocity is at the upflowing lower chromospheric layers, giving observational support to the interpretation that CBGs are manifestations of upward-propagating acoustic shock waves in downflowing atmospheres (Carlsson & Stein 1997).The upflows in the lower chromosphere shift the opacity responsible for the two K 2 features symmetrically positioned about the Ca ii K line core to the blue wing of the Ca ii K line.
These upward-propagating shock waves enhance the gas density at lower chromospheric layers, which couples the Ca ii populations to the local conditions.The local temperature enhancements at the lower chromosphere result in an enhancement in the source function in the blue wing of the Ca ii K line (nominal K 2V position).Since the upper chromosphere is downflowing, the overlying opacity is redshifted (i.e., there is little opacity to absorb this blueshifted radiation), resulting in enhanced emission at the K 2V wavelength position of the Ca ii K line (Carlsson & Stein 1997).The downflows above log τ 500 = −4.2suggest that the K 3 opacity is shifted redward of the line core, which causes an opacity removal effect at the K 2V peak, which in turn enhances the K 2V intensity while suppressing the K 2R peak, an effect termed as "Opacity Window" by Bose et al. (2019aBose et al. ( , 2021) ) who studied the formation of the Ca ii K line of the ondisk spicules.However, the response to perturbations in LOS velocity is only significant between log τ 500 = [−4.5,−3.5] and is negligible at higher atmospheric layers; hence, the amplitudes of these downflows are not reliable.The above results underline the role of velocity gradients in the atmosphere in producing excess emission in one of the two (blue or red) wings of spectral lines.Velocity gradients have also been found to be the cause of asymmetric emission in other physical mechanisms.The effect of wave train of upflows and downflows on the line source function has been studied by Scharmer (1981Scharmer ( , 1984)), who suggested that velocity gradients can enhance one of the K 2 peaks.de la Cruz Rodríguez et al. (2015a) found an enhancement in the red wing intensity of the Ca ii 8542 Å line due to velocity gradients produced by an upflowing magnetic bubble against a downflowing background.Steep velocity gradients due to upflows in the lower chromosphere and downflows in the upper chromosphere can give rise to excess emission in blue wing (blue asymmetry) of the Hα spectral line and vice versa (Kuridze et al. 2015).
A153, page 13 of 21 The Ca ii K profiles in the atmosphere surrounding the core of the CBGs also show a relatively weak enhancement in the K 2V peak intensity and redshifted K 3 features.However, we did not find upflows in LOS velocity in the lower chromosphere.The average temperature enhancement seen in such atmospheres is 0.9 kK with maximum enhancement of 2 kK.The enhancement in the K 2V peaks of such profiles can be explained by the "opacity window" effect where the K 3 opacity is redshifted, causing the removal of opacity at the K 2V peak enhancing the feature (Bose et al. 2019a(Bose et al. , 2021)).
We also report that the temperature at log τ 500 = −4.2during acoustic shocks has a near-constant value, irrespective of the temperature before the onset of the shocks.The result can also be interpreted in a way that the temperature enhancement is larger for the plasma regions with cooler quiescent temperatures, which is similar to what Houston et al. (2020) found through the inversions of Ca ii 8542 Å spectra in the case of umbral flashes.
The temperature of CBGs at log τ 500 = −3 decreases with the evolution of the CBG, with the minimum temperature at the peak phase of the CBG.The decrease in temperature is, on average, 0.75 kK.It could be an expression of adiabatic cooling caused by expansion after a shock wave has passed through the region.The velocity flows are typically directed outward from the region's center, causing adiabatic expansion, and thus cooling (Wedemeyer et al. 2004).However, the temperature decrease could also be an artifact of the inversions; the Tikhonov regularization is used to enforce smoothness in the model parameters, and thus to keep the model smooth, a sudden increase in temperature at node position log τ 500 = −4.5 could be followed by a decrease in temperature at the node position log τ 500 = −3.The signature of this temperature decrease is also reflected in the relatively poor fit of the Ca ii K line wings compared to the K 2V emission feature (see Appendix B).Moreover, the response to perturbations in temperature is much stronger in the higher layers (log τ 500 −4.2) than in the lower atmospheric layers, and hence this decrease in temperature is less reliable.
The temperature enhancement in the core of the CBGs (∼1-4.5 kK), at spatial locations with maximum brightness in the bright grains, is up to five-fold from the enhancement reported by Beck et al. (2013), who studied the grains in the Ca ii H line assuming LTE conditions, and two-fold from the enhancement found by Eklund et al. (2020), who studied such brightness temperature enhancements in the millimeter continuum using the data from ALMA.Our value of temperature enhancement is more reliable than the above works because, a) LTE is not a valid approximation to model the Ca ii H&K lines that form in the upper chromosphere and where PRD effects are significant (Milkey & Mihalas 1974;Uitenbroek 1989), and b) the ALMA observations studied by Eklund et al. (2020) had a lower spatial resolution (2 ), which is larger than the size of the CBGs we have studied in this paper.This may lead to reduced brightness temperatures and, therefore, an underestimation of the temperature enhancement, as pointed out by Eklund et al. (2021b).
Umbral flashes (UFs), which are interpreted as manifestations of magneto-acoustic shock waves, also show single-peaked emission about the nominal K 2V position of the Ca ii K line with redshifted K 3 features and an enhancement in the blue wing of the Ca ii 8542 Å line with a redshifted line core.The Doppler shifts of the emitting flows in the UF atmospheres have been discussed by many authors in the literature with both upflowing (de la Cruz Rodríguez et al. 2013;Joshi & de la Cruz Rodríguez 2018;Anan et al. 2019;Houston et al. 2020) and downflowing (Bose et al. 2019b;Henriques et al. 2017;Athay 1970) chromospheres.The above set of authors found an increase in temper-ature of about 1-2 kK compared to the background, which is comparable to our result.
Using a simulation of wave propagation in a sunspot umbra, Felipe et al. (2014) showed that before the onset of the shocks, the chromosphere is predominantly downflowing with typical absorption Ca ii 8542 Å profiles.As the shock progresses, the middle chromosphere becomes more and more upflowing, whereas the upper chromosphere becomes more and more downflowing, with an emission peak visible in blue wing of the Ca ii 8542 Å spectral line.There is a cospatial temperature enhancement about the heights with upflowing plasma, which is about 1.5 kK at the peak of the shock.Henriques et al. (2020) have reported counter-flowing solutions with weak upflows in the lower chromosphere and strong downflows in the upper chromosphere during UFs, through the inversions of spectropolarimetric observations of the Ca ii 8542 Å line.However, they argued that the magnitude of temperature enhancement can be overestimated due to the hydrostatic equilibrium assumption used in the inversion codes.Using simulations it has been demonstrated that upward-propagating shock waves increase the gas density at lower chromospheric layers, resulting in a strong coupling of the source function to the local conditions (Carlsson & Stein 1997).However, due to the requirement of the hydrostatic equilibrium assumption in inversions, such density perturbations cannot be modeled by any of the current inversion codes.As a result, the density perturbations are ignored, and the temperature and the Doppler velocities are the only contributors to the source function enhancement, potentially leading to an overestimation of the magnitude of the temperature enhancement.
In summary, we report the evolution of temperature, LOS velocity, and microturbulence from the photosphere through the lower to the upper chromosphere, through the simultaneous multi-line non-LTE inversions of spectroscopic observations of CBGs in the Ca ii K, Ca ii 8542 Å, and Fe i 6173 Å lines.Our analysis of temperature and LOS velocity response functions support the interpretations using simulations (Carlsson & Stein 1997;Wedemeyer et al. 2004) that CBGs are manifestations of upward-propagating acoustic shocks against a background of downflowing atmospheres.The k-means clustering (MacQueen 1967) is one of the many unsupervised learning algorithms that are used to find patterns and structures in an unlabeled data set.The algorithm works by partitioning an unlabeled data set (m data points with n features) into k clusters.In our case, the unlabeled data set are the spatially resolved image elements and features are 30 wavelength positions, 29 wavelength of which were sampled in the Ca ii K line and one in the continuum 4000 Å position.The algorithm is iterative and its objective is to minimize inertia (σ k ), that is, within cluster sum of squared distances from the cluster center.A cluster center is nothing but the mean of all data points belonging to the cluster.In our case, a cluster center is the mean of all spectral profiles belonging to the cluster, which we call the RP.The algorithm is initialized by k number of predefined cluster centers and each data point is assigned to a cluster with the nearest (euclidean distance, χ 2 ) center.In subsequent iterations, new cluster centers are calculated from the clusters defined in the previous iteration and the process continues until the algorithm converges.The performance of the kmeans clustering algorithm has a heavy dependence on the seed cluster centers.A random selection of seed cluster centers generally results in a poor clustering.Thus, to improve the quality of clustering, we used the k-means++ (Arthur & Vassilvitskii 2007) algorithm, which ensures that seed cluster centers are as far away from each other as possible.The k-means++ algorithm proceeds in an iterative manner by first selecting a cluster center at random, then each new center (up to k) is chosen such that the distance from the nearest previously chosen center is maximum.
Our main objective in performing the k-means algorithm was to identify CBG-like profiles and use the stratification of atmospheric parameters inferred from inversions of RPs as an initial guess atmosphere while inverting actual observations.Hence, the actual spectral profiles must be close to the corresponding RPs at all wavelength positions.The k-means clustering algorithm calculates the euclidean distance between a data point and the cluster mean before assigning it to a cluster, which implies that while finding the nearest cluster, the χ 2 is more affected for the wavelength positions that have higher variance.Hence, we must normalize the data before feeding it to the k-means algo-rithm, which is done by subtracting the mean and dividing by the standard deviation at each wavelength position before doing k-means clustering, ensuring that all wavelength positions have a variance of one.However, we have more wavelength samples between the Ca ii K far wings and the K 1 features (19) compared to the samples between the two K 1 features (10).In addition, we have one wavelength point in the 4000 Å continuum too.Since, the radiation at wavelength positions between the K 1 features originate from higher atmospheric layers (the lower and upper chromosphere) compared to that of at wavelength positions beyond the K 1 features and the continuum position, it will adversely affect the atmospheric parameters retrieved from inversions of RPs at chromospheric layers, which is detrimental to our study.Therefore, we multiplied every pixel with a normalized weight such that the sum of variances at wavelength positions between the K 1 features equals the of sum of the variances at wavelength positions beyond the K 1 features and the continuum point.
Thereafter, seven time-frames under the best seeing conditions were selected for clustering using the k-means algorithm, which contained a total of 15,988,896 pixels.We used the Elbow method to determine the optimal number of clusters k in which our data set could be grouped, where we analyzed change in σ k with respect to k.We performed clustering of the data with 2 ≤ k ≤ 197 and plotted the σ k , and also the change (σ k − σ k+1 ), as shown in the Elbow plot in Fig. A.1.The k is chosen such that σ k decreases linearly.However, in reality, it is generally a smooth curve, and hence we chose a k large enough that σ k −σ k+1 did not change.According to the Elbow plot, the σ k − σ k+1 does not change after k ∼ 70.However, as our reason for using kmeans was to use the guess atmospheres inferred from inversions of the RPs as an initial guess in inversions of the actual observations, we chose a larger value k = 100, to get a better guess atmospheric model.

Appendix B: Quality of fits
We discuss the match between the synthetic and observed Stokes However, in ROI B the observed Ca ii 8542 Å spectral line is seen in emission, suggesting that the observation coincided with the maximum phase of the CBG seen in the Ca ii K spectra, and qualitatively similar to the synthesized Ca ii 8542 Å spectra.This suggests that although the Ca ii 8542 Å line has one-fourth the weightage in inversions compared to the Ca ii K, if the observation in the Ca ii 8542 Å coincides with the CBG, we are able to get good fits with the utilized weighing scheme.

Fig. 1 .Fig. 2 .
Fig. 1.Overview of the quiet-Sun region observed on 6 June 2019.(a) Interface Region Imaging Spectrograph SJI at the Mg ii k 2796 Å.(b) Solar Dynamics Observatory/Helioseismic and Magnetic Imager LOS magnetogram.The white boxes in panels a and b mark the spatial location of the FOV observed with the CHROMIS instrument.(c) Continuum image at 4000 Å.(d) Line-of-sight magnetic field inferred using Milne-Eddington inversions of the Fe i 6173 Å line observed with the CRISP instrument.Panels e and f show images at wavelength offsets of −108.2 and +22.6 mÅ from the Ca ii K line core, respectively.Panels g and h display images at wavelength offsets of −15.8 and +64.2 mÅ from the Ca ii 8542 Å line core, respectively.A gamma correction with γ = 1.5 and γ = 0.5 was applied on panels a and e-h, respectively.The ROI s (A-H) highlighted by colored squares of 50 × 50 pixels (1 .85×1 .85) indicate the locations of the Ca ii K grains that are analyzed in this paper.Region of interest A shows a CBG in Ca ii K narrowband images (panels e and f), the brightenings in other ROI s occurred at different times.

Fig. 4 .
Fig. 4. Same as for Fig. 2, but for ROI C. corresponding RPs.Out of the 100 RPs, we identified four RPs that show clear spectral signatures of CBGs in the Ca ii K line, a single peak emission around the nominal K 2V , and no emission about the K 2R wavelength position (Sect.3.1).We selected these four RPs based on the criterion that the intensity of the K 2V emission peak must be double the intensity at the K 1V wavelength position.The four selected RPs in the Ca ii K line with the corresponding RPs in the Ca ii 8542 Å and the Fe i 6173 Å lines are shown in Fig. 5.The RPs 78, 18, 36, and 85 from the Ca ii K data in decreasing order show the maximum enhancement in the K 2V peak.A hint of emission in the blue wing of the Ca ii 8542 Å line is seen only in RPs 78 and 18, whereas the RPs 85 and 36 are typical absorption profiles.The corresponding Fe i 6173 Å RPs do not show any peculiar spectral signatures.The density distribution of the Ca ii K spectra within a particular cluster is

Fig. 5 .
Fig. 5. Illustrating the efficiency of the partitioning algorithm (k-means) in CBG clusters.The RPs in the Ca ii K, Ca ii 8542 Å, and Fe i 6173 Å lines are shown using solid blue, maroon, and red lines, respectively.The dashed lines mark the zero position in the wavelength axis.The density distribution plots of the Ca ii K spectra for four RPs classified as CBG RPs are shown in blue.The density (darker meaning a higher concentration of spectra) corresponding to each Ca ii K RP shows the distribution of profiles over the entire time series.The percentage of profiles out of the whole data set belonging to a particular cluster is denoted as n.

Fig. 6 .
Fig. 6.Maps of initial guess atmospheres for ROI A based on inversions of corresponding RPs.The first row shows images of the 4000 Å continuum, the Ca ii K line wing image (−0.6 Å), and the line core blue wing image (−0.1 Å).The second, third, and fourth rows show values at log τ 500 = −1, −3, and −4.2 of T , V LOS , and V turb , respectively.The outer contour (cyan) is made with pixels belonging to RPs 85, 36, 18, and 78, while the middle contour (azure blue) is made with RPs 18 and 78.The inner contour (navy blue) shows the region belonging to RP 78, the strongest CBG RP.

Fig. 7 .Fig. 8 .Fig. 9 .
Fig. 7. Evolution of inferred stratified atmospheric parameters from the inversions of ROI A. The parameter maps for three values of log τ 500 = [−4.2,−3, and −0.1] are shown in the panels a, b, and c, respectively.The first, second, and third rows show the evolution of T , V LOS , and V turb .Panel d shows cross cuts along the depth of the atmospheres through the y-axis for two time steps as indicated.The x positions for these cross cuts are indicated by dashed blue lines in panels a-c.The contours in each panel show the pixels classified as CBGs, as shown in Fig. 2.

Fig. 10 .Fig. 11 .
Fig. 10.Time variation of T and V LOS of a brightest pixel of the CBGs in ROI A-H.For example, the brightest pixel is marked by a "+" in panel (a) of Fig. 2. The evolution of atmospheric parameters are shown at log τ 500 = −4.2,−3, and −0.1.The vertical lines indicate the time of maximum T enhancement.

Fig. 12 .
Fig. 12. Relationship of T enhancement and V LOS at log τ 500 = [−4.2,−3, −0.1].Panel a shows the scatter plot between T when the CBG is at its maximum brightness (maximum enhancement in the K 2V peak intensity) with respect to T before the onset of the CBG.Panel b shows the scatter plot between the T enhancement with respect to the V LOS when the CBG is at its maximum brightness.

Fig
Fig. A.1.Finding the optimum number of clusters k for the k-means clustering of Ca ii K profiles.The dark blue curve denotes inertia (σ k ) for k clusters, whereas σ k − σ k+1 is represented by the cyan curve.The vertical dashed line indicates the used number of clusters k = 100.

I
profiles in the Ca ii K, Ca ii 8542 Å, and Fe i 6173 Å lines for ROI A and ROI B in Figs.B.1 and B.2.The overall morphological structure of the CBG is very well reproduced in the narrowband images of the Ca ii K line.A typical granulation structure can be seen in the synthesized Fe i 6173 Å narrowband images, closely resembling that in observed images.The synthesized profiles of the Ca ii K, Ca ii 8542 Å, and Fe i 6173 Å lines show a good fit with the observed ones with a relatively better fitting of the K 2V emission feature compared to the K 1 features and the line wing wavelength positions.The synthesized Ca ii 8542 Å profile of the CBG in ROI A is in emission compared to the observed profile, which shows typical absorption, because the evolution of the CBG is poorly captured in the Ca ii 8542 Å data due to CRISP having ∼ 4.5 times less cadence than CHROMIS.

Fig
Fig. B.1.Illustrating the quality of fits of synthesized narrowband images and spectral profiles inferred using inversions for ROI A. Comparison of fits synthesized with observed narrowband images at wavelength offsets of −108.2 mÅand +22.6 mÅfrom the Ca ii K line core, and −354.9 mÅ from the Fe i 6173 Å line core at different times.A gamma correction was applied with γ = 0.7 on both the observed and synthesized Ca ii K narrowband images before display.The right panel shows the evolution of the observed (dotted) and synthesized (dashed) profiles of a pixel that appear brightest at t = 49.6 s in the narrowband images at the wavelength offset of −108.2 mÅfrom the Ca ii K line core.The Ca ii 8542 Å and the Fe i 6173 Å spectral cubes were acquired 27 s before the Ca ii K cube.
Fig. B.2. Plots illustrating the quality of fits in ROI B in the same format as for Fig. B.1.The Ca ii 8542 Å and the Fe i 6173 Å spectral cubes were acquired 8 s before the Ca ii K cube.

Fig
Fig. D.4.Inversion results for ROI F, in the same format as for Fig. 7.

Table 1 .
Node positions in log τ 500 scale for the temperature, LOS velocity, and microturbulence used for inversions of different groups of RPs.