RADYN+RH Simulations

• RADYN+RH simulations are a two-stage workflow that couples dynamic radiative-hydrodynamic modeling with detailed NLTE radiative transfer post-processing for spectral synthesis.
• They employ advanced techniques such as multi-level NLTE, partial frequency redistribution, and Stark broadening to diagnose spectral features despite challenges from 1D/1.5D approximations.
• The method is applied to solar and stellar flare spectroscopy, Ellerman-bomb modeling, and photospheric diagnostics, offering insights into energy deposition, line asymmetries, and observational tensions.

Searching arXiv for the cited RADYN/RH papers to ground the article. RADYN+RH simulations are a two-stage radiative-hydrodynamic and radiative-transfer workflow in which RADYN computes a time-dependent flare or impulsively heated atmosphere, and RH or RH 1.5D then synthesizes selected spectra from individual atmospheric snapshots without feeding the spectral solution back into the atmospheric evolution. In this division of labor, RADYN provides the evolving stratification of temperature, density, electron density, velocity, and, in many applications, hydrogen populations, whereas RH provides detailed line formation, source-function diagnostics, contribution functions, formation heights, partial frequency redistribution, and, when configured, polarized transfer. The workflow is used in solar flare spectroscopy, Ellerman-bomb modeling, photospheric line diagnostics, stellar flare continuum synthesis, and multidimensional post-processing of simulation columns [(Zhu et al., 2019); (Pereira et al., 2014)].

1. Core workflow and code division of labor

The defining feature of RADYN+RH simulations is the snapshot-based handoff from a dynamical atmosphere calculation to a separate spectral-synthesis calculation. RADYN solves the time-dependent 1D radiative-hydrodynamic response of an atmosphere to imposed heating, often a non-thermal electron beam. RH then takes those atmospheric states as fixed inputs and solves a more detailed radiative-transfer problem for diagnostics that are not treated adequately in the hydrodynamic stage, or that require different transfer physics. The RH stage is therefore a post-processing step rather than an online coupling.

This structure is stated explicitly in multiple applications. In Na I D$_1$ flare modeling, RADYN snapshots generated at different time steps were used as input atmospheres to RH, and Na I D$_1$ profiles were synthesized for every single time step during a 60 s RADYN run (Kuridze et al., 2016). In Mg II flare-ribbon modeling, RADYN supplied height-dependent temperature, electron density, hydrogen populations, bulk velocity, and microturbulent velocity, while RH recalculated Mg II h, k, and triplet spectra from selected snapshots (Zhu et al., 2019). In flare-induced photospheric velocity studies, the strongest RADYN atmosphere was exported to RH at 1 s cadence over the full 50 s simulation for Fe I 617.3, 630.1, and 630.2 nm synthesis (Monson et al., 2021).

The approach is modular rather than monolithic. RADYN is used for the time-dependent atmospheric state; RH is used as the line-formation engine. This is why the method is often applied to observables whose diagnostic content depends on source-function stratification, contribution functions, opacity asymmetries, or redistribution physics that would be impractical or unavailable in the baseline hydrodynamic calculation. A plausible implication is that the phrase “RADYN+RH” refers less to a single software stack than to a recurrent methodological pattern: dynamic atmosphere first, detailed spectroscopy second.

2. Radiative-transfer physics and numerical formulation

RH and RH 1.5D extend the workflow beyond the native RADYN transfer treatment by emphasizing multi-level NLTE transfer, overlapping radiative transitions, multi-species synthesis, PRD, and, in RH 1.5D, massively parallel column-by-column processing (Pereira et al., 2014). RH 1.5D is based on the RH code of Uitenbroek and specializes in solving the transfer problem independently in each vertical column of a 1D, 2D, or 3D atmosphere. This “1.5D” approximation neglects the effect of inclined rays in a 2D or 3D problem, so only emergent spectra in the vertical direction are computed and horizontal radiative coupling is ignored. The approximation is suitable for continua, weak lines, and many moderate lines, but the main errors occur mostly in the cores of very strong lines in inhomogeneous atmospheres (Pereira et al., 2014).

The numerical architecture is the RH formalism for multi-level NLTE with PRD and blends. RH 1.5D “employs a local approximate $\Psi$ operator and preconditioning of the rate equations, and allows for overlapping radiative transitions” (Pereira et al., 2014). For Mg II flare synthesis, RH is used precisely because RADYN’s CRD-like treatment is inadequate for Mg II h and k, whereas RH includes the hybrid angle-dependent PRD approximation of Leenaarts et al. (2012) (Zhu et al., 2019). The same logic motivates RH use for hydrogen Lyman lines, where PRD is more accurate than the CRD treatment used in RADYN (Brown et al., 2018).

