© ESO 2021

1. Introduction

When observing quiet regions close to the edge of the solar disk (limb), several strong resonance lines, such as H I Ly-α, Ca II K, Ca I 4227, or Na ID₂, show significant scattering polarization signals, characterized by a relatively sharp peak in the line core and broader lobes in the wings (e.g., Gandorfer 2000, 2002; Kano et al. 2017). The observation and modeling of these polarization signals have received increased attention by the scientific community. In fact, their sensitivity to the presence of relatively weak magnetic fields is presently one of the most promising tools to probe the magnetism of the solar chromosphere and transition region (e.g., Trujillo Bueno et al. 2017). This magnetic sensitivity arises from the combined action of the Hanle, Zeeman, and magneto-optical effects (e.g., Alsina Ballester et al. 2017). The interest in these signals is clearly attested by the recent Chromospheric Ly-Alpha SpectroPolarimeter (CLASP) and Chromospheric LAyer SpectroPolarimeter (CLASP2) sounding rocket experiments, which provided unprecedented spectropolarimetric observations of the H I Ly-alpha line (Kano et al. 2017; Trujillo Bueno et al. 2018) and of the Mg II h and k lines (Ishikawa et al. 2021), respectively. In the visible and infrared ranges of the solar spectrum, scattering polarization signals of strong chromospheric lines can be observed from the ground with instruments such as the Zurich IMaging POLarimeter (ZIMPOL, Ramelli et al. 2010) or the Tenerife Infrared Polarimeter II (TIP-II, Collados et al. 2007), and will be among the targets of the new generation of large-aperture solar telescopes (e.g., DKIST and EST).

The scattering polarization signals of strong resonance lines are primarily produced by scattering processes that are coherent in frequency in the atomic reference frame, and they are subject to Doppler redistribution in the observer’s frame. The core of these signals can often be suitably modeled under the limit of complete frequency redistribution (CRD; e.g., Sampoorna et al. 2010; Štěpán et al. 2015). However, this approximation turns out to be completely inadequate to model their extended wing lobes, for which a general approach accounting for partial frequency redistribution (PRD) effects is of the utmost importance. The numerical solution of the radiative transfer (RT) problem for polarized radiation under nonlocal thermodynamic equilibrium (NLTE) conditions, taking scattering polarization and PRD effects into account, is a notoriously challenging problem from a computational point of view. For this reason, the so-called angle-averaged (AA) simplifying assumption, first introduced in the unpolarized case (e.g., Mihalas 1978), is still commonly applied to the modeling of scattering polarization signals (see, e.g., Rees & Saliba 1982; Belluzzi & Trujillo Bueno 2014; Alsina Ballester et al. 2017). By smearing out geometrical aspects of the problem, the AA assumption can, however, introduce significant inaccuracies, and it would be completely inadequate when aimed at modeling the aforementioned signals in realistic atmospheric models, which consider the detailed three-dimensional (3D) structure of the solar atmosphere.

This highlights the need for methods that, exploiting the increased computational resources available today, make it possible to model these signals by solving the NLTE RT problem for polarized radiation, taking PRD effects into account in the general angle-dependent (AD) case. In fact, considerable efforts have already been made in this direction. Dumont et al. (1977) pioneered the modeling of spectral line polarization accounting for AD PRD effects. Subsequently, Nagendra et al. (2002) considered the AD PRD problem in the weak-field Hanle regime for the case of a two-level atom, using the theoretical approach of Stenflo (1994), based on the Kramers-Heisenberg formula. Sampoorna et al. (2008, 2017) generalized that work to the Hanle-Zeeman regime considering arbitrary magnetic fields. Sampoorna et al. (2011) and Supriya et al. (2012) discussed the validity of the AA assumption in the presence of micro-turbulent magnetic fields, while Supriya et al. (2013) extended the analysis to vertical magnetic fields. Sampoorna & Nagendra (2015) developed an AD PRD RT code capable of including vertical bulk velocities. For computational simplicity, the AD PRD calculations presented in the aforementioned works were performed considering one-dimensional (1D) homogeneous atmospheric models. More recently, del Pino Alemán et al. (2020) solved the AD PRD problem for a three-term atomic model in more realistic 1D atmospheric models, applying the theoretical approach of Casini et al. (2014, 2017). Their code can take into account the effect of AD PRD in the presence of magnetic fields of an arbitrary strength and orientation, as well as vertical bulk velocities. However, their AD PRD calculations were limited to problems with axial symmetry because of the extreme computational cost.

