© The Authors 2026

1. Introduction

The chromosphere is a region of the solar atmosphere that is intrinsically dynamic. It is characterised by short-lived finely structured features such as fibrils, jets, and spicules. The complicated physics governing this region makes the chromosphere challenging to both model and interpret from observations. Processes in the chromosphere are driven by photospheric motions, and the chromosphere acts as a funnel of energy from the interior of the Sun to the corona. Therefore, determining energy transfer and dissipation through the chromosphere is necessary to understand how the hot corona is maintained.

In the chromospheric internetwork, steepening of acoustic waves creates rapid heating seen as bright grains in spectral lines such as Ca II K (Rutten & Uitenbroek 1991). Carlsson & Stein (1992) performed 1D numerical simulations of the acoustic wave propagation with non-LTE radiative transfer of the Ca II K line, finding that waves with periods of 180 s reproduce the bright grains seen in observations. Nonetheless, the heating contribution from these acoustic shock waves is still uncertain (da Silva Santos et al. 2024), and accurate heating estimates are limited by several observational challenges such as spatial resolution (noted in Wedemeyer-Böhm et al. 2007), cadence, and non-LTE effects (Carlsson & Stein 1992).

Improvements in solar instruments over the past decades have furthered our knowledge of the solar atmosphere, and brought new attention to the solar chromosphere (see introduction in Carlsson et al. 2019). Better spatially and spectrally resolved solar observations along with UV spectra of the chromosphere (De Pontieu et al. 2014) have provided us with new opportunities to understand the dynamics and energy transfer in the solar atmosphere. Pushing the limits of solar observations requires a theoretical counterpart in order to understand and explain the energy transfer through the solar atmosphere. Often this theoretical framework is built upon 3D radiative magnetohydrodynamic (rMHD) simulations (e.g. MURaM, Vögler et al. 2005; Bifrost, Gudiksen & Nordlund 2005; and CO5BOLD, Freytag et al. 2012). Three-dimensional rMHD simulations have achieved a high amount of realism in the photosphere (Pereira et al. 2013a), and they can also mimic chromospheric dynamics and processes by including recipes for approximating non-LTE radiative losses (Carlsson & Leenaarts 2012) in codes such as Bifrost, or the more recent chromospheric extension to MURaM (Przybylski et al. 2022).

Using spectral observations to understand energy transfer through the chromosphere is complicated by the non-locality of the radiation formed in this layer. Therefore, studies of chromospheric heating often rely on simplified models such as the hydrostatic atmosphere, for example the famous VALC model (Vernazza et al. 1981) whose radiative losses are often taken as a flux to be balanced by heating processes (e.g. the works of Fossum & Carlsson 2005, 2006; Bello González et al. 2009; Abbasvand et al. 2021), or they rely on non-LTE inversions in a 1.5D fashion (e.g. da Silva Santos et al. 2024). Another strategy is to use rMHD simulations to estimate heating of the chromosphere (Molnar et al. 2023; da Silva Santos et al. 2024). Forward modelling of spectral lines from simulations is computationally cheaper than non-LTE inversions, and it does not suffer from the degeneracy of the source function and linear extinction in the formal solution of radiative transfer that inversions do. By using spectra synthesised from 3D rMHD simulations, we can unambiguously link spectral signatures to atmospheric processes.

For this paper, we synthesised chromospheric spectral lines from 3D rMHD Bifrost simulations, obtaining a comprehensive dataset spanning hundreds of simulation snapshots. The spectra were used to identify signatures of chromospheric heating processes in quiet Sun conditions, and we estimate the total heating contribution coming from locations where we observe these signatures. We also identify the atmospheric properties and heating mechanisms behind these spectra, and study their temporal evolution.

2. Methods

2.1. Simulations and synthetic spectra

To study the formation of chromospheric lines, we employed two 3D simulations produced with the Bifrost code. The simulations in our work cover the convection zone, photosphere, chromosphere, and corona, and we used recipes for non-LTE radiative losses in the chromosphere from Carlsson & Leenaarts (2012). Both simulations were magnetically quiet, with different magnetic field configurations mimicking quiet Sun conditions. One simulation, qs024031, was a typical quiet Sun atmosphere, and it is hereafter referred to as simulation A. The other simulation, ch024031_by200bz005, has a unipolar photospheric magnetic field and vertical coronal magnetic field similar to that of a coronal hole, and it is referred to as simulation B. The simulation boxes were of identical size, with 768³ grid cells, a horizontal grid spacing of 32.25 km, and a mean vertical grid spacing of 21.9 km. At z = 0, simulation B has a mean unsigned magnetic field of 5.8 mT, and simulation A has 4.1 mT. Simulation B has been used in several previous studies (Wedemeyer et al. 2022; Silva et al. 2022; Finley et al. 2022; Brchnelova et al. 2025; Udnæs & Pereira 2025). Simulation A is, as far as we know, has only been used in Udnæs & Pereira (2025).

While both simulations are magnetically quiet and can be thought of as quiet Sun simulations, their different magnetic topologies introduce changes in the atmosphere. Simulation B has a hotter corona and a transition region that occurs lower than for simulation A. The average chromospheric heating is also a few kilowatts per square metre higher in simulation B. These simulations do not have very different chromospheres, but the inclusion of two simulations gives us a broader range of temperatures and densities. Using two different simulations was a choice we made to get a better representation of quiet Sun profiles and explore heating signatures from different magnetic configurations.

