© The Authors 2025

1. Introduction

The structure and dynamics of the solar chromosphere is determined by the local magnetic environment. The magnetic field couples the solar interior to the atmosphere, controlling the supply of energy and mass from turbulent near-surface convection. Many mechanisms have been proposed for the heating of the chromosphere, including magnetoacoustic and Alfvénic waves, spicules and swirls, tornados, and reconnection driven eruptive events and jets. A large fraction of the solar surface is covered by the internetwork, i.e. the space between the magnetic network in the quiet Sun. This paper deals with the internetwork chromosphere.

Observations of the internetwork chromosphere show clear signatures of acoustic shocks, seen as the ‘bright grains’ in CaII H&K lines (Rutten & Uitenbroek 1991; Carlsson & Stein 1997). Localised brightenings are also seen where small-scale magnetic loops reach into the chromosphere (Martínez González & Bellot Rubio 2009). Internetwork fields cover a large fraction of the solar surface, and provide a significant amount of the magnetic flux at the solar surface (Gošić et al. 2016).

At the photosphere, magnetic field measurements utilising Zeeman and Hanle diagnostics reveal ubiquitous small-scale fields (de Wijn et al. 2009). In radiation-magnetohydrodynamic (rMHD) simulations of the solar convection zone, these quiet Sun fields are generated by a small-scale dynamo (SSD) mechanism (Vögler & Schüssler 2007; Pietarila Graham et al. 2010; Moll et al. 2011). Simulations including SSD fields produce a close match to magnetic field strengths inferred from observations with the Hinode spectropolarimeter (Lites 2011; Danilovic et al. 2010, 2016) and Hanle measurements (Trujillo Bueno et al. 2004; del Pino Alemán et al. 2018; Zeuner et al. 2024). The temporal variation of the IN fields is also consistent with an SSD process (Buehler et al. 2013). This dynamo-process is not restricted to the near-surface layers, and is present throughout the convection zone (Hotta et al. 2015; Hotta & Kusano 2021). Recirculation in the deeper convection zone is required to match observationally inferred magnetic field strengths (Rempel 2014). The large Reynolds and magnetic Reynolds number of the Sun put simulations with realistic diffusivities beyond current computational capabilities. Brandenburg & Rempel (2019) compared direct numerical simulations (DNSs), with explicit viscous and resistive diffusivities, to large-eddy simulations (LESs) with numerical diffusive terms. They demonstrated a similar behaviour of the SSD between the DNSs and LESs at different Reynolds and magnetic Reynolds numbers. Warnecke et al. (2023) showed that in DNSs approaching the solar parameter regime the SSD is able to operate (for a detailed review of the SSD in the context of quiet Sun magnetism, see Rempel et al. 2023).

Previous simulations of the quiet Sun are often performed with a weak imposed field. Martínez-Sykora et al. (2019) performed a simulation of the quiet Sun at high resolution, using the Bifrost code (Gudiksen et al. 2011), including a chromosphere and corona. This simulation included a weak (2.5 G) guide field, which is amplified by near-surface convection. This simulation includes non-local thermodynamic equilibrium (NLTE) radiative transfer effects; however, it does not include a non-equilibrium (NE) treatment of hydrogen ionisation. They find a chromosphere dominated by shocks, which are able to amplify the magnetic field. In a subsequent paper, Martínez-Sykora et al. (2023) found that including both NE ionisation and a generalised Ohm’s law leads to a hotter chromosphere.

A number of studies have investigated the capability of SSD generated fields to supply an energy flux into the solar atmosphere, sufficient to counteract chromospheric radiative losses. Amari et al. (2015) demonstrated that the SSD model could provide a high enough Poynting-flux to heat the solar chromosphere. They used a Boussinesq SSD simulation of the near-surface convection zone, coupled to a simplified atmosphere model. They found the observationally inferred magnetic field strengths could provide a Poynting-flux of 4 × 10⁶ erg cm⁻² s⁻¹ into the upper atmosphere, and the re-connection of these fields could reproduce chromospheric features. Similar results were found in a LTE MURaM¹ simulation of the quiet Sun in Tilipman et al. (2023). They found that the Poynting flux in the lower atmosphere is significantly higher than the canonical values required to heat the chromosphere 2.28 × 10⁷ erg cm⁻² s⁻².

At the photosphere, the gas pressure dominates and the effect of SSD fields on the convective dynamics is relatively minor. Simulations with SSD fields have slightly smaller granules, reduced convective upflow velocities and form bright points (Bhatia et al. 2024). In the chromosphere, as the plasma-β drops below unity, the presence of the magnetic field comes to dominate the dynamics. Simulations of the corona and chromosphere including a magnetic field self-consistently generated through a SSD mechanism have been performed with the MURaM code (Rempel 2017; Chen et al. 2021, 2022). These simulations are large enough to contain several super-granules and self-consistently generate network-scale fields. The resulting simulation contains a self-consistent corona. These simulations are performed at relatively low resolutions (50 − 120 km) and use an equation of state (EoS) with a LTE treatment of hydrogen ionisation. Despite the success of photospheric SSD simulations in reproducing photospheric observations, the ability of these fields to heat the chromosphere has not yet been self-consistentlymodelled.

In this paper we investigate the structure and dynamics of the magnetic field in the internetwork chromosphere utilising a NE rMHD simulation. We study the dynamics and structure of the internetwork chromosphere generated purely by an SSD. There is no net flux imposed on the simulation, nor any prescribed flux emergence. This work includes for the first time a SSD simulation including an NLTE treatment of the radiative losses and hydrogen ionisation in the chromosphere. The simulation has a high resolution and low diffusion in order to resolve the small-scale features of the quiet Sun. We study the impact of SSD fields in the chromosphere, and the structure of the resulting chromosphere and corona of the model.

In Section 2 we explain the numerical methods and experimental set-up. In Sect. 3.1 we provide an overview of the simulation, including the resulting photospheric field. The chromospheric magnetic field generated by the SSD mechanism is shown in Sect. 3.2, the energy balance in Sect. 3.3, and the resulting Poynting flux in Sect. 3.4. Finally, in Sect. 4 we discuss the results and in Sect. 5 present our conclusions.

2. Numerical approach

We performed a 3D rMHD simulation using the MURaM code (Vögler et al. 2005). The code was extended to include the physics required to model the corona (Rempel 2017) and the NLTE chromosphere (Przybylski et al. 2022). We used the hyperbolic heat conduction and an Alfvén speed limiter that incorporates a semi-relativistic treatment with a reduced speed of light (Boris 1970), implemented through a projection operator in the momentum equation (Gombosi et al. 2002). In order to accurately treat the solar chromosphere we included a 3D scattering multi-group scheme (Skartlien 2000; Hayek et al. 2010), chromospheric line losses (Carlsson & Leenaarts 2012), and a time-dependent treatment of hydrogen ionisation (Sollum 1999; Leenaarts et al. 2007). The simulation code solves the MHD equations in a cartesian geometry for density (ρ), velocity (v), magnetic field (B), and hydrodynamic energy $(E_{\mathrm{hydro}}=E_{\mathrm{int}}+\rho \frac{v_x^2+v_y^2+v_z^2}{2})$, where E_int is the internal energy. We did not solve for the total energy to prevent numerical errors when plasma-β becomes very low in the corona. The simulations do not include a generalised Ohm’s law (Braginskii 1965). We provide a summary of the numerical set-up in Appendix A.