This paper shows the results obtained with a new code capable of solving the NLTE RT problem in static 1D semi-empirical models of the solar atmosphere, taking scattering polarization and AD PRD effects into account. Based on the theoretical approach developed by Bommier (1997a,b), the code allows one to consider two-level atomic models with an unpolarized and infinitely sharp lower level in the presence of arbitrary magnetic fields. The present investigation is specifically focused on the Ca I 4227 Å line, which is an ideal test bench for new approaches to the modeling of scattering polarization including PRD effects. Indeed, the scattering polarization signal of this line clearly shows the impact of such effects (Faurobert-Scholl 1992) and can be suitably modeled by considering a simple two-level atom. Furthermore, the Stokes profiles of this line are of high scientific interest for magnetic field diagnostics, and high-precision spectropolarimetric observations are available (e.g., Bianda et al. 2011; Supriya et al. 2014; Alsina Ballester et al. 2018; Capozzi et al. 2020).

The structure of this paper is as follows: Sect. 2 outlines the formulation of the problem, exposing the relevant continuous equations and their discretization and presenting a suitable iterative numerical approach to solve it. Next, Sect. 3 focuses on the modeling of the Ca I 4227 Å line, presenting a quantitative comparison between the results of AA and AD calculations. Finally, Sect. 4 summarizes the main results of this paper and discusses their implications.

2. Formulation of the problem and methods

In this section, we first present the continuous formulation of the considered NLTE RT problem for polarized radiation, including AD PRD effects. We then expose the discretization and the algebraic formulation of the problem and outline the numerical iterative solution method that was applied. In addition, we finally recall the simplifying AA approximation.

2.1. Continuous formulation of the problem

The physical quantities entering the problem are, in general, functions of the spatial point r, and of the frequency ν and propagation direction Ω of the radiation beam under consideration. A complete description of the polarization state of the radiation field is provided by the four Stokes parameters

$$\begin{array}{c}\hfill I_i(\mathit{r},\mathbf{\Omega },\nu ),\end{array}$$

with i = 1, …, 4 standing for Stokes I, Q, U, and V, respectively. The transfer of partially polarized light along direction Ω at frequency ν is described by the system of coupled first-order inhomogeneous ordinary differential equations given by

$$\begin{array}{c}\hfill \frac{\mathrm{d}}{\mathrm{d}s}{I}_i(\mathit{r},\mathbf{\Omega },\nu )=-{\displaystyle \sum _{j=1}^4}{K}_{\mathit{ij}}(\mathit{r},\mathbf{\Omega },\nu )I_j(\mathit{r},\mathbf{\Omega },\nu )+{\epsilon }_i(\mathit{r},\mathbf{\Omega },\nu ),\end{array}$$(1)

where s is the spatial coordinate measured along the direction Ω, K_ij are the elements of the so-called propagation matrix

$$\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}$$

where η_i and ρ_i are the dichroism coefficients and the anomalous dispersion coefficients, respectively, while ε_i are the emission coefficients in the four Stokes parameters. In general, K_ij and ε_i contain contributions due to both line and continuum processes and can be expressed as

$$\begin{array}{cc}\hfill K_{\mathit{ij}}(\mathit{r},\mathbf{\Omega },\nu )& =K_{\mathit{ij}}^{\ell }(\mathit{r},\mathbf{\Omega },\nu )+K_{\mathit{ij}}^c(\mathit{r},\nu ),\hfill \end{array}$$(2)

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

with the superscripts ℓ and c standing for line and continuum, respectively.

In this work, we consider a two-level atomic model with an unpolarized and infinitely sharp lower level in the presence of magnetic fields of an arbitrary intensity and orientation. For such an atomic model, an analytical solution of the statistical equilibrium equations is available and the line contribution to the emission coefficient can thus be expressed through the redistribution matrix formalism. Stimulated emission is negligible around the Ca I 4227 Å frequency in the solar atmosphere and is consequently neglected. Moreover, bulk velocity fields are not considered; their impact on this problem will be the subject of an upcoming investigation.

2.1.1. Line contributions

The line contributions to η_i and ρ_i for the considered atomic model are given by (see Alsina Ballester et al. 2017)

$$\begin{array}{cc}\hfill {\eta }_i^{\ell }(\mathit{r},\mathbf{\Omega },\nu )& =k_L(\mathit{r}){\displaystyle \sum _{K=0}^2}{\widehat{\mathcal{T}}}_{0,i}^K(\mathit{r},\mathbf{\Omega }){\mathrm{\Phi }}_0^{0K}(\mathit{r},\nu ),\hfill \end{array}$$(4)

$$\begin{array}{cc}\hfill {\rho }_i^{\ell }(\mathit{r},\mathbf{\Omega },\nu )& =k_L(\mathit{r}){\displaystyle \sum _{K=0}^2}{\widehat{\mathcal{T}}}_{0,i}^K(\mathit{r},\mathbf{\Omega }){\mathrm{\Psi }}_0^{0K}(\mathit{r},\nu ),\hfill \end{array}$$(5)

where k_L is the frequency-integrated absorption coefficient, which depends on the population of the lower level, ${\widehat{\mathcal{T}}}_{Q,i}^K$ is the polarization tensor (see Landi Degl’Innocenti & Landolfi 2004, Chapt. 5) evaluated in the magnetic reference system¹, while ${\mathrm{\Phi }}_Q^{K{K}^{\prime }}$ and ${\mathrm{\Psi }}_Q^{K{K}^{\prime }}$ are the generalized profiles defined in Landi Degl’Innocenti & Landolfi (2004, Appendix 13).