The spectra used in this work were synthesised with the RH1.5D code (Uitenbroek 2001; Pereira & Uitenbroek 2015) from a 4 Mm × 4 Mm patch from both simulations. The spectra were calculated with a 10 s cadence for approximately 30 minutes from each simulation: from 133 minutes until 163 minutes for simulation A and from 148 minutes until 181 minutes for simulation B. We performed spectral synthesis for every second grid point in the horizontal directions to reduce the computation time, and because it has a negligible effect on the spatial resolution (in Bifrost, the grid size typically oversamples the real resolution; see discussion in Moe et al. 2022). The spatial sampling of the spectra used in this paper is then 62.5 km × 62.5 km.

We focused on five chromospheric lines: Mg II h&k, Ca II H&K, and Ca II 8542 Å. The lines were synthesised in non-LTE simultaneously using a six-level calcium model atom and a five-level magnesium model atom. The Mg II h&k and Ca II H&K lines were modelled with partial redistribution (PRD), while the infrared Ca II 8542 Å was modelled under the assumption of complete redistribution.

The near-UV lines Mg II h&k are regularly observed by the Interface Region Imaging Spectrograph (IRIS; De Pontieu et al. 2014), and their formation properties are described in a series of IRIS papers (Leenaarts et al. 2013a,b; Pereira et al. 2013b). The lines provide diagnostics of the upper chromosphere, with the full lines tracing heights from the upper photosphere. Upward-propagating shocks are seen in spectrograms of Mg II h (De Pontieu et al. 2015), making this line suitable for our analysis. Synthetic Mg II h&k line profiles from Bifrost simulations are narrower than observations (Pereira et al. 2013b), which Hansteen et al. (2023) suggest could originate from a lack of mass loading into the chromosphere. Ondratschek et al. (2024) used a MURaM-ChE simulation and report synthetic Mg II h&k profiles that are broader than those from Bifrost but still not as broad as observed by IRIS. The authors suggest that insufficient dynamics in the simulations may explain the difference. The specific reason for the narrower synthetic Mg II profiles is of no concern for this work. The mechanisms for line formation and energy deposition are still general, even if our simulations display fewer broad Mg II h&k profiles when compared with observations.

The formation of the Ca II H&K lines is similar to that of the Mg II h&k lines, but Ca is less abundant in the solar atmosphere (Asplund et al. 2009) and thus formed at lower heights. Ca II H&K are in the visible part of the spectrum and therefore observable with ground-based instruments, for example the CHROMIS instrument in the Swedish 1-meter solar telescope (SST; Scharmer et al. 2003). These lines provide good diagnostics for chromospheric heating and show bright grains associated with acoustic shocks (Carlsson & Stein 1992). The Ca II H&K lines are formed at heights where 3D effects could be important for the line cores. Bjørgen et al. (2018) show that PRD effects are more important than 3D effects in the line wings, and 1.5D modelling of the lines under PRD is a feasible option for the wings. Wedemeyer-Böhm & Carlsson (2011) also show that statistical equilibrium is a good approximation for Ca II lines.

The Ca II 8542 Å line is one of calcium’s infrared triplet lines, formed lower than Ca II H&K. This line is a popular choice for inferring the chromospheric magnetic field (e.g. de la Cruz Rodríguez et al. 2012; Quintero Noda et al. 2016; Kuridze et al. 2017), but it also shows the bright grain structure associated with acoustic shock waves (Leenaarts et al. 2009; Martínez González et al. 2023).

2.2. Clustering

We employed clustering to find correlations between synthetic spectra and chromospheric heating in the simulations. This was necessary given the sheer number of spectra and to get a statistical view of chromospheric heating. We used k-means clustering to extract the shapes of the spectra. This algorithm was chosen for its speed and simplicity, making it possible to perform several experiments for different initialisation parameters of the algorithm.

K-means clustering is an unsupervised machine learning technique that is routinely used to analyse spectral data. K-means has been used to analyse both solar observations (e.g. Joshi et al. 2020; Barczynski et al. 2021; Nóbrega-Siverio et al. 2021; Woods et al. 2021; Kleint & Panos 2022) and synthetic spectra from simulations (Moe et al. 2023, 2024; Mathur et al. 2025). The algorithm groups data into k clusters based on similarity, where the similarity is defined by the L₂ norm. Formally, k-means selects clusters by trying to minimise the within-cluster sum of squares,

$$\begin{array}{c}\hfill {\displaystyle \sum _i}||I_{\lambda ,i}-I_{\lambda ,c}{||}^2,\end{array}$$(1)

where I_λ, c is the centroid of the cluster. This centroid is computed as the average spectral profile of a cluster, and is often referred to as the representative profile in spectral line clustering.

In this work, we clustered the Mg II h, Ca II H, and Ca II 8542 Å lines simultaneously, by concatenating the lines. Stacking spectral information from three lines was found to better constrain the underlying atmosphere and resulted in better estimates of chromospheric heating and stratifications. Clustering different lines simultaneously is not a new idea, and has been done before in observational studies (e.g. clustering Stokes components of Fe I lines Khomenko et al. 2003; Pietarila et al. 2007; Viticchié & Sánchez Almeida 2011).

Before clustering the spectra, we performed three important pre-processing steps. We interpolated the spectra into a wavelength grid with constant spacing, normalised the spectra separately for each of the three lines, and performed feature selection in order to favour clusters with emission signatures.

Firstly, we interpolated the spectra into a wavelength grid with constant spacing, to give equal weight to each wavelength bin in the clustering. We restricted the wavelength ranges to ±30 km s⁻¹ around the line core for each line, which discards the continuum and concentrates only on the chromospheric part of the lines.

