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. Reference Dijkstra, Prochaska, Ouchi and Hayes2019). 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 Reference Hummer1962). 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 Reference Ahn and Lee2015; Seon & Kim Reference Seon2020):

(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 (Reference Lao and Smith2020), 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 (Reference Lao and Smith2020), 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 (Reference Harrington1973), Neufeld (Reference Neufeld1990), and Dijkstra et al. (Reference Dijkstra, Haiman and Spaans2006) 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 (Reference Lao and Smith2020) 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 (Reference Lao and Smith2020) 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. Reference Smith, Safranek-Shrader and Bromm2015, Reference Smith, Ma, Bromm, Finkelstein and Hopkins2019, Reference Smith2022; McClymont, Smith, & Tacchella Reference McClymont, Smith and Tacchella2025, 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 (Reference Seon2020) 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. Reference Sansonetti, Kerber, Reader and Rosa2004) 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. Reference Dijkstra, Haiman and Spaans2006), 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 (Reference Seon, Song and Chang2022), 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 Reference Ahn, Lee and Lee2002). 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. (Reference Smith, Safranek-Shrader and Bromm2015). 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 Reference Seon2020). 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 Reference Rybicki2006)

(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. (Reference Nebrin, Smith and Lorinc2025), 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. (Reference Nebrin, Smith and Lorinc2025) 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. (Reference Nebrin, Smith and Lorinc2025) 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. (Reference Nebrin, Smith and Lorinc2025). 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.Footnote ^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. Reference Dijkstra, Haiman and Spaans2006; Michel-Dansac et al. Reference Michel-Dansac, Blaizot, Garel, Verhamme, Kimm and Trebitsch2020), 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 (Reference Burles and Tytler1998). As described in Dijkstra et al. (Reference Dijkstra, Haiman and Spaans2006), 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. (Reference Dijkstra, Haiman and Spaans2006), 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 Reference Neufeld1990; McClellan et al. Reference McClellan, Davis and Arras2022). 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. (Reference Dijkstra, Haiman and Spaans2006) 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 Reference Lao and Smith2020). 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. (Reference Nebrin, Smith and Lorinc2025) and Smith et al. (Reference Smith, Lorinc, Nebrin and Lao2025), 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 Reference Lao and Smith2020):

(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 (Reference Lao and Smith2020). 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. Reference Nebrin, Smith and Lorinc2025).

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.