Within the redistribution matrix formalism, one can distinguish between thermal and scattering contributions to the line emission coefficients, namely,

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

Assuming that inelastic collisions are isotropic, the thermal contribution to the line emission coefficients for the considered atomic system reads as follows (see Alsina Ballester et al. 2017):

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\epsilon }_i^{\ell ,th}(\mathit{r},\mathbf{\Omega },\nu )& =k_L(\mathit{r})ϵ(\mathit{r})W_{\mathrm{T}}(\mathit{r}){\displaystyle \sum _{K=0}^2}{\widehat{\mathcal{T}}}_{0,i}^K(\mathbf{\Omega }){\mathrm{\Phi }}_0^{0K}(\mathit{r},\nu )\hfill \\ \hfill & =ϵ(\mathit{r})W_{\mathrm{T}}(\mathit{r}){\eta }_i^{\ell }(\mathit{r},\mathbf{\Omega },\nu ),\hfill \end{array}\end{array}$$(6)

where ϵ is the thermalization parameter and W_T is the Planck function in the Wien limit at the line-center frequency ν₀. The scattering contribution to the line emission coefficients reads as follows:

$$\begin{array}{c}\hfill {\epsilon }_i^{\ell ,sc}(\mathit{r},\mathbf{\Omega },\nu )=k_L(\mathit{r}){\displaystyle \int \mathrm{d}\nu \prime \oint \frac{\mathrm{d}\mathbf{\Omega }\prime }{4\pi }\sum _{j=1}^4{R}_{\mathit{ij}}(\mathit{r},\mathbf{\Omega },\mathbf{\Omega }\prime ,\nu ,\nu \prime )I_j(\nu \prime ,\mathbf{\Omega }\prime ,\mathit{r}),}\end{array}$$(7)

where R_ij are the elements of the redistribution matrix and I_j are the Stokes parameters of the radiation that pumps the atomic system. Hereafter, the primed and unprimed quantities refer to the incoming and outgoing radiation, respectively. The redistribution matrix derived by Bommier (1997b) is given by the sum of two terms

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

where R^II and R^III describe scattering processes that are coherent and completely uncoherent, respectively, in the atomic reference frame. The transformation of R^II of Bommier (1997b) from the atomic to the observer’s frame (accounting for the Doppler effect due to thermal motions) is described in detail in Alsina Ballester et al. (2017). The assumption of CRD in the observer’s frame was instead made for R^III (e.g., Bommier 1997b; Alsina Ballester et al. 2017). This approximation appears to be suitable for modeling strong chromospheric lines, whose core forms at atmospheric heights where the rate of elastic collisions is relatively low, and the contribution of R^III is therefore much smaller than that of $R^{\mathrm{II}}$. A quantitative analysis of the impact of this approximation on our modeling goes beyond the scope of this work and is left for a future investigation. The reader is referred to Bommier (1997b) and Sampoorna et al. (2017) for more details on the suitability of this assumption. The ensuing redistribution matrices can be written as

$$\begin{array}{cc}\hfill R_{\mathit{ij}}^{\mathrm{II}}(\mathit{r},\mathbf{\Omega },\mathbf{\Omega }\prime ,\nu ,\nu \prime )& ={\displaystyle \sum _{KK\prime Q}}{\mathcal{R}}_Q^{\mathrm{II},KK\prime }(\mathit{r},\nu ,\nu \prime ,\mathrm{\Theta })P_{Q,ij}^{KK\prime }(\mathbf{\Omega },\mathbf{\Omega }\prime ),\hfill \\ \hfill R_{\mathit{ij}}^{\mathrm{III}}(\mathit{r},\mathbf{\Omega },\mathbf{\Omega }\prime ,\nu ,\nu \prime )& ={\displaystyle \sum _{KK\prime Q}}{\mathcal{R}}_Q^{\mathrm{III},KK\prime }(\mathit{r},\nu ,\nu \prime )P_{Q,ij}^{KK\prime }(\mathbf{\Omega },\mathbf{\Omega }\prime ),\hfill \end{array}$$

where Θ represents the scattering angle, that is the angle between the incoming Ω′ and outgoing Ω directions. The redistribution functions ${\mathcal{R}}_Q^{\mathrm{II},KK\prime }$ and ${\mathcal{R}}_Q^{\mathrm{III},KK\prime }$ and the elements of the scattering phase matrix $P_{Q,ij}^{KK\prime }$ are also described in Alsina Ballester et al. (2017).

2.1.2. Continuum contributions

In the visible part of the solar spectrum, continuum processes significantly contribute to the diagonal element η_I of the propagation matrix alone, namely,

$$\begin{array}{cc}\hfill {\eta }_i^c(\mathit{r},\nu )& =k_c(\mathit{r},\nu ){\delta }_{i1},\hfill \\ \hfill {\rho }_i^c& =0,\hfill \end{array}$$