A central diagnostic relation across these studies is the formal solution for emergent intensity,

$$I_{\nu}^{0} = \frac{1}{\mu}\int_z S_\nu e^{-{\tau_\nu}/\mu} \chi_\nu \, dz = \frac{1}{\mu} \int_z C_i \, dz,$$

where $S_\nu$ is the source function, $\tau_\nu$ the monochromatic optical depth, $\chi_\nu$ the opacity, and $C_i$ the contribution function (Costa et al., 2015). RADYN+RH analyses repeatedly rely on this decomposition to identify where line-core and wing photons form, how velocity gradients alter opacity asymmetry, and why line cores may brighten, reverse, or shift.

RH 1.5D also contains explicit convergence-control machinery for difficult PRD problems. Defining

$$\rho(\nu, \mathbf{n}) = \psi(\nu, \mathbf{n}) / \phi(\nu, \mathbf{n}),$$

the code introduces

$\rho' = \rho \epsilon + (1 - \epsilon),$

so that when $_1$0, one has CRD, and as $_1$1, $_1$2 in approximately $_1$3 iterations (Pereira et al., 2014). Together with Ng acceleration, grid optimization, and collisional-radiative switching, this is a major reason RH 1.5D is practical for large chromospheric synthesis jobs.

At the same time, RH often recalculates atomic populations in statistical equilibrium while using non-equilibrium electron densities from RADYN. Multiple papers identify this as a controlled inconsistency. For Mg II and Ca II H/H$_1$4, the approximation is used pragmatically (Zhu et al., 2019, Tamburri et al., 17 Feb 2026). For hydrogen Lyman lines, the discrepancy between RADYN and RH after beam switch-off is interpreted as evidence that non-equilibrium effects can dominate over redistribution differences in the decay phase (Brown et al., 2018).

3. Canonical solar-flare applications

The canonical use case of RADYN+RH simulations is solar flare spectroscopy. A representative example is Na I D$_1$5, modeled with an F11 RADYN atmosphere with beam parameters $_1$6, $_1$7, and $_1$8 keV, heated constantly for 20 s and then allowed to relax for 40 s (Kuridze et al., 2016). In that atmosphere, beam heating increases collisional rates and excites electrons from the ground state $_1$9 to the first excited state $\Psi$0, reducing the $\Psi$1 population ratio by about a factor of 10 at heights around 500–1000 km within the first second. RH then shows that Na I D$\Psi$2 changes from absorption to emission during heating, develops a central reversal only after heating stops, and acquires blue asymmetry from a negative velocity gradient associated with weak upflows of about 2–3 km s$\Psi$3 in the 400–800 km region (Kuridze et al., 2016).

Mg II h, k, and triplet modeling is another defining application. In a 5F11 RADYN atmosphere with a 20 s FWHM heating pulse, RH recalculates Mg II with PRD and shows that non-reversed h and k profiles can appear self-consistently in a compressed upper chromosphere with $\Psi$4, $\Psi$5, and temperature of about $\Psi$6 K, where the source function only slightly deviates from the Planck function (Zhu et al., 2019). The same study showed that RH’s default quadratic Stark broadening for Mg II is about one order of magnitude smaller than STARK-B at $\Psi$7 K, and that even after implementing STARK-B widths an additional empirical factor of 30 was required to match the broadest observed flare-ribbon profiles (Zhu et al., 2019).

RADYN+RH also changed the interpretation of nominally photospheric diagnostics. In Fe I 617.3, 630.1, and 630.2 nm synthesis from a hard-beam RADYN run, RH revealed that the overall line profile can show an apparent red asymmetry by as much as 40 m s$\Psi$8 near maximum beam heating, even though the chromospheric emission contaminating the line is itself blueshifted by as much as 400 m s$\Psi$9 (Monson et al., 2021). The apparent redshift is produced because significant chromospheric emission fills in the blue side of a near-stationary photospheric absorption profile.

Higher-order hydrogen and Ca II diagnostics have been treated similarly. In the first DKIST/ViSP Ca II H and H$$I_{\nu}^{0} = \frac{1}{\mu}\int_z S_\nu e^{-{\tau_\nu}/\mu} \chi_\nu \, dz = \frac{1}{\mu} \int_z C_i \, dz,$$0 flare comparison, state-of-the-art RADYN+RH simulations reproduced certain salient properties such as the width of H$$I_{\nu}^{0} = \frac{1}{\mu}\int_z S_\nu e^{-{\tau_\nu}/\mu} \chi_\nu \, dz = \frac{1}{\mu} \int_z C_i \, dz,$$1, but severely underestimated the width of Ca II H in the red wing and failed to reproduce the exact relative intensity of Ca II H to H$$I_{\nu}^{0} = \frac{1}{\mu}\int_z S_\nu e^{-{\tau_\nu}/\mu} \chi_\nu \, dz = \frac{1}{\mu} \int_z C_i \, dz,$$2 (Tamburri et al., 17 Feb 2026). In hydrogen Lyman-line work, RADYN and RH together showed that upflows in the simulated atmospheres lead to blueshifts in centrally reversed line cores, but after convolution to EVE resolution these profiles can be interpreted as redshifted emission, so apparent downflows do not necessarily correspond to actual downflows (Brown et al., 2018).