Secondly, we performed Z-normalisation for each line separately, to give all three lines the same weight in the distance metric in k-means. Otherwise, the clusters would have poor sensitivity to the Mg II h line, which usually is much fainter than the Ca II lines. Z-normalisation preserved spectral line information with respect to the spectral line mean and standard deviation.

Lastly, we filtered out the most common profiles to make the clustering more sensitive to the spectral signatures of heating. In the quiet Sun chromosphere, most Ca II line profiles are in absorption. For our purposes, absorption profiles are less relevant since the strongest heating events usually show emission in Ca II lines. Therefore, we discarded most of the spectra where the Ca II lines were in absorption before clustering the spectra. This was done through an initial clustering with only eight clusters. We selected the clusters where the centroids of the Ca II lines were in absorption, and removed 90% of the profiles from these these clusters. This left us with 20% of the total spectra, where ca. 60% of the profiles had central emission in the Ca II lines. Reducing the amount of absorption profiles allowed us to capture more variation in clusters of emission profiles while using fewer clusters.

3. Results

3.1. Clustering

We clustered the Mg II h, Ca II H, and Ca II 8542 Å lines from the two Bifrost simulations described in Section 2.1. The number of clusters was chosen to minimise both the variation within each cluster and the number of clusters. A good rule of thumb is to use the inflection point of inertia as the number of clusters via the elbow method. In our analysis, there was no clear elbow in the inertia plot (see Fig. B.1), which is common for complex data such as synthetic spectra. We saw a good correspondence between cluster variation and representative profiles for ≳20 clusters, and after reducing the absorption profiles from the dataset (described in Section 2.2), we clustered the data with 25 clusters. After the k-means algorithm converged, we obtained representative profiles of the reduced dataset that gave more weight to emission profiles. Using the representative profiles of the reduced data, we clustered the discarded absorption profiles in order to include the entire dataset of spectra in our analysis. 87.5% of the discarded profiles were placed in clusters where both calcium lines were in pure absorption. Using different initialisations or a larger number of clusters had a negligible effect on our results, and we discuss this in more detail in Appendix B.

Figure 1 shows the representative profiles along with the spectral distribution of each cluster. The clusters are sorted by size, meaning that cluster 0 is the most common and cluster 24 the least common. Cluster 0 has 42% of the spectra, while cluster 24 has only 0.1% of the spectra. The large amount of samples in cluster 0 and very small amount of samples in the last cluster, is both a consequence and a desired effect of the data filtration we performed before clustering the data. In Fig. 1, the most common clusters show typical absorption profiles in the Ca II lines and weak emission in the Mg II h2V and h2R peaks. Many of the less common clusters show strong emission in the core of the calcium spectral lines, which is typically associated with chromospheric heating.

Fig. 1. Simultaneous k-means clustering for the chromospheric Mg II h, Ca II H, and Ca II 8542 Å spectral lines. The clusters, labelled in bold, are ranked from most common in the top to least common in the bottom. The fraction of profiles in each cluster is labelled in the top right of each panel. The spectra are interpolated to a constant velocity separation for Doppler velocity between between ±30 km s−1. Clusters of interest are divided into two classes: blue grains, plotted with blue; and red grains, plotted with red. In the figure, all clusters are normalised separately, but the scaling between the three lines in each cluster is preserved.

We find that spectra with the strongest emission are often associated with strong Doppler shifts and shock waves. Following the naming of Beckers (1964), we call the profiles with the strongest emission by ‘grains’ (blue or red-shifted). Clusters 14, 15, 17, 19, 22, 23, and 24 therefore correspond to bright blue grains, where the Ca II H peak is above 15 nW m⁻² Hz⁻¹ sr⁻¹ and the Ca II 8542 Å peak is stronger than 10 nW m⁻² Hz⁻¹ sr⁻¹. Clusters 16 and 20 correspond to bright red grains with the same peak intensity thresholds as for the bright blue grains. The peak intensity values are summarised in Table 1, and the velocity shift of the emission peaks are summarised in Table 2. Most of the Ca II K lines in strong emission display only one peak, and therefore we only record the strongest peak for this line.

By averaging the atmospheric quantities per columns mass for each cluster, we calculated the ‘typical’ atmosphere for each cluster, and for each simulation. The clusters we were interested in were the clusters with strong emission, i.e. the blue and red grains. The atmospheric profiles for the representative profiles with brightest blue grains are presented in Fig. 2. The top row of this figure shows the average spectra for each cluster. Because we calculated the average spectra and atmospheric quantities separately for the two simulations, the spectra in the first row of the figure are slightly different than the representative profiles shown in Fig. 1. The bottom three rows in Fig. 2 show the average atmospheric profiles for the clusters with high emission in the blue wing. While Fig. 2 shows some differences in atmospheres between the simulations, the trends are the same. Compared to the average, the temperature minima of bright blue grain atmospheres are hotter, occur at higher column masses, and are followed by temperature rises to more than 6000 K, around 2500 K hotter than the mean temperature at the same column mass in the chromosphere. The line-of-sight velocities show a reversal from up-flowing to down-flowing plasma. This is the shock front created by acoustic waves leaking from the convection zone. The shock front creates an excess of chromospheric heating, seen in the bottom row of the figure. In the same panel, the τ_ν = 1 distribution shows that most of the profiles form near the shock-front, where the chromospheric temperature peak is located.

Fig. 2. Atmospheric profiles for the clusters with high emission in the blue wing of the representative profiles. The columns are split into the two different simulations (left: simulation A, right: simulation B). Top row: Average spectra of the clusters. Bottom three rows: Average atmospheric quantities in the clusters. From the top, they are the gas temperature, line-of-sight velocity, and total dissipative heating (Qviscous + Qjoule). The dashed black lines show the averages for the entire simulations. The distribution in the bottom row is the column mass of the blue emission clusters where τν equals unity for the inner blue wing (−5.4 km s−1) of Ca II H.