For the EoS we followed the method described in Przybylski et al. (2022). We used a pre-tabulated LTE EoS in the interior, merging it to the NE EoS at a pressure of p_cutoff = 3 × 10⁵ dyn cm⁻². The NE EoS includes a time dependent treatment of hydrogen, and H₂ molecules, while all other elements are treated in LTE. We solved the advection of the hydrogen populations with the bulk fluid motions, the system of rates governing the ionisation or recombination, and excitation or de-excitation of atomic hydrogen, and the gas-phase transitions forming or destroying molecular hydrogen. A Newton-Raphson method is used to solve the system of rate equations for the number density n for a 5 level hydrogen atom, protons and molecular hydrogen, as well as temperature (T) and electron number (n_e). The pressure can then be calculated p = nk_BT, where n is the total number density.

We used the FreeEoS package (Irwin 2012) to generate the LTE EoS, with the abundances of Asplund et al. (2009). For the NE EoS, a slight modification is made to the convergence scheme described in Przybylski et al. (2022). Due to the low magnetic field strengths in the corona the time-step can reach values of Δt ≥ 10⁻² s. These large time-steps can cause slow convergence of the NE EoS. To address this we applied the solution of the hydrogen ionisation in two steps. First we solved the system of rates with a time-step of 10⁻⁴ s, which settles the fast rates. A second step was then made using the time-step of the MHD solver (Δt_mhd). The MHD timestep is determined by the maximum fast-wave velocity in the simulation $\mathrm{\Delta }{t}_{\mathrm{mhd}}=\min \left(\frac{\mathrm{\Delta }x}{v_A+c_s+|\mathbf{v}|}\right)$, in terms of the Alfvén velocity (v_A), the sound speed (c_s), and the magnitude of the fluid velocity (|v|). In order to prevent unrealistic temperatures, a floor on the temperature was set to 2200 K, and on internal energy to 1.6 × 10¹² erg g⁻¹. The heating term does not contribute significantly towards the energy balance of the chromosphere, and acts only to improve the convergence of the NE EoS at low temperatures.

For the scattering multi-group radiation scheme we used opacities calculated using the MPS-ATLAS package (Witzke et al. 2021). The multi-group scheme includes 4 bands with the binning determined by the formation height relative to a reference wavelength. We used the continuum at 500 nm as a reference, with bin boundaries located at log₁₀τ₅₀₀ = −2.5, −1.5, −0.5.

The diffusion scheme follows Rempel (2017) with slight modifications. This scheme uses a slope-limiter, acting on the primitive variables, combined with a h-value (see Rempel 2014) which limits diffusion in sufficiently smooth areas. A larger h-value allows a steeper gradient to be reached before limiting occurs, and therefore will give a lower numerical diffusion. In the photosphere we followed Rempel (2014) and use a value of 2 for the mass (h_m), velocity (h_v), energy (h_E), and magnetic field diffusion (h_B). For the chromospheric regions, below a density of 5.0 × 10⁻⁹ g cm⁻³ we increased the diffusion (reduce the h-parameter) to h_m = 1.5, h_v = 1.0, h_E = 1.0 and h_B = 1.0 to minimise the overshoot at shock fronts. The reduced mass diffusion, relative to the energy diffusion, prevents the formation of an instability, which occurs in cold chromospheric shock expansions with strong flows, as described in Przybylski et al. (2022). This method supersedes the use of increased diffusion in regions with low-γ, where γ is the adiabatic index, which was used in the previous work. Finally, when the density is lower than 10⁻¹² g cm⁻³ we use a value of h_m = 1.5, h_v = 1.0, h_E = 1.0, and h_B = 4.0 to achieve a high-Prandtl number regime, as described in Rempel (2017). We applied diffusion separately to the hydrogen populations and the instantaneous internal energy terms, which includes the non-hydrogen components and the microscopic kinetic energy. Additional hyper-diffusion in the z-direction was applied to the vertical momentum, density, internal energy and vertical magnetic field.

The lower boundary condition is the Open Boundary, Symmetric Field (OSb) of Rempel (2014). This boundary condition is open to flows. The entropy in the boundary cells is interpolated from a stellar structure model of the Sun, and the pressure is extrapolated into the ghost cells. The magnetic field is passively advected (set symmetrically in the ghost cells). This allows horizontal field to be advected into the domain, mimicking the effects of recirculation deeper in the convection zone. The upper boundary is open to outflows, and imposes a potential magnetic field (Rempel et al. 2009). The horizontal boundary conditions are periodic.

The simulation was initialised with a horizontally isotropic energy and density profile, based on a G2V stellar model, extending to 1 Mm above the photosphere. The governing parameters of the simulation are summarised in Appendix A. A small random velocity field was added and the simulation was run until convection stabilises. A zero net flux seed field was added, consisting of a vertical (B_z) component, which was set to a random value for each column and is constant in the vertical direction. This gives a zero net-flux magnetic field with an RMS field strength of ≈10⁻³ G. The simulation was run for two days of solar time to allow the field to saturate. The behaviour of the SSD is dependent on properties of the turbulent convection, such as the Mach number, Reynolds number, magnetic Reynolds number, and magnetic Prandtl number. Due to the numerical diffusion scheme used in large eddy simulations, like those performed in this work, it is not possible to unambiguously determine the Reynolds number (based on the electrical conductivity of plasma). Estimates can be made based on the integral scale of turbulence, and the Taylor microscale (Pietarila Graham et al. 2010; Moll et al. 2011). These methods typically underestimate the magnetic Reynolds number, giving a value lower than the threshold for the SSD instability. Comparison between large eddy simulations and simulations including an explicit resistivity show that the SSD displays similar behaviour (Vögler & Schüssler 2007; Brandenburg & Rempel 2019).

The simulation was then extended to 11 Mm above the surface using a constant atmosphere in hydrostatic equilibrium with an arctan temperature distribution and potential magnetic field extrapolation. The simulation was allowed to evolve for another day with an Alfvén speed limit of 100 km s⁻¹. Finally the Alfvén limiter was relaxed to 3000 km s⁻¹ and the simulation was run including the NE treatment of hydrogen for 60 minutes. The simulation was then run for an additional 30 minutes where we saved the state of the atmosphere every 10 s. The final simulation spans a region of 12 × 12 × 18 Mm, ranging from −6.6 Mm below to 11.4 Mm above the solar surface. The simulation includes 512 × 512 × 900 grid points giving a horizontal resolution of 23.4375 km and a vertical resolution of 20 km.

3. Results

3.1. Simulation overview

We first investigated the photosphere and low chromosphere of the simulation in Figure 1. A typical quiet Sun granulation pattern is seen in the continuum radiation (Panel a) and vertical velocity (Panel b). At the photosphere, the magnetic field consists of numerous salt-and-pepper concentrations (Panel d). The vertical velocity and vertical magnetic field are taken on the surface where τ₅₀₀ = 1. A number of flux concentrations over 1000 G are seen, and the average maximum unsigned vertical field strength over the time series is |B_z|_max = 2147 ± 199 G. The mean of the horizontally averaged unsigned vertical field is ⟨|B_z|⟩_havg = 71.6 ± 2.98 G and the RMS magnetic field strength is B_rms = 141.3 ± 7.42 G. The horizontal component of the field at the photosphere, $B_h=⟨\sqrt{B_x^2+B_y^2}{⟩}_{\mathrm{havg}}=103.6\pm 9.07\mathrm{G}$. This is consistent with the inferred field strengths of Hanle measurements of around 130 − 170 G (Trujillo Bueno et al. 2004; Zeuner et al. 2024), and previous simulations with the MURaM code (Rempel 2014; Bhatia et al. 2022). In order to see the upper photosphere we plot in panel c the intensity of the 4th radiation band of the 3D multi-group scheme. This band consists of frequencies (ν) where the lines become optically thick (τ_ν = 1) at heights in the atmosphere with τ₅₀₀ < 10^−2.5, including upper-photospheric and chromospheric spectral lines. The image reveals a shock interference pattern with persistent bright points above the strong magnetic field concentrations.

