© The Authors 2025

1. Introduction

To understand the hot corona, we must look to the more complex dynamics of the chromosphere, a cooler part of the atmosphere situated closer to the surface. With a higher mass density and greater radiative losses than the corona, the chromosphere needs more energy to be heated; it requires a heating rate approximately an order of magnitude higher than the corona. De Moortel & Browning (2015) suggested that coronal heating is only a side effect of chromospheric heating, so the coronal heating problem is really a chromospheric heating problem. Still, heating of the chromosphere remains a major unresolved issue (Carlsson et al. 2019).

The current understanding of chromospheric heating includes two main processes: heating by wave dissipation and by magnetic reconnection. The magnetic field may also provide heating without reconnection, from current dissipation due to ambipolar diffusion, or field line braiding. It is unclear how these processes play out in the chromosphere since we are limited both by spatial and temporal resolution of solar observations, and spectral data can be complex to interpret. Heating of the solar atmosphere by waves generated in the convection zone was predicted already in the 1940s (Biermann 1946; Schwarzschild 1948; Schatzman 1949), but this theory was partly discredited by high-resolution observations with the transition region and coronal explorer (TRACE; Handy et al. 1999) in work done by Fossum & Carlsson (2005, 2006). Shortly after these studies, however, Wedemeyer-Böhm et al. (2007) used a 3D simulation to argue that while waves may not single-handedly balance radiative losses in the solar atmosphere, they play an important role for the energy balance of the solar atmosphere.

Today, the frontier in solar research lies in the interpretation of solar observations with the help of numerical simulations. This combination allows for a detailed description of the plasma, can disentangle the formation of radiation, and helps us study phenomena that are limited by atmospheric seeing conditions (especially for wave analyses) or that go beyond the diffraction limit of telescopes. 3D radiative magnetohydrodynamic (rMHD) simulations of the solar atmosphere – such as those produced with Bifrost (Gudiksen et al. 2011) – are self-contained; they produce atmospheric waves without periodic drivers at the lower boundary, as is common in 1D experiments of waves (e.g. the works done by Stein & Schwartz 1972; Leibacher et al. 1982, and Carlsson & Stein 1992 to name a few).

3D rMHD simulations can contain a convection zone, chromosphere, and corona, and with self-consistent wave generation are a good tool to study wave propagation through the solar atmosphere. In the convection zone of the simulation, p-mode waves are created. In the photosphere, the p-modes occur as the fast MHD wave, since the local sound speed tends to be much larger than the Alfvén speed. The fast mode is subject to mode conversion and transmission when it crosses the β = 1 layer (Bogdan et al. 2003; Vigeesh et al. 2009) and converts to MHD slow waves, Alfvén waves, or surface refraction, which leads to downward-propagating waves. In the chromosphere, p-modes can steepen into shocks that dissipate energy to heat the plasma in a non-linear way.

According to the 1D hydrostatic model of Vernazza et al. (1981), the radiative losses in the chromosphere are 4.6 kW m⁻². This value is often used as a canonical energy budget that has to be filled by either waves or magnetic energy. Different studies have produced different estimates of acoustic-wave fluxes, which do not always compensate for the canonical value of radiative losses. Fossum & Carlsson (2005) estimated acoustic fluxes around 0.4 kW m⁻² by combining TRACE observations with 1D dynamic simulations from RADYN (Carlsson & Stein 1992). Molnar et al. (2021) also estimated acoustic fluxes of 1 kW m⁻² from RADYN simulations, but combined with Dunn Solar Telescope and Atacama Large Millimetre Array observations. On the other hand, Bello González et al. (2009) and Bello González et al. (2010b) used spectropolarimetric observations from the Vacuum Tower Telescope (see Bello González & Kneer 2008, and references therein) to estimate an acoustic flux in the chromosphere of 4.6 kW m⁻². More recently, Molnar et al. (2023) used observations from the Interface Region Imaging Spectrograph (IRIS, De Pontieu et al. 2014) and inferred acoustic fluxes that are much lower than those in 3D rMHD Bifrost simulations. A problem, as noted by Molnar et al. (2023), is that inferred acoustic fluxes are highly model-dependent. Moreover, the acoustic flux can be underestimated by the limited spatial resolution of observations, especially for studies of the chromosphere in radio wavelengths (Molnar et al. 2021).

Even if acoustic fluxes could be perfectly constrained, it remains unclear what portion of these fluxes in either observations or simulations contributes to atmospheric heating. Such analyses are obscured by acoustic-wave reflection, refraction, or transmission in the chromosphere. Complicating these analyses further is the mode conversion of magneto-acoustic waves at the c_a = c_s layer. In this work, we studied how acoustic-wave flux dissipates to chromospheric heating in simulations of quiet Sun chromospheres. We find the characteristic frequencies and wavenumbers of acoustic waves being dissipated, as well as the typical altitudes of acoustic flux dissipation. Our analysis is limited to vertically propagating waves, an assumption often made in solar physics (see e.g. Bello González et al. 2010a).

2. Simulations

