Multi-resonant-line radiative transfer: Lyman-Alpha fine structure and deuterium coupling

Abstract

Resonance lines encode rich information about astrophysical sources and their environments, yet fully analytic treatments of multi-line radiative transfer remain almost entirely unexplored. We present exact, closed-form solutions for steady-state resonant-line radiativeP transfer in “V-shaped” atomic networks, where a single ground state couples to multiple transitions. Starting from the full angle-dependent transfer equation, we generalise absorption and emission coefficients to an arbitrary number of lines, derive a modified Fokker–Planck expansion of the frequency-redistribution integral, and use a judicious change of variables to reduce the problem to a Helmholtz equation with point-like sources in frequency space. This transformation admits analytic solutions for arbitrary sets of lines with fixed frequency offsets and strengths in both slab and spherical geometries. We implement V-shaped line networks in the colt Monte Carlo radiative transfer code and find excellent agreement with the analytic predictions across a wide range of line separations, optical depths, and damping parameters, establishing our solutions as stringent validation benchmarks. For concrete applications related to the Lyman-alpha (Ly$\alpha$) transition of neutral hydrogen, we examine how fine-structure splitting and deuterium injection modify the emergent spectra, internal radiation field, and radiative force multiplier. We show that these effects leave previous conclusions about Ly$\alpha$ feedback in the early universe essentially unchanged. Even when direct observational diagnostics are subtle, our framework provides novel analytic and numerical insights into coupled resonance-line transport and facilitates progress in general modelling of multi-line radiative transfer in diverse astrophysical settings.

1. Introduction

Radiative transfer (RT) plays a central role in astrophysics by describing how light propagates and interacts with astrophysical media. The RT equation (RTE) provides a powerful framework for interpreting observations and modelling multi-scale radiation-induced phenomena, but solving it exactly is often intractable due to the complexity of scattering processes, especially for general three-dimensional geometries (Haworth et al. 2018). As a result, various approximations are often made to obtain analytical solutions. Examples include assuming a steady state, integrating over all photon frequencies, or simplifying the scattering physics involved (Allakhverdian & Naumov 2023). These approximations trade some physical fidelity for mathematical simplifications, and careful choices must be made to capture the dominant effects.

One application with particularly rich physics is the RT of the hydrogen Lyman-alpha (Ly $\alpha$ ) line. Ly $\alpha$ is intrinsically bright due to the pervasive nature of photo-ionisation producing recombination events, each of which has probability $P(\text{Ly}\alpha) \approx 0.68$ of yielding a Ly $\alpha$ photon with a high cross-section for scattering with neutral hydrogen (e.g. Dijkstra et al. 2019). Understanding any mechanisms impacting Ly $\alpha$ RT in the local and high-redshift universe is therefore critical for interpreting observations of star-forming galaxies and the intervening intergalactic medium (IGM). However, Ly $\alpha$ photons undergo multiple resonant scatterings, making their propagation highly diffusive in both physical and frequency space. This process shapes the emergent Ly $\alpha$ spectra and surface brightness profiles, complicating any direct connection between observed signals and their source properties. Analytical models, even if idealised, are invaluable for clarifying these effects and guiding interpretations of Ly $\alpha$ observations and simulations.

Over the years, an increasing number of analytic solutions for Ly $\alpha$ have been extensively developed under various simplifying assumptions. An important milestone was the adoption of the Fokker–Planck approximation for partial frequency redistribution, allowing closed-form solutions for infinite plane-parallel slab geometry with central point and uniform sources by Harrington (1973) and Neufeld (1990), building on work done by Unno (1955). Then, configurations with spherical and cubic symmetry were derived by Dijkstra et al. (2006) and Tasitsiomi (2006), respectively, extending the theory to different geometries. These classical solutions assumed static, homogeneous media and yielded characteristic double-peaked Ly $\alpha$ line profiles in the emergent spectrum. While the basic behaviour is well understood, significant progress is still being made with explicit formulas for the internal radiation field within the medium by Lao & Smith (2020) and Seon & Kim (2020), whereas McClellan, Davis, & Arras (2022) incorporated time-dependent, impulsive emission sources into the analytic framework. Variations in density and source distributions were explored by Lao & Smith (2020) for power-law density and emissivity gradients, while velocity gradients have been included by Nebrin et al. (2025) and Smith et al. (2025) to account for bulk flows, e.g. outflows. Most comprehensively, Nebrin et al. (2025) presented a generalised series solution that combines several additional physical effects, such as macroscopic velocity fields, atomic recoil, Ly $\alpha$ photon destruction via $2\text{p} \to 2\text{s}$ transitions, small-scale density fluctuations induced by turbulence, and continuum absorption/scattering, all within a single model. Each of these advancements has improved the realism of analytic RT models. However, despite this extensive progress, all of the above studies are limited to a single spectral line at a time. In other words, existing analytic solutions describe strong line RT in isolation, without accounting for the possibility of multiple resonant transitions interacting or overlapping in a common medium.

Modelling the transfer of multiple spectral lines simultaneously is a natural next step with broad implications. In the context of Ly $\alpha$ RT, considering multiple lines can capture processes such as fine-structure splitting and the Wouthuysen–Field effect (the coupling of Ly $\alpha$ photons with the 21-cm hyperfine transition). More generally, multi-line RT is essential in environments where different lines overlap or feed into each other. Notable examples include the rich Lyman-series cascades in recombining gas, the Balmer lines in dense nebulae, and line-driven outflows around active galactic nuclei (AGNs), accretion disk winds, and winds around massive stars and supernovae. In such systems, radiation in one transition can alter the population of atoms or ions in excited states, thereby affecting the optical depth and emission in another transition. Numerical simulations have been used to model multi-line RT (Lucy 2002; Matthews et al. 2025; Druett & Zharkova 2018; Chang et al. 2025). To make this feasible, Monte Carlo codes typically employ the ‘macro-atom’ formalism or Sobolev approximation to handle the increased complexity, but often by sacrificing detailed treatment of the frequency redistribution or assuming instant re-emission in a different line. In some contexts, iterative schemes have been used to solve coupled RTEs for all optically thick lines of hydrogen simultaneously. These approaches have provided valuable insights, but often discard the spectral memory of previous scatterings. To date, we are not aware of any closed-form analytic solution that treats more than one resonant line at a time with detailed partial frequency redistribution physics.

