© P. Benedusi et al. 2022

1 Introduction

The magnetic field plays a key role in the physics of the solar chromosphere and transition region, but determining its intensity and orientation in these hot and dynamic layers is a notoriously difficult task. The standard techniques that are routinely applied at the photospheric level are of reduced utility for detecting the elusive magnetic fields in these external regions, and new tailored diagnostic methods are urgently needed. A promising way for investigating the magnetism of the chromosphere and transition region is to exploit the magnetic sensitivity of the linear polarization signals produced by the scattering of anisotropic radiation (scattering polarization) in strong resonance lines (e.g., Trujillo Bueno 2014; Trujillo Bueno et al. 2017). Spectral lines of particular interest are, for instance, HI Ly-α, Mg II k, Ca II K, Ca I 4227, and Na I D₂. The scattering polarization signals of these lines generally show broad profiles with extended wings (e.g., Gandorfer 2000, Gandorfer 2002; Kano et al. 2017) and are sensitive to the presence of weak magnetic fields through the Hanle and magneto-optical effects (e.g. Alsina Ballester et al. 2017). The Hanle effect operates in the core of the line (e.g., Landi Degl’Innocenti & Landolfi 2004), while magneto-optical effects operate in the wings (e.g., del Pino Alemán et al. 2016; Alsina Ballester et al. 2016). Complementary information can be obtained from the circular polarization signals produced by the well-known Zeeman effect. The possibility of exploiting the combined action of the Hanle, magneto-optical, and Zeeman effects in strong resonance lines to investigate the magnetism of the outer solar atmosphere is receiving increasing attention by the scientific community, as testified by the recent efforts from both the observational and theoretical point of view (e.g., Štěpán et al. 2015; Kano et al. 2017; Alsina Ballester et al. 2018; Trujillo Bueno et al. 2018; Ishikawa et al. 2021).

In order to correctly model these scattering polarization signals, it is essential to take frequency correlations between the incoming and outgoing radiation in the scattering processes into account, that is, partial frequency redistribution (PRD) effects (e.g., Faurobert-Scholl 1992; Holzreuter et al. 2005; Belluzzi & Trujillo Bueno 2012; Belluzzi et al. 2012; Supriya et al. 2014). Moreover, deterministic bulk velocity fields in the solar atmosphere must also be considered, because they can significantly modify scattering polarization (e.g., Carlin et al. 2012, 2013). We recall that scattering polarization is intimately related to the geometrical properties (i.e., the anisotropy) of the radiation field that pumps the atoms. It therefore becomes apparent that a general angle-dependent treatment of PRD effects, without introducing simplifying assumptions, such as the angle-averaged one (e.g., Rees & Saliba 1982; Leenaarts et al. 2012; Alsina Ballester et al. 2017), can be essential for an accurate modeling (Janett et al. 2021a). When angle-dependent PRD effects are accounted for, the numerical modeling of scattering polarization signals results in an extremely challenging computational problem with significant memory requirements. Solving the discretized radiative transfer (RT) problem for polarized radiation within a heterogeneous and anisotropic medium is a highly coupled vectorial problem in high dimensions (up to seven for the 3D case). This ambitious task requires the development and application of scalable and efficient computational strategies and the use of last-generation supercomputers.

A remarkable effort has already been made in that direction. Nagendra et al. (2002) and Sampoorna et al. (2008, 2017) proposed solution strategies for the angle-dependent PRD problem for a two-level atom in Hanle and Hanle-Zeeman regimes, respectively. Sampoorna & Nagendra (2015) developed an angle-dependent PRD RT code capable of including vertical bulk velocities. Supriya et al. (2013) solved the angle-dependent PRD problem for an atomic model with hyperfine structure in the nonmagnetic case, while Nagendra et al. (2020) extended that work to the incomplete Paschen–Back effect regime. Finally, del Pino Alemán et al. (2020) modeled the Mg II h & k spectral lines including angle-dependent PRD effects, considering a three-term atom. For computational simplicity, the angle-dependent PRD calculations presented in the aforementioned works only considered homogeneous 1D atmospheric models and/or axial symmetry. It is worth mentioning the work of Anusha & Nagendra (2011, 2012), who pioneered angle-dependent PRD calculations in multidimensional test atmospheric models. A comprehensive review of most of these works is given by Nagendra (2019).

In this article, we first reframe the PRD RT problem for polarized radiation as a linear system by using a suitable physical assumption. Using the 1D case as model problem, we then apply a scalable matrix-free Krylov approach providing a novel preconditioning technique, which is based on the simplifying assumption of complete frequency redistribution (CRD) for the modeling of scattering processes. This strategy provides fast and accurate solutions that allow for an extensive and detailed modeling of the considered scattering polarization signals. The article is organized as follows: in Sect. 2, the RT problem for polarized radiation is introduced and described. In Sect. 3, we present the model problem considered in this work, while in Sect. 4 we present its discretization along with a convenient algebraic formulation. In Sect. 5, we describe the solution strategy and the preconditioning technique. In Sect. 6, scaling and convergence numerical experiments, as well as synthetic emergent Stokes profiles are reported. Finally, Sect. 7 provides remarks and conclusions. This article is the third part of a series: a convenient algebraic formulation of the linear RT problem in the CRD limit was introduced in Janett et al. (2021b), while preconditioned Krylov methods were described and applied to the same problem in Benedusi et al. (2021).