Fig. 1. Snapshot of the simulation. Panel a) Intensity in the continuum (first) band of the multi-group scheme; panel b) Vertical velocity at τ500 = 1; panel c) Intensity in the low-chromospheric (fourth) band of the multi-group scheme; and panel d) Vertical magnetic field strength at τ500 = 1. The maximum field strength has been saturated to better illustrate the salt-and-pepper field concentrations. The dashed line represents the slice shown in Figs. 2 and 8. Animation available online; the image represents t = 0 of the animation.

The physical properties of the chromosphere are illustrated in Fig. 2. The chromosphere is seen to be highly structured with large variations in temperature (panel b) and density (panel d) on small horizontal scales. This variability makes it difficult to meaningfully separate the chromosphere into regions based on a geometric height above the photosphere. We instead followed an approach similar to that of Withbroe & Noyes (1977) to distinguish between the different layers of the solar atmosphere. The near-photosphere extends upwards from the surface where ⟨τ₅₀₀⟩  =  1, to the temperature minimum. We define the base of the chromosphere as the height at which τ₅₀₀ = 10⁻⁵. To simplify analysis, we took a fixed horizontal slice at the height where the horizontally averaged ⟨τ₅₀₀⟩_havg = 10⁻⁵ as the base of the chromosphere. This slice lies approximately 800 km above the average τ₅₀₀ = 1 layer.

Fig. 2. Vertical slice through the simulation taken in the xz plane at y = 6 Mm. Panel a) Temperature; b) Vertical velocity; c) Vertical magnetic field; and d) Plasma density. We used an arcsinh norm for the magnetic field, and limit all panels to highlight the chromosphere. The horizontal dashed lines represent the photosphere τ500 = 0 and base of the chromosphere z = 0.8 Mm. Animation available online; the image represents t = 0 of the animation.

Above the temperature minimum, in the low chromosphere, hydrogen is still largely neutral, plasma-β is predominantly larger than unity, and acoustic shocks dominate the dynamics. In the mid chromosphere the ionisation fraction of hydrogen begins to rise, and the plasma temperature increases to ≈8 kK. In the upper chromosphere, hydrogen is almost fully ionised and the material is heated to near-transition region temperatures, 10–50 kK, as helium is ionised. A large amount of the upper chromospheric material exists at approximately 10 and 30 kK, the preferred temperatures of the first and second ionisation stage of helium in LTE (Golding et al. 2016). Finally, above ≈50 kK a transition region and corona forms.

The tenuous transition region and corona that forms in the simulation is strikingly different from previous simulations which include additional imposed fields (Carlsson et al. 2016; Martínez-Sykora et al. 2019; Breu et al. 2022). When large-scale magnetic fields are present, the plasma-β is always below zero in the upper atmosphere and the dynamics become field aligned, giving the chromosphere its observed spicular nature. Additionally, in smaller simulations a hot plate is often used to force the formation of a transition region and corona (e.g. Leenaarts et al. 2007; Hansteen et al. 2006; Rempel & Przybylski 2021). The atmosphere in the present simulation is more similar in appearance to the flux-emergence experiment of Hansteen et al. (2023), or in large quiet Sun simulations including network fields (Chen et al. 2021; Rempel 2017). In these models, chromospheric material can be seen to reach high into the atmosphere. The structure of the internetwork in these models, in regions with only a small net flux, has not been studied in detail. There is observational evidence of a more corrugated transition region. Trujillo Bueno et al. (2018) state that a more corrugated transition region than present in the Bifrost model is required to match observations from the Chromospheric Lyman-Alpha SpectroPolarimeter (CLASP) rocket experiment (Kano et al. 2017). A study of the Mg II h&k lines in a MURaM enhanced network model (Ondratschek et al. 2024) came to a similar conclusion.

The animation of Fig. 2 reveals a highly dynamic chromosphere. Large regions of the upper atmosphere exist at transition region temperatures and above. The bulk of the low-to-mid chromosphere has velocities in the range of 5 − 10 km s⁻¹. In the upper chromosphere features reaching above 50 km s⁻¹ are seen. Shock fronts propel chromospheric plasma of ≈10 kK into the upper atmosphere, and strong shock expansions frequently reach the minimum allowed temperature of 2200 K. The animation of Fig. 2 reveals regions of localised heating at the interaction of shocks. Additionally small-scale reconnection events and bi-directional jets can be seen, for example at x = 5 Mm, z = 5.5 − 6 Mm, t = 14 − 15 s. The shocks appear to carry chromospheric material, and magnetic field, into the upper atmosphere. The vertical component of the magnetic field (panel c) reveals both persistent loops (for example at x = 10 Mm) and highly turbulent, quickly evolving small-scale field structures. Strong up-flows can be seen to constantly supply magnetic field high into the chromosphere. For example, at about x = 3 Mm, t = 4 minutes a flux emergence event occurs, causing a cool loop to rise up to about z = 5 − 6 Mm.

The averaged thermodynamic properties of the chromosphere are plotted in Fig. 3. The horizontally averaged hydrogen ionisation fraction (top panel) drops sharply at the photosphere, reaching a minimum value around z = 400 km. The temperature minimum in the simulation is at about z = 800 km. This corresponds closely to the height we have defined as the base of the chromosphere, based on the horizontally averaged ⟨τ₅₀₀⟩_havg = 10⁻⁵. The top of the chromosphere is difficult to define due to the large variations of temperature and the presence of chromospheric material high in the atmosphere. We define a chromospheric filling factor as the area fraction at each height with a temperature below 20 kK. The chromospheric filling factor begins to fall at 2 Mm, reaching 50% at 5.5 Mm and falling to 30% at the top of the simulation domain. The average temperature reaches typical low transition region values, around 20 kK at 3.5 Mm. The average temperature plateaus at about 60 kK towards the top of the simulation domain. In the upper atmosphere, approximately 10 pixels in a million reach temperatures greater than 1 MK. The spread of temperature increases in the upper chromosphere, as the NE ionisation of hydrogen prevents recombination, leading to larger temperature variations.

Fig. 3. Horizontally and temporally averaged thermodynamic properties of the simulated chromosphere. Top: Hydrogen ionisation fraction and chromospheric filling factor, which is the area fraction at each height with a temperature below 20 kK. Middle: Temperature (red) and density (black). The shaded areas cover two standard deviations around the mean. The density recovered from hydrostatic equilibrium is included in the middle panel. This was calculated including only the gas pressure (dashed); the gas and turbulent pressures (dash-dotted); and the gas pressure, turbulent pressure, and magnetic term (indistinguishable from the averaged density). Bottom: Error in the density calculated recovered with the hydrostatic equilibrium approximation. The vertical dashed lines represent the photosphere (z = 0 Mm) and the base of the chromosphere (z = 0.8 Mm).