In this paper, we present the first analytic solution of the RTE that self-consistently accounts for multiple resonant lines simultaneously. We focus on the simplest non-trivial case of a ‘V-shaped’ two-level atom network, in which a single ground state is coupled to multiple independent excited states and their corresponding atomic line transitions. This scenario represents the least-complex case of a broader class of multi-line networks, where one would consider higher-order transitions and radiative cascades, greatly increasing the scattering physics complexity. Still, the V-shaped system already captures phenomena like Ly $\alpha$ fine-structure doublet splitting and hydrogen–deuterium line blending. Our analytic solution is obtained under the same assumptions as used in single-line Ly $\alpha$ solutions. We also apply the diffusion approximation for frequency redistribution, which treats each scattering as an effectively instantaneous, small frequency change event. Within this coupled framework, we derive a closed-form expression for the stead-state spectral intensity as a function of position and frequency, for both the internal radiation field and the emergent flux escaping the medium. An important development enabling this solution is the generalisation of the Fokker–Planck formulation of the redistribution function for multiple lines, a powerful result that will be derived more generally in a future paper.

To validate our analytic results and illustrate their utility, we have also developed a complementary numerical approach. We modified the Monte Carlo based Cosmic Ly $\alpha$ Transfer code (COLT; Smith et al. 2015) to handle V-shaped multi-line networks with accurate frequency redistribution. In essence, our code allows photons to scatter into any of the available upper-state transitions, with the choice of the line determined probabilistically by the local relative line optical depths. Using this tool for explicit multi-line resonant scattering physics, we perform simulations of Ly $\alpha$ fine-structure transfer and deuterium–hydrogen mixed media, and we show that the simulation results agree quantitatively with the predictions of our analytic solution. This cross-validation gives us confidence in the correctness of the theory, highlighting the improved physical accuracy as both analytic and numerical methods involve non-trivial generalisations.

The rest of this paper is organised as follows. In Section 2, we introduce the governing RTE and generalisations of the many variables in the RTE, such as the definition of a multi-line dimensionless frequency coordinate and a combined absorption profile. In Section 3, we present the full analytical solution for both internal and emergent spectra in the case of the V-shaped network, with the details of the derivation and other details provided in the appendices. In Section 4, we outline the modifications made to colt to simulate V-shaped scattering, and we demonstrate its accuracy with the analytical solutions developed in Section 3 for various test cases using arbitrary atomic parameters. In Section 5, we apply our new tools to two relevant astrophysical scenarios of Ly $\alpha$ fine-structure and deuterium injection. These examples show how multi-line effects manifest in realistic conditions and confirm that our models reproduce known results from the literature in the appropriate limits. Finally, in Section 6, we summarise our findings, discuss the astrophysical implications of multi-line RT, and outline potential directions for future work, including the incorporation of additional physics and the extension to more complex multi-level atomic networks. Appendix A contains the derivation of the analytical solutions of this text, and Appendix B contains a discussion on the required change of frequency variables. Finally, Appendix C contains a table of the parameters used for each figure for reproducibility purposes and to explain line asymmetries.

2. Governing equations

The general form of the steady-state RTE in static media is:

