© The Authors 2025

1. Introduction

The chromosphere is the region located above the solar photosphere, characterized by strong brightness in the Hα waveband seen in emission during eclipses. In one-dimensional (1D) empirical models of the solar atmosphere, the temperature in the chromosphere gradually increases with altitude, from approximately 4000 K–8000 K at around 2000 km above the bottom of the photosphere (Avrett & Loeser 2008). High-resolution observations and three-dimensional (3D) radiation-magnetohydrodynamics (MHD) simulations show a considerably more complex picture. The chromosphere is structured and dynamic, and the transition region highly corrugated. In the chromosphere, the energy transmitted from the solar interior is redistributed, with the majority radiated into space, and a smaller portion transported to the corona (Jess et al. 2009). Understanding the chromosphere, the distribution of energy within it, and its coupling to the corona is essential for addressing key questions in solar physics, such as coronal heating and the driving of the solar wind (Cranmer & Winebarger 2019; Van Doorsselaere et al. 2020). Additionally, the chromosphere itself has many interesting issues to be solved, such as its heating and the formation of fine structures such as spicules. A recent review of the chromosphere has been provided by Carlsson et al. (2019), offering a comprehensive overview of researchers’ current understanding and outstanding questions in the field.

Simulating the evolution of the lower solar atmosphere is challenging, partly because radiation plays a critical role in energy transfer and interacts strongly with the plasma. Thus the simulation codes must solve the radiative transfer (RT) equation to model the radiative energy exchange processes (Leenaarts 2020). In the photosphere, which is close to local thermodynamic equilibrium (LTE), the level populations follow Saha-Boltzmann statistics. The source function of the main spectral band, the continuum spectrum, can then be approximated by the temperature-dependent Planck function. The chromosphere is particularly challenging due to non-local thermodynamic equilibrium (NLTE) effects. In the chromosphere, strong spectral lines such as the Ca II, Mg II, and H I Lyman-α lines play a key role in radiative energy transfer (Vernazza et al. 1981). Due to long ionization and recombination timescales in the chromosphere, the level populations of the atoms and ions, required for the solution of the RT equation, become decoupled from the local temperature (Carlsson & Stein 2002). Since the evolution of these populations depends in turn on the radiation field, this interdependent relationship makes the modeling of radiative energy exchange complicated and potentially global.

In terms of computational effort, fully solving the RT equation to obtain an accurate, frequency-dependent radiation field for spectral lines is highly expensive – not only due to the large number of frequency points required, but also because of the effects of partially coherent scattering (Leenaarts 2020). Due to the high computational cost, fully solving 3D RT is currently only feasible in spectral synthesis, such as with MULTI3D (Leenaarts & Carlsson 2009) and PORTA (Štěpán & Trujillo Bueno 2013). Time-dependent calculations including a detailed treatment of the radiation transfer are typically restricted to 1D simulations (e.g., Carlsson & Stein 1997; Kašparová et al. 2009).

In chromospheric simulations, the primary emphasis is typically placed on the evolution of the plasma state. In particular, the evolution of hydrogen receives significant attention due to its high abundance and dominant role in plasma dynamics. When modeling the interaction between matter and radiation, using simplified radiation fields is an effective approach to reduce computational complexity. One widely used method is that of Sollum (1999), which approximates the radiation field of hydrogen non-Lyman lines using a Planck function evaluated at a specified local radiation temperature. This method has been implemented in several simulation codes, including CO5BOLD (Leenaarts & Wedemeyer-Böhm 2006), the Oslo Stagger Code (Leenaarts et al. 2007), Bifrost (Golding 2010; Gudiksen et al. 2011), and MURaM (Przybylski et al. 2022), which incorporate NLTE ionization effects in modeling chromospheric plasma. An alternative approach of simplifying the radiation field, utilizing an escape probability method, was suggested by Judge (2017).

The Lyman series, particularly the Lyman-α line, plays a significant role in cooling the chromosphere and thus cannot be ignored in NLTE simulations. One approach to modeling the Lyman series is to assume detailed radiative balance, in which the local upward radiative rate equals the downward radiative rate. Under the assumption of radiative balance, population changes due to radiative excitation and de-excitation are neglected (e.g., Leenaarts et al. 2007; Przybylski et al. 2022). Detailed radiative balance is a good approximation for Lyman-α in the lower chromosphere, but it does not apply in the upper chromosphere (Vernazza et al. 1981; Carlsson & Leenaarts 2012). In the upper chromosphere, radiative de-excitation generally dominates over radiative excitation, with Lyman-α photons carrying a substantial amount of energy into space. Additionally, photons emitted from the transition region are reabsorbed in the upper chromosphere, potentially influencing its thermal and dynamical evolution. Golding et al. (2016) introduced RT of the Lyman-α transition into the Bifrost code and demonstrated that the resulting hydrogen ionization degree near wavefronts differs significantly from that obtained under the assumption of radiative detailed balance. To accurately capture the hydrogen population dynamics in the upper chromosphere, computing the radiation field of the Lyman series by solving the RT equation might be necessary.

In this paper, we present an extension of the MURaM code. Building on the chromospheric extension introduced in Przybylski et al. (2022), we develop a patch to incorporate the RT treatment, initially restricted to modeling the Lyman series. The treatment can easily be expanded to model other species that are important for the chromospheric energy balance, such as He, Ca II, and Mg II. This involves expanding the treatment of the RT to calculate the radiation field utilizing the 3D short-characteristics solver.

We aim to obtain simulated hydrogen population distributions that more closely reflect conditions in the real solar atmosphere, where the heating or cooling due to the radiation field is consistent with that used to treat the atomic populations. In the upper chromosphere, the inclusion of RT enables the simulation of important processes, whereby photons emitted from the transition region and from shock fronts are absorbed by the upper chromosphere. In the lower chromosphere, using radiative rate coefficients from RT can help align the population with the expected detailed radiative balance. In simulations assuming detailed radiative balance, excitation and de-excitation processes are dominated by collisions, with radiative rate coefficients set to zero, thereby imposing no constraints on population evolution. As a result, assuming detailed radiative balance may lead to population simulation results that fall outside of the detailed radiative balance. To reduce computational cost, a single-frequency approach was employed for treating the transitions, following Golding et al. (2016), in which the detailed spectral line profile is neglected, and both opacity and emissivity are assumed to be uniform across the specified passband.

The theoretical model of the single-frequency approach, along with the numerical method used to solve the RT, is described in Sect. 2. In Sect. 3, the accuracy of the treatment in computing radiative rate coefficients and radiative heating or cooling, its impact on the low atmosphere in simulation, the impact of the solver and the Doppler effect, as well as the influence of time step selection on the simulation results are investigated. Finally, a summary is given in Sect. 4.

2. Method