where δ_ij is the Kronecker delta. The continuum contribution to the emission coefficients can be divided into a thermal term, which is unpolarized and isotropic, and a scattering term

$$\begin{array}{c}\hfill {\epsilon }_i^c(\mathit{r},\mathbf{\Omega },\nu )={\epsilon }_I^{c,th}(\mathit{r},\nu ){\delta }_{i1}+{\epsilon }_i^{c,sc}(\mathit{r},\mathbf{\Omega },\nu ).\end{array}$$(8)

The scattering contribution to the continuum emission coefficients reads as follows (see Alsina Ballester et al. 2017):

$$\begin{array}{c}\hfill {\epsilon }_i^{c,sc}(\mathit{r},\mathbf{\Omega },\nu )=\sigma (\mathit{r},\nu ){\displaystyle \sum _{\mathit{KQ}}}{(-1)}^Q{\mathcal{T}}_{Q,i}^K(\mathbf{\Omega })\oint \frac{\mathrm{d}\mathbf{\Omega }\prime }{4\pi }\sum _{j=1}^4{\mathcal{T}}_{-Q,j}^K(\mathbf{\Omega })I_j(\mathit{r},\mathbf{\Omega },\nu ),\end{array}$$(9)

where σ is the continuum opacity for scattering and ${\mathcal{T}}_{Q,i}^K$ is the polarization tensor (see Landi Degl’Innocenti & Landolfi 2004, Chapt. 5) evaluated in the reference system of the problem².

In summary, in order to determine the Stokes parameters I_i of the radiation emerging from the considered atmospheric model, one needs to solve the transfer Eq. (1). This requires knowledge on the emission vector ε_i which, in turn, depends on the Stokes parameters I_i through (7) and (9). The problem thus consists in finding a self-consistent solution for I_i and ε_i.

2.2. Discretization and algebraic formulation

After a suitable discretization of the continuous problem (see, e.g., Janett et al. 2021), the solution of the transfer Eq. (1) and the calculation of the emission coefficients (3) can be expressed in the compact matrix form

$$\begin{array}{ccc}\hfill \mathbf{I}& =\mathrm{\Lambda }\mathit{\epsilon }+\mathbf{t},\hfill & \hfill \text{with}\mathrm{\Lambda }\in {\mathbb{R}}^{4N_s{N}_{\mathrm{\Omega }}{N}_{\nu }\times 4N_s{N}_{\mathrm{\Omega }}{N}_{\nu }},\end{array}$$(10)

$$\begin{array}{ccc}\hfill \mathit{\epsilon }& =\mathrm{\Sigma }\mathbf{I}+\mathbf{c},\hfill & \hfill \text{with}\mathrm{\Sigma }\in {\mathbb{R}}^{4N_s{N}_{\mathrm{\Omega }}{N}_{\nu }\times 4N_s{N}_{\mathrm{\Omega }}{N}_{\nu }},\end{array}$$(11)

respectively, with the vectors I, ε, t, c$\in {\mathbb{R}}^{4N_s{N}_{\mathrm{\Omega }}{N}_{\nu }}$. We note that I collects the discretized Stokes parameters and ε the discretized emission coefficients. Moreover, Λε describe the transfer of the radiation generated inside the atmospheric model, t represents the radiation transmitted from the boundaries, while ΣI and c represent the scattering and thermal contributions to the emission coefficients, respectively.

By choosing I as the problem’s unknown, Eqs. (10)–(11) can then be combined in a single discrete problem, namely,

$$\begin{array}{c}\hfill (Id-\mathrm{\Lambda }\mathrm{\Sigma })\mathbf{I}=\mathrm{\Lambda }\mathbf{c}+\mathbf{t},\end{array}$$(12)

where Id is the identity matrix with a size of 4N_sN_ΩN_ν. The discrete problem (12) is generally nonlinear, because the coefficients of Λ and t depend, in a nonlinear way, on the absorption coefficient η_I, which depends on the lower level atomic population, which in turn is affected by the unknown I. However, if the lower level atomic population is fixed a priori, the absorption coefficient η_I becomes a constant of the problem. In this case, Λ and t become independent from the unknown I and the discrete problem (12) becomes linear.

2.3. Numerical iterative solution

The numerical solution of the full NLTE RT problem is divided into two main steps. In the first step, the nonlinear NLTE RT problem is first solved considering a multilevel atomic model, taking PRD effects into account, but neglecting polarization and the effect of magnetic fields. An accurate estimate of the populations of the atomic levels is thus provided. In the second step, a two-level atomic model with an unpolarized and infinitely sharp lower level is considered. The population of the lower level obtained in the first step is used and kept fixed. The resulting linear problem (12) is solved in the general AD case. Although this problem can in principle be solved in a single step with a direct method, this task is beyond current computational capabilities. A stationary or Krylov iterative method can therefore be applied (see, e.g., Janett et al. 2021; Benedusi et al. 2021), using the unpolarized radiation field obtained in the first step as an initial guess.

