© The Authors 2024

1 Introduction

The numerical integration of the radiative transfer (RT) equation, namely, the formal solution, is a key step in the solution of the RT problem, both at the local thermodynamic equilibrium (LTE), and in non-LTE conditions. From a numerical point of view, the formal solution consists of applying a suitable numerical method (the formal solver). In the polarized case, the formal solver provides an approximate solution of a linear system of first-order coupled inhomogeneous ordinary differential equations (ODEs). Over the years, the effort of the community has produced an extensive literature on the different formal solvers, for both unpolarized (e.g., Auer 1976, 2003; Mihalas et al. 1978; Olson & Kunasz 1987; Kunasz & Auer 1988; Auer & Paletou 1994) and polarized (e.g., Janett et al. 2017, 2018, and references therein) radiation.

In most RT applications, it is common to express the spatial dependency of the transfer equation in terms of the optical depth. This allows for the use of exponential integrators, while also enforcing the numerical stability of the formal solution (e.g., Janett & Paganini 2018). Crucially, it is common to perform this spatial scale conversion by means of numerical methods, thus introducing numerical errors. While considerable efforts have already been exercised in the search for the best-performing formal solver, the numerical conversion to optical depth has often been overlooked. Indeed, a rigorous and clear investigation into this critical numerical step is still lacking. The aim of this work is to analyze the impact of different numerical schemes for the conversion to optical depth, both individually and combined with different common formal solvers. The final aim is to quantify the impact of the numerical conversion to optical depth on the accuracy of the solution of the non-LTE RT problem for polarized radiation and to identify the relevant issues.

The article is organized as follows. Section 2 introduces the non-LTE RT problem for polarized radiation, focusing on the analytical and numerical conversion to optical depth. In Sect. 3, we present our numerical convergence experiments, analyzing the performances of different numerical conversions to optical depth for different problem settings. Finally, Sect. 4 provides our remarks and conclusions. A detailed description of the optical depth conversion schemes considered in this paper can be found in Appendix A. Appendix B presents the synthetic emergent Stokes profiles.

2 Transfer problem for polarized radiation

The intensity and polarization of a beam of radiation are fully described by the four-component Stokes vector: $$I=(I,Q,U,V)=\left(I_1,I_2,I_3,I_4\right)\in {ℝ}_{+}\times {ℝ}^3,$$

with the intensity I being positive and the polarization being encoded in Q, U, and V. The propagation of a radiation beam in a given medium (e.g., the plasma of a stellar atmosphere) is described by the transfer equation for polarized radiation, which consists of a system of four coupled first-order inhomogeneous ODEs. The RT equation is inherently 1D, as it describes how the radiation gets modified while propagating along a straight line, and is decoupled in frequency and direction. In this work, we thus consider the RT problem for the simplest case of a 1D plane-parallel atmosphere, pointing out that the results and conclusions also hold in a 3D setting.

In a plane-parallel atmosphere (see Fig. 1), all physical quantities only depend on the height, z ∊ [z_min, z_max] ⊂ ℝ, and the transfer equation for a beam of radiation of frequency v ∊ [v_min, v_max] ⊂ ℝ₊ propagating along the direction specified by the unit vector Ω = (θ,χ) ∊ [0, π] × [0,2π) can be written as: $$\cos (\theta )\frac{\text{d}}{\text{d}z}I(z,\Omega ,v)=-K(z,\Omega ,v)I(z,\Omega ,v)+\epsilon (z,\Omega ,v).$$(1)

The propagation matrix K ∊ ℝ^4×4 describes how the medium absorbs radiation, coupling the different Stokes parameters. The emission vector ε = (ε₁, ε₂, ε₃, ε₄) ∊ ℝ₊ × ℝ³ describes the radiation emitted by the plasma in the four Stokes parameters.

In general, K and ε include contributions from both line and continuum processes, respectively labeled with the superscripts ℓ and c, which can simply be added together as: $$\begin{array}{c}K(z,\Omega ,v)=K^{\ell }(z,\Omega ,v)+K^c(z,\Omega ,v),\\ \epsilon (z,\Omega ,v)={\epsilon }^{\ell }(z,\Omega ,v)+{\epsilon }^c(z,\Omega ,v).\end{array}$$

The continuum and line contributions to the emissivity contain in general contributions from both thermal and scattering processes. The values of K^ℓ and ε^ℓ depend on the state of the atom giving rise to the considered spectral line. This state has to be determined by solving the statistical equilibrium (SE) equations, which describe the interaction of the atom with the radiation field, other particles in the plasma, and the possible presence of external fields.

The whole non-LTE RT problem consists in finding a self-consistent solution of the RT equation for the radiation field, I, and of the SE equations for the atomic system. The solution of Eq. (1), given knowledge of the initial conditions and the spatial, angular, and frequency dependence of the propagation matrix, as well as the emission vector at a discrete set of points, is the so-called formal solution (Mihalas 1978; Auer 2003; Janett et al. 2017).

2.1 Conversion to optical depth

Formally, the optical depth scale takes the form of the map: $$\tau :\left[z_{\min },z_{\max }\right]\to \left[0,{\tau }_{\max }(\Omega ,v)\right]\subset {ℝ}_{+},$$

defined as the solution of the initial value problem:

$$\text{d}\tau (z,\Omega ,v)=-\frac{\eta (z,\Omega ,v)}{\cos (\theta )}\text{d}z{\text{\hspace{0.17em}with\hspace{0.17em}}}^1\tau \left(z_{\max },\Omega ,v\right)=0,$$(2)

where η is the total absorption coefficient for intensity, which corresponds to the diagonal element of the propagation matrix K, and θ ∊ [0,π/2) ∪ (π/2,π]. In the definition of Eq. (2), we followed the usual convention that the optical depth is measured along the considered line of sight, inwardly in the solar atmosphere (i.e., from the point of view of an observer). Since η > 0, the function τ is strictly monotone decreasing and, thus, there is a differentiable bijection from [z_min, z_max] to [0,τ_max(Ω, v)], with τ_max(Ω, v) = τ(z_min, Ω, v) > 0. By applying the conversion to optical depth defined by EQ. (2), Eq. (1) takes the equivalent form (e.g., Janett & Paganini 2018) as: $$\frac{\text{d}}{\text{d}\tau }I(\tau ,\Omega ,v)=\frac{K(\tau ,\Omega ,v)}{\eta (\tau ,\Omega ,v)}I(\tau ,\Omega ,v)-\frac{\epsilon (\tau ,\Omega ,v)}{\eta (\tau ,\Omega ,v)},$$(3)