We investigate the extent to which the averaged chromosphere can be described in pressure balance. In hydrostatic equilibrium, the density $\rho =-\frac{{\partial }_{z}p}{g}$ in terms of the pressure p, and gravity g. For a discussion of the hydrostatic balance at the photosphere see Bhatia et al. (2022). In the middle panel of Fig. 3 we compare the averaged simulation density to that calculated including only the gas pressure, to that calculated including the turbulent pressure p_turb = ρv_z² and magnetic term $f_{\mathrm{lor}}=\frac{1}{8\pi }(B_h^2-B_z^2)$. This magnetic term is not the magnetic pressure, but results from the horizontal average of the vertical component of the Lorentz force. In the chromosphere, the density can only be recovered when including both the gas and magnetic pressure terms. In the bottom panel we show the error in the density recovered. When only the gas pressure is included this error is 50–70% in the chromosphere. Including the turbulent pressure but not the Lorentz force term, reduces the error to between 20–50%. This illustrates that the gas, kinetic and magnetic terms are all important in the structuring of the quiet Sun chromosphere.

3.2. The chromospheric magnetic field

The simulation shows a number of small kilo-Gauss magnetic field concentrations at the photosphere. Above the photosphere the magnetic field strength falls off quickly. By the base of the chromosphere, the variation of the horizontally averaged, horizontal field strength is ⟨B_h⟩_havg = 22.14 ± 1.34 G. This is a factor of two higher than the unsigned vertical field, ⟨|B_z|⟩_havg = 11.24 ± 0.81 G. The dependence of magnetic field on geometric height in the atmosphere of the simulation is shown in Fig. 4. We plot the horizontally and temporally averaged values of the unsigned vertical field strength ⟨|B_z|⟩_avg and the horizontal field strength ⟨B_h⟩_avg. The averaged horizontal field is higher than the averaged unsigned vertical field throughout the atmosphere. This difference is greatest in the upper photosphere and low chromosphere, which may be a signature of the presence of low-lying loops, as discussed by Schüssler & Vögler (2008). The horizontal field falls to within a factor of two of the vertical field by 2 Mm. The fact that the two averaged field components approach the same magnitude towards the top of the box could be an effect of the upper boundary conditions. There is a distinct peak in the ⟨B_h⟩_havg about 450 km above the photosphere (Rempel 2014). This peak occurs at the minimum of the horizontally averaged RMS velocity above the photosphere. Rempel et al. (2023) interpret this peak as the height to which turbulent convective motions are able to carry the horizontal magnetic field. Comparison to a potential field extrapolated from B_z at the photosphere (dash-dotted lines) reveals that the field strength in the dynamic simulation is approximately two times higher in the low chromosphere. In the upper atmosphere of the simulation the decrease of the magnetic field strength flattens, so that the field is approximately an order of magnitude higher than the potential field by 6 Mm in height.

Fig. 4. Horizontally and temporally averaged values of the unsigned vertical field |Bz| and the horizontal field Bh = (Bx2+By2)0.5 through the atmosphere. The shaded regions cover two standard deviations. The dash-dotted lines illustrate the magnetic energy assuming a potential field extrapolation from z = 0.

The influence of the magnetic field on the plasma dynamics can be quantified by the plasma-β, $\beta =\frac{8\pi p}{B_x^2+B_y^2+B_z^2}$. A high-β plasma is dominated by the gas pressure (log₁₀β > 0) and a low-β plasma has a higher magnetic pressure. Horizontally averaged log₁₀β is shown in Fig. 5 and the shaded region represents one standard deviation. At the photosphere, the plasma-β is high, values of log₁₀β around or below zero are found only in the inter-granular flux concentrations. At approximately 1.1 Mm the average log₁₀β becomes negative, and by 2 Mm almost all the plasma is low-β. In the mid-to-upper chromosphere the plasma-β flattens, as the majority of closed loops turn over and only large-scale loops remain. The limited simulation domain prevents large-scale flux imbalance from occurring, resulting in low field strengths and regions of high-β towards the upper boundary of the simulation. The minimum in the average log₁₀β occurs at about 3 − 4 Mm above the photosphere, and increases slightly to a value of approximately log₁₀β = −1 for the bulk of the atmosphere.

Fig. 5. Horizontally and temporally averaged values of log10β. The shaded region covers two standard deviations. The vertical dashed lines represent the photosphere (z = 0) and the base of the chromosphere (z = 800 km).

3.3. The chromospheric energy balance

The horizontally and temporally averaged energy densities were calculated for the kinetic E_kin = ρv²/2, magnetic E_mag = B²/(8π), and internal E_int components, and plotted in Fig. 6. At the solar surface ⟨E_int⟩_avg is an order of magnitude higher than ⟨E_kin⟩_avg, and two orders of magnitude higher than the ⟨E_mag⟩_avg. As the density drops the magnetic energy becomes larger than the kinetic energy in the low-chromosphere (≈1 Mm), and larger than the internal energy in the mid chromosphere (≈2 Mm). A greater scatter is seen in the internal energy in the chromosphere, as frequent low-temperature post-shock expansions are present in the domain, as well as plasma at chromospheric and coronal temperatures. The scatter at low energies would be larger without the artificial limits on temperature and internal energy.

Fig. 6. Horizontally and temporally averaged values of the kinetic energy (blue), internal energy (yellow), and magnetic energy (green) in the solar atmosphere. The shaded regions cover two standard deviations. We included the magnetic energy of a potential field above z = 0 (dot-dashed line). The vertical dashed lines represent the height of the photosphere (z = 0) and the base of the chromosphere (z = 800 km).

In the low chromosphere (0.5 − 1 Mm) the kinetic and magnetic energies are roughly in equipartition. The kinetic energy drops sharply in the mid chromosphere, where it is typically around an order of magnitude lower than the magnetic and internal energies. The internal and kinetic energies decrease roughly three orders of magnitude between the photosphere and base of the chromosphere (vertical dashed lines in Fig. 6). This sharp decrease is largely caused by the drop in density above the photosphere, which also falls approximately three orders of magnitude (See Fig. 3). Comparison to a potential field, extrapolated from the photosphere (dash-dotted lines) shows magnetic energies from the simulation are around an order of magnitude higher in the low to mid chromosphere. The energy of the potential field remains significantly lower than the kinetic energy aside from a region from 2 Mm to 4.5 Mm. The magnetic energy in the potential field decreases significantly by the upper chromosphere, and is approximately two orders of magnitude lower at 8 Mm above the photosphere. In contrast, the magnetic energy in the MHD simulation decreases more slowly, in rough equipartition with the thermal energy.

In panel a) of Fig. 7 we show histograms of the temperature at three different heights. The histograms show a number of peaks, which shift to higher temperatures at higher heights. These peaks occur at the temperature where dissociation of molecular hydrogen, or the ionisation of hydrogen and helium occur. In LTE the peaks are much sharper than in Fig. 7, as the transitions are instantaneous, while in NE the slow ionisation and recombination rates in the chromosphere lead to a broader spread of temperatures (Carlsson & Stein 2002; Golding et al. 2016). The peak at 2200 K at 1 Mm and at 2 Mm is a consequence of the artificial heating term. At a height of 2 Mm two peaks are seen, at the preferred ionisation temperature of hydrogen, and the first ionisation temperature of helium. At 4 Mm, there is a third peak at the preferred temperature of He II ionisation (Golding et al. 2016). Large regions of the upper atmosphere, at 8 Mm, have temperatures of above 100 kK. However, a persistent million degree corona is unable to form, with only approximately 1% of the plasma above the photosphere reaching a million kelvin.

