<figure><img src="image_1.png" /></figure>

HAL
open science

# Full non-LTE multi-level radiative transfer

T. Lagache, F. Paletou, M. Sampoorna

► To cite this version:

T. Lagache, F. Paletou, M. Sampoorna. Full non-LTE multi-level radiative transfer. Astronomy & Astrophysics - A&A, 2025, 699, pp.A198. 10.1051/0004-6361/202555008. hal-05175595

HAL Id: hal-05175595

https://hal.science/hal-05175595v1

Submitted on 22 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Full non-LTE multi-level radiative transfer

## I. An atom with three bound infinitely sharp levels

T. Lagache<fnref n="1" />,<fnref n="2" />,*, F. Paletou<fnref n="1" />,<fnref n="2" />, and M. Sampoorna<fnref n="3" />

Received 2 April 2025 / Accepted 20 May 2025

## ABSTRACT

**Context.** The standard nonlocal thermodynamic equilibrium (non-LTE) multi-level radiative transfer problem only takes into account the deviation of the radiation field and atomic populations from their equilibrium distribution.

**Aims.** We aim to show how to solve for the full non-LTE (FNLTE) multi-level radiative transfer problem, also accounting for deviation of the velocity distribution of the massive particles from Maxwellian. We considered, as a first step, a three-level atom with zero natural broadening.

**Methods.** In this work, we present a new numerical scheme. Its initialisation relies on the classic, multi-level approximate A-iteration (MALI) method for the standard non-LTE problem. The radiative transfer equations, the kinetic equilibrium equations for atomic populations, and the Boltzmann equations for the velocity distribution functions were simultaneously iterated in order to obtain self-consistent particle distributions. During the process, the observer's frame absorption and emission profiles were re-computed at every iterative step by convolving the atomic frame quantities with the relevant velocity distribution function.

**Results.** We validate our numerical strategy by comparing our results with the standard non-LTE solutions in the limit of a two-level atom with Hummer's partial redistribution in frequency, and with a three-level atom with complete redistribution. In this work, we considered the so-called cross-redistribution problem. We then show new FNLTE results for a simple three-level atom while evaluating the assumptions made for the emission and absorption profiles of the standard non-LTE problem with partial and cross-redistribution.

**Key words.** line: formation – line: profiles – radiative transfer – stars: atmospheres

## 1. Introduction

The standard multi-level non-LTE (non-local thermodynamic equilibrium) radiative transfer problem is at the very heart of the interpretation of spectra from astrophysical objects such as stars, circumstellar discs, molecular clouds, and so on (see Hubený & Mihalas 2014). While the standard approach, which has been continuously developed since the late 1960s, accounts for deviations of the radiation field and the atomic populations from their equilibrium distributions, it assumes velocity distribution of radiating atoms to be Maxwellian.

Hereafter, we consider the problem of full non-LTE (FNLTE) radiative transfer, which mainly consists of the self-consistent calculation of the velocity distribution functions (VDFs) of massive particles with the radiation field; i.e. the distribution of photons in the atmosphere. In that context, the main numerical burden is to work directly with atomic velocities and therefore solve the kinetic equilibrium equations for each velocity.

The kinetic theory of particles and photons (see Oxenius 1986) on which our approach is based also introduces an additional physical process of atomic velocity-changing collisions, which appear naturally in the Boltzmann equations. Hubený & Cooper (1986, see also; e.g. Hubený et al. 1983a,b) studied the consequences of these collisions on the radiative transfer problem in detail (see also the discussion in Sect. 4 of

Sampoorna et al. 2024). The multi-level problem was originally formulated by Hubený et al. (1983a,b, hereafter HOSI and HOSII, respectively; see also Oxenius 1965, 1986) using a semi-classical picture.

Until now, only the two-level FNLTE problem has been solved by Paletou et al. (2023, hereafter PSP23; see also Lagache et al. 2025) for infinitely sharp levels, and by Sampoorna et al. (2024) for naturally broadened excited levels. In the FNLTE multi-level radiative transfer problem, the Boltzmann equation for the VDFs of excited atoms need to be solved simultaneously with the radiative transfer equation for all radiatively allowed transitions. Furthermore, at every iteration we also need to recompute all the macroscopic absorption and emission profiles by convolving the atomic profiles with the relevant VDFs. Clearly, this is a numerically very challenging problem. As a first step, we considered a simple atomic model, namely a three-bound-level atom with zero natural broadening, using most data from Avrett (1968).

The FNLTE formalism for multi-level atoms naturally accounts for both resonance and Raman scattering contributions<fnref n="1" />. To incorporate these scattering mechanisms – including the effects of partial redistribution (PRD) for resonance scattering and cross-redistribution (a.k.a. XRD) for Raman scattering – in the standard (in the sense that VDF deviations from Maxwellian are only partially considered) non-LTE multi-level

<footnote n="1">¹ Université de Toulouse, Observatoire Midi-Pyrénées, Cnrs, Cnes, Irap, Toulouse, France</footnote>

<footnote n="2">² Cnrs, Institut de Recherche en Astrophysique et Planétologie, 14 av. E. Belin, 31400 Toulouse, France</footnote>

<footnote n="3">³ Indian Institute of Astrophysics, Koramangala, Bengaluru 560034, India</footnote>

<footnote n="1">¹ In a scattering transition from level $a \to b \to c$, if $a = c$ we have resonance scattering; if $a \neq c$ and $E_a, E_c < E_b$ it is Raman scattering.</footnote>

<footnote>* Corresponding author: tristan.lagache@univ-tlse3.fr</footnote>
problem an approximate approach was adopted; namely, the absorption profile for all transitions was assumed to be the usual Voigt (or Doppler) profile, and a simplified form of the emission profile derived by HOSII was used (see e.g. Hubený 1985; Uitenbroek 1989; Paletou & Auer 1995; Sampoorna et al. 2013). We remark that the simplifying assumptions behind this approximate approach were detailed in Hubený (1985), and their validity remains to be evaluated. However, none of these assumptions are used in our more complete FNLTE approach.

This paper is organised as follows. In Sect. 2, we present the formalism of the FNLTE transfer problem and apply it to the case of a three-bound-level atom with zero natural broadening in Sect. 3. Our new numerical scheme is described in Sect. 4. It is validated against various benchmark solutions in Sect. 5 (and Appendix C). Supplementary material, useful for understanding our tests, is presented in the Appendix B. New results are presented in Sect. 6. Finally, conclusions and future developments are discussed in Sect. 7.

## 2. The FNLTE multi-level problem

Solving a multi-level radiative transfer problem requires the self-consistent resolution of a set of equations, each of them describing the distribution of particles of different nature: photons and massive particles (atoms, ions, molecules, or electrons). Here, we used the kinetic approach developed by Oxenius (1986, 1965), which was extended by HOSI and HOSII and recently dubbed FNLTE by Paletou & Peymirat (2021). In this formalism, the equations of the transfer problem are all kinetic equations.

### 2.1. The kinetic equations of photons

The distribution of photons associated with the transition from level i of energy $E_i$, to a level j of energy $E_j$, is the specific intensity $I_{ij}$. For an unpolarised, time-independent problem in a plane-parallel 1D geometry, this distribution depends on the frequency $v_{ij}$ of the photon in the observer's frame, its direction cosine $\mu_{ij} = \cos(\theta_{ij})$, and the optical depth $\tau_{ij}(v_{ij})$. Its kinetic equation is the usual radiative transfer equation (RTE), which is written, with the convention $i < j$, as follows:

$$ \mu_{ij} \frac{\partial I_{ij}(v_{ij}, \mu_{ij}, \tau_{ij})}{\partial \tau_{ij}} = I_{ij}(v_{ij}, \mu_{ij}, \tau_{ij}) - \frac{\eta_{ji}(v_{ij}, \tau_{ij})}{\chi_{ij}(v_{ij}, \tau_{ij})}. \quad (1) $$

Hereafter, we will not display dependence of relevant quantities on the optical depth $\tau_{ij}$ without loss of clarity. In the RTE, we introduced the emissivity $\eta_{ji}$ and the absorption coefficient $\chi_{ij}$. They depend on the frequency, the optical depth, the direction of propagation and, implicitly, the VDFs of the massive particles, knowledge of which is indispensable for the computation of the emission and absorption profiles.

In the following, we assume the isotropy of these quantities<fnref n="2" />, and we write

$$ \eta_{ji}(v_{ij}) = \frac{hv_{0,ij}}{4\pi} n_j A_{ji} \psi_{ji}(v_{ij}), \quad (2) $$

and

$$ \chi_{ij}(v_{ij}) = \frac{hv_{0,ij}}{4\pi} n_i B_{ij} \varphi_{ij}(v_{ij}), \quad (3) $$