2 Transfer Problem of Polarized Radiation

The intensity and polarization of a beam of radiation are fully described by the four Stokes parameters I, Q, U, and V ∈ ℝ, which are usually combined in the four-component Stokes vector

with the intensity I being positive and the polarization being encoded in Q, U, and V. The Stokes vector and the other physical quantities entering the RT problem are in general functions of time t, of the spatial point r, and of the frequency ν and propagation direction Ω of the considered radiation beam. We consider the stationary setting, assuming that all the quantities are independent of time.

2.1 Transfer Equation of Polarized Radiation

When a beam of radiation propagates in a medium (e.g., the plasma of a stellar atmosphere), it interacts with the atoms, molecules, and other particles present therein, and its intensity and polarization are modified. This modification is fully described by the RT equation for polarized radiation, which is integro-differential by nature. In a spatial domain D ⊂ ℝ^d with d ∈{1, 2, 3}, this equation can be written as (1)

where i = 1,…, 4, r ∈ D, Ω = (θ,χ) ∈ [0, π] × [0, 2π) is a unit vector, and ν ∈ [v_min, v_max] ⊂ ℝ₊. Equation (1) is equipped with suitable boundary conditions on ∂D, and describes a set of initial values problems (IVPs). The 4 × 4 propagation matrix

describes how the medium absorbs radiation, accounting for the coupling between different Stokes parameters. A detailed derivation of the explicit expression of the elements of K for different atomic models can be found in Landi Degl’Innocenti & Landolfi (2004). The emission vector

describes instead the radiation emitted by the plasma in the four Stokes parameters.

The elements of K and ε generally depend on the state of the atom, which in turn has to be determined by solving a set of rate equations, usually under the assumption of statistical equilibrium. These equations are local and describe the interaction of the atomic system with the radiation field (radiative processes), with the other particles present in the plasma (collisional processes), and with the possible presence of external fields (e.g., a magnetic field). When the statistical equilibrium equations have an analytical solution, the emissivity can be directly expressed as function of the radiation field that illuminates the atom. More precisely, the emission vector can be written as the sum of two terms, (2)

The first term on the r.h.s. describes the contribution to the emissivity brought by atoms that are collisionally excited (collisional or thermal term). This term only depends on the local atmospheric model quantities; its explicit expression for the case of a two-level atom, under the assumption of isotropic collisions (which is considered throughout this work), is reported in Appendix A. The second term describes the contribution brought by radiatively excited atoms (scattering term). This term clearly depends on the radiation field I. By using the redistribution matrix formalism (which is particularly useful for describing PRD phenomena), the scattering term is given by the following integral operator: (3)

where k_L is the frequency-integrated absorption coefficient and R_ij are the elements of the so-called redistribution matrix, which couples all Stokes parameters, all directions, and all frequencies at each spatial point. We follow the convention that primed and unprimed quantities refer to the incident and scattered radiation, respectively.

2.2 Scattering Processes

There are two main reasons why the inclusion of scattering processes makes the numerical approach to the RT problem particularly complex. First, the scattering term (3) locally connects all frequencies ν and directions Ω considered in the RT problem. Therefore, the transfer Eq. (1) must be simultaneously solved for all the frequencies and directions. Second, the evaluation of the integral (3) can be very expensive, depending on the considered scattering model. The redistribution matrix generally couples all the incident (ν, Ω) and scattered (ν′, Ω′) frequencies and directions. However, this strong coupling is mitigated either by the use of light-weight models for the description of scattering (e.g., the assumption of CRD) or by the use of simplifying approximations (e.g., angle-averaging). In this section, we briefly discuss the limit of CRD, and the main simplifying assumptions for treating the general PRD case.

2.2.1 Limit of Complete Frequency Redistribution

Although scattering is intrinsically a second-order process in a perturbative expansion of matter-radiation interaction, under given circumstances, it can be suitably described as a temporal succession of independent first-order absorption and re-emission processes, in which no correlation between the frequencies of the incoming and outgoing radiation can be accounted for (e.g., Casini & Landi Degl’Innocenti 2008). This is generally referred to as the CRD limit. It can be shown that CRD is a suitable assumption either when the rate of elastic collisions (which relax the coherence of scattering) is very high, or when the radiation field that illuminates the atom is spectrally flat (e.g., Casini & Landi Degl’Innocenti 2008). A solid CRD theory for the generation and transfer of polarized radiation, based on the flat-spectrum assumption, is the one described in Landi Degl’Innocenti & Landolfi (2004). Through this theory, it is possible to suitably model the scattering polarization signals of weak photospheric lines, such as the Sr I line at 4607 Å (e.g., del Pino Alemán & Trujillo Bueno 2021), as well as the Doppler core signals of strong chromospheric lines, such as the HI Ly-α line at 1215 Å (e.g., Štěpán et al. 2015).

2.2.2 Coherent Scattering and Collisional Redistribution

A proper description of polarization phenomena that takes frequency correlations between the incoming and outgoing radiation in the scattering processes into account requires the application of suitable theoretical frameworks. Various approaches are available today. We mention here the one based on the Kramers–Heisenberg scattering equation presented in Stenflo (1994), the one based on the idea of metalevels (Weisskopf & Wigner 1930) proposed by Landi Degl’Innocenti et al. (1997), and that of Bommier (1997a,b, 2017), which is based on a high-order perturbative expansion of matter-radiation interaction and self-consistently includes frequency redistribution effects due to elastic collisions. Finally, we recall the approach based on a diagrammatic treatment of atom-photon interaction proposed by Casini et al. (2014, 2017a,b).