In this work, we investigated acoustic-gravity waves in solar atmospheric simulations produced by the stellar atmosphere code Bifrost (Gudiksen et al. 2011). The simulations provide access to the plasma parameters and we can readily analyse the state of the atmosphere. In our work, we studied two simulations with different magnetic field configurations: a quiet Sun simulation, cb24oi; and a coronal hole simulation, ch024031_by200bz005 (used in Finley et al. 2022; that work contains a comprehensive description of the simulation). Both simulations cover the convection zone, photosphere, chromosphere, and corona, which has horizontal domains of 24 Mm × 24 Mm and extends from 2.5 Mm into the convection zone to the corona 14.3 Mm above the surface. Both simulation boxes have 768³ grid points.

The quiet Sun simulation has a mean unsigned magnetic field at z = 0 of 4.1 mT. It was run for a total of 2 hours and 42 minutes after the simulation is relaxed. The coronal hole simulation has a mean unsigned magnetic field at z = 0 of 5.8 mT, and the simulation is run for 3 hours and 14 minutes after the simulation is relaxed. In the upper atmosphere, the magnetic field of the coronal hole is open and aligned in the z-direction. The magnetic field configurations of the simulations are shown in Fig. 1. The figure shows that not only are the field lines more vertical in the coronal hole, the magnetic field is also stronger in its upper atmosphere.

Fig. 1. Ratio of sound speed to Alfvén speeds in xz cuts taken from quiet Sun and coronal hole simulations. The simulation times are taken at 140 minutes and 40 seconds s and at y coordinates 3.10 Mm. Overplotted are the magnetic field lines, with the thickness of the lines being scaled by the logarithm of the field strength, log10‖B‖.

These simulations are suited for wave studies due to their long duration of more than two hours, which gives a good frequency resolution for Fourier analyses. For both simulations we used 679 snapshots with a cadence of 10 s between each snapshot. The quiet Sun simulation was studied from 50 to 163 minutes, and the coronal hole simulation was studied from 113 to 226 minutes. The Nyquist frequency of the simulations is 50 mHZ, and with a 113-minute time span we have a frequency resolution of ≈0.15 mHZ. This frequency resolution is appropriate for, e.g. resolving individual peaks of the p-mode ridges in a power spectrum.

While MHD simulations are a great tool to study the solar atmosphere, there are limitations coming from assumptions that are made in the models. One artefact that is described in Carlsson et al. (2016) is an unnaturally large amplitude of the global p mode with a period of 450 s. This is also discussed in Molnar et al. (2023), where signals were removed from low frequencies by Fourier filtering. Another challenge of Bifrost simulation is a too-dim UV spectrum (Pereira et al. 2013). The too-weak Mg II UV lines are most likely due to too little mass loading of the chromosphere (Hansteen et al. 2023), which can affect wave propagation in the upper chromosphere. For analyses of shockwaves, we also have limitations from the six-point stencil used by Bifrost, which, along with a quenching operator, suppresses small-amplitude perturbations and large gradients. Therefore, we cannot resolve a shockfront with fewer than six points in space, and shocks manifest as steep gradients instead of discontinuities.

3. Analysis of MHD waves

In MHD theory, three forces can act as restoring forces and transport waves through the medium: gas pressure, buoyancy from gravity, and magnetic tension. If one of them dominates, the MHD wave is a pure one-force mode (respectively, acoustic, gravity, or Alfvénic). Analysing the general three-force mode can be complex, and it is often a good approximation to neglect one component and focus on a two-force mode (Stein & Leibacher 1974).

For the quiet solar chromosphere, we can neglect magnetic tension and approximate MHD waves as an acoustic-gravity two-force mode. Its dispersion relation is thus

$$\begin{array}{c}\hfill k_z^2=c_{\mathrm{s}}^{-2}({\omega }^2-{\omega }_{\mathrm{a}}^2)-({\omega }^2-{\omega }_{\mathrm{g}}^2)\frac{k_{\mathrm{h}}^2}{{\omega }^2},\end{array}$$(1)

where k_z and k_h are the vertical and average horizontal wavenumbers, c_s is the sound speed, and ω_a and ω_g are the acoustic and internal gravity cut-off frequencies:

$$\begin{array}{cc}\hfill {\omega }_{\mathrm{a}}& =\frac{\gamma g}{2c_{\mathrm{s}}},\hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill {\omega }_{\mathrm{g}}& =\frac{\sqrt{\gamma -1}g}{c_{\mathrm{s}}},\hfill \end{array}$$(3)

where γ is the adiabatic index.

We show a dispersion relation diagram (k_h–ν) in Fig. 2 for the special case of an isothermal atmosphere with γ = 5/3. The grey regions represent vertically propagating waves (k_z² > 0), while the white regions represent evanescent waves (k_z² < 0). In a non-isothermal atmosphere, the propagation boundaries are reduced (see Fig. 3 of Vigeesh et al. 2017). The diagram is also affected by radiative cooling (see Fig. 102.1 in Mihalas 1978 for the Newtonian cooling case).

Fig. 2. Diagnostic diagram for acoustic-gravity waves in an isothermal atmosphere. Similar to Mihalas & Mihalas (1984, Fig. 53.1), but with γ = 5/3, as used in our simulations.