where, for notational simplicity, the dependence of τ on the direction, Ω, and frequency, v, has not been explicitly indicated. Analyzing the transfer equations for polarized radiation in Eqs. (1) and (3), Janett & Paganini (2018) showed that the conversion from the geometrical spatial scale to the optical depth spatial scale mitigates the variation of the propagation matrix elements along the integration path, providing significant stability enhancements in the formal solution². Moreover, the conversion to optical depth is a key step for applying the Diagonal Element Lambda Operator (DELO) methods, a family of well performing formal solvers belonging to the class of exponential integrators (Guderley & Hsu 1972; Rees et al. 1989; Trujillo Bueno 2003; de la Cruz Rodríguez & Piskunov 2013; Štěpán & Trujillo Bueno 2013).

Crucially, to apply a formal solver to Eq. (3), it is necessary to have a certain knowledge of the bijection, τ. The solution of the initial value problem given by Eq. (2) is expressed as: $$\tau (z,\Omega ,v)=-\frac{1}{\cos (\theta )}{{{\displaystyle \int }}^{\text{}}}_{z_{\max }}^z\eta \left(z^{\prime },\Omega ,v\right)\text{d}{z}^{\prime }.$$(4)

In 1D radiative transfer applications, we are mainly interested in the absolute value of the difference of optical depth between two generic height points z and z₀ (assuming z > z₀) for a given direction, Ω. This information is provided by the integral in Eq. (4) evaluated in the interval [z₀, z], namely, $$\text{Δ}\tau (\Omega ,v)=\tau \left(z_0,\Omega ,v\right)-\tau (z,\Omega ,v)=\frac{1}{\cos (\theta )}{{{\displaystyle \int }}^{\text{}}}_{z_0}^z\eta \left(z^{\prime },\Omega ,v\right)\text{d}{z}^{\prime }\ge 0.$$(5)

Fig. 1 1D plane-parallel atmospheric model. The spatial dependency of atmospheric physical quantities is limited to the vertical spatial coordinate z. For simplicity, the τ dependence on the inclination, θ, and the frequency, v, is not indicated.

2.2 Numerical conversion to optical depth

Discretizing the atmosphere with N_z spatial points along the vertical, and providing the corresponding total absorption coefficients $${\left\{z_k\right\}}_{k=1}^{N_z},{\left\{{\eta }_k(\Omega ,v)\right\}}_{k=1}^{N_z},$$

with η_k(Ω, v) := η(z_k, Ω,v) and z_k < z_k+1, it is common to employ piecewise quadratures to compute the integral Eq. (5), where increments can be computed on each cell [z_k, z_k+1] with a numerical quadrature. Defining τ_k(Ω, v) := τ(z_k, Ω, v), the optical depth between the points z_k and z_k+1 for a direction Ω is thus given by $$\text{Δ}{\tau }_k(\Omega ,v)={\tau }_k(\Omega ,v)-{\tau }_{k+1}(\Omega ,v)=\frac{1}{\cos (\theta )}{{{\displaystyle \int }}^{\text{}}}_{z_k}^{z_{k+1}}\eta (z,\Omega ,v)\text{d}z\ge 0.$$(6)

Moreover, defining the variable $$u=\frac{z-z_k}{z_{k+1}-z_k},\text{\hspace{0.17em}\hspace{0.17em}for\hspace{0.17em}}z\in \left[z_k,z_{k+1}\right]\text{,\hspace{0.17em}}$$(7) the integral Eq. (6) over the k-th cell [z_k, z_k+1] can be equivalently written as an integral in the unit interval [0, 1] $$\text{Δ}{\tau }_k(\Omega ,v)=\frac{z_{k+1}-z_k}{\cos (\theta )}{{{\displaystyle \int }}^{\text{}}}_0^1\tilde{\eta }(u,\Omega ,v)\text{d}u,$$(8) where the function $\tilde{\eta }$ encodes the behavior of the function η, namely $\tilde{\eta }(u,\Omega ,v)=\eta (z,\Omega ,v)$ for u and z satisfying the terms of Eq. (7).

The total optical depth measured from the top boundary (where we set τ(z_max, Ω, v) = 0) until the spatial point z_k for an inclination θ is then given by: $${\tau }_k(\Omega ,v)={\displaystyle \underset{j=k}{\overset{N_z-1}{{{\displaystyle \sum }}^{\text{}}}}}\text{Δ}{\tau }_j(\Omega ,v).$$(9)

Numerical approximations of Δτ_k (for k = 1,…,N_z − 1) can be obtained by replacing the integrals in Eqs. (6) or (8) with a numerical quadrature. Crucially, this numerical approximation must be strictly monotone increasing, because one needs to access the values of the inverse τ⁻¹. The approximation of the integral Eq. (6) in the k-th cell [z_k, z_k+1] at a given inclination θ and frequency v can thus be expressed as $$\text{Δ}{\tau }_k(\Omega ,v)={\displaystyle \underset{j\in S_p^q}{{\sum }^{\text{}}}}{\omega }_j(\theta ){\eta }_j(\Omega ,v),$$(10)

where ω_j are the quadrature weights. The general p-point stencil $$S_p^q=\{k-p+q,\dots ,k+q-1\},$$(11)

with p ≥ 2 and 2 ≤ q ≤ p, indicates the spatial grid points considered for a particular quadrature. This particular stencil always includes the indices k and k + 1 . As the notation suggests, the quadrature on the k-th cell could (in principle) require the knowledge of points lying outside the considered cell. Typically, such points are the immediate neighbors. As a paradigmatic example, we present the very common and simple trapezoidal quadrature. Omitting, for notational simplicity, the dependency on the direction Ω and the frequency, v, the linear interpolant approximating $\tilde{\eta }$ for u ∊ [0, 1] can be expressed as $${\tilde{\eta }}^{\text{trap\hspace{0.17em}}}(u)={\eta }_k(1-u)+{\eta }_{k+1}u.$$

In this case, the conversion to the optical depth via Eq. (8) is given by $$\text{Δ}{\tau }_k^{\text{trap\hspace{0.17em}}}=\frac{z_{k+1}-z_k}{\cos (\theta )}\frac{{\eta }_k+{\eta }_{k+1}}{2}.$$(12)

This is commonly known as the trapezoidal quadrature. In terms of the notation introduced by Eqs. (10) and (11), we have p = 2 and q = 2, that is $S_2^2=\{k,k+1\}$, with the weights as follows: $${\omega }_k={\omega }_{k+1}=\frac{z_{k+1}-z_k}{2\cos (\theta )}.$$

We note that this quadrature is monotone and second-order accurate and does not make use of information from outside the considered k-th cell.