where $A_{ji}$ and $B_{ij}$ are the Einstein coefficients associated, respectively, with spontaneous emission and radiative absorption; $v_{0,ij}$ is the frequency of the $i \leftrightarrow j$ transition, and $n_i$ and $n_j$ are the population densities of the $i$ and $j$ excited levels. Strictly speaking, we should also consider the stimulated emission characterised by the Einstein coefficient $B_{ji}$. However, in most astrophysical cases (with the exception of masers; see e.g. Gray 2012), $n_j B_{ji} \ll n_i B_{ij}$, and stimulated emission is therefore considered as 'negative absorption'. Also, taking this additional phenomenon into account makes it much more difficult to obtain atomic frame absorption and emission profiles (see Sect. 5 of HOSI). Therefore, as a first step, we neglected stimulated emission completely. We also introduced the absorption profile $\varphi_{ij}$ and the emission profile $\psi_{ji}$ in the observer's frame defined as

$$ \varphi_{ij}(v_{ij}) = \oint \frac{d\Omega_{ij}}{4\pi} \int f_i(v) \alpha_{ij}(\xi_{ij}) d^3 v, \quad (4) $$

and

$$ \psi_{ji}(v_{ij}) = \oint \frac{d\Omega_{ji}}{4\pi} \int f_j(v) \beta_{ji}(\xi_{ji}) d^3 v, \quad (5) $$

where $\beta_{ji}^3$ and $\alpha_{ij}$ are the atomic emission and absorption profiles; $f_i(v)$ and $f_j(v)$ are the VDFs, normalised to unity for, respectively, atoms in the excited level $i$ and $j$; $v$ is the atomic velocity; $\Omega_{ij}$ is the direction of a photon of frequency $v_{ij}$ depending on azimuth $\phi_{ij}$ and co-latitude $\theta_{ij}$; and $\xi_{ij}$ is the photon frequency in the atomic frame. The latter is directly linked to the frequency in the observer's frame by the Fizeau-Doppler relationships:

$$ \xi_{ij} = v_{ij} - \frac{v_{0,ij}}{c} v \cdot \Omega_{ij}. \quad (6) $$

Solving the RTE and therefore the radiative transfer problem requires knowledge of the dynamic and radiative properties of the atoms: (i) their velocity distribution, $f_i$, (ii) the atomic densities, $n_i$, and (iii) the atomic profiles, $\alpha_{ij}$ and $\beta_{ji}$. Each of these elements are discussed below.

### 2.2. The kinetic equations of massive particles

Massive particles constituting an atmosphere are also described, as photons, by a distribution (in velocity) governed by a kinetic equation. Apart from their nature, the only properties that distinguish each of these atoms are their speed and direction. Each of them is therefore characterised by its own VDF. The standard radiative transfer, however, deals with the self-consistent determination of the number density of the atoms – i.e. the first moment of these VDFs – together with the radiation field so that nothing is known about these VDFs. On the contrary, our more complete FNLTE approach requires the self-consistent resolution of all the kinetic equations describing all the particles distributions.

In the following, we assume that the temperature of the medium is sufficiently small so that $n_1 \gg n_2, n_3, ...$ In that case, atoms are so numerous on the fundamental level that their population varies very slightly with respect to the LTE value, and we consider that their distribution is Maxwellian. We only considered bound-bound radiative and collisional processes. For atoms

<footnote n="2">² Taking into account their angular dependence would require us to consider angular redistribution.</footnote>

<footnote n="3">³ The atomic emission profile was previously written, $\eta$ (see e.g. Oxenius 1986 or HOSI and HOSII), but we renamed it $\beta$ in order to avoid any confusion with the emissivity.</footnote>
in excited levels, let $F_i = n_i f_i$ be the distribution associated with the excited level $i \ge 2$. The kinetic equation describing their distribution is the Boltzmann equation (we neglected the streaming of particles as argued by Hubený 1981; for more details, please refer to Paletou & Peymirat 2021):

$$C[F_i] = \left(\frac{\delta F_i}{\delta t}\right)_{\text{rad.}} + \left(\frac{\delta F_i}{\delta t}\right)_{\text{inel.}} + \left(\frac{\delta F_i}{\delta t}\right)_{\text{v.c.c.}} = 0. \quad (7)$$

Here, C is the collisional operator composed of three terms of different physical origins: a radiative term (collisions with photons), an inelastic term (collisions with free electrons), and a velocity-changing collision term, i.e. collisions between atoms (this nomenclature was discussed in detail in Sect. 4 of Sampoorna et al. 2024).

Let $N_i = n_i f_i d^3 v dV$ be the number of atoms at level $i$ contained in a volume $d^3 v dV$ of the phase space. Per second, $N_i A_{ij}$ atoms spontaneously emit a photon and deexcite towards a level $j < i$, and $N_i B_{ji} J_{ij}(v)$ atoms absorb a photon and are excited to a level $j > i$. In the phase space, these atoms are destroyed. Here, $J_{ij}(v)$ denotes the partial scattering integral for atoms at level $i$ and velocity in the range of $(v, v + dv)$, which is defined as

$$J_{ij}(v) = \int \alpha_{ij}(\xi_{ij}) I_{ij} d\xi_{ij}. \quad (8)$$

Similarly, the number of atoms created at level $i$ from a level $j > i$ is $N_j A_{ji}$, and that from a level $j < i$ is $N_j B_{ji} J_{ij}(v)$. Finally, the total variation of the number density of atoms at level $i$ in a volume $d^3 v dV$ of the phase space, i.e. the distribution $F_i = N_i/(d^3 v dV)$, resulting from the radiative processes alone, is

$$\left(\frac{\delta F_i}{\delta t}\right)_{\text{rad.}} = \sum_{j<i} (n_j f_j B_{ji} J_{ij}(v) - n_i f_i A_{ij}) + \sum_{j>i} (n_j f_j A_{ji} - n_i f_i B_{ji} J_{ij}(v)). \quad (9)$$

For inelastic collisions, we proceeded in an equivalent way by introducing the collisional excitation and deexcitation coefficients $C_{ij}^4$. We have:

$$\left(\frac{\delta F_i}{\delta t}\right)_{\text{inel.}} = \sum_{j \neq i} (n_j f_j C_{ji} - n_i f_i C_{ij}). \quad (10)$$

Velocity-changing collisions also destroy and create particles in phase space. The rate of these collisions is quantified by the parameter $Q_{v,i}$ for each atomic level (see Oxenius 1986, HOSI, and Paletou & Peymirat 2021). We thus write:

$$\left(\frac{\delta F_i}{\delta t}\right)_{\text{v.c.c.}} = n_i Q_{v,i} [f^M - f_i], \quad (11)$$

where $f^M$ is the Maxwell-Boltzmann velocity distribution defined as

$$f^M(v) = \frac{1}{v_{\text{th}}^3 \pi^{3/2}} e^{-v \cdot v / v_{\text{th}}^2}, \quad (12)$$

with $v_{th}$ being the most probable velocity. The Boltzmann equation is finally written, for level $i$, as

$$\sum_{j<i} (n_j f_j B_{ji} J_{ij}(v) - n_i f_i A_{ij}) + \sum_{j>i} (n_j f_j A_{ji} - n_i f_i B_{ji} J_{ij}(v)) + \sum_{j \neq i} (n_j f_j C_{ji} - n_i f_i C_{ij}) + n_i Q_{v,i} [f^M - f_i] = 0. \quad (13)$$

As we can see from Eq. (13), there is a coupling between each of the VDFs $f_i$ in the Boltzmann equations. This is precisely what the standard non-LTE radiative transfer cannot deal with. We also note that the Boltzmann equations are formally kinetic equilibrium equations (KEE). Integrating over all velocities<fnref n="5" />, we have the equilibrium between all the processes that populate and depopulate the atomic levels in a given space of volume $dV$. Hereafter, we refer to the following as integrated (and more usual) equations of the kinetic equilibrium or IKEE:

$$\sum_{j<i} (n_j B_{ji} \mathcal{J}_{ji} - n_i A_{ij}) + \sum_{j>i} (n_j A_{ji} - n_i B_{ji} \mathcal{J}_{ij}) + \sum_{j \neq i} (n_j C_{ji} - n_i C_{ij}) = 0, \quad (14)$$

where $\mathcal{J}_{ij}$ is the scattering integral, completely equivalent to the standard definition of $\bar{J}_{ij}$, and defined as

$$\mathcal{J}_{ij} = \int_{v} f_i(v) J_{ij}(v) d^3 v. \quad (15)$$

## 2.3. Atomic emission and absorption profiles of a three-bound-level atom

In what follows, we used a model atom with three bound levels. Using phenomenological arguments, HOSI obtained the explicit form of the atomic profiles required to calculate the observer's frame profiles (see Eqs. (4) and (5)) and the partial scattering integral $J_{ij}(v)$ (see Eq. (8)). For the sake of clarity, we briefly recall their main arguments. For more details, we refer the reader to HOSI and HOSII.

First of all, it is important to come back to the notion of natural population. For an ensemble of atoms of the same species, we are interested in a radiative transition from a level $b$ to a level $a$. Moreover, there is a given radiation field permeating the atmosphere, i.e. a given quantity of photons whose characteristics (frequencies and direction of propagation) are well known. We consider a particular physical process populating the level $b$. If the number of atoms created in atomic level $b$ is independent of the number of photons already present, then the $b \to a$ radiative transition will result in a new distribution of photons that is completely independent of the pre-existing photon distribution in the atmosphere. Thus, we can say that the physical process populating level $b$ leads to its natural population. This is the case for spontaneous emission, inelastic collisions, and velocity-changing collisions. However, radiative absorption depends on the radiation field and so does not naturally populate the atomic levels.