In this work, we estimated the energy deposition of waves in the upper atmosphere and considered upward-propagating waves. These could in principle be acoustic or internal gravity waves. From the dispersion diagram in Fig. 2 one sees that only low-frequency gravity waves or higher frequency acoustic waves can propagate. Internal gravity waves are strongly damped by radiation in the lower atmosphere (Mihalas & Toomre 1982) and can be reflected back into the atmosphere beneath β = 1 (Vigeesh et al. 2017) or dissipated by non-linear breaking (Mihalas & Toomre 1981). While internal gravity waves are important for the lower atmosphere, acoustic waves have a greater flux above ∼700 km (Straus et al. 2008). Therefore, we are more interested in the high-frequency acoustic waves’ ability to heat the chromosphere.

Upward-propagating acoustic waves can dissipate energy in the upper atmosphere by shocks or by coupling to other MHD modes, where magnetic tension can no longer be ignored. The critical quantities governing which MHD mode occurs are the sound speed,

$$\begin{array}{c}\hfill c_{\mathrm{s}}=\sqrt{\gamma \frac{p}{\rho }},\end{array}$$

and Alfvén speed,

$$\begin{array}{c}\hfill c_{\mathrm{A}}=\frac{B}{\sqrt{{\mu }_0\rho }},\end{array}$$

where B is the magnetic field strength and μ₀ is the vacuum magnetic permeability. The c_s = c_A layer is the boundary between acoustic waves (c_s > c_A) and other MHD modes. The standard definition of plasma-β,

$$\begin{array}{c}\hfill \beta =\frac{n{k}_{\mathrm{B}}T}{B^2/2{\mu }_0},\end{array}$$(4)

can be written in terms of c_s and c_A:

$$\begin{array}{c}\hfill \beta =\frac{2}{\gamma }\frac{c_{\mathrm{s}}^2}{c_{\mathrm{A}}^2}.\end{array}$$(5)

We see that β = 1 is almost at the same height as c_s = c_A when γ = 5/3. In the lower atmosphere where c_s > c_A, most propagating waves are acoustic. At low plasma-β regions of the chromosphere, when c_A > c_s, the acoustic waves convert to either slow-mode MHD waves that travel with the sound speed, or continue as fast-mode waves traveling faster than the sound speed. Bogdan et al. (2003) showed that it is only at the critical c_s = c_A layer that most coupling between fast and slow MHD waves takes place. When c_A > c_s, most acoustic waves can develop into shocks, which are important for depositing energy in the higher layers.

In Figure 1, we show an illustration of the distinct atmospheric regions by plotting the c_s/c_A ratio from vertical cuts from the numerical simulations we employ. The region of equal sound speed and Alfvén speed is shown in white. The quiet Sun simulation is mainly dominated by pressure in the atmosphere beneath z = 1.5 Mm. The quiet coronal hole simulation has a more magnetically dominated atmosphere, with the c_s = c_A layer closer to 1 Mm on average.

A common tool to analyse waves is the Fourier transform. As part of a discrete Fourier transform (DFT) method, we use the amplitude of the wave with frequency k as the magnitude of the DFT normalised by the number of data points in the input data. To obtain the Fourier power spectrum, we square the amplitude of the DFT and normalise it by the product of the frequency resolution and the wavenumber resolution. This way, we achieve a result that is independent of sampling frequency and spatial resolution.

The Fourier transform is used to calculate the coherence between two signals, which we define as

$$\begin{array}{c}\hfill K_{\mathrm{f},\mathrm{g}}^2=\frac{C_{\mathrm{f},\mathrm{g}}^2+Q_{\mathrm{f},\mathrm{g}}^2}{S_{\mathrm{f},\mathrm{f}}{S}_{\mathrm{g},\mathrm{g}}},\end{array}$$(6)

using the same terminology as Vigeesh et al. (2017). The terms are given as

$$\begin{array}{cc}\hfill C_{\mathrm{f},\mathrm{g}}& =\mathfrak{R}\left\{\mathcal{F}(\mathit{k},\omega )\overline{\mathcal{G}(\mathit{k},\omega )}\right\},\hfill \\ \hfill Q_{\mathrm{f},\mathrm{g}}& =\mathfrak{I}\left\{\mathcal{F}(\mathit{k},\omega )\overline{\mathcal{G}(\mathit{k},\omega )}\right\},\hfill \\ \hfill S_{\mathrm{f},\mathrm{f}}& ={\langle \mathcal{F}(\mathit{k},\omega )\rangle }^2.\hfill \end{array}$$

4. Wave flux

To study the propagation of wave energy, we used the mechanical flux term from Lighthill (1978):

$$\begin{array}{c}\hfill I={\rho }_0{u}^2{c}_{\mathrm{s}},\end{array}$$(7)

where u is the speed of the fluid and ρ is the gas density. By taking the Fourier transform of the velocity, $\widehat{u}=\mathcal{F}(u)$, we obtained the mechanical flux as a function of frequency and wavenumber through the atmosphere:

$$\begin{array}{c}\hfill {\mathcal{F}}_{\mathrm{m}}={\rho }_0{\widehat{u}}^2{c}_{\mathrm{s}}\end{array}$$(8)

where ${\widehat{u}}^2$ is the power spectrum of u. The vertical velocity component u = u_z was used to measure the vertical propagation of the mechanical flux, similarly to Molnar et al. (2023).