Figure 3 shows the typical atmospheres for bright red grains. These clusters have similar temperature profiles to the clusters with blue-shifted emission, but different line-of-sight velocities. Many atmospheric profiles in cluster 20 have a down-flow at a low height in the atmosphere, where the red emission peak is formed. This can be seen as a weak negative line-of-sight velocity in Fig. 3, but the average line-of-sight velocity is not corresponding to the velocity shift of the calcium emission peaks (seen in Table 2). The average line-of-sight velocity for cluster 16 does not show any down-flow around the formation height of the line.

Fig. 3. Atmospheric profiles for the clusters with high emission in the red wing of the representative profiles. See caption of Fig. 2. Here, the distribution in the bottom row is the column mass of the red emission clusters where τν equals unity for the inner red wing (5.4 km s−1) of Ca II H.

Cluster 16 shows broad representative profiles with strong emission peaks that are slightly red-shifted in the Ca II lines. The Mg II h representative profile does however have a stronger h2V than a h2R peak. Analysing the individual atmospheres of these profiles, we could not identify a common pattern in how the lines form. As seen in Fig. 3, the height of formation is spread over a large range, and the average line-of-sight velocity is mostly positive around the height of formation for Ca II H. Upon further inspection of their atmospheric stratifications, we found that the average line-of-sight velocity is poorly constrained for both clusters 16 and 20, unlike for the bright blue grains. For cluster 16, the atmospheres fall into two main distributions: a typical atmosphere where there is a chromospheric up-flow of about 7 km s⁻¹, and a less typical atmosphere with a chromospheric down-flow of about 8 km s⁻¹. A break-down of the two different components is shown in Fig. 4. For the up-flowing atmospheres, the red emission peak is formed lower in the atmosphere, at a column mass of around 0.1 kg m⁻². Above the formation of the red peak, there is a strong up-flow which does not absorb the intensity in the red. The blue part of the line is formed much higher in the chromosphere, between 1 × 10⁻³ kg m⁻² and 1 × 10⁴ kg m⁻². For the down-flowing atmospheres, the formation of the red peak forms mainly in the down-flow. For simulation A, only 17% of the profiles in cluster 16 have a down-flow in the chromosphere, while for simulation B this fraction is 31%.

Fig. 4. Two typical atmospheres for cluster 16: one atmosphere with down-flow and one with up-flow in the chromosphere. The blue lines are average quantities for all pixels with a negative line-of-sight velocity at cm = 10−3 and the red lines are average quantities for all pixels with a positive line-of-sight velocity at the same column mass. We show the τ = 1 distributions for the up-flowing (blue) and down-flowing (red) at a Doppler velocity of +5.4 km s−1 for the Ca II H line.

3.2. Heating signatures and temporal evolution

From the clustering we established that the strongest spectral signatures of chromospheric heating are associated with bright grains. We then focused on two objectives: quantifying the proportion of total chromospheric heating associated with bright grains; and characterizing the temporal evolution of spectra arising from the heating events.

To quantify the energy dissipated in the chromosphere, we used the viscous and Joule heating terms from the simulations, Q_ν + Q_η (we refer to the discussion and Appendix E in Noraz et al. 2025 for a discussion of the terms in the energy equation that are dissipative). The total chromospheric heating was calculated by integrating the sum of the volumetric viscous and Joule dissipation terms from the simulation column-by-column through the chromosphere. The result is an energy flux of total dissipative energy from viscous and Joule heating in the chromosphere. To avoid including photospheric heating, we limited the integration to a range in column mass. The lower limit was 0.1 kg m⁻², which was between 600 km to 800 km in height for most columns. The upper limit was 2 × 10⁻⁵ kg m⁻² for simulation B (average height of 3.1 Mm) and 3 × 10⁻⁶ kg m⁻² for simulation A (average height of 5.2 Mm). We chose a higher column-mass limit for simulation B since its transition region appears at a higher column mass (and lower height) than for simulation A. In the simulations, strong heating events represented a large fraction of the total chromospheric heating, with 15% of the columns responsible for around half of the total dissipative heating. Cumulative distributions of chromospheric heating are shown in Fig. A.1 of Appendix A.

Figure 5 shows the distribution of chromospheric heating per cluster. The figure shows a clear correlation between clusters with line core emission and chromospheric heating. Clusters 22, 23, and 24 have the strongest emission peaks for the three lines and are also the clusters where the chromosphere has the highest total heating. Viscous heating is the highest for all clusters, which is expected in a quiet Sun region with weak magnetic fields.

Fig. 5. Median of the integrated chromospheric viscous and Joule heating per cluster. Vertical black lines give the median absolute deviation, i.e. the spread of values inside each cluster. The dashed horizontal lines are the averages of the integrated chromospheric heating (Qν + Qη) for each simulation.

Until now, we have only looked at the spectra as individual pixels. We wanted to study the temporal evolution of chromospheric heating events, which could give additional context necessary to discern heating processes. Figure 6 shows images in Ca II H from snapshot 931 in simulation A. We highlighted the pixels classified as bright blue grains in blue, and pixels classified as bright red grains in red. The blue pixels account for 3.0% of the total area, and the red pixels for 1.3%. The cyan and magenta crosses mark the locations that are plotted as spectrograms in Fig. 7. The spectrograms in Fig. 7 have cluster information plotted to the left. Here, the dark blue and dark red colours correspond to the same clusters as for the top row (i.e. the bright blue and red grains). The lighter colours indicate clusters with weaker emission in the red or blue part of the calcium lines. A light blue colour is given to columns in clusters 3, 6, 7, 8, 11, 12, 18, and 21; and a light red colour to clusters 4, 5, 9, and 13. White means that the profiles belong to one of the clusters where the calcium lines are in absorption.