An atomic profile characterises a specific radiative transition; say, $b \to a$. To find an expression for it, we need to study all the transitions and series of transitions that populate level $b$ and, ultimately, lead to the $b \to a$ radiative transition. Let $c \to b$ be a transition populating level $b$:

(i) If $c$ naturally populates level $b$, then the resulting distribution of $b \to a$ will be independent of the resulting distribution of $c \to b$. Hence, we need to consider the physical phenomenon involving a single photon, namely the $b^* \to a$ radiative transition (where the notation $l^*$ means that level $l$ is naturally populated). To describe it, we use the generalised redistribution function $r_{ba}(\xi_{ba})$ related to

<footnote n="4">⁴ The $C_{ij}$ depend on the electronic density, which should, in future works, also be determined self-consistently with the kinetic equations.</footnote>

<footnote n="5">⁵ By integrating over all the velocities, we will be missing all the potential effects that may originate from the velocity-changing collisional process.</footnote>
the probability that a photon of frequency between $\xi_{ba}$ and $\xi_{ba} + d\xi_{ba}$ is involved in the $b \to a$ transition.

(ii) If $c$ does not naturally populate level $b$ but is itself naturally populated, we need to consider a two-photon phenomenon ($c^* \to b \to a$). To handle it, we define the conditional probability $j_{cba}(\xi_{ba})$; i.e. the probability that a photon of frequency between $\xi_{ba}$ and $\xi_{ba} + d\xi_{ba}$ is involved in the $b \to a$ transition, knowing that a photon of frequency between $\xi_{cb}$ and $\xi_{cb} + d\xi_{cb}$ was previously involved in the $c^* \to b$ transition. This quantity is defined as (see Eq. (4.4) of HOSI):

$$j_{cba}(\xi_{ba}) = \frac{\int I_{cb} r_{cba}(\xi_{cb}, \xi_{ba}) d\xi_{cb}}{\int I_{cb} r_{cba}(\xi_{cb}) d\xi_{cb}}. \qquad (16)$$

In this equation, $r_{cba}(\xi_{cb}, \xi_{ba})$ is the generalised redistribution function describing the probability that two photons whose frequencies are, respectively, between $\xi_{cb}, \xi_{cb} + d\xi_{cb}$ and $\xi_{ba}, \xi_{ba} + d\xi_{ba}$, are engaged in the $c^* \to b \to a$ transition (see Sect. 3 of HOSI). To understand the meaning of $j_{cba}(\xi_{ba})$, we need to interpret each quantity probabilistically. The specific intensity $I_{cb}$ can be related to the so-called photon number density (see e.g. Hubeny & Mihalas 2014). Thus, we can see the denominator term of Eq. (16) as the total number of photons involved in the $c \to b$ transition and the numerator term as the fraction of photons created by the $b \to a$ transition that actually comes from the $c \to b$ transition.

(iii) If $c$ does not naturally populate the $b$ level and is naturally populated by a third transition $d^* \to c$, we need to consider a phenomenon involving three photons in the series of transitions $d^* \to c \to b \to a$. To describe such a phenomenon, we should use the conditional probability $j_{cba}(\xi_{ba})$ defined analogously to Eq. (16).

(iv) It is clear that the $d$ level itself can be populated non-naturally by a given $e \to d$ transition, and similarly for $e$ with an $f \to e$ transition, and so on. Finally, when a photon is involved in the $b \to a$ transition, there will be a given probability that it originally comes from one of the series of transitions previously described. The atomic profile is the probability density describing all these possibilities simultaneously. It is defined as the average of all the series of possible transitions weighted by their occurrence probabilities.

For atoms with three levels, we therefore need to determine three absorption profiles $-\alpha_{12}, \alpha_{13}$, and $\alpha_{23}$ – as well as three emission profiles: $\beta_{21}, \beta_{31}$, and $\beta_{32}$. The $\alpha_{12}$ and $\alpha_{13}$ profiles are the simplest to determine. Since stimulated emission has been neglected, the fundamental level is naturally populated. Consequently, the probability of occurrence of such a phenomenon is certain; only case (i) needs to be taken into account. We therefore have

$$\alpha_{12} = r_{12}(\xi_{12}) \qquad (17)$$

and

$$\alpha_{13} = r_{13}(\xi_{13}). \qquad (18)$$

For $\beta_{21}$ and $\alpha_{23}$, we need to take into account all the series of transitions leading to level 2. It can be naturally populated, and we denote $prob(\to 2^*)$ as the relevant probability following HOSI and HOSII notations. This is case (i). It can also be populated non-naturally from the fundamental level, and we express

$prob(1^* \to 2)$ as its probability. This is case (ii). Combining these two phenomena, we have (see Eqs. (4.3) and (4.6) of HOSI)

$$\alpha_{23} = prob(\to 2^*)r_{23} + prob(1^* \to 2)j_{123} \qquad (19)$$

and

$$\beta_{21} = prob(\to 2^*)r_{21} + prob(1^* \to 2)j_{121}. \qquad (20)$$

For $\beta_{32}$ and $\beta_{31}$, we are interested in the processes populating level 3. We had to take into account the processes which do not naturally populate level 3 from level 2 and therefore treat the cascade $1^* \to 2 \to 3 \to 1$ or $1^* \to 2 \to 3 \to 2$ (case (iii)). We then have (see Eqs. (4.7) and (4.8) of HOSI)

$$\beta_{31} = prob(\to 3^*)r_{31} + prob(1^* \to 3)j_{131} + prob(\to 2^* \to 3)j_{231} + prob(1^* \to 2 \to 3)j_{1231} \qquad (21)$$

and

$$\beta_{32} = prob(\to 3^*)r_{32} + prob(1^* \to 3)j_{132} + prob(\to 2^* \to 3)j_{232} + prob(1^* \to 2 \to 3)j_{1232}. \qquad (22)$$

For more details on the expressions of ‘prob’ probabilities, we refer the reader to HOSI or Oxenius (1986). Knowledge of these quantities will not be explicitly necessary later on. In this work, we only studied the simplified case of an atom with three infinitely sharp levels. In Sect. 3, we will see how this simplifies the expression of the general atomic profiles given above and, consequently, the whole problem.

### 3. A three-level atom with infinitely sharp levels

With the exception of the fundamental level, atomic levels are generally broadened in energy. In the theoretical developments of HOS presented above, no assumptions on the atomic model were made. Considering the same physical processes, the assumption of infinitely sharp levels does not degrade the problem, although it does simplify it. However, we warn the reader that this hypothesis is part of a heuristic approach and that the realistic study of an astrophysical atmosphere does not leave us the choice of the atomic model, it imposes it. In fact, it is the physical conditions in which the atoms are immersed and the intrinsic nature of the atoms that dictate the width of the atomic levels (broadened or infinitely sharp). For example, if the conditions under investigation lead to collisional broadening, these physical conditions impose an atomic model. Similarly, if we assume an atom such that levels 4 and 5 are naturally broadened, but levels 2 and 3 are not (e.g. metastable levels), this atomic model must be adopted for a realistic study. Hereafter, we restrict ourselves to an atmosphere composed of atoms with three bound and infinitely sharp levels, neglecting stimulated emission and considering that the VDF associated with the fundamental level is Maxwellian ($f_1 \equiv f^M$).

#### 3.1. Atomic profiles

To calculate the atomic profiles of a three-level atom, we distinguish five phenomena that ultimately lead to the radiative transition under investigation: (i) absorption or emission of a single photon in a given $b \to a$ transition ($r_{ab}$ terms); (ii) resonance scattering, i.e. absorption and re-emission of a photon in a single line ($j_{aba}$ terms); (iii) Raman scattering, i.e. absorption and re-emission of a photon in two distinct lines ($j_{cba}$ terms);
(iv) two-photon absorption processes ($j_{123}$ term) and (v) three-photon processes ($j_{dcba}$ terms). Since our levels are assumed to be infinitely sharp, the probability that a photon involved in a given $b \to a$ transition has a frequency $\xi_{ba} \neq v_{0,ba}$ is zero. We then write

$$r_{ba}(\xi_{ba}) = \delta(\xi_{ba} - v_{0,ba}), \quad (23)$$

where $\delta$ is the Dirac distribution. Resonance scattering was studied by Hummer (1962), which obtained analytical expressions for the redistribution functions $r_{aba}(\xi_{ab}, \xi_{ba})$. Later, Hubený (1982) obtained general expressions for all two-photon redistribution functions $r_{cba}(\xi_{cb}, \xi_{ba})$. For atoms without natural broadening, we have, for the Raman or resonance scattering processes (see Eq. (10.135) of Hubený & Mihalas 2014),