This section presents the governing equations for NLTE RT in the MURaM code, along with the numerical methods employed to solve them. For a detailed description of other aspects of the MURaM code, please refer to Vögler et al. (2005), Rempel (2017) and Przybylski et al. (2022).

2.1. Time-dependent treatment of hydrogen ionization

The nonequilibrium (NE) treatment of hydrogen ionization for the chromosphere is incorporated from the version of the MURaM code presented in Przybylski et al. (2022). The NE treatment considers five energy levels of the hydrogen atom, as well as protons and H₂. Additionally, H₂⁺ and H⁻ are treated in chemical equilibrium, and non-hydrogen species are treated in LTE. The hydrogen population updates at each time step account for the effects of advection, bound-bound transitions, bound-free transitions, and chemical reactions, in which the contribution from advection, bound-bound transitions and bound-free transitions is expressed as

$$\begin{array}{c}\hfill \frac{\partial n_i}{\partial t}+\mathrm{\nabla }·(n_i\mathit{v})={\displaystyle \sum _{j\ne i}}{n}_j(C_{\mathit{ji}}+R_{\mathit{ji}})-n_i{\displaystyle \sum _{j\ne i}}(C_{\mathit{ij}}+R_{\mathit{ij}}),\end{array}$$(1)

where n represents number density, i and j denote energy levels, C is the collisional rate coefficient, and R is the radiative rate coefficient. The problem is split, with the population advection calculated using the partial donor cell method presented in Zhang et al. (2019).

The temperature- and electron-dependent collisional rates were calculated using approximate cross-sections for hydrogen-electron collisions from Johnson (1972), with the formula referenced in Leenaarts et al. (2007). The radiative rates for the Balmer, Paschen, and Brackett series were calculated following the prescription presented in Sollum (1999), where the radiation field is expressed as the Planck function of a defined radiation temperature. The Lyman series were treated with a different method, which is introduced later.

The radiative cooling term is given by

$$\begin{array}{c}\hfill Q_{\mathrm{rad}}=Q_{\mathrm{phot}}+Q_{\mathrm{lines}}+Q_{\mathrm{thin}}+Q_{\mathrm{back}},\end{array}$$(2)

where Q_phot represents the contribution from photospheric radiation, Q_lines accounts for the influence of strong chromospheric spectral lines, Q_thin corresponds to coronal optically thin losses and is applicable in regions where temperatures exceed 10 000 K and pressures fall below a user-defined threshold, and Q_back represents the back-heating of chromospheric plasma by coronal radiation. The photospheric radiation was derived using multigroup RT, accounting for the effects of scattering (Ludwig 1992). The radiative heating or cooling contribution from spectral lines Q_lines = Q_H + Q_Mg + Q_Ca includes the contributions from strong spectral lines of hydrogen, magnesium, and calcium. The terms Q_Mg and Q_Ca were computed following the recipe of Carlsson & Leenaarts (2012), in which only radiative losses are considered. The cooling rate in this recipe is expressed as a function of a defined column mass, electron number density, temperature, and element number density. The calculation method for Q_H is presented later. The coronal optically thin losses, Q_thin = −n_en_HΛ(T), is a function of the electron number density, n_e, hydrogen number density, n_H, and a temperature-dependent curve, Λ(T), obtained from the CHIANTI atomic database (Landi et al. 2012). The back-heating by coronal radiation, Q_back, was calculated using RT, following the approach presented in Carlsson & Leenaarts (2012). In this approach, the emissivity, Q_thin/4π, and the opacity at the ionization edge of helium are incorporated into the RT calculations.

In the previous version of the code presented in Przybylski et al. (2022), the Lyman series is treated assuming detailed balance in radiative transitions, with R_ij = R_ji = 0 applied when updating populations. The heating or cooling rate from hydrogen lines (Q_H) is calculated using the recipe of Carlsson & Leenaarts (2012), i.e., the same approach as is used for magnesium (Q_Mg) and calcium (Q_Ca). The radiative detailed balance assumption deviates significantly from real conditions in the upper chromosphere. Therefore, in this extension of the code, the assumption is discarded when treating the Lyman series. Instead NLTE RT is explicitly solved for the two strongest lines in Lyman series: Lyman-α and Lyman-β. The radiative rates obtained from the RT solution are used to update hydrogen populations, and Q_H is the total radiative gain or loss from the two lines. Other Lyman series lines are much weaker than the two strongest lines, and their impact on chromospheric evolution is relatively minor (Carlsson & Leenaarts 2012). To reduce numerical effort, they are still treated under the detailed radiative balance assumption.

2.2. Theoretical model of NLTE radiative transfer in the Lyman series

Radiative transfer solves the time-independent, frequency- and direction-dependent transport equation

$$\begin{array}{c}\hfill \frac{\mathrm{d}{I}_{\nu }}{\mathrm{d}{\tau }_{\nu }}=S_{\nu }-I_{\nu },\end{array}$$(3)

where I_ν, τ_ν, and S_ν represent the ray intensity, optical depth, and source function in a given direction, respectively. For NLTE RT in the Lyman series, the source function depends primarily on the hydrogen level populations. These spectral lines span a finite bandwidth, and their profiles exhibit significant frequency-dependent variation. Ideally, RT should be solved across multiple frequencies to accurately capture this behavior. However, due to the high computational cost of such calculations, we adopted a single-frequency approximation for the lines, following the approach of Golding et al. (2016).

In the single-frequency approach, spectral line-related photons are assumed to be evenly distributed within a given bandwidth, W_ν, centered at ν₀. The frequency-averaged optical depth, ${\tau }_{\overline{\nu }}$, is given by the integral of the frequency-averaged opacity ${\chi }_{\overline{\nu }}$ along the line of sight:

$$\begin{array}{cc}\hfill {\tau }_{\overline{\nu }}(s)& ={\displaystyle {\int }_{-\infty }^s}{\chi }_{\overline{\nu }}(s\prime )\mathrm{d}s\prime ,\hfill \end{array}$$(4)

$$\begin{array}{cc}\hfill \mathrm{d}{\tau }_{\overline{\nu }}& ={\chi }_{\overline{\nu }}\mathrm{d}s,\hfill \end{array}$$(5)

where s represents the distance along the line of sight. The frequency-averaged opacity and source function are defined, respectively, as follows:

$$\begin{array}{cc}\hfill {\chi }_{\overline{\nu }}& =\frac{h{\nu }_0}{4\pi }\frac{n_{\mathrm{l}}{B}_{\mathrm{lu}}}{W_{\nu }},\hfill \end{array}$$(6)

$$\begin{array}{cc}\hfill S_{\overline{\nu }}& =\frac{h{\nu }_0}{4\pi }\frac{n_{\mathrm{u}}{A}_{\mathrm{ul}}}{W_{\nu }{\chi }_{\overline{\nu }}},\hfill \end{array}$$(7)