It can be shown that the most general form of the redistribution matrix for an atomic system with infinitely sharp lower states is given by the linear combination of two terms (e.g., Bommier 1997a,b)

where R^II describes scattering processes that are coherent in frequency in the atomic rest frame, while R^III describes scattering process that are totally uncorrelated in the same reference frame. These two terms contain branching ratios that depend on the rate of elastic collisions (see next section for more details), so that the linear combination of R^II and R^III allows accounting for the partial frequency redistribution effects induced by elastic collisions.

2.2.3 Doppler Redistribution

By means of a standard procedure (e.g., Hummer 1962; Mihalas 1978), the redistribution matrix can be transformed from the atomic reference frame into the observer frame by taking into account the Doppler effect due to both thermal motions and bulk velocities. We assumed that thermal velocities are described by a Maxwellian distribution, and we accounted for arbitrary bulk velocities. The Doppler effect causes frequency redistribution effects in the observer frame, inducing a complex coupling between frequencies and directions. This coupling makes the numerical problem extremely demanding from the computational point of view. For this reason, simplifying approximations have been proposed in the past and are still commonly used, such as the angle-averaged one for the case of (e.g., Mihalas 1978; Rees & Saliba 1982; Leenaarts et al. 2012; Alsina Ballester et al. 2017), and the approximation of CRD in the observer frame for (e.g., Mihalas 1978; Sampoorna et al. 2017; Alsina Ballester et al. 2017).

3 Model Problem

We considered a two-level atom with unpolarized and infinitely sharp lower level in the presence of magnetic and bulk velocity fields. We used the redistribution matrix for this atomic model derived by Bommier (1997a,b). We considered the exact angle-dependent expression of in the observer frame and the approximation of CRD in the observer frame for , which is a suitable choice for modeling strong chromospheric lines (see Janett et al. 2021b). When the formalism of irreducible spherical tensors for polarimetry is used (see Chapt. 5 of Landi Degl’Innocenti & Landolfi 2004) and the impact of an arbitrary bulk velocity is accounted for, the and redistribution matrices can be written as (4) (5)

where are the geometrical tensors defined in Landi Degl’Innocenti & Landolfi (2004, Sect. 5.11), evaluated in a reference system with the quantization axis directed along the magnetic field. The explicit expressions of and can be recovered from Alsina Ballester et al. (2017) for the case without bulk velocity. Fully equivalent expressions were also derived by Sampoorna (2011), working within the Kramers– Heisenberg approach of Stenflo (1994). A bulk velocity field only affects the frequencies of the incident and scattered radiation perceived by the atom. This effect can thus be included by introducing the corresponding Doppler shifts in the and expressions. The quantities α_Q and are the branching ratios for and , respectively, with

where Γ_ℝ, Γ_I, and Γ_E are the rates for radiative de-excitation, inelastic collisions, and elastic collisions, respectively, and the quantity D^(K) is the depolarizing rate of elastic collisions (e.g., Landi Degl’Innocenti & Landolfi 2004, Sect. 7.13). The imaginary term in the denominator describes the Hanle effect, v_L is the Larmor frequency, and g_u is the Landé factor of the upper level.

We note that in the absence of elastic collisions (i.e., Γ_E = D^(K) = 0), the branching ratio for vanishes and the total redistribution matrix reduces to ; this is the limit of coherent scattering in the atomic frame. If the rate of elastic collisions is very high (i.e., Γ_E ≫ Γ_ℝ, Γ_I), then the branching ratio for vanishes and the total redistribution matrix reduces to ; this is the limit of CRD in the atomic frame¹. In this case, the depolarizing rate of elastic collisions becomes very high, that is, D^(K) ≫ Γ_r, Γ_I for K ≠ 0 (we recall that D⁽⁰⁾ = 0), and polarization phenomena thus become negligible. We finally note that the factor

is common to both α_Q and and represents the branching ratio of the scattering contribution to the total emissivity in Eq. (2). In the case Γ_I ≫ Γ_ℝ, the scattering term vanishes and the emis-sivity is dominated by the collisional term; this is the limit of local thermodynamic equilibrium (LTE), which is a light-weight RT problem because the emissivity does not depend on the radiation field I, thus decoupling all frequencies ν and directions Ω considered in the problem.

Fig. 1 Cartesian reference system considered in the problem. The z-axis is directed along the local vertical. The angles θ (inclination) and χ (azimuth) specifying a given propagation direction Ω are indicated in the figure. The considered atmospheric model extends from zmin to zmax. The various physical quantities are constant over the horizontal x-y planes.

3.1 Linearization

The problem presented in Sect. 2 is in general nonlinear. For the considered atomic model, the elements of both K and ε depend on the coefficient k_L. This quantity is proportional to the population of the lower level, which in turn depends on I in a nonlinear way through the statistical equilibrium equations. The redistribution matrix R only depends on local atmospheric quantites, therefore a suitable assumption to retrieve linearity with respect to I is to fix the coefficient k_L a priori (see Belluzzi & Trujillo Bueno 2014; Alsina Ballester et al. 2017; Janett et al. 2021a). In this way, K is independent of I and ε linearly depends on it. The whole problem (1) thus becomes linear in I, since it consists of the set of linear IVPs (1) linearly coupled through the scattering term (3).