$$r_{cba}(\xi_{cb}, \xi_{ba}) = \delta(\xi_{cb} - v_{0,cb})\delta(\xi_{ba} - v_{0,ba}) = r_{cb}(\xi_{cb})r_{ba}(\xi_{ba}). \quad (24)$$

The conditional probabilities $j_{cba}(\xi_{ba})$ are then written as follows:

$$j_{cba}(\xi_{ba}) = \frac{\int I_{cb} r_{cb}(\xi_{cb}) r_{ba}(\xi_{ba}) d\xi_{cb}}{\int I_{cb} r_{cb}(\xi_{cb}) d\xi_{cb}} = r_{ba}(\xi_{ba}). \quad (25)$$

In this study, we did not consider the two-photon absorption process. As discussed in Hubený & Mihalas (2014, pp. 322–323), this process is negligible in most astrophysical atmospheres. We also neglected all three-photon processes for which redistribution functions were obtained by Hubený & Oxenius (1987) from quantum mechanical calculations. For $\alpha_{23}$, for example, the term proportional to $j_{123}$ is assumed to be equal to zero, and therefore $\alpha_{23} = \text{prob}(\to 2^*)r_{23}$. We know that atomic profiles and generalised redistribution functions are normalised to unity because they are probability densities. Hence, to guarantee this condition, we renormalised $\alpha_{23}$, which gives $\alpha_{23} = r_{23}$. We proceeded in the same way when three-photon processes were neglected.

Finally, for a three-level atom with infinitely sharp levels, each of the six atomic profiles established in Eqs. (17)–(22) can then be written simply, $\forall j > i$, after renormalisation, as

$$\alpha_{ij}(\xi_{ij}) = r_{ij}(\xi_{ij}) = \delta(\xi_{ij} - v_{0,ij}) \quad (26)$$

and

$$\beta_{ji}(\xi_{ji}) = r_{ji}(\xi_{ji}) = \delta(\xi_{ji} - v_{0,ji}). \quad (27)$$

## 3.2. The multi-level atom problem

Our new problem therefore involves five coupled distributions: three specific intensities, associated with each of the radiatively allowed transitions, and two velocity distributions associated with the two excited levels. We must now jointly solve all the RTEs (see Eq. (1)), together with the IKEE (see Eq. (14)) and the Boltzmann equations (see Eq. (13)). The set of IKEEs is redundant. So, for atoms with three levels, we chose to replace one of the three equations, namely the one associated with the fundamental density level $n_1$, by the conservation equation $n_1 + n_2 + n_3 = n_{tot}$. Thus, solving the IKEE is equivalent to solving the system:

$$\begin{pmatrix} 1 & 1 & 1 \\ E_{21} & D_2 & E_{23} \\ E_{31} & E_{32} & D_3 \end{pmatrix} \begin{pmatrix} n_1 \\ n_2 \\ n_3 \end{pmatrix} = \begin{pmatrix} n_{tot} \\ 0 \\ 0 \end{pmatrix}, \quad (28)$$

with

$$E_{ij} = A_{ji} + B_{ji} \mathcal{J}_{ji} + C_{ji} \quad (29)$$

and

$$D_i = - \sum_{j \neq i} (A_{ij} + B_{ij} \mathcal{J}_{ij} + C_{ij}). \quad (30)$$

In the above equations, we (numerically) assume that $A_{ij} = 0$ for $i < j$ and $B_{ij} = 0$ for $i > j$.

Solving Boltzmann equations means solving a set of three equations, one of which is trivial ($f_1 = f^M$). We then have

$$\begin{pmatrix} 1 & 0 & 0 \\ e_{21} & d_2 & e_{23} \\ e_{31} & e_{32} & d_3 \end{pmatrix} \begin{pmatrix} f_1(\mathbf{v}) \\ f_2(\mathbf{v}) \\ f_3(\mathbf{v}) \end{pmatrix} = \begin{pmatrix} 1 \\ -n_2 Q_{V,2} \\ -n_3 Q_{V,3} \end{pmatrix} f^M(\mathbf{v}), \quad (31)$$

with

$$e_{ij} = n_j(A_{ji} + B_{ji} J_{ji}(\mathbf{v}) + C_{ji}) \quad (32)$$

and

$$d_i = -n_i \left( Q_{V,i} + \sum_{j \neq i} (A_{ij} + B_{ij} J_{ij}(\mathbf{v}) + C_{ij}) \right). \quad (33)$$

The choice of the atomic and scattering models does not directly affect the form of these equations. However, their ingredients do depend explicitly on them: the absorption and emission profiles in the observer's frame $\varphi_{ij}$ and $\psi_{ji}$ and the partial scattering integrals $J_{ij}(\mathbf{v})$, given, respectively, by Eqs. (4), (5), and (8). For infinitely sharp levels, the emission and absorption profiles take the form (see Appendix A)

$$\psi_{ji}(x_{ji}) = \frac{1}{2\Delta_{ji}} \int_{|x_{ji}|}^{\infty} f_j(u) u du \quad (34)$$

and

$$\varphi_{ij}(x_{ij}) = \frac{1}{2\Delta_{ij}} \int_{|x_{ij}|}^{\infty} f_i(u) u du, \quad (35)$$

together with

$$J_{ij}(\mathbf{u}) = \oint \frac{d\Omega_{ij}}{4\pi} \int \delta(x_{ij} - \mathbf{u} \cdot \mathbf{\Omega}_{ij}) I_{ij} dx_{ij}. \quad (36)$$

In the equations above, we introduce the reduced frequency $x_{ij} = (v_{ij} - v_{0,ij})/\Delta_{ij}$; the normalised velocity $\mathbf{u} = \mathbf{v}/v_{th}^6$; and the Doppler width $\Delta_{ij} = (v_{0,ij}/c)v_{th}$. We also note that $f_i(\mathbf{u})$ is the modulus of the VDF $f_i(\mathbf{u})$, i.e. its integration over all the directions $\mathbf{\Omega}_u$ of the atoms. Finally, it is easy to show that $\varphi_{12}$ and $\varphi_{13}$ are given by the mere Doppler profile,

$$\varphi_{ij}(x_{ij}) = \frac{1}{\Delta_{ij}} \times \frac{1}{\sqrt{\pi}} e^{-x_{ij}^2}, \quad (37)$$

because $f_1 = f^M$. A priori, all the other profiles will be different, and only a self-consistent resolution of the RTE, IKEE, and Boltzmann equations will make it possible to obtain an expression that is physically consistent with the FNLTE approach to radiative transfer.

<footnote n="6">⁶ We now write $f^M(\mathbf{u}) = \pi^{-3/2} e^{-u \cdot u}$. The factor $v_{th}^3$ disappears, in the integrals, with the change of variable $\mathbf{v} \to \mathbf{u}$.</footnote>
## 4. Numerical strategy

In this section, we present our numerical strategy for solving the FNLTE problem for three infinitely sharp bound-level atoms. We first initialised the iterative process to the multi-level CRD solution with the MALI method of Rybicki & Hummer (1991, hereafter the MALI-CRD method; see also Paletou & Lagache 2025). Then, we successively updated each of the radiative quantities using the following very simple strategy:

0: Initialize with a MALI-CRD solution, which gives the population densities $n_1, n_2$, and $n_3$ and therefore the CRD source functions.

1: For each of the transitions, calculate the emission and absorption coefficients, $\eta_{ji}, \chi_{ij}$, using Eqs. (2) and (3).

2: Calculate the source function $S_{ij}$:

$$S_{ij}(x_{ij}) = \frac{\eta_{ji}(x_{ij})}{\chi_{ij}(x_{ij})} \qquad (38)$$

3: Knowing the source function at all optical depths, compute photon distributions $I_{ij}$ for each of the transitions.

4: Compute the partial scattering integrals, $J_{ij}(u)$, given by Eq. (36), using the current VDF (a Maxwellian is adopted for the first iteration; see also Sect. 4 of PSP23).

5: Compute scattering integrals, $J_{ij}$, using Eq. (15).

6: Update all populations, $n_i$, by solving the IKEE system in Eq. (28).

7: Update the velocity distributions, $f_i$, by solving the system in Eq. (31) at each optical depth and atomic velocity $u$.

8: Update $\phi_{ij}$ and $\psi_{ij}$ profiles for each of the radiatively allowed transitions using Eqs. (35) and (34), and return to Step 1.

The output of this iterative process is a self-consistent determination of all the photon and massive particle distributions.

## 5. Validation

Since the multi-level FNLTE radiative transfer problem was never fully addressed before, we proceeded with the validation of our calculations through several indispensable verification steps. Furthermore, since we considered the restricted case of a three-level atom with infinitely sharp levels and neglected velocity-changing collisions (see below), the numerical solutions presented here essentially validate our numerical scheme and do not represent the real spectral lines of hydrogen formed in a stellar atmosphere.