Janett et al. (2018) showed that high-order formal solvers require a corresponding high-order numerical approximation of Eq. (4) to maintain their high-degree of accuracy. High-order monotone quadrature schemes can be obtained by replacing linear interpolation with higher order monotone interpolants. In Appendix A, we present some concrete examples of high-order interpolatory quadratures suitable to approximate the integral Eq. (4), such as parabolic, Hermite, and Bezier quadratures. We finally note that for an exponentially stratified atmospheric model, Mihalas (1978) suggested to approximate η in Eq. (4) with an exponential function. Unfortunately, the accuracy of this exponential quadrature strongly depends on the assumption of exponential stratification of the atmosphere, which is not always suitable (see Janett & Paganini 2018).

3 Numerical experiments

This work is aimed at clearly assessing the performances of different numerical optical depth conversion schemes in the synthesis of emergent Stokes profiles. In particular, we first study the convergence rates of different quadrature schemes for the optical depth conversion as a function of the number of spatial grid points (Sect. 3.3). We then analyze the convergence rates (as a function of the number of spatial grid points) of the conversion schemes combined with different formal solvers in the synthesis of the emergent Stokes profiles (Sect. 3.4). In particular, we analyze the trapezoidal, (backward) parabolic, spline, and cubic Hermite conversion schemes³ presented in Appendix A, combined with the DELO-linear, DELOPAR, and DELO-parabolic formal solvers (e.g., Rees et al. 1989; Trujillo Bueno 2003; Janett et al. 2017).

3.1 Physical setting

We solved the non-LTE RT problem for polarized radiation accounting for both line and continuum processes in a 1D plane-parallel atmospheric model. The continuum processes have been modeled following Landi Degl’Innocenti & Landolfi (2004), assuming frequency coherence in the comoving frame for the scattering contribution to the emissivity. The line processes were modeled by considering the simplest case of a two-level atom with an unpolarized and infinitely sharp lower level, applying the theoretical framework of Bommier (1997). We recall that this framework allows for the partial frequency redistribution (PRD) effects to be considered; namely, the correlations between the frequencies of the incoming and outgoing photons in the scattering processes. In particular, we applied the redistribution matrix formalism, through which the scattering contribution (label “sc”) to the line emissivity is given by: $${\epsilon }^{\ell ,\text{sc}}(z,\Omega ,v)=k_L(z){{\displaystyle \int }}^{\text{}}\text{d}{v}^{\prime }{\displaystyle \oint }\frac{\text{d}{\Omega }^{\prime }}{4\pi }R\left(z,\Omega ,{\Omega }^{\prime },v,v^{\prime }\right)I\left(z,{\Omega }^{\prime },v^{\prime }\right),$$(13)

where R ∊ ℝ^4×4 is the redistribution matrix (which contains an analytic solution of the SE equations), I is the radiation field that locally irradiates the atom, and k_L is the frequency integrated absorption coefficient. The complete explicit expressions of the elements of K^ℓ and ε^ℓ for the considered atomic model, as well as those of K^ℓ and ε^ℓ can be found in Alsina Ballester et al. (2017).

To clearly assess the impact of the different optical depth conversion schemes, we analyzed three different settings of increasing complexity. First, we modeled the continuum only. Then we considered the weak photospheric Sr I 4607 Å line, which we modeled in the limit of complete frequency redistribution (CRD), that is, by neglecting any frequency correlations between incoming and outgoing photons in the scattering processes. The CRD limit is obtained by suitably modifying the redistribution matrix of Bommier (1997). Finally, we considered the stronger chromospheric Ca I 4227 Å line, which we modeled, also taking into account the PRD effects, under the angle-average approximation (e.g., Rees & Saliba 1982; Belluzzi & Trujillo Bueno 2014; Alsina Ballester et al. 2017).

For the sake of simplicity, the impact of bulk velocity fields was neglected. Moreover, when present, we only considered a vertical magnetic field. Under such assumptions, it can be easily verified that the absorption coefficient, η, is essentially isotropic and the whole 1D problem is characterized by axial symmetry around the vertical. Because of this symmetry, both the optical depth and the radiation field only depend on the inclination θ (and not on the azimuth χ). Moreover, taking the reference direction for positive Stokes Q parallel to the atmosphere, it can be shown that in the absence of magnetic fields, the only non-zero Stokes parameters are I and Q.

3.2 Numerical setting

After a suitable discretization of the spatial, angular, and frequency domains, the problem can be solved by applying a two-step approach. In the first step, the full nonlinear non-LTE RT problem is solved neglecting polarization, using the RH code of Uitenbroek (2001). In the second step, the problem is solved including polarization phenomena, keeping the population of the lower level calculated at step 1 fixed. By so doing, the problem of step 2 is linear with respect to I, and it is solved by applying a matrix-free GMRES iterative method equipped with a physics-based preconditioner (see Benedusi et al. 2022). More details on this solution strategy can be found in Janett et al. (2021, 2024), Benedusi et al. (2021).

In all the calculations, the considered spectral interval [v_min, v_max] is discretized with a N_v-nodes grid, which is finer in the line core, where the nodes are equally spaced, and coarser in the wings, where the nodes distance increases logarithmically. The angular discretization, needed to evaluate the scattering integral of Eq. (13), is provided by a tensor product quadrature: we considered two six-node Gauss-Legendre grids for the inclination µ = cos(θ) ∊ [−1,0) ∪ (0,1]; whereas, we considered a grid with 8 equidistant nodes for the azimuth χ ∊ (0, 2π]⁴.

To study the convergence properties of the optical depth conversion schemes as a function of the number of spatial grid points, we need a sequence of discrete spatial grids with an increasingly large number of points. The original spatial grid (with corresponding atmospheric quantities) is provided by the 1D semi-empirical plane-parallel solar atmospheric model C of Fontenla et al. (1993), hereafter FAL-C. This atmospheric model discretizes the height interval [−100 km, 2219 km] with N_z = 70 unevenly distributed spatial nodes. Depending on the application, we extracted a particular subdomain from the original FAL-C grid, that is, [z_min, z_max] ⊂ [−100 km, 2219 km]. Then, the logarithm of the atmospheric quantities in this new domain is fitted with a polynomial of a degree of 10, and the sequence of atmospheric models with an increasing number of uniformly distributed grid-points was generated. The number of cells of the j-th refinement is given by $$N_{\text{cells\hspace{0.17em}}}^j:=2^j\cdot N_{\text{cells\hspace{0.17em}}}^0\text{\hspace{0.17em}\hspace{0.17em}with\hspace{0.17em}}j\in \{1,\dots ,8\},$$(14)

and, accordingly, $$N_z^j=N_{\text{cells\hspace{0.17em}}}^j+1.$$

Throughout this work, we always use $N_{\text{cells}}^0=10$.