Fig. 7. Histograms of the energy densities at three different heights in the atmosphere and temperature at four different heights. Panel a) shows a normalised histogram of temperature in the chromosphere at 1 Mm, 2 Mm, 4 Mm, and 8 Mm. The lower panels show normalised histograms of internal, magnetic, and kinetic energy in the chromosphere, at heights of b) 4 Mm, c) 2 Mm, and d) 1 Mm.

To understand the distribution of energy with height, we plot in Fig. 7 histograms of the different energy components at heights of 1 Mm, 2 Mm, and 4 Mm. Above 4 Mm the different energy components decrease slowly with height, with the relative difference remaining similar. The internal energy histogram keeps a similar distribution at different heights in the atmosphere. The magnetic energy histogram becomes narrower with height, while the kinetic energy histogram broadens with height. The narrowing of the magnetic energy histogram is expected as the magnetic field becomes space filling and plasma-β is lower than one. The broadening of the kinetic energy histogram is also expected due to the larger contribution of shocks with height, with a high density contrast between shock fronts and expansions. The average magnetic energy in the upper chromosphere is substantially higher than the kinetic energy as seen in Fig. 6. However, there is still substantial overlap with the wings of the kinetic energy distribution and the magnetic energy distribution.

The magnetic fields in the chromosphere are substantially higher than expected from a potential extrapolation of the photospheric fields. The distribution of magnetic energy in the upper chromosphere can be better understood by looking at the spatial and temporal variation of the energies through a slice in the simulation domain. In Fig. 8 we show the internal, kinetic, andmagnetic energies as well as the plasma-β for the xz-slice through the atmosphere presented in Fig. 2. The animation of the figure clearly illustrates that shock fronts with high thermal and kinetic energy, where plasma-β ≈ 1, are able to push the magnetic field high into the atmosphere. This results in a relatively smooth distribution of the magnetic field with height and explains much of the excess field strength relative to a potential field in the upper chromosphere. Above ≈1 Mm the simulation is predominately low-β, although regions with β ≥ 1 are present in shock fronts, or in low-field regions of the upper atmosphere.

Fig. 8. Spatial variation of energy density in the xz plane at y = 6 Mm, including a) internal energy, b) kinetic energy, c) magnetic energy, and d) plasma-β. The horizontal dashed lines represent the height of the photosphere (z = 0) and base of the chromosphere z = 0.8 Mm. Animation available online; the image represents t = 0 of the animation.

3.4. Poynting flux

The magnetic energy at the photosphere represents most of the free energy, making it a strong candidate for heating the chromosphere and corona. We now investigate the magnetic energy (Poynting) flux into the chromosphere. Figure 9 shows the horizontally and temporally averaged vertical component of the Poynting flux P_z. The vertical Poynting flux is further split into two parts, the advection of horizontal field $P_{z,v_z}=\frac{v_z}{4\pi }(B_y^2+B_x^2)$, and a term proportional to the vertical field $P_{z,B_z}=-\frac{B_z}{4\pi }(B_x{v}_x+v_y{B}_y)$. The latter term represents horizontal motions around a vertical field, including motions with a shear term, braiding, and vortical or wave motions. Vigorous near-surface convection leads to amplification of magnetic energy near the surface, and recirculation of this field gives a net negative Poynting flux below the surface. Above the photosphere, the horizontally and temporally averaged vertical Poynting flux P_z = P_z, vz + P_z, Bz is positive. The majority of this Poynting flux is provided through the P_z, Bz term. At the base of the chromosphere, which is the height at which ⟨τ₅₀₀⟩ = 10⁻⁵, the Poynting flux is ⟨P_z⟩_havg = (5.06 ± 1.08)×10⁶ erg cm⁻² s⁻¹. This is split into the advective term ⟨P_z, vz⟩_havg = (1.86 ± 1.07)×10⁶ erg cm⁻² s⁻¹, and the ⟨P_z, Bz⟩_havg = (3.21 ± 0.12)×10⁶ erg cm⁻² s⁻¹. The P_z, vz term is negative at the photosphere, increasing to a maximum in the low chromosphere before becoming negative through the mid-to-upper chromosphere up to about 4 Mm, when it turns positive again.

Fig. 9. Horizontally and temporally averaged values of the vertical Poynting flux in the solar atmosphere. The top panel shows its stratification in the photosphere and low chromosphere and the bottom panel displays its behaviour in the chromosphere and above on an expanded vertical scale. The blue line shows the shear (wave) −Bz(vxBx + vyBy), the red line the flux advection vzBh2, and the black line the total vertical Poynting flux. The solid lines show the mean value and the shaded area shows two standard deviations.

The time-dependence of the horizontally averaged Poynting flux is shown in Fig. 10. The P_z, Bz component of the Poynting flux does not vary strongly with time. In contrast, the P_z, vz component varies. Although largely negative in the chromosphere it periodically supplies magnetic energy into the upper chromosphere. The Poynting flux term representing the vertical advection of horizontal-field P_z, vz varies approximately with the 5-minute turnover time of solar convection. In a simulation domain with limited size and depth, only a limited range of modes can be excited, and these modes are excited with excess power (see Appendix B). The oscillations seen in Fig. 10 are significantly shorter than the fundamental box mode, which has a period of 12 minutes. Instead, a clear signature is seen with a 5 minute periodicity, which are the acoustic p-modes that propagate throughout the solar interior. In the upper atmosphere of the simulation the shear(wave) term decreases quicker with height than the advective term. By 9 Mm the advective term is positive and carries 60% of the total Poynting flux ⟨P_z, vz⟩_{havg, 9 Mm} = (4.71 ± 1.19)×10⁴ erg cm⁻² s⁻¹. The reduction in the shear Poynting flux term could be affected by the height of the upper boundary. A gradient can be seen in the ⟨P_z, Bz⟩ component of the Poynting flux with height. This is caused by a gradual increase of B_tot over the time series studied. The average field strength in the chromospheres evolves slowly on a timescale of ≈1 − 2 hours, which is not captured completely in the limited time series of this study. A sharper increase of B_tot in the low chromosphere occurs due to a large flux emergence event which occurs around 18 minutes into the simulation.

Fig. 10. Time-distance image of the of the horizontal averages of a) the shear (wave) −Bz(vxBx + vyBy) and b) the flux advection vzBh2 components of the vertical Poynting flux. The horizontal dashed lines represent the photosphere z = 0 and base of the chromosphere ⟨τ500⟩havg = 10−5.

Finally, we aim to understand the spatial distribution of the Poynting flux in the internetwork chromosphere. Horizontal slices through the chromosphere are shown in Fig. 11 for three heights, at 1.2 Mm, 2.2 Mm, and 4.2 Mm above z = 0 Mm. An animation of the figure is available. The top row shows the temperature at these different heights. In the lower chromosphere (panel a), the temperature shows a shock interference pattern, with small regions of sustained high temperature present above strong field concentrations. By 2.2 Mm the plasma-β is low and the dynamics begin to follow the magnetic field, becoming highly structured perpendicular to the magnetic field. Most of the plasma is around 20 kK, with large pockets of low-chromospheric temperatures, and small regions reaching transition region temperatures. The plasma is highly structured in the horizontal direction with sharp variations from swirling motions and shock buffeting. In the highest slice most of the plasma remains at upper chromospheric temperatures of around 20 kK, with some regions being heated to reach coronal temperatures, nearly 1 MK. A large swirling event is seen in panel c), at around x = 4 Mm, y = 6 Mm, with significant low chromospheric material reaching high in the atmosphere. It is also seen at 2.2 Mm, although less clearly. The swirl is easier to detect at all heights in other physical parameters, such as B_z or Poynting flux.

