© The Authors 2023

1. Introduction

Significant scattering polarization signals are observed in several strong resonance lines of the solar spectrum, such as H I Ly-α (Kano et al. 2017), Mg II k (Rachmeler et al. 2022), Ca II K, Ca I 4227 Å, and Na I D₂ (e.g., Stenflo et al. 1980; Stenflo & Keller 1997; Gandorfer 2000, 2002). These signals, which are characterized by broad profiles with large amplitudes in the line wings, encode a variety of information on the thermodynamic and magnetic properties of the upper layers of the solar atmosphere (e.g., Trujillo Bueno 2014; Trujillo Bueno et al. 2017). A correct modeling of these profiles requires solving the radiative transfer (RT) problem for polarized radiation in nonlocal thermodynamic equilibrium (non-LTE) conditions, taking partial frequency redistribution (PRD) effects (i.e., correlations between the frequencies of incoming and outgoing photons in scattering processes) into account (e.g., Faurobert-Scholl 1992; Stenflo 1994; Holzreuter et al. 2005; Belluzzi & Trujillo Bueno 2012).

A powerful formalism to describe PRD phenomena is that of the redistribution function (e.g., Hummer 1962; Mihalas et al. 1978), which is generalized to the redistribution matrix in the polarized case (e.g., Domke & Hubeny 1988; Stenflo 1994; Bommier 1997a,b). For resonance lines, the redistribution matrix is given by the sum of two terms, commonly labeled as R^II and R^III according to the notation introduced by Hummer (1962). The R^II matrix describes scattering processes that are coherent in frequency in the atomic reference frame, while R^III describes scattering processes that are totally incoherent in the same frame (e.g., Bommier 1997a,b)¹. The linear combination of R^II and R^III allows for frequency redistribution effects due to elastic collisions with neutral perturbers to be taken into account (e.g., Bommier 1997b). In the observer’s frame, the Doppler shifts due to the thermal motions of the atoms are responsible for further frequency redistribution effects. The Doppler effect actually induces a complex coupling between the frequencies and propagation directions of the incident and scattered radiation, which makes the evaluation of both R^II and R^III, as well as the solution of the whole non-LTE RT problem, notoriously challenging from the computational standpoint. For this reason, approximate expressions of the redistribution matrices in the observer’s frame, in which such coupling is loosened, have been proposed and extensively used. In this context, a common choice is to use the angle-averaged (AA) expression of R^II, hereafter R^II−AA (e.g., Mihalas et al. 1978; Rees & Saliba 1982; Bommier 1997b; Leenaarts et al. 2012), and the expression of R^III obtained under the assumption that the scattering processes described by this matrix are also totally incoherent in the observer’s frame, hereafter R^III−CRD (e.g., Mihalas et al. 1978; Bommier 1997b; Alsina Ballester et al. 2017). The latter is also referred to as the assumption of complete frequency redistribution (CRD) in the observer’s frame.