where h is Planck’s constant, l and u represent the lower and upper levels of the transition, n_l the population density of the lower level, n_u the population density of the upper level, B_lu the Einstein coefficient for radiative excitation, and A_ul the Einstein coefficient for spontaneous de-excitation. Stimulated emission is significantly weaker than spontaneous emission for these lines and is therefore neglected in the model. The direction-dependent, frequency-averaged ray intensity, $I_{\overline{\nu }}$, is obtained by substituting ${\chi }_{\overline{\nu }}$ and $S_{\overline{\nu }}$ into the transport equation (Eq. (3)) and solving it. Once the radiation field is known, the heating or cooling associated with this spectral line can be determined. The heating rate is given by

$$\begin{array}{c}\hfill Q=h{\nu }_0(n_{\mathrm{l}}{R}_{\mathrm{lu}}-n_{\mathrm{u}}{R}_{\mathrm{ul}}),\end{array}$$(8)

where the downward radiative rate coefficient is R_ul = A_ul, and the upward radiative rate coefficient depends on the frequency-averaged mean intensity, $J_{\overline{\nu }}$:

$$\begin{array}{c}\hfill R_{\mathrm{lu}}=B_{\mathrm{lu}}{J}_{\overline{\nu }}=\frac{B_{\mathrm{lu}}}{4\pi }{\displaystyle \oint I_{\overline{\nu }}}\mathrm{d}\mathrm{\Omega },\end{array}$$(9)

where Ω is a solid angle.

We adopted line center frequencies of ν₀ = c/(121.5 nm) for the Lyman-α line and ν₀ = c/(102.6 nm) for the Lyman-β line. The bandwidths of the transitions are user-defined parameters, with default values of W_ν = 1.0 × 10¹² Hz for the Lyman-α line and W_ν = 1.25 × 10¹² Hz for the Lyman-β line. The single-frequency approach with the default bandwidths gives accurate radiative rate coefficients in the test, as demonstrated in Section 3.2. The regions of the solar atmosphere in which the NLTE treatment was applied were also user-defined. This was implemented in order to avoid an overlap with cooling from the multigroup losses in the photosphere. Additionally, the solution must be smoothly joined with the LTE photosphere. In order to achieve this, the non-LTE RT is solved in regions where p < P_NLTE-RT, and the heating or cooling rates and radiative rate coefficients are applied in regions where p < P_RT-USE. By default, P_NLTE-RT = 1 × 10⁵ dyn cm⁻² and P_RT-USE = 1 × 10⁴ dyn cm⁻², corresponding to typical pressures in the photosphere (Avrett & Loeser 2008).

2.3. Numerical radiative transfer

The implementation of numerical RT for the Lyman lines was guided by the principles of simplicity, computational efficiency, and robustness. For typical MURaM radiative MHD simulations, extreme events including strong shocks and eruptive events up to the scale of small flares can occur. The radiation field and populations must be 100% convergent, and radiative instabilities that can feed back into the energy through heating or cooling must be avoided.

The MURaM code utilizes two sets of meshes: one for numerically solving the MHD equations (MHD mesh) and the other for RT equations (RT mesh). Both meshes are uniform and Cartesian. The grid points of the RT mesh are positioned at the center of the cuboids formed by adjacent grid points of the MHD mesh (left panel of Fig. 1). The populations updated and stored in the MHD mesh are used to solve NLTE RT, with the goal of obtaining the radiative heating or cooling rate (Q) and radiative rate coefficients (R_lu and R_ul) at the MHD mesh. The steps are as follows:

Fig. 1. Schematic illustration of the mesh concerning RT. The left panel shows the computational meshes employed for solving the governing equations in the MURaM code, while the right panel illustrates the relationship between the radiation ray and the mesh structure. The MHD mesh is employed for solving the MHD equations, whereas the RT mesh is used for solving the RT equation. Both meshes are uniform and Cartesian, with the RT grid points positioned at the centers of the cells formed by adjacent MHD grid points. Information exchange between the two meshes is accomplished through bi-linear or tri-linear interpolation. Rays are involved in solving the RT equation using the short-characteristics method, in which the ray intensity at the grid points (o) is the quantity being solved for in MURaM. The solution involves the source functions at the upwind (u), downwind (d), and grid points, along with the ray intensity at the upwind point. The rays typically do not intersect with the grid vertices. Linear or bi-linear interpolation was employed to compute the source function, opacity, and ray intensity at the upwind and downwind points.

1)

• Obtain populations (n_l and n_u) at the RT mesh from the MHD mesh via bi-linear or tri-linear interpolation.

• 2)

  Calculate the opacity (χ) and source function (S) at the RT mesh.

• 3)

  Solve the RT equation to obtain J at RT mesh.

• 4)

  Compute Q at the RT mesh.

• 5)

  Interpolate Q and J to the MHD mesh.

• 6)

  Calculate R_lu and R_ul from J, n_l, and n_u at the MHD mesh.

The RT mesh was introduced to enhance robustness. When Q and J are interpolated from the RT mesh to the MHD mesh, their spatial distributions become smoother, thereby reducing the likelihood of numerical instabilities, which could compromise the stability of the simulation. In addition, the existing multigroup RT for photospheric radiation is solved on the RT mesh. Solving the Lyman RT on the same mesh allows us to reuse many of the existing functions, thereby preventing unnecessary code expansion and keeping the implementation streamlined. A four-step Runge-Kutta-like method was used for time integration in solving the MHD equations (Jameson 2017). The complete form of the NE equation of state, including an update of the rate equation, was solved only for the first step. A simplified solution, which conserves charge, hydrogen nucleus number, and energy was solved for the other steps.

The RT problem was solved using the short characteristics method, employing the type A quadrature scheme described in Carlson et al. (1963) for angular integration. As was noted earlier, the RT is solved on the RT mesh, where our objective is to determine J at each grid point. For this purpose, we selected 24 ray directions (three angular bins per quadrant) following Carlson et al. (1963), sequentially computed the ray intensity spatial distribution for each direction, and finally averaged these intensity distributions to obtain J. The computation of the intensity distribution of a direction requires an iterative solution: at each iteration, the intensity at every grid point is updated based on the current values at adjacent points, and the process continues until convergence is reached. The short characteristics method was employed to perform these intensity updates. The details of the short characteristics method and the iterative scheme are described below.

In the short characteristics method, the local ray intensity is updated based on the intensities at nearby grid points. The equation of ray intensity can be expressed as

$$\begin{array}{c}\hfill I_i({\tau }_{\mathrm{o},i})=I_i({\tau }_{\mathrm{u},i}){\mathrm{e}}^{-({\tau }_{\mathrm{o},i}-{\tau }_{\mathrm{u},i})}+{\displaystyle {\int }_{{\tau }_{\mathrm{u},i}}^{{\tau }_{\mathrm{o},i}}}{S}_i(t){\mathrm{e}}^{t-{\tau }_{\mathrm{o},i}}\mathrm{d}t,\end{array}$$(10)