Prior to taking the Fourier transform of the vertical velocity signal, we normalised and apodised it. The normalisation was done by subtracting the mean, which removes the zero-frequency component of the Fourier transform. The apodisation was done with a Tukey window (also called a tapered cosine), which mitigates discontinuities at the edges of the signal. The coordinates in time and space were transformed to frequency f with units millihertz, and wavenumber k_x and k_y with units per megametre. The 3D Fourier power of velocity was azimuthally averaged, transforming k_x, k_y → k_h. This is a common technique used to more easily analyse Fourier powers in a k_h − ν diagram (Jess et al. 2023).

In this work, we focused on the interface region around the β = 1 layer. In geometric height, this layer is dynamic and corrugated (see Fig. 1, where c_s = c_A), and therefore we limited the analysis to the heights from the photosphere up to 2 Mm.

We calculated the mechanical flux for horizontal planes of the simulations using Eq. (8) and show it in Fig. 3. The vertical velocity power was multiplied by the average gas density and sound speed at each height. The lines in the figure are the dispersion relation of acoustic-gravity waves described by Eq. (1). In the 1D linear theory outlined previously, the regions above and under these lines are where pressure waves and internal gravity waves propagate. The region between the lines is the evanescent region, where waves do not propagate. We calculated the dispersion relation using the average quantities at each layer. Due to the temperature profile of the chromosphere, the average sound speed increases towards the temperature minimum, before declining again towards the hotter temperatures of the transition region. Therefore, the dispersion relation of acoustic waves has a maximal cut-off frequency at the temperature minimum.

Fig. 3. kh − ν diagram of mechanical flux in quiet Sun and coronal hole simulations at five different heights in the atmosphere. The solid black lines are the dispersion relations calculated from time-averaged quantities at each height. The x-axes are divided by 2π and are therefore not in angular units.

In the k_h − ν diagram, we see strong wave power in the internal gravity wave region in the photosphere. This is expected, and it is the signature of convective motions. Higher up, the p-mode ridges are the only signature present. The velocity power at 0.5 Mm resembles the Fe I power calculated from observations (see, e.g. Fig. 5 of Kneer & Bello González 2011, which shows both p-mode ridges and a convective signature beneath the dispersion relation). The p modes are the box modes of the simulation, and they lie mainly inside the acoustic and internal gravity cut-off frequency (for simpler atmospheric models, this should be completely inside; shown in Pintér & Goossens 1999). The overshooting ridges above the cut-off frequency appear because these 3D rMHD simulations are not an idealised isothermal plane-parallel atmosphere. They are leaking out of the convection zone and propagating through the chromosphere. These leaking p modes above 5.5 mHz are in fact believed to drive dynamic fibrils and some types of spicules (Hansteen et al. 2006; De Pontieu et al. 2007; Heggland et al. 2007).

To obtain the total vertical flux, we integrated the mechanical flux over frequency and wavenumber. We show the total, acoustic, and internal gravity fluxes for both simulations in Fig. 4. To separate acoustic and internal gravity flux, we used the theoretical dispersion relations as a bounding filter. We calculated the acoustic flux by integrating above the dispersion relation, only including the propagating flux, and removing standing waves. This integration also removes the large global oscillation discussed in Carlsson et al. (2016). The internal gravity flux was obtained by integrating in the region under the dispersion relation (labelled internal gravity waves in Fig. 2), and it should be taken as an upper bound since the group velocity of internal gravity waves is lower than the sound speed. Since our 3D simulations do not follow all assumptions used for the dispersion relation in Eq. (1), these fluxes are not totally accurate and should be taken as proxies for the propagating and internal gravity-wave flux.

Fig. 4. Total integrated flux in quiet Sun and coronal hole simulations. Solid lines show the total mechanical flux; dashed lines show integrated acoustic flux; dotted lines represent integrated internal gravity flux. Contributions are determined from the dispersion relation that is calculated in Eq. (1) and shown in Fig. 2.

The flux terms in Fig. 4 are relatively high compared to studies with 1D simulations (e.g. Fossum & Carlsson (2005)). At 0.5 Mm, the acoustic flux is 6 kW m⁻² for both simulations. This is high enough to balance the radiative losses of the chromosphere, at least the 4.6 kW m⁻² estimate of Vernazza et al. (1981). Still, this does not tell us that waves are solely responsible for heating a quiet Sun or coronal hole chromosphere. First, the acoustic flux presented here is a proxy. Second, different physical mechanisms can reduce the flux without leading to dissipation into heat as follows: wave reflection of acoustic waves back to the photosphere, mode conversion into different kinds of waves that travel into the corona, and wave damping by radiation. To heat the chromosphere, pressure waves have to develop into shocks that dissipate wave energy into heat. We analyse the energy deposition by looking at the heating terms in the simulations.

5. Energy dissipation from waves

In our simulations, the wave flux decreases with height (Fig. 4). This can be caused by many different processes. Waves can be reflected (e.g. by the transition region), or they may be refracted and change direction. However, the most relevant behaviour for this study is that waves can dissipate energy into heating of the chromosphere. To investigate the heat deposition by acoustic waves in the chromosphere, we analysed the heating terms in the simulations.

Our Bifrost simulations have two main dissipative processes that heat the atmosphere: viscous and Ohmic (or Joule) dissipation (Gudiksen & Nordlund 2005). Viscous dissipation is the mechanism converting wave flux to heat when acoustic waves develop into shocks, and it dominates Joule heating by about an order of magnitude in the mid chromospheres of the simulations. Our goal is to identify when viscous heating is brought on by waves and analyse what wave characteristics contribute to the heating.