Although it is not fully consistent, this two-step approach permits one to obtain an accurate estimate of the lower level population from a multilevel atomic model, while retaining the simplicity of a two-level atom when accounting for polarization. Moreover, by fixing the lower level population, the line absorption coefficient becomes a constant, and the problem at step 2 becomes linear (see Sect. 2.2). It is worth noting that the population of the upper level is not kept fixed at step 2. This allows one to take into account the impact of atomic polarization (neglected at step 1) on the ratio between the population of the upper and lower level. In this respect, it must be observed that the degree of atomic polarization is in general quite low and, therefore, this impact is in most cases very small. The validity of the approach is corroborated by previous works (e.g., Belluzzi & Trujillo Bueno 2012; Belluzzi et al. 2012; Alsina Ballester et al. 2016, 2018).

2.4. Angle-averaged approximation

In the general AD approach, one considers the inherent coupling between incoming and outgoing directions and frequencies, which manifests in the presence of the scattering angle Θ in the redistribution function ${\mathcal{R}}_Q^{\mathrm{II},KK\prime }$. Accurate AD calculations are computationally very expensive, and the AA approximation (e.g., Mihalas 1978; Rees & Saliba 1982) was thus proposed to simplify the numerical problem. According to the formulation outlined in Sect. 2.1, this approximation consists in averaging the expression of ${\mathcal{R}}_Q^{\mathrm{II},KK\prime }$ over the scattering angle Θ, namely,

$$\begin{array}{c}\hfill {\overline{\mathcal{R}}}_Q^{\mathrm{II},KK\prime }(\mathit{r},\nu ,\nu \prime )=\frac{1}{2}{\displaystyle {\int }_0^{\pi }}\mathrm{d}\mathrm{\Theta }sin(\mathrm{\Theta }){\mathcal{R}}_Q^{\mathrm{II},KK\prime }(\mathit{r},\mathrm{\Theta },\nu ,\nu \prime ).\end{array}$$

The AA redistribution function ${\overline{\mathcal{R}}}_Q^{\mathrm{II},KK\prime }$ is thus independent of the scattering angle, partially decoupling (in the absence of bulk velocities) frequencies and directions in the calculation of the line emissivity (see Eq. (7)). This allows for a drastic reduction in the computational cost (e.g., Alsina Ballester et al. 2017).

This approximation has been proved suitable to model the intensity profiles (e.g., Leenaarts et al. 2012; Sukhorukov & Leenaarts 2017). In the modeling of polarization signals, the range of validity of this approximation was analyzed by Faurobert (1987, 1988) for the nonmagnetic case, by Sampoorna et al. (2011) for weak magnetic fields, and by Sampoorna et al. (2017) for arbitrary magnetic fields. The overall conclusion is that the AA approximation does not allow for a satisfactory modeling of the linear polarization signals in the core of spectral lines. These analyses have been, however, limited to academic scenarios, with hypothetical lines and homogeneous atmospheric models. A rigorous comparison between the AA and the AD approaches for realistic atomic and atmospheric models is still lacking and is the main topic of the next section.

3. Numerical results

This section presents the numerical modeling of the linear scattering polarization profiles of the Ca I line at 4227 Å, as an illustrative application of the problem’s formulation, discretization, and solution method described in Sect. 2. The calculations were carried out in 1D semi-empirical models of the solar atmosphere, in the presence of arbitrary magnetic fields, and taking AD PRD effects into account.

3.1. The Ca I 4227 Å spectral line

The Ca I 4227 Å line is produced by the transition between the ground level of neutral calcium (4s^2 1S₀) and the excited level $4\text{s}4\text{p}{1}^{}{\text{P}}_1^{\text{o}}$. Since both levels pertain to singlets and the upper one is not radiatively connected to other levels with lower energy, the scattering polarization signal of this line can be suitably modeled by considering a two-level atomic model.

The intensity spectrum of this line shows extended damping wings and a nearly saturated line core. For a disk center line of sight (LOS) and considering the semi-empirical model C of Fontenla et al. (1993), the line core forms at a height of about 900 km above the τ = 1 surface, with τ being the continuum optical depth at 5000 Å along the vertical. When observed in quiet regions close to the limb, this line exhibits the largest scattering polarization signal in the visible range of the so-called second solar spectrum (e.g., Stenflo et al. 1980; Stenflo 1982; Gandorfer 2002). This signal presents a clear triple-peak structure, typical of strong resonance lines, with a sharp peak in the line core and broad lobes in the wings. The wing lobes form at photospheric heights (at around 180 km for an observation at μ = cosθ = 0.1, with θ being the heliocentric angle), and they can only be suitably modeled by taking PRD effects into account. The scattering polarization signal of this line is of particular interest because of its sensitivity to the magnetic fields present across a wide range of atmospheric heights. Indeed, the line-core peak is sensitive to low-chromospheric magnetic fields through the Hanle effect (Stenflo 1982; Faurobert-Scholl 1992), while the wing lobes are sensitive to photospheric fields via magneto-optical effects (Alsina Ballester et al. 2018). Furthermore, the well-known Zeeman effect can be exploited as well in this line.