where i represents the direction index, τ_o, i is the optical depth of the i-th ray at the grid point, and τ_u, i is the optical depth at the upwind point. The upwind point, u, is different for different directions. Here, the second-order scheme BESSER is used as the default solver to numerically integrate the source function, assuming a monotonic quadratic Bézier variation, given by

$$\begin{array}{c}\hfill {\displaystyle {\int }_{{\tau }_{\mathrm{u}}}^{{\tau }_{\mathrm{o}}}}S(t){\mathrm{e}}^{t-{\tau }_{\mathrm{o}}}\mathrm{d}t=\mathrm{\Delta }I={\mathrm{\Psi }}_{\mathrm{u}}S({\tau }_{\mathrm{u}})+{\mathrm{\Psi }}_{\mathrm{o}}S({\tau }_{\mathrm{o}})+{\mathrm{\Psi }}_{\mathrm{c}}{c}_{\mathrm{u}},\end{array}$$(11)

where c_u denotes the source function at a specified control point, and Ψ_u, Ψ_o and Ψ_c are interpolation weights. As a second-order scheme, the value of c_u depends on S(τ_u), S(τ_o), and downwind point source function S(τ_d), with the methods for computing c_u and the weights are detailed in Štěpán & Trujillo Bueno (2013). The spatial relationship between rays and points is depicted in Fig. 1. As is shown in the figure, the ray does not coincide with the mesh lines, so the upwind and downwind points of the ray generally do not coincide with the grid points. This requires interpolation to determine the source function and ray intensity at these points. Alternatively, a first-order linear solver (Kunasz & Auer 1988) is available for source function integration.

The solution of ray intensity was obtained using the Gauss–Seidel iterative method, in which the updated values are used within each iteration, based on the following expression:

$$\begin{array}{c}\hfill I^{\mathrm{new}}({\tau }_{\mathrm{o}})=\mathrm{\Delta }I+I^{\mathrm{old}/\mathrm{new}}({\tau }_{\mathrm{u}}){\mathrm{e}}^{-\mathrm{\Delta }{\tau }_{\mathrm{ou}}},\end{array}$$(12)

where ΔI is given by Eq. (11) and is computed in advance prior to the iteration, and Δτ_ou is computed in advance as well. In each iteration step, we updated I at all grid points in the region of interest according to Eq. (12). The code employs MPI parallel computing, in which the simulation box is divided into multiple regions, with each processor responsible for managing the data within one region. When computing the ray intensity at grid points located at the boundary of a processor’s domain, the upwind point may fall outside its domain, depending on the ray direction, leaving the processor without the necessary intensity data. Consequently, in each iteration, cross-processor communication is required to obtain the ray intensities at points outside the domain. The grid points are swept from the upwind side. If the upwind point falls outside the domain, I(τ_u) is assigned the value from the preceding iteration (I^old(τ_u)); otherwise, the updated value is employed (I^new(τ_u)). For further details on data exchange related to ray intensity, refers to Vögler (2003). The iteration process continues until the error falls below a predefined threshold at all grid points in the region of interest:

$$\begin{array}{c}\hfill max\left(\frac{I^{\mathrm{new}}-I^{\mathrm{old}}}{max(I_{\min },I^{\mathrm{old}})}\right)<5\times {10}^{-5},\end{array}$$(13)

where I_min is a defined small value used to prevent division by zero. Radiative transfer is solved in the regions where Δτ_ou ≥ Δτ_min, where Δτ_min is a user-specified parameter that allows control over whether RT is also solved in the optically thin corona. In regions where Δτ_ou < Δτ_min, the ray intensity I is set to zero, as radiation absorption can be neglected due to the low opacity.

The mean intensity required in Eq. (9) was computed as the average of the ray intensities over 24 directions:

$$\begin{array}{c}\hfill J={\displaystyle \sum _i}{\omega }_i{I}_i,\end{array}$$(14)

where ω_i is the weight assigned to direction i given by the quadrature, in this case Carlson et al. (1963). Radiative heating contributed by the transition was calculated from J in regions where the optical thickness is relatively small, and from the divergence of the radiative energy flux (∇ ⋅ F) in optically thick regions, similar to the approach adopted for multigroup RT (Bruls et al. 1999). The heating rate is given by

$$\begin{array}{c}\hfill Q=h{\nu }_0(n_{\mathrm{l}}{B}_{\mathrm{lu}}J-n_{\mathrm{u}}{A}_{\mathit{ul}}){\mathrm{e}}^{-\mathrm{\Delta }\tau /{\tau }_0}-\mathrm{\nabla }·F(1-{\mathrm{e}}^{-\mathrm{\Delta }\tau /{\tau }_0}),\end{array}$$(15)

where τ₀ = 0.1, Δτ = χmin(Δx, Δy, Δz), χ is the opacity, and Δx, Δy, and Δz are the grid cell sizes in the three spatial directions. In optically thick regions, the terms n_lB_luJ and n_uA_ul are nearly equal, making round-off errors significant. In such cases, computing the heating rate from the radiative energy flux yields more accurate results.

3. Results

3.1. The influence of time step selection on the accuracy of simulation results

In this subsection, we investigate the influence of time step selection on the evolution of the level populations in simulations. Radiative excitation and de-excitation are processes that occur on very short timescales. For example, the radiative de-excitation timescale for the Lyman-α transition is given by 1/R₂₁ ≈ 1/A₂₁ ≈ 2 ns. Resolving the radiative de-excitation process in simulations would require a time step of Δt < 1 ns, which is computationally prohibitive for a 3D MHD simulation, where timesteps are typically on the order of 10⁻³–10⁻¹ seconds. However, due to this extremely short timescale, the radiative relaxation of the level populations occurs rapidly, such that the population distribution can be considered to be in balance, dependent on the local radiation field, at each moment. To capture the evolution of the level populations, a time step on the timescale of radiation field variations is required, rather than one based on the nanosecond timescales of radiative excitation and de-excitation. The timescale of radiation field variations should be comparable to that of MHD processes, such as convection, as population evolution dominated by collisions occurs over longer timescales. According to Carlsson & Stein (2002), the ionization/relaxation timescale associated with chromospheric shocks range from 10 to 10³ seconds, increasing dramatically in the lower, cooler layers of the chromosphere (10³–10⁵ seconds). Consequently, the evolution of the level populations can, in principle, be captured using the MHD time step in the simulation. To verify this, we evaluated the convergence of the simulation results by varying the time step size.