All the solutions presented in this sections are calculated with a common set of parameters taken from Avrett (1968, Sect. 9.A; see also the code shared by Paletou & Lagache 2025)<fnref n="7" />. We considered an isothermal atmosphere composed of a three-bound-level hydrogen atom in a 1D, semi-infinite, plane-parallel geometry having a maximum total optical depth of $\tau_{max} = 10^{14}$ for the transition $1 \leftrightarrow 2$. We neglected the effects of velocity-changing collisions, i.e. $Q_V = 0$ (however, $Q_V \neq 0$ adds no numerical difficulties). We sample the optical depth scale with a logarithmic grid using four points per decade extending from 0 to $\tau_{max}$, with an initial step of $10^{-3}$. Integration over $\mu$ is performed using a six-point Gauss-Legendre quadrature. Integration over the azimuths (required to calculate $J_{ij}(u)$ – see PSP23) was performed using a ten-point rectangular quadrature. Finally, grids of reduced frequency, $x$, and normalised velocity, $u$, are identical and extend from 0 to $x_{max} = u_{max} = 4$ with a step of

<figure><img src="image_3.png" /><figcaption>Fig. 1. Frequency variation of normalised source function $S_{12}$ of $1 \leftrightarrow 2$ line at different optical depths $\tau$ for a two-level atom. Horizontal dashed black lines show CRD solutions for $\tau = 0, 1, 10, 100, 10^3, 10^4$. We compare FNLTE results (coloured lines) with Hummer’s PRD results (black open circle). CRD and PRD solutions were computed using the ALI and FBF methods, respectively. Clearly, we accurately recover the results presented by Hummer (1969) and PSP23.</figcaption></figure>

0.1. Velocity and frequency integrations were performed with trapezoidal quadrature. Each solution computed with our FNLTE code was obtained after 75 iterations, and the process is always initialised by the CRD solution, obtained after 100 iterations of the MALI-CRD method. Finally, for quantitative comparison purposes, we introduce the following relative error term, defined for each optical depth and frequency,

$$Err_{ij}^{rel} = \left| \frac{S_{ij}^{\text{FNLTE}} - S_{ij}^{\text{ref}}}{S_{ij}^{\text{ref}}} \right|, \qquad (39)$$

between our FNLTE source functions $S_{ij}$ and reference solutions described hereafter.

### 5.1. Benchmarking against a two-level atom in PRD

The first step in validating our multi-level code is to reduce it to a two-level atom problem with coherent scattering in the atom’s frame. PSP23 showed that when $Q_V = 0$ it is equivalent to Hummer’s results for PRD $R_{I-A}$ (Hummer 1962). Fig. 1 displays the source functions for $1 \leftrightarrow 2$ line, obtained with our FNLTE method (coloured lines) and compared to the reference $R_{I-A}$ PRD solution (open black circles) calculated using the FBF method (Paletou & Auer 1995). The mean of the relative error between these solutions is $\approx 0.31\%$. We are also able to reproduce the results obtained by PSP23 for the $f_2$ distribution function.

### 5.2. Benchmarking against multi-level CRD

Next, by forcing, in our new three-level FNLTE code, that $f_1 = f_2 = f_3 = f^M$, we recover the well-known results of MALI-CRD for Avrett’s hydrogen case (Avrett 1968, see also Paletou & Léger 2007). In Fig. 2, we show the optical depth variations of the

<footnote n="7">⁷ In particular, the temperature and collision rates $C_{ji}$ are fixed and independent of the optical depth.</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Optical depth variations of normalised source functions for each line of a three-level hydrogen atom in CRD. Our three-level FNLTE code (coloured lines) accurately reproduces reference solutions (open black circle) computed with MALI-CRD.</figcaption></figure>

source function calculated for each of the three lines using our multi-level FNLTE code and the MALI-CRD method. Mean relative errors of about 0.12% for 1 ↔ 2, 0.33% for 1 ↔ 3, and 0.32% for 2 ↔ 3 lines are reported, and they more likely come from completely different numerical schemes.

### 5.3. Benchmarking against XRD

On the basis of the theoretical work of HOSI & HOSII, the FNLTE approach of radiative transfer presented in Sect. 2 appears, at present, as the most complete approach (however, the QED approach of Bommier 2016 goes beyond the semi-classical picture of HOS formalism used here). In the past, indeed, numerous attempts have been made to incorporate it in more advanced radiative models, but with several simplifying assumptions. As a result of this, none of them could solve all the kinetic equations self-consistently with the RTE. The case where only resonance scattering is taken into account, namely the case of standard PRD, was studied by Hubený (1985). When Raman scattering is taken into account, we speak of 'cross-redistribution' (this term was used by Milkey et al. 1975 to describe this scattering process) or XRD (e.g. Hubený & Lites 1995; Sampoorna et al. 2013; Uitenbroek 1989). In all cases, these approaches were incomplete, with persisting approximations: (i) all absorption profiles and (ii) scattering integrals $J_{ij}$ were calculated using Maxwellian VDFs, (iii) $d_i \approx n_i(D_i - Q_{V,i})^8$ in the Boltzmann equations, and (iv) there is no coupling between the VDFs; in other words, to calculate $f_i$, we assumed that $f_j = f^M \forall j \neq i$. Considering these assumptions, the Boltzmann equations presented in Eq. (31) reduce to (assuming $Q_V = 0$ again), for all atomic excited levels,

$$f_2 \approx \frac{f^M}{n_2 P_2} \left\{ [n_3(A_{32} + C_{32}) + n_1 C_{12}] + n_1 B_{12} J_{12}(\mathbf{u}) \right\} \quad (40)$$

<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Optical depth variations of normalised source functions for the 1 ↔ 3 line of a three-level hydrogen atom in the frame of XRD approximation. We compare FNLTE results (coloured lines) for various frequencies $x = 0, 1, 2, 3, 4$ with XRD reference solutions (open black circle) computed using the MALI-XRD method of Sampoorna et al. (2013). The dashed black line shows the CRD solution computed using the MALI-CRD method.</figcaption></figure>

and

$$f_3 \approx \frac{f^M}{n_3 P_3} \left\{ [n_1 C_{13} + n_2 C_{23}] + n_1 B_{13} J_{13}(\mathbf{u}) + n_2 B_{23} J_{23}(\mathbf{u}) \right\}, \quad (41)$$

where $P_i$ are defined as

$$P_i = \sum_{j \neq i} (A_{ij} + B_{ij} J_{ij} + C_{ij}). \quad (42)$$

From these analytical expressions, we can deduce a general expression for the emission profiles (see Appendix B for more details):

$$\psi_{ji} = \frac{\varphi_{ij}^{*M}}{n_j P_j} \left\{ n_i B_{ij} \bar{J}_{iji}^M + \sum_{\substack{l \neq i,j \\ l < j}} n_l B_{lj} \bar{P}_{lji}^M \right. \\ \left. + \left[ n_i C_{ij} + \sum_{k \neq i,j} n_k (A_{kj} + C_{kj}) \right] \right\}, \quad (43)$$

where the notation $X^M$ means that quantities were calculated using a Maxwellian VDF, in particular for the generalised redistribution functions, $R_{ijk}$, contained in the definition of the diffusion integrals $\bar{J}_{iji}$ and $\bar{P}_{lji}$. This last equation is completely equivalent to Eq. (11) of Sampoorna et al. (2013) obtained in XRD; i.e. considering resonance and Raman scattering effects. Whether one works directly with the emission profiles (XRD) or solves the uncoupled kinetic equations (namely Eqs. (40) and (41)) self-consistently with the transfer equations, respective results should be identical. We show that, within the framework of the approximations presented above, there is an equivalence between our FNLTE formalism and the XRD formalism used by Sampoorna et al. (2013), whose numerical implementation is inspired by the MALI-PRD method developed by Paletou (1995).

The source functions obtained with these two methods are shown in Figs. 3 and 4, respectively, for the 1 ↔ 3 and 2 ↔ 3

<footnote n="8">⁸ Without v.c.c., we have $n_i D_i = \int d_i(\mathbf{u}) d^3\mathbf{u}$. Assuming $d_i \approx n_i D_i$ is equivalent to assuming that $\int d_i(\mathbf{u}) f(\mathbf{u}) d^3\mathbf{u} \approx \int d_i(\mathbf{u}) d^3\mathbf{u} \times \int f(\mathbf{u}) d^3\mathbf{u}$.</footnote>
<figure><img src="image_6.png" /><figcaption>Fig. 4. Same as Fig. 3, but for 2 ↔ 3 line.</figcaption></figure>

lines. For the 1 ↔ 3 line, we measure a mean relative error between our degraded FNLTE code solution and the XRD solution of about 0.23%. For 2 ↔ 3, it is about 0.37%. We chose not to show the source function associated with the 1 ↔ 2 line, as it differs only slightly from the solution shown in Fig. 1 (a relative maximum difference of 4.6% with respect to the two-level atom solution). However, XRD results are very well reproduced for this line with our FNLTE code, with a mean relative error between these two solutions of about 0.23%.

In Fig. 4, at the surface and at high frequency, we observe an overshoot; i.e. $S_{23}/B_W > 1$. The possibility that we are facing a numerical problem has been explored, in particular because we explored regions where the optical depths become very tiny. However, this feature is reproduced identically by the two different numerical methods. Furthermore, as the numerical problems remain identical, no similar behaviour has been observed in any of the other cases studied (such as MALI-CRD and MALI-PRD; the latter is discussed in Appendix C). This characteristic seems to be specific to XRD, which suggests that the consistency of some of its assumptions would need to be discussed further. However, as we see in Sect. 6, this artefact disappears when one finally deals with the complete problem. The above-mentioned – and very satisfactory – tests now lead us with confidence towards new results, using the full description.