To properly study the accuracy of a particular optical depth conversion scheme, we evaluated its impact on all the $N_{\text{cells}}^0=10$ “macrocells” of the initial grid. The approximation of the total optical depth of the ĸ-th macrocell (with $\kappa =1,\dots ,N_{\text{cells\hspace{0.17em}}}^0$) for the j-th refinement is then given by $${\mathcal{T}}_{\kappa ,j}(\mu ,v):={\displaystyle \underset{\ell =2j^{j_{(k-1)+1}}}{\overset{2^{j}k}{{\sum }^{\text{}}}}}\text{Δ}{\tau }_{\ell }(\mu ,v).$$

The relative error is defined by $$e_{\tau }(\mu ,v;\kappa ,j):=\left|\frac{{\mathcal{T}}_{\kappa ,j}-{\mathcal{T}}_{\kappa ,\text{\hspace{0.17em}ref\hspace{0.17em}}}}{{\mathcal{T}}_{\kappa ,\text{\hspace{0.17em}ref\hspace{0.17em}}}}\right|,$$(15)

where 𝒯_κ,ref is the reference solution, obtained with the cubic Hermite optical depth conversion and using the most refined spatial grid with $N_{\text{cells\hspace{0.17em}}}^8=2560$. The convergence rate, k_κ, of an optical depth conversion scheme is then provided through a linear regression model, which fits the quantity $$\log \left(\underset{k\in I_k,v}{\max }{e}_{\tau }(\mu ,v;\kappa ,j)\right),$$

as a function of log ($N_z^j$), using the nonlinear least squares approach. By defining $$I_{\kappa }:=\left\{2,\dots ,N_{\text{cells\hspace{0.17em}}}^0-1\right\},$$(16)

we intentionally exclude the first (κ = 1) and last ($\kappa =N_{\text{cells\hspace{0.17em}}}^0$) macrocells from the max_κ, because many higher order schemes cannot be fully implemented in those macrocells, and their inclusion would strongly alter the numerical analysis, especially the value of the convergence rate, k_τ.

Besides the convergence properties of the individual optical depth conversion schemes, it is also important to investigate how these schemes combine with formal solvers for the synthesis of emergent Stokes profiles. The relative error of the emergent profiles using the j-th refined spatial grid is defined by: $$e_X(\mu ,v;j):=\frac{X^j(\mu ,v)-X^{\text{ref}}(\mu ,v)}{X^{\text{ref}}(\mu ,v)+\mathrm{sgn}\left(X^{\text{ref}}(\mu ,v)\right)\cdot {\delta }_X},$$(17)

where X = I, Q, V is the considered Stokes parameter, sgn is the sign function, and δ_X > 0 is an offset used to prevent Eq. (17) from diverging in the case of X^ref → 0. The considered values for δ_X are given along with the results. The reference profile X^ref is obtained with the cubic Hermite optical depth conversion, combined with DELO-parabolic formal solver, in the most refined spatial grid. The convergence rate, k_X, of a formal solver combined with an optical depth conversion scheme is provided through a linear regression model, fitting the quantity of $$\log \left(\underset{\mu ,v}{\max }{e}_X(\mu ,v;j)\right),$$

as a function of log ($N_z^j$), using the nonlinear least squares approach. We note that the emergent Stokes profiles are computed from the converged radiation field by performing a single formal solution on the set of lines-of-sight with µ ∊ {0.1,0.15,0.2,0.25,0.33,0.41,0.5,0.58,0.66,0.75,0.9}.

For the first application, where we modeled the continuum only, we considered a single wavelength at λ = 4605.8 Å, and the spatial domain [z_min, z_max] = [−100 km, 1378 km], extracted from FAL-C and discretized according to Eq. (14). In a second application, we modeled the Sr I 4607 Å line, while also accounting for continuum contributions. We considered the wavelength interval [λ_min, λ_max] = [4605.8 Å, 4611.5 Å] discretized with N_v = 71 frequency nodes, and the spatial domain [z_min, z_max] = [−100 km, 2062 km]. Finally, to model the Ca I 4227 Å line, while also accounting for continuum contributions, we considered the wavelength interval [λ_min, λ_max] = [4220.2 Å, 4235.7 Å], discretized with N_v = 101 frequency nodes, and the spatial domain [z_min, z_max] = [−100 km, 2168 km].

3.3 Optical depth calculations

In this section, we analyze the numerical conversion to optical depth. The magnetic fields are neglected, noting that their impact on the absorption coefficient is negligible for the typical strengths found in quiet or moderately active regions. In the first, third, and fifth rows of Fig. 2, we show the optical depth Eq. (9) as a function of height calculated with four different conversion schemes (trapezoidal, parabolic, spline, and cubic Hermite) for different values of N_z for the continuum, Sr I 4607, and Ca I 4227 settings, respectively. In the calculations, we considered the line of sight µ = 0.17, and a single frequency, v_c, which (for the Sr I 4607 and Ca I 4227 settings) is the line-center frequency. In the second, fourth, and sixth rows of Fig. 2, we also report the corresponding relative error, $$\underset{\kappa \in I_k}{\max }{e}_{\tau }\left(\mu ,v_{\text{c}};\kappa ,j\right),$$

as a function of $N_z^j$, with e_τ given by Eq. (15) and I_κ by Eq. (16)⁵. We also provide the corresponding convergence rate, k_τ.

In all settings, the trapezoidal conversion is second-order accurate, the parabolic scheme is third-order accurate, whereas the spline and cubic Hermite conversions are fourth-order accurate, as summarized in Table 1. In general, the fourth-order conversion schemes outperform lower order methods. Indeed, the trapezoidal and parabolic conversions struggle in providing highly accurate results even for N_z = 81. We also note that the parabolic conversion completely fails in approximating the optical depth for N_z ≤ 21, even providing negative values of τ. The spline scheme also shows a problematic pre-asymptotic behavior, providing very inaccurate results for coarse spatial grids (N_z = 11 , 21 ). For all the considered settings, the cubic Her-mite conversion is the best performing scheme, always showing fourth-order accuracy and providing reliable results even for the coarsest spatial grids.

3.4 Synthesis of emergent Stokes profiles

In Fig. 3, we show, for all combinations of the considered optical depth conversions and formal solvers, the relative error: $$\underset{\mu ,v}{\max }{e}_X(\mu ,v;j),$$