The impact and the range of validity of the R^II−AA approximation in the modeling of scattering polarization has been discussed by several authors (e.g., Faurobert 1987, 1988; Nagendra et al. 2002, 2020; Nagendra & Sampoorna 2011; Sampoorna et al. 2011, 2017; Anusha & Nagendra 2012; Sampoorna & Nagendra 2015; del Pino Alemán et al. 2020; Janett et al. 2021a). These studies have shown that the use of R^II−AA can introduce significant (and hardly predictable) inaccuracies in the modeling of the line-core signals, while it seems to be suited for modeling the wing lobes and their magnetic sensitivity through the magneto-optical (MO) effects. By contrast, little effort has been directed to determine the suitability of the R^III−CRD assumption when modeling scattering polarization. Although this approximation has no true physical justification, it proved to be suitable for modeling the intensity profiles of spectral lines (e.g., Chap. 13 of Mihalas et al. 1978, and references therein). However, Bommier (1997b) pointed out that it may lead to appreciable errors when polarization phenomena are taken into account. In that work, the author considered a 90° scattering process of an unpolarized beam of radiation in the presence of a weak magnetic field, and compared the polarization of the scattered radiation calculated considering the exact expression of R^III and the R^III−CRD approximation, finding appreciable differences between the two cases². The exact expression of R^III has been used in the non-LTE RT calculations with isothermal one-dimensional (1D) atmospheric models by Sampoorna et al. (2011) and Supriya et al. (2012) for the nonmagnetic case, Nagendra et al. (2002) and Supriya et al. (2013) in the weak field Hanle regime, and, more recently, Sampoorna et al. (2017) in the more general Hanle-Zeeman regime. In this last paper, the authors analyzed the suitability of the R^III−CRD approximation for modeling scattering polarization, considering an ideal spectral line. They first concluded that the use of R^III−CRD introduces some inaccuracies, especially in the line wings, when considering optically thin self-emitting slabs in the presence of weak magnetic fields (i.e., fields for which the ratio #x0393;_B between the magnetic splitting of the Zeeman sublevels and the natural width of the upper level of the considered transition is on the order of unity). Moreover, these discrepancies are much smaller when stronger magnetic fields (#x0393;_B ≈ 100) are considered (see Fig. 1 and Sect. 4 in Sampoorna et al. 2017). However, when considering an atmospheric model with a greater optical depth, they found that the use of R^III−CRD does not seem to produce any noticeable effect on the emergent Stokes profiles already for weak fields (see Fig. 2 in Sampoorna et al. 2017). To the best of the authors’ knowledge, the aforementioned works are the only ones reporting PRD calculations of scattering polarization performed using the exact expression of R^III.

Taking advantage of a new solution strategy for the non-LTE RT problem for polarized radiation, tailored for including PRD effects (see Benedusi et al. 2022), we provide a quantitative analysis of the suitability of the R^III−CRD approximation, considering more general and realistic settings than in previous studies. In particular, we show the results of non-LTE RT calculations of the scattering polarization signals of two different spectral lines (i.e., Ca I 4227 Å and Sr I 4607 Å) in 1D models of the solar atmosphere, both semi-empirical and extracted from 3D magneto-hydrodynamic (MHD) simulations, and accounting for the impact of realistic magnetic and bulk velocity fields.

The paper is structured as follows. Section 2 presents the basic equations for the RT problem for polarized radiation, as well as the redistribution matrix formalism for describing PRD effects. Section 3 briefly presents the discretization and algebraic formulation of the problem, together with the considered numerical solution strategy. Section 4 exposes analytical and numerical comparisons between R^III and R^III−CRD. Section 5 provides some general considerations on the role of R^III in the RT modeling of scattering polarization, and on the expected impact of R^III−CRD. Sections 6 and 7 report comparisons of non-LTE RT calculations of scattering polarization signals in realistic 1D settings, performed with R^III and R^III−CRD. Finally, Sect. 8 provides remarks and conclusions.

2. Transfer problem for polarized radiation

The intensity and polarization of a radiation beam are fully described by the four Stokes parameters I, Q, U, and V. The Stokes parameter I is the intensity, Q and U quantify the linear polarization, while V quantifies the circular polarization (e.g., Landi Degl’Innocenti & Landolfi 2004, hereafter LL04). The Stokes parameters and other physical quantities involved in the RT problem are generally functions of the frequency and propagation direction of the radiation beam under consideration, as well as of the spatial point and time. Hereafter, we assume stationary conditions so that all quantities are time-independent.

A beam of radiation propagating in a medium (e.g., the plasma of a stellar atmosphere) is modified as a result of the interaction with the particles therein and the possible presence of magnetic, electric, and bulk velocity fields. This modification is fully described by the RT equation for polarized radiation, consisting of coupled systems of first-order, inhomogeneous, ordinary differential equations. Defining the Stokes parameters as the four components of the Stokes vector I = (I, Q, U, V)^T = (I₁, I₂, I₃, I₄)^T ∈ ℝ⁴, the RT equation for a beam of radiation of frequency ν, propagating along the direction specified by the unit vector Ω = (θ, χ)∈[0, π]×[0, 2π) can be written as³

$$\begin{array}{c}\hfill {\mathbf{\nabla }}_{\mathbf{\Omega }}\mathit{I}(\mathit{r},\mathbf{\Omega },\nu )=-K(\mathit{r},\mathbf{\Omega },\nu )\mathit{I}(\mathit{r},\mathbf{\Omega },\nu )+\mathit{\epsilon }(\mathit{r},\mathbf{\Omega },\nu ),\end{array}$$(1)

where ∇_Ω denotes the directional derivative along Ω, and r ∈ D is the spatial point in the domain D ⊂ ℝ³. The propagation matrix K ∈ ℝ^4 × 4 is given by

$$\begin{array}{c}\hfill K=\left(\begin{array}{cccc}{\eta }_1& {\eta }_2& {\eta }_3& {\eta }_4\\ {\eta }_2& {\eta }_1& {\rho }_4& -{\rho }_3\\ {\eta }_3& -{\rho }_4& {\eta }_1& {\rho }_2\\ {\eta }_4& {\rho }_3& -{\rho }_2& {\eta }_1\\ \end{array}\right).\end{array}$$(2)

The elements of K describe absorption (η₁), dichroism (η₂, η₃, and η₄), and anomalous dispersion (ρ₂, ρ₃, and ρ₄) phenomena. The emission vector ε = (ε₁, ε₂, ε₃, ε₄)^T ∈ ℝ⁴ describes the radiation emitted by the plasma in the four Stokes parameters.

In this work, we consider the contribution to K and ε brought by both line and continuum processes. These contributions, labeled with the superscripts ℓ and c, respectively, simply add to each other:

$$\begin{array}{cc}\hfill K(\mathit{r},\mathbf{\Omega },\nu )& =K^{\ell }(\mathit{r},\mathbf{\Omega },\nu )+K^c(\mathit{r},\mathbf{\Omega },\nu ),\hfill \\ \hfill \mathit{\epsilon }(\mathit{r},\mathbf{\Omega },\nu )& ={\mathit{\epsilon }}^{\ell }(\mathit{r},\mathbf{\Omega },\nu )+{\mathit{\epsilon }}^c(\mathit{r},\mathbf{\Omega },\nu ).\hfill \end{array}$$

In the frequency interval of a given spectral line, the line contribution to the elements of K and ε depends on the state of the atom (or molecule) giving rise to that line. In general, this state has to be determined by solving a set of rate equations (statistical equilibrium equations), which describe the interaction of the atom with the radiation field (radiative processes), other particles present in the plasma (collisional processes), and external magnetic and electric fields. When the statistical equilibrium equations have an analytic solution, the line contribution to the emission vector can be directly related to the radiation field that illuminates the atom (incident radiation) through the redistribution matrix formalism. More precisely, the line contribution to the emission vector can be written as the sum of two terms, namely,

$$\begin{array}{c}\hfill {\mathit{\epsilon }}^{\ell }(\mathit{r},\mathbf{\Omega },\nu )={\mathit{\epsilon }}^{\ell ,\mathrm{sc}}(\mathit{r},\mathbf{\Omega },\nu )+{\mathit{\epsilon }}^{\ell ,\mathrm{th}}(\mathit{r},\mathbf{\Omega },\nu ),\end{array}$$

where ε^ℓ, sc describes the contribution from atoms that are radiatively excited (scattering term), and ε^ℓ, th describes the contribution from atoms that are collisionally excited (thermal term). Following the convention that primed and unprimed quantities refer to the incident and scattered radiation, respectively, the line scattering term can be written through the scattering integral

$$\begin{array}{c}\hfill {\mathit{\epsilon }}^{\ell ,\mathrm{sc}}(\mathit{r},\mathbf{\Omega },\nu )=k_L(\mathit{r}){\displaystyle {\int }_{{\mathbb{R}}_{+}}}\mathrm{d}{\nu }^{\prime }{\displaystyle \oint \frac{\mathrm{d}{\mathbf{\Omega }}^{\prime }}{4\pi }}R(\mathit{r},\mathbf{\Omega },{\mathbf{\Omega }}^{\prime },\nu ,{\nu }^{\prime })\mathit{I}(\mathit{r},{\mathbf{\Omega }}^{\prime },{\nu }^{\prime }),\end{array}$$(3)

where the factor k_L is the frequency-integrated absorption coefficient, and R ∈ ℝ^4 × 4 the redistribution matrix. The redistribution matrix element R_ij relates the emissivity in the ith Stokes parameter in direction Ω and frequency ν to the jth Stokes parameter of the incident radiation with direction Ω′ and frequency ν′.

In this work, we consider an atomic system composed of two-levels (two-level atom) with an unpolarized and infinitely sharp lower level (see Appendix for more details), and we apply the corresponding redistribution matrix as derived in Bommier (1997b):

$$\begin{array}{c}\hfill R(\mathit{r},\mathbf{\Omega },{\mathbf{\Omega }}^{\prime },\nu ,{\nu }^{\prime })=R^{\mathrm{II}}(\mathit{r},\mathbf{\Omega },{\mathbf{\Omega }}^{\prime },\nu ,{\nu }^{\prime })+R^{\mathrm{III}}(\mathit{r},\mathbf{\Omega },{\mathbf{\Omega }}^{\prime },\nu ,{\nu }^{\prime }).\end{array}$$

The explicit expressions of R^II and R^III in both the atomic and observer’s reference frames are given in Appendix B. The line contribution to the elements of the propagation matrix K, and the line thermal emissivity ε^ℓ, th in the observer’s frame can be found in Appendix C. The continuum contribution to the emissivity can also be written as the sum of a thermal and a scattering term (see Eq. (D.1)), where the latter can be expressed as a scattering integral fully analogous to Eq. (3). A detailed discussion of the continuum terms is provided in Appendix D. For simplicity, hereafter we use the notation ε^th and ε^sc to refer to the thermal and scattering contributions to the emissivity, including both line and continuum.

The RT problem consists in finding a self-consistent solution for the RT equation (1) and the equation for the scattering contribution to the emissivity (3). This problem is in general nonlinear because of the factor k_L appearing in the line contribution to both the elements of the propagation matrix and the emission coefficients. This factor is proportional to the population of the lower level, which in turn depends nonlinearly on the incident radiation field through the statistical equilibrium equations.

We linearize the problem with respect to I, by fixing a priori the population of the lower level, and thus the factor k_L. In such scenario, whose suitability is discussed in Janett et al. (2021b) and Benedusi et al. (2021), the propagation matrix K and the thermal contribution to the emissivity ε^th is independent of I, while the scattering term ε^sc depends on it linearly through the scattering integral. The population of the lower level can be taken either from the atmospheric model (if provided) or from independent calculations. The latter can be carried out with available non-LTE RT codes that possibly neglect polarization (which is expected to have a minor impact on the population of ground or metastable levels), but allow considering multilevel atomic models. In this way, accurate estimates of the lower level population can be used, and reliable results can be obtained in spite of the simplicity of the considered two-level atomic model (e.g., Janett et al. 2021a; Alsina Ballester et al. 2021).

3. Numerical solution strategy

Following the works by Janett et al. (2021b) and Benedusi et al. (2021, 2022), we first present an algebraic formulation of the considered linearized RT problem for polarized radiation. Starting from this formulation, we then apply a parallel solution strategy, based on Krylov iterative methods with physics-based preconditioning. This strategy allows us to solve the problem in semi-empirical 1D models of the solar atmosphere, considering the exact expression of both R^II and R^III redistribution matrices. The same approach, coupled with a new domain decomposition technique, has been recently generalized to the 3D case (Benedusi et al. 2023).

The problem presented in Sect. 2 is discretized by introducing suitable grids for the continuous variables r, Ω, and ν. The angular grid ${\left\{{\mathbf{\Omega }}_i\right\}}_{i=1}^{N_{\mathrm{\Omega }}}$ is determined by the quadrature rule chosen to solve the angular integral in Eq. (3). We consider a right-handed Cartesian reference system with z-axis directed along the vertical. A given direction Ω is specified by the inclination θ and the azimuth χ, defined as shown in Fig. 1.

Fig. 1. Right-handed Cartesian reference system considered in the problem. The z-axis is directed along the local vertical. Any vector is specified through its polar angles θ (inclination) and χ (azimuth).

We discretize the inclination through two Gauss-Legendre grids with N_θ/2 points each, one for μ = cos(θ)∈(−1, 0), and one for μ ∈ (0, 1). For the azimuth, we use a grid of N_χ equidistant points. The angular quadrature is then the spherical Cartesian product (e.g., Davis & Rabinowitz 2007) of the Gauss-Legendre rule for the inclination and the trapezoidal rule for the azimuth. This approach is commonly used in RT applications and allows for the implementation of fast algorithms. The frequency grid ${\left\{{\nu }_j\right\}}_{j=1}^{N_{\nu }}$ is chosen to adequately sample the spectral line under investigation.

In a 1D (plane-parallel) setting, the spatial coordinate r = (x, y, z) can be replaced by the vertical coordinate z ∈ [z_min, z_max], and the RT equation (1) can be rewritten as

$$\begin{array}{c}\hfill cos(\theta )\frac{\mathrm{d}}{\mathrm{d}z}\mathit{I}(z,\theta ,\chi ,\nu )=-K(z,\theta ,\chi ,\nu )\mathit{I}(z,\theta ,\chi ,\nu )+\mathit{\epsilon }(z,\theta ,\chi ,\nu ).\end{array}$$(4)

The spatial grid ${\left\{z_k\right\}}_{k=1}^{N_z}$ is provided by the considered atmospheric model. We assume that the radiation entering the atmosphere from the lower boundary is isotropic, unpolarized, and equal to the Planck function, and that no radiation is entering from the upper boundary. Equation (4) is thus equipped with the following boundary conditions:

$$\begin{array}{cc}\hfill \mathit{I}(z_{\min },\theta ,\chi ,\nu )& ={[B_P(\nu ,T(z_{\min })),0,0,0]}^T\theta \in [0,\pi /2),\forall \chi ,\forall \nu ,\hfill \\ \hfill \mathit{I}(z_{\max },\theta ,\chi ,\nu )& =\mathbf{0}\theta \in (\pi /2,\pi ],\forall \chi ,\forall \nu ,\hfill \end{array}$$

where B_P is the Planck function and T the effective temperature. Given the propagation matrices and the emission vectors at all height points ${\left\{z_k\right\}}_{k=1}^{N_z}$, for a given direction (θ, χ) and frequency ν, the RT equation (1) can be numerically solved along that direction and at that frequency by applying a suitable integrator. This process is generally known as a formal solution. In this work, we use the L-stable DELO-linear method combined with a linear conversion to optical depth (e.g., Janett et al. 2017, 2018). An analysis of the stability properties of this method can be found in Janett & Paganini (2018).

With N = 4N_zN_θN_χN_ν being the total number of unknowns, we introduce the collocation vectors $\mathbf{I}\in {\mathbb{R}}^N$ and $\mathit{ϵ}\in {\mathbb{R}}^N$, which contain the numerical approximations of the Stokes vector and the emission vector, respectively, at all nodes. Given ϵ, the solution of all discretized RT equations (4) can be written as

$$\begin{array}{c}\hfill \mathbf{I}=\mathrm{\Lambda }\mathit{ϵ}+\mathbf{t},\end{array}$$(5)

where Λ : ℝ^N → ℝ^N is the transfer operator, which encodes the formal solver and the propagation matrix, and $\mathbf{t}\in {\mathbb{R}}^N$ is a vector encoding the boundary conditions.

Similarly, given I, the discrete computation of the emission vector can be written as

$$\begin{array}{c}\hfill \mathit{ϵ}={\mathit{ϵ}}^{\mathrm{sc}}+{\mathit{ϵ}}^{\mathrm{th}}=\mathrm{\Sigma }\mathbf{I}+{\mathit{ϵ}}^{\mathrm{th}},\end{array}$$(6)

where ϵ^sc and ϵ^th encode the scattering and thermal contributions (including both line and continuum processes), respectively, as described in Sect. 2. The scattering operator Σ : ℝ^N → ℝ^N encodes the numerical evaluation of the scattering integral (3) and thus depends on the chosen numerical quadratures. In general, the operator Σ is given by the sum of different components, namely,

$$\begin{array}{c}\hfill \mathrm{\Sigma }={\mathrm{\Sigma }}^{\mathrm{I}I}+{\mathrm{\Sigma }}^{\mathrm{I}II}+{\mathrm{\Sigma }}^{\mathrm{c}},\end{array}$$

where Σ^II and Σ^III encode the contributions from R^II and R^III, respectively, and Σ^c the scattering contribution from the continuum. The vector ϵ^th ∈ ℝ^N encodes the thermal emissivity.

Under the assumption that k_L is known a priori (see end of Sect. 2), the operators Λ and Σ are linear with respect to I. Combining Eq. (5) and Eq. (6), the whole discrete RT problem can be formulated as a linear system of size N with unknown I, namely

$$\begin{array}{c}\hfill (Id-\mathrm{\Lambda }\mathrm{\Sigma })\mathbf{I}=\mathrm{\Lambda }{\mathit{ϵ}}^{\mathrm{th}}+\mathbf{t},\end{array}$$(7)

where $Id\in {\mathbb{R}}^{N\times N}$ is the identity matrix. The right-hand-side vector Λϵ^th + t can be computed a priori by solving Eq. (4) with thermal contributions only (i.e., by performing a single formal solution with ϵ^sc = 0). The action of the matrices Λ, Σ, and Id − ΛΣ can be encoded in a matrix-free form (see Benedusi et al. 2022).

We solve the linear system (7) by applying a matrix-free, preconditioned GMRES method, with a relative tolerance of 10⁻⁷. Preconditioning is performed by describing scattering processes in the limit of CRD. This corresponds to substituting the operator Σ^II + Σ^III with a new operator Σ^CRD that is much cheaper to evaluate. The explicit expression of Σ^CRD for the considered atomic model can be found in Chap. 10 of LL04. The reader is referred to Benedusi et al. (2022) for more details on this solution strategy.

We conclude this section by briefly describing the numerical strategy used to handle bulk velocities. These velocities introduce Doppler shifts, which depend on the propagation direction of the considered radiation, in the expressions of the elements of K and ε (see Appendix B). We compute the emission vector in a reference frame in which the bulk velocity is zero (comoving reference frame). In this reference frame, there are no Doppler shifts, and this allows us to significantly reduce the number of evaluations of the redistribution matrix when performing the scattering integral (3). The drawback of this approach is that it requires transforming the incident radiation field from the observer’s frame to the comoving frame and then transforming the emission vector back to the observer’s frame. These steps are performed by means of high-order interpolations (e.g., cubic splines) on the frequency axis. The additional computational cost introduced by such interpolations is more than compensated by the reduced number of evaluations of the redistribution matrices.

4. Comparison of R^III and R^III−CRD

In this section, we first present an analytical comparison between the general angle-dependent expression of the R^III matrix (hereafter denoted as R^III) and its R^III−CRD approximation. We then review the reasons why the R^III−CRD approximation allows for a significant simplification of the problem from a computational point of view, and also provide a presentation of the challenges faced when considering R^III, focusing on algorithmic aspects.

4.1. Analytic considerations

In the formalism of the irreducible spherical tensors for polarimetry (see Chap. 5 of LL04), the R^III redistribution matrix in the observer’s reference frame (Eq. (B.29)) is given by the product between the scattering phase matrix ${\mathcal{P}}_Q^{K{K}^{\prime }}\in {\mathbb{C}}^{4\times 4}$ (Eq. (B.12)), which only depends on r, Ω, and Ω′, and the redistribution function ${\mathcal{R}}_Q^{\mathrm{III},K{K}^{\prime }}\in \mathbb{C}$ (Eq. (B.30)), which also depends on ν and ν′. This factorization also holds for the R^III−CRD matrix, leaving the ${\mathcal{P}}_Q^{K{K}^{\prime }}$ matrix unchanged, while replacing the redistribution function ${\mathcal{R}}_Q^{\mathrm{III},K{K}^{\prime }}$ with the approximation ${\mathcal{R}}_Q^{\mathrm{III}-\mathrm{CRD},K{K}^{\prime }}$ given by Eq. (B.35). A fundamental difference between ${\mathcal{R}}_Q^{\mathrm{III},K{K}^{\prime }}$ and ${\mathcal{R}}_Q^{\mathrm{III}-\mathrm{CRD},K{K}^{\prime }}$ is that the former depends on the scattering angle Θ = arccos(Ω ⋅ Ω′) (i.e., the angle between the directions Ω and Ω′), while the latter does not.

To analyze the dependence on Θ, we start considering the simpler expressions that ${\mathcal{R}}_Q^{\mathrm{III},K{K}^{\prime }}$ and ${\mathcal{R}}_Q^{\mathrm{III}-\mathrm{CRD},K{K}^{\prime }}$ assume in the absence of magnetic fields (B = 0). In this case, the magnetic shifts u_MuMℓ (see Eq. (B.20)) are zero, and the sums over the magnetic quantum numbers M appearing in Eqs. (B.30) and (B.36) can be performed analytically. If we additionally assume that there are no bulk velocities (v_b = 0) or, without loss of generality, we evaluate the redistribution matrix in the comoving frame (see Sect. 3), the Doppler shifts u_b (see Eq. (B.20)) also vanish, and the dependence on the propagation directions Ω and Ω′ is only through the scattering angle Θ. Expressing the functions in terms of the reduced frequencies u and u′ (see Eq. (B.19)), in the limit of B = 0 and v_b = 0, one finds:

$$\begin{array}{c}\hfill {\mathcal{R}}_Q^{\mathrm{III},K{K}^{\prime }}(\mathit{r},\mathbf{\Omega },{\mathbf{\Omega }}^{\prime },u,u\prime ){\underset{B=0,v_b=0}{\overset{\delta }{\to }}}_{K{K}^{\prime }}{\tilde{\mathcal{R}}}^{\mathrm{III},K}(\mathit{r},\mathrm{\Theta },u,u^{\prime }),\end{array}$$

with

$$\begin{array}{c}\hfill {\tilde{\mathcal{R}}}^{\mathrm{III},K}(\mathit{r},\mathrm{\Theta },u,u^{\prime })=({\beta }^K(\mathit{r})-\alpha (\mathit{r}))W^K(J_{\ell },J_u)\tilde{\mathcal{I}}(\mathit{r},\mathrm{\Theta },u,u^{\prime }).\end{array}$$(8)

The quantities β^K and α are given by Eqs. (B.3) and (B.4), respectively, in the limit of no magnetic fields, while W^K is defined as (see Eq. (10.17) of LL04)

$$\begin{array}{c}\hfill W^K(J_{\ell },J_u)=3(2J_u+1){\left\{\begin{array}{ccc}1& 1& K\\ J_u& J_u& J_{\ell }\end{array}\right\}}^2,\end{array}$$

where J_ℓ and J_u are the total angular momenta of the lower and upper level, respectively, and the operator in curly brackets is the Wigner’s 6-j symbol. The quantity $\stackrel{\sim }{\mathcal{I}}$ takes different analytical forms depending on the value of Θ:

• if Θ ∈ (0, π) (i.e., if Ω′≠Ω, −Ω)

  $$\begin{array}{cc}\hfill \tilde{\mathcal{I}}(\mathit{r},\mathrm{\Theta },u,u^{\prime })& =\frac{1}{{\pi }^{2}sin(\mathrm{\Theta })}{\displaystyle \int \mathrm{d}}y\exp (-y^2)\hfill \\ \hfill & \times H(\frac{a(\mathit{r})}{sin(\mathrm{\Theta })},\frac{u+ycos(\mathrm{\Theta })}{sin(\mathrm{\Theta })})\varphi (a(\mathit{r}),u^{\prime }+y),\hfill \end{array}$$(9)

• if Θ = 0 (i.e., forward scattering Ω′=Ω)

  $$\begin{array}{cc}\hfill \tilde{\mathcal{I}}(\mathit{r},\mathrm{\Theta }=0,u,u^{\prime })& =\frac{1}{{\pi }^{5/2}}{\displaystyle \int \mathrm{d}}y\exp (-y^2)\hfill \\ \hfill & \times \varphi (a(\mathit{r}),u+y)\varphi (a(\mathit{r}),u^{\prime }+y),\hfill \end{array}$$(10)

  and

• if Θ = π (i.e., backward scattering Ω′= − Ω)

  $$\begin{array}{cc}\hfill \tilde{\mathcal{I}}(\mathit{r},\mathrm{\Theta }=\pi ,u,u^{\prime })& =\frac{1}{{\pi }^{5/2}}{\displaystyle \int \mathrm{d}}y\exp (-y^2)\hfill \\ \hfill & \times \varphi (a(\mathit{r}),u-y)\varphi (a(\mathit{r}),u^{\prime }+y).\hfill \end{array}$$(11)

In the previous equations, H is the Voigt profile (i.e., the real part of the Faddeeva function, see Chap. 5 of LL04) and ϕ the Lorentzian profile (i.e., the real part of the function φ defined in Eq. (B.28)). Similarly, it can be shown that

$$\begin{array}{c}\hfill {\mathcal{R}}_Q^{\mathrm{III}-\mathrm{CRD},K{K}^{\prime }}(\mathit{r},\mathbf{\Omega },{\mathbf{\Omega }}^{\prime },u,u^{\prime }){\underset{B=0,v_b=0}{\overset{\delta }{\to }}}_{K{K}^{\prime }}{\tilde{\mathcal{R}}}^{\mathrm{III}-\mathrm{CRD},K}(\mathit{r},u,u^{\prime }),\end{array}$$

with

$$\begin{array}{cc}\hfill {\tilde{\mathcal{R}}}^{\mathrm{III}-\mathrm{CRD},K}(\mathit{r},u,u^{\prime })& =({\beta }^K(\mathit{r})-\alpha (\mathit{r}))W^K(J_{\ell },J_u)\hfill \\ \hfill & \times \frac{1}{\pi }H(a(\mathit{r}),u)H(a(\mathit{r}),u^{\prime }).\hfill \end{array}$$

Recalling that the Voigt profile is, by definition, a convolution between a Gaussian and a Lorentzian distribution (e.g., Chap. 5 of LL04), from Eq. (9), it can be easily seen that for Θ = π/2

$$\begin{array}{cc}\hfill \tilde{\mathcal{I}}(\mathit{r},\mathrm{\Theta }=\pi /2,u,u^{\prime })& =\frac{1}{{\pi }^2}H(a(\mathit{r}),u){\displaystyle \int \mathrm{d}}y\exp (-y^2)\varphi (a(\mathit{r}),u^{\prime }+y)\hfill \\ \hfill & =\frac{1}{\pi }H(a(\mathit{r}),u)H(a(\mathit{r}),u^{\prime }).\hfill \end{array}$$

This shows that the approximate ${\stackrel{\sim }{\mathcal{R}}}^{\mathrm{III}-\mathrm{CRD},K}$ function corresponds to the exact ${\stackrel{\sim }{\mathcal{R}}}^{\mathrm{III},K}$ for Θ = π/2, namely,

$$\begin{array}{c}\hfill {\tilde{\mathcal{R}}}^{\mathrm{III}-\mathrm{CRD},K}(\mathit{r},u,u^{\prime })={\tilde{\mathcal{R}}}^{\mathrm{III},K}(\mathit{r},\mathrm{\Theta }=\pi /2,u,u^{\prime }).\end{array}$$(12)

For scattering angles Θ ≠ π/2, the functions ${\stackrel{\sim }{\mathcal{R}}}^{\mathrm{III},K}$ and ${\stackrel{\sim }{\mathcal{R}}}^{\mathrm{III}-\mathrm{CRD},K}$ are generally different, and this difference increases as Θ approaches the two limiting cases of forward and backward scattering (i.e., Θ = 0 and Θ = π, respectively). This can be clearly seen in Fig. 2, where ${\stackrel{\sim }{\mathcal{R}}}^{\mathrm{III},0}$ is plotted as a function of u′ for different values of u and Θ. In particular, we compare the profiles for Θ = 0 and π, to that of Θ = π/2, which, as shown above, corresponds to ${\stackrel{\sim }{\mathcal{R}}}^{\mathrm{III}-\mathrm{CRD},0}$. For any value of u, the function ${\stackrel{\sim }{\mathcal{R}}}^{\mathrm{III}-\mathrm{CRD},0}$ shows a relatively broad profile (full width at half maximum of about 2), centered at u′ = 0. The amplitude of the peak is maximum at u = 0 (left panel) and quickly decreases by various orders of magnitude already at u ≈ 3 (middle panel). The behavior of the function ${\stackrel{\sim }{\mathcal{R}}}^{\mathrm{III},0}$ is much more complex. For u = 0 (left panel) and Θ = 0 or π, the profile is centered at u′ = 0 and it is much sharper and with a larger amplitude than that of ${\stackrel{\sim }{\mathcal{R}}}^{\mathrm{III}-\mathrm{CRD},0}$. For u ≈ 3 (middle panel) and Θ = 0 (resp. Θ = π), the function is characterized by a broad profile, similar in amplitude and width to that of Θ = π/2 (which is equivalent to ${\stackrel{\sim }{\mathcal{R}}}^{\mathrm{III}-\mathrm{CRD},0}$, see Eq. (12)), but with its maximum slightly shifted to positive (resp. negative) values of u′. Additionally, it shows a secondary sharp peak centered at u′=u (resp. u′= − u). For u ≈ 9 (right panel), in the case of Θ = 0 (resp. Θ = π), the secondary peak becomes negligible, while the main one is very similar to that of ${\stackrel{\sim }{\mathcal{R}}}^{\mathrm{III}-\mathrm{CRD},0}$. Noting that the dependence on K is limited to the factors β^K and W^K (see Eq. (8)), it is possible to state that ${\stackrel{\sim }{\mathcal{R}}}^{\mathrm{III}-\mathrm{CRD},K}$ and ${\stackrel{\sim }{\mathcal{R}}}^{\mathrm{III},K}$ are de facto equivalent in the wings, independently of the value of Θ. By contrast, they differ in the core and near wings, where the magnitude of the discrepancies strongly depends on Θ.

Fig. 2. Comparison of ${\stackrel{\sim }{\mathcal{R}}}^{\mathrm{III},0}$ as a function of u′ for three different scattering angles Θ (see legends), for u = 0 (left panel), u = 3.4 (middle panel), and u = 8.9 (right panel). The function is calculated considering a transition between levels with total angular momenta Jℓ = 0 and Ju = 1, and a damping parameter a = 0.01. The factor (β0 − α) is calculated setting α = 0 and using the value of β0 corresponding to the Ca I 4227 Å line at a height of 300 km in the FAL-C atmospheric model.

In the presence of a magnetic field, the redistribution function ${\mathcal{R}}_Q^{\mathrm{III},K{K}^{\prime }}$ (Eq. (B.30)) is given by a linear combination of the functions ${\mathcal{I}}_{(M_u{M}_{\ell }),(M_u^{\prime }{M}_{\ell }^{\prime })}$ (Eqs. (B.31)–(B.33)). These functions are fully analogous to $\stackrel{\sim }{\mathcal{I}}$ (Eqs. (9)–(11)), the only difference being that the second and third functions in the integrand (i.e., the profiles depending on the scattered and incident radiation, respectively) are shifted by u_MuMℓ and $u_{M_u^{\prime }{M}_{\ell }^{\prime }}$, respectively. It can be shown that also in the presence of magnetic fields (but still neglecting bulk velocities), the following relation holds:

$$\begin{array}{c}\hfill {\mathcal{R}}_Q^{\mathrm{III}-\mathrm{CRD},K{K}^{\prime }}(\mathit{r},u,u^{\prime })={\mathcal{R}}_Q^{\mathrm{III},K{K}^{\prime }}(\mathit{r},\mathrm{\Theta }=\pi /2,u,u^{\prime }).\end{array}$$

When Θ ≠ π/2, ${\mathcal{R}}_Q^{\mathrm{III},K{K}^{\prime }}$ (as function of u′) differs in general from ${\mathcal{R}}_Q^{\mathrm{III}-\mathrm{CRD},K{K}^{\prime }}$. As in the unmagnetized case, the difference is marginal for large values of u (i.e., u > 6) independently of the magnetic field strength, while it can be very significant in the core and near wings, especially when Θ is close to 0 or π. When approaching these limit cases, the curves for ${\mathcal{R}}_Q^{\mathrm{III},K{K}^{\prime }}$ can show very sharp peaks and the contribution of the various Zeeman components, split in frequency by the magnetic field, can eventually be resolved. As an illustrative example, Fig. 3 shows the component ${\mathcal{R}}_{-2}^{\mathrm{III},22}$ plotted as a function of u′, for u = 0.76, B = 30 G, and different values of Θ. We note that this component has a nonzero imaginary part since Q ≠ 0. As for the unmagnetized case, the curve for Θ = π/2, which corresponds to the CRD approximation, shows a broad profile centered at u′ = 0. As Θ departs from π/2, the corresponding profiles show increasingly large differences from the previous one. In particular, for Θ approaching 0 (resp. π), the position of the maximum moves from u′ = 0 toward u′=u (resp. u′= − u), while the width of the profile becomes sharper and the amplitude larger, in both the real and imaginary parts. As the profiles become sharper, the Zeeman components become visible, producing small lobes of opposite signs with respect to the central peak in the real part (left panels) and small substructures in the central peak in the imaginary part (right panels).

Fig. 3. Real (left column) and imaginary (right column) parts of ${\mathcal{R}}_Q^{\mathrm{III},K{K}^{\prime }}$ as a function of u′, for different scattering angles (see legend). We consider the component with K = K′ = 2 and Q = −2. The function is evaluated at u = 0.76, including a magnetic field of 30 G. The other parameters are the same as in Fig. 2.

In summary, the R^III−CRD approximation should be suitable in the far wings of the spectral lines (see right panel of Fig. 2), while it can introduce inaccuracies in the core and near wings, mainly due to scattering processes with Θ close to 0 or π, for which R^III−CRD and R^III differ significantly.

4.2. Computational considerations on R^III−CRD

The R^III−CRD approximation is based on the assumption that scattering processes are totally incoherent also in the observer’s frame, which in turn implies that in a scattering event, absorption and reemission can be treated as completely independent processes. Consistently with this picture, Eq. (B.35) shows that in R^III−CRD the joint probability of absorbing radiation with frequency ν′ and reemitting radiation with frequency ν is simply given by the product of two generalized profiles ${\mathrm{\Phi }}_Q^{K{K}^{\prime }}$ (Eq. (B.36)). Following the approach discussed at the end of Sect. 3, we evaluate the emission vector in the comoving frame, where no bulk velocities are present. In this case, the generalized profiles do not depend on the propagation directions Ω and Ω′, and from Eqs. (B.29), (B.35), and (B.12), one finds that the emission coefficients corresponding to R^III−CRD are given by

$$\begin{array}{cc}\hfill {\epsilon }_i^{\mathrm{III}-\mathrm{CRD}}(\mathit{r},& \mathbf{\Omega },\nu )=k_L(\mathit{r})\hfill \\ \hfill & \times {{\displaystyle \int {\mathbb{R}}_{+}}}_{}\mathrm{d}{\nu }^{\prime }{\displaystyle \oint \frac{\mathrm{d}{\mathbf{\Omega }}^{\prime }}{4\pi }\sum _{j=1}^4}{R}_{\mathit{ij}}^{\mathrm{III}-\mathrm{CRD}}(\mathit{r},\mathbf{\Omega },{\mathbf{\Omega }}^{\prime },\nu ,{\nu }^{\prime })I_j(\mathit{r},{\mathbf{\Omega }}^{\prime },{\nu }^{\prime })\hfill \\ \hfill =& \frac{k_L(\mathit{r})}{\mathrm{\Delta }{\nu }_D^2(\mathit{r})}{\displaystyle \sum _{K,K^{\prime }=0}^2}{\displaystyle \sum _{Q=-K_{min}}^{K_{min}}}{\displaystyle \sum _{K^{″}=|Q|}^{2J_u}}({\beta }_Q^{K^{″}}(\mathit{r})-{\alpha }_Q(\mathit{r}))\hfill \\ \hfill & \times {\displaystyle \sum _{Q^{\prime }=-K}^K}{\displaystyle \sum _{Q″=-K^{\prime }}^{K^{\prime }}}{(-1)}^{Q^{\prime }}{\mathcal{T}}_{Q″,i}^{K^{\prime }}(\mathbf{\Omega }){\overline{\mathcal{D}}}_{QQ″}^{K^{\prime }}(\mathit{r}){\mathcal{D}}_{Q{Q}^{\prime }}^K(\mathit{r})\hfill \\ \hfill & \times {\mathrm{\Phi }}_Q^{K^{″}{K}^{\prime }}(\mathit{r},\nu ){\displaystyle {\int }_{{\mathbb{R}}_{+}}}\mathrm{d}{\nu }^{\prime }{\mathrm{\Phi }}_Q^{K^{″}K}(\mathit{r},{\nu }^{\prime })J_{-Q^{\prime }}^K(\mathit{r},{\nu }^{\prime }),\hfill \end{array}$$(13)

where i = 1, …, 4, K_min = min(K, K′), and #x0394;ν_D, ${\mathcal{T}}_{Q,i}^K$, and ${\mathcal{D}}_{Q{Q}^{\prime }}^K$ are the Doppler width (in frequency units), the polarization tensor, and the rotation matrices, respectively (see Appendix B for more details). Finally, $J_Q^K$ is the radiation field tensor, defined as (see Eq. (5.157) of LL04)

$$\begin{array}{c}\hfill J_Q^K(\mathit{r},\nu )={\displaystyle \oint \frac{\mathrm{d}\mathbf{\Omega }}{4\pi }\sum _{j=1}^4}{\mathcal{T}}_{Q,j}^K(\mathit{r},\mathbf{\Omega })I_j(\mathit{r},\mathbf{\Omega },\nu ).\end{array}$$

Equation (13) shows that the dependencies on the frequency and propagation direction of the incoming and outgoing radiation are completely decoupled. This allows the implementation of simple and fast computational algorithms: once the values of $J_Q^K$ are obtained from the formal solution of the RT equation, it is possible to independently compute the values of ${\mathrm{\Phi }}_Q^{K^{\prime }K}$ for the frequencies of the incoming and outgoing radiation and combine them only during the final calculation of the emission coefficient. Using the big-O notation, the time complexity (e.g., Sipser 1996) of the algorithm scales as O$\left(N_{\nu }^d\right)$, where d ∈ (1, 2) grows monotonically with N_ν. For the typical size of the frequency grids needed to synthesize one (or a few) spectral lines (i.e., N_ν ≈ 100), d can be considered close to 1. This kind of time complexity is justified because the calculation of the generalized profile is slow (complex algorithms are required), while the subsequent combination of them is fast as it only implies products of complex numbers.

In numerical applications considering the R^III−CRD approximation, it is a common practice to perform the integration over the frequencies of the incoming radiation in Eq. (13) using the frequency grid from the problem discretization (see Sect. 3) without further refinements, and applying the trapezoidal rule. This methodology has the advantages of a decreased time-to-solution (TTS) and implementation simplicity. On the other hand, it should be pointed out that the accuracy of the results can be improved by using denser and more specific grids for the frequencies of the incoming radiation without any major loss of overall performance. This is because, even if a finer frequency grid is used in Eq. (13), in a PRD setting the TTS for evaluating the total emissivity remains mainly dominated by the contribution of the angle-dependent R^II redistribution matrix.

4.3. Computational considerations on R^III

The evaluation of ε^III considering the exact expression of R^III (see equations in Appendix B.2.2) is computationally challenging for the following main reasons: (a) it involves a 4-dimensional integration, leading to a very large number of evaluations of the integrands in Eqs. (B.31)–(B.33). This issue is otherwise known as the curse of dimensionality; (b) the integration variables (ν′, Ω′, and y) cannot be algebraically decoupled; (c) already for simple atomic transitions, the total number of combinations of magnetic quantum numbers coupled with the tensorial polarimetric indices is high, which adds another layer of complexity to the 4 dimensions of the overall problem; (d) the integrands include the Faddeeva function (Faddeeva & Terent’ev 1961) whose evaluation requires a large TTS (e.g., Oeftiger et al. 2016; e) the integrands show a very complex behavior with respect to the various integration variables and parameters, thus requiring the use of very fine, unstructured, and case-dependent grids.

Our overall approach for calculating ε^III is to first evaluate the integral over y, followed by that over the frequencies ν′, and finally, the one over the propagation direction Ω′. The analytic expressions of ℐ (Eqs. (B.31)–(B.33)) show that in the absence of bulk velocities (or if the redistribution matrix is evaluated in the comoving frame), the coupling between the propagation directions of the incoming and outgoing radiation occurs through the scattering angle Θ. The integrand in ℐ shows a complex behavior with respect to the integration variable y, especially when the scattering angle approaches 0 and π. In these cases, the integrand becomes close to a Lorentz distribution or its square. These functions are not easy to numerically integrate due to the presence of sharp peaks and extended wings. Indeed, the convergence rate for the numerical quadrature of the Lorentz distribution (and its square) is generally slow. Furthermore, the computation of a single ℐ (see Eq. (B.31)) can require up to thousands of evaluations of the Faddeeva function (accounting for more than 70% of the TTS). In order to perform the quadratures over y, we applied an adaptive quadrature method based on the Gauss-Kronrod approach (e.g., Kronrod 1965; Piessens et al. 2012). The advantage of this method is that it is capable of automatically inferring the behavior of the integrand by achieving very high accuracy with a relatively low number of function evaluations.

In general, the time complexity for the computation of ε^III is O$\left(N_{\mathrm{\Omega }}^2{N}_{\nu }^3\right)$. The cubic contribution given by the number of frequency grid points is due to the fourth dimension induced by y. It must be noted that the number of grid points needed to adequately perform the quadrature of the integral over y is generally larger than N_ν.

For the angular integral in Eq. (3), it is convenient to apply a quadrature rule characterized by a regular angular grid, because in this case, the effective number of different scattering angles is significantly lower than the total pairs of directions. The results reported in the next sections were obtained with the quadrature method described in Sect. 3, with 12 Gauss-Legendre inclinations and 8 equally spaced azimuths. In this case, the total number of scattering angles is limited to 200. In our algorithm, the quantities ℐ are precomputed for the whole set of different scattering angles corresponding to the chosen angular grid, thus avoiding repeating the calculation of ℐ, which is rather expensive, for different pairs of directions having the same scattering angle.

The quantity ℐ depends on two pairs of magnetic quantum numbers, (M_u, M_ℓ) and $(M_u^{\prime },M_{\mathcal{l}}^{\prime })$, where the labels u and ℓ indicate the upper and lower level, respectively. Equation (B.30) shows that the expression of R^III actually includes four quantities ℐ, which differ for the values of the magnetic quantum numbers and are coupled inside six nested loops over such quantum numbers. It can be easily realized that in such loops several tuples of magnetic sublevels are repeated, which suggests the opportunity to increase the efficiency of the algorithm by precomputing the quantities ℐ for all possible combinations of magnetic quantum numbers, thus avoiding to recalculate the same quantity ℐ various times. The precomputation of ℐ as described above leads to a drastic reduction of the total number of calculations needed to compute ε^III, which otherwise would be very significant due to the extra dimension in the integration. The precomputed values of ℐ are stored out-of-core (e.g., on disk) because they require a large footprint, and are accessed during the calculation of ε^III through a system of look-up tables. It must be observed that the quantity ℐ and, consequently, the grids used for performing the numerical integration, also depends on the spatial point through the damping parameter a, the magnetic field, and the Doppler width #x0394;ν_D. Our strategy thus calls for relatively large data storage capabilities that, however, remain manageable in the case of 1D applications.

5. Impact of R^III on spectral lines formation

The R^III redistribution matrix describes scattering processes during which the atom undergoes elastic collisions with neutral perturbers (mainly hydrogen and helium atoms in the solar atmosphere) that completely relax any correlation between the frequencies of the incident and scattered radiation, thus making the scattering totally incoherent. Informally speaking, due to such collisions, the atom does not keep any memory of the frequency of the incident photon. On the contrary, the R^II redistribution matrix describes scattering processes in which the atom is not subject to any elastic collisions so that the frequencies of the incident and scattered radiation remain fully correlated (coherent scattering). A quantitative estimate of the relative weight of R^III with respect to R^II is provided by the branching ratio for R^II (see Eq. (B.9)) in the absence of magnetic fields (a.k.a. coherence fraction):

$$\begin{array}{c}\hfill \tilde{\alpha }(\mathit{r})=\frac{A_{u\ell }+C_{u\ell }(\mathit{r})}{A_{u\ell }+C_{u\ell }(\mathit{r})+Q_{\mathrm{el}}(\mathit{r})},\end{array}$$

where A_uℓ is the Einstein coefficient for spontaneous emission, C_uℓ is the rate of inelastic de-exciting collisions, and Q_el is the rate of elastic collisions with neutral perturbers. A value of $\stackrel{\sim }{\alpha }$ close to unity (corresponding to a very low rate of elastic collisions compared to the rates for spontaneous emission and collisional de-excitation) means that R^II dominates over R^III, while a value of $\stackrel{\sim }{\alpha }$ close to zero (corresponding to a very high value of Q_el compared to A_uℓ and C_uℓ) means instead that R^III dominates with respect to R^II.

The impact of R^III is thus expected to be marginal in the core of strong spectral lines (i.e., lines showing broad intensity profiles with extended wings). Indeed, their line-core radiation generally originates from the upper layers of the solar atmosphere, where the number density of neutral perturbers (and thus the rate of elastic collisions) is relatively low. The relative weight of R^III can, however, become significant in the wings of such lines, as they usually form much lower in the atmosphere. On the other hand, it must be observed that the profiles entering the definition of R^III are all centered around the line-center frequency, and therefore the net contribution of this redistribution matrix to the emissivity in the line wings is generally marginal with respect to that of R^II. This can be seen as an analytical proof of the well-known fact that scattering processes are mainly coherent in the line wings (e.g., Mihalas et al. 1978). Consistently with this picture, Alsina Ballester et al. (2022) showed that in the wings of strong resonance lines, the contribution of R^III needs to be taken into account in order not to overestimate the weight of R^II, but its net contribution to scattering polarization is fully negligible. These considerations suggest that the exact analytical form of R^III should not be crucial for modeling the core and far wings of strong lines, and that the R^III−CRD approximation should therefore be suitable in these spectral ranges. The most critical regime is that of the near wings, where strong resonance lines may show very significant scattering polarization signals. There, R^III may bring a non-negligible contribution, and it is harder to estimate a priori the suitability of the R^III−CRD approximation.

The relative weight of R^III is also non-negligible in the case of spectral lines forming in the deeper layers of the solar atmosphere (photosphere), where the number density of colliders is significant. On the other hand, these lines are generally weak in the intensity spectrum, showing narrow absorption profiles with a Doppler core and no wings. Since Doppler redistribution is generally very efficient in the line-core, the limit of CRD (i.e., to assume that all scattering processes are totally incoherent) has always been considered a good approximation for modeling both the intensity (e.g., Mihalas et al. 1978) and scattering polarization profiles of these lines (e.g., LL04).

In order to quantitatively verify these considerations, and assess the suitability of the R^III−CRD approximation for modeling scattering polarization, we model the intensity and polarization of two different spectral lines, namely a strong spectral line with extended wings forming in the upper layers of the solar atmosphere, and a weaker line forming deeper in the atmosphere. Excellent examples for these two typologies of spectral lines are, respectively, the Ca I line at 4227 Å and the Sr I line at 4607 Å. Both lines result from a resonant transition between the ground level of the considered atomic species, which in both cases has total angular momentum J_ℓ = 0, and an excited level with J_u = 1. Both of them show conspicuous scattering polarization signals, which can be suitably modeled considering a two-level atom with an unpolarized and infinitely-sharp lower level.

Figure 4 shows the variation of the coherence fraction $\stackrel{\sim }{\alpha }$ with height in the 1D semi-empirical model C of Fontenla et al. (1993, hereafter FAL-C) for the two considered spectral lines. In the same figure, we also plot the height at which the optical depth τ, in the frequency intervals of the considered spectral lines, is unity. It can be shown (e.g., Mihalas et al. 1978) that this height provides an approximate estimate of the atmospheric region from which the emergent radiation originates (formation height). We recall that the optical depth at the frequency ν along direction Ω is defined as

$$\begin{array}{c}\hfill \tau (s,\mathbf{\Omega },\nu )={\displaystyle {\int }_0^s}{\eta }_1(x,\mathbf{\Omega },\nu )\mathrm{d}x,\end{array}$$(14)

Fig. 4. Height variation of $\stackrel{\sim }{\alpha }$ (blue line) for Sr I 4607 (top panel) and Ca I 4227 (bottom panel) in the FAL-C model. The red and dashed red isolines show the height where the optical depth as a function of wavelength is unity (i.e., τ = 1) for μ = 0.034 and μ = 0.966, respectively.

where s is the spatial coordinate along direction Ω (measured with respect to an arbitrary initial point), and η₁ is the absorption coefficient for the intensity (see also Janett et al. 2017). For calculating the formation height, the direction Ω must point inward in the atmosphere, and the initial point s = 0 is taken at the upper boundary. Since the absorption coefficient is largest at the line center and decreases moving to the wings, it can be immediately seen that the formation height is highest in the line core and decreases moving to the wings. From Eq. (14) it is also clear that, for a given frequency, the formation height is higher for a line of sight (LOS) close to the edge of the solar disk (limb) than for one at the disk center. Clearly, the formation height is higher for strong lines which have a large absorption coefficient (because the number density of atoms in the lower level of the transition is particularly large in the solar atmosphere) than for weak ones.

5.1. Ca I 4227 Å line

In the intensity spectrum, the Ca I 4227 Å line shows a very broad and deep absorption profile (e.g., Gandorfer 2002) with an equivalent width of 1426 mÅ (Moore et al. 1966), that is, one of the strongest spectral lines in the visible part of the solar spectrum. When observed in quiet regions close to the limb, this line shows a large scattering polarization signal with a sharp peak in the line core and broad lobes in the wings (e.g., Gandorfer 2002). This signal, with its peculiar triplet-peak structure, has been extensively observed and modeled in the past (e.g., Stenflo et al. 1980; Faurobert-Scholl 1992; Bianda et al. 2003; Sampoorna et al. 2009; Anusha et al. 2011; Supriya et al. 2014; Carlin & Bianda 2017; Alsina Ballester et al. 2018; Janett et al. 2021a; Jaume Bestard et al. 2021, and references therein). In particular, it was clearly established that its broad wing lobes are produced by coherent scattering processes with PRD effects. We recall that the line-core peak is sensitive to the presence of magnetic fields via the Hanle effect. The Hanle critical field (i.e., the magnetic field strength for which the sensitivity to the Hanle effect is maximum) is B_c ≈ 25 G. The wing lobes are sensitive to longitudinal magnetic fields of similar strength via the magneto-optical (MO) effects (Alsina Ballester et al. 2018).

The lower panel of Fig. 4 shows that in the FAL-C atmospheric model, the core of the Ca I 4227 Å line forms above 800 km (low chromosphere). At these heights, the number density of neutral perturbers is very low, the coherence fraction is thus very close to unity, and R^II dominates with respect to R^III. On the other hand, as we move from the core to the wings, the formation height quickly decreases to photospheric levels. At these heights, the coherence fraction is much lower than one and the weight of R^III is significant. On the other hand, for the reasons discussed above, its net impact in the line wings is expected to be relatively low. The most critical region is that of the near wings, where the Ca I 4227 Å line shows strong scattering polarization signals, and the net contribution from R^III can be non-negligible. The suitability of the R^III−CRD approximation in this spectral region can only be assessed numerically, and it will be analyzed in Sect. 6.

5.2. Sr I 4607 Å line

The Sr I 4607 Å line is a rather weak absorption line in the intensity spectrum, with an equivalent width of 36 mÅ only (see Moore et al. 1966). Nonetheless, this line shows a prominent scattering polarization signal at the solar limb, characterized by a sharp profile (e.g., Gandorfer 2002). This signal has been also extensively observed and modeled in the past, especially in order to investigate, via the Hanle effect, the small-scale, unresolved magnetic fields that permeate the quiet solar photosphere (e.g., Stenflo & Keller 1997; Trujillo Bueno et al. 2004; del Pino Alemán et al. 2018; del Pino Alemán & Trujillo Bueno 2021; Dhara et al. 2019; Zeuner et al. 2020, 2022). The Hanle critical field for this line is B_c ≈ 23 G. The limit of CRD has always been considered suited for modeling both the intensity and scattering polarization profiles of this line.

The upper panel of Fig. 4 shows that in the FAL-C atmospheric model, the core of this line forms in the photosphere, below 500 km. The curve for the coherence fraction shows that the weight of R^III is, as expected, significant at these heights. In the next section, we explore the suitability of the R^III−CRD approximation in the modeling of the scattering polarization signal of this line through full PRD RT calculations in the presence of both deterministic and nondeterministic (e.g., Stenflo 1982, 1994; Landi Degl’Innocenti & Landolfi 2004) magnetic fields. For the latter, we consider unimodal micro-structured isotropic (MSI) magnetic fields, namely magnetic fields with a given strength and an orientation that changes on scales below the mean free path of photons, uniformly distributed over all directions. The expressions of RT quantities in the presence of MSI magnetic fields are exposed in Appendix E (see also Sect. 4 of Alsina Ballester et al. 2017).

6. Numerical results: FAL-C atmospheric model

In this section and in the next one, we present the numerical results⁴ of non-LTE RT calculations of the scattering polarization profiles of the Ca I 4227 Å and Sr I 4607 Å lines, performed with the numerical solution strategy described in Sect. 3. All calculations are performed using the general (angle-dependent) expression of the R^II redistribution matrix, while considering both the exact form of the R^III matrix and its R^III−CRD approximation. The lower level population, which we keep fixed in the problem (see the end of Sect. 2), is precomputed with the RH code (Uitenbroek 2001), which solves the nonlinear unpolarized non-LTE RT problem. For the Ca I 4227 Å line, we run RH using a model for calcium composed of 25 levels, including 5 levels of Ca II and the ground level of Ca III. For the Sr I 4607 Å line, we considered a model composed of 34 levels, including the ground level of Sr II. These preliminary computations also provided the rates for elastic and inelastic collisions, as well as the continuum quantities.

For a quantitative comparison between the emergent Stokes profiles obtained using R^III and R^III−CRD, we plot the error defined by

$$\begin{array}{c}\hfill \text{Error}(\mathit{a},\mathit{b})=\frac{|\mathit{a}{(max\mathit{a})}^{-1}-\mathit{b}{(max\mathit{a})}^{-1}|}{1+\left|\mathit{a}{(max\mathit{a})}^{-1}\right|},\end{array}$$(15)

where a and b represent the values of a given Stokes parameter of the emergent radiation, for a given direction and all considered wavelengths, obtained considering R^III and R^III−CRD, respectively. The maximum is calculated with respect to wavelength, over the considered interval. The error definition in Eq. (15) does not correspond to the standard relative error and it was introduced to prevent amplifying the discrepancies where a and b are close to zero and consequently to the numerical noise. Where the signals a and b are relevant, this definition provides a value that is smaller, by a factor of two in the worst case (i.e., where a takes its maximum value), than the usual relative error. This definition is thus justified as it provides the correct order of magnitude of the error where the signal is relevant while damping it when the signal becomes negligible.

In this section, we show calculations performed in the FAL-C atmospheric model (N_z = 70 height nodes), in the presence of constant (i.e., height-independent) deterministic magnetic fields (i.e., magnetic fields having a well-defined strength and orientation at each spatial point), in the absence of bulk velocity fields. The LOS toward the observer is taken on the x − z plane of the considered reference system (see Fig. 1) and is specified by μ = cos θ, with θ the inclination with respect to the vertical. Typically, we present the emergent Stokes profiles for two specific directions: μ = 0.034, which represents radiation coming from the solar limb (nearly horizontal LOS), and μ = 0.996, which represents radiation coming from near the center of the solar disk (nearly vertical LOS).

6.1. Ca I 4227 Å line

We first consider the modeling of the Stokes profiles of the Ca I 4227 Å line, which is discretized with N_ν = 99 unevenly spaced nodes. The characteristics of this spectral line are adequately reproduced by our calculations. The triple-peak structure in Q/I is evident in the non-magnetized model shown in Fig. 5, while the magnetized case shown in Fig. 6 displays a clear depolarization of the core as a result of the Hanle effect and a mild depolarization of the lobes due to the magneto-optical effects.

Fig. 5. Emergent Stokes I (top panels) and Q/I (bottom panels) profiles of the Ca I 4227 Å line at μ = 0.034 (left panels) and μ = 0.966 (right panels) calculated in the FAL-C atmospheric model in the absence of magnetic fields. Calculations take into account PRD effects considering the exact expression of RIII (blue lines) and the RIII−CRD approximation (green marked lines). The reference direction for positive Q is taken parallel to the limb. The red dotted lines report the error between RIII and RIII−CRD calculations, given by Eq. (15).

Fig. 6. Emergent Stokes I (first row), Q/I (second row), U/I (third row), and V/I (fourth row) profiles for the Ca I 4227 Å line at μ = 0.034 (left column) and μ = 0.966 (right column) calculated in the FAL-C atmospheric model in the presence of a horizontal (θB = π/2,  χB = 0) magnetic field with B = 20 G. Calculations take into account PRD effects considering the exact expression of RIII (blue lines) and the RIII−CRD approximation (green marked lines). The reference direction for positive Q is taken parallel to the limb. The red dotted lines report the error between RIII and RIII−CRD calculations, given by Eq. (15).

Figure 5 shows the Stokes I and Q/I profiles for a non-magnetic scenario, where we observe no significant difference between R^III and R^III−CRD calculations, with an error that is always smaller than 7.5 × 10⁻³. In the absence of magnetic fields, the U/I and V/I profiles vanish and are consequently not shown.

Figure 7 shows the contributions from R^II, R^III, and R^III−CRD to |ε_I| and |ε_Q| for the non-magnetized case as a function of height. For the sake of completeness, the figure also shows the contributions from continuum and line thermal emissivities to Stokes I and from continuum to Stokes Q (see blue lines). In this figure, we consider μ = 0.17 and λ = 4227.1 Å, which corresponds to the wavelength at the maximum of the red Q/I lobe in Fig. 5 (left bottom panel). The limit of the gray region represents the height where the optical depth at the considered wing wavelength and LOS is unity. The top panel of Fig. 7 shows that the contributions to ε_I from R^III and R^III−CRD practically coincide at all heights, and that at the height where the optical depth is unity, they are very similar to that from R^II. By contrast, in the bottom panel, we see that the contribution to ε_Q of R^II dominates over that of R^III at all heights, and it can be thus considered the only relevant contributor to the formation of the Q/I wing lobe. In this case, the contribution from R^III−CRD is different from that of R^III, but it remains negligible with respect to that of R^II. This explains why the computations of Fig. 5 do not show any appreciable differences between R^III and R^III−CRD in the line wings.

Fig. 7. Various contributions (see legend) to the emission coefficients for the Stokes parameters I (top panel) and Q (bottom panel) as a function of height in the FAL-C model, in the absence of magnetic fields. The emission coefficients are evaluated at the wavelength λ = 4227.1 Å, corresponding to the maximum of the Q/I lobe in the red wing of the line, for a LOS with μ = 0.17. The shaded area in the panels highlights the atmospheric region where the optical depth at this wavelength and LOS is greater than 1. εX (where X= I or Q) is the total emissivity, ${\epsilon }_X^{\mathrm{II}}$, ${\epsilon }_X^{\mathrm{III}}$, and ${\epsilon }_X^{\mathrm{III}-\mathrm{CRD}}$ are the contributions from RII, RIII, and RIII−CRD, respectively, while ${\epsilon }_X^{\ell ,\mathrm{th}}$ and ${\epsilon }_X^c$ are the contributions from the line thermal emissivity and continuum, respectively.

Finally, Fig. 6 displays all the Stokes profiles in the presence of a height-independent, horizontal magnetic field with B = 20 G. An overall good agreement between R^III and R^III−CRD settings is observed in all profiles. We note that the error is generally larger in the core and near wings of the Q/I, U/I, and V/I profiles. This conclusion remains valid for magnetic fields ranging from 10 G to 200 G. Indeed, these results show no significant differences and are therefore not reported.

6.2. Sr I 4607 Å line

We now consider the modeling of the Sr I 4607 Å line, which is discretized with N_ν = 130 unevenly spaced nodes. Our non-LTE RT calculations adequately reproduce this weak line, showing a small absorption profile in the intensity spectrum without wings, and a prominent and sharp Q/I scattering polarization peak (see Fig. 8).

Fig. 8. Same as Fig. 5, but for the Sr I 4607 Å line.

In the I and Q/I profiles shown in Fig. 8 for a non-magnetic scenario, we only note minor differences between R^III and R^III−CRD calculations. By contrast, we note some relevant discrepancies between R^III and R^III−CRD cases when a deterministic magnetic field is considered. Figures 9–11 display the Q/I, U/I, and V/I profiles in the presence of a height-independent, horizontal magnetic field with B = 20 G, B = 50 G, and B = 100 G, respectively. We omit to report the intensity I profiles, because they are essentially identical to those exposed in Fig. 8 for the non-magnetized case. We first note that, in all cases, there are no observable discrepancies between R^III and R^III−CRD calculations in the U/I and V/I profiles. By contrast, some relevant differences appear in the Q/I profiles, such as the one shown in the top-left panel of Fig. 10. Moreover, R^III−CRD calculations generally present larger Q/I line-core signals at μ = 0.966 with respect to R^III ones. For completeness, Fig. 12 displays the Q/I profile where we observed the maximal error, corresponding to the B = 20 G case at μ = 0.83.

Fig. 9. Same as Fig. 6, but for the Sr I 4607 Å line.

Fig. 10. Same as Fig. 9, but for B = 50 G.

Fig. 11. Same as Fig. 9, but for B = 100 G.

Fig. 12. Same as Fig. 9, but for μ = 0.83. In this Q/I profile, we observed the maximal error between RIII and RIII−CRD calculations.

It is worth observing that for μ = 0.966 and a magnetic field with B = 20 G (i.e., a forward scattering Hanle effect geometry), the Ca I 4227 Å line shows a positive Q/I peak (see Fig. 6), while the Sr I 4607 Å line a negative one (see Fig. 9). Bearing in mind that both spectral lines originate from a J_ℓ = 0 → J_u = 1 transition and have similar Hanle critical fields, we may suggest that this inversion is due to their different formation heights and properties. An in-depth analysis of this result goes beyond the scope of this paper and will be the object of a future investigation.

As anticipated, the Sr I 4607 Å line has been extensively exploited to investigate the small-scale unresolved magnetic fields that fill the quiet solar photosphere. For this reason, we have also analyzed the case of unimodal MSI magnetic fields. We recall that in the presence of such magnetic fields, the signatures of Zeeman and magneto-optical effects vanish due to cancellation effects, and the only impact of the magnetic field is the depolarization of Q/I due to the Hanle effect, which depends on the field strength. Figure 13 shows the I and Q/I profiles for a height-independent MSI magnetic field with B = 20 G. No significant differences between R^III and R^III−CRD calculations are visible, and the error is always smaller than 5 × 10⁻³. This result suggests that the R^III−CRD approximation provides accurate results when the problem is characterized by cylindrical symmetry, while it can introduce appreciable discrepancies when the direction of a non-vertical, deterministic magnetic field breaks this symmetry.

Fig. 13. Same as Fig. 8, but in the presence of a uniform micro-structured isotropic magnetic field with B = 20 G.

7. Numerical results: 1D atmospheric model from 3D MHD simulation

In Sect. 6, we limited our calculations to the semi-empirical FAL-C atmospheric model, possibly including spatially uniform magnetic fields, in the absence of bulk velocities. In this section, we compare the impact of R^III and R^III−CRD in PRD calculations in a 1D atmospheric model extracted from a 3D MHD simulation, which includes height-dependent magnetic and bulk velocity fields. As we have seen in previous sections, the magnetic field impacts the polarization profiles through the Hanle, magneto-optical, and Zeeman effects. On the other hand, bulk velocities introduce Doppler shifts as well as amplitude enhancements and asymmetries in the polarization profiles (e.g., Sampoorna & Nagendra 2015, 2016; Carlin et al. 2012, 2013; Jaume Bestard et al. 2021). The joint impact of height-dependent magnetic fields and bulk velocities, together with the inherent thermodynamic structure of the atmospheric model, results in complex emergent Stokes profiles, in which it is generally difficult to distinguish the contribution of each factor.

In this work, we adopted a 1D atmospheric model extracted from the 3D magnetohydrodynamic (MHD) simulation of Carlsson et al. (2016), performed with the Bifrost code (Gudiksen et al. 2011). In practice, we took a vertical column of such 3D simulation, clipped in the height interval between −100 and 1400 km, which includes the region where the considered spectral lines form, discretized with N_z = 79 nodes. The corresponding vertical resolution ranges between 19 km and 100 km. Figure 14 shows the variation of the magnetic field, temperature, and vertical component of the bulk velocity as a function of height in the considered model (hereafter, Bifrost model), which shows relatively quiet conditions. Since the 1D module of the RH code (used to calculate the lower-level population) can only handle vertical bulk velocities, in this study we only considered the vertical component of the model’s velocity, although our code can in principle take into account velocities of arbitrary direction. As additional information, Fig. 15 shows the variation of the coherence fraction $\stackrel{\sim }{\alpha }$ with height in the Bifrost model for the two considered spectral lines, as well as the height at which the optical depth τ, in the frequency intervals of the two spectral lines, is unity (see Sect. 5).

Fig. 14. Physical quantities of the considered 1D atmospheric model extracted from the 3D MHD Bifrost simulation en024048_hion (snapshot 385, column 120 × 120). Upper panel: strength (blue solid line), inclination (gray dashed line), and azimuth (light-blue dashed-dotted line) of the magnetic field as a function of height. Lower panel: vertical component of bulk velocity (green solid line) and temperature (red solid line) as a function of height. We adopt the convention that the velocity is positive if pointing outward in the atmosphere and negative if pointing inward. For clarity, the horizontal green dotted line indicates zero velocity.

Fig. 15. Same as Fig. 4, but for the Bifrost atmospheric model.

7.1. Ca I 4227 Å line

We first consider the modeling of the Ca I 4227 Å line, which is now discretized with N_ν = 141 unevenly spaced nodes. Figure 16 displays all the Stokes profiles, showing an overall good agreement between R^III and R^III−CRD calculations for a LOS with μ = 0.034. By contrast, for a LOS with μ = 0.996, we observe some appreciable discrepancies in the Q/I and U/I profiles, with an error up to 3 × 10⁻¹. We recall that these linear scattering polarization signals, obtained for a LOS close to the disk center, are due to the forward-scattering Hanle effect (e.g., Trujillo Bueno 2001). Interestingly, such discrepancies disappear in the absence of bulk velocities (profiles not reported here). This suggests that the differences between R^III and R^III−CRD calculations are amplified in the presence of bulk velocities, and they are larger for a LOS close to the vertical because in this case the Doppler shifts are more pronounced. We finally note that the error affects the amplitudes of the main peaks of the profiles but not their shape.

Fig. 16. Emergent Stokes profiles for the Ca I 4227 Å line calculated in the Bifrost model (see Sect. 7), which includes height-dependent magnetic and (vertical) bulk velocity fields. The vertical gray lines show the line-center wavelength in the absence of bulk velocities.

7.2. Sr I 4607 Å line

We now consider the modeling of the Sr I 4607 Å line, which is discretized with N_ν = 141 unevenly spaced nodes. First, we note that in this Bifrost model, this line shows an emission profile in intensity for a LOS of μ = 0.034. We verified that this is due to the thermodynamic structure of this atmospheric model, and it can be explained from the behavior of the source function at the formation height of the line for this limb LOS (see Fig. 15). Figure 17 shows that appreciable discrepancies between R^III and R^III−CRD calculations are found in all Stokes profiles for μ = 0.034, and in Q/I and U/I for μ = 0.966. The maximal error is found in the U/I profile at μ = 0.966, where however the polarization signal is very weak and thus of limited practical interest. The discrepancies between R^III and R^III−CRD computations become slightly larger in the Bifrost model, where also a bulk velocity field is included. On the other hand, observing that the most significant errors appear in very weak polarization signals, we can conclude that for practical applications the R^III−CRD approximation can be safely used to obtain reliable results.

Fig. 17. Same as Fig. 16, but for the Sr I 4607 Å line.

8. Conclusions

In this work, we assessed the suitability and range of validity of a widely used approximation for R^III, that is, its expression under the assumption of CRD in the observer’s frame, labeled with R^III−CRD. To this aim, we solved the full non-LTE RT problem for polarized radiation in 1D models of the solar atmosphere, considering both the exact expression of R^III and its R^III−CRD approximation. With respect to the previous work of Sampoorna et al. (2017), we considered semi-empirical models, as well as 1D models extracted from 3D MHD simulations, which provide a reliable approximation of the solar atmosphere. The analysis was focused on the chromospheric Ca I line at 4227 Å and the photospheric Sr I line at 4607 Å, accounting for the impact of magnetic fields (both deterministic and micro-structured) and bulk velocities.

We first compared the analytical forms of R^III and R^III−CRD, showing that they correspond when the scattering angle Θ is equal to π/2, while their difference increases as Θ approaches 0 (forward scattering) or π (backward scattering). The numerical results for the Ca I 4227 Å line show that the R^III−CRD approximation is suited to model the scattering polarization signals of strong chromospheric lines. This is not surprising considering that in these lines the contribution of R^II dominates with respect to that of R^III, both in the core (which forms high in the atmosphere) and in the far wings (where scattering is basically coherent). On the other hand, we verified that the R^III−CRD approximation is also adequate in the near wings, where the contribution of R^III cannot be neglected a priori. The numerical results for the Sr I 4607 Å line show that the R^III−CRD approximation also provides reliable results in photospheric lines, for which the contribution of R^III is significant. Observing that the scattering polarization signal of the Sr I 4607 Å line is extensively used to investigate the small-scale magnetic fields that fill the quiet solar photosphere, we verified that the R^III−CRD approximation is also suitable in the presence of micro-turbulent magnetic fields. On the other hand, some appreciable discrepancies (i.e., relative errors larger than 0.01 according to the error definition (15)) are found when deterministic magnetic fields or bulk velocities are included, and the polarization signals are below 0.4%. However, the qualitative agreement between the two settings remains satisfactory in general; differences in the overall shapes and signs are exceptions and appear when the signals are very weak. These discrepancies may suggest that the limit of CRD (i.e., to fully neglect PRD effects), which is generally adopted to model the intensity and polarization of these weak lines, can possibly introduce some appreciable inaccuracies with respect to a PRD treatment. A quantitative comparison between PRD and CRD calculations for this line is ongoing, and it will be presented in a separate publication.

In conclusion, we can state that in optically thick media, the use of the lightweight R^III−CRD approximation guarantees reliable results in the modeling of scattering polarization in strong and weak resonance lines.

Acknowledgments

We thank the anonymous referee for carefully reading the manuscript and for many useful comments, suggestions, and corrections. The financial support by the Swiss National Science Foundation (SNSF) through grant CRSII5_180238 is gratefully acknowledged.

References

Appendix A: Atomic model and data

We consider an atomic system composed of two levels (two-level atom). Each level is characterized by the energy E, the quantum number for the total angular momentum J (positive integer or semi-integer values only), and the Landé factor g. Hereafter, the physical quantities referring to the upper and lower levels will be labeled with subscripts u and ℓ, respectively. Each level is composed of 2J + 1 magnetic sublevels characterized by the magnetic quantum number M (M = −J, −J + 1, ..., J). The magnetic sublevels are degenerate in the absence of magnetic fields, while they split in energy when a magnetic field B is present (Zeeman effect). Their energies are given by E(M) = E_J + h ν_L g M, where E_J is the energy of the considered J-level, h the Planck constant, and ν_L = eB(4πm_ec)⁻¹ (with e the elementary charge, m_e the electron mass, and c the speed of light) the Larmor frequency⁵. The line-center frequency is ν₀ = (E_u − E_ℓ)/h. The frequency of the Zeeman component corresponding to the transition between the upper magnetic sublevel M_u and the lower magnetic sublevel M_ℓ is

$$\begin{array}{c}\hfill {\nu }_{M_u{M}_{\ell }}=\frac{E(M_u)-E(M_{\ell })}{h}={\nu }_0+{\nu }_L(g_u{M}_u-g_{\ell }{M}_{\ell }).\end{array}$$

The energies, angular momenta, Landé factors, and Einstein coefficients for spontaneous emission A_uℓ for the two-level atomic models considered in this work to synthesize the intensity and polarization of the Ca I 4227 Å and Sr I 4607 Å lines are summarized in Table A.1. The Landé factors of levels with J_ℓ = 0 are formally taken to be zero.

Appendix B: Redistribution matrix

All the results presented in this work are obtained considering the redistribution matrix for a two-level atom with an unpolarized and infinitely sharp lower level, in the presence of magnetic fields, as derived by Bommier (1997b). This redistribution matrix is given by the sum of two terms, commonly referred to as R^II and R^III (e.g., Hummer 1962). In this appendix, we provide their analytic expressions, written with a slightly different notation than the one used in the original paper or in subsequent works in which the same redistribution matrices were reported and applied (e.g., Alsina Ballester et al. 2017). We first present their expressions in the atomic reference frame, taking the quantization axis along the magnetic field (Sect. B.1). We briefly comment on the branching ratios (Sect. B.1.1) and we derive their expressions in a new reference system with the quantization axis along the vertical (Sect. B.1.2). Subsequently, we transform them in the observer’s reference frame, including bulk velocities (Sect. B.2). Finally, we provide their expressions under a series of simplifying approximations, including the one for R^III analyzed in this work (Sect. B.3).

B.1. Expression in the atomic reference frame

In the atomic reference frame, taking the quantization axis parallel to the magnetic field (magnetic reference system), the R^II and R^III redistribution matrices are respectively given by Eqs. (51) and (49) in Bommier (1997b). We refer to these redistribution matrices with the symbol ${\widehat{R}}^{\mathrm{X}}$ (with X = II,III) to distinguish them from those expressed taking the quantization axis along a different direction. After some variable and index renaming,⁶ we rewrite them in the following equivalent form

$$\begin{array}{c}\hfill {\widehat{R}}_{\mathit{ij}}^{\mathrm{X}}(\mathit{r},\mathbf{\Omega },{\mathbf{\Omega }}^{\prime },\xi ,{\xi }^{\prime })={\displaystyle \sum _{K,K^{\prime }=0}^2}{\displaystyle \sum _{Q=-K_{\min }}^{K_{\min }}}{\mathcal{R}}_Q^{\mathrm{X},K{K}^{\prime }}(\mathit{r},\xi ,{\xi }^{\prime }){(-1)}^Q{\widehat{\mathcal{T}}}_{Q,i}^{K^{\prime }}(\mathit{r},\mathbf{\Omega }){\widehat{\mathcal{T}}}_{-Q,j}^K(\mathit{r},{\mathbf{\Omega }}^{\prime }),\end{array}$$(B.1)

where r is the spatial point, Ω the propagation direction, and ξ the radiation frequency in the atomic reference frame. The convention that primed and unprimed quantities refer to the incident and scattered radiation, respectively, is used. The indices i and j can take values 1, 2, 3, and 4, while K_min = min(K, K′). The quantity ${\widehat{\mathcal{T}}}_{Q,i}^K$ (with K = 0, ..., 2 and Q = −K, ..., K) is the geometrical tensor defined in Sect. 5.11 of LL04, evaluated in the magnetic reference system. Given that the direction of the magnetic field may vary with the spatial point, this tensor depends in general on r.

The function ${\mathcal{R}}_Q^{\mathrm{III},K{K}^{\prime }}$ is given by (see Eq. (49) of Bommier 1997b)

$$\begin{array}{c}\hfill {\mathcal{R}}_Q^{\mathrm{III},K{K}^{\prime }}(\mathit{r},\xi ,{\xi }^{\prime })={\displaystyle \sum _{K^{″}=|Q|}^{2J_u}}({\beta }_Q^K(\mathit{r})-{\alpha }_Q(\mathit{r})){\mathrm{\Phi }}_Q^{K^{″}{K}^{\prime }}(\mathit{r},\xi ){\mathrm{\Phi }}_Q^{K^{″}K}(\mathit{r},{\xi }^{\prime }),\end{array}$$(B.2)

where

$$\begin{array}{c}\hfill {\beta }_Q^K(\mathit{r})=\frac{{\mathrm{\Gamma }}_R}{{\mathrm{\Gamma }}_R+{\mathrm{\Gamma }}_I(\mathit{r})+D^{(K)}(\mathit{r})+2\pi \mathrm{i}{\nu }_L(\mathit{r})g_{u}Q},\end{array}$$(B.3)

and

$$\begin{array}{c}\hfill {\alpha }_Q(\mathit{r})=\frac{{\mathrm{\Gamma }}_R}{{\mathrm{\Gamma }}_R+{\mathrm{\Gamma }}_I(\mathit{r})+{\mathrm{\Gamma }}_E(\mathit{r})+2\pi \mathrm{i}{\nu }_L(\mathit{r})g_{u}Q}.\end{array}$$(B.4)

The quantities #x0393;_R, #x0393;_I, and #x0393;_E are the line broadening constants due to radiative decays, inelastic collisions, and elastic collisions, respectively:

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_R=A_{u\ell },{\mathrm{\Gamma }}_I=C_{u\ell },{\mathrm{\Gamma }}_E=Q_{\mathrm{el}},\end{array}$$

where C_uℓ is the rate of inelastic collisions inducing transitions from the upper to the lower level and Q_el is the rate of elastic collisions with neutral perturbers (mainly hydrogen and helium atoms). The quantity D^(K) is the depolarizing rate due to elastic collisions. In the absence of experimental data or detailed theoretical calculations for the rates D^(K), the approximate relation D⁽²⁾ = 0.5 Q_el is generally used. Under the assumption that the interaction between the atom and the perturber is described by a single tensor operator of rank K′, the rates with K ≠ 2 can be obtained from D⁽²⁾ as discussed in Sect. 7.13 of LL04 (see equation after (7.108) and Eq. (7.109), which relates D⁽¹⁾ and D⁽²⁾ for the case of K′ = 2). The generalized profile ${\mathrm{\Phi }}_Q^{K{K}^{\prime }}$ is given by Eq. (10.40) of LL04 and can be equivalently written as

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_Q^{K{K}^{\prime }}(\mathit{r},\xi )={\displaystyle \sum _{M_u,M_u^{\prime }=-J_u}^{J_u}}{\displaystyle \sum _{M_{\ell }=-J_{\ell }}^{J_{\ell }}}{\displaystyle \sum _{q,q^{\prime }=-1}^1}{\mathcal{B}}_{K{K}^{\prime }Q{M}_u{M}_u^{\prime }{M}_{\ell }q{q}^{\prime }}\frac{1}{2}[{\mathrm{\Phi }}_{M_u{M}_{\ell }}(\mathit{r},\xi )+{\overline{\mathrm{\Phi }}}_{M_u^{\prime }{M}_{\ell }}(\mathit{r},\xi )],\end{array}$$(B.5)

where K = 0, ..., 2J_u, K′ = 0, 1, 2, and Q = −K_min, ..., K_min with K_min = min(K, K′). The notation $\overline{f}(·)$ refers to the complex conjugate function. The quantity ℬ_{KK′QMuMu′Mℓqq′} is given by

$$\begin{array}{cc}\hfill {\mathcal{B}}_{K{K}^{\prime }Q{M}_u{M}_u^{\prime }{M}_{\ell }q{q}^{\prime }}& ={(-1)}^{1+J_u-M_u+q^{\prime }}\sqrt{3(2J_u+1)(2K+1)(2K^{\prime }+1)}\hfill \\ \hfill & \times \left(\begin{array}{ccc}{J}_u& J_{\ell }& 1\\ -M_u& M_{\ell }& -q\end{array}\right)\left(\begin{array}{ccc}{J}_u& J_{\ell }& 1\\ -M_u^{\prime }& M_{\ell }& -q^{\prime }\end{array}\right)\left(\begin{array}{ccc}{J}_u& J_u& K\\ M_u^{\prime }& -M_u& -Q\end{array}\right)\left(\begin{array}{ccc}1& 1& K^{\prime }\\ q& -q^{\prime }& -Q\end{array}\right),\hfill \end{array}$$

where the quantities in parentheses are the Wigner’s 3-j symbols (e.g., Sect. 2.2 of LL04). It must be observed that the 3-j symbols are nonzero only if the sum of the arguments of the lower raw is zero. Because of this, the values of q and q′ are uniquely determined once M_u, $M_u^{\prime }$, and M_ℓ are assigned. Thus, the sums over these indices in Eq. (B.5) are in practice redundant. The complex profile Φ_MuMℓ is defined as

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{M_u{M}_{\ell }}(\mathit{r},\xi )=\frac{1}{\pi }\frac{1}{\mathrm{\Gamma }(\mathit{r})-\mathrm{i}({\nu }_{M_u{M}_{\ell }}(\mathit{r})-\xi )},\end{array}$$(B.6)

with #x0393; = (#x0393;_R + #x0393;_I + #x0393;_E)/4π.

The function ${\mathcal{R}}_Q^{\mathrm{II},K{K}^{\prime }}$ is given by (see Eq. (51) of Bommier 1997b)

$$\begin{array}{cc}\hfill {\mathcal{R}}_Q^{\mathrm{II},K{K}^{\prime }}(\mathit{r},\xi ,{\xi }^{\prime })& ={\alpha }_Q(\mathit{r}){\displaystyle \sum _{M_u,M_u^{\prime }=-J_u}^{J_u}}{\displaystyle \sum _{M_{\ell },M_{\ell }^{\prime }=-J_{\ell }}^{J_{\ell }}}{\displaystyle \sum _{p,p^{\prime },p^{″},p^{″\prime }=-1}^1}{\mathcal{C}}_{K{K}^{\prime }Q{M}_u{M}_u^{\prime }{M}_{\ell }{M}_{\ell }^{\prime }p{p}^{\prime }{p}^{″}{p}^{″\prime }}\hfill \\ \hfill & \times \delta (\xi -{\xi }^{\prime }-{\nu }_{M_{\ell }{M}_{\ell }^{\prime }}(\mathit{r}))\frac{1}{2}[{\mathrm{\Phi }}_{M_u^{\prime }{M}_{\ell }}(\mathit{r},{\xi }^{\prime })+{\overline{\mathrm{\Phi }}}_{M_u{M}_{\ell }}(\mathit{r},{\xi }^{\prime })],\hfill \end{array}$$(B.7)

where α_Q is given by Eq. (B.4), δ(⋅) is the Dirac delta, and Φ_MuMℓ is the complex profile defined in Eq. (B.6). The quantity ${\nu }_{M_{\mathcal{l}}{M}_{\mathcal{l}}^{\prime }}$ is given by

$$\begin{array}{c}\hfill {\nu }_{M_{\ell }{M}_{\ell }^{\prime }}(\mathit{r})={\nu }_L(\mathit{r})(g_{\ell }{M}_{\ell }-g_{\ell }{M}_{\ell }^{\prime }),\end{array}$$

and ${\mathcal{C}}_{K{K}^{\prime }Q{M}_u{M}_u^{\prime }{M}_{\ell }{M}_{\ell }^{\prime }pp\prime p″p^{″\prime }}$ is defined as (see Eq. (12) of Bommier 1997b)

$$\begin{array}{cc}\hfill {\mathcal{C}}_{K{K}^{\prime }Q{M}_u{M}_u^{\prime }{M}_{\ell }{M}_{\ell }^{\prime }p{p}^{\prime }{p}^{″}{p}^{″\prime }}& =3(2J_u+1)\sqrt{2K+1}\sqrt{2K^{\prime }+1}{(-1)}^{2J_u-M_{\ell }-M_{\ell }^{\prime }}\hfill \\ \hfill & \times \left(\begin{array}{ccc}{J}_u& J_{\ell }& 1\\ M_u& -M_{\ell }& -p\end{array}\right)\left(\begin{array}{ccc}{J}_u& J_{\ell }& 1\\ M_u^{\prime }& -M_{\ell }& -p^{\prime }\end{array}\right)\left(\begin{array}{ccc}{J}_u& J_{\ell }& 1\\ M_u& -M_{\ell }^{\prime }& -p^{″}\end{array}\right)\left(\begin{array}{ccc}{J}_u& J_{\ell }& 1\\ M_u^{\prime }& -M_{\ell }^{\prime }& -p^{″\prime }\end{array}\right)\hfill \\ \hfill & \times \left(\begin{array}{ccc}1& 1& K\\ -p& p^{\prime }& Q\end{array}\right)\left(\begin{array}{ccc}1& 1& K^{\prime }\\ -p^{″}& p^{″\prime }& Q\end{array}\right).\hfill \end{array}$$

As for q and q′ in Eq. (B.5), the sums over p, p′, p″, and p‴ in Eq. (B.7) are in practice redundant.

B.1.1. Branching ratios

The quantities $({\beta }_Q^K-{\alpha }_Q)$ and α_Q appearing in Eqs. (B.2) and (B.7), respectively, are the branching ratios between R^II and R^III. These terms also contain the branching ratio for the scattering contribution to the total emissivity, 1 − ϵ, with

$$\begin{array}{c}\hfill ϵ(\mathit{r})=\frac{{\mathrm{\Gamma }}_I(\mathit{r})}{{\mathrm{\Gamma }}_R+{\mathrm{\Gamma }}_I(\mathit{r})}\end{array}$$(B.8)

the photon destruction probability. Factorizing this term (i.e., writing ${\alpha }_Q=(1-ϵ){\stackrel{\sim }{\alpha }}_Q$ and ${\beta }_Q^K=(1-ϵ){\stackrel{\sim }{\beta }}_Q^K$), the net branching ratios for R^II and R^III are therefore