# An alternative iterative scheme to the full non-LTE transfer problem with coherent scattering in the atom's frame

T Lagache, M Sampoorna, F Paletou

► To cite this version:

T Lagache, M Sampoorna, F Paletou. An alternative iterative scheme to the full non-LTE transfer problem with coherent scattering in the atom's frame. RAS Techniques and Instruments, 2025, 4, pp.rzaf023. 10.1093/rasti/rzaf023. hal-05148611

HAL Id: hal-05148611

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

Submitted on 7 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# An alternative iterative scheme to the full non-LTE transfer problem with coherent scattering in the atom's frame

T. Lagache<sup>1,2*</sup> M. Sampoorna<sup>1,3</sup> and F. Paletou<sup>1,2</sup>

Accepted 2025 June 6. Received 2025 May 23; in original form 2025 January 14

## ABSTRACT

We propose an alternative scheme for the computation of the so-called 'full non-LTE' (non-local thermodynamic equilibrium) radiative transfer problem, assuming coherent scattering in the atom's frame. This generalized problem should explicitly deal with the coupling between atomic velocities and the photon paths, thereby requiring the implementation of new numerical strategies to solve it. Recently, Paletou et al., presented a numerical scheme, based on the Λ-iteration, to solve this problem; however it needs to be initialized using the standard non-LTE solution with complete redistribution to achieve convergence. Our new scheme is based on accelerated Λ-iteration (ALI), which is robust and insensitive to the choice of the initial guess solution. After bench-marking our new iterative scheme against the previously developed Λ-iteration method, we demonstrate its robustness by studying its convergence behaviour. This new iterative scheme has been coded in the *Julia* language, and its main characteristics are hereafter described with some details.

**Key words:** Radiative transfer – Line formation – Numerical methods.

## 1 INTRODUCTION

Although the non-local thermodynamic equilibrium (hereafter non-LTE) radiative transfer problem was first solved numerically during the 60's (see e.g. Hubený & Mihalas 2014), most of this endeavour, until recently, have avoided to consider possible departures of the velocity distribution functions (VDFs) of the massive particles constituting the astrophysical plasma, from Maxwellian. Paletou & Peymirat (2021, hereafter PP21) opened again this case, which was left aside in the astrophysical community since the mid-eighties (see references in PP21). As stated by Hubený, Oxenius & Simonneau (1983), '[...] both the line profile coefficients of an atom and those of a volume element of a gas can be determined only if one knows the distribution functions of photons and particles, that is, if one solves the entire problem of spectral line formation in a self-consistent manner'. Hereafter we shall therefore deal with the so-called 'full non-LTE' (FNLTE) problem, which aims at computing self-consistently such VDFs with the radiation field in a stellar atmosphere. Recent developments of the FNLTE approach to radiative transfer for the case of naturally broadened upper level atoms (Sampoorna et al. 2024) and for multilevel atoms (Lagache, Paletou & Sampoorna 2025) make it possible to describe in greater detail the dynamical state of massive particles in astrophysical atmospheres. Extending this formalism to the calculations of electrons distribution would also give access to a physical quantity that is critical for plasma processes. FNLTE is also relevant for the transfer of polarized light (see e.g. Landi Degl'Innocenti & Landolfi 2004, §13). In general,

precise knowledge of the VDF's would improve the coupling of the non-LTE radiation transfer problem with (magneto-)hydrodynamic simulations that could potentially affect any atmosphere that requires constraints on these VDF's.

The main difficulty of solving the standard non-LTE radiative transfer come from the coupling between the radiation field and matter properties described respectively by the specific intensity and the source function. To calculate the source function, we need to know the radiation field and vice versa. The radiative transfer is indeed a 'difficult problem' (see e.g. Rutily & Chevallier 2006). An iterative process such as the accelerated Λ-iteration (hereafter ALI) with a diagonal operator, represents a method of choice since the advent of Olson, Auer & Buchler (1986) and, opened-up a wide range of applications, from multidimensional problems (e.g. Auer & Paletou 1994), multi-level ones (Rybicki & Hummer 1991, 1992) to complete and partial redistribution on frequency (hereafter CRD and PRD) ones (Paletou 1995; Paletou & Auer 1995) or polarization (e.g. Sampoorna, Nagendra & Stenflo 2013, and references therein) etc... (see also Hubený & Mihalas 2014, for a more comprehensive list of applications). The numerical difficulties involved in solving the full non-LTE problem are further escalated, as the scattering integral now depends not only on the radiation field but also on the velocity of the radiating atom. In other words the VDF of the excited level of the radiating atom needs to be determined consistently with the radiation field. Thus the ALI based methods need to be suitably generalized to treat such problems.

