Issue 
A&A
Volume 664, August 2022



Article Number  A197  
Number of page(s)  12  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/202243059  
Published online  31 August 2022 
Numerical solutions to linear transfer problems of polarized radiation
III. Parallel preconditioned Krylov solver tailored for modeling PRD effects
^{1}
Euler Institute, Università della Svizzera italiana (USI),
6900
Lugano, Switzerland
^{2}
IRSOL Istituto Ricerche Solari “Aldo e Cele Daccò”, Università della Svizzera italiana,
6605
LocarnoMonti, Switzerland
email: gioele.janett@irsol.usi.ch
^{3}
LeibnizInstitut für Sonnenphysik (KIS),
79104
Freiburg i. Br., Germany
Received:
7
January
2022
Accepted:
1
May
2022
Context. The polarization signals produced by the scattering of anistropic radiation in strong resonance lines encode important information about the elusive magnetic fields in the outer layers of the solar atmosphere. An accurate modeling of these signals is a very challenging problem from the computational point of view, in particular when partial frequency redistribution (PRD) effects in scattering processes are accounted for with a general angledependent treatment.
Aims. We aim at solving the radiative transfer problem for polarized radiation in nonlocal thermodynamic equilibrium conditions, taking angledependent PRD effects into account. The problem is formulated for a twolevel atomic model in the presence of arbitrary magnetic and bulk velocity fields. The polarization produced by scattering processes and the Zeeman effect is considered.
Methods. The proposed solution strategy is based on an algebraic formulation of the problem and relies on a convenient physical assumption, which allows its linearization. We applied a nested matrixfree GMRES iterative method. Effective preconditioning is obtained in a multifidelity framework by considering the lightweight description of scattering processes in the limit of complete frequency redistribution (CRD).
Results. Numerical experiments for a onedimensional (1D) atmospheric model show near optimal strong and weak scaling of the proposed CRDpreconditioned GMRES method, which converges in few iterations, independently of the discretization parameters. A suitable parallelization strategy and highperformance computing tools lead to competitive run times, providing accurate solutions in a few minutes.
Conclusions. The proposed solution strategy allows the fast systematic modeling of the scattering polarization signals of strong resonance lines, taking angledependent PRD effects into account together with the impact of arbitrary magnetic and bulk velocity fields. Almost optimal strong and weak scaling results suggest that this strategy is applicable to realistic 3D models. Moreover, the proposed strategy is general, and applications to more complex atomic models are possible.
Key words: radiative transfer / polarization / scattering / stars: atmospheres / Sun: atmosphere
© P. Benedusi et al. 2022
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the SubscribetoOpen model. Subscribe to A&A to support open access publication.
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_{2}. 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 magnetooptical 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 magnetooptical 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 wellknown Zeeman effect. The possibility of exploiting the combined action of the Hanle, magnetooptical, 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., FaurobertScholl 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 angledependent treatment of PRD effects, without introducing simplifying assumptions, such as the angleaveraged 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 angledependent 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 lastgeneration 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 angledependent PRD problem for a twolevel atom in Hanle and HanleZeeman regimes, respectively. Sampoorna & Nagendra (2015) developed an angledependent PRD RT code capable of including vertical bulk velocities. Supriya et al. (2013) solved the angledependent 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 angledependent PRD effects, considering a threeterm atom. For computational simplicity, the angledependent 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 angledependent 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 matrixfree 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 fourcomponent 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 integrodifferential 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 twolevel 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 frequencyintegrated absorption coefficient and R_{ij} are the elements of the socalled 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 lightweight models for the description of scattering (e.g., the assumption of CRD) or by the use of simplifying approximations (e.g., angleaveraging). 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 secondorder process in a perturbative expansion of matterradiation interaction, under given circumstances, it can be suitably described as a temporal succession of independent firstorder absorption and reemission 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 flatspectrum 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 highorder perturbative expansion of matterradiation interaction and selfconsistently includes frequency redistribution effects due to elastic collisions. Finally, we recall the approach based on a diagrammatic treatment of atomphoton 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 angleaveraged 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 twolevel 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 angledependent 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 deexcitation, 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^{1}. 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)} = 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 emissivity is dominated by the collisional term; this is the limit of local thermodynamic equilibrium (LTE), which is a lightweight 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 zaxis 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 z_{min} to z_{max}. The various physical quantities are constant over the horizontal xy 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 PlaneParallel Atmospheric Model
We considered a semiempirical planeparallel 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 ingoing 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 semidiscrete 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 GaussLegendre 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 GaussLegendre 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 semiempirical 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 semidiscrete 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 shortcharacteristic 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 η_{1} is the diagonal element of the propagation matrix and
For a given ray (θ, χ, ν), the exact integration of Eq. (9) in the interval [z_{0}, z] yields the variation of constants formula (10)
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 DELOlinear method, which is an Lstable 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 DELOlinear 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_{1}. 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., DELOlinear) 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 matrixfree 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 singlesize N linear system with unknown I, namely (15)
The righthand 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 Mmatrix 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 fixedpoint 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 blockJacobi 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^{6} for 1D models and reaches 10^{10} for full 3D models. Crucially, Krylov methods are well suited for matrixfree approaches and have proved to be highly effective solution strategies for large linear systems.
Fig. 2 Sparsity patterns of Λ, Σ, and Id – ΛΣ using N_{s} = 10, N_{θ} = 2, N_{χ} = 1, and Ν_{ν} = 10, i.e., N^{2} = 800^{2} = 6.4 × 10^{5} 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 N_{p} ≤ N_{s} with N_{s} = 70, N_{θ} = 12, N_{χ} = 8, and N_{v} = 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^{−1} (Id – ΛΣ) should have a smaller condition number than Id – ΛΣ, with P^{−1} 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 matrixfree 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 matrixfree application of P^{−1} is described in algorithm 4.
5.2 Parallel MatrixFree Krylov Solver
For the solution of the linear system in Eq. (15), we propose a matrixfree preconditioned GMRES (PGMRES) method, where the matrixfree actions of Id – ΛΣ and P^{−1} 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^{−1}). 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 lightweight application of Λ, which is sequential by nature. Parallelization on N_{p} ≤ N_{s} processors is achieved with a rowwise 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 speedup 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 rowwise partition of Σ for N_{s} = 4 and N_{p} = 3, i.e., P_{1}, P_{2}, and P_{3}, 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 CrayMPICH compiler. We stress that the PETSc default rowwise partition follows a loadbalancing principle, and except in trivial cases (i.e., when N_{s}/N_{p} ∈ ℕ), it does not correspond to the rowwise 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 largerscale 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)^{2}. The used partition features computing nodes with two 18core Intel Xeon E52695v4 (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^{−5} 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 N_{s} = 70, N_{θ} = 12, N_{χ} = 8, and N_{v} = 99; the run time of preconditioned methods for N_{p} = 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 twolevel atom (e.g., FaurobertScholl 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 semiempirical 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 righthand side of Eq. (15), which is negligible, however. We remark that the matrixfree 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. Strongscaling (left panel) and weakscaling (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 linecore peak of the Q/I profile while producing a peak in the core of U/I. The impact of magnetooptical 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 wellknown 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^{−1} . 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 linecore peak signal of the Q/I profiles is significantly enhanced with respect to the zerovelocity 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 angledependent 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 matrixfree preconditioned GMRES (PGMRES) solver, where the preconditioner is designed in a multifidelity framework by considering scattering in the lightweight CRD limit.
For a twolevel 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 matrixfree 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 angledependent 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 apriori), 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., twoterm atoms), and scattering models (e.g., limit of CRD or angleaveraged PRD).
Fig. 7 Emergent Q/I (left panel), U/I (middle panel), and V/I (right panel) profiles calculated in the FALC atmospheric model (N_{s} = 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 xy 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 FALC atmospheric model (N_{s} = 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 xy 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 HighPerformance Computing Joint Undertaking (JU) under grant agreement No 955701 (project TIMEX). 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)
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
 Alsina Ballester, E., Belluzzi, L., & Trujillo Bueno, J. 2016, ApJ, 831, L15 [NASA ADS] [CrossRef] [Google Scholar]
 Alsina Ballester, E., Belluzzi, L., & Trujillo Bueno, J. 2017, ApJ, 836, 6 [NASA ADS] [CrossRef] [Google Scholar]
 Alsina Ballester, E., Belluzzi, L., & Trujillo Bueno, J. 2018, ApJ, 854, 150 [NASA ADS] [CrossRef] [Google Scholar]
 Anusha, L.S., & Nagendra, K.N. 2011, ApJ, 739, 40 [NASA ADS] [CrossRef] [Google Scholar]
 Anusha, L.S., & Nagendra, K.N. 2012, ApJ, 746, 84 [NASA ADS] [CrossRef] [Google Scholar]
 Anusha, L., Nagendra, K., Paletou, F., & Léger, L. 2009, ApJ, 704, 661 [CrossRef] [Google Scholar]
 Badri, M., Jolivet, P., Rousseau, B., & Favennec, Y. 2018, J. Comput. Phys., 360, 74 [NASA ADS] [CrossRef] [Google Scholar]
 Badri, M., Jolivet, P., Rousseau, B., & Favennec, Y. 2019, Comput. Math. Applic., 77, 1453 [CrossRef] [Google Scholar]
 Balay, S., Gropp, W.D., McInnes, L.C., & Smith, B.F. 1997, in Modern Software Tools in Scientific Computing, eds. E. Arge, A.M. Bruaset, & H.P. Langtangen (Birkhäuser Press), 163 [CrossRef] [Google Scholar]
 Balay, S., Abhyankar, S., Adams, M.F., et al. 2015a, PETSc Users Manual, Tech. Rep. ANL95/11  Revision 3.6, Argonne National Laboratory [Google Scholar]
 Balay, S., Abhyankar, S., Adams, M.F., et al. 2015b, PETSc Web page, http://www.mcs.anl.gov/petsc [Google Scholar]
 Belluzzi, L., & Trujillo Bueno, J. 2012, ApJ, 750, L11 [NASA ADS] [CrossRef] [Google Scholar]
 Belluzzi, L., & Trujillo Bueno, J. 2014, A&A, 564, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Belluzzi, L., Trujillo Bueno, J., & Stepân, J. 2012, ApJ, 755, L2 [NASA ADS] [CrossRef] [Google Scholar]
 Benedusi, P., Janett, G., Belluzzi, L., & Krause, R. 2021, A&A, 655, A88 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bommier, V. 1997a, A&A, 328, 706 [NASA ADS] [Google Scholar]
 Bommier, V. 1997b, A&A, 328, 726 [NASA ADS] [Google Scholar]
 Bommier, V. 2017, A&A, 607, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Carlin, E.S., Manso Sainz, R., Asensio Ramos, A., & Trujillo Bueno, J. 2012, ApJ, 751, 5 [NASA ADS] [CrossRef] [Google Scholar]
 Carlin, E.S., Asensio Ramos, A., & Trujillo Bueno, J. 2013, ApJ, 764, 40 [NASA ADS] [CrossRef] [Google Scholar]
 Casini, R., & Landi Degl'Innocenti, E. 2008, Astrophysical Plasmas, 44, in eds. T. Fujimoto, & A. Iwamae, Atsushi, 247 [NASA ADS] [Google Scholar]
 Casini, R., Landi Degl'Innocenti, M., Manso Sainz, R., Land i Degl'Innocenti, E., & Landolfi, M. 2014, ApJ, 791, 94 [NASA ADS] [CrossRef] [Google Scholar]
 Casini, R., del Pino Alemân, T., & Manso Sainz, R. 2017a, ApJ, 835, 114 [NASA ADS] [CrossRef] [Google Scholar]
 Casini, R., del Pino Alemân, T., & Manso Sainz, R. 2017b, ApJ, 848, 99 [NASA ADS] [CrossRef] [Google Scholar]
 del Pino Alemân, T., & Trujillo Bueno, J. 2021, ApJ, 909, 180 [CrossRef] [Google Scholar]
 del Pino Alemân, T., Manso Sainz, R., & Trujillo Bueno, J. 2014, ApJ, 784, 46 [CrossRef] [Google Scholar]
 del Pino Alemân, T., Casini, R., & Manso Sainz, R. 2016, ApJ, 830, L24 [CrossRef] [Google Scholar]
 del Pino Alemân, T., Trujillo Bueno, J., Casini, R., & Manso Sainz, R. 2020, ApJ, 891, 91 [CrossRef] [Google Scholar]
 Fabiani Bendicho, P., Trujillo Bueno, J., & Auer, L. 1997, A&A, 324, 161 [NASA ADS] [Google Scholar]
 FaurobertScholl, M. 1992, A&A, 258, 521 [NASA ADS] [Google Scholar]
 Fontenla, J.M., Avrett, E.H., & Loeser, R. 1993, ApJ, 406, 319 [NASA ADS] [CrossRef] [Google Scholar]
 Gandorfer, A. 2000, The Second Solar Spectrum: A High Spectral Resolution Polarimetric Survey of Scattering Polarization at the Solar Limb in Graphical Representation, Volume I: 4625 A6995 A (Zurich: vdf ETH) [NASA ADS] [Google Scholar]
 Gandorfer, A. 2002, The Second Solar Spectrum: A High Spectral Resolution Polarimetric Survey of Scattering Polarization at the Solar Limb in Graphical Representation. Volume II: 3910 A4630 A (Zurich: vdf ETH) [NASA ADS] [Google Scholar]
 Guderley, G., & Hsu, C.C. 1972, Math. Comp., 26, 51 [CrossRef] [Google Scholar]
 Holzreuter, R., Fluri, D.M., & Stenflo, J.O. 2005, A&A, 434, 713 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hummer, D.G. 1962, MNRAS, 125, 21 [NASA ADS] [CrossRef] [Google Scholar]
 Ishikawa, R., Bueno, J.T., del Pino Alemân, T., et al. 2021, Sci. Adv., 7, eabe8406 [NASA ADS] [CrossRef] [Google Scholar]
 Janett, G., & Paganini, A. 2018, ApJ, 857, 91 [NASA ADS] [CrossRef] [Google Scholar]
 Janett, G., Carlin, E.S., Steiner, O., & Belluzzi, L. 2017a, ApJ, 840, 107 [NASA ADS] [CrossRef] [Google Scholar]
 Janett, G., Steiner, O., & Belluzzi, L. 2017b, ApJ, 845, 104 [NASA ADS] [CrossRef] [Google Scholar]
 Janett, G., Steiner, O., & Belluzzi, L. 2018, ApJ, 865, 16 [NASA ADS] [CrossRef] [Google Scholar]
 Janett, G., Alsina Ballester, E., Guerreiro, N., et al. 2021a, A&A, 655, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Janett, G., Benedusi, P., Belluzzi, L., & Krause, R. 2021b, A&A, 655, A87 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kano, R., Trujillo Bueno, J., Winebarger, A., et al. 2017, ApJ, 839, L10 [NASA ADS] [CrossRef] [Google Scholar]
 Landi Degl'Innocenti, E., & Landolfi, M. 2004, Polarization in Spectral Lines (Dordrecht: Kluwer Academic Publishers) Astrophysics and Space Science Library, 307 [Google Scholar]
 Landi Degl'Innocenti, E., Landi Degl'Innocenti, M., & Landolfi, M. 1997, in Forum THÉMIS, Science with THÉMIS, eds. N. Mein, & S. SahalBréchot (Paris: Obs. de Paris), 59 [Google Scholar]
 Leenaarts, J., Pereira, T., & Uitenbroek, H. 2012, A&A, 543, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mihalas, D. 1978, Stellar Atmospheres, 2nd edn. (San Francisco: W.H. Freeman and Company) [Google Scholar]
 Nagendra, K.N. 2019, in Astronomical Society of the Pacific Conference Series, eds. L. Belluzzi, R. Casini, M. Romoli, & J. Trujillo Bueno, 526, 99 [NASA ADS] [Google Scholar]
 Nagendra, K.N., Frisch, H., & Faurobert, M. 2002, A&A, 395, 305 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Nagendra, K.N., Sowmya, K., Sampoorna, M., Stenflo, J.O., & Anusha, L.S. 2020, ApJ, 898, 49 [NASA ADS] [CrossRef] [Google Scholar]
 Rees, D.E., & Saliba, G.J. 1982, A&A, 115, 1 [NASA ADS] [Google Scholar]
 Rees, D., Murphy, G., & Durrant, C. 1989, ApJ, 339, 1093 [NASA ADS] [CrossRef] [Google Scholar]
 Ren, K., Zhang, R., & Zhong, Y. 2019, J. Comput. Phys., 399, 108958 [NASA ADS] [CrossRef] [Google Scholar]
 Sampoorna, M. 2011, ApJ, 731, 114 [NASA ADS] [CrossRef] [Google Scholar]
 Sampoorna, M., & Nagendra, K.N. 2015, ApJ, 812, 28 [NASA ADS] [CrossRef] [Google Scholar]
 Sampoorna, M., Nagendra, K.N., & Stenflo, J.O. 2008, ApJ, 679, 889 [NASA ADS] [CrossRef] [Google Scholar]
 Sampoorna, M., Stenflo, J.O., Nagendra, K.N., et al. 2009, ApJ, 699, 1650 [NASA ADS] [CrossRef] [Google Scholar]
 Sampoorna, M., Nagendra, K.N., & Stenflo, J.O. 2017, ApJ, 844, 97 [NASA ADS] [CrossRef] [Google Scholar]
 Stenflo, J. 1994, Solar Magnetic Fields: Polarized Radiation Diagnostics, 189 (Springer) [Google Scholar]
 Štěpán, J., Trujillo Bueno, J., Leenaarts, J., & Carlsson, M. 2015, ApJ, 803, 65 [CrossRef] [Google Scholar]
 Supriya, H.D., Smitha, H.N., Nagendra, K.N., Ravindra, B., & Sampoorna, M. 2013, MNRAS, 429, 275 [NASA ADS] [CrossRef] [Google Scholar]
 Supriya, H.D., Smitha, H.N., Nagendra, K.N., et al. 2014, ApJ, 793, 42 [NASA ADS] [CrossRef] [Google Scholar]
 Trujillo Bueno, J. 2014, in Solar Polarization 7, eds. K.N. Nagendra, J.O. Stenflo, Z.Q. Qu, & M. Sampoorna, Astronomical Society of the Pacific Conference Series, 489, 137 [Google Scholar]
 Trujillo Bueno, J., & Manso Sainz, R. 1999, ApJ, 516, 436 [CrossRef] [Google Scholar]
 Trujillo Bueno, J., Landi Degl'Innocenti, E., & Belluzzi, L. 2017, Space Sci. Rev., 210, 183 [NASA ADS] [CrossRef] [Google Scholar]
 Trujillo Bueno, J., Stepân, J., Belluzzi, L., et al. 2018, ApJ, 866, L15 [NASA ADS] [CrossRef] [Google Scholar]
 Uitenbroek, H. 2001, ApJ, 557, 389 [NASA ADS] [CrossRef] [Google Scholar]
 Warsa, J., Thompson, K., & Morel, J. 2003, Trans. Am. Nuclear Soci., 449 [Google Scholar]
 Weisskopf, V., & Wigner, E. 1930, Z. Phys., 63, 54 [CrossRef] [Google Scholar]
It is interesting to observe that if we artificially set α_{Q} = 0 while keeping the nominal value of , and when we assume that the incident field is spectrally flat, the ensuing emission coefficient coincides with the one derived in the CRD approach of Landi Degl’Innocenti & Landolfi (2004).
All Tables
All Figures
Fig. 1 Cartesian reference system considered in the problem. The zaxis 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 z_{min} to z_{max}. The various physical quantities are constant over the horizontal xy planes. 

In the text 
Fig. 2 Sparsity patterns of Λ, Σ, and Id – ΛΣ using N_{s} = 10, N_{θ} = 2, N_{χ} = 1, and Ν_{ν} = 10, i.e., N^{2} = 800^{2} = 6.4 × 10^{5} 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 N_{p} ≤ N_{s} with N_{s} = 70, N_{θ} = 12, N_{χ} = 8, and N_{v} = 99. 

In the text 
Fig. 4 Example of rowwise partition of Σ for N_{s} = 4 and N_{p} = 3, i.e., P_{1}, P_{2}, and P_{3}, aware of the block structure of Σ. 

In the text 
Fig. 5 Convergence history for the solution of Eq. (15) with N_{s} = 70, N_{θ} = 12, N_{χ} = 8, and N_{v} = 99; the run time of preconditioned methods for N_{p} = 70 is reported in square brackets. 

In the text 
Fig. 6 Scaling experiments. Strongscaling (left panel) and weakscaling (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 FALC atmospheric model (N_{s} = 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 xy 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 FALC atmospheric model (N_{s} = 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 xy plane (see Fig. 1). The PGMRES solver converges in five iterations. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.