Fig. 6. Images of the Ca II H blue wing at −5 km s−1 and red wing at +5 km s−1 from simulation A. Overplotted on the images are locations where the pixels belong to a cluster with strong blue emission (cluster 14, 15, 17, 19, 22, 23, and 24) in blue or strong red emission (cluster 16 and 20) in red. The cyan and magenta crosses are placed on locations with especially strong blue or red emission, respectively. The pixel marked with cyan belongs to cluster 24, and the pixel marked with magenta belongs to cluster 20. Spectrograms from these pixels can be seen in Fig. 7.

Fig. 7. Spectrograms of two heating events from simulation A. Top row: bright blue grain from the pixel that is marked with a cyan cross in Fig. 6. Bottom row: bright red grain from the pixel that is marked with a magenta cross in Fig. 6. Both events display an excess of chromospheric heating, shown in the column to the right. The spectrograms are scaled non-linearly with a γ factor to enhance contrast. For Mg II h and Ca II H, γ = 1/4; for Ca II 8542 Å, γ = 1/2. The left inset in the spectrograms shows the cluster class of the spectra at a given time-step: blue and red colours correspond to blue and red grains; light blue and light red colours correspond to clusters where the Ca II lines show weaker emission.

In the Ca II H blue wing of Fig. 6, we see the typical bright blue grains marked with blue. The spectrogram of the bright blue grain that is marked in cyan (first row of Fig. 7) shows an oscillating line core with the typical saw-tooth emission feature of upward-propagating shock waves. To the right, we see the integrated chromospheric heating versus time for the same pixel, which rises to ≈30 kW m⁻² in total.

The spectrogram of the magenta pixel is not associated with a blue grain, neither spatially nor temporally. This is shown in the bottom row of the figure, where we can see the time evolution of the pixel marked with the magenta cross. The Ca II spectra are formed in a region where there is a shear in horizontal flow that creates strong heating that rises to ≈60 kW m⁻² in total.

For the pixel marked with a cyan cross in Fig. 6, the τ = 1 height is lower in the blue part of the line core of Ca II 8542 Å, and higher in the red part of the line core of Ca II 8542 Å. The emission in the blue wing is formed in an upward-propagating shock front, similar to the line-of-sight velocity seen in the cyan line for cluster 24 in Fig. 2. For the spectra marked with the magenta cross, the formation of the red emission peak happens at 0.7 Mm. At this height, there is a down-flow in the line-of-sight velocity, and the heating that creates the emission is produced by a gradient in the horizontal velocity component v_y (shear flow).

To further investigate the spectral signatures of blue and red grains, we looked at a bigger sample of spectrograms from both of the simulations. We show a subset of spectrograms from clusters labelled as blue and red grains in Fig. 8. The figure is representative for the dataset as a whole: pixels labelled as bright blue grains usually show a temporal evolution of an oscillating line core with emission in the blue wing as a shock front propagates into the chromosphere. This signature is usually, but not always, present in all three lines. The figure does not show a temporal coherence for bright red grains, something we did not find in the rest of the spectra either. Some red grains show oscillating line cores, while others have no movement in the line core.

Fig. 8. Evolution of spectrograms in the clusters classified as blue grains (left column) and red grains (right column). The horizontal dashed lines mark the snapshot that correspond to the cluster label given in the upper-left corner of each panel, and the vertical dashed lines mark the line centre. The spectra are scaled non-linearly with a γ factor to enhance contrast. For Mg II h and Ca II H, γ = 1/3; for Ca II 8542 Å, γ = 2/3. All spectra are plotted with a time interval of Δt = 340 s. The colours in the left inset of the spectrograms show the cluster classes similarly to Fig. 7.

Often, pixels were grouped into clusters that were the same or similar as neighbouring pixels. In the colour-coded maps to the left of the spectrograms in Fig. 8, we plotted the classification of the spectra at each timestep. From these maps, we can see that most blue grains are temporally connected to other blue grains, while red grains usually appear alone. We extended this analysis to all spectra in the dataset by grouping the pixels of bright blue or bright red grains together, and calculated the typical lifetimes and areas of these two groups of clusters. The dimensions of these groups were computed by connecting them in both time and space with bi-connectivity, and discarding cluster groups with fewer than two neighbouring pixels in space. The average lifetime and area of these groups are reported in Table 3. The findings support the trend we see in Fig. 8, where blue grains have longer lifetimes than red grains. Blue grains also have bigger areas, which is seen in Fig. 6 where the red grains seem more scattered.

4. Discussion

Using synthetic spectra from two quiet Sun 3D rMHD simulations, we looked for signatures of chromospheric heating in the Mg II h, Ca II H, and Ca II 8542 Å lines. We clustered the spectra together and find that profiles associated with the strongest chromospheric heating have strong emission and are mostly of the type associated with chromospheric bright grains, produced by upward-propagating acoustic shocks, as shown in Fig. 2.

Typical quiet Sun spectral observations show three-minute oscillations from p-mode waves. Transient brightenings in the Ca II lines are signs of heating caused by the steepening of p-mode waves. This internetwork region is sometimes called the clapotissphere (Rutten & Uitenbroek 1991; Rutten 1995) due to the steep shock waves that characterise this area. Chromospheric shock waves also create oscillations in the spectra with sawtooth shaped emission in the blue wing. Therefore, Ca II bright grains have long been a signature of chromospheric wave heating (e.g. Carlsson & Stein 1992).