If we take the height derivative of the acoustic flux, F_ac, we end up with the amount of acoustic-wave energy lost in the upward direction of the atmosphere. Energy is lost from reflection or conversion from wave energy to heat, so the acoustic-wave heating should be smaller than this gradient:

$$\begin{array}{c}\hfill \frac{\partial F_{\mathit{ac}}}{\partial z}>\text{acoustic wave heating}.\end{array}$$

If high-frequency acoustic waves are the dominant source of viscous heating in the simulated chromospheres, then the equation above must be true.

We calculated the derivative of the acoustic flux shown in Fig. 4, and compared it with the average viscous heating per height in the simulations in Fig. 5. For intermediate heights (approx. 0.7–1.7 Mm), the acoustic-wave dissipation is greater than the viscous dissipation for both of the simulations. Around 0.6 Mm, the gradient of acoustic flux is almost zero, which means the acoustic flux is not damped around the temperature minimum. The coronal hole simulation has a smaller loss of flux above 1.5 Mm, before a spike appears at 2 Mm, which is caused by a large gradient in sound speed in the transition region. This simulation has greater viscous dissipation higher up, which leads to larger flux. However, both the acoustic flux and its gradient decline rapidly above 2 Mm, to the point where the acoustic flux is dropping below the viscous dissipation, meaning that other processes are driving the viscous dissipation. The quiet Sun simulation has a smaller acoustic flux gradient, but it remains above the viscous dissipation at higher layers.

Fig. 5. Gradient of vertical acoustic flux (solid lines) and time-averaged viscous dissipation (dashed lines) in simulations of the quiet Sun (blue) and coronal hole (red).

From Fig. 5, we can see that it is plausible that acoustic waves evolving into shockwaves are responsible for most of the viscous dissipation heating in the mid chromospheres of the quiet Sun and coronal hole simulations. A Fourier analysis of the viscous heating term, ℱ_q, also shows high power in the p-mode ridges observed in ℱ_m (see Fig. 3). While we do not show ℱ_q here, we related the viscous heating to the mechanical flux by computing the coherence between ℱ_m and ℱ_q at each height using Eq. (6). Calculating the 80% coherence threshold, which is common in solar observations (Jess et al. 2023), we obtain the characteristic frequencies and wavenumbers of wave heating, which we show in the k_h − ν diagram of Fig. 6. The figure shows the coherence above the 80% confidence limit, which happens for intermediate heights only (0.75 Mm to 1.25 Mm).

Fig. 6. kh − ν diagram of coherence Km, q between mechanical flux ℱm and viscous dissipation ℱq. Km, q, given by Eq. (6), is calculated at five different heights in the simulated chromospheres. The solid black lines are the dispersion relations calculated from time-averaged quantities at each height. The red contour line outlines the 80% confidence interval.

Acoustic-gravity waves that contribute to viscous heating have relatively high frequencies above 4 mHz, but there are no signs of acoustic-wave heating above 8 mHz. None of the coherence is inside the region of internal gravity waves, which is expected. However, there are regions with high coherence below the acoustic cut-off frequency. This is likely because we are using a dispersion relation for an isothermal atmosphere, while our 3D model has local variations in temperature (e.g. cool pockets of gas in the chromosphere with significantly lower sound speeds).

6. A qualitative view of heating in a wavefront

So far, we have looked at average wave energy deposition across hundreds of simulation snapshots. We also wanted to know, however, what happens at the level of individual events. Here, we take a detailed look at a single event where a vertically propagating acoustic wave evolved into a strong shockfront.

We selected the pixel at (x = 3.25 Mm, y = 3.06 Mm) from the coronal hole simulation and identified several vertically propagating wavefronts that developed into shocks in the chromosphere. In Fig. 7, we show the column’s vertical velocity, temperature, and dissipation coefficients for 15 mins of simulated time. Wavefronts can be seen in all panels of the plot, but they are most identifiable in the vertical velocity. In the vertical velocity plot, we see gradients that propagate from the photosphere into the chromosphere and often reflect downwards off the transition region. These waves have typical periods of 3–4 min and come from the p modes of the simulation.

Fig. 7. Vertical velocity, temperature, and dissipation coefficients in a column of the coronal hole simulation at x = 3.25 Mm and y = 3.06 Mm.

We analysed a shock that developed at t = 149 min in depth; this is visualised in Fig. 8. The shockfront is shown in vertical velocity in panel a of Fig. 8, with the location of the shockfront depicted by the solid black line. The dashed lines trace a particle propagating upwards with the local fast-mode speed or with the local sound speed. When the shock passes through the chromosphere, the gas density decreases, which leads to two things. First, the wave increases in amplitude due to energy conservation. The increased amplitude makes the wave develop into a strong shock and leads to wave heating, which can be seen in panels c and e. Second, the magnetic pressure becomes comparable to the plasma pressure (β ≈ 1). Therefore, the magnetic field becomes important in wave propagation. We see this in panel a: the wavefront is supersonic; its speed is very close to the fast-mode speed from 0.9 Mm until 1.4 Mm.