The survey comprises four 2D simulation cases, all of which begin with the same initial condition. In two of these cases, the time step is determined by the Courant–Friedrichs–Lewy (CFL) condition, with CFL numbers set to 0.6 and 0.1, respectively. The resulting average time steps are approximately 10⁻² s and 2 × 10⁻³ s. In the remaining two cases, the time step is explicitly specified as 10⁻⁴ s and 10⁻⁵ s, respectively. The initial condition originates from the output generated by the MURaM chromospheric extension. The simulation domain spans 8.192 Mm × 8.192 Mm, which is covered by 256 × 256 grid points. This corresponds to a resolution of 32 km both horizontally and vertically. The lower boundary is located in the convection zone, the upper boundary is in the corona. A Lyman-α bandwidth of W_ν = 1 × 10¹² Hz and a Lyman-β bandwidth of W_ν = 1.25 × 10¹² Hz are adopted, as described in Sec. 2.2. Radiative transfer is solved for both the corona and the chromosphere, with the parameter Δτ_min = 10⁻⁸ (see Sect. 2.3) used in MURaM.

The temperature and n₂ population number density distributions at t = 60 s are compared in Fig. 2. The four cases exhibit similar temperature and population distributions, as is shown in the figure. The profiles in Fig. 2 (e)-(f) show a slight change in the results when the CFL number is reduced from 0.6 to 0.1, with no noticeable effect from further reductions in time step size. As expected, the population evolution is well captured when using the MHD time step. A balance solution for n₂ in the Δt = 10⁻⁵ case is shown in Fig. 2f, which is given by

Fig. 2. Temperature and n2 population distributions at t = 60 s for the four simulation cases. The average time steps of the CFL = 0.6 and CFL = 0.1 cases are approximately 10−2 s and 2 × 10−3 s, respectively. Panels (a) and (b) display the temperature and n2 distributions for the case with CFL = 0.6. Panels (c) and (d) show the corresponding regions (outlined in panels (a) and (b)) for all four cases, respectively. Panels (e) and (f) present 1D profiles of temperature and n2, respectively, extracted along the dashed lines in panels (a) and (b), for all four cases. In panel (f), an additional n2 profile computed from the balance solution for the Δt = 10−5 s case is also shown for reference.

$$\begin{array}{c}\hfill n_2=\frac{{\sum }_{j\ne 2}{n}_j(C_{j2}+R_{j2})}{{\sum }_{j\ne 2}(C_{2j}+R_{2j})}.\end{array}$$(16)

The simulated n₂ distribution is in good agreement with the balance solution, as expected.

The temporal evolution of level populations can be classified into fast and slow evolution (Carlsson & Stein 2002; Judge 2005). The slow evolution corresponds to changes occurring under approximate balance, a condition typical in the chromosphere. In contrast, the fast evolution describes the rapid establishment of approximate balance from an out-of-balance state, driven by the radiative transitions. In the code, the time evolution is solved using the Newton–Raphson method, which captures the faster rates, and allows for accurate modeling of the slow rates, playing a crucial role in the accurate modeling of population dynamics (Przybylski et al. 2022). The iterative form of the rate equation Eq. (1) is given by

$$\begin{array}{c}\hfill f_i=\frac{n_i^{t_0+\mathrm{\Delta }t}}{n_i^{t_0}}-\frac{\mathrm{\Delta }t}{n_i^{t_0}}[{\displaystyle \sum _{j\ne i}}{n}_j(C_{\mathit{ji}}+R_{\mathit{ji}})-n_i{\displaystyle \sum _{j\ne i}}(C_{\mathit{ij}}+R_{\mathit{ij}})]-1=0,\end{array}$$(17)

where the advection term is neglected. In the case of slow evolution – that is, when the timescale of change is much larger than the time step, Δt – the behavior of this method is consistent with that of the explicit method. For fast evolution, where the timescale is much shorter than Δt, the populations will reach an approximate balance given by

$$\begin{array}{c}\hfill {\displaystyle \sum _{j\ne i}}{n}_j(C_{\mathit{ji}}+R_{\mathit{ji}})-n_i{\displaystyle \sum _{j\ne i}}(C_{\mathit{ij}}+R_{\mathit{ij}})\approx 0,\end{array}$$(18)

within one time step, in agreement with the actual physical process. For more detailed study and discussion of population evolution, refer to Carlsson & Stein (2002) and Judge (2005).

3.2. Accuracy of the single-frequency approach for the Lyman-α and β transitions

To assess the accuracy of the single-frequency approach applied to the Lyman lines, we compare its results with reference solutions obtained using a more detailed multifrequency approach. Using the output of the CFL = 0.6 case presented in Sect. 3.1, we computed the radiative heating and cooling rates, as well as the radiative rate coefficients, with both methods.

The single-frequency results were obtained from the MURaM code, while the multifrequency reference results were generated using the Lightweaver framework (Osborne & Milić 2021). In the single-frequency approach, only the hydrogen level population data were used in the computation with a specified line width. In the multifrequency calculations, temperature data were additionally used to model spectral line broadening effects. By setting the velocity to zero, the Doppler effect is neglected. Sect. 3.4 demonstrates that it has an insignificant impact on the radiative rates. Therefore this simplification has a negligible impact on our results. Lightweaver solves the RT equation using the short-characteristics approach for 1D and 2D problems. It accounts for partial frequency redistribution (PRD) in spectral lines and computes the specific intensity at each frequency, spatial point, and direction. The Lyman-α (line center at 121.5684 nm) multifrequency results are computed using 99 frequency points within a wavelength range of 120.8384 – 122.2983 nm, while the Lyman-β (line center at 102.5733 nm) results are obtained from 48 frequency points within a wavelength range of 102.3167 – 102.8299 nm. For consistency, the RT equation was solved using the second-order solver BESSER, the same solver used in the MURaM code. A 24-bin spherical space discretization was employed in solving the equation, as in the MURaM code, but with a different quadrature provided by Štěpán et al. (2020). The multifrequency radiative heating or cooling rates and radiative rate coefficients were first computed on the RT mesh, then interpolated onto the MHD mesh, and subsequently compared to the single-frequency results at the MHD mesh, which were computed using the method described in Sec. 2.3.

The results of the single-frequency approach and the multifrequency reference calculations are presented and compared in Figs. 3 and 4. As is shown in Figs. 3c–f, the spatial distribution of upward radiative rate coefficients (R_lu) of two Lyman lines are well reproduced by the single-frequency approach. The point-to-point comparison in Figs. 3g–h show that the rate coefficient amplitudes are in excellent agreement across most of the locations, with only minor deviations observed. These minor discrepancies are located at the upper chromosphere and the transition region, where the medium transitions from optically thick to optically thin. Most regions of the chromosphere are optically thick in Lyman-α, allowing the single-frequency approach to yield accurate upward rate coefficients in these areas. The amplitude and spatial distribution of the heating and cooling effects are also well captured, as is shown in Fig. 4. In terms of correlation with the reference, the heating or cooling effect is not captured as accurately as the radiative rate coefficient. Throughout most of the chromosphere, the rates of photon emission and absorption are nearly balanced in the transitions, such that small deviations in the radiative rate coefficient can lead to proportionally larger variations in radiative heating or cooling. The strongest radiative heating and cooling effects occur at the transition region boundary. The detailed radiative balance of the lines breaks down there, refer to Fig. 6 and Sect. 3.3.