Considerable effort was indeed directed to observing and modeling the scattering polarization signal of the Ca I 4227 Å line. It is worthwhile to mention the observations of Bianda et al. (2003), which showed unexpected U/I wing signals, as well as clear spatial variations in the amplitude of both the Q/I and U/I wing lobes. These observations could be explained in terms of magneto-optical effects by Alsina Ballester et al. (2018). This theoretical interpretation has gained further support from the observational work by Capozzi et al. (2020). In the core of this line, scattering polarization signals can also be observed at the disk center, the interpretation of which requires 3D radiative transfer modeling (see Jaume Bestard et al. 2021). Observations of such disk-center signals were presented by Bianda et al. (2011) and subsequently analyzed by Anusha et al. (2011) to infer the magnetic fields present in the low chromosphere under the assumption that the observed signals are due to the Hanle effect in forward-scattering geometry (see Trujillo Bueno 2001). More recently, Carlin & Bianda (2017) theoretically studied the temporal evolution of these polarization signals using the CRD and 1.5D approximations, emphasizing their sensitivity to gradients in the vertical component of the bulk velocities, while Jaume Bestard et al. (2021) investigated the linear and circular polarization signals of the Ca I 4227 Å line in full 3D, pointing out the important symmetry breaking produced by the spatial gradients of the horizontal component of the macroscopic velocities. Moreover, Supriya et al. (2014) modeled the center-to-limb variation in the Q/I profile of this line, and they compared it to the observations recorded with the ZIMPOL-III instrument (e.g., Ramelli et al. 2010).

3.2. Numerical setting and methods

The Stokes profiles of the Ca I 4227 Å line presented in the following sections were calculated in the 1D semi-empirical models of Fontenla et al. (1993, hereafter FAL models), whose spatial grids contain between 60 and 70 nodes. A Cartesian reference system with the z-axis (quantization axis for total angular momentum) directed along the vertical is considered. A given propagation direction Ω is specified through its inclination θ ∈ [0, π] (i.e., μ = cosθ ∈ [ − 1, 1]), with respect to the vertical, and azimuth χ ∈ [0, 2π) measured counterclockwise from the x-axis. The angular grid is given by the tensor product of eight equally spaced nodes in the azimuthal interval [0, 2π], and six Gauss-Legendre quadrature nodes in each inclination subinterval [ − 1, 0] and [0, 1], for a total of 96 nodes. The wavelength interval [4071 Å, 4261 Å] is discretized onto a spectral grid with 155 nodes, equally spaced in the line core and logarithmically distributed in the wings. In all calculations, the reference direction for positive Stokes Q is considered to be parallel to the limb.

The whole problem was solved through the two-step approach presented in Sect. 2.3. The first step was carried out by means of the RH code (Uitenbroek 2001), using the parabolic short-characteristic formal solver. In step 1, an atomic model for calcium with 25 levels (19 Ca I levels, five Ca II levels, and the ground level of Ca III) was considered. This model accounts for 21 line transitions and 24 continuum transitions. All spectral lines were computed in the CRD limit, except for the Ca I 4227 Å line and the Ca II H & K lines, which received a PRD treatment. The RH code was also used to calculate the continuum total opacity k_c, scattering opacity σ, and thermal emissivity ${\epsilon }_I^{\text{c},\text{th}}$, as well as the collisional rates. The rates for inelastic collisions with electrons, inducing transitions between the upper and lower levels of the Ca I 4227 Å line, were computed according to the formula given in Seaton (1962). The rates of elastic collisions were calculated including the contributions from the Van der Waals interaction with neutral hydrogen and helium atoms, computed following Unsöld (1955), and from the quadratic Stark effect due to electrons and singly charged ions (Traving 1960; Gray 1992).

Step 2 was carried out by means of a Matlab code that solves the linear problem (12) in the general AD case with a Richardson iterative method (see Janett et al. 2021), using the unpolarized radiation field calculated at step 1 as an initial guess. In step 2, a two-level atomic model was considered. The population of the lower level was fixed to the value calculated in the first step and left unchanged. The DELO-parabolic formal solver (Janett et al. 2017, 2018) was used to solve the transfer Eq. (1).

Concerning the choice of the formal solver, it must be recalled that high-order schemes generally outperform low-order methods if the atmospheric model guarantees sufficient smoothness in the solution. This turns out to be the case for the FAL models, and for this reason, the third-order accurate DELO-parabolic method of Janett et al. (2017) was used to perform the calculations presented in this paper. Even though Janett et al. (2018) did not recommend multistep schemes for practical applications, the smoothness of the FAL atmospheric models allows for the safe use of the DELO-parabolic method, guaranteeing its proximity to L-stability (see Janett & Paganini 2018).