In the simulations, chromospheric bright grains were relatively rare, but associated with an excess of chromospheric heating. The strong blue emission clusters accounted for only 3.3% of the pixels in simulation A and 2.6% of the pixels in simulation B. However, these pixels had heating rates much larger than the average: four times higher in simulation A and three times higher in simulation B. Bright blue grains dominated the heating locally, and their fraction of the total chromospheric heating balance was as high as 13.4% in simulation A, and 8.1% in simulation B. Using a similar simulation but different methods, Noraz et al. (2025) found that 15% of chromospheric energy dissipation comes from shocks, which is comparable, but slightly higher, than our results. Not all shock-front-heating in our simulations was related to a bright blue grain spectral signature. For instance, some shock front heating happened beneath the height of formation of the Ca II 8542 Å line. In these cases, as the hot plasma was advected upwards in the chromosphere, emission in the Ca II and Mg II lines was seen after one time-step (10 s). Therefore, our estimation of shock-front heating from blue grains is likely to undershoot by a few percent if we compare with the results of Noraz et al. (2025).

Simulation A had about 40% more bright grains than simulation B. However, in both simulations the strong heating was confined to a small area (50% of the heating came from 15% of the area as seen in Fig. A.1). Therefore, a larger proportion of the chromospheric heating in simulation B does not show clear spectral signatures in the Ca II and Mg II lines. Since the main difference between these simulations is the magnetic configuration, this suggests that the inclination of the chromospheric magnetic field determines the occurrence of bright grains.

Other spectral profiles associated with large chromospheric heating have strong emission in the red wing. These were included in two groups of clusters (16 and 20). Unlike the bright blue grain profiles, these profiles did not share a common atmospheric profile nor temporal evolution. Instead, they appear to be formed in a variety of conditions. The profiles in cluster 16 were produced by two different groups of atmospheres. Most of the emission profiles were formed in the lower chromosphere in regions with very low velocities, with an up-flow above that makes the atmosphere transparent to brightenings in the red wing. A smaller fraction of the profiles were actually formed in chromospheric down-flows. Emission in the red wing of the calcium lines can also be caused by magnetic flux emergence (magnetic bubbles, Ortiz et al. 2014; de la Cruz Rodríguez et al. 2015), but in this cluster we saw no conclusive trends in the magnetic field. Ortiz et al. (2014) also explain red emission with enhanced absorption in the blue wing, which is also the case for the up-flowing atmospheres we saw in cluster 16, where the blue wing was formed higher in the chromosphere. The degeneracy in the formation of the red emission peak in calcium lines in clusters 16 and 20 makes it important to use caution when analysing such profiles from observations, as seen by the ambiguity in velocities inferred from Ca II 8542 Å alone (Kuckein et al. 2025). A possible way to resolve the degeneracy of different heating profiles would be to use additional spectral lines, or use different Stokes components to resolve the degeneracy of different heating processes (e.g. as Zhou et al. 2025 have done to identify signatures of reconnection in the 8542 Å line).

The bright red grains comprised only 0.7% of all profiles, but similar to the bright blue grains, their chromospheric heating was much larger than the average (four times higher for simulation A and 3.6 times higher for simulation B). Their contribution to the total chromospheric heating is 3.3% for simulation A and 2.2% for simulation B. This is lower than the contribution of the blue grains because red grains are not as common. If we include all clusters where Ca II H has emission in the line core, we obtain a much larger proportion of the total heating. In simulation A these pixels make up 47% of the total chromospheric heating, and in simulation B the same proportion is 38%.

In absolute values, the spatially and temporally averaged chromospheric heating in our simulation was 3.8 kW m⁻² for simulation A and 6.2 kW m⁻² for simulation B. These values are comparable to the canonical value of radiative losses in the VAL chromosphere of 4.6 kW m⁻² (Vernazza et al. 1981). The average heating from blue and bright grains was substantially higher: 15.2 kW m⁻² and 19.4 kW m⁻² for the blue grains; and 15.2 kW m⁻² and 22.4 kW m⁻² for the red grains. Given that bright grains are rare and their heating is far higher than the average values, it is best to compare them to spatially resolved heating estimates. Díaz Baso et al. (2021) make use of inversions of solar observations to derive spatially resolved estimates of radiative losses in an active region. Averaging over small patches, they find radiative losses of 4.5 kW m⁻² for the quiet Sun, 32 kW m⁻² for a weaker plage, and 108 kW m⁻² for the most active patches. Morosin et al. (2022) use similar methods and find average radiative losses of 28 kW m⁻² for a plage. These plage regions have magnetic fields and conditions much more extreme than our simulated setup, but both Díaz Baso et al. (2021) and Morosin et al. (2022) find many locations farther from the most active parts that locally have radiative losses in the range of 10−20 kW m⁻², fully consistent with our values for the bright grains.

The temperature excursions in shock fronts that produce the bright blue grain signatures were found to be around 2500 K on average. This value is comparable to those inferred from the inversion of SST observations (1000 K to 2000 K in de la Cruz Rodríguez et al. 2013; Joshi & de la Cruz Rodríguez 2018; Mathur et al. 2022), but somewhat larger than values inferred from brightness temperature at millimetre wavelengths (450 K to 750 K, see Eklund et al. 2020). The ALMA observations of Eklund et al. (2020) have a spatial resolution of about $2\stackrel{″}{.}0$, which can very well average small-scale temperature peaks. From Fig. 6 we see that strong emission regions appear at scales smaller than ≈250 km, or about $0\stackrel{″}{.}3$.