Fig. 3. Comparison of the upward radiative rate coefficients obtained by MURaM (using the single-frequency approach) and reference (using the multifrequency approach by Lightweaver). Panels (a) and (b) show the atmospheric density and temperature distributions, respectively. Panels (c) and (d) display the Lyman-α and Lyman-β rate coefficients in the cool region (T <  = 0.2 MK) computed by MURaM, while panels (e) and (f) present the corresponding coefficients from the reference. Panels (g) and (h) display the results of the point-to-point amplitude comparison, with the horizontal axis representing the values from MURaM and the vertical axis corresponding to the amplitudes from the reference. The points are colored according to their corresponding color in the MURaM result panels (c, d), so the area where the points are located can be inferred from the color. Since the highly ionized corona is not the region of interest for NLTE RT, the rate coefficients at high temperatures (T > 0.2 MK) were manually set to zero to emphasize the comparison results in the lower atmosphere. The time corresponding to the simulated data is t = 1 s.

Fig. 4. Comparison of the upward radiative rate coefficients obtained by MURaM and Lightweaver (reference). Panels (a) and (b) show the heating effects contributed by the Lyman-α and Lyman-β lines from MURaM, while panels (c) and (d) present the corresponding heating effects in the reference. Panels (e) and (f) give the point-to-point amplitude comparison, where the horizontal axis corresponds to the values obtained from MURaM and the vertical axis to the amplitudes from the reference. The points are colored based on their sign in both the MURaM and reference results, with red and pink indicating heating, and blue and cyan indicating cooling. The time corresponding to the simulated data is t = 1 s.

The spatial distribution of the heating and cooling effects of the Lyman-β line is similar to that of the Lyman-α line, but its amplitude is two orders of magnitude weaker. The downward radiative rate coefficients, R_ul, are very close to the constants A_ul in the lower atmosphere, as the radiative de-excitation of the Lyman lines in this region primarily governed by the spontaneous de-excitation process. Therefore, the results are not compared here.

3.3. Effects of the NLTE radiative transfer treatment of the Lyman series

We now evaluate the impact of employing NLTE RT in the treatment of the Lyman series on chromospheric simulations. We conducted simulations using both the updated version of the MURaM code and the previous version, as it is described in Przybylski et al. (2022), and compared the results. As is described in Sect. 2, the previous version assumes detailed radiative balance for the Lyman series by setting the radiative rates to zero in the implementation, and computes radiative cooling from hydrogen using the recipe from Carlsson & Leenaarts (2012). In contrast, the updated version solves the NLTE RT equations to calculate the radiative rate coefficients for the Lyman-α and Lyman-β transitions, as well as the associated radiative heating and cooling rates. More details are provided in Sect. 2.1. Both simulations use the same initial condition, as used in Sect. 3.1. A CFL number of 0.6 is used in both simulations to limit the timestep.

The chromosphere produced by the two versions of the code differ significantly, particularly in terms of temperature and hydrogen level populations. Fig. 5 presents a comparison of the temperature distributions and hydrogen level population fraction distributions from the two simulations at t = 20 s. In terms of temperature, the updated version of the code produces a hotter upper chromosphere. Given that the two versions employ different methods to compute radiative cooling, differences in temperature are expected. Additionally, differences in the population states will change the amount of internal energy held in ionization and excitation, changing the temperature. From a physical perspective, the previous version accounts only for radiative cooling and neglects radiative heating related to the absorption of Lyman-α at the wavefront (Carlsson & Leenaarts 2012). Therefore, it is reasonable that the updated version, which includes both processes, yields a higher chromospheric temperature. The radiative cooling due to the Lyman continuum and H-α have not been included in the update version, which may also contribute to the observed temperature differences.

Fig. 5. Comparison of the temperature and population distributions obtained by two versions of the MURaM code. The left column shows results from the previous version, while the right column presents results from the updated version. The first row displays the temperature distribution, while the second, third, fourth, and fifth rows show the population fractions of level 1, level 2, level 3, and protons, respectively, where np is the proton number density and nH is the total hydrogen number density. To compare the extended height of the incompletely ionized region, a reference line is provided in panels (c), (d), (i) and (j). The previous version assumes Rlu = Rul = 0 for the Lyman series and computes radiative cooling using an empirical recipe, whereas the updated version calculates the Lyman-α and β lines rate coefficients and the associated radiative heating and cooling through solving the RT equation. The time corresponding to the simulated data is t = 20 s

In terms of populations, noticeable differences between the two versions are observed. The ground level population fraction and ionization fraction distributions (Figs. 5 c, d, i, j) show that the incompletely ionized region reaches higher altitudes in the version incorporating Lyman-α radiation transfer (updated version), which aligns with the conclusions of Golding et al. (2016). The spatial distributions of the populations of the upper level of the transitions, n₂ and n₃, are further decoupled from the temperature and exhibit smoother profiles in the updated version. In the previous version, collision processes dominate excitation and de-excitation, since the radiative rate coefficients for the lines are set to zero. Since the collisional rate coefficients are temperature-dependent, some structures in the temperature distribution also appear in the distributions of n₂ and n₃. In the updated version, the inclusion of radiative de-excitation processes suppresses the occurrence of high n₂ and n₃ fractions in the chromosphere, resulting in reduced peak values. The inclusion of RT in the updated version introduces a new channel for energy exchange between different regions, potentially contributing to the differences in level populations. The populations obtained in the updated version are closer to satisfying the detailed radiative balance condition, as demonstrated in Fig. 6. In the results from the updated version, only the region near the wavefront deviates from detailed radiative balance, while the rest of the chromosphere can be approximated as being in detailed radiative balance. Despite adopting the assumption of detailed radiative balance, the previous version produces simulation outcomes that deviate markedly from detailed radiative balance condition.

Fig. 6. Comparison of the degree to which the population distributions provided by the two versions of the MURaM code conform to the radiative balance. Radiative upward and downward transition rates are computed from the populations shown in Fig. 5 using the Lightweaver code and are illustrated in panels (a)-(d). Panels (e) and (f) show the degree of deviation from the detailed radiative balance condition in the results from the two versions, with the calculation formula provided in the title. The time corresponding to the simulated data is t = 20 s. The rate coefficients at high temperatures (T > 0.2 MK) were manually set to zero to emphasize the results in the lower atmosphere.