4. Cross-diagnostic tests and observational tensions

A major strength of RADYN+RH simulations is that they permit simultaneous testing of dynamical and spectroscopic consistency across multiple diagnostics. A major weakness is that those diagnostics often impose competing constraints.

Ellerman-bomb modeling provides a clear example. In 1D RADYN atmospheres post-processed with MULTI for H$$I_{\nu}^{0} = \frac{1}{\mu}\int_z S_\nu e^{-{\tau_\nu}/\mu} \chi_\nu \, dz = \frac{1}{\mu} \int_z C_i \, dz,$$3 and Ca II 8542 Å and with RH in PRD for Mg II h&k, heating around 500–700 km at about $$I_{\nu}^{0} = \frac{1}{\mu}\int_z S_\nu e^{-{\tau_\nu}/\mu} \chi_\nu \, dz = \frac{1}{\mu} \int_z C_i \, dz,$$4 produced the best classical EB-like H$$I_{\nu}^{0} = \frac{1}{\mu}\int_z S_\nu e^{-{\tau_\nu}/\mu} \chi_\nu \, dz = \frac{1}{\mu} \int_z C_i \, dz,$$5 and Ca II 8542 Å wing brightenings, while Mg II h&k favored heating around 700–900 km. But the Mg II-favorable models also made Ca II 8542 Å show strong core emission, so the same simple 1D heating scenario could not reproduce all diagnostics simultaneously (Reid et al., 2017).

An analogous tension appears in impulsive-phase chromospheric evaporation studies that combined F-CHROMA RADYN atmospheres with RH 1.5D for C II 1334.5 Å and CHIANTI for Fe XXI 1354.1 Å (Sadykov et al., 2018). The observations gave a gentle-to-explosive transition energy flux of $$I_{\nu}^{0} = \frac{1}{\mu}\int_z S_\nu e^{-{\tau_\nu}/\mu} \chi_\nu \, dz = \frac{1}{\mu} \int_z C_i \, dz,$$6, while the models gave $$I_{\nu}^{0} = \frac{1}{\mu}\int_z S_\nu e^{-{\tau_\nu}/\mu} \chi_\nu \, dz = \frac{1}{\mu} \int_z C_i \, dz,$$7, so the thresholds were comparable. Yet significant discrepancies remained: Fe XXI Doppler shifts were substantially stronger in the models than in the data, and the C II mean blueshifts predicted by the models were absent in the observations (Sadykov et al., 2018).

The DKIST Ca II H/H$$I_{\nu}^{0} = \frac{1}{\mu}\int_z S_\nu e^{-{\tau_\nu}/\mu} \chi_\nu \, dz = \frac{1}{\mu} \int_z C_i \, dz,$$8 study sharpened a related point: the width of H$$I_{\nu}^{0} = \frac{1}{\mu}\int_z S_\nu e^{-{\tau_\nu}/\mu} \chi_\nu \, dz = \frac{1}{\mu} \int_z C_i \, dz,$$9 is not solely related to condensation properties. The modeled condensation electron densities at the height of maximum H$S_\nu$0 contribution span from $S_\nu$1 to $S_\nu$2, more than an order of magnitude, yet the resulting H$S_\nu$3 widths are broadly similar (Tamburri et al., 17 Feb 2026). This was interpreted to mean that deeper flare layers shape both H$S_\nu$4 width and intensity, not just the peak condensation density.

These tensions have methodological significance. They show that RADYN+RH simulations are most informative when they are used as multi-diagnostic consistency checks rather than as single-line fitting engines. This suggests that disagreements among H$S_\nu$5, Ca II, Mg II, Fe I, and hydrogen continua are not merely nuisances; they are direct evidence of missing structure, missing broadening physics, or overly restrictive heating assumptions.

5. Multidimensional and stellar extensions

Although RADYN itself is 1D, the RH side of the workflow has been extended naturally to multidimensional atmospheres. RH 1.5D can read 1D, 2D, or 3D atmospheres and solve the transfer independently in each column, which makes it suitable for post-processing large simulation snapshots when the 1.5D approximation is acceptable (Pereira et al., 2014). This makes possible a multidimensional analogue of the familiar RADYN→RH strategy.