3.2 Continuum Processes

So far, we have only considered the contribution to the emission vector due to line processes, that is, due the transition between the upper and lower level of the considered atomic model. However, the proposed formalism allows for a straightforward inclusion of continuum processes. Continuum processes contribute to the emissivity through a scattering term and a thermal term that are formally analogous to these brought by line processes (see Appendix B). It is therefore possible to consider both line and continuum processes via expressions that are formally identical to those introduced in the previous section. Under the typical conditions found in the solar atmosphere, continuum processes only contribute to the diagonal elements for the propagation matrix (see Appendix B).

3.3 Plane-Parallel Atmospheric Model

We considered a semi-empirical plane-parallel model of the solar atmosphere, in which all physical quantities are constant over horizontal planes. In this setting, r = (x, y, z) can be replaced by the vertical coordinate z ∈ [z_mm, z_max] ⊂ ℝ, see Fig. 1. With this simplification, the RT Eq. (1) can be rewritten as (6)

For in-going rays at the domain boundary ∂D, Eq. (6) is equipped with the following boundary conditions: (7)

At the lower boundary, we considered an isotropic and unpolarized incident radiation I_in(v) = [B(v), 0, 0, 0]^T, where B(v) is the Planck function at the effective temperature in z_min.

4 Discretization and Algebraic Formulation

In astrophysical applications, the discretization of the RT problem described by Eqs. (1) and (3) is usually performed in terms of the discrete ordinates method (DOM). Within this method, the continuous problem is discretized on a predefined set of directions determined by the quadrature rule chosen to approximate the angular integral appearing in Eq. (3), leading to a set of semi-discrete RT equations. This angular discretization is then coupled with the spatial one, usually provided by the considered atmospheric model, and the spectral one, tailored to the considered spectral line. The DOM method is easy to implement, usually performs with good accuracy and flexibility, and leads to algorithms of reasonable computational efficiency. We considered a double Gauss-Legendre angular discretization coupled with an exponential integrator method in space; these quadrature methods are briefly described in Sects. 4.2 and 4.3, respectively. Finally, the considered spectral interval was discretized through a grid, which provides a suitable sampling of the considered spectral line.

4.1 Frequency Discretization

The considered spectral interval [v_min, v_max] is discretized with a unevenly spaced grid with N_v nodes, namely

The frequency grid is usually finer in the line core, where the nodes are equally spaced, and coarser in the wings, where the node distances increase logarithmically.

4.2 Angular Discretization

For the angular discretization of Ω = (θ, χ), we used a tensor product quadrature. For the inclination µ = cos(θ) ∈ [−1, 1], we considered two Gauss-Legendre grids (and corresponding weights) with N_θ/2 nodes each, namely

for µ ∈ (−1, 0) and µ ∈ (0, 1), respectively. This grid corresponds to

with θ_j = arccos(µ_j) for j = 1,…, N_θ, with Ν_θ even. For the azimuth χ ∈ (0, 2π], we considered an equidistant grid (and corresponding trapezoidal weights) with N_χ nodes, namely

4.3 Spatial Discretization and Formal Solution

The discretization of the spatial domain D is usually given by the considered discrete atmospheric model. We considered the 1D semi-empirical model C of Fontenla et al. (1993), in which the spatial domain is discretized with an unevenly spaced grid with N_s nodes, namely (8)

For (θ_j, χ_k, v_n) and a corresponding emissivity vector ε at all , the solution of the semi-discrete ODEs system (6) is numerically approximated in applying a suitable integrator: this process is known as the formal solution. The exponential integrators class (e.g., Guderley & Hsu 1972) provides different suitable formal solvers, known as short-characteristic methods in the unpolarized case and as DELO methods in the polarized case (e.g., Rees et al. 1989; Janett et al. 2017a).

For the sake of clarity, we briefly recall the idea that lies behind DELO methods. Equation (6) is first rewritten as (9)

where η₁ is the diagonal element of the propagation matrix and

For a given ray (θ, χ, ν), the exact integration of Eq. (9) in the interval [z₀, z] yields the variation of constants formula (10)

where (11)

A large variety of numerical schemes are available to approximate the integral in Eq. (10). In particular, a linear approximation of S_i in the integration interval produces the DELO-linear method, which is an L-stable method particularly suited for the solution of Eq. (6) since it guarantees stability by exactly resolving the exponential decay of the solution due to absorption processes (e.g., Janett & Paganini 2018). For the numerical experiments presented in Sect. 6, the DELO-linear method was used as formal solver.

We finally remark that DELO methods requires the spatial scale conversion from the geometrical scale to the optical depth scale described by Eq. (11); details about this conversion are given in Appendix C. We refer to Janett et al. (2017a,b, Janett et al. 2018) and Janett & Paganini (2018), for a detailed review of DELO methods and their convergence and stability properties in the context of polarized radiation transfer.

4.4 Algebraic Formulation

As presented in Janett et al. (2021b) and Benedusi et al. (2021) for a simplified setting, we considered the algebraic formulation of Eq. (6), which is a key step for the setup of an efficient solution strategy. For the total number of unknowns N = 4N_sΝ_θΝ_χΝ_ν, we considered the collocation vectors of the Stokes and emission vectors, I and ϵ ∈ ℝ^N, respectively, ordered with a lexicographic criterion, namely