(1) \begin{align} \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\nabla} I_\nu = j_\nu - k_\nu I_\nu + \iint I_{\nu'} (kR)_{\nu', \boldsymbol{n}' \rightarrow \nu, \boldsymbol{n}}\,\text{d}\Omega' \text{d}\nu' ,\end{align}

where $I_\nu(\boldsymbol{r}, \boldsymbol{n})$ is the specific intensity (for frequency $\nu$ at position $\boldsymbol{r}$ travelling in direction $\boldsymbol{n}$ ), $j_\nu(\boldsymbol{r})$ and $k_\nu(\boldsymbol{r})$ denote the emission and absorption coefficients, and the last term accounts for frequency redistribution due to partially coherent scattering (Dijkstra et al. 2019). The term (kR) relates to the redistribution function $R_{\nu', \boldsymbol{n}' \rightarrow \nu, \boldsymbol{n}}$ , which is the differential probability per unit initial photon frequency $\nu'$ and per unit initial directional solid angle $\Omega'$ that the scattering of such a photon travelling in direction $\boldsymbol{n}'$ would place the scattered photon at frequency $\nu$ and directional unit vector $\boldsymbol{n}$ (Hummer 1962). The notational convention of (kR) was chosen for convenience with multi-line extensions and has the same units as $k_\nu$ . Physically, R encodes the quantum mechanical scattering processes (e.g. Doppler shifts, natural broadening, recoil, etc.), and couples the RTE across different frequencies.

The frequency dependence of the absorption coefficient is given by the Voigt profile $\phi_{\text{Voigt}}$ centred on the line resonance. For convenience we define the Hjerting–Voigt function $H(x, a) = \sqrt{\pi} \Delta \nu_{D} \phi_{\text{Voigt}}(\nu)$ as the dimensionless convolution of Lorentzian and Maxwellian distributions,

(2) \begin{align} H(x, a) = \frac{a}{\pi} \int_{-\infty}^\infty \frac{e^{-y^2}\text{d}y}{a^2+(y-x)^2} \approx \begin{cases} e^{-x^2} & \quad \text{`core'} \\[2pt] {\displaystyle \frac{a}{\sqrt{\pi} x^2} } & \quad \text{`wing'} \end{cases}\end{align}

where the ‘damping parameter’ is $a \equiv \Delta \nu_{L} /2 \Delta \nu_{D}$ , $\Delta \nu_{L}$ is the natural line width of the line, $\Delta \nu_{D} \equiv \nu_0\,v_{\text{th}} / c$ is the Doppler width, and $v_{\text{th}} = (2\,k_{{B}}T/m)^{1/2}$ is the thermal velocity at temperature T for carrier mass m, such that $x = (\nu-\nu_0)/\Delta\nu_{D}$ measures the frequency offset from line centre in Doppler units.

For a multi-line system, we generalise the notation to arrive at a compact modified RTE. With N lines, we have multiple distinct resonance frequencies $\nu_{0, i}$ for $i \in \{1,\, ...,\, N\}$ . We denote the dimensionless frequency variables for each individual line as $x_i$ and redefine x as the unified dimensionless frequency variable as

(3) \begin{align} x_i \equiv \frac{\nu - \nu_{0,i}}{\Delta\nu_{{D},i}} = \mu_i(x-\Delta_i) \quad \text{and} \quad x \equiv \frac{\nu - \nu_0}{\Delta\nu_{D}} = \frac{1}{N} \sum_{i = 1}^{N}x_i , \end{align}

respectively, where we defined $\mu_i = \Delta\nu_{D} / \Delta\nu_{{D},i}$ . Here $\nu_0$ and $\Delta\nu_{D}$ can be arbitrary reference values but for concreteness we may adopt the mean line frequency $\nu_0 \equiv \sum_{i=1}^N\nu_{0,i}/N$ and Doppler width $\Delta\nu_{D} = \sum_{i = 1}^N \Delta\nu_{{D}, i} / N$ . The dimensionless offset of each line in this common x coordinate is $\Delta_i \equiv (\nu_{0, i} - \nu_0) / \Delta\nu_{D}$ . In practice, $\mu_i$ departs from unity when lines have relatively large separations or different carrier masses as $\Delta\nu_{{D}, i} \propto \nu_{0,i} / \sqrt{m_i}$ . Throughout this work we will use x as the master frequency variable and treat each line’s profile in shifted coordinates $x_i$ .

The total absorption coefficient is the sum of contributions from all lines (Ahn & Lee 2015; Seon & Kim 2020):

(4) \begin{align} k_x &\equiv \sum_{i=1}^N k_{x,i} = \sum_{i=1}^N n_i\,\sigma_{0, i}\,H(a_i, x_i) \nonumber \\[2pt] &= \sum_{i = 1}^N n_i\,f_{lu,i}\,\frac{\pi\,e^2}{m_e\,c\,\Delta\nu_{{D},i}}\,H(a_i, x_i) \equiv k(\boldsymbol{r})\bar{H}(x) ,\end{align}

where $n_i$ is the number density of atoms in the lower level of a given transition and $\sigma_{0, i}$ is the line centre cross-section defined in terms of the oscillator strength $f_{lu,i}$ . In the last equality, we have factored out the position dependence into $k(\boldsymbol{r}) = \sum_{i = 1}^Nn_i\sigma_{0, i}$ , the total line-centre absorption coefficient of all lines, and defined a multi-lined Voigt profile that encodes the relative strength of each line within dimensionless weights to represent each contribution:

(5) \begin{align} \bar{H}(x) \equiv \sum_{i=1}^N \omega_i H(a_i,x_i) \quad \text{where} \quad \omega_i \equiv \frac{n_i\,f_{lu, i}} {\Delta\nu_{{D}, i}} \Big/ \sum_{i=1}^N \frac{n_i\,f_{lu, i}} {\Delta\nu_{{D}, i}} \, .\end{align}

By construction, the weights are normalised to unity, $\sum_{i=1}^N $ $\omega_i = 1$ , and reduce to the simple abundance ratio of level populations in the case of lines with similar oscillator strengths and thermal widths. We later simplify each of these coefficients further for the specific applications explored in this paper.

We must extend the scattering term in Eq. (1) to account for multiple lines. For V-shaped networks, the redistribution function can be written as a sum of contributions from each line scattering back into the same line. Adopting the common x coordinate for the integrand but with line-specific parameters gives:

(6) \begin{align} &\iint I_{x'} (kR)_{x', \boldsymbol{n}' \rightarrow x, \boldsymbol{n}}\,\text{d}x' \text{d}\Omega' \longrightarrow \notag \\[2pt] & \, k(\boldsymbol{r}) \sum_{i = 1}^N \omega_i \iint I_{x'} H(a_i, x'_i) R_{x' - \Delta_i, \boldsymbol{n}' \rightarrow x - \Delta_i, \boldsymbol{n}}\,\text{d}x'\text{d}\Omega' \, .\end{align}

We prefer this form over an equivalent sum of shifted frequency integrals to unify Fokker–Planck terms.

Figure 1. Multi-line Voigt profile $\bar{H}(x)$ from Eq. (5) for a two-line system with asymmetric line strengths. The stronger line dominates except in the immediate vicinity of the weaker line. The normalisation also affects the weak-line heights.

3. Analytical solutions

Using the formulations above, we now derive analytical solutions to the RTE for multi-line V-shaped networks. The full derivation can be found in Appendix A; here we summarise the results and provide the final expressions for the specific intensity. We focus primarily on the simplest case of two resonant ground state transitions ( $N=2$ ), generalising afterwards. For clarity, we follow the notation of Lao & Smith (2020), and write the general multi-line angle-averaged intensity J as governed by

(7) \begin{align} \frac{1}{k(\boldsymbol{r})} \boldsymbol{\nabla}\!\boldsymbol{\cdot}\!\left(\frac{\boldsymbol{\nabla} J}{k(\boldsymbol{r})} \right)\!+\!\frac{3}{2} \bar{H}(x) \frac{\partial}{\partial x}\!\left(\!\hat{H}(x) \frac{\partial J}{\partial x} \right) = -\frac{3\mathcal{L}}{4\pi}\!\frac{\eta(\boldsymbol{r})}{k(\boldsymbol{r})}\!\frac{\bar{H}\tilde{H}}{\sqrt{\pi}} ,\end{align}

where we defined further combinations of individual Voigt profiles

(8) \begin{align} \tilde{H}(x) = \sum_{i=1}^{N}\omega_i\mu_iH(x_i) \quad \text{and} \quad \hat{H}(x) = \sum_{i=1}^{N}\frac{\omega_i}{\mu_i^2}H(x_i) ,\end{align}

which are used to make the equations analytically tractable, and are justified by the following frequency change of variables:

(9) \begin{align} \text{d}\tilde{x} = \sqrt{\frac{2}{3}} \frac{\text{d}x}{\tau_0(\bar{H}\hat{H})^{1/2}} \quad \implies \quad \tilde{x} = \int_0^{x}\sqrt{\frac{2}{3}}\frac{\text{d}x'}{\tau_0(\bar{H}\hat{H})^{1/2}} \, .\end{align}

This yields a solution for any general N-line V-shaped network:

(10) \begin{align} J(\tilde{\boldsymbol{r}}, \tilde{x}) = \sum_{i = 1}^N \omega_i\mu_i J_i(\tilde{\boldsymbol{r}}, \tilde{x}-\tilde{\Delta}_i) \, ,\quad \text{where} \\[-10pt] \nonumber \end{align}

(11) \begin{align} J_i(\tilde{\boldsymbol{r}},\tilde{x}-\tilde{\Delta}_i) = \frac{\mathcal{L} \sqrt{6}}{8 \pi} \sum_{n=1}^\infty \frac{Q_{n}}{\lambda_{n}} \vartheta_n(\tilde{\boldsymbol{r}}) e^{-\lambda_{n}|\tilde{x}-\tilde{\Delta}_i|} , \\[8pt] \nonumber \end{align}

where $\mathcal{L}$ is the total luminosity of the source embedded in the medium, $\vartheta_n(\boldsymbol{r})$ and $\lambda_n$ are the spatial eigenfunctions and eigenvalues obtained from the diffusion equation depending on the geometry and boundary conditions, and $Q_n$ is a geometrical coupling coefficient defined in Appendix A. The tildes on $\tilde{\boldsymbol{r}}$ and $\tilde{x}$ indicate that these coordinates have been transformed into dimensionless rescaled coordinates, i.e. the position normalised to the optical depth and a line-integrated frequency defined and calculated in Appendix B that greatly simplifies the differential equation. When compared with the single line solutions of Lao & Smith (2020), the multi-line solution is simply a linear combination of single-line solutions centred on each line, i.e. a superposition weighted by the relative line strength $\omega_i$ and scaled by $\mu_i$ to account for differences in thermal width. This is a gratifyingly simple result, as the complicated coupling between lines can be encapsulated by known results, provided that it is possible to work out the common transformed frequency variable $\tilde{x}$ .

To make these results concrete, we now specialise to commonly assumed geometries and source configurations. We consider (i) an infinite homogeneous slab of height Z with a symmetric source at the mid-plane, and (ii) a uniform sphere of radius R with a central source, both characterised by centre-to-edge optical depths of $\tau_0$ . These are the same setups originally considered by Harrington (1973), Neufeld (1990), and Dijkstra et al. (2006) for single-line transfer, and therefore provide benchmark cases to understand how the presence of a second line alters classic solutions.

For slab geometry, we adopt the dimensionless spatial coordinate $\tilde{z} = z / Z$ and place a delta-function source at $z=0$ . The solution for the multi-line internal angle-averaged intensity is

(12) \begin{align} J(\tilde{z},\tilde{x}) = \frac{\mathcal{L}\sqrt{6}}{8\pi^2} \sum_{i = 1}^N \omega_i\mu_i\,\text{tanh}^{-1}\!\left\{\text{cos}\!\left(\frac{\pi\tilde{z}}{2}\right)\text{sech}\!\left[\frac{\pi(\tilde{x}-\tilde{\Delta}_i)}{2}\right]\right\} ,\end{align}

and the emergent (surface) spectrum from the slab is obtained by evaluating the boundary condition at $\tilde{z} = 1$ , which yields

(13) \begin{align} J(\tilde{x}) = \frac{\mathcal{L}\sqrt{6}}{16\pi f\tau_0\bar{H}(\tilde{x})} \sum_{i = 1}^N \omega_i\mu_i\,\text{sech}\!\left[\frac{\pi(\tilde{x}-\tilde{\Delta}_i)}{2}\right] ,\end{align}

where f is an order unity constant determined by the choice of boundary conditions. These results reduce to the standard single-line solutions when one line dominates, i.e. $\omega_1 \approx 1$ and $\tilde{\Delta}_1 \approx 0$ . The general behaviour of the $N=2$ slab solution is illustrated in the top panel of Figure 2, showing that as the optical depth increases, the intensity profile broadens and the two lines begin to blend together into a single, wider profile once the line separation is within the reach of the resonant broadening, i.e. $\Delta \lesssim (a \tau_0)^{1/3}$ . For high enough optical depth, photons collectively scatter far enough in frequency to bridge the gap between the two resonances, making the RT resemble that of a single effective line with combined optical depth.

Figure 2. A two-line V-shaped network with various optical depths and line separations. The asymmetry in chosen line parameters ( $\omega_1 = 1/3$ and $\omega_2 = 2/3$ ) causes a difference between the peaks, with the second line being much stronger than the first. Upper Panel: Internal spectrum for a plane-parallel slab at $\tilde{r} = 0.5$ , changing the optical depth over $\tau_0 \in \{10^4, 10^5, 10^6, 10^7\}$ at a temperature of $T = 10^4$ K, based on Eq. (12). Higher optical depth enables more diffusion across frequency space, since photons must typically reach a critical frequency of order $x_{\text{esc}} \sim (a \tau_0)^{1/3}$ before escape, increasing the likelihood of coupling between nearby lines. Lower Panel: Internal spectrum for a spherical cloud geometry with different line separations from the centre as parametrised by $\Delta \in \pm \{0.5, 1.25, 2, 2.75\}\,(a\tau_0)^{1/3}$ at $\tilde{r} = 0.5$ , based on Eq. (14). As the spacing increases the solution separates into two distinct profiles, decoupling from each other. While there is still significant overlap the relative heights are similar due to mixing between lines.

For spherical geometry, we adopt $\tilde{r} = r / R$ and place a delta-function source at $r=0$ . The multi-line internal solution is

(14) \begin{align} J(\tilde{r},\tilde{x}) = \frac{\mathcal{L} \sqrt{6}}{32 \pi^2 R^2 \tilde{r}} \sum_{i = 1}^N \frac{\omega_i\mu_i\,\sin(\pi \tilde{r})}{\cosh\![\pi(\tilde{x}-\tilde{\Delta}_i)] - \cos\!(\pi \tilde{r})} ,\end{align}

and the emergent intensity is

(15) \begin{align} J(\tilde{x}) = \frac{\mathcal{L}\sqrt{6}}{64 \pi R^2 f \tau_0 \bar{H}(\tilde{x})} \sum_{i = 1}^N \omega_i\mu_i\,\text{sech}^2\!\left[\frac{\pi(\tilde{x}-\tilde{\Delta}_i)}{2}\right] \, .\end{align}

In the bottom panel of Figure 2, we show the internal spectrum at an arbitrary radius $\tilde{r} = 0.5$ for a two-line system. Again, when the separation between the lines is small compared to the diffusion width of each line, the two peaks overlap and blend into each other. Conversely, when the lines are well separated (relative to their widths), even the internal spectrum shows two distinct diffusion peaks, quantifiable by the separability condition $|\Delta_2 - \Delta_1| \gg (a \tau_0)^{1/3}$ . More precisely, the full-width at half-maximum (FWHM) of a single-line profile can be found using equation (90) in Lao & Smith (2020) in $\tilde{x}$ space, converting to x space using Eq. (A9) and the approximation $\bar{H}(x) \approx \omega_i H_i(x_i, a_i)$ valid around the neighbourhood $x \approx \Delta_i$ , which gives $\text{FWHM}_i \approx 1.67\,(\omega_i a_i \tau_0)^{1/3}$ . Overlap occurs when $|\Delta_i - \Delta_j| \lesssim \text{FWHM}_i + \text{FWHM}_j$ . For example, in Figure 2, we have $(\text{FWHM}_1 + \text{FWHM}_2) / (a\tau_0)^{1/3} = 2.57$ , indicating a significant overlap for $\Delta \leq 1.3\,(a\tau_0)^{1/3}$ . A more subtle variable that also impacts line separation is temperature, as we can see from $\Delta_i \propto T^{-1/2}$ and $(a_i\tau_0)^{1/3} \propto T^{-1/3}$ , and so at higher temperatures the distributions merge together more easily.

We stress that while Eq. (10) suggests a simple superposition of single-line solutions, one must be careful in interpreting the individual terms. The transformation $\tilde{x} = \tilde{x}(x)$ depends on the combined profile $\bar{H}(x)$ (see Appendix B), which means that in regions of frequency space dominated by one line, the coordinate mapping for the other line can become highly non-linear. An illustration of this phenomenon is given in Figure 3, where we present the contributions of each of the lines from Eq. (14). Individual terms do not match the standard single-line solutions derived in Lao & Smith (2020) centred on the same frequency. In fact, the $J_2$ term deviates strongly because $\text{d}\tilde{x}/\text{d}x \rightarrow 0$ around $x \approx \Delta_1$ . Therefore, these individual terms are not to be used in isolation, only the sum J is physically meaningful across all frequencies. To better understand the behaviour of $\tilde{x}$ , see Appendix B for more details. Finally, we note that many other quantities typically derived also have similar superposition relations, i.e. energy density.

Figure 3. Internal spectrum $(\tilde{r} = 0.5)$ contributions of individual lines, i.e. $J_1$ and $J_2$ from Eq. (14), in comparison to standard single-line solutions. We emphasise that the terms exhibit qualitatively different behaviour as only the sum J is physically meaningful across the full spectrum, rather than $J_1$ and $J_2$ in isolation.

4. Numerical simulations

To check the validity of our analytic solutions, we implemented a new flexible multi-line resonant scattering Monte Carlo radiative transfer (MCRT) module within the Cosmic Ly $\alpha$ Transfer code (COLT; Smith et al. 2015, 2019, 2022; McClymont, Smith, & Tacchella 2025, for public code access and documentation see colt.readthedocs.io ). The colt code is well-tested for Ly $\alpha$ and we extend this functionality to handle N-line V-shaped atomic networks as in Section 2. The modifications are similar to those described by Seon & Kim (2020) for including the Wouthuysen–Field effect, i.e. including fine-structure for Ly $\alpha$ –21 cm coupling in MCRT frameworks. In essence, the algorithm proceeds as follows: to trigger scattering events we draw a random optical depth $\tau_{\text{scat}} = -\ln\xi$ , where $\xi$ is a random univariate, from the exponential distribution $e^{-\tau}$ and calculate the scattering distance by equating this to the total traversed optical depth from a minimal set of nearby lines. We then draw another random number to determine which line the photon actually interacts with. This is done by comparing the relative contributions of each line to the total optical depth, as the probability to scatter into line i is $P_i = \tau_i / \sum \tau_i$ , where the $\tau_i$ are partial optical depths for each line. After selecting the line, we treat the event as a standard resonant scattering in that transition, drawing the outgoing direction and frequency based on the atomic properties of the line (either artificial or real data, e.g. Sansonetti et al. 2004) and the corresponding phase and redistribution functions. This sampling captures the true microphysics, where photons interact with one atom at a time. In this paper, we focus on V-shaped networks, assuming all atoms are in the ground state, which is valid where the radiation field is not intense enough to significantly populate excited states. Future work will explore non-LTE resonant pumping effects by iterating between transport and level population statistics. We emphasise that our MCRT implementation fully retains frequency-dependent scattering physics without Sobolev-like or escape probability approximations; thus, it serves as a robust numerical check on our analytic diffusion solutions. We refer to a specific set of levels, allowed transitions, and their corresponding atomic data as a ‘network’. We feature three networks in this paper: (i) an arbitrary two-line V-shaped network with variable line separation, (ii) the hydrogen Ly $\alpha$ fine-structure doublet, and (iii) a Deuterium injected Ly $\alpha$ network.

For simplicity, in all simulations we adopt the same phase functions as most Ly $\alpha$ codes (Dijkstra et al. 2006), i.e. dipole scattering in the wings with modifications in the core. While emergent spectra are in fact nearly independent of the phase function because of multiple scatterings, in principle, one could extend the networks to have dynamic phase functions that vary with incoming frequency for a self-consistent quantum mechanical treatment or polarisation using a Stokes vector approach as in Seon, Song, & Chang (2022), but this is beyond the scope of the current study.

For the applications of this paper we use a static, uniform-density configuration with either slab or spherical geometry. We inject photons at the centre of the medium (point source) with an initial frequency chosen such that the photon is equally likely to be in any of the lines, i.e. we match line i with probability weight $\omega_i$ . For convenience, we calculate $a\tau_0$ using the $a_i$ from the line with lower frequency, unless otherwise specified. In many cases, the $a_i$ used can be arbitrary, as $a\tau_0$ is simply a frequency diffusion scaling factor. For concreteness we list the specific parameter values used in each figure in Appendix C for reproducibility.

For optically-thick tests, we have confirmed that the standard constant core-skipping algorithm works well as long as $x_{\text{crit}}$ for each independent line is not too large, where $x_{\text{crit}}$ is the minimum perpendicular velocity drawn from a volcano–Gaussian using the Box–Muller method (Ahn, Lee, & Lee 2002). In general, there is a negligible impact on the emergent spectra for constant values $x_{\text{crit}} \ll (a\tau_0)^{1/3}$ or dynamical core skipping as in Smith et al. (2015). However, for certain tests where the line separation decreases, such as the Ly $\alpha$ and deuterium cases, core skipping causes an unphysical bias towards one of the lines if they have different optical depths. Therefore, to show the accuracy of the analytical solutions, we turned off core skipping for such runs. To remedy this, one must generalise the dynamic core-skipping algorithm to use the minimum of $0.2\,(a\tau_0)_i^{1/3}$ for the line and the frequency offset placing the photon in the influence of nearby neighbouring lines, e.g. $k_{x,i} \lesssim k_x/2$ , but we leave this for a future study.

In Figure 4, we present results for the two-line V-shaped network with varying frequency separations between the lines, comparing these with our analytical solutions from Section 2. As we decrease the line separation, the two distinct double-peaked solutions merge together like superimposed solutions. The physical reason for this is that once the line separation is $\Delta \lesssim (a\tau_0)^{1/3}$ the resonant broadening is enough to diffuse into the range of the other line, leading to mixing during random walk transport. We note that the flux approaches zero at each line centre because of the $\bar{H}(x)$ in the denominator has orders of magnitude higher cross-sections line centre, so frequency diffusion is required for escape. We also see that when the line separation becomes very small, the solution starts to recover the standard double-peaked structure of a singular line. The analytic solution is indistinguishable from the converged simulation output, confirming the validity of the Fokker–Planck approximation in high optical depth scenarios.

Figure 4. A toy network with two arbitrary lines of equal oscillator strength, illustrating the transition from well-separated lines with distinct peaks to overlapping lines with merged profiles as the frequency separation is reduced. We compare emergent spectra for slab geometry with different separations between the lines (top: widely separated; middle: intermediate; bottom: nearly overlapping). The analytic solution from Eq. (13) and numerical colt simulations are in excellent agreement even as the lines become highly coupled to produce nontrivial spectra.

Finally, while we do have an analytical form of $\tilde{x} = \tilde{x}(x)$ as solved for in Appendix B, for this paper we converted between $\tilde{x}$ and x space by numerically integrating Eq. (B7). Despite the very strong agreement between the closed form solution and and numerical calculation, we chose to use the latter to avoid assumptions or approximations breaking down, and since it was computationally inexpensive.

5. Applications to real-world systems

We now turn to two specific astrophysical applications to demonstrate V-shaped multi-line RT: Ly $\alpha$ fine-structure and deuterium injection. In both cases, multi-line effects are known to play a role, and our analytic solutions provide a convenient way to isolate these. The methods are also adaptable to other such systems.

5.1. Ly $\boldsymbol{\alpha}$ fine structure

Relativistic corrections to the Schrödinger equation and spin–orbit coupling causes the 2p energy level of the hydrogen atom to split into two distinct energy states. We ignore hyper-fine splitting of the ground state. This modifies the standard Ly $\alpha$ network, creating two closely spaced lines corresponding to the following atomic transitions: $1{}^2S_{1/2} \longleftrightarrow 2{}^2P^0_{1/2}$ , and $1{}^2S_{1/2} \longleftrightarrow 2{}^2P^0_{3/2}$ , which we denote as $1/2$ and $3/2$ , respectively. The latter has twice the oscillator strength of the former, due to the J degeneracy, but their wavelengths differ by $1.34\,\text{km s}^{-1}$ . While fine-structure splitting is a more physically accurate picture of Ly $\alpha$ RT, most models ignore it because the separation becomes smaller than the Doppler broadening at gas temperatures higher than 100 K (Seon & Kim 2020). Since there is one lower state and two upper states, we can apply our solutions to this V-shaped two-line network. Inputting the relevant atomic data, we have $n_{1/2} = n_{3/2} = n_{H}$ , and because both lines stem from hydrogen and $\nu_{1/2} \approx \nu_{3/2}$ we have $\Delta\nu_{{D}, 1/2} \approx \Delta\nu_{{D}, 3/2}$ , such that $\mu_1 \approx \mu_2$ ; similarly since the Einstein coefficients are approximately equal this leads to damping parameters $a_{1/2} \approx a_{3/2}$ . However, the weights differ according to $\omega_{1/2} \approx \frac{1}{2}\omega_{3/2} = 1/3$ due to the ratio of the oscillator strengths.

In Figure 5, we present results from idealised cold clouds in slab geometry at three temperatures, $T = \{1, 10, 100\}\,\text{K}$ , demonstrating excellent agreement between our analytic solution and numerical simulations. As temperature decreases, the thermal Doppler width decreases, exposing the fine-structure splitting of the Ly $\alpha$ line. However, even at 1 K, the emergent spectra strongly resembles non-split Ly $\alpha$ because the peak locations are still much larger than the line separations $\Delta \ll (a\tau_0)^{1/3}$ , blending the line interactions. It is important to note that our base solution becomes physically inaccurate once we consider the effect of recoil, as Ly $\alpha$ photons systematically lose energy to absorbing atoms. This would add a term to the redistribution function yielding (Rybicki 2006)

(16) \begin{align} &\int J_{x'} H(x'_i)R_{x'-\Delta_i \rightarrow x-\Delta_i} \text{d}x' \nonumber \\[4pt] &= H(x_i)J_{x} + \frac{1}{2}\frac{\partial}{\partial x}\left[H(x_i)\left( \frac{\partial J_{x}}{\partial x} +\frac{\Delta\nu_{D} hJ_x}{k_{{B}}T}\right)\right] \, .\end{align}

Because of the inverse temperature term in the correction, recoil becomes more significant in cold gas, which matches the result of our simulations. To our knowledge, the only existing analytical solution including the effects of recoil is derived in Nebrin et al. (2025), but this is only for a single line treatment and produces a series solution with coefficients defined recursively. We now outline how recoil might be incorporated into our multi-line framework. The full RTE becomes

(17) \begin{align} &\frac{\boldsymbol{\nabla} }{k(\boldsymbol{r})} \boldsymbol{\cdot} \left(\frac{\boldsymbol{\nabla} J}{k(\boldsymbol{r})} \right) + \frac{3}{2}\sum_{i=0}^N\frac{\Delta\nu_{{D},i}^2}{\Delta\nu_{D}^2} H(x_i) \frac{\partial}{\partial x}\!\left[H(x_i)\!\left(\frac{\partial J}{\partial x} + \frac{h\Delta\nu_{{D},i}}{k_{B}T}J\right)\right] \nonumber \\[4pt] &= \frac{3\bar{H}(x)}{k(\boldsymbol{r})}\int \frac{j_x}{4\pi} \, \text{d} \Omega \, . \end{align}

We note that Nebrin et al. (2025) also includes other effects such as gas velocity gradients and continuum absorption/scattering, but as a proof of concept we neglect these effects and save them for a more comprehensive future study.

Figure 5. A physical Ly $\alpha$ fine-structure setup where we vary the temperature of the system. Even for very low gas temperatures the two lines do not cleanly separate in the emergent spectrum. We compare emergent Ly $\alpha$ spectra for spherical geometry at $T = \{1, 10, 100\}\,\text{K}$ . Decreasing the temperature reduces the Doppler width and causes the line separation to increase. Even at 1 K, however, the two components overlap strongly and the profiles remain double-peaked rather than fully splitting. Atomic recoil becomes important at extremely low temperatures.

Although we do not incorporate recoil with full self-consistency, as an experiment in the regimes where a fine-structure treatment is needed, we combine our principle of superposition with the solution in Nebrin et al. (2025) to get

(18) \begin{align} J(x) = \frac{\sqrt{18}\mathcal{L}}{(4\pi)^2\tau_0\Delta\nu_{D}R\bar{H}(x)} \sum_{n=1}^\infty \frac{2(-1)^n n\pi}{9(\mathcal{R}_+ - \mathcal{R}_-)} \sum_{i=0}^N \mathcal{F}\!\left(\frac{x-\Delta_i}{\beta_n}\right) ,\end{align}

where $\mathcal{R}_+, \mathcal{R}_-$ , and $\mathcal{F}(x/\beta_n)$ are defined in appendix A2 in Nebrin et al. (2025). The results are illustrated in Figure 6, where we used the solutions for a fine-structure network (top) and a toy network with larger line separation (bottom) in the low temperature regime $T = 10$ K. ^a Naively applying the superposition yields results better than the single line case, but clearly the solutions are not sufficient, as evident in the large line separation case where the solutions produce the wrong peak locations and heights. A proper treatment would require generalising the single-line wing approximation $H(x,a) \approx a/\sqrt{\pi} x^2$ used in the derivation of this solution, and most likely require a double-wing or Wentzel–Kramers–Brillouin (WKB) approximation to obtain an analytically tractable differential equation. This is a formidable task with relatively little physical benefit, so we leave this for future studies.

Figure 6. We compared the superimposed recoil solution with our non-recoil solution, along with colt simulations. Our high optical depth ( $a\tau_0 = 10^5$ ), low temperature setup ( $T = 10$ K) ensures convergence of analytical solutions and strong recoil effects. Top: Eq. (18) is used for fine-structure, and we see strong agreement. Note that the peaks are not perfectly aligned, an issue that is exacerbated for larger line separations. Bottom: Eq. (18) is used for an arbitrary setup with larger line separation, and we observe misaligned peaks as well as incorrect heights, likely indicating that the line-specific parameters like optical depth or damping parameters are poorly modelled. In both panels, the qualitative features match, as we see the photons biased towards the red side of the plot, indicating that the true solution likely resembles this one with some non-linear correction.

In summary, our analytic solution is accurate and fine-structure splitting does not significantly alter emergent Ly $\alpha$ spectra unless the gas is extremely cold (although there may be other consequences for internal properties, such as the number of scatterings), at which point recoil must also be included.

5.2. Deuterium injection

We also consider the effect that a small deuterium (D) fraction on Ly $\alpha$ RT. ${D}\,$ i has a Ly $\alpha$ transition nearly identical to that of hydrogen but slightly shifted in frequency due to the greater mass. Deuterium injection has been implemented into Monte Carlo codes (Dijkstra et al. 2006; Michel-Dansac et al. 2020), and while this is not exactly a V-shaped network since the lower energy levels are not the same energy, the lines are independent, so they can still be modelled with the same solution as before. Qualitatively, deuterium produces a narrow absorption dip on the blue side of the Ly $\alpha$ line. We use the ratio of deuterium atoms to hydrogen atoms as $n_{D}/n_{H} = 3.39 \times 10^{-5}$ , obtained from Burles & Tytler (1998). As described in Dijkstra et al. (2006), the deuterium line centre frequency is shifted blueward by 82 $\text{km s}^{-1}$ . We choose to centre x at Ly $\alpha$ line centre, so we have $x = (\nu - \nu_{0,{H}})/\Delta\nu_{{D},{H}}$ . This corresponds to $\Delta_{H} = 0$ and $\Delta_{D} = +6.3\,T_4^{-1/2}$ . Due to the doubled mass relative to hydrogen, factors of $\sqrt{2}$ must be accounted for, yielding $\mu_1 = 1$ and $\mu_2 = \sqrt{2}$ . Deuterium, like hydrogen, will also encounter fine-structure splitting. This detail was ignored in previous analyses, but we have included it in our numerical and analytical implementation. Therefore, our network consists of four lines, with two fine structure Ly $\alpha$ lines each. Adopting the values of ${D}/{H}$ , $f_{lu,i}$ , and $\Delta\nu_{{D},i}$ for Ly $\alpha$ and deuterium into Eq. (5), we get $\omega_{{H},1/2} = 0.333318$ , $\omega_{{H},3/2} = 0.666634$ , $\omega_{{D},1/2} = 1.59726 \times 10^{-5}$ , and $\omega_{{D},3/2} = 3.1944\times 10^{-5}$ . These weights indicate that the emergent spectra consists almost entirely of H Ly $\alpha$ emission.

In the top panel of Figure 7, we reproduce Figure 3 in Dijkstra et al. (2006), and compare our analytical solution to the emergent spectra. The new solution qualitatively predicts the absorption feature caused by deuterium injection, albeit with a different peak height blueward of the deuterium line. However, $a_{H}\tau_{{H},0} = 343$ , far less than the $a_{H}\tau_{{H},0} \gtrsim 10^3$ needed for analytical convergence of our solutions (Neufeld 1990; McClellan et al. 2022). In the bottom panel of Figure 7, we increase the optical depth to $a_{H}\tau_{{H},0} = 10^4$ and lower the temperature to 987 K, which improves convergence in both lines, since $a_{D}\tau_{{D},0}$ is also around 30 times higher than the previous test. However, we still have $a_{D}\tau_{{D},0} \ll 10^3$ so the deuterium line itself is not in the diffusion regime. Nevertheless, the overall agreement is good, and the analytic model captures the multi-line physics in sufficiently optically-thick media.

Figure 7. Emergent spectra of Ly $\alpha$ with deuterium injection, showing the analytical solutions give strong qualitative agreement with simulations, and converge in the regime $a\tau_0 \gtrsim 10^3$ . The values of $a\tau_0$ and T were specially chosen to align the deuterium absorption feature with the peak of the Ly $\alpha$ spectrum. Top: We reproduce the setup in Figure 3 in Dijkstra et al. (2006) with $T = 10^4$ K and moderate optical depth $a\tau_0 = 343$ . The differences are due to the relatively low optical depth. Bottom: With a higher optical depth $a\tau_0 = 10^4$ and lower temperature of $T = 987$ K to align the absorption feature (which is also more narrow), we see much better agreement between our analytical solution and the emergent flux.

To determine when the deuterium line will have an observable impact in the Ly $\alpha$ spectrum, we demand that $x_{p} \approx 0.93099\,(a_{H}\tau_{{H},0})^{1/3} \approx \Delta_{D} = +6.3(T_4)^{-1/2}$ , where $x_{p}$ is the peak of the Ly $\alpha$ solution (Lao & Smith 2020). Therefore, increasing the optical depth requires decreasing the temperature to align the absorption feature with the emergent flux. This relationship is shown in Figure 8, where one can see which values of temperature (T) and hydrogen column density ( $N_{{H}_{I}}$ ) cause the deuterium line centre to lie near the peak of Ly $\alpha$ solution. For a given optical depth, the alignment condition $\Delta_{D} = x_{p}$ yields $T = 2.77\times 10^7 \, (N_{{H}_{I}}/10^{20}\,\text{cm}^{-2})/(a\tau_0)$ . Around this ‘sweet spot’ region in parameter space, one must account for the effects of deuterium injection, while outside this range the absorption is either too close to the Ly $\alpha$ core or beyond the reach of resonant scattering. We caution that in realistic astrophysical environments, clouds have bulk motions (infall, outflow, and turbulence) that broaden and shift the Ly $\alpha$ line, likely smearing out deuterium absorption features across frequency space, making deuterium-induced signatures extremely challenging to observe. Our model of a static uniform cloud can be regarded as an optimistic upper limit on the impact of deuterium. Incorporating velocity gradients into the multi-line framework, building on Nebrin et al. (2025) and Smith et al. (2025), would allow one to quantify the interplay of bulk motions and deuterium injection. We can also compare how deuterium affects the force multiplier with (Lao & Smith 2020):

(19) \begin{align} M_{{F},{D}} &= \frac{-c}{3\mathcal{L}}\int \nabla u(r)\,\text{d}V = \frac{-c}{3\mathcal{L}}\int \frac{4\pi}{c}\int \nabla J(x) \,\text{d}x\text{d}V \nonumber \\[4pt] &= \frac{2\sqrt{6}}{3\pi}\sum_{n=1}^{\infty} \frac{[1-(-1)^n]}{n}\int \sum_{i=1}^4\left[\omega_i\mu_ie^{-\lambda_n |\tilde{x}-\tilde{\Delta}_i|} \right] \, \text{d}x(\text{d}\tilde{x}) \, . \nonumber \\[-12pt] \end{align}

Unfortunately, the non-trivial change of variables means that our expression for $M_{F}$ is not closed form, but we can still compute this numerically. We also calculate the Ly $\alpha$ force multiplier using Eq. (19) but with the single line treatment $N = 1, \,\bar{H}(x) \rightarrow H(x,a), \, \text{d}\tilde{x} = \sqrt{2/3}\,\text{d}x/[\tau_0H(x,a)]$ as in Lao & Smith (2020). We show these quantities in Figure 9, and find that the difference is below two percent for most physical scenarios, except for low temperatures, where the low- $a\tau_0$ difference can be as high as eight percent. Therefore, deuterium injection does not change any prior conclusions about Ly $\alpha$ feedback, and another potential source of uncertainty in such modelling is removed (Nebrin et al. 2025).

Figure 8. Deuterium absorption in column density and temperature parameter space. The black dotted line shows where the deuterium line core lies exactly on the peak of the Ly $\alpha$ flux, and the grey dotted lines indicate $\pm0.5\,(a\tau_0)^{1/3}$ away. The Ly $\alpha$ flux as a fraction of its peak is heat-mapped in the background.

Figure 9. We numerically evaluate the force multiplier for hydrogen $(M_{{F},\,{H}})$ and for deuterium $(M_{{F},\,{D}})$ using Eq. (19), and show that for most relevant astrophysical scenarios, the difference between the two is at the percent level. As temperature or optical depth increases, the force multipliers become equal.

6. Discussion

We have developed and tested the first (to our knowledge) analytic framework for resonant RT involving multiple spectral lines simultaneously. Focusing on V-shaped atomic systems with one ground state and multiple excited states, i.e. coupled transitions, we derive closed-form solutions for the angle-averaged intensity in both infinite slab and spherical geometries. These solutions generalise previous single-line treatments (e.g. Harrington 1973; Neufeld 1990; Lao & Smith 2020; Nebrin et al. 2025) by showing that, under the same diffusion and high optical-depth approximations, the multi-line problem can be expressed as a weighted linear superposition of single-line solutions in an appropriately transformed frequency coordinate. This represents a powerful and somewhat unexpected conceptual advance: the complicated coupling between lines is encoded entirely in the composite absorption profile $\bar{H}$ , frequency variable transformation $\tilde{x}(x)$ , and weighting coefficients $\omega_i$ and $\mu_i$ , making multi-line transfer analytically tractable for the first time.

Our comparisons with MCRT simulations reinforce the validity of the analytic solutions. By extending the colt code to handle multi-line scattering, we were able to simulate idealised and physical scenarios with overlapping lines and directly verify that the intensity profiles match the numerical ones to within the sampling noise. The agreement holds across a broad range of line separations, optical depths, and temperatures, as long as we remain in the regime where the diffusion approximation is valid, i.e. $a\tau_0 \gtrsim 10^3$ for each line. In cases where one line has a much lower optical depth than this, such as with trace deuterium injection, the solution can still predict qualitative effects but may not capture peak heights exactly. Still, considering the deuterium optical depths are often orders of magnitude below the expected criterion for applicability, we consider the multi-line solutions to be successful. This strong agreement confirms that the generalised absorption coefficient, modified Fokker–Planck approximation, and generalised change of variables accurately capture the underlying physics of resonant scattering and that the newly implemented multi-line extension to colt behaves as expected. Furthermore, the Monte Carlo method and diffusion solvers can now be applied to more complex problems beyond what analytic solutions can handle, providing a numerical laboratory for future studies.

The Ly $\alpha$ fine-structure application highlights both the strengths and limitations of the formalism. On one hand, the model correctly recovers the merging of the $2{}^2P^0_{1/2}$ and $2{}^2P^0_{3/2}$ components with increasing temperature. On the other hand, at very low temperatures the neglect of atomic recoil introduces a systematic discrepancy, since recoil adds an additional drift term in frequency space that becomes significant, producing asymmetric red and blue peaks when Doppler broadening is minimal ( $T \lesssim 100$ K). In an initial attempt to account for recoil, we also attempted to use the superposition principle on the solution of Nebrin et al. (2025), turning their single-line results into a multi-line solution. However, the agreement with MCRT simulations is poor, especially as the line separation increases. Future work could involve relaxing the approximations used in the single-line derivation Nebrin et al. (2025) to generate generalised solutions that include more complex effects like recoil, velocity gradients, or dust absorption. One could also consider incorporating polarisation, which may have interesting multi-line effects.

The deuterium injection case demonstrates the versatility of the approach for realistic astrophysical multi-species mixtures. Despite being outside the strict V-shaped configuration, the photons see two independent resonances, so the method successfully predicts the deuterium absorption feature and its evolution with optical depth, reproducing the qualitative behaviour seen in Dijkstra et al. (2006). By surveying parameter space, we were able to identify where deuterium is most likely to matter along a narrow region of column density and temperature. The improved agreement at higher optical depths indicates that the analytic formalism converges to the exact solution in the same diffusion-dominated limit as single-line transfer. This opens the possibility of applying the model to other chemically or isotopically mixed media where multiple resonant transitions overlap, such as dense nebular regions, metal-enriched outflows, or partially ionised circumgalactic gas. We also find a negligible impact on the force multiplier, which means that recent results indicating the importance of Ly $\alpha$ radiation pressure feedback, but neglecting deuterium scattering, remain unchanged (e.g. Smith, Bromm, & Loeb 2017; Tomaselli & Ferrara 2021; Nebrin et al. 2025; Manzoni & Ferrara 2025). ^b However, once again the model reaches its limits, as bulk velocity motions and turbulence can smear the deuterium absorption feature across frequency space, an effect not accounted for in this paper. Future studies could straightforwardly incorporate details from Nebrin et al. (2025) and Smith et al. (2015). Each additional effect adds complexity, but our work provides a baseline solution upon which others can explore further.

Beyond the specific applications we examined, the present work provides the foundation for a more general theory of multi-line RT networks. The recent work of Chang et al. (2025), who performed MCRT modelling of Rayleigh and Raman scattering, and branching transitions in between, underscores the importance of developing accurate multi-line RT methods. Their simulations show that multi-line scattering can produce non-intuitive line ratios and broadening that challenge simple interpretations of observed spectra. Extending our framework to include branching transitions would enable analytic exploration of systems where photons can scatter between multiple excited states before de-excitation, such as stellar outflows, accreting compact binaries, and Little Red Dots (Harikane et al. 2023; Kocevski et al. 2023; Greene et al. 2024; Matthee et al. 2024; Chang et al. 2025). For instance, $\Lambda$ -shaped networks with two lower levels and one upper level or a cascade where an excited state can decay via multiple paths (as in more complete hydrogen treatments) would be natural next targets. Initial explorations are promising, but require a more careful treatment of generalised frequency redistribution functions. Bridging the gap between detailed numerical RT and analytic understanding will enable us to more effectively interpret the ever-increasing wealth of spectral data revealing unexpected phenomena. Such results establish a benchmark for future studies of multi-line resonant transfer and supply an analytic counterpart to numerical codes, facilitating rapid exploration of parameter space and deeper physical insight into the coupling of radiation and matter across multiple spectral lines.