Fig. 11. Horizontal cuts in the chromosphere at heights of (from left to right) z = 1.2, 2.2, and 4.2 Mm above z = 0 Mm. The panels (top to bottom) include the temperature, vertical magnetic field Bz, and two components of the vertical Poynting flux, Pz, Bz = −Bz(vxBx + vyBy) and Pz, vz = Bh. The temperature panels are limited to compare the structures at different layers in the chromosphere. In panel c) approximately 8 pixels in a million are above 1 MK. The field strength and Poynting flux panels are limited to highlight the chromospheric structures. Animation available online; the image represents t = 0 of the animation.

The second row shows the vertical magnetic field strength. In the low chromosphere a number of small, sub-granular loops are seen, with field strengths above 50 G. These loops can be identified in the horizontal slices by a smooth transition of polarity along the structure. At 2 Mm the field is dominated by loops on the granular scale, and the field has dropped to a maximum value of around 10 G. At 4 Mm the field shows a large-scale flux imbalance of around 4 G over horizontal scales of a few granules (4 − 8 Mm). Even in the upper chromosphere significant small-scale structure is seen in the magnetic field, down to ≈100 km with opposite polarity loops frequently interacting.

The bottom two rows of Fig. 11 show the P_z, vz and P_z, Bz components of the Poynting flux. The animation of Fig. 11 reveals the dynamics that lead to strong Poynting flux in the chromosphere. A swirl at x = 4 Mm, y = 7 Mm expands into the upper chromosphere and carries significant Poynting flux upwards. The swirl is present at the beginning of the simulation and remains visible for approximately 20 min before starting to decay. A number of smaller-scale swirls can be seen in the simulation, but none has as large an impact on the upper chromosphere. At x = 8 Mm, y = 4 Mm the interaction of a number of small loops carries significant Poynting flux into the upper chromosphere. Most of the other regions of enhanced Poynting flux are more transient, lasting only minutes.

4. Discussion

Whereas previous work included an imposed vertical field, which dominates the field geometry of at least the higher layers of the chromosphere, in the present work, we produce a model of the quiet Sun internetwork generated through a SSD mechanism. The SSD, acting in the upper convection zone of the simulation, generates magnetic fields consistent with photospheric observations of the internetwork field, as has been showed in a number of papers (e.g. Danilovic et al. 2010, 2016; Lites et al. 2008; Lites 2011; Trujillo Bueno et al. 2004; del Pino Alemán et al. 2018; Zeuner et al. 2024) We studied the strength of fields in the chromosphere, and investigate their contribution to the chromospheric energy balance. We demonstrated that a zero net-flux SSD simulation is capable of generating magnetic fields that reach high into the chromosphere. The small simulation domain (12 × 12 Mm) restricts the formation of larger-scale magnetic loops reaching high into the atmosphere and likely dominating the magnetic structure there. Even so, the magnetic field in the upper chromosphere is strong enough that the magnetic energy is more than an order of magnitude higher than the kinetic energy, and around a factor of 2 higher than the thermal energy.

We found that the magnetic energy falls off far more slowly with height than expected from a potential field extrapolation. A mechanism for the enhanced field strengths in the chromosphere is shown to be the interaction between shocks and the magnetic field. The shock fronts are seen to push the turbulent internetwork fields higher into the atmosphere (see for example the animation of Fig. 8). Although the shocks in the chromosphere act to locally enhance the magnetic field, this does not constitute a dynamo process. It is possible for local enhancement of the field to occur without a dynamo process. An additional source of energy is contained in the braiding of the field, caused by small-scale vortical motions, which is ubiquitous in simulations of the solar atmosphere.

These changes will cause extrapolations of quiet inter-network regions to significantly underestimate the magnetic field strength in the chromosphere. This will further impact the inferred pressure balance, density, and free magnetic energy present in the chromosphere. Finally the height of the network canopy will be increased due to the higher pressure in the internetwork regions.

At the photosphere, the behaviour of quiet Sun internetwork fields have been well studied in both observations (Buehler et al. 2013; Gošić et al. 2014, 2016; Lagg et al. 2016; Smitha et al. 2017) and through numerical simulations (Danilovic et al. 2010; Rempel 2014; Bhatia et al. 2022, 2024). Due to the difficulties in measuring chromospheric magnetic fields, and the complexity of detailed chromospheric simulations, the extent to which these fields impact the chromosphere has remained poorly understood. Khomenko et al. (2018) performed a simulation including the low chromosphere, with an upper boundary at 1.4 Mm. These simulations study the propagation of magnetic energy including non-ideal effects such as ambipolar diffusion and the Hall effect. They find a 90% reduction of Poynting flux before the upper boundary of the simulation. However, the upper boundary is too low, so that it is expected to significantly affect the field topology in the chromosphere. Previous numerical studies of the internetwork chromosphere, including high resolution and an upper boundary above the transition region, are limited to the model of Martínez-Sykora et al. (2019) (MS19).

Observational inference of weak quiet Sun fields, in the range of 1–100 G in the chromosphere, typically requires Hanle depolarisation measurements (Stenflo 1982; Trujillo Bueno & del Pino Alemán 2022). Faurobert-Scholl (1992) investigated depolarisation of the Ca I 4227 Å line, which forms in the mid chromosphere, roughly 700 − 1200 km above the solar surface in a 1D plane parallel solar atmosphere model. Assuming a height-dependent formation due to the chromospheric canopy they found a 20 − 100 G field strength. The observations of Bianda et al. (1998), of the same line, measure an average of 5 − 15 G value over the formation height of the line. Lagg et al. (2009) investigated Hanle depolarisation of the He I 10830 Å line and found a field strength of 20 − 50 G in a weak field region. These results are consistent with those in our simulation, with an average field strength of ≈15–30 G found over the formation region 700 − 1200 km.

A previous study of the quiet Sun chromospheric field by MS19 found similar properties of the magnetic field in the low atmosphere. They find the magnetic field is strongly horizontal in the low-chromosphere. Their work included a 2.5 G vertical guide field, leading to the field becoming vertical in the upper chromosphere. The internetwork SSD simulation presented in this work does not include any imposed field, and the limited domain size prevents the formation of large-scale network fields. This causes the chromospheric field to remain highly inclined through the bulk of the chromosphere, in contrast to MS19. We see significantly stronger magnetic energy in the chromosphere, around a factor of 2 higher than MS19. This is likely caused by the larger, deeper convection zone in combination with the reduced diffusivity of the slope-limited numerical diffusion scheme, leading to a stronger dynamo-generated magnetic field in the photosphere. In the internetwork SSD simulation, there is a close equipartition of kinetic and magnetic energy in the low chromosphere, and magnetic energy begins to dominate in the mid chromosphere. The kinetic energies in the two simulations are of a similar magnitude at 1 Mm, which is a factor of ≈4 higher than the magnetic energy in MS19. By the upper chromosphere the 2.5 G imposed field begins to dominate and the two simulations can no longer be directly compared.