For example, the first 4Ν_θN_χΝ_ν elements of I and ϵ correspond to the first spatial point z₁. We refer to algorithms 1–2 for the precise indices correspondence relation. For notational simplicity, an abuse of notation is employed: the symbol I_i is used to denote both the continuous solution of Eq. (1) and its numerical approximation.

The solution of the RT problem is divided into two components: the solution of the transfer equation (formal solution), and the calculation of the emissivity. We recall that both IVP (1) and the integral (3) are linear with respect to I and ε. In a discrete emissivity field ϵ, the solution of all discretized transfer IVPs in Eq. (6) can be written as (12)

where Λ : ℝ^N → ℝ^N is the linear transfer operator that encodes the numerical scheme (e.g., DELO-linear) and the propagation matrix K. An example of the explicit assembly of Λ is given in Janett et al. (2021b). The action of Λ can be encoded in matrixfree form (see algorithm 1). The vector t ∈ ℝ^N encodes boundary conditions and can be computed using a modification of algorithm 1, setting ϵ= 0 and using the prescribed I_in instead of the zero boundary conditions.

Similarly, the discrete computation of the emissivity of Eq. (2) can be written as (13)

where the linear operator Σ : ℝ^N → ℝ^N encodes the discretized scattering integral (3) and depends on the chosen numerical quadratures. When accounting for PRD effects, we can write the action of Σ as the sum of two components, namely (14)

where Σ^II and Σ^III encode the contributes of R^II and R^III , respectively. As for the operator Λ, the action of Σ can be encoded in a matrix-free form (see algorithm 2). The vector ϵ^th ∈ ℝ^N encodes the thermal emissivity.

Combining Eqs. (12) and (13), we can formulate the whole discrete RT problem as a single-size N linear system with unknown I, namely (15)

The right-hand side vector Λϵ^th + t can be computed a priori by solving IVP (6)–(7) with thermal contributions alone, that is, by performing a single formal solution with ɛ^SC = 0. The matrixfree action of the operator Id – ΛΣ is described by algorithm 3.

Since anisotropic and heterogeneous coefficients as well as unevenly spaced grids are present, the matrices Λ, Σ, and Id – ΛΣ are usually not symmetric. In particular, the matrix Id – ΛΣ is an M-matrix and is nonsingular if ρ(ΛΣ) < l, where ρ is the spectral radius. Figure 2 shows an example of the sparsity pattern of these three matrices. We remark that the block diagonal structure of Σ is well suited for a parallel evaluation on each spatial node,

as the structure of algorithm 2 suggests. Moreover, the action of A couples different spatial nodes, but each ray (θ_j,χ_k,v_n) can be evaluated independently. As sparsity patterns suggest, the application of Λ and Σ comes at a very different computational costs; algorithm 1 indicates that the application of Λ requires O(N) floating point operations, while algorithm 2 shows that O(N N_θ N_v N_χ) operations are required to apply Σ. Since the redistribution matrix R^III is separable, the operator Σ^III has complexity of O(N) and the cost of Σ is thus dominated by the Σ^II term. Figure 3 reports an example of run times for the application of the three operators, showing that the cost of applying Λ and Σ^III is indeed negligible when compared to the cost of Σ^II.

5 Scalable Solution Strategy

The structure of the linear system in Eq. (15) suggests the use of a fixed-point iteration (also known as Richardson method or Λ-iteration), that is, iterating between Eqs. (12) and (13) until convergence. However, this iteration converges too slowly for practical applications, and more effective strategies are often employed, such as Jacobi or block-Jacobi preconditioners (e.g., Trujillo Bueno & Manso Sainz 1999; Janett et al. 2021b) and multigrid techniques (e.g., Fabiani Bendicho et al. 1997).

More recently, preconditioned Krylov methods gained popularity in the RT context (e.g., Warsa et al. 2003; Anusha et al. 2009; Ren et al. 2019; Badri et al. 2018, 2019; Benedusi et al. 2021). We recall that the explicit assembly of the dense matrix Id – ΛΣ in Eq. (15) is indeed not suitable for practical applications because N typically exceeds 10⁶ for 1D models and reaches 10¹⁰ for full 3D models. Crucially, Krylov methods are well suited for matrix-free approaches and have proved to be highly effective solution strategies for large linear systems.

Fig. 2 Sparsity patterns of Λ, Σ, and Id – ΛΣ using Ns = 10, Nθ = 2, Nχ = 1, and Νν = 10, i.e., N2 = 8002 = 6.4 × 105 elements. The corresponding number of nonzero elements (nz) is reported. The elements in the lower and upper triangular parts of Id – ΛΣ correspond to downward (negative µ) and upward (positive µ) directions, respectively.

Fig. 3 Run times for a single application of Σ, ΣIII, and Λ using algorithms 1–2 varying the number of processors Np ≤ Ns with Ns = 70, Nθ = 12, Nχ = 8, and Nv = 99.

5.1 Preconditioning