A direct example is the 3D StellarBox + RH study of particle-beam heating of the solar atmosphere. There the dynamical atmosphere is evolved with the 3D radiative-MHD code StellarBox, and RH 1.5D is then used to synthesize Fe I 6173 Å Stokes profiles from extracted vertical columns. The beam flux density is fixed at $S_\nu$6, the spectral index is $S_\nu$7, and the low-energy cutoffs are $S_\nu$8, $S_\nu$9, $\tau_\nu$0, and $\tau_\nu$1 keV (Granovsky et al., 30 Dec 2025). The simulations produce strong upper-chromospheric heating, multiple shock fronts, and continuum enhancements up to a factor of 2.5 relative to pre-flare levels, yet still do not produce full Fe I 6173 Å emission. The comparison with 1D RADYN was interpreted to show that realistic 3D fine structure changes shock complexity, continuum structuring, and cooling, while not overturning the standard result that these electron beams deposit their energy mostly above the low chromosphere (Granovsky et al., 30 Dec 2025).

RADYN+RH methods have also been exported to stellar flare work. In high-cadence ULTRACAM observations of dMe flares, RADYN supplied F13 atmospheres and RH post-processing was used to synthesize spectra with higher hydrogen levels and Landau-Zener opacity modifications (Kowalski et al., 2016). The favored model used a beam flux of $\tau_\nu$2 and, in the improved case, a single power law with $\tau_\nu$3 evaluated at $\tau_\nu$4 s. This RH-processed atmosphere reproduced $\tau_\nu$5 and $\tau_\nu$6, close to the observed impulsive flare values, and supported an interpretation in which the hot optical color temperature and small Balmer jump arise from optically thick hydrogen recombination in a dense chromospheric condensation plus deeper heated stationary layers rather than from a literal blackbody (Kowalski et al., 2016).

6. Model bases, assumptions, and unresolved issues

RADYN+RH simulations depend strongly on the atmospheric archive and beam prescription that feed the RH stage. The F-CHROMA grid provides a standardized public set of 96 1D RADYN flare models spanning $\tau_\nu$7–8, $\tau_\nu$8 to $\tau_\nu$9, and $\chi_\nu$0–25 keV, with a triangular 20 s beam and 30 s relaxation (Carlsson et al., 2023). Those models do not use RH internally, but they are explicitly presented as atmospheres that can be post-processed with RH, Lightweaver, or CHIANTI. They also come with explicit caveats: all bound-bound transitions are treated with CRD; PRD is not treated properly for hydrogen Lyman lines or Ca II H&K; the models are 1D; and additional heating from return currents is not included (Carlsson et al., 2023).

Several recurrent assumptions delimit the scope of the method. First, RH post-processing usually has no feedback on the RADYN atmosphere. Second, RH often assumes statistical equilibrium even when RADYN evolved the atmosphere with time-dependent non-equilibrium populations. Third, RH 1.5D neglects horizontal radiative transfer, so only vertical emergent spectra are computed in multidimensional post-processing (Pereira et al., 2014). Fourth, broadening prescriptions remain a live issue: Mg II required STARK-B plus an additional factor of 30 in one flare-ribbon study (Zhu et al., 2019), while Ca II H red-wing broadening remains severely underpredicted in DKIST flare spectra (Tamburri et al., 17 Feb 2026).

Upstream flare physics also matters. Work coupling the Stanford unified Fokker–Planck code FLARE to RADYN showed that the shape of the non-thermal electron distribution, not only the total beam energy flux, substantially changes where energy is deposited, how strongly the atmosphere evaporates, and what line profiles emerge (Costa et al., 2015). That study was not itself an RH paper, but it implies that downstream RH synthesis is sensitive to whether RADYN was driven by an ad hoc power law or by a quasi-thermal-plus-power-law spectrum. A plausible implication is that some apparent disagreements attributed to transfer physics may in fact originate in the beam distribution feeding the hydrodynamic solution.

For these reasons, RADYN+RH simulations are best understood as a high-fidelity but not fully closed forward-modeling framework. They excel when the target problem requires dynamic atmospheres together with detailed NLTE/PRD line formation, contribution-function analysis, or synthetic observables beyond the native RADYN setup. They remain limited by 1D or 1.5D geometry, statistical-equilibrium assumptions in the post-processing stage, incomplete broadening physics, and the dependence of spectral results on the upstream heating prescription. The continuing literature treats those limitations not as minor technical details, but as the central agenda for the next generation of flare and transient-atmosphere modeling (Tamburri et al., 17 Feb 2026, Granovsky et al., 30 Dec 2025, Carlsson et al., 2023).