In this paper, we propose a new ALI based iterative scheme which provides an alternative and a more robust solution method, as compared to the one first proposed by Paletou, Sampoorna &

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

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

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

<footnote>*E-mail: tristan.lagache@univ-tlse3.fr</footnote>
## 2 T. Lagache, M. Sampoorna and F. Paletou

Peymirat (2023, hereafter PP23). Section 2 explains this new scheme, for the case of coherent scattering in the atom's frame. Then Section 3 provides classic benchmark results, validating the new approach. Most important technical details are presented in Appendices A and B. Conclusions are given in Section 4.

## 2 VELOCITY-BY-VELOCITY METHOD

In this section we develop an ALI based velocity-by-velocity method to solve the problem at hand. The concerned basic equations are presented in detail in PP21 and PSP23. We refer the reader to these articles. Here we focus directly on the new numerical method of solution. We remark in passing that like in PP21 and PSP23, the intensity and the source function are normalized to the Planck function (specifically the Wien function, as we neglect stimulated emission).

### 2.1 The iterative problem

When we examine each of the quantities presented in previous works (PP21, PSP23; Sampoorna et al. 2024), the most critical quantity that needs to be calculated is the partial scattering integral $\tilde{J}_{12}(u, \tau)$. It is directly related to the radiation field, which allows us to implement an ALI-based method (see e.g. Hubený & Mihalas 2014, and references therein). Numerically, we compute the radiation field by performing successive formal solutions of the radiative transfer equation. This operation is expressed, using the so-called $\Lambda_{x\mu}$-operator<fnref n="1" /> as:

$$I_{x\mu} = \Lambda_{x\mu}[S_x], \quad (1)$$

with $S_x$ the isotropic source function.

By injecting equation (1) into the definition of $\tilde{J}_{12}(u, \tau)$ given in equations (6) and (1) of PSP23, we then have a problem expressed in term of a velocity-dependent quantity:

$$\tilde{J}_{12}(u, \tau) = \oint \frac{d\Omega_u}{4\pi} \oint \frac{d\Omega}{4\pi} \int_{-\infty}^{+\infty} \delta(x - \vec{u} \cdot \vec{\Omega}) \Lambda_{x\mu}[S_x] dx, \quad (2)$$

where $\vec{\Omega}_u$ and $\vec{\Omega}$ are respectively the direction of the velocity of the atoms and the direction of the light ray (note that here we use the conventional 'reduced frequency', $x$, as well as the normalized velocity of the atoms, $u$ which was introduced by PP21).

Our aim is therefore to solve the FNLTE radiative transfer problem working directly on velocities $u$, and avoiding the explicit use of standard redistribution functions (Hummer 1962). Here, we propose to write an iterative scheme working directly on $\tilde{J}_{12}(u, \tau)$. First, we replace the source function $S_x$ by its expression, where $\varphi(x)$ and $\psi(x)$ are respectively the absorption and emission profiles in the observer's frame:

$$S(x, \tau) = [\varepsilon + (1 - \varepsilon)\mathcal{J}_{12}] \frac{\psi(x, \tau)}{\varphi(x)}, \quad (3)$$

with $\mathcal{J}_{12}$ the velocity-average of $\tilde{J}_{12}$ defined as:

$$\mathcal{J}_{12} = 4 \int_{0}^{\infty} \tilde{J}_{12}(u, \tau) u^2 f^M(u) du. \quad (4)$$

Here, we define the 1D Maxwellian distribution $f^M(u)$ as:

$$f^M(u) = \frac{1}{\sqrt{\pi}} e^{-u^2}. \quad (5)$$

Using equation (5) of PSP23, we can rewrite equation (3) as:

$$S(x, \tau) = \frac{1}{1 + \zeta} \cdot \frac{2}{\varphi(x)} \int_{|x|}^{\infty} [\varepsilon + (1 - \varepsilon)\tilde{J}_{12}(u, \tau)] uf^M(u) du \\ + \frac{\zeta}{1 + \zeta} [\varepsilon + (1 - \varepsilon)\mathcal{J}_{12}] \cdot \frac{2}{\varphi(x)} \int_{|x|}^{\infty} uf^M(u) du. \quad (6)$$