## 6. New results

Hereafter, we used the same model and parameters as the ones used for our previous tests; in particular, velocity-changing collisions continue to be neglected (although their treatment does not present any numerical additional difficulty). Firstly, the $S_{12}$ source function for the 1 ↔ 2 line changes very slightly compared to the results obtained by PSP23 with a two-level atom with a mean relative error between these solutions of the order of 1.92%. This can be explained by the fact that atomic levels are much less populated when their energy is higher. In other words, we have $n_1 \gg n_2 \gg n_3 \gg \dots$. Thus, the source function for the 2 ↔ 1 transition is largely dominated by radiative processes populating level 2 from level 1 and vice versa. The impact of multi-level modelling on this line is therefore marginal. On the other hand, the 1 ↔ 3 and 2 ↔ 3 lines are substantially affected by FNLTE effects. In Figs. 5 and 6, we show the $S_{13}$ and $S_{23}$ source functions obtained for these lines. Differences can be seen at the surface with XRD, comparing Figs. 5 and 6, respectively, with Figs. 3 and 4. We know that it is at the surface that non-LTE

<figure><img src="image_7.png" /><figcaption>Fig. 5. Optical depth variations of normalised source functions for 1 ↔ 3 line of a three-bound-level hydrogen atom. It shows FNLTE results (coloured lines) obtained for various frequencies: $x = 0, 1, 2, 3, 4$. The dashed black line represents the MALI-CRD solution.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 6. Same as Fig. 5, but for 2 ↔ 3 line.</figcaption></figure>

effects are more pronounced. It is therefore not surprising that any change in the approach, describing all deviations from LTE, modifies the characteristics of the atmosphere in this region. As we explain in the previous section, the overshoot observed using XRD for the $S_{23}$ source function now disappears. This is because the frequency variations of the source function are mainly due to the ratio between the emission and the absorption profile:

$$ \rho_{ij}(x) = \frac{\psi_{ji}(x)}{\varphi_{ij}(x)}. \qquad (44) $$

With XRD, we have $\rho_{23} \gg 1$ precisely in this region close to the surface and far from the line centre where the overshoot is observed, as can be seen from the left panel of Fig. 7. According to XRD assumptions, all absorption profiles are calculated using a Maxwellian. However, we see in Sects. 2 and 3 that the absorption profile $\varphi_{23}$ should also be determined in a self-consistent way with all the other radiative quantities; i.e. calculated with $f_2 \neq f^M$. It therefore seems contradictory to use $f_2$ to determine the emission profile, $\psi_{21}$, but ignore this fact for the absorption profile, $\varphi_{23}$. This suggests that XRD’s assumption about absorption profiles is largely responsible for the overshoot. This characteristic is no longer present when $\varphi_{23}$ is self-consistently calculated, and, quantitatively, we can see that $\rho_{23}$ decreases significantly in the very same regions (see the right panel of Fig. 7).
<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Frequency variations of $\rho_{23}$ ratio for 2 ↔ 3 line of a three-bound-level hydrogen atom. We show the XRD (left plot) and FNLTE (right plot) results obtained for different optical depths: $\tau = 10^{-12}, 10^{-10}, 10^{-8}, 10^{-5}, 1$.</figcaption></figure>

A more detailed examination of the influence of each of the four XRD assumptions (see Sect. 5.3) taken separately shows that the assumption on absorption profiles has a major impact on the $S_{23}$ source function. Quantitatively, we measure a mean relative difference between the solution thus obtained and the FNLTE solution of the order of 65%. Qualitatively, and as discussed above, this assumption considerably increases the $S_{23}$ source function in the areas where the overshoot is observed in XRD. However, the amplitude of the overshoot is increased by a factor ≈2. This hypothesis alone does not fully explain this behaviour. Another XRD hypothesis must therefore have the opposite effect, by reducing the $S_{23}$ source function near the surface and at high frequency. This is precisely the qualitative effect of the hypothesis on the coupling between the VDFs in the Boltzmann equations. The inclusion of these two assumptions alone explains the appearance of this overshoot, with a maximum relative deviation between the solution thus obtained and the XRD solution of the order of ≈3%. The other assumptions have a much smaller impact, of the order of a few percent on average. However, what is more important is that the FNLTE approach allows us to get rid of all these assumptions and solves the XRD overshooting issue. It is important to note that the XRD overshooting problem seen in the present work (see Fig. 4) is also due to the use of a three-level model atom with infinitely sharp levels. More specifically, the overshoot in $S_{23}$ from the XRD hypothesis is a pathological feature of $R_I$ (Hummer 1962) and $P_I$ (Hubený 1982) and does not occur with $R_{II}$ and $P_{II}$, which include the natural broadening of energy levels. Indeed, the latter has been verified using the MALI-XRD method of Sampoorna et al. (2013). Moreover, in a more realistic case including more atomic levels, one would expect to observe a strong influence of the higher levels on level 3, modifying the behaviour near the surface.

Finally, it has been emphasised throughout this paper that the FNLTE approach for radiative modelling deals not only with the possible deviations of the radiation field and atomic populations from their equilibrium distributions, but also for the deviation of the VDFs from Maxwellian. Therefore, in addition to the usual source functions, we are now able to obtain all VDFs, $f_i \neq f^M$, of the excited levels. As $S_{21}$, we find that $f_2$ varies only slightly from the two-level results previously obtained by PSP23. Indeed,

<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Normalised atomic velocities variations of $f_3$ velocity distribution of the second excited level of a three-bound-level hydrogen atom. Results (coloured lines) were obtained for different optical depths: $\tau_{31} = 0.5, 1, 10, 10^2, 10^8, 5 \cdot 10^8, 10^{12}$. The dashed black line shows the Maxwellian limit $f_3(u) = f^M(u)$.</figcaption></figure>

we measure a mean relative deviation of 0.084% between these two solutions. However, the velocity distribution associated with level 3 is a new feature introduced by the FNLTE approach, as displayed in Fig. 8. The higher the atomic velocity modulus *u* is, the more $f_3$ deviates from the Maxwellian. In deeper layers, as $f_3 \rightarrow f^M$ the atmosphere tends to be thermalised. We also note that the variations in $f_3$ are very similar to those in $\rho_{23}$ (XRD), which is perfectly understandable given the explicit dependence of the profiles on the distribution functions (especially because $\psi_{32}$ is calculated using $f_3$).

## 7. Conclusions

Although they were thoroughly formulated more than forty years ago by Hubený et al. (1983a,b) (see also Oxenius 1986), FNLTE effects for a multi-level atom had never been computed before the present study. As a first step, we restricted ourselves to a simple three-bound-level model atom, based on parameters initially
provided by Avrett (1968), and infinitely sharp levels for several
benchmark purposes.

As our preliminary solutions to the multi-level FNLTE prob-
lem cannot be directly compared with other studies, we verified
our ability to reproduce several ‘degraded’ problems with accu-
racy. This also led us to point out some inconsistencies of the
quite advanced still XRD modelling.

The level of realism of the atomic and scattering models will
need to be improved in order to take into account the natural
broadening of excited atomic levels. However, the major con-
ceptual jump from two to three or more levels has already been
made, and the Sampoorna et al. (2024) step will be generalised
further.

In the longer term, we could also consider non-Maxwellian
free electrons in a self-consistent way as other massive parti-
cles. Cases where stimulated emission is not negligible could
also be explored, although this would involve major changes in
the calculation of atomic profiles.

Finally, we are aware that realism comes at a high computa-
tional price. For this purpose, we need to develop a more efficient
numerical scheme. However, another study of ours (Lagache
et al. 2025) shows that it is possible to adapt approximate opera-
tor methods (see e.g. Hubený 2003; Paletou & Auer 1995) to the
two infinitely-sharp-level FNLTE problem. Therefore, in the near
future, a robust and convergent method similar to these ones will
be adapted to our new multi-level, multi-distribution problem.