The measured Poynting fluxes in this simulation are consistent with those presented in Tilipman et al. (2023). They present simulations, performed with the MURaM code, include the upper convection zone and low chromosphere of a SSD-driven quiet Sun simulation. Their treatment of radiation and hydrogen ionisation in the chromosphere is limited to LTE. They observe a turnover in the Poynting flux just below the photosphere and a Poynting flux peak of 2.28 × 10⁷ erg cm⁻² s⁻¹ at 128 km. Comparing MURaM simulations to observational data they find that interpretation of the observations (see also Silva et al. 2022) is complicated by the measurement of the field on an optical depth surface. Additionally, it is difficult to measure both the vertical and horizontal field components that are required to determine the shear (wave) component (P_z, Bz). We significantly extend the study of the photospheric Poynting flux by investigating the propagation of energy into the solar chromosphere, where we treat NE hydrogen ionisation and NLTE radiative transfer.

In order to understand the energy flux propagating from the turbulent convection zone into the atmosphere we took spatial and temporal averages at the base of the chromosphere. The spatially and temporally averaged vertical Poynting flux into the chromosphere is ⟨P_z⟩ = 5.01 × 10⁶ erg cm⁻² s⁻¹, approximately 25% higher than the canonical energy flux required to counteract radiative losses in the quiet Sun chromosphere (4.3 × 10⁶ erg cm⁻² s⁻¹) and corona (3 × 10⁵ erg cm⁻² s⁻¹) (Withbroe & Noyes 1977). This is consistent with large simulations, such as Rempel (2017) and Chen et al. (2021), that supply an average vertical Poynting flux of approximately 5 × 10⁵ erg cm⁻² s⁻¹ to the corona. Therefore, for the action of fields self-consistently generated by an SSD, convective turbulence easily provides a large enough energy flux into the chromosphere. We have demonstrated that this is the case even for zero net magnetic flux.

In order to determine the complete energy balance of the modelled chromosphere we additionally calculated the vertical kinetic energy flux, F_z, Ekin = 1/2ρv_z(v_x²+v_y²+v_z²) and enthalpy flux F_z, H = v_z(E_int+p). The averaged kinetic energy flux into the chromosphere, ⟨F_z, Ekin⟩_avg = −4.8 × 10⁵ erg cm⁻² s⁻¹, is an order of magnitude lower than the Poynting flux and negative. The averaged vertical enthalpy flux at the base of the chromosphere is high ⟨F_z, H⟩_havg = 3.63 × 10⁶ erg cm⁻² s⁻¹. This value is highly sensitive to the height at which the energy flux is measured. It is negative below 450 km, with a peak at 600 km, falling an order of magnitude by 1 Mm. This strong peak is presumably related to the height of the canopy in the low chromosphere (e.g. Solanki & Steiner 1990), the sharp decrease in Poynting flux in the upper photosphere leads to an upwards enthalpy flux in the low chromosphere. To complete the energy balance, we include radiative losses in the chromosphere $\sum _{z=0.8\mathrm{Mm}}^{10\mathrm{Mm}}⟨Q_{\mathit{rad}}{⟩}_{\mathrm{havg}}dz=-8.44\times {10}^6\mathrm{erg}{\mathrm{cm}}^{-2}{\mathrm{s}}^{-1}$. Here dz is the vertical resolution and Q_rad is the total radiative losses/heating rate from the multi-group scheme, chromospheric lines, optically thin losses, and EUV back heating. These losses are nearly double those expected from the values of Withbroe & Noyes (1977). Possible causes of these differences could be the difficulty in defining the ‘base of the chromosphere’ in a dynamic, 3D model, long term variations of the chromosphere, or the lack of a canopy and network in the simulation.

Attributing chromospheric heating to a particular energy flux at the base of the chromosphere is difficult, and would require a more detailed study than presented in this work. The RMS values of both components of the Poynting flux and the kinetic energy flux are of a similar magnitude, but the mean values differ significantly. The shear (P_z, Bz) component of the Poynting flux and the enthalpy flux are the dominant contributors of energy into the chromosphere. An approach such as the Lagrangian ’corks’ described by Leenaarts (2018) could help to understand the mass transfer between the photosphere, chromosphere and corona. However, tracing the transport of energy carried by Alfvénic waves is complicated by the complex topology and frequent restructuring of the chromospheric magnetic fields. The dissipation of these energies, and the region of the chromosphere in which they are thermalised will be dependent on the resolution and numerical diffusion scheme. Additional sources of energy dissipation, for example through ion-neutral or multi-fluid effects (Khomenko 2017; Martínez-Sykora et al. 2012, 2020) will change the dissipation of kinetic and magnetic energy.

We do not include a generalised Ohm’s law to model the effect of ambipolar diffusion. Khomenko et al. (2018) perform a LTE simulation of SSD generated fields, including a generalised Ohm’s law, with an upper boundary in the low chromosphere (1.4 Mm). In their simulation the inclusion of ambipolar diffusion leads to a significant dissipation of Poynting flux in the solar chromosphere. The inclusion of NE ionisation could significantly change these results as hydrogen is the major electron donor in the chromosphere. The study of MS19 was extended to include NLTE effects by Martínez-Sykora et al. (2023), where the inclusion of NE ionisation and a generalised Ohm’s law leads to increased temperatures in the upper chromosphere. However, they do not include results from a simulation including NE ionisation without a generalised Ohm’s law, to allow a direct comparison to this study. In the presented simulation the energy into the chromosphere is a conserved quantity, with the largest diffusivity dissipating the most energy. Only the details of the distribution of energy dissipation between the different terms will change as resolution is increased, or if non-ideal effects are included. The results presented in this paper are therefore expected to be robust.

The extent of the numerical domain is insufficient in depth and width to self-consistently model supergranulation. The horizontal motions at scales of supergranulation (≈20 Mm) advect field to edges of the supergranules. This produces a net flux imbalance on these scales, and network loops will reach higher and likely contribute to heating the plasma in higher layers to coronal temperatures. As seen in Fig. 6, the scale height of the vertical magnetic field in the upper atmosphere is roughly equal to k = 2π/6 Mm. In a larger SSD simulation, such as Chen et al. (2021), the magnetic fields emerging in the internetwork interact with these large-scale fields. The horizontally averaged RMS field strength is ≈3 G at 4 Mm above the photosphere (see Fig. 4). A typical super-granular flux imbalance is approximately ±2.8 G (Smitha et al. 2017). In the upper chromosphere, the dynamics seen in the current simulation will become dominated by the magnetic canopy provided by the magnetic network. The emerging internetwork fields, which currently propel material to the top of the simulation domain, will interact with these overlying fields. The interaction between the highly turbulent internetwork fields and canopy will drive reconnection events, further heating the chromosphere and corona.

5. Conclusion

In this work we demonstrate that even in the quietest regions of the Sun the magnetic field is strong enough to dominate the energy balance above the mid chromosphere. The kinetic and magnetic energies are in equipartition throughout the low chromosphere, while in the higher layers the magnetic energy density is an order of magnitude higher. In the higher layers the magnetic energy is also a factor of 2 higher than the thermal energy density. This equipartition is caused by the advection and compression of field by shocks in the lower atmosphere, and the anisotropy of the field and hot plasma in the upper atmosphere. The resulting field provides a Poynting flux strong enough to heat the chromosphere. Additionally the magnetic canopy in the upper photosphere produces a significant enthalpy flux into the chromosphere.

However, the resulting simulation is unable to form a steady million degree corona. Only small, fluctuating pockets of million-degree gas are formed. This is likely due to the too low magnetic energy at greater heights in our simulation. Within the limited spatial domain of the simulation, a magnetic network is not formed. Therefore, no strong magnetic loops which reach high into the atmosphere are present in the simulation. Longer loops will provide a mechanism to allow a larger Poynting flux to reach higher into the atmosphere.