Since:

$$\int_{|x|}^{\infty} f^M(u) u du = \frac{1}{2} \varphi(x) = \frac{1}{2} \cdot \frac{1}{\sqrt{\pi}} e^{-x^2}, \quad (7)$$

we can also write:

$$S(x, \tau) = \frac{1}{1 + \zeta} \left[ \varepsilon + (1 - \varepsilon) \cdot \frac{2}{\varphi(x)} \int_{|x|}^{\infty} \tilde{J}_{12}(u, \tau) uf^M(u) du \right] \\ + \frac{\zeta}{1 + \zeta} [\varepsilon + (1 - \varepsilon) \mathcal{J}_{12}]. \quad (8)$$

By inserting this expression for the source function into equation (2) together with equation (4), we can express $\tilde{J}_{12}$ as a function of itself, namely,

$$\tilde{J}_{12}(u, \tau) = \frac{1}{1+\zeta} \mathcal{I} \left\{ \Lambda_{x\mu} \left[ \varepsilon + (1-\varepsilon) \cdot \frac{2}{\varphi(x)} \times \int_{|x|}^{\infty} \tilde{J}_{12}(u', \tau) u' f^M(u') du' \right] \right\} \\ + \frac{\zeta}{1+\zeta} \mathcal{I} \left\{ \Lambda_{x\mu} \left[ \varepsilon + (1-\varepsilon) \cdot 4 \times \int_{0}^{\infty} \tilde{J}_{12}(u', \tau) u'^2 f^M(u') du' \right] \right\}, \quad (9)$$

where $\mathcal{I}$ describes the operation:

$$\mathcal{I}\{f(x)\} = \oint \frac{d\Omega_u}{4\pi} \oint \frac{d\Omega}{4\pi} \int_{-\infty}^{+\infty} \delta(x - \vec{u} \cdot \vec{\Omega}) f(x) dx. \quad (10)$$

The expression in equation (9) provides the mathematical statement of the numerical problem that we intend to solve iteratively.

### 2.2 An approximate operator approach

As usual, we decompose the $\Lambda_{x\mu}$-operator:

$$\Lambda_{x\mu} = \Lambda_{x\mu}^* + (\Lambda_{x\mu} - \Lambda_{x\mu}^*), \quad (11)$$

where $\Lambda_{x\mu}^*$ is an approximate (diagonal) operator of $\Lambda_{x\mu}$. At a given iteration, we update $\tilde{J}_{12}$ by adding a small correction $\delta\tilde{J}_{12}$ to the previously calculated value $\tilde{J}_{12}^\dagger$:

$$\tilde{J}_{12} = \tilde{J}_{12}^\dagger + \delta\tilde{J}_{12}. \quad (12)$$

The key point is then to have an expression for the correction $\delta\tilde{J}_{12}$ at every iteration. Combining these last two equations with equation (9), we have, to first order:

$$\tilde{J}_{12}^\dagger + \delta\tilde{J}_{12} \\ \approx \frac{1}{1+\zeta} \mathcal{I} \left\{ \Lambda_{x\mu} \left[ \varepsilon + (1-\varepsilon) \cdot \frac{2}{\varphi(x)} \int_{|x|}^{\infty} \tilde{J}_{12}^\dagger(u', \tau) u' f^M(u') du' \right] \right\} \\ + \frac{\zeta}{1+\zeta} \mathcal{I} \left\{ \Lambda_{x\mu} \left[ \varepsilon + (1-\varepsilon) \cdot 4 \int_{0}^{\infty} \tilde{J}_{12}^\dagger(u', \tau) u'^2 f^M(u') du' \right] \right\} \\ + \frac{1}{1+\zeta} \mathcal{I} \left\{ \Lambda_{x\mu}^*(1-\varepsilon) \cdot \frac{2}{\varphi(x)} \int_{|x|}^{\infty} \delta\tilde{J}_{12}(u', \tau) u' f^M(u') du' \right\} \\ + \frac{\zeta}{1+\zeta} \mathcal{I} \left\{ \Lambda_{x\mu}^*(1-\varepsilon) \cdot 4 \int_{0}^{\infty} \delta\tilde{J}_{12}(u', \tau) u'^2 f^M(u') du' \right\}. \quad (13)$$

If we look at these terms in detail, we can see that the first two terms only depend on old quantities. Their sum can be interpreted

<footnote n="1">¹The radiation field is axisymmetric in a planar atmosphere, so that it and other quantities such as $\Lambda$ depend only on the cosine of the angle made by the ray with the atmospheric normal, namely $\mu$.</footnote>
as a formal solution i.e. the calculation of $\tilde{J}_{12}$ with a source function that depends only on $S_x^\dagger$:

$$S_x^\dagger = \frac{1}{1+\zeta} [\varepsilon + (1-\varepsilon)\tilde{J}_x^\dagger] + \frac{\zeta}{1+\zeta} [\varepsilon + (1-\varepsilon)\mathcal{J}_{12}^\dagger], \quad (14)$$

with:

$$\tilde{J}_x^\dagger = \frac{2}{\varphi(x)} \int_{|x|}^{\infty} \tilde{J}_{12}^\dagger(u', \tau) u' f^M(u') du'. \quad (15)$$

Defining:

$$I_1 = \mathcal{I}\{\Lambda_{x\mu}^* \cdot \frac{2}{\varphi(x)} \int_{|x|}^{\infty} \delta \tilde{J}_{12}(u', \tau) u' f^M(u') du'\}, \quad (16)$$

and,

$$I_2 = \mathcal{I}\{\Lambda_{x\mu}^* \cdot 4 \int_0^\infty \delta \tilde{J}_{12}(u', \tau) u'^2 f^M(u') du'\}, \quad (17)$$

Equation (13) can be summarized, by writing

$$\mathcal{I}\{\Lambda_{x\mu}[S_x^\dagger]\} = \Lambda_u[S_x^\dagger],$$

as:

$$\delta \tilde{J}_{12} - \frac{1-\varepsilon}{1+\zeta} (I_1 + \zeta I_2) = \Lambda_u [S_x^\dagger] - \tilde{J}_{12}^\dagger = r_u, \quad (18)$$

where $r_u$ can be interpreted as the residual between the formal solution $\Lambda_u[S_x^\dagger]$ and $\tilde{J}_{12}^\dagger$, calculated at velocity $u$. Because $I_1$ and $I_2$ are both integrals on $\delta \tilde{J}_{12}(u')$, the previous equation can be rewritten in a matrix form. We therefore introduce the matrix $A(u, u', \tau)$ defined at all velocities $u$ and $u'$, and all optical depths $\tau$:

$$A \cdot \delta \tilde{J}_{12} = \vec{r}. \quad (19)$$

By inverting matrix $A$, we can compute, at a given iteration $\delta \tilde{J}_{12}$ and thus update the current value of $\tilde{J}_{12}$. Indeed integrals $I_1$ and $I_2$ can be rewritten, by discretizing (i for $u$; j for $u'$), as:

$$I_1 = \sum_j C_{ij} \delta \tilde{J}_{12}^j(\tau), \quad (20)$$

and,

$$I_2 = \sum_j D_{ij} \delta \tilde{J}_{12}^j(\tau), \quad (21)$$

where:

$$C_{ij} = \begin{cases} 2 \frac{\Lambda_x^*(x=0)}{\varphi(x=0)} W_j f_j^M u_j & \text{if } i=1 \\ \frac{1}{u_i} W_j f_j^M u_j \int_{-\vec{u}}^{\vec{u}} \frac{\Lambda_x^*(x)}{\varphi(x)} dx & \text{elsewhere} \end{cases} \quad (22)$$

$$D_{ij} = \begin{cases} 4 \Lambda_x^*(x=0) W_j f_j^M u_j^2 & \text{if } i=1 \\ \frac{2}{u_i} W_j f_j^M u_j^2 \int_{-\vec{u}}^{\vec{u}} \Lambda_x^*(x) dx & \text{elsewhere} \end{cases} \quad (23)$$

(detailed calculation are presented in Appendix A).

In the above equations, $\vec{u} = \min(u_i, u_j)$, and $W_j$ denotes the j-th integration weights; $f_j^M$ is the Maxwellian velocity distribution corresponding to $u'$, and $\Lambda_x^*$ results from the angular integration of $\Lambda_{x\mu}^*$:

$$\Lambda_x^* = \frac{1}{4\pi} \oint \Lambda_{x\mu}^* d\Omega. \quad (24)$$

Now, we can discretize equation (18) as:

$$\delta \tilde{J}_{12}^i - (1 - \varepsilon) \sum_j B_{ij} \delta \tilde{J}_{12}^j = r_i, \quad (25)$$

where:

$$B_{ij} = \frac{1}{1+\zeta} C_{ij} + \frac{\zeta}{1+\zeta} D_{ij}. \quad (26)$$

Comparing equation (19) with equation (18), we can identify the left hand-side of equation (25) as the sum:

$$\sum_j A_{ij} \delta \tilde{J}_{12}^j,$$

wherein the elements of $A$ at each depth are:

$$A_{ij} = \delta_{ij} - (1 - \varepsilon) B_{ij}, \quad (27)$$

with $\delta_{ij}$ the Kronecker delta function.

The form of matrix $A$ given above is very similar to that obtained by Paletou & Auer (1995) for the 'frequency by frequency' (FBF) scheme they developed in order to solve the two-level atom with standard PRD problem. We propose to call this new method as UBU (velocity by velocity). It can be described as:

(i) Step 0: at each optical depth, compute $A_{ij}$ coefficients (only once) and initialize a first value of the partial scattering integral $\tilde{J}_{12}(u, \tau)$;

(ii) Step 1: compute scattering integral $\mathcal{J}_{12}^\dagger(\tau)$ from equation (4);

(iii) Step 2: compute excited atoms velocity distribution function $f_2^\dagger(u, \tau)$ from equation (3) of PSP23;

(iv) Step 3: compute emission profile $\psi^\dagger(x, \tau)$ from equation (4) of PSP23;

(v) Step 4: compute source function $S_x^\dagger(\tau)$ from equation (3);

(vi) Step 5: compute a formal solution $\Lambda_u[S_x^\dagger]$;

(vii) Step 6: compute the residual $r_u = \Lambda_u[S_x^\dagger] - \tilde{J}_{12}^\dagger$;

(viii) Step 7: solve $A \cdot \delta \tilde{J}_{12} = \vec{r}$ to obtain the correction $\delta \tilde{J}_{12}$;

(ix) Step 8: update $\tilde{J}_{12} = \tilde{J}_{12}^\dagger + \delta \tilde{J}_{12}$ and go back to Step (1).

Our UBU method has been coded in Julia because of the significant gain in computation time (≈ 13 times faster than Python to solve the FNLTE problem).

## 2.3 Numerical implementation

In practice, the computation of quantities such as the partial scattering integral, the velocity distribution of excited atoms, the emission profile, and source function remain unchanged, as compared to PSP23. To implement an approximated operator method such as UBU, we only need to compute matrix $A$, and solve the system in equation (19). Care must be taken when computing $A_{ij}$ coefficients to ensure normalization of critical quantities such as emission profile or velocity distributions. With the FBF scheme, we use normalized frequency redistribution functions given a priori. The quantities to be normalized are therefore obvious. In contrast, looking at equations (22) and (23), normalization is a problem by itself. Here, we propose to check the following relation (for more details, see Appendix B):

$$\sum_j C_{ij} = \sum_j D_{ij} = \frac{1}{2u} \int_{-u}^{u} \Lambda_x^* dx \triangleq N_i. \quad (28)$$

Therefore the normalization procedure is:

(i) Step 0: calculate $C_{ij}$ and $D_{ij}$ using equations (22) and (23), then calculate $\sum_j C_{ij}$ and $\sum_j D_{ij}$;

(ii) Step 1: calculate $N_i$ and deduce the ratio $\alpha = N_i / \sum_j C_{ij} = N_i / \sum_j D_{ij}$;

(iii) Step 2: at each depth, normalize by writing $C_{ij}^{\text{new}} = \alpha C_{ij}$ and similarly for $D_{ij}$.
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Frequency variation of normalized source function at different optical depth $\tau$. Horizontal black dashed lines show CRD solutions for $\tau = 0, 1, 10, 100, 10^3, 10^4$. We compare UBU results with 'standard' PRD results (black open circle). CRD and PRD solutions are computed using the ALI and FBF methods, respectively. Clearly we recover the results presented in fig. 1(c) of Hummer (1969) and fig. 1 of PSP23.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Maximum relative error to FBF reference solution in three cases: Initialization of source function to CRD solution for PSP23 method (dashed line) and our new UBU method (full line); initialization of source function to Wien function $B_W$ for PSP23 method (dash-dotted line). These solutions are calculated under the same conditions as Fig. 1.</figcaption></figure>

## 3 RESULTS

PSP23 have shown that, in the case where velocity-changing collisions are negligible ($\zeta = 0$), we have an equivalence between FNLTE and standard PRD, thereby providing a valuable benchmark. For our new scheme, we then compare our results with standard PRD results obtained using the FBF method. We compute the source function $S(x, \tau)$ for an 1D semi-infinite and plan parallel atmosphere of total optical thickness $\tau_{max} = 10^6$, with $\zeta = 0$ and $\epsilon = 10^{-4}$. We sample the optical thickness with a logarithmic grid using 5 points per decade. We perform angular integration with (i) a 6 points Gauss-Legendre quadrature for $\mu$ integration and (ii) a 10 points rectangular quadrature for the azimuths $\varphi \in [0, 2\pi]$. Reduced frequencies, $x$, and normalized velocities $u$ grids are identical; they span from 0 to $x_{max}$, $u_{max} = 6$ with a step size of 0.1. Frequency and velocity integrations are performed using a trapezoidal quadrature. We use the Wien function $B_W$ as the initial guess for the partial scattering integral $\tilde{J}_{12}$. We run the UBU method for 75 iterations [for this case, there is no need to go beyond as the error term defined in equation (29) becomes constant at around 60 iterations].

The UBU solution for the source function is shown in Fig. 1, where we also compare it to the standard PRD results obtained with FBF. Our UBU method clearly reproduces accurately the benchmark result of Hummer (1969), and also that of PSP23 thereby validating our alternative scheme. For a more quantitative comparison we show in Fig. 2 the maximum relative error $E_{rel}$ between the UBU source function, $S$, and from the standard PRD source function, $S^{FBF}$, where $E_{rel}$ is defined as:

$$E_{rel} = \max_{\tau,x} \left| \frac{S - S^{FBF}}{S^{FBF}} \right| \qquad (29)$$

Both the UBU and PSP23 method were initialized with an ALI-CRD solution. In Fig. 2, we also show as dash-dotted black line the maximum relative error to FBF for the PSP23 method, when it was initialized with the Wien function $S = B_W$. Clearly, we see that a Wien initialization is very unsatisfactory for the PSP23 scheme, while the UBU method presented in this paper converges well, under the same conditions, as we have already seen in Fig. 1. This is expected, as the PSP23 scheme is based on a mere lambda iteration-like numerical scheme, while our UBU method is another ALI-based iterative technique. Also, under the same conditions, our new UBU method is significantly faster than PSP23's method, with a gain of about 20 iterations to achieve the solution with the same precision ($E_{rel} \approx 3$ per cent).

However, regardless of the scheme used, it seems difficult to obtain a solution that reproduces the FBF solution with an accuracy better than $E_{rel} \approx 3$ per cent. This can be explained by the radically different methods used in our FNLTE approach.

As already noted above, when compared to the PSP23 method, our new UBU method is less sensitive to initial conditions. This behaviour is illustrated in Fig. 3 where we compute a normalized source function using the UBU method initialized using $S = yB_W$ with $y$ between 0 and 1 (namely, a fraction of the Wien function). Here we chose to show the source function at $\tau = 1$, and at the line centre $x = 0$. To compare these solutions, we also showed the FBF reference solution $S^{FBF}(x = 0, \tau = 1)$. From this figure, we see that our new method is also quite insensitive to the choice of the initial source function.
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Normalized source function $S(x, \tau)/B_W$ computed via UBU method and displayed at $x = 0$ and $\tau = 1$ for an initialization at $S_{\text{init}}/B_W = 0, 0.25, 0.5, 0.75, 1$ (coloured full lines). As a reference for the convergence evaluation, we use the solution computed using the FBF method (dashed black line).</figcaption></figure>

## 4 CONCLUSIONS

We have shown that our new UBU scheme is significantly more robust than the method previously developed by Paletou et al. (2023). In particular, we looked at the convergence of these two schemes and showed that the UBU scheme was quite insensitive to the choice of the initial guess. We also show that fewer iterations were needed to achieve the same level of accuracy. These two elements are important advantages when it comes, in our future work, to dealing with more complex and more realistic cases (natural broadening, multi-level atoms) where the computation time will be more demanding and for which a proper initialization may be an issue.

## ACKNOWLEDGEMENTS

TL is kindly supported by a ‘PhD booster’ granted by the Toulouse Graduate School of Earth and Space Science (TESS; https://tess.omp.eu/). MS would like to thank IRAP, OMP, for providing the financial support to visit IRAP for 15 d in March 2024, when some part of this work was carried out.

## CONFLICT OF INTEREST

Authors declare no conflict of interest.

## DATA AVAILABILITY

The validation of the above-described results uses numerical solutions obtained with the FBF scheme originally described in Paletou & Auer (1995), and using the standard redistribution functions of Hummer (1962). The latter results were not numerized, to the best of our knowledge. However, a *Julia* version of the FBF numerical code will be available to the public via the HAL (https://hal.science/) platform in the coming months, enabling to reproduce easily the main reference results used for Fig. 1.

## REFERENCES

* Auer L., Paletou F., 1994, A&A, 285, 675

* Hubený I., Mihalas D., 2014, Theory of Stellar Atmospheres: An Introduction to Astrophysical Non-Equilibrium Quantitative Spectroscopic Analysis. Vol. 26, Princeton Univ. Press, Princeton

* Hubený I., Oxenius J., Simonneau E., 1983, J. Quant. Spectrosc. Radiat. Transfer, 29, 477

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

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

* Lagache T., Paletou F., Sampoorna M., 2025, A&A, in press

* Landi Degl'Innocenti E., Landolfi M., 2004, Astrophysics and Space Science Library Vol. 307, Polarization in Spectral Lines. Kluwer Academic Publishers, Dordrecht

* Olson G. L., Auer L., Buchler J. R., 1986, J. Quant. Spectrosc. Radiat. Transfer, 35, 431

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

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

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

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

* Rutily B., Chevallier L., 2006, Eur. Astron. Soc. Publ. Ser., 18, 1

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

* Rybicki G., Hummer D., 1992, A&A, 262, 209

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

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

## APPENDIX A: DERIVATION OF THE MATRIX COEFFICIENTS $C_{ij}$ AND $D_{ij}$

To obtain an expression for matrix A appearing in equation (19), in particular for the $C_{ij}$ and $D_{ij}$ coefficients, we need to start from the definition of integrals $I_1$ and $I_2$ expressed in equations (16) and (17). Using for $u \neq 0$,

$$ \oint \delta(x - \vec{u} \cdot \vec{\Omega}) d\Omega_u = \frac{2\pi}{u} H(u - |x|), \quad (A1) $$

where $H$ is the Heaviside function, we can rewrite the operator $\mathcal{I}$ as:

$$ \mathcal{I}\{f(x)\} = \frac{1}{2u} \oint \frac{d\Omega}{4\pi} \int_{-\infty}^{+\infty} f(x)H(u - |x|)dx. \quad (A2) $$

In the definitions of $I_1$ and $I_2$, we have $I_1 = \mathcal{I}\{f_1(x)\}$ and $I_2 = \mathcal{I}\{f_2(x)\}$ with:

$$ f_1(x) = \Lambda_{x\mu}^* \cdot \frac{2}{\varphi(x)} \int_{|x|}^{\infty} \delta \tilde{J}_{12}(u', \tau) f^M(u')u'du', \quad (A3) $$

and

$$ f_2(x) = \Lambda_{x\mu}^* \cdot 4 \int_0^\infty \delta \tilde{J}_{12}(u', \tau) f^M(u')u'^2 du'. \quad (A4) $$

Since only the approximate operator $\Lambda_{x\mu}^*$ is dependent on the direction of the ray (via $\mu$), using:

$$ \frac{1}{4\pi} \oint \Lambda_{x\mu}^* d\Omega = \Lambda_x^*, \quad (A5) $$

we obtain:

$$ I_1 = \frac{1}{u} \int_{-\infty}^{+\infty} H(u - |x|) \frac{\Lambda_x^*}{\varphi(x)} \int_{|x|}^{\infty} \delta \tilde{J}_{12}(u', \tau) f^M(u')u'dx dx, \quad (A6) $$

and,

$$ I_2 = \frac{2}{u} \int_{-\infty}^{+\infty} H(u - |x|) \Lambda_x^* \int_0^\infty \delta \tilde{J}_{12}(u', \tau) f^M(u') u'^2 du' dx. \quad (A7) $$

By writing:

$$ \int_{|x|}^{\infty} g(u') du' = \int_{0}^{\infty} H(u' - |x|) g(u') du', \quad (A8) $$
it can be shown that:

$$ \int_{-\infty}^{+\infty} H(u - |x|)H(u' - |x'|) \frac{\Lambda_x^*}{\varphi(x)} dx = \int_{-\bar{u}}^{\bar{u}} \frac{\Lambda_x^*}{\varphi(x)} dx, \quad (A9) $$

where $\bar{u} = \min(u, u')$.

Finally, by reversing the order of integration between $u'$ and $x$ and using the previous equation, we get:

$$ I_1 = \frac{1}{u} \int_0^\infty \delta \tilde{J}_{12}(u', \tau) f^M(u') u' \int_{-\bar{u}}^{\bar{u}} \frac{\Lambda_x^*}{\varphi(x)} dx du', \quad (A10) $$

and,

$$ I_2 = \frac{2}{u} \int_0^\infty \delta \tilde{J}_{12}(u', \tau) f^M(u') u'^2 \int_{-u}^{u} \Lambda_x^* dx du'. \quad (A11) $$

These last result are only valid if $u \neq 0$. When $u = 0$, the Dirac distribution forces $x = 0$, which in this case is equivalent to:

$$ I_1 = 2 \frac{\Lambda_x^*(0)}{\varphi(0)} \int_0^\infty \delta \tilde{J}_{12}(u', \tau) f^M(u') u' du'; \quad (A12) $$

$$ I_2 = 4 \Lambda_x^*(0) \int_0^\infty \delta \tilde{J}_{12}(u', \tau) f^M(u') u'^2 du'. \quad (A13) $$

By discretizing the integration over $u'$ with the index $j$ and introducing an analogous index $i$ for $u$, we have:

$$ \begin{aligned} I_1 &= \sum_j C_{ij} \delta \tilde{J}_{12}^j \\ &= \begin{cases} 2 \frac{\Lambda_x^*(0)}{\varphi(0)} \sum_j W_j f_j^M u_j \delta \tilde{J}_{12}^j & \text{if } i=1 \\ \frac{1}{u_i} \sum_j W_j f_j^M u_j \int_{-\bar{u}}^{\bar{u}} \frac{\Lambda_x^*}{\varphi(x)} dx \cdot \delta \tilde{J}_{12}^j & \text{else.} \end{cases} \end{aligned} \quad (A14) $$

and,

$$ \begin{aligned} I_2 &= \sum_j D_{ij} \delta \tilde{J}_{12}^j \\ &= \begin{cases} 4 \Lambda_x^*(0) \sum_j W_j f_j^M u_j^2 \delta \tilde{J}_{12}^j & \text{if } i=1 \\ \frac{2}{u_i} \sum_j W_j f_j^M u_j^2 \int_{-u_i}^{u_i} \Lambda_x^* dx \cdot \delta \tilde{J}_{12}^j & \text{else.} \end{cases} \end{aligned} \quad (A15) $$

By identifying equations (A14) and (A15) with, respectively, equations (A12) and (A13), we find the expression for $C_{ij}$ and $D_{ij}$ given in equations (22) and (23).

## APPENDIX B: DERIVATION OF THE NORMALIZATION RELATIONSHIP

The aim here is to obtain an analytical expression for the sums:

$$ \sum_j C_{ij} \quad \text{and} \quad \sum_j D_{ij}. \quad (B1) $$

These sums can be calculated by injecting $\delta \tilde{J}_{12} = 1$ into the integrals $I_1$ and $I_2$. In $u \neq 0$, we then have, using equations (A6) and (A11):

$$ \sum_j C_{ij} = \frac{1}{u} \int_{-\infty}^{+\infty} H(u - |x|) \frac{\Lambda_x^*}{\varphi(x)} \int_{|x|}^{\infty} f^M(u') u' du' dx, \quad (B2) $$

and,

$$ \sum_j D_{ij} = \frac{2}{u} \int_{-\infty}^{+\infty} f^M(u') u'^2 \int_{-u}^{u} \Lambda_x^* dx du'. \quad (B3) $$

Using:

$$ \int_{|x|}^{\infty} f^M(u') u' du' = \frac{1}{2} \varphi(x) \quad \text{and} \quad \int_{0}^{\infty} f^M(u') u'^2 du' = \frac{1}{4}, $$

we have:

$$ \begin{aligned} \sum_j C_{ij} &= \sum_j D_{ij} = \frac{1}{2u} \int_{-\infty}^{+\infty} H(u - |x|) \Lambda_x^* dx \\ &= \frac{1}{2u} \int_{-u}^{u} \Lambda_x^* dx \triangleq N_i. \end{aligned} \quad (B4) $$

This paper has been typeset from a TeX/Latex file prepared by the author.