Fig. 8. A column of the coronal hole simulation where a wave front develops into a shock. x-axes show simulation time in minutes, and y-axes show height above the surface in Mm. Panel (a) shows vertical velocity, (b) shows temperature, and (e) and (f) show the viscous and Ohmic (or Joule) volumetric dissipation coefficients. The lines of panels (c) and (d) are the volumetric dissipation coefficients integrated over the 5 grid points constituting the computational stencil of the shock front. The solid lines in panels (a), (b), (c), and (d) show the location of the shock front. The dashed and dash-dotted magenta lines in panel (a) traces particles travelling with the local fast-mode and sound speed cs respectively.

From the time when the shock starts to form at 149 min and 10 s, the dissipation coefficients increase as the shock travels upwards and becomes stronger. After 40 s, the integrated viscous dissipation reaches its peak (Fig. 8 panel c) with a heating rate of 4 kW m⁻², before the shock decreases. When the shockfront passes z = 1.4 Mm, it reflects downwards towards the photosphere.

In this example, as seen in panels c and d of Fig. 8, both viscous and Ohmic dissipation are enhanced in the shock, but viscous dissipation is much higher. This is expected from friction arising in the shockfront. Integrating over the shockfront, we find that the maximum viscous dissipation is about 6.3 times higher than the Ohmic dissipation. We find a similar scenario by looking at other individual events (not shown), where viscous dissipation often culminates above 10 kW m⁻². The heating in these events is of the same order of magnitude as radiative losses inferred from quiet Sun observations: 4.5 kW m⁻² in Díaz Baso et al. (2021) or 2.8 kW m⁻² in da Silva Santos et al. (2024). The detailed view from individual events is therefore consistent with the average flux results of the previous section, in particular those shown in Fig. 5, and it shows that shocks can deposit significant amounts of energy in the chromosphere. To find out exactly how much energy is deposited by all shocks, we would require a systematic accounting of shocks over all locations, which is beyond the scope of this work.

7. Discussion

We studied wave fluxes and energy deposition in 3D rMHD simulations, both from a global perspective and a detailed perspective of individual shocks. Estimating acoustic fluxes and wave dissipation from these complex numerical experiments is fraught with uncertainty, and so we adopted several assumptions to make the problem more tractable.

We adopted the dispersion relation of an isothermal atmosphere, while the simulations are dynamic and inhomogeneous. We approximated the wave fluxes in Eq. (8) by using the horizontal averages of the gas density and sound speed. Using our dispersion relation implies adiabatic wave propagation, an assumption that can be troublesome in the presence of vertical temperature gradients, which increase the cut-off frequency (Vigeesh et al. 2017). However, the effects of radiation counteract vertical temperature gradients, pushing the propagation closer to the adiabatic case (Mihalas 1978; Khomenko et al. 2008). For example, using the Newtonian radiation field described in Souffrin (1966) and Mihalas (1978), with solar conditions at the temperature minimum (T ≈ 2700 K and γ = 5/3), the acoustic cut-off frequency for an instantly relaxed radiation field is 3.6 mHz, while an adiabatic atmosphere has an acoustic cut-off frequency of 4.6 mHz (a Newtonian radiation field mimics optically thin radiation). With optically thick losses, radiation is not effectively altering the dispersion relation of Heggland et al. (2011). However, inclined magnetic fields allow the propagation of low-frequency acoustic waves in the chromosphere (De Pontieu et al. 2004; Hansteen et al. 2006; Jefferies et al. 2006). Therefore, an accurate dispersion relation for wave propagation is not feasible, and the adiabatic propagation is a fair approximation for our purposes. Propagation at lower frequencies will happen above inclined magnetic field concentrations, and we reiterate that our flux estimates are a lower bound.

Wave flux at high frequencies can differ significantly between 3D rMHD simulations run with different codes by up to two orders of magnitude (Fleck et al. 2021). Therefore, high-frequency flux should be approached with care. However, in our case its contribution to the total flux is small because most of the chromospheric flux is concentrated around 5 mHz.

We find that internal gravity waves have very low power in the upper photosphere and chromosphere and are effectively damped in these higher layers, which is consistent with results of earlier works (e.g. Mihalas & Toomre 1981, 1982; Vigeesh et al. 2017). Non-linear wave breaking of internal gravity waves (Mihalas & Toomre 1981; Vigeesh et al. 2017) is dominant in the lower and middle chromospheres, turning propagating internal waves into turbulence. At the height of wave breaking, the maximal vertical velocity of internal waves is a few kilometres per second; that is, much lower than the group velocity of acoustic waves (Mihalas & Toomre 1981). Damping at lower heights is shown in Mihalas & Toomre (1982), which identified radiative transfer effects that smooth the temperature variations from internal gravity waves, damping them from the low photosphere to the temperature minimum. The suppression of internal gravity waves shifts the peak of the acoustic-gravity flux to higher frequencies, making high-frequency waves more relevant for chromospheric heating.

In Fig. 5, we show that internal gravity waves are not as important as acoustic waves for chromospheric heating. The acoustic and gravity-wave-heating regions correspond to the high-frequency part of the p-mode ridges seen in Fig. 3. These leaked p modes come from the convection zones of the simulations, travel vertically, and develop into shocks that dissipate wave energy in the chromosphere. Since most of the coherence is located at heights of around 1 Mm (Fig. 6), acoustic-wave heating takes place in the middle of the chromosphere. Above this region, there may still be wave heating, but it will be from waves that have the magnetic field as a restoring force; that is, from mode conversion.

