Open Access
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/0004-6361/202243059
Published online 31 August 2022

© P. Benedusi et al. 2022

Licence Creative CommonsOpen 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 Subscribe-to-Open 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 D2. 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 I=(I,Q,U,V)T=(I1,I2,I3,I4)T,$ {\bf{I}} = {\left( {I,Q,U,V} \right)^T} = {\left( {{I_1},{I_2},{I_3},{I_4}} \right)^T}, $

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 ΩIi(r,Ω,v)=j=14Kij(r,Ω,v)+εi(r,Ω,v),$ {\bf{\Omega }} \cdot \nabla {I_i}\left( {{\bf{r}},{\bf{\Omega }},v} \right) = - \sum\limits_{j = 1}^4 {{K_{ij}}\left( {{\bf{r}},{\bf{\Omega }},v} \right) + {\varepsilon _i}\left( {r,{\bf{\Omega }},v} \right)} , $(1)

where i = 1,…, 4, rD, Ω = (θ,χ) ∈ [0, π] × [0, 2π) is a unit vector, and ν ∈ [vmin, vmax] ⊂ ℝ+. Equation (1) is equipped with suitable boundary conditions on ∂D, and describes a set of initial values problems (IVPs). The 4 × 4 propagation matrix K=(η1η2η3η4η2η1ρ4ρ3η3ρ4η1ρ2η4ρ3ρ2η1)$ K = \left( {\matrix{ {{\eta _1}} \hfill & {{\eta _2}} \hfill & {{\eta _3}} \hfill & {{\eta _4}} \hfill \cr {{\eta _2}} \hfill & {{\eta _1}} \hfill & {{\rho _4}} \hfill & { - {\rho _3}} \hfill \cr {{\eta _3}} \hfill & { - {\rho _4}} \hfill & {{\eta _1}} \hfill & {{\rho _2}} \hfill \cr {{\eta _4}} \hfill & {{\rho _3}} \hfill & { - {\rho _2}} \hfill & {{\eta _1}} \hfill \cr } } \right) $

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 ε=(ε1,ε2,ε3,ε4)T$ \varepsilon = {\left( {{\varepsilon _1},{\varepsilon _2},{\varepsilon _3},{\varepsilon _4}} \right)^T} $

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, ε(r,Ω,v)+εth(r,Ω,v)+εsc(r,Ω,v).$ {\bf{\varepsilon }}\left( {{\bf{r}},{\bf{\Omega }},v} \right) + {{\bf{\varepsilon }}^{{\rm{th}}}}\left( {{\bf{r}},{\bf{\Omega }},v} \right) + {{\bf{\varepsilon }}^{{\rm{sc}}}}\left( {{\bf{r}},{\bf{\Omega }},v} \right). $(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: εisc(r,Ω,v)=kL(r)=dvdΩ4πj=14Rij(r,Ω,Ω,v,v)Ij(r,Ω,v),$ \varepsilon _i^{{\rm{sc}}}\left( {{\bf{r}},{\bf{\Omega }},v} \right) = {k_L}\left( {\bf{r}} \right) = \smallint {\rm{d}}v'\oint {{{{\rm{d}}\Omega '} \over {4\pi }}\sum\limits_{j = 1}^4 {{R_{ij}}\left( {{\bf{r}},{\bf{\Omega '}},{\bf{\Omega }},v',v} \right){I_j}} \left( {{\bf{r}},{\bf{\Omega '}},v'} \right)} , $(3)

where kL is the frequency-integrated absorption coefficient and Rij 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) Rij(r,Ω,Ω,v,v)=RijII(r,Ω,Ω,v,v)+RijII(r,Ω,Ω,v,v)$ {R_{ij}}\left( {{\bf{r}},{\bf{\Omega '}},{\bf{\Omega }},v',v} \right) = R_{ij}^{{\rm{II}}}\left( {{\bf{r}},{\bf{\Omega '}},{\bf{\Omega }},v',v} \right) + R_{ij}^{{\rm{II}}}\left( {{\bf{r}},{\bf{\Omega '}},{\bf{\Omega }},v',v} \right) $

where RII describes scattering processes that are coherent in frequency in the atomic rest frame, while RIII 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 RII and RIII 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 RijII$R_{ij}^{{\rm{II}}}$ (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 RijIII$R_{ij}^{{\rm{III}}}$ (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 RijII$R_{ij}^{{\rm{II}}}$ in the observer frame and the approximation of CRD in the observer frame for RijIII$R_{ij}^{{\rm{III}}}$, 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 RijII$R_{ij}^{{\rm{II}}}$ and RijIII$R_{ij}^{{\rm{III}}}$ redistribution matrices can be written as Rij(r,Ω,Ω,v,v)=KKQαQ(r)rQII,KK(r,Ω,Ω,v,v)^Q,iK(Ω)^Q,iK(Ω),$ \matrix{ {{R_{ij}}\left( {{\bf{r}},{\bf{\Omega '}},{\bf{\Omega }},v',v} \right) = } \hfill & {\sum\limits_{KK'Q} {{\alpha _Q}\left( {\bf{r}} \right)r_Q^{{\rm{II,}}KK'}\left( {{\bf{r}},{\bf{\Omega '}},{\bf{\Omega }},v',v} \right)} } \hfill \cr {} \hfill & { \cdot \hat \Im _{Q,i}^{K'}\left( {\bf{\Omega }} \right)\hat \Im _{ - Q,i}^{K'}\left( {\bf{\Omega }} \right),} \hfill \cr } $(4) Rijm(r,Ω,Ω,v,v)=KKKQ(βQK((r)αQ(r))rQIII,KKK(r,Ω,Ω,v,v)^Q,iK(Ω)^Q,iK(Ω),$ \matrix{ {R_{ij}^m\left( {{\bf{r}},{\bf{\Omega '}},{\bf{\Omega }},v',v} \right) = } \hfill & {\sum\limits_{KK'K''Q} {(\beta _Q^{K''}\left( {\left( {\bf{r}} \right) - {\alpha _Q}\left( {\bf{r}} \right)} \right)r_Q^{{\rm{III,}}KK'K''}\left( {{\bf{r}},{\bf{\Omega '}},{\bf{\Omega }},v',v} \right)} } \hfill \cr {} \hfill & { \cdot \hat \Im _{Q,i}^{K'}\left( {\bf{\Omega }} \right)\hat \Im _{ - Q,i}^K\left( {{\bf{\Omega '}}} \right),} \hfill \cr } $(5)

where JQ,iK${\cal J}_{Q,i}^K$ 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 rQII,KK$r_Q^{{\rm{II,}}KK\prime }$ and rQIII,KKK$r_Q^{{\rm{III,}}KK\prime K\prime \prime }$ 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 rQII,KK$r_Q^{{\rm{II,}}KK\prime }$ and rQIII,KKK$r_Q^{{\rm{III,}}KK\prime K\prime \prime }$ expressions. The quantities αQ and (βQKαQ)$\left( {\beta _Q^K - {\alpha _Q}} \right)$ are the branching ratios for RijII$R_{ij}^{{\rm{II}}}$ and RijIII$R_{ij}^{{\rm{III}}}$, respectively, with αQ(r)=ΓRΓR+Γl(r)ΓR+Γl(r)ΓR+Γl(r)+ΓE(r)+2πivL(r)guQ,βQK(r)=ΓRΓR+Γl(r)ΓR+Γl(r)ΓR+Γl(r)+D(K)(r)+2πivL(r)guQ,$ \eqalign{ & {\alpha _Q}\left( {\bf{r}} \right) = {{{\Gamma _R}} \over {{\Gamma _R} + {\Gamma _l}\left( {\bf{r}} \right)}}{{{\Gamma _R} + {\Gamma _l}\left( {\bf{r}} \right)} \over {{\Gamma _R} + {\Gamma _l}\left( {\bf{r}} \right) + {\Gamma _E}\left( {\bf{r}} \right) + 2\pi {\rm{i}}{v_L}\left( {\bf{r}} \right){g_u}Q}}, \cr & \beta _Q^K\left( {\bf{r}} \right) = {{{\Gamma _R}} \over {{\Gamma _R} + {\Gamma _l}\left( {\bf{r}} \right)}}{{{\Gamma _R} + {\Gamma _l}\left( {\bf{r}} \right)} \over {{\Gamma _R} + {\Gamma _l}\left( {\bf{r}} \right) + {D^{\left( K \right)}}\left( {\bf{r}} \right) + 2\pi {\rm{i}}{v_L}\left( {\bf{r}} \right){g_u}Q}}, \cr} $

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, vL is the Larmor frequency, and gu 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 RijIII$R_{ij}^{{\rm{III}}}$ vanishes and the total redistribution matrix reduces to RijII$R_{ij}^{{\rm{II}}}$; 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 RijII$R_{ij}^{{\rm{II}}}$ vanishes and the total redistribution matrix reduces to RijIII$R_{ij}^{{\rm{III}}}$; this is the limit of CRD in the atomic frame1. 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 ΓRΓR+Γl(r)$ {{{\Gamma _R}} \over {{\Gamma _R} + {\Gamma _l}\left( {\bf{r}} \right)}} $

is common to both αQ and βQK${\beta _Q^K}$ 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.

thumbnail 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 kL. 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 kL 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 ∈ [zmm, zmax] ⊂ ℝ, see Fig. 1. With this simplification, the RT Eq. (1) can be rewritten as cos(θ)ddzI(z,θ,χ,v)=K(z,θ,χ,v)I(z,θ,χ,v)+ε(z,θ,χ,v).$ \cos \left( \theta \right){{\rm{d}} \over {{\rm{d}}z}}{\bf{I}}\left( {z,\theta ,\chi ,v} \right) = - {\bf{K}}\left( {z,\theta ,\chi ,v} \right){\bf{I}}\left( {z,\theta ,\chi ,v} \right) + \varepsilon \left( {z,\theta ,\chi ,v} \right). $(6)

For in-going rays at the domain boundary ∂D, Eq. (6) is equipped with the following boundary conditions: I(zmin,θ,χ,v)=Iin(v)forθ[0,π/2),χ,v,I(zmax,θ,χ,v)=0forθ(π/2,π),χ,v.$ \matrix{ {{\bf{I}}\left( {{z_{\min }},\theta ,\chi ,v} \right) = {{\bf{I}}_{{\rm{in}}}}\left( v \right)} \hfill & {} \hfill & {} \hfill & {{\rm{for}}} \hfill & {\theta \in [0,\pi /2),{\forall _\chi },{\forall _v},} \hfill \cr {{\bf{I}}\left( {{z_{\max }},\theta ,\chi ,v} \right) = 0} \hfill & {} \hfill & {} \hfill & {{\rm{for}}} \hfill & {\theta \in (\pi /2,\pi ),{\forall _\chi },{\forall _v}.} \hfill \cr } $(7)

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

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 [vmin, vmax] is discretized with a unevenly spaced grid with Nv nodes, namely vmin=v1<v2<<vNv=vmax.$ {v_{\min }} = {v_1} < {v_2} < \cdots < {v_{{N_v}}} = {v_{\max }}. $

The frequency grid { vn }n=1Nv$\left\{ {{v_n}} \right\}_{n = 1}^{{N_v}}$ 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 1<μ1<μ2<<μNθ/2<0<μNθ/2+1<<μNθ<1,$ - 1 < {\mu _1} < {\mu _2} < \cdots < {\mu _{{N_\theta }/2}} < 0 < {\mu _{{N_\theta }/2 + 1}} < \cdots < {\mu _{{N_\theta }}} < 1, $

for µ ∈ (−1, 0) and µ ∈ (0, 1), respectively. This grid corresponds to 0<θ1<θ2<<θNθ/2<π/2<θNθ/2+1<<θNθ<π,$ 0 < {\theta _1} < {\theta _2} < \cdots < {\theta _{{N_\theta }/2}} < \pi /2 < {\theta _{{N_\theta }/2 + 1}} < \cdots < {\theta _{{N_\theta }}} < \pi , $

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 χk=k2π/Nχ    for   k=1,,Nχ.$ {\chi _k} = k \cdot 2\pi /{N_\chi }\,\,\,\,{\rm{for}}\,\,\,k = 1, \ldots ,{N_\chi }. $

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 Ns nodes, namely zmin=z1<z2<<zNs=zmax.$ {z_{\min }} = {z_1} < {z_2} < \cdots < {z_{{N_{\rm{s}}}}} = {z_{\max }}. $(8)

For (θj, χk, vn) and a corresponding emissivity vector ε at all { zl }l=1Ns$\left\{ {{z_l}} \right\}_{l = 1}^{{N_{\rm{s}}}}$, the solution of the semi-discrete ODEs system (6) is numerically approximated in { zl }l=1Ns$\left\{ {{z_l}} \right\}_{l = 1}^{{N_{\rm{s}}}}$ 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 cos(θ)ddzI(z,θ,χ,v)=η1(z,θ,χ,v)[ S(z,θ,χ,v)I(z,θ,χ,v) ],$ \cos \left( \theta \right){{\rm{d}} \over {{\rm{d}}z}}{\bf{I}}\left( {z,\theta ,\chi ,v} \right) = {\eta _1}\left( {z,\theta ,\chi ,v} \right)\left[ {{\bf{S}}\left( {z,\theta ,\chi ,v} \right) - {\bf{I}}\left( {z,\theta ,\chi ,v} \right)} \right], $(9)

where η1 is the diagonal element of the propagation matrix and S(z,θ,χ,v)=(IDK(z,θ,χ,v)η1(z,θ,χ,v))I(z,θ,χ,v)+ε(z,θ,χ,v)η1(z,θ,χ,v).$ {\bf{S}}\left( {z,\theta ,\chi ,v} \right) = \left( {ID - {{K\left( {z,\theta ,\chi ,v} \right)} \over {{\eta _1}\left( {z,\theta ,\chi ,v} \right)}}} \right){\bf{I}}\left( {z,\theta ,\chi ,v} \right) + {{\varepsilon \left( {z,\theta ,\chi ,v} \right)} \over {{\eta _1}\left( {z,\theta ,\chi ,v} \right)}}. $

For a given ray (θ, χ, ν), the exact integration of Eq. (9) in the interval [z0, z] yields the variation of constants formula I(z,θ,χ,v)=eτ(z,θ,v)I(z0,θ,χ,v)+z0zdxeτ(x,θ,χ,v)S(x,θ,χ,v),$ {\bf{I}}\left( {z,\theta ,\chi ,v} \right) = {e^{ - \tau \left( {z,\theta ,v} \right)}}{\bf{I}}\left( {{z_0},\theta ,\chi ,v} \right) + \int_{{z_0}}^z {{\rm{d}}x\,{e^{ - \tau \left( {x,\theta ,\chi ,v} \right)}}} {\bf{S}}\left( {x,\theta ,\chi ,v} \right), $(10)

where τ(z,θ,χ,v)=1| cos(θ) |z0zdxη1(x,θ,χ,v).$ \tau \left( {z,\theta ,\chi ,v} \right) = {1 \over {\left| {\cos \left( \theta \right)} \right|}}\int_{{z_0}}^z {{\rm{d}}x\,{\eta _1}\left( {x,\theta ,\chi ,v} \right)} . $(11)

A large variety of numerical schemes are available to approximate the integral in Eq. (10). In particular, a linear approximation of Si 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 = 4NsΝθΝχΝν, we considered the collocation vectors of the Stokes and emission vectors, I and ϵ ∈ ℝN, respectively, ordered with a lexicographic criterion, namely I=[I1(z1,θ1,χ1,v1),,I4(z1,θ1,χ1,v1),I1(z1,θ1,χ1,v2),,,I3(zNs,θNθ,χNv),I4(zNs,θNθ,χNχ,vNv)]T=[ε1(z1,θ1,χ1,v1),,ε4(z1,θ1,χ1,v1),ε1(z1,θ1,χ1,v2),,,ε3(zNs,θNθ,χNχ,vNv),ε4(zNs,θNθ,χNχ,vNv)]T.$ \matrix{ \hfill {{\bf{I}} = [{I_1}\left( {{z_1},{\theta _1},{\chi _1},{v_1}} \right), \ldots ,{I_4}\left( {{z_1},{\theta _1},{\chi _1},{v_1}} \right),{I_1}\left( {{z_1},{\theta _1},{\chi _1},{v_2}} \right), \ldots ,} \cr \hfill { \ldots ,{I_3}\left( {{z_{{N_{\rm{s}}}}},{\theta _{{N_\theta }}},{\chi _{{N_v}}}} \right),{I_4}\left( {{z_{{N_{\rm{s}}}}},{\theta _{{N_\theta }}},{\chi _{{N_\chi }}},{v_{{N_v}}}} \right){]^T}} \cr \hfill { \in = [{\varepsilon _1}\left( {{z_1},{\theta _1},{\chi _1},{v_1}} \right), \ldots ,{\varepsilon _4}\left( {{z_1},{\theta _1},{\chi _1},{v_1}} \right),{\varepsilon _1}\left( {{z_1},{\theta _1},{\chi _1},{v_2}} \right), \ldots ,} \cr \hfill { \ldots ,{\varepsilon _3}\left( {{z_{{N_{\rm{s}}}}},{\theta _{{N_\theta }}},{\chi _N}_{_\chi },{v_{{N_v}}}} \right),{\varepsilon _4}\left( {{z_{{N_{\rm{s}}}}},{\theta _{{N_\theta }}},{\chi _{{N_\chi }}},{v_{{N_v}}}} \right){]^T}.} \cr } $

For example, the first 4ΝθNχΝν elements of I and ϵ correspond to the first spatial point z1. We refer to algorithms 1–2 for the precise indices correspondence relation. For notational simplicity, an abuse of notation is employed: the symbol Ii 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 I=Λ+t$ {\bf{I}} = \Lambda + {\bf{t}} $(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 Iin instead of the zero boundary conditions.

Similarly, the discrete computation of the emissivity of Eq. (2) can be written as ϵ=ϵsc+ϵth=ΣI+ϵth,$ = {^{{\rm{sc}}}} + {^{{\rm{th}}}} = \Sigma {\bf{I}} + {^{{\rm{th}}}}, $(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 Σ=ΣII+ΣIII,$ \Sigma = {\Sigma ^{{\rm{II}}}} + {\Sigma ^{{\rm{III}}}}, $(14)

where ΣII and ΣIII encode the contributes of RII and RIII , 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 (IdΛΣ)I=Λϵth+t.$ \left( {Id - \Lambda \Sigma } \right){\bf{I}} = \Lambda {^{{\rm{th}}}} + t. $(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,vn) 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θ Nv Nχ) operations are required to apply Σ. Since the redistribution matrix RIII 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 106 for 1D models and reaches 1010 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.

thumbnail 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.

thumbnail Fig. 3

Run times for a single application of Σ, ΣIII, and Λ using algorithms 1–2 varying the number of processors NpNs 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−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 P=IdΛΣIII,$ P = Id - \Lambda {\Sigma ^{{\rm{III}}}}, $(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−1 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−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 light-weight application of Λ, which is sequential by nature. Parallelization on NpNs 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 ⌈Ns/Np⌉ or to ⌈Ns/Np⌉ 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 ⌈Ns/Np⌉ blocks), and the run time with Np processors, denoted by T(Np), can be bounded by T(Np) ≤ T(Ns) • ⌈Ns/Np⌉.

thumbnail 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 Ns/Np ∈ ℕ), 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χNv.

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)2. 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−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.

Table 1

Number of iterations until convergence of GMRES and PGM-RES for the solution of Eq. (15), varying Nθ = Nχ and Nv with Ns = 70.

thumbnail 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 Nv = 99 nodes. The problem was solved in the 1D semi-empirical model C of Fontenla et al. (1993), with [zmin, zmax] = [−100, 2219] km, discretized with Ns = 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 Nv, 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 Np 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 Ns. This is a key finding in view of 3D atmospheric models, where Ns is expected to grow considerably.

thumbnail 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.

Table 2

Run time (in minutes) to convergence for the solution of Eq. (15) using Ns = 128, Νθ = 12, Nχ = 8, and Νν = 99, i.e., Ν ≃ 4.9 × 106 total unknowns.

Table 3

Run time (in minutes) to convergence for the solution of Eq. (15) varying Ns = Np and with Νθ = 12, Nχ = 8, and Νν = 99; corresponding iterations are reported in square brackets.

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: υb(z)=0.018s1.z,$ {\upsilon _{\rm{b}}}\left( z \right) = 0.018\,{{\rm{s}}^{ - 1}}.z, $

with the atmospheric height z given in km and vb 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 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).

thumbnail 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 (μNθ/2+2,χ1)=(0.17,0)$\left( {{\mu _{{{{N_\theta }} \mathord{\left/{\vphantom {{{N_\theta }} {2 + 2}}} \right. \kern-\nulldelimiterspace} {2 + 2}}}},{\chi _1}} \right) = \left( {0.17,0} \right)$. 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.

thumbnail 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 (μNθ/2+2,χ1)=(0.17,0)$\left( {{\mu _{{{{N_\theta }} \mathord{\left/{\vphantom {{{N_\theta }} {2 + 3}}} \right.\kern-\nulldelimiterspace} {2 + 3}}}},{\chi _1}} \right) = \left( {0.38,0} \right)$ (leftpanels) and (μNθ/2+3,χ1)=(0.38,0)$\left( {{\mu _{{{{N_\theta }} \mathord{\left/{\vphantom {{{N_\theta }} {2 + 3}}} \right.\kern-\nulldelimiterspace} {2 + 3}}}},{\chi _1}} \right) = \left( {0.38,0} \right)$ (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) εi,th(r,Ω,v)=kL(r)ΓI(r)ΓR+ΓI(r)B(r)K=02Φ00K(r,Ω,v)^0,iK(r,Ω),$ \varepsilon _i^{\ell ,{\rm{th}}}\left( {r,\Omega ,v} \right) = {k_L}\left( r \right){{{\Gamma _I}\left( r \right)} \over {{\Gamma _R} + {\Gamma _I}\left( r \right)}}B\left( r \right)\sum\limits_{K = 0}^2 {\Phi _0^{0K}} \left( {r,\Omega ,v} \right)\hat \Im _{0,i}^K\left( {r,\Omega } \right), $

where i = 1 ,…, 4, B is the Planck function (if stimulated emission is neglected, the Wien limit must be considered), and Φ00K{\Phi _0^{0K}} is the generalized profile defined in Landi Degl’Innocenti & Landolfi (2004, App. 13). The factor ΓI(r)ΓR+ΓI(r)$ {{{\Gamma _I}\left( r \right)} \over {{\Gamma _R} + {\Gamma _I}\left( r \right)}} $

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 εisc(r,Ω,v)=εi,sc(r,Ω,v)+εic,sc(r,Ω,v),εith(r,Ω,v)=εi,sc(r,Ω,v)+εic,th(r,Ω,v).$ \eqalign{ &amp; \varepsilon _i^{{\rm{sc}}}\left( {r,\Omega ,v} \right) = \varepsilon _i^{\ell {\rm{,sc}}}\left( {r,\Omega ,v} \right) + \varepsilon _i^{c{\rm{,sc}}}\left( {r,\Omega ,v} \right), \cr &amp; \varepsilon _i^{{\rm{th}}}\left( {r,\Omega ,v} \right) = \varepsilon _i^{\ell {\rm{,sc}}}\left( {r,\Omega ,v} \right) + \varepsilon _i^{c{\rm{,th}}}\left( {r,\Omega ,v} \right). \cr} $

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) εisc(r,Ω,v)=σ(r,v)KQ(1)QJQ,iK(Ω)JQK(r,v),$ \varepsilon _i^{{\rm{sc}}}\left( {r,\Omega ,v} \right) = \sigma \left( {r,v} \right)\sum\limits_{KQ} {{{\left( { - 1} \right)}^Q}} _{Q,i}^K\left( \Omega \right)J_{ - Q}^K\left( {r,v} \right), $

where σ is the continuum opacity for scattering, and JQ,iK${\cal J}_{Q,i}^K$ is the geometrical tensor introduced in Landi Degl’Innocenti & Landolfi (2004, Sect. 5.11), and defined in the vertical reference system. Finally, the quantity JQK$J_Q^K$ is the radiation field tensor, defined as JQK(r,v)=dΩ4πj=14JQ,jK(Ω)Ij(r,Ω,v).$ J_Q^K\left( {r,v} \right) = \oint {{{{\rm{d}}\Omega } \over {4\pi }}} \sum\limits_{j = 1}^4 {_{Q,j}^K\left( \Omega \right){I_j}\left( {r,\Omega ,v} \right)} . $

The continuum term can thus be written in an integral form analogous to Eq. (3) in terms of the following redistribution matrix Rijc(r,Ω,Ω,v,v)=δ(vv)σ(r,v)KQ(1)QJQ,iK(Ω)JQ,jK(Ω).$ R_{ij}^c\left( {r,\Omega ',\Omega ,v',v} \right) = \delta \left( {v - v'} \right)\sigma \left( {{\bf{r}},v} \right)\sum\limits_{KQ} {{{\left( { - 1} \right)}^Q}} _{Q,i}^K\left( \Omega \right)_{ - Q,j}^K\left( {\Omega '} \right). $

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 εic,th(r,Ω,v)=εIc,th(r,v)δi1.$ \varepsilon _i^{c,{\rm{th}}}\left( {{\bf{r}},\Omega ,v} \right) = \varepsilon _I^{c,{\rm{th}}}\left( {{\bf{r}},v} \right){\delta _{i1}}. $(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 ηI(r,Ω,v)=ηI(r,Ω,v)+ηIc(r,v).$ {\eta _I}\left( {{\bf{r}},\Omega ,v} \right) = \eta _I^\ell \left( {{\bf{r}},\Omega ,v} \right) + \eta _I^c\left( {{\bf{r}},v} \right). $(B.2)

Routines for the calculation of σ, εIc,th$\varepsilon _I^{c,{\rm{th}}}$, and ηIc$\eta _I^c$ are publicly available (e.g., Uitenbroek 2001).

Appendix C Conversion to Optical Depth

To convert the spatial grid { zl }l=1Ns$\left\{ {{z_l}} \right\}_{l = 1}^{{N_s}}$ to optical depth, the map g:[ zmin,zmax ][ τmin,τmax ]+$ g:\left[ {{z_{\min }},{z_{\max }}} \right] \to \left[ {{\tau _{\min }},{\tau _{\max }}} \right] \subset {_ + } $

can be considered, which in a 1D geometry is defined as the solution of the IVP g(z)=1| cos(θ) |η1(z,θ,χ,v),$ g'\left( z \right) = {1 \over {\left| {\cos \left( \theta \right)} \right|}}{\eta _1}\left( {z,\theta ,\chi ,v} \right), $(C.1)

with g(zmin)=τmin   and  g(zmax)=τmax forθ[ 0,π/2 ],g(zmin)=τmin   and  g(zmax)=τmax forθ[ π/2,π ].$ \eqalign{ &amp; g\left( {{z_{\min }}} \right) = {\tau _{\min }}\,\,\,{\rm{and}}\,\,g\left( {{z_{\max }}} \right) = {\tau _{\max }}\,{\rm{for}}\theta \in \left[ {0,\pi /2} \right], \cr &amp; g\left( {{z_{\min }}} \right) = {\tau _{\min }}\,\,\,{\rm{and}}\,\,g\left( {{z_{\max }}} \right) = {\tau _{\max }}\,{\rm{for}}\theta \in \left[ {\pi /2,\pi } \right]. \cr} $

From the IVP (C.1), we obtain g(z)=τmin+1| cos(θ) |zminzdxη1(x,θ,χ,v) for θ[ 0,π/2 ],$ g\left( z \right) = {\tau _{\min }} + {1 \over {\left| {\cos \left( \theta \right)} \right|}}\int_{{z_{\min }}}^z {{\rm{d}}x{\eta _1}\left( {x,\theta ,\chi ,v} \right)} \,{\rm{for}}\,\theta \in \left[ {0,\pi /2} \right], $(C.2) g(z)=τmin1| cos(θ) |zminzdxη1(x,θ,χ,v) for θ[ π/2,π ],$ g\left( z \right) = {\tau _{\min }} - {1 \over {\left| {\cos \left( \theta \right)} \right|}}\int_{{z_{\min }}}^z {{\rm{d}}x{\eta _1}\left( {x,\theta ,\chi ,v} \right)} \,{\rm{for}}\,\theta \in \left[ {\pi /2,\pi } \right], $(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

  1. Alsina Ballester, E., Belluzzi, L., & Trujillo Bueno, J. 2016, ApJ, 831, L15 [NASA ADS] [CrossRef] [Google Scholar]
  2. Alsina Ballester, E., Belluzzi, L., & Trujillo Bueno, J. 2017, ApJ, 836, 6 [NASA ADS] [CrossRef] [Google Scholar]
  3. Alsina Ballester, E., Belluzzi, L., & Trujillo Bueno, J. 2018, ApJ, 854, 150 [NASA ADS] [CrossRef] [Google Scholar]
  4. Anusha, L.S., & Nagendra, K.N. 2011, ApJ, 739, 40 [NASA ADS] [CrossRef] [Google Scholar]
  5. Anusha, L.S., & Nagendra, K.N. 2012, ApJ, 746, 84 [NASA ADS] [CrossRef] [Google Scholar]
  6. Anusha, L., Nagendra, K., Paletou, F., & Léger, L. 2009, ApJ, 704, 661 [CrossRef] [Google Scholar]
  7. Badri, M., Jolivet, P., Rousseau, B., & Favennec, Y. 2018, J. Comput. Phys., 360, 74 [NASA ADS] [CrossRef] [Google Scholar]
  8. Badri, M., Jolivet, P., Rousseau, B., & Favennec, Y. 2019, Comput. Math. Applic., 77, 1453 [CrossRef] [Google Scholar]
  9. 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]
  10. Balay, S., Abhyankar, S., Adams, M.F., et al. 2015a, PETSc Users Manual, Tech. Rep. ANL-95/11 - Revision 3.6, Argonne National Laboratory [Google Scholar]
  11. Balay, S., Abhyankar, S., Adams, M.F., et al. 2015b, PETSc Web page, http://www.mcs.anl.gov/petsc [Google Scholar]
  12. Belluzzi, L., & Trujillo Bueno, J. 2012, ApJ, 750, L11 [NASA ADS] [CrossRef] [Google Scholar]
  13. Belluzzi, L., & Trujillo Bueno, J. 2014, A&A, 564, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  14. Belluzzi, L., Trujillo Bueno, J., & Stepân, J. 2012, ApJ, 755, L2 [NASA ADS] [CrossRef] [Google Scholar]
  15. Benedusi, P., Janett, G., Belluzzi, L., & Krause, R. 2021, A&A, 655, A88 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. Bommier, V. 1997a, A&A, 328, 706 [NASA ADS] [Google Scholar]
  17. Bommier, V. 1997b, A&A, 328, 726 [NASA ADS] [Google Scholar]
  18. Bommier, V. 2017, A&A, 607, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  19. Carlin, E.S., Manso Sainz, R., Asensio Ramos, A., & Trujillo Bueno, J. 2012, ApJ, 751, 5 [NASA ADS] [CrossRef] [Google Scholar]
  20. Carlin, E.S., Asensio Ramos, A., & Trujillo Bueno, J. 2013, ApJ, 764, 40 [NASA ADS] [CrossRef] [Google Scholar]
  21. Casini, R., & Landi Degl'Innocenti, E. 2008, Astrophysical Plasmas, 44, in eds. T. Fujimoto, & A. Iwamae, Atsushi, 247 [NASA ADS] [Google Scholar]
  22. 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]
  23. Casini, R., del Pino Alemân, T., & Manso Sainz, R. 2017a, ApJ, 835, 114 [NASA ADS] [CrossRef] [Google Scholar]
  24. Casini, R., del Pino Alemân, T., & Manso Sainz, R. 2017b, ApJ, 848, 99 [NASA ADS] [CrossRef] [Google Scholar]
  25. del Pino Alemân, T., & Trujillo Bueno, J. 2021, ApJ, 909, 180 [CrossRef] [Google Scholar]
  26. del Pino Alemân, T., Manso Sainz, R., & Trujillo Bueno, J. 2014, ApJ, 784, 46 [CrossRef] [Google Scholar]
  27. del Pino Alemân, T., Casini, R., & Manso Sainz, R. 2016, ApJ, 830, L24 [CrossRef] [Google Scholar]
  28. del Pino Alemân, T., Trujillo Bueno, J., Casini, R., & Manso Sainz, R. 2020, ApJ, 891, 91 [CrossRef] [Google Scholar]
  29. Fabiani Bendicho, P., Trujillo Bueno, J., & Auer, L. 1997, A&A, 324, 161 [NASA ADS] [Google Scholar]
  30. Faurobert-Scholl, M. 1992, A&A, 258, 521 [NASA ADS] [Google Scholar]
  31. Fontenla, J.M., Avrett, E.H., & Loeser, R. 1993, ApJ, 406, 319 [NASA ADS] [CrossRef] [Google Scholar]
  32. 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 A-6995 A (Zurich: vdf ETH) [NASA ADS] [Google Scholar]
  33. 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 A-4630 A (Zurich: vdf ETH) [NASA ADS] [Google Scholar]
  34. Guderley, G., & Hsu, C.-C. 1972, Math. Comp., 26, 51 [CrossRef] [Google Scholar]
  35. Holzreuter, R., Fluri, D.M., & Stenflo, J.O. 2005, A&A, 434, 713 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  36. Hummer, D.G. 1962, MNRAS, 125, 21 [NASA ADS] [CrossRef] [Google Scholar]
  37. Ishikawa, R., Bueno, J.T., del Pino Alemân, T., et al. 2021, Sci. Adv., 7, eabe8406 [NASA ADS] [CrossRef] [Google Scholar]
  38. Janett, G., & Paganini, A. 2018, ApJ, 857, 91 [NASA ADS] [CrossRef] [Google Scholar]
  39. Janett, G., Carlin, E.S., Steiner, O., & Belluzzi, L. 2017a, ApJ, 840, 107 [NASA ADS] [CrossRef] [Google Scholar]
  40. Janett, G., Steiner, O., & Belluzzi, L. 2017b, ApJ, 845, 104 [NASA ADS] [CrossRef] [Google Scholar]
  41. Janett, G., Steiner, O., & Belluzzi, L. 2018, ApJ, 865, 16 [NASA ADS] [CrossRef] [Google Scholar]
  42. Janett, G., Alsina Ballester, E., Guerreiro, N., et al. 2021a, A&A, 655, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  43. Janett, G., Benedusi, P., Belluzzi, L., & Krause, R. 2021b, A&A, 655, A87 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  44. Kano, R., Trujillo Bueno, J., Winebarger, A., et al. 2017, ApJ, 839, L10 [NASA ADS] [CrossRef] [Google Scholar]
  45. Landi Degl'Innocenti, E., & Landolfi, M. 2004, Polarization in Spectral Lines (Dordrecht: Kluwer Academic Publishers) Astrophysics and Space Science Library, 307 [Google Scholar]
  46. Landi Degl'Innocenti, E., Landi Degl'Innocenti, M., & Landolfi, M. 1997, in Forum THÉMIS, Science with THÉMIS, eds. N. Mein, & S. Sahal-Bréchot (Paris: Obs. de Paris), 59 [Google Scholar]
  47. Leenaarts, J., Pereira, T., & Uitenbroek, H. 2012, A&A, 543, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  48. Mihalas, D. 1978, Stellar Atmospheres, 2nd edn. (San Francisco: W.H. Freeman and Company) [Google Scholar]
  49. 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]
  50. Nagendra, K.N., Frisch, H., & Faurobert, M. 2002, A&A, 395, 305 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  51. Nagendra, K.N., Sowmya, K., Sampoorna, M., Stenflo, J.O., & Anusha, L.S. 2020, ApJ, 898, 49 [NASA ADS] [CrossRef] [Google Scholar]
  52. Rees, D.E., & Saliba, G.J. 1982, A&A, 115, 1 [NASA ADS] [Google Scholar]
  53. Rees, D., Murphy, G., & Durrant, C. 1989, ApJ, 339, 1093 [NASA ADS] [CrossRef] [Google Scholar]
  54. Ren, K., Zhang, R., & Zhong, Y. 2019, J. Comput. Phys., 399, 108958 [NASA ADS] [CrossRef] [Google Scholar]
  55. Sampoorna, M. 2011, ApJ, 731, 114 [NASA ADS] [CrossRef] [Google Scholar]
  56. Sampoorna, M., & Nagendra, K.N. 2015, ApJ, 812, 28 [NASA ADS] [CrossRef] [Google Scholar]
  57. Sampoorna, M., Nagendra, K.N., & Stenflo, J.O. 2008, ApJ, 679, 889 [NASA ADS] [CrossRef] [Google Scholar]
  58. Sampoorna, M., Stenflo, J.O., Nagendra, K.N., et al. 2009, ApJ, 699, 1650 [NASA ADS] [CrossRef] [Google Scholar]
  59. Sampoorna, M., Nagendra, K.N., & Stenflo, J.O. 2017, ApJ, 844, 97 [NASA ADS] [CrossRef] [Google Scholar]
  60. Stenflo, J. 1994, Solar Magnetic Fields: Polarized Radiation Diagnostics, 189 (Springer) [Google Scholar]
  61. Štěpán, J., Trujillo Bueno, J., Leenaarts, J., & Carlsson, M. 2015, ApJ, 803, 65 [CrossRef] [Google Scholar]
  62. Supriya, H.D., Smitha, H.N., Nagendra, K.N., Ravindra, B., & Sampoorna, M. 2013, MNRAS, 429, 275 [NASA ADS] [CrossRef] [Google Scholar]
  63. Supriya, H.D., Smitha, H.N., Nagendra, K.N., et al. 2014, ApJ, 793, 42 [NASA ADS] [CrossRef] [Google Scholar]
  64. 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]
  65. Trujillo Bueno, J., & Manso Sainz, R. 1999, ApJ, 516, 436 [CrossRef] [Google Scholar]
  66. Trujillo Bueno, J., Landi Degl'Innocenti, E., & Belluzzi, L. 2017, Space Sci. Rev., 210, 183 [NASA ADS] [CrossRef] [Google Scholar]
  67. Trujillo Bueno, J., Stepân, J., Belluzzi, L., et al. 2018, ApJ, 866, L15 [NASA ADS] [CrossRef] [Google Scholar]
  68. Uitenbroek, H. 2001, ApJ, 557, 389 [NASA ADS] [CrossRef] [Google Scholar]
  69. Warsa, J., Thompson, K., & Morel, J. 2003, Trans. Am. Nuclear Soci., 449 [Google Scholar]
  70. Weisskopf, V., & Wigner, E. 1930, Z. Phys., 63, 54 [CrossRef] [Google Scholar]

1

It is interesting to observe that if we artificially set αQ = 0 while keeping the nominal value of βQK$\beta _Q^{K}$, 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

Table 1

Number of iterations until convergence of GMRES and PGM-RES for the solution of Eq. (15), varying Nθ = Nχ and Nv with Ns = 70.

Table 2

Run time (in minutes) to convergence for the solution of Eq. (15) using Ns = 128, Νθ = 12, Nχ = 8, and Νν = 99, i.e., Ν ≃ 4.9 × 106 total unknowns.

Table 3

Run time (in minutes) to convergence for the solution of Eq. (15) varying Ns = Np and with Νθ = 12, Nχ = 8, and Νν = 99; corresponding iterations are reported in square brackets.

All Figures

thumbnail 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
thumbnail 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
thumbnail Fig. 3

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

In the text
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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 (μNθ/2+2,χ1)=(0.17,0)$\left( {{\mu _{{{{N_\theta }} \mathord{\left/{\vphantom {{{N_\theta }} {2 + 2}}} \right. \kern-\nulldelimiterspace} {2 + 2}}}},{\chi _1}} \right) = \left( {0.17,0} \right)$. 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
thumbnail 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 (μNθ/2+2,χ1)=(0.17,0)$\left( {{\mu _{{{{N_\theta }} \mathord{\left/{\vphantom {{{N_\theta }} {2 + 3}}} \right.\kern-\nulldelimiterspace} {2 + 3}}}},{\chi _1}} \right) = \left( {0.38,0} \right)$ (leftpanels) and (μNθ/2+3,χ1)=(0.38,0)$\left( {{\mu _{{{{N_\theta }} \mathord{\left/{\vphantom {{{N_\theta }} {2 + 3}}} \right.\kern-\nulldelimiterspace} {2 + 3}}}},{\chi _1}} \right) = \left( {0.38,0} \right)$ (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

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.