Data availability

Movies associated with Figs. 1, 2, 8, and 11 are available at https://www.aanda.org

Acknowledgments

D.P. would like to thank T. Bhatia for discussions on the SSD. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 101097844 ). The work of DP was funded by the Federal Ministry for Economic Affairs and Climate Action (BMWK) through the German Space Agency at DLR based on a decision of the German Bundestag (Funding code: 50OU2201). We gratefully acknowledge the computational resources provided by the Raven supercomputer of the Max Planck Computing and Data Facility (MPCDF) in Garching, Germany. D.P. would like to thank A. Irwin (Free-EoS) and V. Witzke (MPS-ATLAS). This material is based upon work supported by the NSF National Center for Atmospheric Research, which is a major facility sponsored by the U.S. National Science Foundation under Cooperative Agreement No. 1852977. This project has received funding from the Swedish Research Council (2021-05613) and the Swedish National Space Agency (2021-00116).

References

Appendix A: Simulation set-up

In this section we provide a summary of the numerical set-up used in the simulation. Table A.1 includes information on the numerical scheme, boundary conditions and physics modules included. In Table A.2 we provide the governing parameters of the run. These include the simulation size and resolution, and the net magnetic flux. The type of star which will be simulated is set by the inflowing entropy and pressure at the lower boundary, the gravity and elemental abundances.

Appendix B: Simulation box mode

Turbulent near-surface convection excites modes with a range of frequencies and wavenumbers. In a simulation of limited width and depth, with a partially absorbing lower boundary, the power is spread over a limited set of modes. In particular, a global box mode can be excited with excess power. To visualise this effect, we plot in Figure B.1 the horizontally averaged vertical velocity in the simulation, see also Fig. 8 of Carlsson et al. (2016). At photospheric heights a periodicity of approximately 5-6 minutes is seen. The period of the fundamental mode can be approximated by calculating the travel time of an acoustic wave between the photosphere to the lower boundary of the simulation $T=2{\int }_{z=-6.8}^0\frac{\mathit{dz}}{C_s}$, where C_s is the horizontally averaged sound speed. This gives a value of 12.8 minutes, which is significantly longer than the dominant oscillation observed in Fig. B.1.

Fig. B.1. Variation in the horizontally averaged vertical velocity in the simulation. The figure shows different heights in the interior and near-photosphere (top panel), and the chromosphere (bottom panel).

All Tables

All Figures

Fig. 1. Snapshot of the simulation. Panel a) Intensity in the continuum (first) band of the multi-group scheme; panel b) Vertical velocity at τ500 = 1; panel c) Intensity in the low-chromospheric (fourth) band of the multi-group scheme; and panel d) Vertical magnetic field strength at τ500 = 1. The maximum field strength has been saturated to better illustrate the salt-and-pepper field concentrations. The dashed line represents the slice shown in Figs. 2 and 8. Animation available online; the image represents t = 0 of the animation.

In the text

Fig. 2. Vertical slice through the simulation taken in the xz plane at y = 6 Mm. Panel a) Temperature; b) Vertical velocity; c) Vertical magnetic field; and d) Plasma density. We used an arcsinh norm for the magnetic field, and limit all panels to highlight the chromosphere. The horizontal dashed lines represent the photosphere τ500 = 0 and base of the chromosphere z = 0.8 Mm. Animation available online; the image represents t = 0 of the animation.

In the text

Fig. 3. Horizontally and temporally averaged thermodynamic properties of the simulated chromosphere. Top: Hydrogen ionisation fraction and chromospheric filling factor, which is the area fraction at each height with a temperature below 20 kK. Middle: Temperature (red) and density (black). The shaded areas cover two standard deviations around the mean. The density recovered from hydrostatic equilibrium is included in the middle panel. This was calculated including only the gas pressure (dashed); the gas and turbulent pressures (dash-dotted); and the gas pressure, turbulent pressure, and magnetic term (indistinguishable from the averaged density). Bottom: Error in the density calculated recovered with the hydrostatic equilibrium approximation. The vertical dashed lines represent the photosphere (z = 0 Mm) and the base of the chromosphere (z = 0.8 Mm).

In the text

	Fig. 4. Horizontally and temporally averaged values of the unsigned vertical field |Bz| and the horizontal field Bh = (Bx2+By2)0.5 through the atmosphere. The shaded regions cover two standard deviations. The dash-dotted lines illustrate the magnetic energy assuming a potential field extrapolation from z = 0.
In the text

	Fig. 5. Horizontally and temporally averaged values of log10β. The shaded region covers two standard deviations. The vertical dashed lines represent the photosphere (z = 0) and the base of the chromosphere (z = 800 km).
In the text

Fig. 6. Horizontally and temporally averaged values of the kinetic energy (blue), internal energy (yellow), and magnetic energy (green) in the solar atmosphere. The shaded regions cover two standard deviations. We included the magnetic energy of a potential field above z = 0 (dot-dashed line). The vertical dashed lines represent the height of the photosphere (z = 0) and the base of the chromosphere (z = 800 km).

In the text

	Fig. 7. Histograms of the energy densities at three different heights in the atmosphere and temperature at four different heights. Panel a) shows a normalised histogram of temperature in the chromosphere at 1 Mm, 2 Mm, 4 Mm, and 8 Mm. The lower panels show normalised histograms of internal, magnetic, and kinetic energy in the chromosphere, at heights of b) 4 Mm, c) 2 Mm, and d) 1 Mm.
In the text

	Fig. 8. Spatial variation of energy density in the xz plane at y = 6 Mm, including a) internal energy, b) kinetic energy, c) magnetic energy, and d) plasma-β. The horizontal dashed lines represent the height of the photosphere (z = 0) and base of the chromosphere z = 0.8 Mm. Animation available online; the image represents t = 0 of the animation.
In the text

Fig. 9. Horizontally and temporally averaged values of the vertical Poynting flux in the solar atmosphere. The top panel shows its stratification in the photosphere and low chromosphere and the bottom panel displays its behaviour in the chromosphere and above on an expanded vertical scale. The blue line shows the shear (wave) −Bz(vxBx + vyBy), the red line the flux advection vzBh2, and the black line the total vertical Poynting flux. The solid lines show the mean value and the shaded area shows two standard deviations.

In the text

	Fig. 10. Time-distance image of the of the horizontal averages of a) the shear (wave) −Bz(vxBx + vyBy) and b) the flux advection vzBh2 components of the vertical Poynting flux. The horizontal dashed lines represent the photosphere z = 0 and base of the chromosphere ⟨τ500⟩havg = 10−5.
In the text

Fig. 11. Horizontal cuts in the chromosphere at heights of (from left to right) z = 1.2, 2.2, and 4.2 Mm above z = 0 Mm. The panels (top to bottom) include the temperature, vertical magnetic field Bz, and two components of the vertical Poynting flux, Pz, Bz = −Bz(vxBx + vyBy) and Pz, vz = Bh. The temperature panels are limited to compare the structures at different layers in the chromosphere. In panel c) approximately 8 pixels in a million are above 1 MK. The field strength and Poynting flux panels are limited to highlight the chromospheric structures. Animation available online; the image represents t = 0 of the animation.

In the text

	Fig. B.1. Variation in the horizontally averaged vertical velocity in the simulation. The figure shows different heights in the interior and near-photosphere (top panel), and the chromosphere (bottom panel).
In the text