as a function of $N_z^j$, with e_X given by Eq. (17) for X = I, Q. In particular, we show the results for the continuum (top two panels), Sr I 4607 (middle two panels), and Ca I 4227 (bottom two panels) settings in the absence of magnetic fields. For the error expressed in Eq. (17), we used the offset δ_X = 0 for both continuum and Sr I 4607 cases, and the offsets δ_I = 0 and δ_Q = 10⁻¹² for the Ca I 4227 case; this is because, in this last case, the Q profile shows a sign reversal for certain µ. In the figure, we also report the corresponding convergence rates, that are summarized in Table 2.

In the continuum and Sr I 4607 settings, DELO-linear shows second-order accuracy and reasonable relative errors both in I and Q for all the conversion schemes. By contrast, its order of accuracy drops to less than 1.5 for the Ca I 4227 case, producing very large relative errors both in I and Q. As shown in Fig. B.1, the DELO-linear formal solver is very inaccurate in the Ca I 4227 Å line core for both I and Q and this unsatisfying behavior is responsible for the very large errors and the bad convergence rate. In all settings, both the DELOPAR and DELO-parabolic formal solvers are third-order accurate in I and Q for all the considered optical depth conversion schemes, except for the trapezoidal scheme, which drops the convergence rate to second-order. When equipped with a high-order optical depth conversion scheme, these high-order formal solvers clearly out-perform DELO-linear both in accuracy and rate of convergence. We note that the third-order accuracy shown by DELOPAR is in agreement with Janett et al. (2017). This study demonstrated that in the absence of magnetic fields (i.e., for a vanishing off-diagonal elements of K), the DELOPAR effectively performs as a third-order numerical scheme. We additionally note that the use of the parabolic conversion scheme leads to unacceptable errors for N_z ≤ 21.

In Fig. 4, we repeat the same convergence error analysis for the Sr I 4607 setting for the Stokes I, Q, and V, but in the presence of a height-independent vertical magnetic field with B = 100 G. The numerical analysis was carried out using the offsets δ_I = 0 and δ_Q = δ_V = 10⁻¹². Also, in this setting, the whole 1D problem is characterized by axial symmetry around the vertical. We note that the convergence behavior of the different numerical schemes for Stokes I corresponds to the one in the absence of magnetic fields (see the third row of Fig. 3). This confirms that in quiet or moderately active solar conditions, magnetic fields have a negligible impact on the emergent intensity profiles. Both DELO-linear and DELOPAR show a second-order accuracy in Stokes Q and V for all the conversion schemes. We note that the second-order accuracy shown by DELOPAR in the presence of magnetic fields is in agreement with Janett et al. (2017). By contrast, DELO-parabolic is third-order accurate for all the considered optical depth conversion schemes, except for the trapezoidal one, which drops the convergence rate to second-order. As shown in Fig. B.2, the trapezoidal and parabolic conversion schemes yield relevant errors also for the Stokes V.

In conclusion, in order to exploit high-order formal solvers, it is necessary to avoid low-order conversion schemes and to guarantee a sufficiently high number of spatial points.

Fig. 2 Optical depth Eq. (9) as a function of height (first, third, and fifth rows) calculated with the trapezoidal (first column), parabolic (second column), spline (third column), and cubic Hermite (fourth column) conversion schemes for a line of sight with µ = 0.17 and different Nz values (see legends on first column), for the continuum (first row), Sr i 4607 (third row), and Ca I 4227 (fifth row) settings. The second, fourth, and sixth rows show the corresponding relative error as a function of Nz, also reporting the convergence rate, kτ, obtained fitting the encircled data. The reference solution is obtained with the cubic Hermite conversion scheme and using Nz = 2561. The relative error of the parabolic and spline conversions for the coarsest spatial grid is out of the plot range.

Fig. 3 Relative error in the emergent Stokes I (first, third, and fifth rows) and Stokes Q (second, fourth, and sixth rows) as a function of the number of spatial grid points Nz, for the DELO-linear (left), DELOPAR (center), and DELO-parabolic (right) formal solvers combined with the trapezoidal, parabolic, spline and cubic Hermite optical depth conversion schemes in the absence of magnetic fields. The corresponding convergence rates kX are reported in Table 2. The considered settings are: continuum (first and second rows), Sr I 4607 (third and fourth rows), and Ca I 4227 (fifth and sixth rows). The reference profiles are obtained with the cubic Hermite conversion scheme and DELO-parabolic, using Nz = 2561. The stars represent the relative error Eq. (17) obtained using the original domain extracted from the FAL-C atmospheric model.

4 Conclusions

In this paper, we analyze the impact of the numerical conversion from the geometrical to the optical depth spatial scale in the solution of the non-LTE RT problem for polarized radiation. We provide a systematic investigation of this numerical step, considering both the optical depth conversion scheme on its own and how it performs when combined with different formal solvers for the synthesis of emergent Stokes profiles. As a test bench, we considered a 1D semi-empirical atmospheric model in the absence of bulk velocity fields, along with three different RT settings: the continuum, the Sr I 4607 Å line, modeled in the CRD limit, and the Ca I 4227 Å line. In the latter setting, angle-averaged PRD effects were taken into account during the modeling.

For the optical depth conversion, in all cases, the trapezoidal scheme is second-order accurate, the parabolic scheme is third-order accurate, whereas the spline and cubic Hermite schemes show fourth-order accuracy. Notably, the trapezoidal and parabolic conversions struggle to provide highly accurate results, yielding relative errors on the order of 1–10% for the common FAL-C spatial grid. Moreover, the parabolic and spline schemes present an erratic pre-asymptotic behavior, leading to very large errors for coarse spatial grids.

In the synthesis of the emergent Stokes profiles, the convergence rate is impacted by the choice of both the conversion scheme and formal solver. The use of low-order conversion schemes reduces the convergence rate of high-order formal solvers. In particular, high-order accuracy in both the optical depth conversion and the formal solver is necessary to guarantee low relative errors (<10⁻²) in the synthetic emergent Stokes profiles; however, this is valid provided that the number of spatial grid points is sufficiently large to adequately sample the atmospheric structure. To obtain accurate results in practical forward-modeling applications, we advise always equipping high-order formal solvers with a high-order optical depth conversion scheme, such as the cubic Hermite conversion.

Finally, we note that standard high-order quadratures rely on smoothness assumptions on the functions to be integrated. Standard quadratures may indeed poorly perform in the presence of discontinuities, which may drastically increase local errors, reducing the accuracy of the solution and thwarting high-order convergence (e.g., Janett 2019). The performance of the numerical conversion to optical depth in 3D atmospheric models from MHD simulations, showing non-smooth variations of the physical quantities, will be the subject of a future investigation. This future study will also include the analysis of non-oscillatory interpolatory quadratures (see Appendix A).

Fig. 4 Same as Fig. 3, but for the emergent Stokes I, Q, and V for the Sr I 4607 setting, in the presence of a height-independent vertical magnetic field with B = 100 G.