Table 1 shows the impact of activating this NLTE RT module on simulation speed. The computational time increases by approximately 9%–21%, depending on the number of spectral lines included and whether RT is resolved in the optically thin corona. Compared with Case 1, the increased computation time in the other cases is primarily attributable to the solution of the Lyman RT problem. There might remain potential to improve the computational efficiency of this module.

3.4. Impacts of the solver and the Doppler effect

Two schemes, the BESSER and the linear solver, were implemented for the module, as is detailed in Section 2.3. To evaluate their differences, we conducted a test comparing the simulation results produced by the two solvers. In the test, simulations were performed using the BESSER and linear solvers, respectively. Both the Lyman-α and Lyman-β lines were activated in the simulations, with Δτ_min set to 10⁻⁸. Fig. 7 presents the radiative rate coefficient (R₁₂) and population distribution (n₂) at t = 100 s for both simulations. As is shown in panel e, slight differences can be observed in the radiative rates obtained by the two solvers. However, this difference in radiative rate does not lead to significant long-term differences in the population distribution, as indicated by the nearly identical population profiles in panel f. There may be an underlying mechanism that drives the convergence of the population evolution.

Fig. 7. Comparison of the simulation results generated by two solvers. The radiative rate-coefficient distributions (R12) produced by the BESSER and linear solvers are shown in panels (a) and (c), respectively, while panel (e) compares the distributions extracted along the dashed line depicted in panels (a) and (c). Similarly, the level-population distributions (n2) obtained from the two solvers are presented in panels (b) and (d), respectively, with panel (f) showing the corresponding comparison along the dashed line. All results correspond to t = 100 s. The rate coefficients at high temperatures (T > 0.2 MK) were manually set to zero to emphasize the results in the lower atmosphere.

In the single-frequency treatment of the Lyman lines, the Doppler effect is neglected. To assess its impact on the Lyman-α line, we conducted a test using Lightweaver to calculate the upward radiative rate coefficient (R₁₂) both with and without the Doppler effect, and subsequently compared the results. The calculations are based on the data produced by the BESSER solver, as described earlier in this chapter, using the hydrogen level populations, temperature and velocity values in the computation. For the calculations that neglect the Doppler effect, the plasma velocity was set to zero, whereas the calculations accounting for the Doppler effect employed the velocity values obtained from the simulation output. Fig. 8 presents the simulation results, including temperature, vertical velocity, and the computed rate coefficients from both approaches. The results show that the difference in rate coefficients between the two methods is on the order of 10⁻³, indicating that the Doppler effect can be safely neglected in this context.

Fig. 8. Comparison of the Lyman-α upward radiative rate coefficient (R12) computed with and without the Doppler effect. Panels (a) and (b) present the simulation outputs for temperature and vertical velocity, respectively. Panels (c) and (d) show the rate coefficients calculated with and without accounting for the Doppler effect. Panel (e) compares the rate coefficients along the dashed line indicated in panels (c) and (d), while panel (f) presents the point-to-point comparison, with the evaluation formula provided above the panel. The rate coefficients at high temperatures (T > 0.2 MK) were manually set to zero to emphasize the results in the lower atmosphere.

4. Discussion and conclusion

To more accurately simulate chromospheric dynamics, we incorporated an NLTE RT module into the MURaM code, and used it to compute radiative rate coefficients and radiative heating or cooling associated with the Lyman-α and β lines. In this module, the time-dependent RT equation is solved based on the evolving population densities within the simulation. To reduce computational cost, a single-frequency approach was employed to treat these spectral lines, assuming uniform light intensity and opacity within a specified bandwidth for each line. To evaluate the accuracy of the single-frequency approach, the radiative rate coefficients and associated heating and cooling effects computed by MURaM were compared with reference solutions provided by the Lightweaver framework. The comparison results demonstrate that the radiative rate coefficients are accurately reproduced, and the radiative heating and cooling effects are also well captured. After incorporating this module, the chromospheric simulation results from MURaM show significant changes, including a hotter upper chromosphere, an upward extension of incomplete ionization region, and hydrogen level populations more consistent with detailed radiative balance across the chromosphere. The module enables simulations to accurately capture population evolution under MHD time steps, demonstrating its practicality for high-resolution multidimensional applications.

In future work, 3D simulations will be employed to further validate the module. Simulations will be performed using both the previous and updated versions, and the resulting atmospheric outputs (e.g., temperature, density, and velocity), as well as the structure of the jets, will be compared to evaluate the impact of the new module. The inclusion of RT for the Lyman-α and β lines aims to enhance the accuracy of hydrogen population modeling in the chromosphere. To assess whether this objective has been successfully achieved, the forward-modeled results will be compared with imaging and spectroscopic observations in the Lyman-α, Lyman-β, and Hα bands. The MULTI3D tool will be used for imaging and spectroscopic synthesis (Leenaarts & Carlsson 2009). On the observational side, data from a variety of instruments can be used for comparison; for example, data from the Solar Ultraviolet Measurements of Emitted Radiation on board the Solar and Heliospheric Observatory (SOHO/SUMER, Wilhelm et al. 1995), the Extreme Ultraviolet Variability Experiment on board the Solar Dynamics Observatory (SDO/EVE, Woods et al. 2012), the Very high angular resolution ultraviolet telescope (VAULT2.0, Vourlidas et al. 2016), the Spectral Imaging of the Coronal Environment on board the Solar Orbiter (SolO/SPICE, SPICE Consortium 2020), the Lyman-alpha Solar Telescope on board the Advanced Space-based Solar Observatory (ASO-S/LST, Li et al. 2019), and the Hα Imaging Spectrograph on board the Chinese Hα Solar Explorer (CHASE/HIS, Liu et al. 2022).

In the single-frequency approach to RT, the bandwidth parameter (W_ν) is a key factor. The default bandwidths for the Lyman-α and Lyman-β lines were selected based on a limited parameter survey, based on how well the resulting upward radiative rate coefficients and heating or cooling effect matched those from Lightweaver. The test case is described in Sect. 3.2. Our choice of Lyman-α line bandwidth differs from that of Golding et al. (2016), who employed a bandwidth of W_ν = 5 × 10¹¹ Hz. This value is interpreted as the frequency range corresponding to half the maximum value of the Doppler-broadened line profile at a temperature of 10 000 K. In our test, a Lyman-α bandwidth of around W_ν = 1 × 10¹² Hz yields a good agreement with reference results. This value may physically correspond to the thermal broadening at a specific location within the wavefront or to the average thermal broadening across the wavefront region. Bandwidth selection remains an area for improvement in our work. With a more extensive parameter survey, it may be possible to identify bandwidths that simultaneously capture accurate radiative rate coefficients and radiative heating and cooling effects in the upper chromosphere. The influence of simulation dimensionality and resolution on the optimal bandwidth also requires further investigation.