Preconditioning can significantly increase robustness and convergence speed of iterative techniques, such as Krylov methods. We recall that an effective left preconditioner P ∈ ℝ^N×N for the linear system of Eq. (15) should be a good and cheap approximation of the matrix Id – ΛΣ. More precisely, P⁻¹ (Id – ΛΣ) should have a smaller condition number than Id – ΛΣ, with P⁻¹ being computationally cheap to apply. As pointed out in Sect. 4.4, the application of Σ^II is by far the most expensive part within the application of the whole operator Id – ΛΣ. For this reason, the preconditioner is designed by neglecting the action of Σ^II , namely (16)

where Σ^III is obtained setting α_Q = 0 in Eqs. (4)-(5) when running algorithm 2. The matrix-free version of P is obtained by replacing Σ with Σ^III in algorithm 3. The operator Σ^III is indeed a cheap approximation of Σ and P corresponds to the RT operator in the CRD limit. The matrix-free application of P⁻¹ is described in algorithm 4.

5.2 Parallel Matrix-Free Krylov Solver

For the solution of the linear system in Eq. (15), we propose a matrix-free preconditioned GMRES (PGMRES) method, where the matrix-free actions of Id – ΛΣ and P⁻¹ are given by algorithms 3–4, respectively. In this sense, the proposed solver consists of two nested GMRES iterations; one for the solution of Eq. (15) and one for the solution of the corresponding CRD problem (i.e., for the evaluation of P⁻¹). Numerical tests revealed that other choices of the iterative scheme (e.g., BiCGSTAB method) and of the preconditioner (e.g., Jacobi) are not competitive in terms of run time with the proposed strategy.

We naturally achieve spatial parallelization given the block diagonal structure of Σ, while global communication in space is achieved by the light-weight application of Λ, which is sequential by nature. Parallelization on N_p ≤ N_s processors is achieved with a row-wise partition and distribution of Σ and Λ with a load balancing principle. Spatial blocks are not split between multiple processors; each processor exactly corresponds either to ⌈N_s/N_p⌉ or to ⌈N_s/N_p⌉ spatial nodes, see Fig. 4. The solution vector I is distributed with the same parallel layout. Parallel speed-up is limited by the slowest processor (corresponding to ⌈N_s/N_p⌉ blocks), and the run time with N_p processors, denoted by T(N_p), can be bounded by T(N_p) ≤ T(N_s) • ⌈N_s/N_p⌉.

Fig. 4 Example of row-wise partition of Σ for Ns = 4 and Np = 3, i.e., P1, P2, and P3, aware of the block structure of Σ.

6 Numerical Experiments

6.1 Implementation

For the distributed data structures and Krylov solvers, we used the C++ framework PETSc (see Balay et al. 2015a,b, 1997) with the Cray-MPICH compiler. We stress that the PETSc default row-wise partition follows a load-balancing principle, and except in trivial cases (i.e., when N_s/N_p ∈ ℕ), it does not correspond to the row-wise partition described in Fig. 4. Therefore, the partition must be adjusted by the user, for example, using the PETSc functions VecSetBlockSize( ) and MatSetBlockSize( ) and specifying the block size 4N_θN_χN_v.

The nature of the application of Λ requires sequential calls to the MPI_Send( ) and MPI_Recv( ) routines to propagate the formal solution through the different processors. Since the complexity of the application of Λ is negligible compared to the one of Σ, this sequential bottleneck is not relevant in the presented numerical experiment; see Fig. 3. For larger-scale experiments, where Amdahl’s law can play a role, a multithreaded implementation of Λ with one thread for each ray (θ, χ, ν) is possible.

Numerical experiments have been performed on the Cray XC40 nodes of the Piz Daint supercomputer of the Swiss national supercomputing centre (CSCS)². The used partition features computing nodes with two 18-core Intel Xeon E5-2695v4 (2.10GHz) processors. We used the default GMRES settings of PETSc for the preconditioning described in algorithm 4 and for the solution of the linear system of Eq. (15), using a zero initial guess, a restart after 30 iterations, and a threshold of 10⁻⁵ for the relative residual as a stopping criterion. A finer tuning of the GMRES in algorithm 4, for instance, by increasing its tolerance, is possible.

Fig. 5 Convergence history for the solution of Eq. (15) with Ns = 70, Nθ = 12, Nχ = 8, and Nv = 99; the run time of preconditioned methods for Np = 70 is reported in square brackets.

6.2 Atomic and Atmospheric Models

We considered the Ca I line at 4227 Å, which is an ideal benchmark for new approaches to the modeling of scattering polarization including PRD effects and can be suitably modeled by considering a simple two-level atom (e.g., Faurobert-Scholl 1992; Sampoorna et al. 2009; Supriya et al. 2014; Alsina Ballester et al. 2018; Janett et al. 2021a). We considered the wavelength interval [λ_min, λ_max] = [4220, 4234] Å, discretized with N_v = 99 nodes. The problem was solved in the 1D semi-empirical model C of Fontenla et al. (1993), with [z_min, z_max] = [−100, 2219] km, discretized with N_s = 70 nodes. In the applications shown in Sect. 6.4, the approach was also tested, including the impact of deterministic magnetic and bulk velocity fields.

6.3 Convergence and Scaling

In Table 1, we report the number of iterations to convergence for GMRES and PGMRES, for which we varied the discretization parameters N_θ, N_χ, and N_v, that is, the number of discrete rays. In this regard, both methods prove to be robust; increasing discrete rays does not correspond to an increase of iterations. In Fig. 5, we present an example of the convergence history for preconditioned and unpreconditioned iterative methods. We recall that each BiCGSTAB iteration requires two evaluations of the RT operator, making it globally slower than PGMRES.