*Acknowledgements.* T. Lagache is kindly supported by a 'PhD booster' grant by the Toulouse Graduate School of Earth and Space Sciences (TESS; https://tess.omp.eu/). Authors thank Prof. Ivan Hubeny for very critical and insightful comments.

References

Avrett, E. H. 1968, in Resonance Lines in Astrophysics, eds. R. G. Athay, J. Mathis, & A. Skumanich (Boulder: National Center for Atmospheric Research), 27

Gray, M. 2012, Maser Sources in Astrophysics, 50 (Cambridge University Press)

Heinzel, P., Gouttebroze, P., & Vial, J.-C. 1987, A&A, 183, 351

Hubený, I. 1981, A&A, 100, 314

Hubený, I. 1982, JQSRT, 27, 593

Hubený, I. 1985, Bull. Astron. Inst. Czech., 36, 1

Hubený, I. 2003, Stellar Atmos. Model., 288, 17

Hubený, I., & Cooper, J. 1986, ApJ, 305, 852

Hubený, I., & Lites, B. W. 1995, ApJ, 455, 376

Hubený, I., & Mihalas, D. 2014, Theory of Stellar Atmospheres: An Introduction to Astrophysical Non-equilibrium Quantitative Spectroscopic Analysis, 26 (Princeton University Press)

Hubený, I., & Oxenius, J. 1987, JQSRT, 37, 65

Hubený, I., Oxenius, J., & Simonneau, E. 1983a, JQSRT, 29, 477

Hubený, I., Oxenius, J., & Simonneau, E. 1983b, JQSRT, 29, 495

Hummer, D. G. 1962, MNRAS, 125, 21

Hummer, D. G. 1969, MNRAS, 145, 95

Lagache, T., Sampoorna, M., & Paletou, F. 2025, RAS Tech. Instrum., rzaf023

Milkey, R. W., Shine, R. A., & Mihalas, D. 1975, ApJ, 199, 718

Oxenius, J. 1965, JQSRT, 5, 771

Oxenius, J. 1986, Kinetic Theory of Particles and Photons: Theoretical Foundations of Non-LTE Plasma Spectroscopy (Berlin: Springer)

Paletou, F. 1995, A&A, 302, 587

Paletou, F., & Auer, L. H. 1995, A&A, 297, 771

Paletou, F., & Lagache, T. 2025, https://hal.science/hal-05051883

Paletou, F., & Léger, L. 2007, JQSRT, 103, 57

Paletou, F., & Peymirat, C. 2021, A&A, 649, A165

Paletou, F., Sampoorna, M., & Peymirat, C. 2023, A&A, 671, A93

Rybicki, G. B., & Hummer, D. G. 1991, A&A, 245, 171

Sampoorna, M., Nagendra, K. N., & Stenflo, J. O. 2013, ApJ, 770, 92

Sampoorna, M., Paletou, F., Bommier, V., & Lagache, T. 2024, A&A, 690, A213

Uitenbroek, H. 1989, A&A, 213, 360
## Appendix A: Derivation of observer's frame profiles

The absorption and emission profiles in the observer's frame are given by the Eqs. (4) and (5). For infinitely sharp levels, we have given the form of the atomic profiles in the Eqs. (26) and (27). Thus, with $j > i$, we write (for example for $\psi$ but the calculation is completely analogous for $\varphi$):

$$ \psi_{ji}(v_{ji}) = \oint \frac{d\Omega_{ji}}{4\pi} \int f_j(\mathbf{u})\delta(\xi_{ji} - v_{0,ij})d^3\mathbf{u}. \quad (A.1) $$

Using the definitions of the reduced frequency $x_{ji}$, the normalised velocity $u$ and the Doppler-Fizeau effect given in Eq. (6), we obtain:

$$ \psi_{ji} = \frac{1}{\Delta_{ji}} \oint \frac{d\Omega_{ji}}{4\pi} \int f_j(\mathbf{u}) \delta(x_{ji} - \mathbf{u} \cdot \mathbf{\Omega}_{ji}) u^2 du d\Omega_u. \quad (A.2) $$

Noticing that:

$$ \oint \frac{d\Omega_{ji}}{4\pi} \delta(x_{ji} - \mathbf{u} \cdot \mathbf{\Omega}_{ji}) = \frac{1}{2u} H(u - |x_{ji}|), \quad (A.3) $$

with $H$ the Heaviside function, then:

$$ \psi_{ji}(x_{ji}) = \frac{1}{2\Delta_{ji}} \int_{|x_{ji}|}^{\infty} f_j(u) u du. \quad (A.4) $$

## Appendix B: Derivation of a general expression for emission profiles in the XRD approximations

The aim here is to obtain a general expression for the emission profiles within the framework of the approximations of the so-called XRD approach of the non-LTE radiative transfer. These approximations are described in detail in Sect. 5.3. Before applying this set of assumptions, we propose a more general formulation based on the kinetic equations given in Eq. (13). For this purpose, we rewrite these equations as:

$$ n_i \Pi_i(\mathbf{u}) f_i = L_i, \quad (B.1) $$

with,

$$ \Pi_i(\mathbf{u}) = Q_{v,i} + \sum_{j \neq i} [\tilde{R}_{ij}(\mathbf{u}) + C_{ij}], \quad (B.2) $$

and, following HOSI, we introduce $L_i$, the density corresponding to the ensemble of processes populating naturally or not the atomic level $i$ in the phase space:

$$ L_i = n_i Q_{v,i} f^M + \sum_{j \neq i} n_j f_j [\tilde{R}_{ji}(\mathbf{u}) + C_{ji}]. \quad (B.3) $$

Analogously, $\Pi_i(\mathbf{u})$ is related to the density corresponding to the set of processes depopulating the level $i$ in phase space<fnref n="9" />. In the previous equations, $\tilde{R}_{ij}(\mathbf{u})$ is the radiative rate dependent on

the velocity associated with the transition $i \to j$. It is defined as:

$$ \tilde{R}_{ij}(\mathbf{u}) = \begin{cases} B_{ij} J_{ij}(\mathbf{u}) & \text{if } i < j, \\ A_{ij} & \text{if } i > j. \end{cases} \quad (B.4) $$

Integration over all velocities of this quantity, i.e.:

$$ R_{ab} = \int d^3 \mathbf{u} f_a(\mathbf{u}) \tilde{R}_{ab}(\mathbf{u}), \quad (B.5) $$

leads to the standard radiative rate which appear in the IKEE defined as:

$$ R_{ij} = \begin{cases} B_{ij} J_{ij} & \text{if } i < j, \\ A_{ij} & \text{if } i > j. \end{cases} \quad (B.6) $$

Also, with velocity-changing collisions, we define $P_i$ as:

$$ P_i = Q_{v,i} + \sum_{j \neq i} (R_{ij} + C_{ij}), \quad (B.7) $$

where stimulated emission is neglected. Then, we can write the $\beta_{ji}$ atomic emission profiles describing the $j \to i$ transition (this formula was adapted from HOSI; see also our Sect. 2.3) as:

$$ \begin{aligned} \beta_{ji} = & \text{prob}(\to j^*) r_{ji} + \sum_{l<j} \text{prob}(\to l^* \Rightarrow j) j_{lji} \\ & + \sum_{k<l<j} \text{prob}(\to k^* \Rightarrow l \Rightarrow j) j_{klji} \\ & + \sum_{m<k<l<j} \text{prob}(\to m^* \Rightarrow k \Rightarrow l \Rightarrow j) j_{mklji} + \dots \end{aligned} \quad (B.8) $$

Next, according to HOSI, we can write:

$$ \text{prob}(\to j^*) = \frac{1}{L_j} \left[ n_j f^M Q_{v,j} + \sum_{p \neq j} n_p f_p (A_{pj} + C_{pj}) \right]. \quad (B.9) $$

The term in the numerator relates to all processes that naturally populate level $j$. From the denominator, we retrieve the term $L_j$, the meaning of which we have already discussed. Also, the probability that level $j$ is non naturally populated from a level $l < j$ is<fnref n="10" />:

$$ \text{prob}(\to l^* \Rightarrow j) = \text{prob}(\to l^*) \times \frac{n_l f_l B_{lj} J_{lj}(\mathbf{u})}{L_j}. \quad (B.10) $$

We also have:

$$ \text{prob}(\to k^* \Rightarrow l \Rightarrow j) = \text{prob}(\to k^* \Rightarrow l) \times \frac{n_l f_l B_{lj} J_{lj}(\mathbf{u})}{L_j}, \quad (B.11) $$

and,

$$ \text{prob}(\to m^* \Rightarrow k \Rightarrow l \Rightarrow j) = \text{prob}(\to m^* \Rightarrow k \Rightarrow l) \times \frac{n_l f_l B_{lj} J_{lj}(\mathbf{u})}{L_j}. \quad (B.12) $$

In Sect. 3, we have neglected all processes involving three or more photons. In other words, we have implicitly made the assumption that $\text{prob}(\to k^* \Rightarrow l \Rightarrow j) = 0, \forall k, l, j$ with $k < l < j$. We must therefore have $n_l f_l B_{lj} J_{lj}(\mathbf{u}) = 0$ which is absurd as none

<footnote n="9">⁹ Velocity-changing collisions cannot populate a given level $i$ from a level $j$. However, they change the velocity of the atom inducing a transition from one state (position and velocity) to another in phase space. These collisions therefore populate a level by creating or destroying particles in phase space.</footnote>

<footnote n="10">¹⁰ For $l = 1, 2$, we get Eqs (7.4), (7.7) and (7.8) of HOSI.</footnote>
of these quantities can be zero or, $prob(\to k^* \Rightarrow l) = 0$, implying<fnref n="11" /> $prob(\to l^*) = 1$. Thus, to satisfy this last condition, we must assume that all levels $l$ such that $l < j$ are naturally populated. This assumption only applies to the calculation of a given emission profile. For example, at three levels, the calculation of $\beta_{31}$ and $\beta_{32}$ is made by assuming that level 2 is naturally populated (which is factually false because radiative absorption $1 \to 2$ does not lead to natural population of level 2). On the other hand, for the calculation of $\beta_{21}$, level 2 will not be considered to be naturally populated in any circumstances. With this clarification, we write the atomic emission profiles as:

$$ \beta_{ji} = \text{prob}(\to j^*) r_{ji} + \sum_{l<j} \text{prob}(\to l^* \Rightarrow j) j_{lji}. \quad (B.13) $$

Finally, the emission profile in the observer's frame $\psi_{ji}$ now writes, using Eq. (5) and neglecting velocity-changing collisions, as:

$$ \psi_{ji} = \oint \frac{d\Omega_{ji}}{4\pi} \int d^3\mathbf{u} \frac{f_j(\mathbf{u})}{L_j} \left\{ r_{ji} \sum_{p \neq j} n_p f_p(A_{pj} + C_{pj}) \right. \\ \left. + \sum_{l<j} \text{prob}(\to l^*) \times n_l f_l B_{lj} J_{lj}(\mathbf{u}) j_{lji} \right\}. \quad (B.14) $$

To continue, we use the Boltzmann equations i.e. Eq. (B.1) and our assumption on the natural population of atomic levels $l < j$, namely, $prob(\to l^*) = 1$; which is a consequence of the assumption on processes with three or more photons. We then have:

$$ \psi_{ji} = \oint \frac{d\Omega_{ji}}{4\pi} \int \frac{d^3\mathbf{u}}{n_j \Pi_j(\mathbf{u})} \left\{ r_{ji} \sum_{p \neq j} n_p f_p(A_{pj} + C_{pj}) \right. \\ \left. + \sum_{l<j} n_l f_l B_{lj} J_{lj}(\mathbf{u}) j_{lji} \right\}. \quad (B.15) $$

In order to retrieve usual non-LTE transfer approximations (including XRD), we use $\Pi_i \approx P_i$ (equivalent to assumption (iii) – see Sect. 5.3; see also HOSI for more details), and the above equation reduces to:

$$ \psi_{ji} = \frac{1}{n_j P_j} \oint \frac{d\Omega_{ji}}{4\pi} \int d^3\mathbf{u} \left\{ r_{ji} \sum_{p \neq j} n_p f_p(A_{pj} + C_{pj}) \right. \\ \left. + \sum_{l<j} n_l f_l B_{lj} J_{lj}(\mathbf{u}) j_{lji} \right\}. \quad (B.16) $$

We recall that $j_{lji} = \int I_{lj} r_{lji} d\xi_{lj} / J_{lj}(\mathbf{u})$ and the generalised scattering integrals denoted $\mathcal{R}_{lji}$ are defined as:

$$ \mathcal{R}_{lji}(x_{ij}) = \oint \frac{d\Omega_{li}}{4\pi} \int I_{lj} \frac{R_{lji}^l(x_{ij}, x_{lj})}{\varphi_{ij}^{M*}} dx_{lj}, \quad (B.17) $$

with $\varphi_{ij}^{k*}$ the absorption profile in the observer's reference frame calculated with the VDFs $f_k$ and $R_{lji}^l$ the generalised redistribution function integrated in angle and velocity (and calculated with $f_j$). Following HOSII, they are defined respectively as:

$$ \varphi_{ij}^{k*} = \oint \frac{d\Omega_{ji}}{4\pi} \int d^3\mathbf{u} f_k(\mathbf{u}) r_{ji}, \quad (B.18) $$

and,

$$ R_{lji}^l(x_{ij}, x_{lj}) = \oint \frac{d\Omega_{ji}}{4\pi} \int d^3\mathbf{u} f_l(\mathbf{u}) r_{lji}. \quad (B.19) $$

The emission profile $\psi_{ji}$ is then finally written as:

$$ \psi_{ji} = \frac{\varphi_{ij}^{M*}}{n_j P_j} \left\{ \sum_{l<j} n_l B_{lj} \mathcal{R}_{lji}^l + \sum_{p \neq j} n_p (A_{pj} + C_{pj}) \frac{\varphi_{ij}^{p*}}{\varphi_{ij}^{M*}} \right\}. \quad (B.20) $$

By distinguishing Raman and resonance scattering, and noting $\mathcal{R}_{lji}^l = \bar{P}_{lji}^l$ and $\mathcal{R}_{lji}^l = \bar{J}_{lji}^l$ respectively their generalised scattering integrals, we have:

$$ \psi_{ji} = \frac{\varphi_{ij}^{M*}}{n_j P_j} \left\{ n_i B_{ij} \bar{J}_{iji}^M + \sum_{\substack{l \neq i, j \\ l < j}} n_l B_{lj} \bar{P}_{lji}^l + \sum_{p \neq j} n_p (A_{pj} + C_{pj}) \frac{\varphi_{ij}^{p*}}{\varphi_{ij}^{M*}} \right\}. \quad (B.21) $$

Finally, we make the last assumption of the XRD approximation: to calculate $f_i$, we assume that $f_j = f^M \forall j \neq i$. Thus, the quantities defined in Eqs. (B.18) and (B.19) are calculated using a Maxwellian VDF and we simply obtain, after manipulating sums and indices:

$$ \psi_{ji} = \frac{\varphi_{ij}^{M*}}{n_j P_j} \left\{ n_i B_{ij} \bar{J}_{iji}^M + \sum_{\substack{l \neq i, j \\ l < j}} n_l B_{lj} \bar{P}_{lji}^M \right. \\ \left. + \left[ n_i C_{ij} + \sum_{k \neq i, j} n_k (A_{kj} + C_{kj}) \right] \right\}, \quad (B.22) $$

which is completely equivalent to Eq. (43).

## Appendix C: Benchmarking against standard PRD

We are fully aware that cross-redistribution (XRD) is a mechanism still rarely invoked in radiative modelling. As an extra validation test, we propose here to study the case of standard Hummer's PRD, which is absolutely central to explaining the observed properties of the $1 \leftrightarrow 2$ and $1 \leftrightarrow 3$ lines of the Sun (see Hubeny & Lites 1995; Heinzel et al. 1987). In the present case, we consider a hydrogen atom with three infinitely sharp levels for which only Hummer's PRD effects will be considered for the two resonance lines ($1 \leftrightarrow 2$ and $1 \leftrightarrow 3$ lines), and where the CRD approximation will be used for $2 \leftrightarrow 3$ line i.e. $\psi_{32} = \varphi_{23}$. For $1 \leftrightarrow 2$ and $1 \leftrightarrow 3$ lines, an emission profile is calculated without taking into account the effects of cross-redistribution. This is equivalent to considering the second term of Eq. (43) as zero, i.e. $\bar{P}_{lji}^M = 0$. For $1 \leftrightarrow 2$ line, this consideration has no impact and so $f_2$ will always be given by Eq. (40). On the other hand, for $1 \leftrightarrow 3$ line, we must now calculate $f_3$ with $J_{23}(\mathbf{u}) = \mathcal{J}_{23} = 0$. Note, moreover, that according to the IKEE, we have:

$$ n_1 C_{13} + n_2 C_{23} = n_3 P_3 - n_2 B_{23} \mathcal{J}_{23} - n_1 B_{13} \mathcal{J}_{13}, \quad (C.1) $$

and we finally write:

$$ f_3 = f^M + \frac{f^M}{n_3 P_3} \left[ n_1 B_{13} (J_{13}(\mathbf{u}) - \mathcal{J}_{13}) \right]. \quad (C.2) $$

<footnote n="11">¹¹ Because $prob(\to k^* \Rightarrow l) + prob(\to l^*) = 1$</footnote>
<figure><img src="image_11.png" /></figure>

**Fig. C.1.** Optical depth variations of normalised source functions for 1 ↔ 3 line of a three-level hydrogen atom. We compare FNLTE results (coloured lines) for various frequencies $x = 0, 1, 2, 3, 4$ with PRD reference solutions (black open circle) computed using the MALI-PRD method of Paletou (1995). The black dashed line shows the CRD solution computed using MALI-CRD method.

<figure><img src="image_12.png" /></figure>

**Fig. C.2.** Optical depth variations of normalised source functions for 2 ↔ 3 line of a three-level hydrogen atom. We compare FNLTE results (red line) with PRD reference solutions (black open circle) computed using the MALI-PRD method of Paletou (1995). The black dashed line shows the CRD solution computed using MALI-CRD method.

The combination of Eqs. (40), (C.2) then allows us to compare
the results given by our FNLTE multi-level code with the refer-
ence results calculated using the MALI-PRD method developed
by Paletou (1995). The results are shown in Fig. C.1 for the 1 ↔ 3
line and in Fig. C.2 for 2 ↔ 3 line. We do not show the results for
1 ↔ 2 line for the same reasons explained in Sect. 5.3. Finally,
we see that the expected results are very well reproduced, with a
mean relative error between the two solutions equal, for 1 ↔ 2,
1 ↔ 3 and 2 ↔ 3 lines respectively, to about 0.22%, 0.37% and
0.34%.