In order to facilitate the convergence in step 2, an intermediate step can be added. In this additional step, the linear problem (12) was solved with a stationary iterative method under the AA approximation (see Alsina Ballester et al. 2017). This calculation provides a temporary solution for the polarized radiation field, which acts as a better initial guess for step 2. At the same time, the Stokes profiles resulting from the intermediate step also provides a comparison suitable to identify the impact and the range of validity of the AA approximation with respect to the general AD treatment of PRD effects. These AA calculations were carried out using a modified version of the Fortran code developed by Alsina Ballester et al. (2017), which solves the linear problem (12) taking PRD effects into account under the AA approximation. Unlike the original version of the code, this modified version uses a damped Jacobi iterative method³ and the DELO-parabolic formal solver. The damped Jacobi method is used because the Jacobi iterative method does not guarantee convergence when used in combination with the DELO-parabolic formal solver (see Janett et al. 2021). A detailed analysis on the convergence properties of various iterative methods applied to transfer problems of polarized radiation can be found in Janett et al. (2021).

3.3. Nonmagnetic case

The impact of the AD treatment on the modeling of the emergent Stokes profiles of the Ca I 4227 Å line was first analyzed in the absence of magnetic fields, considering the four different FAL models. These atmospheric models represent a faint inter-network region of the quiet Sun (FAL-A), an average region of the quiet Sun (FAL-C), a bright network region of the quiet Sun (FAL-F), and a typical plage area (FAL-P).

Figure 1 shows the emergent Q/I profiles for the aforementioned atmospheric models and a LOS with μ = 0.1. In the absence of magnetic fields, the U/I profiles vanish and are consequently not shown. An interesting result is the following: the small trough in the line-core peak of the Q/I profile, obtained when the AA approximation is considered, completely disappears when a general AD calculation is carried out. The trough found in the AA case can thus be attributed to averaging over the scattering angle explained in Section 2.4. This clearly highlights the scientific value of AD calculations, as the sharper Q/I peak attained in this case is in much better agreement with the observations (e.g., Gandorfer 2002). By contrast, a good agreement between AA and AD calculations is found in the wings of the Q/I profile.

Fig. 1. Emergent Q/I profiles calculated in different atmospheric models for a LOS with μ = 0.1, taking into account PRD effects both under the AA approximation (solid lines) and in the general AD case (dashed lines). No magnetic field is considered. The considered atmospheric models are: FAL-C (red and orange lines), FAL-A (blue and light blue lines), FAL-F (green and pale green lines), and FAL-P (purple and violet lines). The reference direction for positive Q is considered to be parallel to the limb.

AA and AD calculations are now compared considering the FAL-C atmospheric model in the absence of magnetic fields for different LOSs. The results shown in Fig. 2 are consistent with those of Fig. 1. For all the considered LOSs, AA calculations yield a trough in the line-core peak of the Q/I profile, whereas the AD counterparts show a peak. Additionally, AA and AD calculations substantially agree in the wings of the Q/I profile. As expected, the amplitude of the Q/I profile, both in the core and in the wings, decreases when moving toward the disk center (i.e., for increasing μ values).

Fig. 2. Emergent Q/I profiles calculated in the FAL-C atmospheric model for different LOSs, taking into account PRD effects both under the AA approximation (solid lines) and in the general AD case (dashed lines). No magnetic field is considered. The considered LOS are: μ = 0.1 (red and orange lines), μ = 0.3 (blue and light blue lines), μ = 0.5 (green and pale green lines), and μ = 0.9 (purple and violet lines). The reference direction for positive Q is considered to be parallel to the limb.

3.4. Magnetic case

The linear scattering polarization signal of the Ca I 4227 Å line is primarily sensitive to the presence of magnetic fields through the combined action of the Hanle and magneto-optical effects. The Hanle effect operates in the line core, causing a reduction in the degree of linear polarization and a rotation in the plane of linear polarization. The Hanle critical field for this line is of approximately 25 G. Magneto-optical effects mainly operate in the line wings, producing a rotation in the plane of linear polarization, and a reduction in the linear polarization degree. The impact of this mechanism becomes appreciable in the presence of longitudinal magnetic fields with strengths similar to those required for the Hanle effect to have a substantial impact (more details in Alsina Ballester et al. 2018, 2019). The polarization of the Ca I 4227 Å line is also sensitive to the well-known Zeeman effect, which produces measurable circular polarization signals in the presence of magnetic fields with strengths larger than 10 G in the upper solar photosphere.