The average lifetime of blue emission signatures was almost one minute, consistent with bright grains observed in Ca II 8542 Å (Martínez González et al. 2023) and Ca II H2V (Mathur et al. 2022). The average area of blue grains corresponds to an angular diameter of $0\stackrel{″}{.}3$, and of $0\stackrel{″}{.}2$ for the red grains. The bright grains in our simulations are smaller than previously reported from observations, with Ca II 8542 Å grain sizes seen with IBIS ranging from $0\stackrel{″}{.}7-1\stackrel{″}{.}0$ (Vecchio et al. 2009). This can be due to the higher spatial resolution of the synthetic data.

Our analysis is limited to the heating that occurs in these simulations. Bifrost uses a hyper-diffusion scheme, which makes the code stable for a broad range of different plasma regimes. However, this scheme smooths out gradients and spreads out heating, which should make heating from shock fronts a lower bound (Udnæs & Pereira 2025). Viscosity and resistivity have a large impact on the plasma dynamics. Færder et al. (2024) study the impact of different resistivity models in Bifrost, which is important for the nature of reconnection and heating in the simulations. They find that the hyper-diffusive model currently employed in Bifrost best reproduced the total heating at lower resolutions. Therefore, while our simulations have different resistivities and viscosities than the solar atmosphere, their approximations to chromospheric heating is as close as we can get with current methods.

5. Conclusions

Chromospheric heating is complex and cannot be measured directly. The key spectral diagnostics we can observe form under optically thick non-LTE conditions, whose non-locality can lead to loss of information. In this study, we use 3D rMHD simulations combined with spatio-temporal non-LTE spectral synthesis to understand how to quantify heating (which we can measure directly in the simulations) from spectral signatures.

We calculated chromospheric spectra in non-LTE from two quiet Sun rMHD simulations with different magnetic configurations: one typical quiet Sun simulation (A), and one simulation with a coronal hole-like magnetic field (B). By using k-means clustering, we identified heating signatures in the chromospheric spectra. The spectral signatures associated with the strongest heating showed strong emission on either the blue or red wings of all lines, and we refer to them as blue or red grains, respectively.

We find that blue grains were produced by upward-propagating shock waves, while red grains does not show a typical atmospheric stratification. The blue grains have a coherent temporal evolution with an oscillating line core and sawtooth-like emission in the blue wing. Red grains have many different temporal signatures, including oscillating line cores with red wing emission and stationary line cores with sudden emission in the inner red wing.

The clusters with strong emission accounted for fewer than 4% of all the profiles in the dataset, but these profiles represented 17% of the chromospheric energy dissipation in simulation A, and 10% of the chromospheric energy dissipation in simulation B. There were 40% more bright grains in simulation A compared to simulation B, reflecting a dependence on the magnetic-field orientation, which was unipolar and vertical in simulation B.

Including all clusters with emission in the Ca II H line accounted for 47% and 38% of the heating in simulation A and B, respectively. We find strong emission to be relatively rare in our simulations, but dominating energy dissipation locally. In the simulations, a substantial portion of the chromospheric heating is dynamic with lifetimes under a minute, and localised in regions of a few 1 × 10⁴ km².

Acknowledgments

This work has been supported by the Research Council of Norway through its Centres of Excellence scheme, project number 262622. We kindly acknowledge the computational resources provided by UNINETT Sigma2 – the National Infrastructure for High Performance Computing and Data Storage in Norway.

References

Appendix A: Heating in the simulations

In the simulations we studied, the distribution of dissipative chromospheric heating was highly skewed. We show the cumulative distribution and the contribution to dissipative heating in Fig. A.1. Small areas of the simulations contribute to a large fraction of the total chromospheric heating. E.g. 50% of the total heating comes from only 15% of the total area.

Fig. A.1. Cumulative distributions of chromospheric heating and the contribution to the total chromospheric heating. The dashed black lines show percentiles of the two distributions.

Appendix B: K-means clustering

In this section, we describe the effects of the pruning of clusters. In general, more clusters are favourable in order to capture more variation in line profiles. But, increasing the number of clusters also increases the model’s complexity, bias, and interpretability.

A textbook example to determine the optimal number of clusters is the elbow method. We show inertia as a function of clusters in Fig. B.1, and observe that there does not exist an elbow for our data. Therefore, there is no straightforward way to determine the perfect number of clusters for our data.

Figure B.2 shows how the total fraction of red and blue grains depends on the number of clusters. The figure shows that the fraction of red blue grains stabilise at a large number of clusters, n > 10². The dashed line in the figure indicates 25 clusters, which was the model we use in our results. This model gives similar values for blue grain heating as the models using a much higher number of clusters. Using 300 clusters, we obtain a total heating from red and blue grains at 18.2% for simulation A and 10.9% for simulation B. This is close to the values we got from using 25 clusters: 16.7% for simulation A and 10.4% for simulation B.

We also investigate the atmospheric stratifications with using a higher number of clusters. The degeneracy in the formation of the spectra in cluster 16 shown in Fig. 4 also occurred when we clustered the data with 100 clusters. Figure B.3 shows the line-of-sight velocities from the new cluster that was most similar to cluster 16 in our results. Therefore, the degeneracy in the formation of the slightly red-shifted emission spectra was not an artefact of choosing a small number of clusters.

Fig. B.1. Inertia versus number of clusters of the k-means clustering on the spectra.

Fig. B.2. Results from the k-means algorithm versus number of clusters. Top row: percentage of pixels classified as red or blue grains. Bottom row: percentage of the total heating coming from pixels classified as red and blue grains. The vertical dashed lines indicate the number of clusters in the model used in our results.