In Table 2, we present a strong scaling experiment, that is, we varied the number of processors N_p for a fixed problem size. In Table 3, we present a weak scaling experiment, that is, we varied the number of processors for a fixed problem size per processor. In particular, we report run times and iterations to convergence for the solution of the linear system of Eq. (15) using GMRES and PGMRES. Run times also include the time spent to assemble the right-hand side of Eq. (15), which is negligible, however. We remark that the matrix-free setting does not require other prior assembly costs. The run times from Tables 2 and 3 are also reported in Fig. 6. Both GMRES and PGMRES show an near optimal strong scaling and the proposed PGMRES outperforms the unprecondiotioned GMRES in terms of iterations and run time. Moreover, PGMRES shows an near optimal weak scaling and robustness with respect to the number of spatial points N_s. This is a key finding in view of 3D atmospheric models, where N_s is expected to grow considerably.

Fig. 6 Scaling experiments. Strong-scaling (left panel) and weak-scaling (right panel) experiments, corresponding to the data of Tables 2 and 3, respectively. The dashed lines represent ideal strong and weak scaling.

6.4 Emergent Stokes Profiles

In Fig. 7, we show the impact of the magnetic field on the emergent Stokes profiles of the Ca I 4227 Å line for a line of sight (LOS) with µ = 0.l7. The magnetic field is horizontal and contained in the plane defined by the local vertical and the LOS. The figure clearly shows the Hanle effect, which decreases the amplitude of the line-core peak of the Q/I profile while producing a peak in the core of U/I. The impact of magneto-optical effects is likewise visible, which decreases the wing lobes of Q/I while producing similar wing lobes in U/I. Furthermore, antisymmetric V/I signals are evidently produced by the well-known Zeeman effect. These results excellently agree with those reported in Janett et al. (2021b).

In order to show the impact of bulk velocities on the emergent Stokes profiles, we considered the following vertical bulk velocities:

with the atmospheric height z given in km and v_b in km s⁻¹ . In Fig. 8 we show the emergent Stokes profiles of the Ca I 4227 Å line for two different LOSs, with µ = 0.l7 and µ = 0.38. We first note that the emergent I and Q/I profiles are shifted in frequency due to the Doppler effect; this shift increases for increasing µ values because of the larger bulk velocity component along the LOS. Moreover, the line-core peak signal of the Q/I profiles is significantly enhanced with respect to the zero-velocity case. The radiation field anisotropy is indeed significantly enhanced in the presence of velocity fields with gradients, thus increasing the linear polarization degree (see Carlin et al. 2012). As expected, the amplitude of the whole Q/I signal in the core and in the wings decreases for increasing µ values (i.e., when moving toward the center of the solar disk). We finally note that the Q/I profiles obtained in the presence of bulk velocities are clearly asymmetric, showing a higher blue wing lobe.

7 Conclusions

We considered the NLTE RT problem for polarized radiation, including scattering polarization and angle-dependent PRD effects, which is currently one of the most challenging settings in the field of numerical RT. In this context, after linearization under a suitable physical assumption, we proposed an innovative, efficient, parallel, and scalable solution strategy.

We first presented a suitable discretization of the linear RT problem and its corresponding algebraic formulation in terms of transfer and scattering operators. This formulation was then used to design a matrix-free preconditioned GMRES (PGMRES) solver, where the preconditioner is designed in a multifidelity framework by considering scattering in the light-weight CRD limit.

For a two-level atom and a 1D atmospheric model, numerical experiments show that the proposed parallel PGMRES solver is robust with respect to all the discretization parameters, shows near optimal strong and weak scaling, and is highly competitive in terms of run time. The proposed approach produces accurate solutions in few iterations, with no need to provide a suitable initial guess. In terms of memory footprint, the matrix-free approach only requires the storage of the two vectors I and ϵ, which are distributed among the available processors. This approach thus promises to provide fast and accurate solutions for the modeling of scattering polarization with angle-dependent PRD effects, considering realistic 3D atmospheric models.

We finally remark that algorithms 1–2–3 are general. As long as the emission vector ε is linear with respect to the incident radiation (and the elements of the propagation matrix are known a-priori), the proposed solution strategy is indeed applicable to different discretization techniques (e.g., finite elements), geometries (e.g., 3D spatial domains), atomic models (e.g., two-term atoms), and scattering models (e.g., limit of CRD or angle-averaged PRD).

Fig. 7 Emergent Q/I (left panel), U/I (middle panel), and V/I (right panel) profiles calculated in the FAL-C atmospheric model (Ns = 70), with Νθ = 12, Nχ = 8, and Νν = 99, for the LOS with . A magnetic field of various intensities (0 G, 10 G, and 50 G) and direction (θΒ, χΒ) = (π/2, 0) is considered. The reference direction for positive Q is taken parallel to the x-y plane (see Fig. 1). The PGMRES solver converges in five iterations.

Fig. 8 Emergent I (upper panels) and Q/I (lowerpanels) profiles calculated in the FAL-C atmospheric model (Ns = 70), with Νθ = 12, Nχ = 8, and Νν = 99, for the LOS with (leftpanels) and (rightpanels). The line line center frequency ν0 is reported. A vertical bulk velocity field is considered (see Sect. 6.4). The reference direction for positive Q is taken parallel to the x-y plane (see Fig. 1). The PGMRES solver converges in five iterations.