We see maximal coherence between acoustic flux and viscous heating around 1 Mm, where the simulations typically have the c_s = c_A divide that is seen in Figure 1. This layer is also where fast and slow magneto-acoustic-gravity waves decouple (Bogdan et al. 2003). Under this layer, where pressure dominates, the fast and slow magneto-acoustic-gravity waves are essentially the fast-mode acoustic wave. The reason for the high coherence between acoustic flux and heating at this layer is two-fold. First, the acoustic waves develop into shocks at this height, which is essential to producing viscous heating. Second, above this height, we have mode conversion of the acoustic waves. The coronal hole atmosphere has a stronger and more uniform magnetic field, which introduces magnetic field dependency to the acoustic waves lower down in the atmosphere. Therefore, the acoustic heating coherence shown in Fig. 6 is smaller for the coronal hole simulation. Other reasons for lower acoustic-wave heating above the c_s = c_A layer include the dissipation and reflection of acoustic shock waves.

In Fig. 4, we see a rapid decline in acoustic flux around 1 Mm, which is also seen as a peak in the gradient of acoustic flux at 1 Mm in Fig. 5. Some of the decline in acoustic flux is due to mode conversion, but from Fig. 5 we also see viscous dissipation that is comparable to the loss of acoustic flux. Looking at an individual shockfront, we find mode conversion from an acoustic shock to a fast magneto-acoustic shock above 1 Mm, while the shockfront produced high heating rates in the chromosphere before reflecting off the transition region (Fig. 8). Reflection of the fast shock can be due to the steep gradient in the Alfvén speed (Khomenko & Cally 2012).

It is difficult to study the viscous dissipation power ℱ_q quantitatively. Because viscous dissipation is not a smooth function, its Fourier transform shows Gibbs phenomenon (Wilbraham 1848) artefacts. However, we saw large regions with more than 80% coherence between acoustic flux and viscous dissipation, which we show in Fig. 6. We then used the coherence between acoustic-gravity flux and viscous dissipation as a proxy to identify the locations and frequencies where acoustic-wave heating takes place. We find that acoustic-wave heating is most dominant in the region where the plasma transitions from being pressure-dominated to magnetically dominated. The simulations also show acoustic-wave heating centred around the acoustic cut-off frequency at spatial scales bigger than 1 Mm. The coronal hole simulation displays a lower acoustic-wave-heating coherence than the quiet Sun simulation, which indicates that wave propagation already depends on the magnetic field at heights of 1 Mm. Magnetic field effects on a shockfront can also be seen in Fig. 8. Here, the shock propagated faster than the sound speed above a height of 1 Mm, where it propagated as a fast magneto-acoustic shock. Above 1.4 Mm, the shock was reflected from the transition region towards the lower atmosphere.

Our analysis rests on two main limitations: how accurately we can estimate shock heating from simulations, and the fact that we analyse only vertically propagating waves.

Sharp shockfronts are predicted to have a width of only a few mean-free paths (Priest 2014). A 3D simulation is not realistically able to resolve the shockfront, because grid sizes are not as small and the numerical stencils need to span several grid points. In a 3D simulation, we see shocks as steep gradients with a width comparable to the numerical stencil instead of being a jump over an infinitesimal thickness. In Bifrost, the numerical diffusion (diffusive and quench operator; see Gudiksen et al. 2011) spreads the immediate heating from the shockfront out over five-to-six grid points. When we recovered the heating from the shockfront, we therefore counted contributions from the grid points that constituted the numerical stencil around the shockfront. If anything, this leads to a smoothing of the estimated heat deposition, and our values should be a lower bound, which has a minimal effect on the results.

While we chose to only study vertically propagating acoustic waves to make the analysis feasible, it is also a reasonable assumption to make when studying energy transport through the chromosphere. Acoustic waves have smaller horizontal amplitudes than vertical ones in the lower atmosphere (Bello González et al. 2010a), and horizontal amplitudes are more important for internal gravity waves (Mihalas & Toomre 1982), which do not transport much flux into the chromosphere. Fossum & Carlsson (2006) discussed that a spherical wavefront that is excited in the deeper atmosphere will be refracted to a more planar wave due to the change in the sound speed, and it will travel in the direction of the gradient of the sound speed, which for most cases is in the vertical direction. Therefore, we believe that only accounting for vertically propagating waves means that our values should be a lower bound, and so this does not affect our results.

Studies using 1D simulations (Fossum & Carlsson 2005, 2006) and 1D hydrostatic models (Sobotka et al. 2016) combined with observations have found that acoustic-wave heating is not high enough to balance radiative losses. Still, these studies discuss that limitations in the spatial resolution can negatively affect flux calculations, a point that is reinforced by findings in Yadav et al. (2021), which found that current observations give flux estimates roughly a factor of three too small. Sobotka et al. (2016) concluded that acoustic waves may still may become the main contributor to chromospheric heating in certain chromospheric regions. da Silva Santos et al. (2024) also found that non-LTE inversions can underestimate acoustic fluxes in the chromosphere and overestimate radiative losses. Here, we did not calculate the proportion of acoustic flux that is dissipated to heat in the chromosphere, and therefore we make no claim about the proportion of heating that comes from acoustic waves in the quiet chromosphere. However, we see that individual events provide substantial heating where acoustic-wave heating dominates. This suggests that acoustic-wave heating must be an important component in the energy balance of the quiet chromosphere.