With regard to radiative cooling from hydrogen, only the contributions from the Lyman-α and β lines are considered in the updated version of the MURaM code. According to Carlsson & Leenaarts (2012), these two spectral lines dominate the cooling across most of the chromospheric temperature range, with Lyman-α being the primary contributor. However, below temperatures of 7000 K, Hα becomes the dominant contributor to radiative cooling. In addition, the Lyman continuum also contributes significantly to chromospheric radiative losses. The impact of the absence of the Hα and Lyman continuum on the evolution of chromospheric temperature remains to be investigated. Including RT for these lines is a potential direction for the future development of the MURaM code. The current scheme could easily be extended to helium (Golding et al. 2016), as well as to Mg II and Ca II, to further improve the models.

Acknowledgments

We would like to thank Narayanamurthy Smitha for helpful discussions related to RT and the solar chromosphere. 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 – project WINSUN). 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 and Viper supercomputer systems of the Max Planck Computing and Data Facility (MPCDF) in Garching, Germany. We thank C. Osborne (Lightweaver) and A. Irwin (FreeEoS).

References

All Tables

All Figures

Fig. 1. Schematic illustration of the mesh concerning RT. The left panel shows the computational meshes employed for solving the governing equations in the MURaM code, while the right panel illustrates the relationship between the radiation ray and the mesh structure. The MHD mesh is employed for solving the MHD equations, whereas the RT mesh is used for solving the RT equation. Both meshes are uniform and Cartesian, with the RT grid points positioned at the centers of the cells formed by adjacent MHD grid points. Information exchange between the two meshes is accomplished through bi-linear or tri-linear interpolation. Rays are involved in solving the RT equation using the short-characteristics method, in which the ray intensity at the grid points (o) is the quantity being solved for in MURaM. The solution involves the source functions at the upwind (u), downwind (d), and grid points, along with the ray intensity at the upwind point. The rays typically do not intersect with the grid vertices. Linear or bi-linear interpolation was employed to compute the source function, opacity, and ray intensity at the upwind and downwind points.

In the text

Fig. 2. Temperature and n2 population distributions at t = 60 s for the four simulation cases. The average time steps of the CFL = 0.6 and CFL = 0.1 cases are approximately 10−2 s and 2 × 10−3 s, respectively. Panels (a) and (b) display the temperature and n2 distributions for the case with CFL = 0.6. Panels (c) and (d) show the corresponding regions (outlined in panels (a) and (b)) for all four cases, respectively. Panels (e) and (f) present 1D profiles of temperature and n2, respectively, extracted along the dashed lines in panels (a) and (b), for all four cases. In panel (f), an additional n2 profile computed from the balance solution for the Δt = 10−5 s case is also shown for reference.

In the text

Fig. 3. Comparison of the upward radiative rate coefficients obtained by MURaM (using the single-frequency approach) and reference (using the multifrequency approach by Lightweaver). Panels (a) and (b) show the atmospheric density and temperature distributions, respectively. Panels (c) and (d) display the Lyman-α and Lyman-β rate coefficients in the cool region (T <  = 0.2 MK) computed by MURaM, while panels (e) and (f) present the corresponding coefficients from the reference. Panels (g) and (h) display the results of the point-to-point amplitude comparison, with the horizontal axis representing the values from MURaM and the vertical axis corresponding to the amplitudes from the reference. The points are colored according to their corresponding color in the MURaM result panels (c, d), so the area where the points are located can be inferred from the color. Since the highly ionized corona is not the region of interest for NLTE RT, the rate coefficients at high temperatures (T > 0.2 MK) were manually set to zero to emphasize the comparison results in the lower atmosphere. The time corresponding to the simulated data is t = 1 s.

In the text

Fig. 4. Comparison of the upward radiative rate coefficients obtained by MURaM and Lightweaver (reference). Panels (a) and (b) show the heating effects contributed by the Lyman-α and Lyman-β lines from MURaM, while panels (c) and (d) present the corresponding heating effects in the reference. Panels (e) and (f) give the point-to-point amplitude comparison, where the horizontal axis corresponds to the values obtained from MURaM and the vertical axis to the amplitudes from the reference. The points are colored based on their sign in both the MURaM and reference results, with red and pink indicating heating, and blue and cyan indicating cooling. The time corresponding to the simulated data is t = 1 s.

In the text

Fig. 5. Comparison of the temperature and population distributions obtained by two versions of the MURaM code. The left column shows results from the previous version, while the right column presents results from the updated version. The first row displays the temperature distribution, while the second, third, fourth, and fifth rows show the population fractions of level 1, level 2, level 3, and protons, respectively, where np is the proton number density and nH is the total hydrogen number density. To compare the extended height of the incompletely ionized region, a reference line is provided in panels (c), (d), (i) and (j). The previous version assumes Rlu = Rul = 0 for the Lyman series and computes radiative cooling using an empirical recipe, whereas the updated version calculates the Lyman-α and β lines rate coefficients and the associated radiative heating and cooling through solving the RT equation. The time corresponding to the simulated data is t = 20 s

In the text

Fig. 6. Comparison of the degree to which the population distributions provided by the two versions of the MURaM code conform to the radiative balance. Radiative upward and downward transition rates are computed from the populations shown in Fig. 5 using the Lightweaver code and are illustrated in panels (a)-(d). Panels (e) and (f) show the degree of deviation from the detailed radiative balance condition in the results from the two versions, with the calculation formula provided in the title. The time corresponding to the simulated data is t = 20 s. The rate coefficients at high temperatures (T > 0.2 MK) were manually set to zero to emphasize the results in the lower atmosphere.

In the text

Fig. 7. Comparison of the simulation results generated by two solvers. The radiative rate-coefficient distributions (R12) produced by the BESSER and linear solvers are shown in panels (a) and (c), respectively, while panel (e) compares the distributions extracted along the dashed line depicted in panels (a) and (c). Similarly, the level-population distributions (n2) obtained from the two solvers are presented in panels (b) and (d), respectively, with panel (f) showing the corresponding comparison along the dashed line. All results correspond to t = 100 s. The rate coefficients at high temperatures (T > 0.2 MK) were manually set to zero to emphasize the results in the lower atmosphere.

In the text

Fig. 8. Comparison of the Lyman-α upward radiative rate coefficient (R12) computed with and without the Doppler effect. Panels (a) and (b) present the simulation outputs for temperature and vertical velocity, respectively. Panels (c) and (d) show the rate coefficients calculated with and without accounting for the Doppler effect. Panel (e) compares the rate coefficients along the dashed line indicated in panels (c) and (d), while panel (f) presents the point-to-point comparison, with the evaluation formula provided above the panel. The rate coefficients at high temperatures (T > 0.2 MK) were manually set to zero to emphasize the results in the lower atmosphere.

In the text