Acknowledgements

The financial support by the Swiss National Science Foundation (SNSF) through grant CRSII5_180238 is gratefully acknowledged. The authors are particularly grateful to the anonymous referee for providing valuable comments that helped to improve the article.

Appendix A Optical depth conversion schemes

We present a series of numerical quadratures for evaluating the optical depth scale obtained by replacing the discrete total opacity η (i.e., the integrand in Eq. (4)) with different interpolants, which are then integrated exactly. For notational simplicity, the dependence of η on the direction and frequency is always omitted.

A.1 Lagrange interpolatory quadratures

The p-order Lagrange polynomial approximating η in the interval z ∊ [z_k, z_k+1] reads $${\eta }^{\text{Lagr\hspace{0.17em}},p,q}(z)={\displaystyle \underset{i\in S_p^q}{{\sum }^{\text{}}}}{\eta }_i{\ell }_i(z),$$(A.1)

where $S_p^q=\{k-p+q,\dots ,k+q-1\}$ is a p-point stencil with p ≥ 2 and 2 ≤ q ≤ p, and l_i are the Lagrange basis polynomials, $${\ell }_i(z)={\displaystyle \underset{m\in S_p^q,m\ne i}{{\prod }^{\text{}}}}\frac{z-z_m}{z_i-z_m},$$

satisfying the relation ℓ_i(z_j) = δ_ij. The interpolation (A.1) exactly matches the data values η_i at positions z_i inside the whole stencil $S_p^q$. The conversion to the optical depth (6) given by the Lagrange interpolant (A.1) is then $$\text{Δ}{\tau }_k^{\text{Lagr},p,q}\left(z_{k-p+q},\dots ,z_{k+q+1},{\eta }_{k-p+q},\dots ,{\eta }_{k+q+1}\right)={\displaystyle \underset{i\in S_p^q}{{\sum }^{\text{}}}}{\eta }_i{L}_i,$$(A.2)

with $$L_i={{{\displaystyle \int }}^{\text{}}}_{z_k}^{z_{k+1}}{\ell }_i(z)\text{d}z.$$

The right-hand side of (A.2) has the same form of (10), and the L_i play the role of the weights, hence the task of computing (6) is reduced to choosing the stencil and computing the weights.

We note that by choosing p = 2 and q = 2, we recover the trapezoidal conversion given by (12). In the following, we provide the explicit expressions for the weights of the backward and forward parabolic quadratures, and of the cubic quadrature.

A.1.1 Backward parabolic quadrature

The backward parabolic interpolation is obtained by choosing p = 3 and q = 2, and consists of determining the quadratic Lagrange polynomial interpolating the data {η_k−1,η_k ,η_k+1} at {z_k−1, z_k, z_k+1}. From (6) we have $$\text{Δ}{\tau }_k^{\text{b-parab\hspace{0.17em}}}={{{\displaystyle \int }}^{\text{}}}_{z_k}^{z_{k+1}}{\eta }^{\text{Lagr},3,2}(z)\text{d}z,$$

and one finds $$\begin{array}{c}\text{Δ}{\tau }_k^{\text{b\hspace{0.17em}parab\hspace{0.17em}}}=-{\eta }_{k-1}\frac{h_k^3}{6h_{k-1}\left(h_{k-1}+h_k\right)}+{\eta }_k\frac{h_k\left(3h_{k-1}+h_k\right)}{6h_{k-1}}\\ +{\eta }_{k+1}\frac{h_k\left(3h_{k-1}+2h_k\right)}{6\left(h_{k-1}+h_k\right)},\end{array}$$(A.3)

with h_k = z_k+1 − z_k. We note that $\text{Δ}{\tau }_1^{\text{b-parab\hspace{0.17em}}}$ is not defined by the quadrature (A.3), and the quadrature on the first cell [z₁, z₂] can be performed using, for instance, the trapezoidal quadrature (12).

A.1.2 Forward parabolic quadrature

The forward parabolic interpolation is obtained by choosing p = 3 and q = 3, and consists of determining the quadratic Lagrange polynomial interpolating the data {η_k, η_k+1, η_k+2} at {z_k, z_k+1, z_k+2}. From (6) we have $$\text{Δ}{\tau }_k^{\text{f-parab\hspace{0.17em}}}={{{\displaystyle \int }}^{\text{}}}_{z_k}^{z_{k+1}}{\eta }^{\text{Lagr},3,3}(z)\text{d}z,$$

and we find: $$\begin{array}{c}\text{Δ}{\tau }_k^{\text{f\hspace{0.17em}parab\hspace{0.17em}}}={\eta }_k\frac{h_k\left(2h_k+3h_{k+1}\right)}{6\left(h_k+h_{k+1}\right)}+{\eta }_{k+1}\frac{h_k\left(h_k+3h_{k+1}\right)}{6h_{k+1}}\\ -{\eta }_{k+2}\frac{h_k^3}{6h_{k+1}\left(h_k+h_{k+1}\right)}.\end{array}$$(A.4)

We note that $\text{Δ}{\tau }_{N_z}^{\text{f-parab\hspace{0.17em}}}$ is not defined by the quadrature (A.4), and the quadrature on the last cell $\left[z_{N_z-1},z_{N_z}\right]$ can be performed using, for instance, the trapezoidal quadrature (12).

A.1.3 Cubic quadrature

The cubic interpolation can be obtained by choosing p = 4 and q = 3, and consists of determining the cubic Lagrange polynomial interpolating the data {η_k−1, η_k, η_k+1, η_k+2} at {z_k−1, z_k, z_k+1, z_k+2}. From Eq. (6) we have $$\text{Δ}{\tau }_k^{\text{cubic\hspace{0.17em}}}={{{\displaystyle \int }}^{\text{}}}_{z_k}^{z_{k+1}}{\eta }^{\text{Lagr,\hspace{0.17em}}4,3}(z)\text{d}z,$$

and we find $$\begin{array}{l}\text{Δ}{\tau }_k^{\text{cubic\hspace{0.17em}}}=-{\eta }_{k-1}\frac{h_k^3\left(2h_{k+1}+h_k\right)}{12h_{k-1}\left(h_{k-1}+h_k\right)\left(h_{k-1}+h_k+h_{k+1}\right)}\hfill \\ +{\eta }_k{h}_k\frac{h_{k-1}\left(4h_k+5h_{k+1}\right)+h_k\left(h_k+h_{k+1}\right)+h_{k+1}\left(h_{k-1}+h_k\right)}{12h_{k-1}\left(h_k+h_{k+1}\right)}\hfill \\ +{\eta }_{k+1}{h}_k\frac{h_{k-1}\left(h_k+h_{k+1}\right)+h_k\left(h_{k-1}+h_k\right)+h_{k+1}\left(5h_{k-1}+4h_k\right)}{12h_{k+1}\left(h_{k-1}+h_k\right)}\hfill \\ -{\eta }_{k+2}\frac{h_k^3\left(2h_{k-1}+h_k\right)}{12\left(h_{k-1}+h_k+h_{k+1}\right)\left(h_k+h_{k+1}\right)h_{k+1}}.\hfill \end{array}$$(A.5)