Acknowledgements

The authors gratefully acknowledge the Swiss National Science Foundation (SNSF) for financial support through grant CRSII5_180238. Rolf Krause acknowledges the funding from the European High-Performance Computing Joint Undertaking (JU) under grant agreement No 955701 (project TIME-X). The JU receives support from the European Union’s Horizon 2020 research and innovation programme and Belgium, France, Germany, Switzerland. We thank the anonymous referee for useful suggestions that helped improving the manuscript.

Appendix A Thermal Contribution to the Line Emissivity

Under the assumption of isotropic inelastic collisions, the thermal contribution to the line emissivity is given by (e.g., Alsina Ballester et al. 2017)

where i = 1 ,…, 4, B is the Planck function (if stimulated emission is neglected, the Wien limit must be considered), and is the generalized profile defined in Landi Degl’Innocenti & Landolfi (2004, App. 13). The factor

represents the branching ratio of the thermal contribution to the total emissivity in Equation (2).

Appendix B Continuum Contributions

Continuum processes contribute to the emissivity through a thermal and a scattering term. These contributions are formally analogous to those due to line processes. Labeling the line and continuum contributions with the indices ℓ and c, respectively, we can write

As far as continuum scattering processes are concerned, it is generally a very good approximation to assume that they are coherent in frequency in the observer frame (e.g., Mihalas 1978; del Pino Alemán et al. 2014). Under this assumption, the continuum scattering term is given by (e.g., Alsina Ballester et al. 2017)

where σ is the continuum opacity for scattering, and is the geometrical tensor introduced in Landi Degl’Innocenti & Landolfi (2004, Sect. 5.11), and defined in the vertical reference system. Finally, the quantity is the radiation field tensor, defined as

The continuum term can thus be written in an integral form analogous to Eq. (3) in terms of the following redistribution matrix

For the conditions typically found in the solar atmosphere, the thermal part of the continuum emissivity is isotropic and only contributes to Stokes I, namely (B.1)

Moreover, continuum processes only contribute to the total opacity (i.e., to the diagonal elements of the propagation matrix), with an isotropic term (B.2)

Routines for the calculation of σ, , and are publicly available (e.g., Uitenbroek 2001).

Appendix C Conversion to Optical Depth

To convert the spatial grid to optical depth, the map

can be considered, which in a 1D geometry is defined as the solution of the IVP (C.1)

with

From the IVP (C.1), we obtain (C.2) (C.3)

and the conversion into the optical depth scale can thus be performed numerically by replacing the integral in Eqs. (C.2) and (C.3) with a suitable numerical quadrature (e.g., Janett & Paganini 2018).

References

All Tables

All Figures

	Fig. 1 Cartesian reference system considered in the problem. The z-axis is directed along the local vertical. The angles θ (inclination) and χ (azimuth) specifying a given propagation direction Ω are indicated in the figure. The considered atmospheric model extends from zmin to zmax. The various physical quantities are constant over the horizontal x-y planes.
In the text

	Fig. 2 Sparsity patterns of Λ, Σ, and Id – ΛΣ using Ns = 10, Nθ = 2, Nχ = 1, and Νν = 10, i.e., N2 = 8002 = 6.4 × 105 elements. The corresponding number of nonzero elements (nz) is reported. The elements in the lower and upper triangular parts of Id – ΛΣ correspond to downward (negative µ) and upward (positive µ) directions, respectively.
In the text

	Fig. 3 Run times for a single application of Σ, ΣIII, and Λ using algorithms 1–2 varying the number of processors Np ≤ Ns with Ns = 70, Nθ = 12, Nχ = 8, and Nv = 99.
In the text

	Fig. 4 Example of row-wise partition of Σ for Ns = 4 and Np = 3, i.e., P1, P2, and P3, aware of the block structure of Σ.
In the text

	Fig. 5 Convergence history for the solution of Eq. (15) with Ns = 70, Nθ = 12, Nχ = 8, and Nv = 99; the run time of preconditioned methods for Np = 70 is reported in square brackets.
In the text

	Fig. 6 Scaling experiments. Strong-scaling (left panel) and weak-scaling (right panel) experiments, corresponding to the data of Tables 2 and 3, respectively. The dashed lines represent ideal strong and weak scaling.
In the text

Fig. 7 Emergent Q/I (left panel), U/I (middle panel), and V/I (right panel) profiles calculated in the FAL-C atmospheric model (Ns = 70), with Νθ = 12, Nχ = 8, and Νν = 99, for the LOS with . A magnetic field of various intensities (0 G, 10 G, and 50 G) and direction (θΒ, χΒ) = (π/2, 0) is considered. The reference direction for positive Q is taken parallel to the x-y plane (see Fig. 1). The PGMRES solver converges in five iterations.

In the text

Fig. 8 Emergent I (upper panels) and Q/I (lowerpanels) profiles calculated in the FAL-C atmospheric model (Ns = 70), with Νθ = 12, Nχ = 8, and Νν = 99, for the LOS with (leftpanels) and (rightpanels). The line line center frequency ν0 is reported. A vertical bulk velocity field is considered (see Sect. 6.4). The reference direction for positive Q is taken parallel to the x-y plane (see Fig. 1). The PGMRES solver converges in five iterations.

In the text