8. Conclusions

In this work, we looked into the role of acoustic-wave heating in simulations of quiet chromospheres. Analysing two rMHD simulations with quiet magnetic fields, we obtained estimates of acoustic flux and wave-heating characteristics in the photosphere and chromosphere of the simulations. One simulation had a magnetic field configuration typical of the quiet Sun, and the other had that of a coronal hole. The open magnetic field of the coronal hole gave slightly different wave heating characteristics than for the quiet Sun.

Using a simple acoustic-gravity dispersion relation, we separated the acoustic flux into two components: acoustic and internal gravity waves. This approach worked well for the lower atmosphere of the simulations, where the plasma is pressure-dominated, and the sound speed is higher than the Alfvén speed. In the chromospheres we find that most wave flux is acoustic; the internal-gravity flux is quickly damped in the photosphere and is weaker in the chromosphere (above 0.5 Mm). k_h − ν diagrams show strong chromospheric flux around the p-mode ridges of the simulations, with frequencies from 3 mHz to 6 mHz and horizontal wavenumbers from 0 Mm to 1 Mm⁻¹. While most of these power ridges are underneath the acoustic cut-off frequency, in the evanescent region, we find that a large fraction of the flux is above the acoustic cut-off frequency in the simulated chromospheres.

The heating from acoustic-gravity waves is centred mainly at heights of around 1 Mm for both simulations, which is the region where, on average, plasma-β reaches unity. Acoustic heating is more pronounced in the quiet Sun simulation. In the coronal hole simulation, acoustic heating is less distinct due to the influence of the magnetic field on the chromospheric waves. A qualitative analysis showed a vertically propagating shockfront converted into a fast magneto-acoustic shock in the coronal hole chromosphere, before being reflected off the transition region. Together, these analyses highlight the importance of the magnetic field configuration on the propagation of acoustic-gravity waves across the chromosphere: a weak magnetic field facilitates the propagation of acoustic waves, and these waves dissipate most of their energy in the region where plasma pressure balances the magnetic pressure.

Finally, we compared the acoustic flux to the average viscous dissipation at different heights. At heights above 0.75 Mm, we find that the energy lost by acoustic waves is more than enough to account for the heat generated by viscous dissipation of waves. This view is complemented by the analysis of single shock events, which shows that acoustic waves are the largest contributors to viscous dissipation in the simulated chromospheres. We conclude that acoustic-wave heating is an important mechanism in the chromospheric energy balance for the quiet Sun and coronal holes. Further work is needed to identify and measure the prevalence of shocks, so one can estimate the total energy deposited by acoustic waves and relate it to radiative losses and ultimately the heating of the chromosphere.

Acknowledgments

This work has been supported by the Research Council of Norway through its Centers 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

All Figures

	Fig. 1. Ratio of sound speed to Alfvén speeds in xz cuts taken from quiet Sun and coronal hole simulations. The simulation times are taken at 140 minutes and 40 seconds s and at y coordinates 3.10 Mm. Overplotted are the magnetic field lines, with the thickness of the lines being scaled by the logarithm of the field strength, log10‖B‖.
In the text

	Fig. 2. Diagnostic diagram for acoustic-gravity waves in an isothermal atmosphere. Similar to Mihalas & Mihalas (1984, Fig. 53.1), but with γ = 5/3, as used in our simulations.
In the text

	Fig. 3. kh − ν diagram of mechanical flux in quiet Sun and coronal hole simulations at five different heights in the atmosphere. The solid black lines are the dispersion relations calculated from time-averaged quantities at each height. The x-axes are divided by 2π and are therefore not in angular units.
In the text

	Fig. 4. Total integrated flux in quiet Sun and coronal hole simulations. Solid lines show the total mechanical flux; dashed lines show integrated acoustic flux; dotted lines represent integrated internal gravity flux. Contributions are determined from the dispersion relation that is calculated in Eq. (1) and shown in Fig. 2.
In the text

	Fig. 5. Gradient of vertical acoustic flux (solid lines) and time-averaged viscous dissipation (dashed lines) in simulations of the quiet Sun (blue) and coronal hole (red).
In the text

	Fig. 6. kh − ν diagram of coherence Km, q between mechanical flux ℱm and viscous dissipation ℱq. Km, q, given by Eq. (6), is calculated at five different heights in the simulated chromospheres. The solid black lines are the dispersion relations calculated from time-averaged quantities at each height. The red contour line outlines the 80% confidence interval.
In the text

	Fig. 7. Vertical velocity, temperature, and dissipation coefficients in a column of the coronal hole simulation at x = 3.25 Mm and y = 3.06 Mm.
In the text

Fig. 8. A column of the coronal hole simulation where a wave front develops into a shock. x-axes show simulation time in minutes, and y-axes show height above the surface in Mm. Panel (a) shows vertical velocity, (b) shows temperature, and (e) and (f) show the viscous and Ohmic (or Joule) volumetric dissipation coefficients. The lines of panels (c) and (d) are the volumetric dissipation coefficients integrated over the 5 grid points constituting the computational stencil of the shock front. The solid lines in panels (a), (b), (c), and (d) show the location of the shock front. The dashed and dash-dotted magenta lines in panel (a) traces particles travelling with the local fast-mode and sound speed cs respectively.

In the text