We note that $\text{Δ}{\tau }_1^{\text{cubic\hspace{0.17em}}}$ and $\text{Δ}{\tau }_{N_z}^{\text{cubic\hspace{0.17em}}}$ are not defined by the quadrature (A.5), and the quadrature on these cells can be performed using, for instance, the trapezoidal quadrature in Eq. (12).

A.2 Cubic Hermite quadratures

The cubic Hermite interpolant approximating $\tilde{\eta }$ for u ∊ [0,1] makes use of both the values of η and its derivative η′, and it is expressed as $$\begin{array}{l}{\tilde{\eta }}^{\text{Herm\hspace{0.17em}}}(u)={\eta }_k\left(1-3u^2+2u^3\right)+{\eta }_k^{\prime }\left(z_{k+1}-z_k\right)\left(u-2u^2+u^3\right)\hfill \\ +{\eta }_{k+1}\left(3u^2-2u^3\right)+{\eta }_{k+1}^{\prime }\left(z_{k+1}-z_k\right)\left(-u^2+u^3\right).\hfill \end{array}$$(A.6)

The conversion to the optical depth in Eq. (8) is given by $$\begin{array}{l}\text{Δ}{\tau }_k^{\text{Herm\hspace{0.17em}}}={{{\displaystyle \int }}^{\text{}}}_0^1{\tilde{\eta }}^{\text{Herm\hspace{0.17em}}}(u)\text{d}u\hfill \\ =\left({\eta }_k+{\eta }_{k+1}\right)\frac{\left(z_{k+1}-z_k\right)}{2}+\left({\eta }_k^{\prime }-{\eta }_{k+1}^{\prime }\right)\frac{{\left(z_{k+1}-z_k\right)}^2}{12}.\hfill \end{array}$$(A.7)

We note that a dataset can be interpolated by applying the above procedure on each interval also without the explicit information on the first derivatives η′_k and η′_k+1 One can indeed employ discrete first derivatives to approximate η′_k and η′_k+1, and this choice is not unique. For instance, approximating η′_k and η′_k+1 with the first-order four-point (nodes z_k−1, z_k, z_k+1 and z_k+2 ) derivatives for nonuniform grids (e.g., Singh & Bhadauria 2009; Janett et al. 2019),⁶ we can recover the Lagrange cubic quadrature using Eq. (A.5).

Another example is given by the cubic spline method, which constructs η by choosing the discrete approximation of η′_k and η′_k+1, making even the second derivative η″ continuous at those nodes. This produces a smoother result, which is also more accurate if the dataset consists of values of a smooth function, but it can also introduce undesirable oscillations. In any case, the quadrature (A.7) is fourth-order accurate provided that the approximation ${\tilde{\eta }}^{\prime }\approx {\eta }^{\prime }$ is at least of third-order (e.g., Dougherty et al. 1989; Janett 2019).

A.3 Non-oscillatory interpolatory quadratures

In 3D radiative transfer applications, the discrete data to be interpolated could present at the same time complex smooth structures and various types of discontinuities. It is well known that standard high-order interpolations (e.g., Lagrange interpolatory quadratures) tend to misrepresent the nonsmooth behavior of a function, introducing oscillations near discontinuities. It is thus necessary to apply high-order well-behaved techniques that are able to interpolate both smooth and discontinuous data.

Based on the cubic Hermite interpolant (A.7), the Shape-Preserving Piecewise Cubic Interpolation (the one applied in this paper) interpolates η, making the first derivative η′ continuous. Moreover, the slopes are chosen in such a way that η preserves the shape of the data and respects monotonicity (Fritsch & Carlson 1980). The fourth-order weighted essentially non-oscillatory (WENO) interpolation technique (Janett et al. 2019) has been devised for use with both complex smooth structures and discontinuities. The fourth-order WENO scheme uses the same four-point stencil as the cubic Lagrange and monotonic cubic interpolations. Another example is given by the Bezier interpolant. To obtain accurate results, de la Cruz Rodríguez & Piskunov (2013) proposed using monotonic quadratic Bezier splines to compute the quadrature in Eq. (6). In the following, we present the quadrature based on the quadratic Bezier interpolations.

A.3.1 Quadratic Bezier quadrature

The quadratic Bezier interpolant approximating $\tilde{\eta }$ for u ∊ [0,1] is expressed as: $${\tilde{\eta }}^{\text{Bez}}(u)={\eta }_k{(1-u)}^2+{\eta }_{k+1}{u}^2+C_{k}2u(1-u),$$(A.8)

with C_k as the so-called Bezier control point of the k-th cell. The conversion to the optical depth for Eq. (8) defined through the quadratic Bezier interpolant is given by $$\text{Δ}{\tau }_k^{\text{Bez}}={{{\displaystyle \int }}^{\text{}}}_0^1{\tilde{\eta }}^{\text{Bez}}(u)\text{d}u=\left({\eta }_k+{\eta }_{k+1}+C_k\right)\frac{z_{k+1}-z_k}{3}.$$(A.9)

The Bezier control point can be set as symmetric via $$C_k=\frac{C_k^0+C_k^1}{2},$$(A.10)

where $$\begin{array}{l}{C}_k^0:={\eta }_k+{\eta }_k^{\prime }\frac{z_{k+1}-z_k}{2},\hfill \\ C_k^1:={\eta }_{k+1}-{\eta }_{k+1}^{\prime }\frac{z_{k+1}-z_k}{2},\hfill \end{array}$$(A.11)