Fig. B.3. Line-of-sight velocities from a redshifted cluster in emission, from a k-means model with 100 clusters. Top row: representative profile of the new cluster (violet lines) compared to cluster 16 (green lines) in the results. The solid black lines are the mean spectra. Bottom row: Distribution of line-of-sight velocities from each pixel in the new cluster. In the line-of-sight velocity, there are two common distributions: one with down-flows and one with up-flows. The blue and red dashed lines are averages of the up-flows and down-flows, respectively. The solid black line is the average velocity of the cluster.

All Tables

All Figures

Fig. 1. Simultaneous k-means clustering for the chromospheric Mg II h, Ca II H, and Ca II 8542 Å spectral lines. The clusters, labelled in bold, are ranked from most common in the top to least common in the bottom. The fraction of profiles in each cluster is labelled in the top right of each panel. The spectra are interpolated to a constant velocity separation for Doppler velocity between between ±30 km s−1. Clusters of interest are divided into two classes: blue grains, plotted with blue; and red grains, plotted with red. In the figure, all clusters are normalised separately, but the scaling between the three lines in each cluster is preserved.

In the text

Fig. 2. Atmospheric profiles for the clusters with high emission in the blue wing of the representative profiles. The columns are split into the two different simulations (left: simulation A, right: simulation B). Top row: Average spectra of the clusters. Bottom three rows: Average atmospheric quantities in the clusters. From the top, they are the gas temperature, line-of-sight velocity, and total dissipative heating (Qviscous + Qjoule). The dashed black lines show the averages for the entire simulations. The distribution in the bottom row is the column mass of the blue emission clusters where τν equals unity for the inner blue wing (−5.4 km s−1) of Ca II H.

In the text

	Fig. 3. Atmospheric profiles for the clusters with high emission in the red wing of the representative profiles. See caption of Fig. 2. Here, the distribution in the bottom row is the column mass of the red emission clusters where τν equals unity for the inner red wing (5.4 km s−1) of Ca II H.
In the text

Fig. 4. Two typical atmospheres for cluster 16: one atmosphere with down-flow and one with up-flow in the chromosphere. The blue lines are average quantities for all pixels with a negative line-of-sight velocity at cm = 10−3 and the red lines are average quantities for all pixels with a positive line-of-sight velocity at the same column mass. We show the τ = 1 distributions for the up-flowing (blue) and down-flowing (red) at a Doppler velocity of +5.4 km s−1 for the Ca II H line.

In the text

	Fig. 5. Median of the integrated chromospheric viscous and Joule heating per cluster. Vertical black lines give the median absolute deviation, i.e. the spread of values inside each cluster. The dashed horizontal lines are the averages of the integrated chromospheric heating (Qν + Qη) for each simulation.
In the text

Fig. 6. Images of the Ca II H blue wing at −5 km s−1 and red wing at +5 km s−1 from simulation A. Overplotted on the images are locations where the pixels belong to a cluster with strong blue emission (cluster 14, 15, 17, 19, 22, 23, and 24) in blue or strong red emission (cluster 16 and 20) in red. The cyan and magenta crosses are placed on locations with especially strong blue or red emission, respectively. The pixel marked with cyan belongs to cluster 24, and the pixel marked with magenta belongs to cluster 20. Spectrograms from these pixels can be seen in Fig. 7.

In the text

Fig. 7. Spectrograms of two heating events from simulation A. Top row: bright blue grain from the pixel that is marked with a cyan cross in Fig. 6. Bottom row: bright red grain from the pixel that is marked with a magenta cross in Fig. 6. Both events display an excess of chromospheric heating, shown in the column to the right. The spectrograms are scaled non-linearly with a γ factor to enhance contrast. For Mg II h and Ca II H, γ = 1/4; for Ca II 8542 Å, γ = 1/2. The left inset in the spectrograms shows the cluster class of the spectra at a given time-step: blue and red colours correspond to blue and red grains; light blue and light red colours correspond to clusters where the Ca II lines show weaker emission.

In the text

Fig. 8. Evolution of spectrograms in the clusters classified as blue grains (left column) and red grains (right column). The horizontal dashed lines mark the snapshot that correspond to the cluster label given in the upper-left corner of each panel, and the vertical dashed lines mark the line centre. The spectra are scaled non-linearly with a γ factor to enhance contrast. For Mg II h and Ca II H, γ = 1/3; for Ca II 8542 Å, γ = 2/3. All spectra are plotted with a time interval of Δt = 340 s. The colours in the left inset of the spectrograms show the cluster classes similarly to Fig. 7.

In the text

	Fig. A.1. Cumulative distributions of chromospheric heating and the contribution to the total chromospheric heating. The dashed black lines show percentiles of the two distributions.
In the text

	Fig. B.1. Inertia versus number of clusters of the k-means clustering on the spectra.
In the text

	Fig. B.2. Results from the k-means algorithm versus number of clusters. Top row: percentage of pixels classified as red or blue grains. Bottom row: percentage of the total heating coming from pixels classified as red and blue grains. The vertical dashed lines indicate the number of clusters in the model used in our results.
In the text

Fig. B.3. Line-of-sight velocities from a redshifted cluster in emission, from a k-means model with 100 clusters. Top row: representative profile of the new cluster (violet lines) compared to cluster 16 (green lines) in the results. The solid black lines are the mean spectra. Bottom row: Distribution of line-of-sight velocities from each pixel in the new cluster. In the line-of-sight velocity, there are two common distributions: one with down-flows and one with up-flows. The blue and red dashed lines are averages of the up-flows and down-flows, respectively. The solid black line is the average velocity of the cluster.

In the text