Figure 3 shows the impact of the AD treatment on the emergent Stokes profiles of the Ca I 4227 Å line, for a LOS with μ = 0.1, when including magnetic fields. The direction of the magnetic field is horizontal and contained in the plane defined by the local vertical and the LOS, while different field strengths have been considered. For the considered range of magnetic field strengths (between 0 and 50 G), the Hanle effect produces a decrease in the amplitude of the line-core peak of the Q/I profile (until the disappearance of the peak), while giving rise to a peak in the core of the U/I profile. A trough is found in both peaks when performing AA calculations, while the AD counterparts show a peak. The same figure allows one to perceive the sensitivity of the line wings to the same magnetic fields through magneto-optical effects: as the magnetic field increases, the wing lobes of Q/I decrease, while similar wing lobes appear in U/I. At such wing wavelengths, AA and AD calculations show a good agreement, thus suggesting that the AA approximation could be suitable to model and exploit magneto-optical effects in the scattering polarization wings of the Ca I 4227 Å line.

Fig. 3. Emergent Q/I (left panel) and U/I (right panel) calculated in the FAL-C atmospheric model for a LOS with μ = 0.1, taking into account PRD effects both under the AA approximation (solid lines) and in the general AD case (dashed lines). The magnetic field is horizontal and contained in the plane defined by the local vertical and the LOS. Different field strengths have been applied: 0 G (red and orange lines), 10 G (blue and light blue line), 20 G (green and pale green lines), and 50 G (purple and violet lines). The reference direction for positive Q is considered to be parallel to the limb.

4. Conclusions

This paper presents a quantitative comparison between the AA and AD approaches in the numerical solution of the NLTE RT problem of polarized radiation, taking scattering polarization and PRD effects into account. In particular, the Q/I and U/I Stokes profiles of the Ca I 4227 Å spectral line have been modeled and discussed, considering realistic 1D semi-empirical atmospheric models in the presence of magnetic fields.

The most relevant feature of the polarization signals obtained when performing AD calculations is the disappearance of the artificial trough in the line-core peak of the Q/I and U/I scattering polarization profiles obtained when using the AA approximation. This difference between AA and AD calculations is found for all LOS, both in the absence and in the presence of magnetic fields. This finding highlights the scientific value of AD calculations, as the corresponding line-core peak (which is a key observable for Hanle effect diagnostics) is in much better agreement with the observations. Previous investigations have shown that other strong resonance lines, such as Na I D₂, show a trough in the line core peak of the Q/I profile when modeled under the AA approximation (e.g., Belluzzi et al. 2015; Alsina Ballester et al. 2021). It is expected that also in these cases, the trough will disappear when AD calculations are carried out.

Additionally, the numerical results support the validity of the AA approximation for modeling the wings of the Q/I and U/I profiles in the static case. This validates the investigation on the Ca I 4227 Å line performed by Alsina Ballester et al. (2018) under the AA approximation, strengthening the diagnostic potential of magneto-optical effects for exploring the magnetism of stellar atmospheres.

The numerical approach presented in this paper can be applied to model other resonance lines, and it can be generalized to a wider range of scenarios, such as two-term atomic models, and 3D atmospheric models with arbitrary magnetic and bulk velocity fields. In particular, the impact of bulk velocity motions on the modeling of the scattering polarization of resonance lines taking PRD effects into account in the general AD case will be subject to an upcoming investigation.

Acknowledgments

The financial support by the Swiss National Science Foundation (SNSF) through grant CRSII5_180238 is gratefully acknowledged. E.A.B. and L.B. gratefully acknowledge financial support by SNSF through grant 200021_175997. J.T.B. acknowledges the funding received from the European Research Counsil (ERC) under the European Union’s Horizon 2020 research and innovation programme (ERC Advanced Grant agreement No 742265).

References

All Figures

Fig. 1. Emergent Q/I profiles calculated in different atmospheric models for a LOS with μ = 0.1, taking into account PRD effects both under the AA approximation (solid lines) and in the general AD case (dashed lines). No magnetic field is considered. The considered atmospheric models are: FAL-C (red and orange lines), FAL-A (blue and light blue lines), FAL-F (green and pale green lines), and FAL-P (purple and violet lines). The reference direction for positive Q is considered to be parallel to the limb.

In the text

Fig. 2. Emergent Q/I profiles calculated in the FAL-C atmospheric model for different LOSs, taking into account PRD effects both under the AA approximation (solid lines) and in the general AD case (dashed lines). No magnetic field is considered. The considered LOS are: μ = 0.1 (red and orange lines), μ = 0.3 (blue and light blue lines), μ = 0.5 (green and pale green lines), and μ = 0.9 (purple and violet lines). The reference direction for positive Q is considered to be parallel to the limb.

In the text

Fig. 3. Emergent Q/I (left panel) and U/I (right panel) calculated in the FAL-C atmospheric model for a LOS with μ = 0.1, taking into account PRD effects both under the AA approximation (solid lines) and in the general AD case (dashed lines). The magnetic field is horizontal and contained in the plane defined by the local vertical and the LOS. Different field strengths have been applied: 0 G (red and orange lines), 10 G (blue and light blue line), 20 G (green and pale green lines), and 50 G (purple and violet lines). The reference direction for positive Q is considered to be parallel to the limb.

In the text