are discrete estimates of η in the middle of the k-th cell. To guarantee monotonicity, the derivative is numerically estimated as $${\eta }_k^{\prime }=\{\begin{array}{ll}\frac{d_{k-1/2}\cdot d_{k+1/2}}{{\alpha }_k\cdot d_{k+1/2}+\left(1-{\alpha }_k\right)\cdot d_{k-1/2}},\hfill & \text{\hspace{0.17em}if\hspace{0.17em}}{d}_{k-1/2}\cdot d_{k+1/2}>0\hfill \\ 0,\hfill & \text{else}\hfill \end{array}$$(A.12)

with $$\begin{array}{l}{\alpha }_k=\frac{1}{3}\left(1+\frac{z_{k+1}-z_k}{z_{k+1}-z_{k-1}}\right),\hfill \\ d_{k+1/2}=\frac{{\eta }_{k+1}-{\eta }_k}{z_{k+1}-z_k},\hfill \\ d_{k-1/2}=\frac{{\eta }_k-{\eta }_{k-1}}{z_k-z_{k-1}}.\hfill \end{array}$$

We note that $C_k^0$ and $C_k^1$ are not defined in the first and in the last cell, respectively.

Appendix B Emergent profiles

The top row of Fig. B.1 shows the emergent I and Q/I profiles as a function of wavelength at µ = 0.1 for the Sr I 4607 Å and Ca I 4227 Å lines in the absence of magnetic fields, for N_z = 41 , which is a typical spatial grid size for semi-empirical atmospheric models. In the bottom row, we also report the corresponding relative errors e_I and e_Q defined by (17). The relative error e_I shows that, independently of the formal solver, the use of the trapezoidal and parabolic conversion schemes yields to relevant errors (almost 10⁻¹ ), that is one order of magnitude worse than the results obtained with the high-order conversion schemes. In particular, the use of the trapezoidal scheme usually leads to an underestimation of the intensity, whereas the use of the parabolic scheme leads to an overestimation. For the high-order conversion schemes, the accuracy in the line core is clearly determined by the formal solver. We finally note that, for the Ca I 4227 Å line, the DELO-linear formal solver is very inaccurate in the line core for both I and Q. This unsatisfying behavior is responsible for the very large errors and the bad convergence rate already noticed in Sect. 3.4. We also note that the more satisfying results in the Q/I profiles are due to a sort of compensation of the errors in I and Q. Indeed, these profiles are usually both either under- or over-estimated. For this reason, fractional polarization signals (i.e., Q/I, U/I, and V/I) are usually less sensitive to the conversion scheme than the I, Q, U, and V signals.

To also analyze Zeeman polarization signals, in Fig. B.2 we show the emergent I and V/I profiles at µ = 1 (solar disk-center) for the Sr I 4607 Å line, in the presence of a height-independent vertical magnetic field with B = 100 G. For this setting and this line of sight, Q/I and U/I signals are zero and, therefore, they are not shown. In the bottom row, we also report the corresponding relative errors e_I and e_V defined by (17), with offsets δ_I = 0 and δ_V = 10⁻¹². We note that the trapezoidal and parabolic conversion schemes yield relevant errors also for the Stokes V, while the more satisfying results in the V/I profiles are due to a sort of compensation of the errors in I and V.

Fig. B.1 Emergent I and Q/I proflies at µ = 0.1 as a function of wavelength (top row) for the Sr I 4607 Å (first and second columns) and Ca I 4227 Å (third and fourth columns) lines (in the absence of magnetic fields), modeled considering Nz = 41 and all 12 combinations of 3 formal solvers, namely DELO-linear (red), DELOPAR (green), and DELOPARABOLIC (blue), and 4 optical depth conversion schemes, namely trapezoidal (dot-dashed), parabolic (dot-dot-dashed), spline (dotted), and cubic Hermite (dashed). The reference solution (black solid line) is obtained with the cubic Hermite conversion, the DELO-parabolic formal solver, and Nz = 2561. The bottom row shows the corresponding relative errors as a function of wavelength.

Fig. B.2 Same as Fig. B.1, but for the emergent I and V/I profiles of the Sr I 4607 A line at µ = 1, in the presence of a height-independent vertical magnetic field with B = 100 G. The grey vertical lines are meant as guide to the eye and indicate the wavelength position of the two maxima of |v/I|.

References

All Tables

All Figures

	Fig. 1 1D plane-parallel atmospheric model. The spatial dependency of atmospheric physical quantities is limited to the vertical spatial coordinate z. For simplicity, the τ dependence on the inclination, θ, and the frequency, v, is not indicated.
In the text

Fig. 2 Optical depth Eq. (9) as a function of height (first, third, and fifth rows) calculated with the trapezoidal (first column), parabolic (second column), spline (third column), and cubic Hermite (fourth column) conversion schemes for a line of sight with µ = 0.17 and different Nz values (see legends on first column), for the continuum (first row), Sr i 4607 (third row), and Ca I 4227 (fifth row) settings. The second, fourth, and sixth rows show the corresponding relative error as a function of Nz, also reporting the convergence rate, kτ, obtained fitting the encircled data. The reference solution is obtained with the cubic Hermite conversion scheme and using Nz = 2561. The relative error of the parabolic and spline conversions for the coarsest spatial grid is out of the plot range.

In the text

Fig. 3 Relative error in the emergent Stokes I (first, third, and fifth rows) and Stokes Q (second, fourth, and sixth rows) as a function of the number of spatial grid points Nz, for the DELO-linear (left), DELOPAR (center), and DELO-parabolic (right) formal solvers combined with the trapezoidal, parabolic, spline and cubic Hermite optical depth conversion schemes in the absence of magnetic fields. The corresponding convergence rates kX are reported in Table 2. The considered settings are: continuum (first and second rows), Sr I 4607 (third and fourth rows), and Ca I 4227 (fifth and sixth rows). The reference profiles are obtained with the cubic Hermite conversion scheme and DELO-parabolic, using Nz = 2561. The stars represent the relative error Eq. (17) obtained using the original domain extracted from the FAL-C atmospheric model.

In the text

	Fig. 4 Same as Fig. 3, but for the emergent Stokes I, Q, and V for the Sr I 4607 setting, in the presence of a height-independent vertical magnetic field with B = 100 G.
In the text

Fig. B.1 Emergent I and Q/I proflies at µ = 0.1 as a function of wavelength (top row) for the Sr I 4607 Å (first and second columns) and Ca I 4227 Å (third and fourth columns) lines (in the absence of magnetic fields), modeled considering Nz = 41 and all 12 combinations of 3 formal solvers, namely DELO-linear (red), DELOPAR (green), and DELOPARABOLIC (blue), and 4 optical depth conversion schemes, namely trapezoidal (dot-dashed), parabolic (dot-dot-dashed), spline (dotted), and cubic Hermite (dashed). The reference solution (black solid line) is obtained with the cubic Hermite conversion, the DELO-parabolic formal solver, and Nz = 2561. The bottom row shows the corresponding relative errors as a function of wavelength.

In the text

	Fig. B.2 Same as Fig. B.1, but for the emergent I and V/I profiles of the Sr I 4607 A line at µ = 1, in the presence of a height-independent vertical magnetic field with B = 100 G. The grey vertical lines are meant as guide to the eye and indicate the wavelength position of the two maxima of |v/I|.
In the text