{'fraction_non_alphanumeric': 0.040845329894465356, 'fraction_numerical': 0.015185409374424296, 'mean_word_length': 4.546343599474529, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Numerical solutions to linear transfer problems of polarized radiation I. Algebraic formulation and stationary iterative methods

Context. The numerical modeling of the generation and transfer of polarized radiation is a key task in solar and stellar physics research and has led to a relevant class of discrete problems that can be reframed as linear systems. In order to solve such problems, it is common to rely on efficient stationary iterative methods. However, the convergence properties of these methods are problem-dependent, and a rigorous investigation of their convergence conditions, when applied to transfer problems of polarized radiation, is still lacking. Aims. After summarizing the most widely employed iterative methods used in the numerical transfer of polarized radiation, this article aims to clarify how the convergence of these methods depends on different design elements, such as the choice of the formal solver, the discretization of the problem, or the use of damping factors. The main goal is to highlight advantages and disadvantages of the different iterative methods in terms of stability and rate of convergence. Methods. We first introduce an algebraic formulation of the radiative transfer problem. This formulation allows us to explicitly assemble the iteration matrices arising from different stationary iterative methods, compute their spectral radii and derive their convergence rates, and test the impact of different discretization settings, problem parameters, and damping factors. Conclusions. The general methodology used in this article, based on a fully algebraic formulation of linear transfer problems of polarized radiation, provides useful estimates of the convergence rates of various iterative schemes. Additionally, it can lead to novel solution approaches as well as analyses for a wider range of settings, including the unpolarized case.


Introduction
The theory of the generation and transfer of polarized radiation is of crucial importance for solar and stellar physics. Currently, one of the main challenges in the field is performing detailed numerical simulations of the transfer of polarized radiation in magnetized and dynamic atmospheres. In particular, it is essential to properly account for scattering processes in complex atomic models and multidimensional geometries, which lead to nonlinear and nonlocal multidimensional problems. In applied science, it is often profitable to convert nonlinear problems into linear ones, allowing the use of the powerful tools from numerical linear algebra, such as stationary or Krylov iterative methods. This is also the case in the field of numerical radiative transfer, where a fair number of problems of practical interest are accurately described through a linear system. For instance, the modeling of the transfer of resonance line polarization considering two-level (or two-term) atomic models can be linearized by assuming that the lower level (or term) is not polarized and that its population is known and fixed a priori (see, e.g., Belluzzi & Trujillo Bueno 2014;Sampoorna et al. 2017;Alsina Ballester et al. 2017).
In principle, linear problems can be solved in a single step with a direct method. However, this task can be highly inefficient or even unfeasible in practice, especially for large problems. Therefore, it is common to rely on efficient stationary iterative methods that, combined with suitable parallelization strategies, make the whole problem computationally tractable.
Trujillo Bueno & Landi Degl'Innocenti (1996) first applied the Richardson (a.k.a. fixed point iteration), Jacobi, and block-Jacobi iterative methods to the numerical transfer of polarized radiation 1 , while Trujillo Bueno & Manso Sainz (1999) applied Gauss-Seidel and successive over-relaxation (SOR) iterative methods to the same problem. In particular, the block-Jacobi iterative method is still frequently used in different scattering polarization problems (e.g., Belluzzi & Trujillo Bueno 2014;Alsina Ballester et al. 2017;Sampoorna et al. 2017;Alsina Ballester et al. 2018). By contrast, Gauss-Seidel and SOR iterative methods are not of common practice (Lambert et al. 2016), because they are nontrivial to parallelize. Despite of its slow con-vergence, the Richardson iterative method is still used in computationally complex transfer problems, where the other iterative methods are not exploitable (e.g., del Pino Alemán et al. 2020;Janett et al. 2021).
A rigorous investigation of the convergence conditions of such iterative methods applied to transfer problems of polarized radiation is still lacking. The specific stability requirements are problem-dependent and often difficult to determine. This paper gives a deeper analysis of these convergence conditions, analyzing their dependence on different design elements, such as the formal solver, the spectral, angular, and spatial numerical grids, and the damping factor.
Section 2 briefly introduces iterative methods, exposes their convergence conditions, and presents the most common stationary methods. Section 3 presents the benchmark analytical problem employed in this paper and exposes its discretization and algebraic formulation in the context of various iterative methods. Section 4 numerically analyzes the impact of different problem elements on the iterative solutions, focusing on the choice of the formal solver, on the various discrete grids, and on damping parameters. Finally, Section 5 provides remarks and conclusions, which are also generalized to more complex problems.

Iterative methods
We consider the linear system where the nonsingular matrix A ∈ R N×N and the vector b ∈ R N are given, and the solution x ∈ R N is to be found. In absence of rounding errors, direct methods give the exact solution to the linear system (1) by a finite sequence of steps. They provide the solution using, for example, a suitable factorization of A, such as Gaussian elimination. Due to their high computational complexity and storage costs, direct methods are often prohibitive in practice, especially for large problems. For this reason, iterative methods are usually preferred for the solution of large linear systems. These methods are convenient in terms of parallel implementation and therefore particularly suitable for today's high performance computing systems. Iterative methods start with an initial guess x 0 ∈ R N and generate a sequence of improving approximate solutions x 1 , x 2 , . . . , x n , x n+1 ∈ R N given by A one-step iterative method computes the approximate solution x n+1 by considering only the previous approximate solution x n . Hence, the transition from x n to x n+1 does not rely on the previous history of the iterative process. Moreover, a stationary iterative method generates the sequence of approximate solutions by the recursive application of the same operator φ. Therefore, stationary one-step iterative methods can be expressed in the form If the operator φ is also linear, the linear stationary one-step iterative method can be expressed in the form where neither the iteration matrix G ∈ R N×N nor the vector g depend upon the iteration n. Hence, a linear stationary one-step iterative method starts from an initial guess x 0 and repeatedly invokes the recursive formula (2) until a chosen termination condition matches (cf. Section 2.4).
In order to study and design linear stationary one-step iterative methods, it is convenient to consider the linear system (1) in the form We recall that multistep iterative methods have been also used in the radiative transfer context. Relaying on multiple previous iterates, multistep iterative methods can accelerate the convergence of one-step iterations. In this regard, we mention Ng (e.g., Olson et al. 1986), ORTHOMIN (e.g., Klein et al. 1989), and state-of-the-art Krylov methods. In particular, Krylov acceleration strategies are the subject of the second part of this paper series (Benedusi et al. 2021).

Convergence conditions
An iterative method is convergent if the approximate solution x n converges to the exact solution x as the iteration process progresses, that is, if x n → x for n → ∞. This implies that the error, defined as e n = x n − x, must approach zero as n increases. Subtracting Equation (3) from Equation (2), one obtains e n+1 = Ge n , indicating that the convergence of the iterative method depends entirely on the iteration matrix G. In fact, the iterative method (2) is convergent for any initial guess x 0 if and only if where ρ(G) is the spectral radius of G (see, e.g., Hackbusch 2016). The quantity ρ(G) is also known as the convergence rate, because it describes the speed of convergence of an iterative process. When employing an iterative method, a rigorous convergence analysis based on the condition (4) is advisable. However, since an analytical expression of ρ(G) is often hard to obtain in practice, heuristic-based analyses are commonly employed to study how the convergence of iterative methods depends on problem and discretization parameters (Hackbusch 2016).

Examples
Linear stationary one-step iterative methods for the solution of (1) usually arise from the splitting with P nonsingular. Rewriting (1) in the form of (3) and introducing iteration indices as in (2), one obtains with G = P −1 R = P −1 (P − A) = Id − P −1 A and g = P −1 b. Crucially, P has to be easier to invert then A, while being a good approximation of it to ensure fast convergence. Sometimes, it is convenient to express iterative methods in terms of the correction In the correction form, the iterative method (6) reads or simply where the correction is given by δx n = P −1 r n and r n = b − Ax n is the iteration residual. Many iterative methods are defined through the following decomposition of the matrix A where D is the diagonal of A, while L and U are its strict lower and strict upper triangular parts, respectively. Different choices of the matrix P yield different iterative methods: Richardson method, P = D, Jacobi method, P = blockdiag(A 1 , . . . , A m ), block-Jacobi method, P = D + L or P = D + U, Gauss-Seidel method, with A 1 , . . . , A m being m square diagonal blocks of A. All these choices of P allow for a convenient, and possibly parallel, computation of P −1 . We note that the iterative method corresponding to P = A always converges in one iteration, resulting in a direct method.
In another perspective, known as preconditioning, all these methods correspond to a Richardson iteration applied to the preconditioned system which is equivalent to the linear system (1). The nonsingular matrix P is known as preconditioner and is usually chosen to reduce the condition number of P −1 A with respect to the one of A, such that the preconditioned system is "easier" to solve with respect to the original one.
Multiplying the preconditioner P (or part of it) by a positive scalar, leads to damped iterative schemes, which are discussed in the following section.

Damped iterative methods
The classical idea of damping (or relaxation) consists in scaling the correction δx n in (7) by a damping factor ω ∈ R >0 , namely The term damped properly holds only for 0 < ω < 1, a.k.a. under-relaxation. In fact, the choice ω = 1 yields the original method, while ω > 1 gives an extrapolated version, a.k.a. overrelaxation. For simplicity, the term damped will be used for all choices of ω. Damped iterations can also be read as that is, the weighted mean of the old iterate x n and one step of the original methodx n+1 = x n + δx n . In order to guarantee the convergence of damped methods, it is possible to obtain bounds on ω as well as optimal choices in terms of convergence rate. In particular, the preconditioned iteration (9) is convergent if and only if λ i being the ith eigenvalue of the matrix P −1 A. If P −1 A has real positive eigenvalues, the optimal convergence rate for the iterative methods presented in Section 2.2 is attained for where λ min and λ max are the minimum and maximum eigenvalues of P −1 A. The SOR method is a damped variant of the Gauss-Seidel method, with P = ω −1 D + L. If A is consistently ordered, G Jac = Id − D −1 A has real eigenvalues, and ρ(G Jac ) < 1, the optimal convergence rate for SOR is attained for A deeper and comprehensive analysis of damped stationary iterative methods is given, for instance, by Hageman & Young (1981); Saad (2003); Quarteroni et al. (2010); Hackbusch (2016).

Termination condition
Typically, iterative solution strategies terminate at the nth step if the norm of the relative residual is smaller then a desired tolerance, that is, if r n 2 / b 2 < tol. However, this stopping criteria becomes less and less reliable as the condition number of A, denoted with κ(A), increases. In fact, the following relation between the relative error and the relative residual holds: Since by construction κ(P −1 A) < κ(A), the same termination condition with respect to the preconditioned system (8), namely, P −1 r n 2 / P −1 b 2 < tol, gives a better approximation of the true relative error.

Benchmark problem
In this section, we present the continuous formulation of a linear benchmark transfer problem for polarized radiation. Secondly, we consider its discretization and derive its algebraic formulation. Finally, we remark how the stationary iterative methods presented in Sections 2.2 and 2.3 are applied in this context.
The problem is formulated within the framework of the theoretical approach described in Landi Degl' Innocenti & Landolfi (2004). Within this theory, a scattering process is treated as a succession of independent absorption and reemission processes. This is the so-called limit of complete frequency redistribution (CRD), in which no frequency correlations between the incoming and the outgoing photons in a scattering process are taken into account. We consider a two-level atom with an unpolarized lower level. Since stimulated emission is completely negligible in solar applications, it is not considered. For the sake of simplicity, the contribution of continuum processes is neglected, as well as that of magnetic and bulk velocity fields. We note, however, that the methodology presented in this section can also be applied to a wider range of settings, such as the unpolarized case, two-term atomic models including continuum contributions, atmospheric models with arbitrary magnetic and bulk velocity fields, and theoretical frameworks accounting for partial frequency redistribution (PRD) effects.

Continuous problem
The physical quantities entering the radiative transfer problem are, in general, functions of the spatial point r, the frequency ν, and the propagation direction Ω of the radiation ray under consideration. In the absence of polarization phenomena, the specific intensity and the atomic level population fully describe the radiation field and the atomic system, respectively. When polarization is taken into account, a more detailed description is instead required. The radiation field is fully described by the four Stokes parameters I i (r, Ω, ν), with i = 1, . . . , 4, standing for Stokes I, Q, U, and V, respectively. The nth energy level of the atomic system, with total angular momentum J n , is fully described by the multipolar components of the density matrix, a.k.a. spherical statistical tensors, ρ K Q,n (r), with K = 0, . . . , 2J n , and Q = −K, . . . , K. We recall that the spherical statistical tensors are in general complex quantities, and that the 0-rank component, ρ 0 0,n , is proportional to the population of the nth level (see Landi Degl'Innocenti & Landolfi 2004, Chapter 3). A two-level atomic system at position r is thus described by the following quantities ρ K Q, (r), and ρ K Q,u (r), with the subscripts and u indicating the lower and upper levels, respectively. The assumption of unpolarized lower level 2 implies that In the absence of magnetic fields, the statistical equilibrium equations of a two-level atom with an unpolarized lower level have the following analytical solution (see Landi Degl'Innocenti & Landolfi 2004, Chapter 14) 3 : where c and h have their usual meaning of light speed and Planck constant, respectively, is the so-called thermalization parameter, W is the Planck function in the Wien limit at the line-center frequency ν 0 , and .
The coefficient w (K) J u J is the polarizability factor 4 , δ (K) u describes the depolarizing rate of the upper level due to elastic collisions, and the frequency-averaged radiation field tensor is given bȳ where T K Q,i (Ω) is the polarization tensor (see Landi Degl'Innocenti & Landolfi 2004, Chapter 5) and φ is the absorption profile. By definition, the maximum rank of the radiation field tensorJ K Q is K = 2 (see Landi Degl'Innocenti & Landolfi 2004, Chapter 5). Consequently, the only nonzero spherical statistical tensors of the upper level ρ K Q,u are those with K ≤ 2. Exploiting the conjugation properties of the spherical statistical tensors (see Landi Degl'Innocenti & Landolfi 2004, Chapter 3), at each position r, the considered two-level atomic system with an unpolarized lower level is thus described by 9 real quantities.
The propagation of the Stokes parameters at frequency ν along Ω is described by the following differential equation where s is the spatial coordinate along the direction Ω, ε i is the emission coefficient of the ith Stokes parameter, and the entries K i j form the 4 × 4 propagation matrix.
In the absence of a magnetic field and lower level polarization, the propagation matrix is diagonal with diagonal element η. The transfer equations for the four Stokes parameters are consequently decoupled and read d ds We note that η does not depend on Ω because bulk velocity fields are neglected, while the source function S i does not depend on ν because of the CRD assumption. Considering dipole scattering and neglecting the continuum contribution, the source function for the ith Stokes parameter is given by (see Landi Degl'Innocenti & Landolfi 2004, Chapter 14) Using (13), the quantities σ K Q , which can be interpreted as the multipolar components of the source function, are given by 4 This quantity is given by where the quantity in curly parentheses is the 6j-symbol (see Landi Degl'Innocenti & Landolfi 2004, Chapter 2).

Discretization and algebraic formulation
In order to replace the continuous problem by a system of algebraic equations, we discretize the continuous variables ν, r, and Ω using, respectively, grids with N ν , N s and N Ω nodes. The quantities of the problem are now evaluated at the nodes at {r k } N s k=1 , {Ω m } N Ω m=1 , and {ν p } N ν p=1 only. For k = 1, . . . , N s , the quantity σ K Q (r k ) is computed through the discrete versions of (18) and (14), namely, where J KQ,i (r k , Ω m , ν p ) depends on the choice of spectral and angular quadratures used in (14) Similarly, the discrete version of (16) reads with Finally, provided initial conditions, the transfer equation (15) can be solved along each propagation direction Ω m and at each frequency ν p to provide the Stokes component I i (r k , Ω m , ν p ) for all k. This numerical step can be expressed as where t i (r k , Ω m , ν p ) represents the radiation transmitted from the boundaries and the coefficients Λ i (r k , r l , Ω m , ν p ) encode the numerical method used to propagate the ith Stokes parameter along the ray Ω m (cf. (C.9)). For notational simplicity, the explicit dependence on variables and indexes is omitted, and Equations (19)-(21) are expressed, respectively, in the following compact matrix form where the vector I ∈ R 4N s N Ω N ν collects the discretized Stokes parameters, S ∈ R 4N s N Ω the discretized source function, and σ σ σ ∈ R 9N s its discretized multipolar components. The explicit expressions of vectors and matrices appearing in (22)-(24) are given in Appendices A-C for a one-dimensional (1D) spatial grid. In general, the entries of the matrix J depend on the type of numerical integration used in Equation (14), while the entries of T depend on the coefficients of Equation (16). Moreover, the entries of Λ depend on the numerical method used to solve the transfer Equation (15) (a.k.a. formal solver), on the spatial grid and, eventually, on the numerical conversion to the optical depth scale (e.g., Janett & Paganini 2018).
By choosing σ σ σ as the unknown vector, Equations (22)-(24) can then be combined into a single discrete problem, namely, with Id − JΛT being a square matrix of size 9N s .
The operator Λ and the vector t depend on the absorption coefficient η, which linearly depends on ρ 0 0, , which in turn enters the definition of the unknown σ σ σ. The discrete problem (25) is consequently nonlinear. However, the nonlinear problem (25) becomes linear if we further assume either that (i) the transfer equation (15) is solved by introducing the optical depth scale defined by dτ(ν) = −η(r, ν)ds, and the optical depth grid is fixed a priori, or that (ii) the density matrix element ρ 0 0, (or equivalently the lower level atomic population) is fixed a priori and, consequently, the absorption coefficient η becomes a constant of the problem. We note that in both cases the operator Λ and the vector t become independent of the unknown σ σ σ.
In the literature, the so-called Λ-operator usually describes the mapping of the discretized source function vector S either to the specific intensity vector I or to the corresponding mean radiation field vector (e.g., Olson et al. 1986;Rybicki & Hummer 1991;Paletou & Auer 1995). This operator was introduced to highlight the linear nature of this mapping inside the whole nonlinear radiative transfer problem. The traditional Λ-operator corresponds either to (24) or to the combined action of (24) and the integral (14).
We note that the linearized system (25) has the same form as (1) with A = Id − JΛT , x = σ σ σ, and b = Jt + c. The matrix A can be assembled by constructing J, Λ, and T , computing the products JΛT , and finally calculating Id − JΛT . Alternatively, the matrix JΛT can be assembled column-by-column, constructing the jth column by applying (20)-(21)-(19) (setting c K q and t i to zero) to a point-like multipolar component of the source function σ σ σ i = δ i j , for i, j = 1, . . . , 9N s .
Sections (3.3)-(3.6) describe stationary iterative methods in the context of radiative transfer problems.

Richardson method
The form of the linear system (25) suggests the operator splitting (5) with P = Id and R = JΛT . The iterative method based on this splitting, that is, is an example of the Richardson method, which is commonly known as Λ-iteration in the context of linear radiative transfer problems. Equation (27) describes the following iterative process: (i) given an initial estimate of σ σ σ, calculate S at each discrete position and direction through (16); (ii) compute the numerical solution of the transfer equation (15) to obtain I at each discrete position, frequency, and direction; (iii) integrate I at each position according to (14) to obtain the frequency averaged radiation field tensor, which is used to update σ σ σ using (18). The damped Richardson method is simply obtained by multiplying the correction δσ σ σ n of the Richardson method by a damping factor ω.

Jacobi method
Considering Equation (25), an usual operator splitting (5) is given by choosing the preconditioner P to be diagonal, with the same diagonal of Id − JΛT , resulting in a minimal cost for the application and computation of P −1 , possibly in parallel. This splitting yields the Jacobi method, while its damped version is obtained by multiplying the correction δσ σ σ n by a damping factor ω.
Trujillo Bueno & Manso Sainz (1999) applied the Jacobi method when dealing with CRD scattering line polarization linear problems. In this context, the Jacobi method is sometimes termed as ALI method or local ALI method (see Hubeny 2003). However, an increasing problem size results in a lower impact of the Jacobi acceleration. For instance, when PRD effects are included, the size of the resulting linear system scales with N s N ν under the angle-averaged approximation (see, e.g., Alsina Ballester et al. 2017), while it scales with N s N ν N Ω in the general angle-dependent treatment (see, e.g., Janett et al. 2021). In both cases, the impact of the Jacobi acceleration is barely perceptible.

Block-Jacobi method
In radiative transfer problems, linear systems often have a natural block structure, which offers different possibilities in defining the blocks. Appendices A-C explicitly illustrate the structure of the matrices that build up the linear system (25) for the 1D case. Considering (25), the preconditioner P can be chosen as the block-diagonal of Id − JΛT . This choice yields the block-Jacobi method, while its damped version is simply obtained by multiplying the correction δσ σ σ n by a damping factor ω. We remark that, within this method, the computation of P −1 requires solving the multiple smaller systems resulting from the diagonal blocks.
Trujillo Bueno & Manso Sainz (1999) first applied the block-Jacobi method to the CRD radiative transfer problem, choosing square blocks that account for the coupling of the different multipolar components of the source function. Moreover, the block-Jacobi method is a common choice when dealing with PRD radiative transfer problems under the angle-averaged approximation (e.g., Belluzzi & Trujillo Bueno 2014;Alsina Ballester et al. 2017;Sampoorna et al. 2017;Alsina Ballester et al. 2018). In this case, the preconditioner is usually built by square blocks of size N ν that account for the coupling of all the frequencies (see the frequency-by-frequency method of Paletou & Auer 1995). In such applications, the block-Jacobi method is referred to as ALI, Jacobi, or Jacobi-based method. Janett et al. (2021) first applied the damped block-Jacobi method to the same problem.

Gauss-Seidel method
Another choice of the operator splitting (5) is given by choosing the preconditioner P as the lower (or upper) triangular part of Id − JΛT . This splitting yields the Gauss-Seidel method. The SOR method is simply obtained by modifying the correction δσ σ σ n with a suitable damping factor ω as explained in Section 2.3. Trujillo Bueno & Manso Sainz (1999) first proposed and applied the Gauss-Seidel and SOR methods to CRD scattering line polarization linear problems, while Sampoorna & Trujillo Bueno (2010) generalized these methods to the PRD case. Trujillo Bueno & Fabiani Bendicho (1995) confirmed with demonstrative results that the optimal parameter ω sor given by (12) is indeed effective in the context of numerical radiative transfer. However, Gauss-Seidel and SOR methods are not common in radiative transfer applications (Lambert et al. 2016), mainly because they have limited and nontrivial parallelization capabilities (see, e.g., the red-black Gauss-Seidel algorithm Saad 2003).

Numerical analysis
It is worth mentioning that realistic radiative transfer problems can be substantially more complex than the one considered here. However, we present a suitable benchmark to better understand how the convergence of various iterative methods depends on different design elements, such as the choice of the formal solver, the discretization of the problem, or the use of damping factors. The conclusions could be then generalized to more complex problems.
In the absence of magnetic and bulk velocity fields, 1D plane-parallel atmospheres satisfy cylindrical symmetry with respect to the vertical. Due to the symmetry of the problem, the only nonzero radiation field tensor components areJ 0 0 andJ 2 0 and, from (18), the only nonzero multipolar components of the source function are σ 0 0 and σ 2 0 . Setting the angle γ = 0 in the explicit expression of the polarization tensor T 2 0,i entering (16), the only nonzero source functions are S 1 and S 2 (i.e., S I and S Q ). Consequently, the only nonzero Stokes parameters are I 1 and I 2 (i.e., I and Q). Moreover, the angular dependence of the problem variables is fully described by the inclination θ ∈ [0, π] with respect to the vertical, or equivalently by µ = cos(θ). The atmosphere is assumed to be homogeneous and isothermal. The thermalization parameter entering (18) is chosen as = 10 −4 , while the damping constant entering the absorption profile φ as a = 10 −3 . A two-level atom with J u = 1 and J = 0 is considered and the depolarizing effect of elastic collisions is neglected. All the problem parameters are summarized below J u = 1, J l = 0, w (0) 10 = w (2) 10 = 1, δ (K) u = 0, a (0,0) 10 = a (2,0) 10 = 1 − , W = 1, = 10 −4 , a = 10 −3 , In 1D geometries, the spatial vector r can be replaced by the scalar height coordinate z. Unless otherwise stated, hereafter, the atmosphere is discretized along the vertical direction through a grid in frequency-integrated optical depth (τ). In particular, a logarithmically spaced grid given by 10 −5 = τ 1 < τ 2 < · · · < τ N s = 10 4 is considered in the following sections. The spectral line is sampled with N ν frequency nodes equally spaced in the reduced frequency interval [−5, 5], and N Ω Gauss-Legendre nodes are used to discretize µ ∈ [−1, 1]. For the sake of clarity, a limited number of formal solvers are analyzed: the first-order accurate implicit Euler method, the second-order accurate DELO-linear and DELOPAR methods, and the third-order accurate DELO-parabolic method 5 . However, the numerical analysis performed here can be unconditionally applied to any formal solver. Further details on the specific properties of the aforementioned formal solvers are given by Janett et al. (2017a,b). Similarly to Trujillo Bueno & Manso Sainz (1999), the block Jacobi preconditioner is used with a block size of two that corresponds to the σ 0 0 and σ 2 0 components local coupling.
In the following sections, we take advantage of the algebraic formulation of the transfer problem (25) to study the convergence of iterative methods presented in Section 2. In particular, we focus on the spectral radius of the iteration matrix ρ(G), with G defined in (26), for various preconditioners P and formal solvers. We remark that the line thermal contributions (encoded in c) and boundary conditions (encoded in t) enter only in the right-hand side of (25) and, consequently, they do not affect the results of this study. Figure 1 displays the magnitude of the entries of the matrix Id−JΛT for various formal solvers. Informally speaking, given the solution vector (cf. Appendix A) σ σ σ = [σ 0 0 (τ 1 ), σ 2 0 (τ 1 ), σ 0 0 (τ 2 ), σ 2 0 (τ 2 ), ..., σ 0 0 (τ N s )σ 2 0 (τ N s )] T , a large value of the entry |Id − JΛT | i j , with i, j ∈ {1, ..., 2N s }, indicates that σ σ σ i depends strongly on σ σ σ j . The asymmetry of the matrices is mainly due to the usage of a nonuniform spatial grid. The matrices entries show that, as expected, for τ ≤ 1 (i.e., i < 90) the solution σ σ σ depends mostly on a region in the neighborhood of τ = 1. As τ increases this effect becomes less evident and the solution depends predominantly on neighboring nodes. Figure 2 illustrates that, once a minimal resolution is guaranteed, varying the number of nodes in the spectral and angular grids N ν and N Ω , respectively, has a negligible effect on ρ(G). Therefore, the setting N ν = 21 and N Ω = 20 remains fixed in the following numerical analysis.

Impact of the formal solver
We now present the convergence rates of various undamped (i.e., with ω = 1) iterative solvers, as a function of the number of spatial nodes for various formal solvers. The numerical results are summarized in Figure 3. For the Richardson method (a.k.a. Λ−iteration) we observe that ρ(G) ≈ 1, resulting in a very slow convergence. Indeed, this method is seldom used in practical applications. Moreover, we generally observe a lower convergence rates as N s increases for all formal solvers. Among the considered formal solvers, DELOPAR shows the best convergence when combined with the Jacobi and Gauss-Seidel iterative methods. For the same iterative methods, DELOPAR is not convergent for N s = 10, because in this case the matrix P −1 (Id − JΛT ) has one negative eigenvalue, and thus (10) is not satisfied. A particular case is given by the DELO-parabolic, for which the Jacobi and block-Jacobi iterative methods never converge.

Impact of the damping parameter
We now discuss the impact of the damping parameter ω introduced in Section 2.3. Figure 4 reports ρ(G) for the damped meth- ods as a function of the damping parameter ω and for multiple formal solvers. The top left panel shows that the impact of the damping parameter on the Richardson method is barely perceivable. Similarly, top right and bottom left panels indicates that the impact of the damping parameter on the Jacobi and block-Jacobi methods is not particularly pronounced. In fact, the spectral radius provided by the optimal parameter does not differ substantially from the value obtained with ω = 1. However, a quite interesting result is revealed: the Jacobi and block-Jacobi methods (at ω = 1) combined with the DELO-parabolic formal solver yield a spectral radius bigger than unity, preventing convergence. However, the use of a suitable damping parameter (ω < 0.81) enforces stability, guaranteeing convergence. An illustrative application of the damped block-Jacobi iterative method to transfer of polarized radiation is given by Janett et al. (2021). The bottom right panel confirms that the impact of the damping parameter in the SOR method is particularly effective. Indeed, a suitable choice of the damping parameter significantly reduces the spectral radius of the iteration matrix, increasing the convergence rate of the SOR method. This is consistent with the results of the previous investigations by Trujillo Bueno & Fabiani Bendicho (1995) and Trujillo Bueno & Manso Sainz (1999). We notice that the optimal damping parameter cannot be calculated with (12) for the DELO-parabolic case, since κ = ρ(G Jac ) is larger then one. Moreover, using the implicit Euler and DE-LOPAR formal solvers, the matrix G Jac has few complex eigenvalues and therefore the estimate (12) is not sharp. Figure 5 collects the convergence rates of the damped Jacobi and SOR methods, using the optimal damping parameters ω opt and ω sor defined in Section 2.3.  (11) and ω sor from (12), respectively. Since (12) cannot be applied if κ < 1, the ω sor corresponding to DE-LOPAR is used for the DELO-parabolic formal solver.

Conclusions
This paper presents a fully algebraic formulation of a linear transfer problem for polarized radiation allowing for a formal analysis that aims to better understand the convergence properties of stationary iterative methods. In particular, we investigate how the convergence conditions of these methods depend on different problem and discretization settings, identifying cases where instability issues appear and exposing how to deal with them.
Numerical experiments confirm that the number of nodes in the spectral and angular grids has a negligible effect on the spectral radius of iteration matrices of various methods and, hence, on their convergence properties. By contrast, a larger number of spatial nodes leads to a reduction of convergence rates for all iterative methods and formal solvers under investigation.
In general, the use of damping parameters can both enforce stability and increase the convergence rate. Additionally, the optimal parameters described in Section 2.3 are indeed effective, especially in the SOR case. In the particular case, where the Jacobi and block-Jacobi methods are combined with the DELOparabolic formal solver, the spectral radius of the iteration matrix becomes larger than unity. Hence, unless a suitable damping parameter is used, this combination is not convergent.
We remark that in this paper, we considered a specific benchmark linear radiative transfer problem. However, the presented methodology can also be applied to a wider range of settings, such as the unpolarized case, two-term atomic models including continuum contributions, 3D atmospheric models with arbitrary magnetic and bulk velocity fields, and theoretical frameworks accounting for PRD effects.
Moreover, the algebraic formulation of the transfer problem is used here exclusively as a practical tool for the convergence analysis. In the setting we have considered, the whole transfer problem is encoded in a dense linear system of size 2N s . Different choices of the unknown vector (I, σ σ σ, or S) result in systems with different size, algebraic properties, and sparsity patterns.
However, this algebraic formulation paves the way to the application of advanced solution methods for linear (or nonlinear) systems, arising from radiative transfer applications. In particular, many state-of-the-art parallel numerical techniques can be employed, such as Krylov methods (see, e.g., Lambert et al. 2015), parallel preconditioners and multigrid acceleration (Hackbusch 2016;Saad 2003). Stationary iterative methods can also be applied as smoothers for multigrid techniques. These methods can also be applied in parallel, using available numerical libraries, which are tailored for high performance computing. In particular, Krylov methods are the subject of study of the second part of this paper series (Benedusi et al. 2021).
The matrix W ∈ R 4N Ω ×9 is generally dense and the matrix T is also given by the Kronecker product T = Id N s ⊗ W, where Id N s is the identity matrix of size N s . j l=k 1 + ∆τ l (µ, ν) and h k (µ, ν) = 1 N s −1 l=k 1 + ∆τ l (µ, ν) .
that is using the same ordering of (A.1). For the sake of clarity, we only considered the implicit Euler method for the explicit assembly of the matrix Λ. However, this analysis can be unconditionally applied to any formal solver, such as exponential integrators or multistep methods. In practice, it would be sufficient to modify (C.6) and (C.7) according to the selected formal solver.