© ESO, 2011

1. Introduction

This paper is concerned with a study of the linear polarization of spectral lines due to the resonance scattering of anisotropic radiation. We deal only with two-level atoms with an unpolarized ground-level. The theory of resonance scattering shows that in general correlations exist between the directions and frequencies of the incident and scattered photons, a situation referred to as “partial frequency redistribution” (PRD). In some cases, these correlations can be ignored and one has the so-called “complete frequency redistribution” (CRD). The line polarization is very sensitive to the nature of the frequency redistribution mechanism. Polarized radiative transfer problems with PRD are much harder to solve than those with CRD. The complexity is particularly significant because the scattering process is in general described by a (4 × 4) redistribution matrix. In the special case of the non-magnetic resonance (Rayleigh) scattering in planar atmospheres with axisymmetric boundary conditions, it is sufficient to use (2 × 2) redistribution matrices for the Stokes vector (I,Q).

We are particularly interested in radiative transfer with an angle-dependent PRD. Dumont et al. (1977) were the first to consider line polarization with angle-dependent PRD, using the Hummer (1962)’s type I redistribution function. McKenna (1985) used Hummer’s types I, II, and III, and Faurobert (1987, 1988) the type II¹. Nagendra et al. (2002, 2003) considered a more general problem of angle-dependent PRD for the weak-field Hanle effect using the PRD matrices derived by Bommier (1997b). An even more general redistribution matrix for the Hanle-Zeeman scattering was proposed in Sampoorna et al. (2007a,b) and used in Sampoorna et al. (2008a) to calculate linear and circular polarizations of spectral lines in magnetic fields of arbitrary strength.

The difficulty in using angle-dependent PRD matrices comes mainly from the evaluation of the scattering terms because the redistribution matrices depend in an intricate way on the directions and frequencies of the incident and scattered beams. For this reason, most of the linear polarization investigations taking into account PRD effects have been carried out with an approximate angle-averaged PRD, as suggested by Rees & Saliba (1982). In this approximation, frequency redistribution is independent of the scattering angle and also of the polarization phase matrix. This so-called “hybrid” approximation has been employed in works by Rees & Saliba (1982), McKenna (1985), Faurobert (1987, 1988), Nagendra (1988, 1994), Faurobert-Scholl (1991), Nagendra et al. (1999), Fluri et al. (2003), Sampoorna et al. (2008b), Sampoorna & Trujillo Bueno (2010), among others. For the Rayleigh scattering, the use of an angle-averaged PRD function is not unjustified as the quantitative differences between the Q/I profiles calculated with angle-dependent (AD) and angle-averaged (AA) frequency redistribution remain around 15% or less as shown by Faurobert (1987) and Nagendra et al. (2002). Improvements in observational techniques is a strong incentive to develop more elaborate diagnostic tools by taking into account the angle-dependent frequency redistribution.

Here we describe three new numerical methods to handle the Rayleigh scattering with an angle-dependent PRD. Two of them are of the accelerated lambda iteration (ALI) type and the other is based on a Neumann series expansion, which amounts to an expansion in the mean number of scattering events. They differ rather strongly from previous methods used for angle-dependent PRD. In Dumont et al. (1977), Faurobert (1987, 1988), Nagendra et al. (2002), Sampoorna et al. (2008a), the radiation field is described by the two Stokes parameters I and Q and the transfer equations for I and Q are solved by Feautrier-type methods, sometimes associated with a perturbation method, or fully perturbative methods, which is based on the linear polarization created by resonance scattering being a few percent only. In McKenna (1985), the radiation field is also described by I and Q, but the radiative transfer equations are solved by a moment equation method.

ALI methods introduced for scalar radiative transfer problems have been generalized to handle Rayleigh scattering and the weak-field Hanle effect (see the review by Nagendra 2003; and also Nagendra & Sampoorna 2009). These methods have the advantage of being much faster while remaining as accurate as the traditional exact or perturbative methods (see Nagendra et al. 1999). They have been so far applied to the angle-averaged PRD only. They make use of a decomposition of the Stokes parameters into a set of new fields, referred to as “reduced intensities” in Faurobert-Scholl (1991) and Nagendra et al. (1998) and as “spherical irreducible components” in Frisch (2007). These decompositions were introduced to study the Hanle effect. The decomposition method in Frisch (2007) is based on the decomposition of the polarization phase matrix in terms of the spherical irreducible tensors for polarimetry ${\mathrm{𝒯}}_{\mathit{Q}}^{\mathit{K}}\mathrm{\left(}\mathit{i,}\mathbf{\Omega }\mathrm{\right)}$ introduced by Landi Degl’Innocenti (1984, see also Bommier 1997a; Landi Degl’Innocenti & Landolfi 2004). In Faurobert-Scholl (1991, see also Nagendra et al. 1998), the decomposition relies on an azimuthal Fourier expansion method. It is a generalization of Chandrasekhar’s (1950) azimuthal Fourier expansion method for Rayleigh scattering, to the case of Hanle scattering. The ${\mathrm{𝒯}}_{\mathit{Q}}^{\mathit{K}}$ based expansion technique turns out to be simpler than the azimuthal Fourier expansion method.

The irreducible spherical components satisfy transfer equations that are simpler than the equations for I and Q. For angle-averaged PRD functions, the irreducible components of the source terms become independent of the ray direction, even in the presence of a magnetic field, and satisfy fairly standard integral equations that can be used to construct ALI numerical methods of solution. It has been shown in Frisch (2009, 2010) that a similar decomposition of the Stokes parameters can be performed with angle-dependent PRD functions. The Hanle effect is considered in Frisch (2009) and the Rayleigh scattering in Frisch (2010, hereafter HF10). Here we show how the decomposition described in HF10 can be used to construct ALI and a Neumann series expansion methods.

The outline of the paper is as follows. In Sect. 2, we present transfer equations for the irreducible components of the radiation field and the main steps of the decomposition technique. In Sect. 3, we present two different ALI methods, one of which generalizes the frequency-by-frequency (FBF) method of Paletou & Auer (1995) into a frequency-angle by frequency-angle (FABFA) method and the other that generalizes the core-wing separation method, also of Paletou & Auer (1995). In this section, we present a third iterative method, which relies on a Neumann series expansion of the irreducible components of the source terms contributing to the polarization. We refer to this approach here as the “scattering expansion method” because it is equivalent to an expansion in the mean number of scattering events. The three methods described in Sect. 3 make use of the azimuthal Fourier coefficients of order 0, 1, and 2 of the Hummer’s PRD functions of types II and III. In Sect. 4, we show how these Fourier coefficients can be calculated and discuss their main properties. Numerical validation and convergence properties of the iterative methods are presented in Sect. 5. Results are presented in Sect. 6, where we discuss in detail the angle-dependent PRD effects on the Q/I profiles. Some concluding remarks are presented in Sect. 7.

2. Governing equations

The atmosphere is assumed to be plane-parallel, the primary source of photons to be of thermal origin, the incident radiation to be either zero or axisymmetric, and the magnetic field to be zero or micro-turbulent. These assumptions imply that the radiation field is cylindrically symmetric and can be described by the two Stokes parameters I and Q, if one chooses the reference direction for the measurement of Q such that positive Q is perpendicular to the surface of the atmosphere. Because of the cylindrical symmetry Stokes U = 0.

The polarized transfer equation for the Stokes parameters I and Q can be written in a component form as $\mathit{\mu }\frac{\mathit{\partial }{\mathit{I}}_{\mathit{i}}}{\mathit{\partial \tau }}\mathrm{=}\left[\mathit{\varphi }\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{+}\mathit{r}\right]\left[{\mathit{I}}_{\mathit{i}}\mathrm{\left(}\mathit{\tau ,x,\mu }\mathrm{\right)}\mathrm{-}{\mathit{S}}_{\mathit{i}}\mathrm{\left(}\mathit{\tau ,x,\mu }\mathrm{\right)}\right]\mathit{,} \mathit{i}\mathrm{=}\mathrm{0}\mathit{,}\mathrm{1}\mathit{,}$(1)where μ = cosθ, with θ being the co-latitude with respect to the atmospheric normal, τ the line optical depth defined by dτ =  −k_ldz, with k_l the frequency-averaged line absorption coefficient, and ϕ(x) the normalized Voigt function. The frequency x is measured in units of the Doppler width, assumed to be constant, with x = 0 at line center. The ratio of continuum to line absorption coefficient is denoted by r. The total source vector is given by ${\mathit{S}}_{\mathit{i}}\mathrm{\left(}\mathit{\tau ,x,\mu }\mathrm{\right)}\mathrm{=}\frac{\mathit{\varphi }\mathrm{\left(}\mathit{x}\mathrm{\right)}{\mathit{S}}_{\mathrm{l}\mathit{,i}}\mathrm{\left(}\mathit{\tau ,x,\mu }\mathrm{\right)}\mathrm{+}\mathit{r}{\mathit{S}}_{\mathrm{c}\mathit{,i}}}{\mathit{\varphi }\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{+}\mathit{r}}\mathit{,}$(2)where S_c,i are the components of the unpolarized continuum source vector. We assume that S_c,0 = B, where B is the Planck function at line center, and S_c,1 = 0. The line source vector can be written as $\begin{array}{ccc}& & {\mathit{S}}_{\mathrm{l}\mathit{,i}}\mathrm{\left(}\mathit{\tau ,x,\mu }\mathrm{\right)}\mathrm{=}{\mathit{G}}_{\mathit{i}}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\\ & & \mathrm{+}\mathrm{\int }\mathrm{􏽉}\sum _{\mathit{j}\mathrm{=}\mathrm{0}\mathit{,}\mathrm{1}}\frac{{\mathit{R}}_{\mathit{ij}}\mathrm{\left(}\mathit{x,}\mathbf{\Omega }\mathit{,}{\mathit{x}}^{\mathrm{\prime }}\mathit{,}{\mathbf{\Omega }}^{\mathrm{\prime }}\mathrm{\right)}}{\mathit{\varphi }\mathrm{\left(}\mathit{x}\mathrm{\right)}}{\mathit{I}}_{\mathit{j}}\mathrm{\left(}\mathit{\tau ,}{\mathit{x}}^{\mathrm{\prime }}\mathit{,}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{\right)}\frac{\hspace{0.17em}\mathrm{d}{\mathrm{\Omega }}^{\mathrm{\prime }}}{\mathrm{4}\mathit{\pi }}\hspace{0.17em}\mathrm{d}{\mathit{x}}^{\mathrm{\prime }}\mathit{,}\end{array}$(3)where dΩ′ = sinθ′   dθ′   dχ′. The outgoing and incoming ray directions Ω and Ω′ are defined, respectively, by their polar angles (θ,χ) and (θ′,χ′). For simplicity, we assume that the primary source is unpolarized, namely that only G₀(τ) is non-zero. It is proportional to the Planck function at line center. The term R_ij(x,Ω,x′,Ω′) denote the elements of the redistribution matrix for Rayleigh scattering (Domke & Hubeny 1988; Bommier 1997a).

According to the decomposition technique described in HF10, we can write I_i and S_l,i as ${\mathit{I}}_{\mathit{i}}\mathrm{\left(}\mathit{\tau ,x,\mu }\mathrm{\right)}\mathrm{=}\sum _{\mathit{K,Q}\mathrm{\ge }\mathrm{0}}\stackrel{}{{{\mathrm{𝒯}}_{\mathit{Q}}^{\mathit{K}}}_{\mathrm{˜}}}\mathrm{\left(}\mathit{i,\mu }\mathrm{\right)}{\mathrm{ℐ}}_{\mathit{Q}}^{\mathit{K}}\mathrm{\left(}\mathit{\tau ,x,\mu }\mathrm{\right)}\mathit{,} \mathit{i}\mathrm{=}\mathrm{0}\mathit{,}\mathrm{1}\mathit{,}$(4)${\mathit{S}}_{\mathrm{l}\mathit{,i}}\mathrm{\left(}\mathit{\tau ,x,\mu }\mathrm{\right)}\mathrm{=}\sum _{\mathit{K,Q}\mathrm{\ge }\mathrm{0}}\stackrel{}{{{\mathrm{𝒯}}_{\mathit{Q}}^{\mathit{K}}}_{\mathrm{˜}}}\mathrm{\left(}\mathit{i,\mu }\mathrm{\right)}{\mathrm{𝒮}}_{\mathrm{l}\mathit{,Q}}^{\mathit{K}}\mathrm{\left(}\mathit{\tau ,x,\mu }\mathrm{\right)}\mathit{,} \mathit{i}\mathrm{=}\mathrm{0}\mathit{,}\mathrm{1.}$(5)We have four terms in the summation over K and Q corresponding to K = Q = 0, K = 2 with Q = 0,1,2. The irreducible tensors $\begin{array}{}{}_{\mathrm{˜}}\\ {\mathrm{𝒯}}_{\mathit{Q}}^{\mathit{K}}\end{array}\mathrm{\left(}\mathit{i,\mu }\mathrm{\right)}$ are defined in HF10. An explicit expression of Eq. (4) can be found in Eq. (42). The total source vector S_i has a decomposition similar to Eq. (5). The irreducible line source vector components ${\mathrm{𝒮}}_{\mathrm{l}\mathit{,Q}}^{\mathit{K}}\mathrm{\left(}\mathit{\tau ,x,\mu }\mathrm{\right)}$ may be written as $\begin{array}{ccc}& & {\mathrm{𝒮}}_{\mathrm{l}\mathit{,Q}}^{\mathit{K}}\mathrm{\left(}\mathit{\tau ,x,\mu }\mathrm{\right)}\mathrm{=}{\mathit{\delta }}_{\mathit{K}\mathrm{0}}{\mathit{\delta }}_{\mathit{Q}\mathrm{0}}{\mathit{G}}_{\mathrm{0}}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\mathrm{+}{\mathrm{\int }}_{\mathrm{-}\mathrm{\infty }}^{\mathrm{+}\mathrm{\infty }}\hspace{0.17em}\mathrm{d}{\mathit{x}}^{\mathrm{\prime }}{\mathrm{\int }}_{-1}^{\mathrm{+}\mathrm{1}}\frac{\hspace{0.17em}\mathrm{d}{\mathit{\mu }}^{\mathrm{\prime }}}{\mathrm{2}}\\ & & \mathrm{\times }\frac{\begin{array}{}\mathrm{˜}\\ {\mathrm{ℛ}}_{\mathit{Q}}^{\mathit{K}}\end{array}\mathrm{\left(}\mathit{x,\mu ,}{\mathit{x}}^{\mathrm{\prime }}\mathit{,}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{\right)}}{\mathit{\varphi }\mathrm{\left(}\mathit{x}\mathrm{\right)}}\sum _{{\mathit{K}}^{\mathrm{\prime }}\mathit{,}{\mathit{Q}}^{\mathrm{\prime }}\mathrm{\ge }\mathrm{0}}\stackrel{}{{{\mathrm{\Gamma }}_{\mathit{Q}{\mathit{Q}}^{\mathrm{\prime }}}^{\mathit{K}{\mathit{K}}^{\mathrm{\prime }}}}_{\mathrm{˜}}}\mathrm{\left(}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{\right)}{\mathrm{ℐ}}_{{\mathit{Q}}^{\mathrm{\prime }}}^{{\mathit{K}}^{\mathrm{\prime }}}\mathrm{\left(}\mathit{\tau ,}{\mathit{x}}^{\mathrm{\prime }}\mathit{,}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{\right)}\mathit{,}\end{array}$(6)where $\stackrel{}{{{\mathrm{\Gamma }}_{\mathit{Q}{\mathit{Q}}^{\mathrm{\prime }}}^{\mathit{K}{\mathit{K}}^{\mathrm{\prime }}}}_{\mathrm{˜}}}\mathrm{\left(}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{=}\sum _{\mathit{j}\mathrm{=}\mathrm{0}\mathit{,}\mathrm{1}}\stackrel{}{{{\mathrm{𝒯}}_{\mathit{Q}}^{\mathit{K}}}_{\mathrm{˜}}}\mathrm{\left(}\mathit{j,}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{\right)}\stackrel{}{{{\mathrm{𝒯}}_{{\mathit{Q}}^{\mathrm{\prime }}}^{{\mathit{K}}^{\mathrm{\prime }}}}_{\mathrm{˜}}}\mathrm{\left(}\mathit{j,}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{\right)}\mathit{.}$(7)The coefficients $\stackrel{}{{{\mathrm{\Gamma }}_{\mathit{Q}{\mathit{Q}}^{\mathrm{\prime }}}^{\mathit{K}{\mathit{K}}^{\mathrm{\prime }}}}_{\mathrm{˜}}}\mathrm{\left(}\mathit{\mu }\mathrm{\right)}$ are given in the Appendix of HF10. The functions $\begin{array}{}{\mathrm{˜}}_{}\\ {\mathrm{ℛ}}_{\mathit{Q}}^{\mathit{K}}\end{array}$ in Eq. (6) take the form ${\begin{array}{}\mathrm{˜}\\ \mathrm{ℛ}\end{array}}_{\mathrm{0}}^{\mathrm{0}}\mathrm{=}\mathit{\alpha }\stackrel{}{{{\mathit{r}}_{\mathrm{II}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}}_{\mathrm{˜}}}\mathrm{+}\mathrm{[}{\mathit{\beta }}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{-}{\mathit{\alpha }}^{\mathrm{]}}\stackrel{}{{{\mathit{r}}_{\mathrm{III}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}}_{\mathrm{˜}}}\mathit{,}$(8)${\begin{array}{}\mathrm{˜}\\ \mathrm{ℛ}\end{array}}_{\mathit{Q}}^{\mathrm{2}}\mathrm{=}{\mathit{W}}_{\mathrm{2}}{\mathit{\mu }}_{\mathrm{2}}\left\{\mathit{\alpha }\stackrel{}{{{\mathit{r}}_{\mathrm{II}}^{\mathrm{\left(}\mathit{Q}\mathrm{\right)}}}_{\mathrm{˜}}}\mathrm{+}\mathrm{[}{\mathit{\beta }}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}\mathrm{-}{\mathit{\alpha }}^{\mathrm{]}}\stackrel{}{{{\mathit{r}}_{\mathrm{III}}^{\mathrm{\left(}\mathit{Q}\mathrm{\right)}}}_{\mathrm{˜}}}\right\}\mathit{,} \mathit{Q}\mathrm{=}\mathrm{0}\mathit{,}\mathrm{1}\mathit{,}\mathrm{2}\mathit{,}$(9)where W₂(J_l,J_u) is an atomic depolarization factor depending on the angular momentum of the lower and upper levels of the transition. The coefficients α and β^(K) (Bommier 1997b) are branching ratios $\mathit{\alpha }\mathrm{=}\frac{{\mathrm{\Gamma }}_{\mathrm{R}}}{{\mathrm{\Gamma }}_{\mathrm{R}}\mathrm{+}{\mathrm{\Gamma }}_{\mathrm{I}}\mathrm{+}{\mathrm{\Gamma }}_{\mathrm{E}}}\mathit{,}$(10)and ${\mathit{\beta }}^{\mathrm{\left(}\mathit{K}\mathrm{\right)}}\mathrm{=}\frac{{\mathrm{\Gamma }}_{\mathrm{R}}}{{\mathrm{\Gamma }}_{\mathrm{R}}\mathrm{+}{\mathrm{\Gamma }}_{\mathrm{I}}\mathrm{+}{\mathit{D}}^{\mathrm{\left(}\mathit{K}\mathrm{\right)}}}\mathit{,}$(11)where Γ_R is the radiative rate, Γ_I and Γ_E the inelastic and elastic collisional rates, and D^(K) the collisional depolarization rate such that D⁽⁰⁾ = 0. The coefficient μ₂ takes the effects of a micro-turbulent magnetic field into account. It depends on the magnetic field probability density function (see e.g., Landi Degl’Innocenti & Landolfi 2004, p. 215) and is unity in the absence of magnetic fields. Since the Hanle effect acts only in the line core, the coefficient μ₂ should be set to unity outside the line core. The $\stackrel{}{{{\mathit{r}}_{\mathrm{X}}^{\mathrm{\left(}\mathit{Q}\mathrm{\right)}}}_{\mathrm{˜}}}$ (with X = II  or  III) are the azimuthal Fourier coefficients of order 0, 1, and 2 of Hummer (1962)’s PRD functions r_II and r_III. They are defined by $\begin{array}{ccc}& & \stackrel{}{{{\mathit{r}}_{\mathrm{X}}^{\mathrm{\left(}\mathit{Q}\mathrm{\right)}}}_{\mathrm{˜}}}\mathrm{\left(}\mathit{x,\mu ,}{\mathit{x}}^{\mathrm{\prime }}\mathit{,}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{=}\frac{\mathrm{2}\mathrm{-}{\mathit{\delta }}_{\mathrm{0}\mathit{Q}}}{\mathrm{2}\mathit{\pi }}\\ & & \mathrm{\times }{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{2}\mathit{\pi }}{\mathit{r}}_{\mathrm{X}}\mathrm{(}\mathit{x,\mu ,}{\mathit{x}}^{\mathrm{\prime }}\mathit{,}{\mathit{\mu }}^{\mathrm{\prime }}\mathit{,\chi }\mathrm{-}{\mathit{\chi }}^{\mathrm{\prime }}\mathrm{)}\cos \mathrm{\left[}\mathit{Q}\mathrm{\right(}\mathit{\chi }\mathrm{-}{\mathit{\chi }}^{\mathrm{\prime }}\mathrm{\left)}\mathrm{\right]}\hspace{0.17em}\mathrm{d}\mathrm{(}\mathit{\chi }\mathrm{-}{\mathit{\chi }}^{\mathrm{\prime }}\mathrm{)}\mathit{.}\end{array}$(12)The components ${\mathrm{ℐ}}_{\mathit{Q}}^{\mathit{K}}$ satisfy a transfer equation similar to Eq. (1). The source term is given by Eq. (2), where S_l,i and S_c,i are replaced by ${\mathrm{𝒮}}_{\mathrm{l}\mathit{,Q}}^{\mathit{K}}$ and ${\mathrm{𝒮}}_{\mathrm{c}\mathit{,Q}}^{\mathit{K}}\mathrm{=}{\mathit{\delta }}_{\mathit{K}\mathrm{0}}{\mathit{\delta }}_{\mathit{Q}\mathrm{0}}\mathit{B}$. Introducing the four-component vectors ${\mathbf{𝒮}}_{\mathrm{l}}\mathrm{\left(}\mathit{\tau ,x,\mu }\mathrm{\right)}\mathrm{=}\mathrm{\{}{\mathrm{𝒮}}_{\mathrm{0}}^{\mathrm{0}}\mathit{,}{\mathrm{𝒮}}_{\mathrm{0}}^{\mathrm{2}}\mathit{,}{\mathrm{𝒮}}_{\mathrm{1}}^{\mathrm{2}}\mathit{,}{\mathrm{𝒮}}_{\mathrm{2}}^{\mathrm{2}}{\mathrm{\}}}^{\mathrm{T}}$ and $\mathbf{ℐ}\mathrm{\left(}\mathit{\tau ,x,\mu }\mathrm{\right)}\mathrm{=}\mathrm{\{}{\mathrm{ℐ}}_{\mathrm{0}}^{\mathrm{0}}\mathit{,}{\mathrm{ℐ}}_{\mathrm{0}}^{\mathrm{2}}\mathit{,}{\mathrm{ℐ}}_{\mathrm{1}}^{\mathrm{2}}\mathit{,}{\mathrm{ℐ}}_{\mathrm{2}}^{\mathrm{2}}{\mathrm{\}}}^{\mathrm{T}}$, we can re-write Eq. (6) in vector form as $\begin{array}{ccc}& & {\mathbf{𝒮}}_{\mathrm{l}}\mathrm{\left(}\mathit{\tau ,x,\mu }\mathrm{\right)}\mathrm{=}\mathbf{𝒢}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\\ & & \mathrm{+}{\mathrm{\int }}_{\mathrm{-}\mathrm{\infty }}^{\mathrm{+}\mathrm{\infty }}{\mathrm{\int }}_{-1}^{\mathrm{+}\mathrm{1}}\frac{\begin{array}{}\mathrm{˜}\\ \mathbf{ℛ}\end{array}\mathrm{\left(}\mathit{x,\mu ,}{\mathit{x}}^{\mathrm{\prime }}\mathit{,}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{\right)}}{\mathit{\varphi }\mathrm{\left(}\mathit{x}\mathrm{\right)}}\mathbf{\Gamma }\mathrm{\left(}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{\right)}\mathbf{ℐ}\mathrm{\left(}\mathit{\tau ,}{\mathit{x}}^{\mathrm{\prime }}\mathit{,}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{\right)}\frac{\hspace{0.17em}\mathrm{d}{\mathit{\mu }}^{\mathrm{\prime }}}{\mathrm{2}}\hspace{0.17em}\mathrm{d}{\mathit{x}}^{\mathrm{\prime }}\mathit{.}\end{array}$(13)The primary source vector is , where G₀(τ) = ϵB with ϵ = Γ_I/(Γ_I + Γ_R). The (4 × 4) matrix $\begin{array}{}{\mathrm{˜}}_{}\\ \mathbf{ℛ}\end{array}$ is diagonal, i.e., $\begin{array}{}{}_{\mathrm{˜}}\\ \mathbf{ℛ}\end{array}\mathrm{=}\mathrm{diag}{}^{\mathrm{[}}\begin{array}{}{}_{\mathrm{˜}}\\ {\mathrm{ℛ}}_{\mathrm{0}}^{\mathrm{0}}\end{array}\mathit{,}\begin{array}{}{}_{\mathrm{˜}}\\ {\mathrm{ℛ}}_{\mathrm{0}}^{\mathrm{2}}\end{array}\mathit{,}\begin{array}{}{}_{\mathrm{˜}}\\ {\mathrm{ℛ}}_{\mathrm{1}}^{\mathrm{2}}\end{array}\mathit{,}{\begin{array}{}{}_{\mathrm{˜}}\\ {\mathrm{ℛ}}_{\mathrm{2}}^{\mathrm{2}}\end{array}}^{\mathrm{]}}$. Its elements are defined in Eqs. (8) and (9). The (4 × 4) matrix Γ is a full matrix with elements $\stackrel{}{{{\mathrm{\Gamma }}_{\mathit{Q}{\mathit{Q}}^{\mathrm{\prime }}}^{\mathit{K}{\mathit{K}}^{\mathrm{\prime }}}}_{\mathrm{˜}}}$. Owing to its symmetry, it has only ten independent elements (see Eq. (7) and also HF10). When $\begin{array}{}{\mathrm{˜}}_{}\\ \mathbf{ℛ}\end{array}$ is independent of μ and μ′ (CRD or angle-averaged PRD), only the two components of the source vector corresponding to the index Q = 0 are non-zero. In this case, one recovers the usual Rayleigh scattering whereby the radiation field and source vector are fully described by two-component vectors.

3. Numerical methods of solution

We present three iterative methods to solve the problem of angle-dependent PRD. The first two are ALI type methods and the third is a scattering expansion method. We assume a plane-parallel slab geometry with a given total optical thickness.

3.1. The polarized accelerated lambda iteration approaches

For notational simplicity, we neglect the explicit dependence of ℐ and on τ. The dependence on x and μ appear as subscripts. The formal solution of the transfer equation for the four-component irreducible vector ℐ can be written as ${\mathbf{ℐ}}_{\mathit{x\mu }}\mathrm{=}{\mathbf{\Lambda }}_{\mathit{x\mu }}\mathrm{\left[}{\mathbf{𝒮}}_{\mathit{x\mu }}\mathrm{\right]}\mathrm{+}{{T}}_{\mathit{x\mu }}\mathit{,}$(14)where T_xμ is the directly transmitted part of the intensity vector and Λ_xμ is the frequency- and angle-dependent (4 × 4) integral operator. We introduce the operator splitting written as ${\mathbf{\Lambda }}_{\mathit{x\mu }}\mathrm{=}{\mathbf{\Lambda }}_{\mathit{x\mu }}^{\mathrm{\ast }}\mathrm{+}\mathrm{(}{\mathbf{\Lambda }}_{\mathit{x\mu }}\mathrm{-}{{\mathbf{\Lambda }}_{\mathit{x\mu }}^{\mathrm{\ast }}}^{\mathrm{)}}\mathit{,}$(15)where ${\mathbf{\Lambda }}_{\mathit{x\mu }}^{\mathrm{\ast }}$ is the diagonal approximate operator (see Olson et al. 1986). We now write the total source vector and the line source vector as $\begin{array}{ccc}& & {\mathbf{𝒮}}_{\mathit{x\mu }}^{\mathit{n}\mathrm{+}\mathrm{1}}\mathrm{=}{\mathbf{𝒮}}_{\mathit{x\mu }}^{\mathit{n}}\mathrm{+}\mathit{\delta }{\mathbf{𝒮}}_{\mathit{x\mu }}^{\mathit{n}}\mathit{,}\\ & & {\mathbf{𝒮}}_{\mathrm{l}\mathit{,x\mu }}^{\mathit{n}\mathrm{+}\mathrm{1}}\mathrm{=}{\mathbf{𝒮}}_{\mathrm{l}\mathit{,x\mu }}^{\mathit{n}}\mathrm{+}\mathit{\delta }{\mathbf{𝒮}}_{\mathrm{l}\mathit{,x\mu }}^{\mathit{n}}\mathit{,}\end{array}$(16)where n is the iteration index. Combining Eqs. (15) and (16) with Eq. (13), we derive an equation for the line source vector corrections that can be written as $\begin{array}{ccc}& & \mathit{\delta }{\mathbf{𝒮}}_{\mathrm{l}\mathit{,x\mu }}^{\mathit{n}}\mathrm{-}{\mathrm{\int }}_{\mathrm{-}\mathrm{\infty }}^{\mathrm{+}\mathrm{\infty }}{\mathrm{\int }}_{-1}^{\mathrm{+}\mathrm{1}}\frac{{\begin{array}{}\mathrm{˜}\\ \mathbf{ℛ}\end{array}}_{\mathit{x\mu ,}{\mathit{x}}^{\mathrm{\prime }}{\mathit{\mu }}^{\mathrm{\prime }}}}{{\mathit{\varphi }}_{\mathit{x}}}{\mathbf{\Gamma }}_{{\mathit{\mu }}^{\mathrm{\prime }}}{\mathit{p}}_{{\mathit{x}}^{\mathrm{\prime }}}\\ & & \mathrm{\times }{\mathbf{\Lambda }}_{{\mathit{x}}^{\mathrm{\prime }}{\mathit{\mu }}^{\mathrm{\prime }}}^{\mathrm{\ast }}\mathrm{[}\mathit{\delta }{{\mathbf{𝒮}}_{\mathrm{l}\mathit{,}{\mathit{x}}^{\mathrm{\prime }}{\mathit{\mu }}^{\mathrm{\prime }}}^{\mathit{n}}}^{\mathrm{]}}\frac{\hspace{0.17em}\mathrm{d}{\mathit{\mu }}^{\mathrm{\prime }}}{\mathrm{2}}\hspace{0.17em}\mathrm{d}{\mathit{x}}^{\mathrm{\prime }}\mathrm{=}\mathbf{𝒢}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\mathrm{+}\overline{\mathbf{𝒥}}\begin{array}{}\hspace{0.17em}\mathit{n}\\ \mathit{x\mu }\end{array}\mathrm{-}{\mathbf{𝒮}}_{\mathrm{l}\mathit{,x\mu }}^{\mathit{n}}\mathit{,}\end{array}$(17)where p_x = ϕ_x/(ϕ_x + r) and $\overline{\mathbf{𝒥}}\begin{array}{}\hspace{0.17em}\mathit{n}\\ \mathit{x\mu }\end{array}\mathrm{=}{\mathrm{\int }}_{\mathrm{-}\mathrm{\infty }}^{\mathrm{+}\mathrm{\infty }}{\mathrm{\int }}_{-1}^{\mathrm{+}\mathrm{1}}\frac{{\begin{array}{}\mathrm{˜}\\ \mathbf{ℛ}\end{array}}_{\mathit{x\mu ,}{\mathit{x}}^{\mathrm{\prime }}{\mathit{\mu }}^{\mathrm{\prime }}}}{{\mathit{\varphi }}_{\mathit{x}}}{\mathbf{\Gamma }}_{{\mathit{\mu }}^{\mathrm{\prime }}}{\mathbf{\Lambda }}_{{\mathit{x}}^{\mathrm{\prime }}{\mathit{\mu }}^{\mathrm{\prime }}}\mathrm{[}{{\mathbf{𝒮}}_{{\mathit{x}}^{\mathrm{\prime }}{\mathit{\mu }}^{\mathrm{\prime }}}^{\mathit{n}}}^{\mathrm{]}}\frac{\hspace{0.17em}\mathrm{d}{\mathit{\mu }}^{\mathrm{\prime }}}{\mathrm{2}}\hspace{0.17em}\mathrm{d}{\mathit{x}}^{\mathrm{\prime }}\mathit{.}$(18)For this investigation, we chose the Jacobi decomposition of the Λ operator. Superior iterative methods such as the Gauss-Seidel (GS) and the successive overrelaxation (SOR) technique were introduced into scalar radiative transfer theory by Trujillo Bueno & Fabiani Bendicho (1995). These methods provide faster convergence rates than the Jacobi method. Trujillo Bueno & Manso Sainz (1999) extended these methods to the case of resonance scattering. The generalization to the Hanle effect was developed by Manso Sainz & Trujillo Bueno (1999). In these references, CRD is assumed for the frequency redistribution. The GS and SOR methods were extended to the case of angle-averaged PRD in Sampoorna & Trujillo Bueno (2010). In this first attempt at developing ALI methods for angle-dependent PRD functions, we chose the simpler Jacobi decomposition. We now describe two different methods of solution for Eq. (17).

3.1.1. Frequency-angle by frequency-angle (FABFA) method

This FABFA method is a generalization of the FBF method introduced by Paletou & Auer (1995). In matrix form, Eq. (17) can be written as $\mathbf{A}\mathit{\delta }{\mathbf{𝒮}}_{\mathrm{l}}^{\mathit{n}}\mathrm{=}{{r}}^{\mathit{n}}\mathit{,}$(19)where the residual vector rⁿ is given by the right-hand side of Eq. (17). At each depth point, rⁿ and $\mathit{\delta }{\mathbf{𝒮}}_{\mathrm{l}}^{\mathit{n}}$ are vectors of length 4N_x   2N_μ, where N_x is the number of frequency points in the range  [0,x_max] ² and N_μ is the number of angle points in the range  [0 < μ ≤ 1] . The matrix A thus has dimensions (4N_x   2N_μ × 4N_x   2N_μ). For a given x, x′, μ, and μ′, the matrix A can be decomposed into (N_x   2N_μ × N_x   2N_μ) blocks of 4 × 4 elements. In each block, denoted by , the elements may be written as $\mathbf{𝒜}\mathrm{=}{\mathit{\delta }}_{\mathit{mn}}{\mathit{\delta }}_{\mathit{\alpha \beta }}\mathbf{E}\mathrm{-}{\mathit{w}}_{\mathit{\beta }}{\mathbf{g}}_{\mathit{m\alpha ,n\beta }}{\mathbf{\Gamma }}_{\mathit{\beta }}{\mathbf{\Lambda }}_{\mathit{n\beta }}^{\mathrm{\ast }}\mathit{,}$(20)where m = 1,...,N_x, n = 1,...,N_x′, α = 1,...,2N_μ, and β = 1,...,2N_μ′, and E is the identity operator. The coefficients w_β denote the μ′ integration weights and g_mα,nβ are defined by ${\mathbf{g}}_{\mathit{m\alpha ,n\beta }}\mathrm{=}\frac{{\begin{array}{}\mathrm{˜}\\ \mathbf{ℛ}\end{array}}_{\mathit{m\alpha ,n\beta }}}{{\mathit{\varphi }}_{\mathit{m}}}\mathit{w̅}\mathit{n}\mathit{,}$(21)where are the frequency (x′) integration weights. The FABFA method requires the calculation of the matrix A^-1 before the iteration cycle.

3.1.2. The core-wing separation method for angle-dependent PRD

Here we describe a generalization of the core-wing method of Paletou & Auer (1995). According to the ALI methods, the right-hand side of Eq. (17) must be calculated as accurately as possible but in the left-hand side there is some flexibility in the choice of the operator acting on the corrections $\mathit{\delta }{\mathbf{𝒮}}_{\mathrm{l}\mathit{,}{\mathit{x}}^{\mathrm{\prime }}{\mathit{\mu }}^{\mathrm{\prime }}}^{\mathit{n}}$. The choice of this operator affects the speed of convergence of the process but not the final solution. For Rayleigh scattering, we already know that the angle-dependence of the PRD functions has only a mild effect on the Q/I profiles. Hence, in the matrix ${\begin{array}{}{\mathrm{˜}}_{}\\ \mathbf{ℛ}\end{array}}_{\mathit{x\mu ,}{\mathit{x}}^{\mathrm{\prime }}{\mathit{\mu }}^{\mathrm{\prime }}}$ in the left-hand side, we retain the two terms corresponding to the index Q = 0 and set to zero the two other terms corresponding to Q = 1 and Q = 2. This is equivalent to making the approximation ${\begin{array}{}\mathrm{˜}\\ \mathbf{ℛ}\end{array}}_{\mathit{x\mu ,}{\mathit{x}}^{\mathrm{\prime }}{\mathit{\mu }}^{\mathrm{\prime }}}\mathrm{\simeq }{\mathbf{ℛ}}_{\mathit{x}{\mathit{x}}^{\mathrm{\prime }}}^{\mathrm{AA}}\mathrm{ℰ}\mathit{,}$(22)where ℰ and ${\mathbf{ℛ}}_{\mathit{x}{\mathit{x}}^{\mathrm{\prime }}}^{\mathrm{AA}}$ are (4 × 4) diagonal matrices defined by ℰ = diag [1,1,0,0]  and ${\mathbf{ℛ}}_{\mathit{x}{\mathit{x}}^{\mathrm{\prime }}}^{\mathrm{AA}}\mathrm{=}\mathrm{diag}\left[{{}^{\mathrm{\right(}}{\mathrm{ℛ}}_{\mathit{x}{\mathit{x}}^{\mathrm{\prime }}}^{\mathrm{AA}}^{\mathrm{\left)}}}_{\mathrm{0}}^{\mathrm{0}}\mathit{,}{{}^{\mathrm{\right(}}{\mathrm{ℛ}}_{\mathit{x}{\mathit{x}}^{\mathrm{\prime }}}^{\mathrm{AA}}^{\mathrm{\left)}}}_{\mathrm{0}}^{\mathrm{2}}\mathit{,}\mathrm{0}\mathit{,}\mathrm{0}\right]$, with ${\mathrm{(}{{\mathrm{ℛ}}_{\mathit{x}{\mathit{x}}^{\mathrm{\prime }}}^{\mathrm{AA}}}^{\mathrm{)}}}_{\mathrm{0}}^{\mathrm{0}}\mathrm{=}\mathit{\alpha }{\mathit{r}}_{\mathit{x}{\mathit{x}}^{\mathrm{\prime }}}^{\mathrm{II}\mathrm{-}\mathrm{AA}}\mathrm{+}\mathrm{[}{\mathit{\beta }}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{-}{\mathit{\alpha }}^{\mathrm{]}}{\mathit{r}}_{\mathit{x}{\mathit{x}}^{\mathrm{\prime }}}^{\mathrm{III}\mathrm{-}\mathrm{AA}}\mathit{,}$(23)${\mathrm{(}{{\mathrm{ℛ}}_{\mathit{x}{\mathit{x}}^{\mathrm{\prime }}}^{\mathrm{AA}}}^{\mathrm{)}}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{=}{\mathit{W}}_{\mathrm{2}}{\mathit{\mu }}_{\mathrm{2}}\left\{\mathit{\alpha }{\mathit{r}}_{\mathit{x}{\mathit{x}}^{\mathrm{\prime }}}^{\mathrm{II}\mathrm{-}\mathrm{AA}}\mathrm{+}\mathrm{[}{\mathit{\beta }}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}\mathrm{-}{\mathit{\alpha }}^{\mathrm{]}}{\mathit{r}}_{\mathit{x}{\mathit{x}}^{\mathrm{\prime }}}^{\mathrm{III}\mathrm{-}\mathrm{AA}}\right\}\mathit{,}$(24)and ${\mathit{r}}_{\mathit{x}{\mathit{x}}^{\mathrm{\prime }}}^{\mathrm{II}\mathrm{-}\mathrm{AA}}$ and ${\mathit{r}}_{\mathit{x}{\mathit{x}}^{\mathrm{\prime }}}^{\mathrm{III}\mathrm{-}\mathrm{AA}}$ are the type-II and type-III angle-averaged (AA) redistribution functions. By comparing with the FABFA method, we find that this approximation leads to a correct converged solution, but at a much less computational cost (see Sect. 5).

Substituting Eq. (22) in Eq. (17), we obtain $\begin{array}{ccc}& & \mathit{\delta }{\mathbf{𝒮}}_{\mathrm{l}\mathit{,x\mu }}^{\mathit{n}}\mathrm{-}\\ & & {\mathrm{\int }}_{\mathrm{-}\mathrm{\infty }}^{\mathrm{+}\mathrm{\infty }}\frac{{\mathbf{ℛ}}_{\mathit{x}{\mathit{x}}^{\mathrm{\prime }}}^{\mathrm{AA}}}{{\mathit{\varphi }}_{\mathit{x}}}{\mathrm{\int }}_{-1}^{\mathrm{+}\mathrm{1}}\mathrm{ℰ}{\mathbf{\Gamma }}_{{\mathit{\mu }}^{\mathrm{\prime }}}{\mathit{p}}_{{\mathit{x}}^{\mathrm{\prime }}}{\mathbf{\Lambda }}_{{\mathit{x}}^{\mathrm{\prime }}{\mathit{\mu }}^{\mathrm{\prime }}}^{\mathrm{\ast }}\mathrm{[}\mathit{\delta }{{\mathbf{𝒮}}_{\mathrm{l}\mathit{,}{\mathit{x}}^{\mathrm{\prime }}{\mathit{\mu }}^{\mathrm{\prime }}}^{\mathit{n}}}^{\mathrm{]}}\frac{\hspace{0.17em}\mathrm{d}{\mathit{\mu }}^{\mathrm{\prime }}}{\mathrm{2}}\hspace{0.17em}\mathrm{d}{\mathit{x}}^{\mathrm{\prime }}\mathrm{=}{{r}}_{\mathit{x\mu }}^{\mathit{n}}\mathit{.}\end{array}$(25)To the above equation, we now apply the core-wing method (Paletou & Auer 1995; Nagendra et al. 1999; Fluri et al. 2003), namely $\begin{array}{ccc}\frac{{\mathit{r}}_{\mathit{x}{\mathit{x}}^{\mathrm{\prime }}}^{\mathrm{II}\mathrm{-}\mathrm{AA}}}{{\mathit{\varphi }}_{\mathit{x}}}\mathrm{\simeq }\{\begin{array}{c}\\ & \\ & \end{array}& & \end{array}$(26)$\begin{array}{ccc}\frac{{\mathit{r}}_{\mathit{x}{\mathit{x}}^{\mathrm{\prime }}}^{\mathrm{III}\mathrm{-}\mathrm{AA}}}{{\mathit{\varphi }}_{\mathit{x}}}\mathrm{\simeq }\{\begin{array}{c}\\ & \\ & \end{array}& & \end{array}$(27)where x_c is the frequency that distinguishes the line core and the wing. A value of x_c = 3.5 Doppler widths for the separation frequency is a reasonable choice. Combining Eqs. (26) and (27) with Eq. (25), we obtain the expression for the line source vector corrections $\mathit{\delta }{\mathbf{𝒮}}_{\mathrm{l}\mathit{,x\mu }}^{\mathit{n}}\mathrm{=}{{r}}_{\mathit{x\mu }}^{\mathit{n}}\mathrm{+}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\alpha }}_{\mathit{x}}\mathrm{)}\mathrm{𝒲ℬℰ}{\mathbf{\Delta }{T}}^{\mathrm{core}}\mathrm{+}{\mathit{\alpha }}_{\mathit{x}}\mathit{\alpha }\mathrm{𝒲ℰ}{\mathbf{\Delta }{T}}_{\mathit{x}}^{\mathrm{wing}}\mathit{,}$(28)where and $\mathrm{ℬ}\mathrm{=}\mathrm{diag}{}^{\mathrm{[}}\mathit{\beta }^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathit{,}{\mathit{\beta }}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}\mathit{,}{\mathit{\beta }}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}\mathit{,}{\mathit{\beta }}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}{}^{\mathrm{]}}$. The frequency and angle-independent four-dimensional vector ΔT^core is given by ${\mathbf{\Delta }{T}}^{\mathrm{core}}\mathrm{=}{\mathrm{\int }}_{\mathrm{-}{\mathit{x}}_{\mathrm{c}}}^{\mathrm{+}{\mathit{x}}_{\mathrm{c}}}{\mathit{\varphi }}_{{\mathit{x}}^{\mathrm{\prime }}}{\mathrm{\int }}_{-1}^{\mathrm{+}\mathrm{1}}{\mathit{p}}_{{\mathit{x}}^{\mathrm{\prime }}}{\mathbf{\Gamma }}_{{\mathit{\mu }}^{\mathrm{\prime }}}{\mathbf{\Lambda }}_{{\mathit{x}}^{\mathrm{\prime }}{\mathit{\mu }}^{\mathrm{\prime }}}^{\mathrm{\ast }}\mathrm{[}\mathit{\delta }{{\mathbf{𝒮}}_{\mathrm{l}\mathit{,}{\mathit{x}}^{\mathrm{\prime }}{\mathit{\mu }}^{\mathrm{\prime }}}^{\mathit{n}}}^{\mathrm{]}}\frac{\hspace{0.17em}\mathrm{d}{\mathit{\mu }}^{\mathrm{\prime }}}{\mathrm{2}}\hspace{0.17em}\mathrm{d}{\mathit{x}}^{\mathrm{\prime }}\mathit{,}$(29)and the frequency-dependent but angle-independent four-dimensional vector ${\mathbf{\Delta }{T}}_{\mathit{x}}^{\mathrm{wing}}$ is given by ${\mathbf{\Delta }{T}}_{\mathit{x}}^{\mathrm{wing}}\mathrm{=}{\mathrm{\int }}_{-1}^{\mathrm{+}\mathrm{1}}{\mathit{p}}_{\mathit{x}}{\mathbf{\Gamma }}_{{\mathit{\mu }}^{\mathrm{\prime }}}{\mathbf{\Lambda }}_{\mathit{x}{\mathit{\mu }}^{\mathrm{\prime }}}^{\mathrm{\ast }}\mathrm{[}\mathit{\delta }{{\mathbf{𝒮}}_{\mathrm{l}\mathit{,x}{\mathit{\mu }}^{\mathrm{\prime }}}^{\mathit{n}}}^{\mathrm{]}}\frac{\hspace{0.17em}\mathrm{d}{\mathit{\mu }}^{\mathrm{\prime }}}{\mathrm{2}}\mathrm{·}$(30)In Eq. (28), α_x is the core-wing separation coefficient. In the core α_x = 0 and in the wings ${\mathit{\alpha }}_{\mathit{x}}\mathrm{=}{\mathit{r}}_{\mathit{xx}}^{\mathrm{II}\mathrm{-}\mathrm{AA}}\mathit{/}{\mathit{\varphi }}_{\mathit{x}}$.

To evaluate ΔT^core, we consider Eq. (28) with α_x = 0 and apply the operator ${{}^{\mathrm{\int }}}_{\mathrm{-}{\mathit{x}}_{\mathrm{c}}}^{\mathrm{+}{\mathit{x}}_{\mathrm{c}}}{\mathit{\varphi }}_{\mathit{x}}{{}^{\mathrm{\int }}}_{-1}^{\mathrm{+}\mathrm{1}}{\mathit{p}}_{\mathit{x}}{\mathbf{\Gamma }}_{\mathit{\mu }}{\mathbf{\Lambda }}_{\mathit{x\mu }}^{\mathrm{\ast }}\mathrm{[}\hspace{0.17em}\mathrm{]}\frac{\hspace{0.17em}\mathrm{d}\mathit{\mu }}{\mathrm{2}}\hspace{0.17em}\mathrm{d}\mathit{x}$ from left. After some algebra, we obtain $\begin{array}{ccc}& & {\mathbf{\Delta }{T}}^{\mathrm{core}}\mathrm{=}\\ & & {\left[\mathbf{E}\mathrm{-}\left({\mathrm{\int }}_{\mathrm{-}{\mathit{x}}_{\mathrm{c}}}^{\mathrm{+}{\mathit{x}}_{\mathrm{c}}}{\mathit{\varphi }}_{\mathit{x}}{\mathrm{\int }}_{-1}^{\mathrm{+}\mathrm{1}}{\mathit{p}}_{\mathit{x}}{\mathbf{\Gamma }}_{\mathit{\mu }}{\mathbf{\Lambda }}_{\mathit{x\mu }}^{\mathrm{\ast }}\frac{\mathrm{d}\mathit{\mu }}{\mathrm{2}}\mathrm{d}\mathit{x}\right)\hspace{0.17em}\mathrm{𝒲ℬℰ}\right]}^{-1}{\overline{{r}}}^{\hspace{0.17em}\mathit{n}}\mathit{,}\end{array}$(31)where ${\overline{{r}}}^{\mathit{n}}\mathrm{=}{\mathrm{\int }}_{\mathrm{-}{\mathit{x}}_{\mathrm{c}}}^{\mathrm{+}{\mathit{x}}_{\mathrm{c}}}{\mathit{\varphi }}_{\mathit{x}}{\mathrm{\int }}_{-1}^{\mathrm{+}\mathrm{1}}{\mathit{p}}_{\mathit{x}}{\mathbf{\Gamma }}_{\mathit{\mu }}{\mathbf{\Lambda }}_{\mathit{x\mu }}^{\mathrm{\ast }}\mathrm{[}{{{r}}_{\mathit{x\mu }}^{\mathit{n}}}^{\mathrm{]}}\frac{\hspace{0.17em}\mathrm{d}\mathit{\mu }}{\mathrm{2}}\hspace{0.17em}\mathrm{d}\mathit{x}\mathit{.}$(32)Similarly, to evaluate ${\mathbf{\Delta }{T}}_{\mathit{x}}^{\mathrm{wing}}$ we apply the operator ${{}^{\mathrm{\int }}}_{-1}^{\mathrm{+}\mathrm{1}}{\mathit{p}}_{\mathit{x}}{\mathbf{\Gamma }}_{\mathit{\mu }}{\mathbf{\Lambda }}_{\mathit{x\mu }}^{\mathrm{\ast }}\mathrm{[}\hspace{0.17em}\mathrm{]}\frac{\hspace{0.17em}\mathrm{d}\mathit{\mu }}{\mathrm{2}}$ from the left to Eq. (28), to obtain $\begin{array}{ccc}& & {\mathbf{\Delta }{T}}_{\mathit{x}}^{\mathrm{wing}}\mathrm{=}{\left[\mathbf{E}\mathrm{-}{\mathit{\alpha }}_{\mathit{x}}\mathit{\alpha }\left({\mathrm{\int }}_{-1}^{\mathrm{+}\mathrm{1}}{\mathit{p}}_{\mathit{x}}{\mathbf{\Gamma }}_{\mathit{\mu }}{\mathbf{\Lambda }}_{\mathit{x\mu }}^{\mathrm{\ast }}\frac{\hspace{0.17em}\mathrm{d}\mathit{\mu }}{\mathrm{2}}\right)\hspace{0.17em}\mathrm{𝒲ℰ}\right]}^{-1}\\ & & \end{array}$(33)where ${\begin{array}{}\mathrm{˜}\\ {r}\end{array}}_{\mathit{x}}^{\mathit{n}}\mathrm{=}{\mathrm{\int }}_{-1}^{\mathrm{+}\mathrm{1}}{\mathit{p}}_{\mathit{x}}{\mathbf{\Gamma }}_{\mathit{\mu }}{\mathbf{\Lambda }}_{\mathit{x\mu }}^{\mathrm{\ast }}\mathrm{[}{{{r}}_{\mathit{x\mu }}^{\mathit{n}}}^{\mathrm{]}}\frac{\hspace{0.17em}\mathrm{d}\mathit{\mu }}{\mathrm{2}}\mathrm{·}$(34)The equations for ${\mathbf{\Delta }{T}}_{\mathit{x}}^{\mathrm{wing}}$ and ΔT^core are then incorporated into Eq. (28) to calculate the corrections to the source vector.

3.2. Scattering expansion approach

We present an iterative method of a different type that can also be used with an angle-dependent PRD. It is based on a Neumann series expansion of the components of the source vector contributing to the polarization. This Neumann series amounts to an expansion in the mean number of scattering events (see Frisch et al. 2009). Its first term yields the so-called single scattering solution. For Rayleigh scattering with an angle-dependent PRD, the single scattering solution is given in Eq. (25) of HF10. Here, following Frisch et al. (2009), we include higher order terms. The main steps of this method are given below.

We first neglect polarization in the calculation of Stokes I, i.e., we assume that Stokes I is given by the component ${\mathrm{ℐ}}_{\mathrm{0}}^{\mathrm{0}}$. This component is the solution of a non-LTE unpolarized radiative transfer equation. We calculate it with a scalar version of the core-wing ALI method described in Sect. 3.1.2, using ${{}^{\mathrm{\right(}}{\mathrm{ℛ}}_{\mathit{x}{\mathit{x}}^{\mathrm{\prime }}}^{\mathrm{AA}}^{\mathrm{\left)}}}_{\mathrm{0}}^{\mathrm{0}}$ as a PRD function. Knowing ${\mathrm{ℐ}}_{\mathrm{0}}^{\mathrm{0}}$, we can calculate the single scattering approximation for each component ${\mathrm{𝒮}}_{\mathrm{l}\mathit{,Q}}^{\mathrm{2}}$ (Q = 0,1,2). It may be written as $\begin{array}{ccc}& & {\mathrm{[}{\begin{array}{}\mathrm{˜}\\ {\mathrm{𝒮}}_{\mathrm{l}\mathit{,Q}}^{\mathrm{2}}\end{array}}^{\mathrm{]}}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}\mathrm{\left(}\mathit{\tau ,x,\mu }\mathrm{\right)}\mathrm{=}\\ & & {\mathrm{\int }}_{\mathrm{-}\mathrm{\infty }}^{\mathrm{+}\mathrm{\infty }}{\mathrm{\int }}_{-1}^{\mathrm{+}\mathrm{1}}\frac{\begin{array}{}\mathrm{˜}\\ {\mathit{R}}_{\mathit{Q}}^{\mathrm{2}}\end{array}\mathrm{\left(}\mathit{x,\mu ,}{\mathit{x}}^{\mathrm{\prime }}\mathit{,}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{\right)}}{\mathit{\varphi }\mathrm{\left(}\mathit{x}\mathrm{\right)}}\stackrel{}{{{\mathrm{\Gamma }}_{\mathit{Q}\mathrm{0}}^{\mathrm{20}}}_{\mathrm{˜}}}\mathrm{\left(}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{\right)}{\mathrm{ℐ}}_{\mathrm{0}}^{\mathrm{0}}\mathrm{\left(}\mathit{\tau ,}{\mathit{x}}^{\mathrm{\prime }}\mathit{,}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\frac{\hspace{0.17em}\mathrm{d}{\mathit{\mu }}^{\mathrm{\prime }}}{\mathrm{2}}\hspace{0.17em}\mathrm{d}{\mathit{x}}^{\mathrm{\prime }}\mathit{.}\end{array}$(35)The superscript 1 stands for single scattering. The corresponding radiation field ${{}^{\mathrm{\right[}}\begin{array}{}{\mathrm{˜}}_{}\\ {\mathrm{ℐ}}_{\mathit{Q}}^{\mathrm{2}}\end{array}^{\mathrm{\left]}}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}$ is calculated with a formal solver and it serves as a starting point for calculating the higher order terms. For order (n), $\begin{array}{ccc}& & {\mathrm{[}{\begin{array}{}\mathrm{˜}\\ {\mathrm{𝒮}}_{\mathrm{l}\mathit{,Q}}^{\mathrm{2}}\end{array}}^{\mathrm{]}}}^{\mathrm{\left(}\mathit{n}\mathrm{\right)}}\mathrm{=}{\mathrm{[}{\begin{array}{}\mathrm{˜}\\ {\mathrm{𝒮}}_{\mathrm{l}\mathit{,Q}}^{\mathrm{2}}\end{array}}^{\mathrm{]}}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}\mathrm{+}{\mathrm{\int }}_{\mathrm{-}\mathrm{\infty }}^{\mathrm{+}\mathrm{\infty }}\hspace{0.17em}\mathrm{d}{\mathit{x}}^{\mathrm{\prime }}{\mathrm{\int }}_{-1}^{\mathrm{+}\mathrm{1}}\frac{\hspace{0.17em}\mathrm{d}{\mathit{\mu }}^{\mathrm{\prime }}}{\mathrm{2}}\\ & & \mathrm{\times }\frac{\begin{array}{}\mathrm{˜}\\ {\mathrm{ℛ}}_{\mathit{Q}}^{\mathrm{2}}\end{array}\mathrm{\left(}\mathit{x,\mu ,}{\mathit{x}}^{\mathrm{\prime }}\mathit{,}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{\right)}}{\mathit{\varphi }\mathrm{\left(}\mathit{x}\mathrm{\right)}}\sum _{{\mathit{Q}}^{\mathrm{\prime }}\mathrm{\ge }\mathrm{0}}\stackrel{}{{{\mathrm{\Gamma }}_{\mathit{Q}{\mathit{Q}}^{\mathrm{\prime }}}^{\mathrm{22}}}_{\mathrm{˜}}}\mathrm{\left(}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{\right)}{\mathrm{[}{\begin{array}{}\mathrm{˜}\\ {\mathrm{ℐ}}_{{\mathit{Q}}^{\mathrm{\prime }}}^{\mathrm{2}}\end{array}}^{\mathrm{]}}}^{\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}}\mathrm{\left(}\mathit{\tau ,}{\mathit{x}}^{\mathrm{\prime }}\mathit{,}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{\right)}\mathit{.}\end{array}$(36)The iteration is continued until a convergence criterion defined in Sect. 5 is satisfied. The component ${\mathrm{ℐ}}_{\mathrm{0}}^{\mathrm{0}}$ is calculated only once. We emphasize that this method will be reliable only if the polarization rate remains small, say smaller than 20%.

4. The Fourier azimuthal averages and coefficients of the PRD functions

The numerical methods described in the preceding section involve the azimuthal Fourier coefficients $\stackrel{}{{{\mathit{r}}_{\mathrm{II}}^{\mathrm{\left(}\mathit{Q}\mathrm{\right)}}}_{\mathrm{˜}}}$ and $\stackrel{}{{{\mathit{r}}_{\mathrm{III}}^{\mathrm{\left(}\mathit{Q}\mathrm{\right)}}}_{\mathrm{˜}}}$ for Q = 0,1,2 defined in Eq. (12). The azimuthal averages with Q = 0 were considered in the classical work by Milkey et al. (1975). Subsequently Vardavas (1976), Faurobert (1987), and Wallace & Yelle (1989) devised methods to evaluate various azimuthal moments of the angle-dependent frequency redistribution functions. We stress that the moments of order Q = 0 are normalized to the absorption profile when integrated over all the incoming frequencies and angles, while the moments of orders Q = 1 and Q = 2 are normalized to zero. For reasons that we explain below, these normalizations should be satisfied to great accuracy.

Here, we use a 31-point Gauss-Legendre quadrature to perform the integration over (χ − χ′). Since r_X(x,μ,x′,μ′,χ − χ′) (X = II and III) are even functions of (χ − χ′), the integration can be limited to the range 0 ≤ (χ − χ′) ≤ π.

In Fig. 1, we show surface plots of $\stackrel{}{{{\mathit{e}}_{\mathrm{II}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}}_{\mathrm{˜}}}\mathrm{=}\stackrel{}{{{\mathit{r}}_{\mathrm{II}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}}_{\mathrm{˜}}}\mathit{/}\mathit{\varphi }\mathrm{\left(}{\mathit{x}}^{\mathrm{\prime }}\mathrm{\right)}$ and $\stackrel{}{{{\mathit{e}}_{\mathrm{III}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}}_{\mathrm{˜}}}\mathrm{=}\stackrel{}{{{\mathit{r}}_{\mathrm{III}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}}_{\mathrm{˜}}}\mathit{/}\mathit{\varphi }\mathrm{\left(}{\mathit{x}}^{\mathrm{\prime }}\mathrm{\right)}$ (left panels and right panels, respectively) as a function of the scattered frequency x and the scattered direction μ for the incoming frequencies x′ = 3 and x′ = 4, and the incoming direction μ′ = 0.3. The damping parameter a of the Voigt profile ϕ_x is taken to be equal to 0.001. The behaviors at x′ = 3 and x′ = 4 are typical of line core and line wings, respectively. We can observe that the functions $\stackrel{}{{{\mathit{e}}_{\mathrm{II}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}}_{\mathrm{˜}}}$ and $\stackrel{}{{{\mathit{e}}_{\mathrm{III}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}}_{\mathrm{˜}}}$ have similar behaviors for x′ = 3 but quite different ones for x′ = 4. For the wing frequency x′ = 4, one recovers features that are typical of the angle-averaged PRD functions, namely $\stackrel{}{{{\mathit{e}}_{\mathrm{II}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}}_{\mathrm{˜}}}$ as a function of x is peaked at the incident frequency x′, while $\stackrel{}{{{\mathit{e}}_{\mathrm{III}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}}_{\mathrm{˜}}}$ shows a CRD–type behavior with a single peak at x = 0 for all values of μ. The results shown for μ′ = 0.3 hold for all incoming directions 0 < μ′ < 1.

Figure 2 shows $\stackrel{}{{{\mathit{e}}_{\mathrm{II}}^{\mathrm{\left(}\mathit{Q}\mathrm{\right)}}}_{\mathrm{˜}}}$ and $\stackrel{}{{{\mathit{e}}_{\mathrm{III}}^{\mathrm{\left(}\mathit{Q}\mathrm{\right)}}}_{\mathrm{˜}}}$ for all the three values Q = 0,1,2 as a function of the outgoing frequency x for different choices of (μ, μ′) and different incoming frequencies x′. Our choice of (μ,μ′) values for Fig. 2 are identical to those of Wallace & Yelle (1989) in their analysis of the azimuth-averaged type II redistribution function. The frequencies x′ = 1 and x′ = 4 are representative of the line center and wing behaviors, respectively. Figure 2 clearly shows that the sharp peaks appearing in Fig. 1 for Q = 0 also exist for Q = 1,2 for forward and backward scattering (see panels b, c, e, f). Actually, surface plots for Q = 1,2 are very similar to the surface plots shown in Fig. 1. However, because they are normalized to zero, $\stackrel{}{{{\mathit{e}}_{\mathrm{II}}^{\mathrm{\left(}\mathit{Q}\mathrm{\right)}}}_{\mathrm{˜}}}$ and $\stackrel{}{{{\mathit{e}}_{\mathrm{III}}^{\mathrm{\left(}\mathit{Q}\mathrm{\right)}}}_{\mathrm{˜}}}$ for Q = 1,2 can take negative values as can be observed in Fig. 2. In the upper panels of this figure, one can observe that the coefficients $\stackrel{}{{{\mathit{e}}_{\mathrm{II}}^{\mathrm{\left(}\mathit{Q}\mathrm{\right)}}}_{\mathrm{˜}}}$ decrease in magnitude with increasing Q. This result agrees with the curves shown by Domke & Hubeny (1988) for μ = 0.8 and μ′ = 0.2 (see their Fig. 3). However, for the special cases of forward and backward scattering (panels b, c, e, f), the Fourier coefficients for Q = 1 and 2 can become larger than the Q = 0 coefficients.

A significant difference between angle-averaged and angle-dependent PRD functions is the presence of very sharp and narrow peaks for x ≃ x′ and μ ≃ μ′ for the angle-dependent PRD functions. This has important implications for the numerical evaluation of the scattering integrals (see Eq. (6)). It is well known in scalar non-LTE radiative transfer calculations that an accurate normalization is required for the profile function and the PRD functions, especially for lines with a large optical thickness and a very small thermalization parameter. Any error will indeed act as a spurious sink or source of photons. To achieve the normalization to a high accuracy, a very fine frequency grid is needed. However, the use of such fine frequency grids in the transfer calculations require large computing resources. Hence, strategies have been developed in the past (see for e.g., Wallace & Yelle 1989) to handle frequency quadratures in angle-dependent PRD problems. For unpolarized transfer problems, Adams et al. (1971) proposed a natural cubic spline representation for the radiation field. In this paper, we have developed our own strategy to maintain the computational cost of the radiative transfer calculations at a reasonable level while ensuring a high accuracy: $\stackrel{}{{{\mathit{r}}_{\mathrm{X}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}}_{\mathrm{˜}}}$ is normalized to ϕ_x with an accuracy of 99% and the integrals of $\stackrel{}{{{\mathit{r}}_{\mathrm{X}}^{\mathrm{\left(}\mathit{Q}\mathrm{\right)}}}_{\mathrm{˜}}}$, Q = 1,2, over incoming frequencies and directions are around 10^-7.

We first start with a frequency grid typical of PRD line transfer problems on which the transfer equation will be solved. The choice of the actual frequency grid is given below in Sect. 5. We then sub-divide each frequency interval into a fine mesh of Simpson quadrature points (e.g., 41-point Simpson formula) on which we calculate the coefficients $\stackrel{}{{{\mathit{r}}_{\mathrm{X}}^{\mathrm{\left(}\mathit{Q}\mathrm{\right)}}}_{\mathrm{˜}}}$. A seven-point Gaussian quadrature formula with μ in the range 0 < μ ≤ 1 is used for the angular grid.

To handle the peaks that occur in the forward and backward scattering situations (see e.g. Fig. 1) we proceed in the following way. We introduce a cut-off scattering angle Θ_{cut − off} = 10^-6 radians and assume that the PRD functions keep a constant value, given by the values at the cut-off, when Θ < Θ_{cut − off} or (π − Θ) < Θ_{cut − off}. This practical trick has implications for the normalization accuracy of the redistribution functions, but we have verified that Θ_{cut − off} = 10^-6 radians is a reasonable choice. To approach the two extreme values of Θ as close as possible, it is sufficient to employ a seven-point Gaussian quadrature in  [0 < μ ≤ 1] .

5. Validation and convergence properties of the iterative methods

We consider isothermal, self-emitting plane-parallel slab atmospheres with no incident radiation at the boundaries. These slab models are characterized by a set of input parameters (T,a,ϵ,r,Γ_E/Γ_R), where T is the optical thickness of the slab. The Planck function is set to unity. The depolarizing collisional rate D⁽²⁾ is assumed to be 0.5 × Γ_E. The thermalization parameter ϵ, which is actually a photon destruction probability, is defined by ϵ = Γ_I/(Γ_R + Γ_I) = 1 − β⁽⁰⁾. The PRD function defined in Eq. (8) can also be written as $\stackrel{}{{{\mathit{R}}_{\mathrm{0}}^{\mathrm{0}}}_{\mathrm{˜}}}\mathrm{=}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{ϵ}\mathrm{)}\mathrm{[}{\mathit{\gamma }}_{\mathrm{coh}}\stackrel{}{{{\mathit{r}}_{\mathrm{II}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}}_{\mathrm{˜}}}\mathrm{+}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\gamma }}_{\mathrm{coh}}\mathrm{\left)}\stackrel{}{{{\mathit{r}}_{\mathrm{III}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}}_{\mathrm{˜}}}\mathrm{\right]}\mathit{,}$(37)where γ_coh = (1 + Γ_E/(Γ_R + Γ_I))^-1. For Γ_E = 0, one has pure r_II since γ_coh = 1. For small values of ϵ, γ_coh is about 1/(1 + Γ_E/Γ_R). In all our calculations, W₂ = μ₂ = 1.

For all the figures presented in this paper, we use a logarithmically spaced τ-grid with 5 points per decade, with the first depth point at τ₁ = 10^-2. For the frequency grid, we use equally spaced points in the line core and logarithmically spaced ones in the wings. Furthermore, the maximum frequency x_max is chosen such that the condition ϕ(x_max)T ≪ 1 is satisfied. We have typically 70 points in the interval  [0,x_max] .

To test the correctness of our three iterative methods, we compared the emergent solutions with those of the perturbative type method developed in Nagendra et al. (2002), for an optically thin (T = 10) slab and a relatively thicker (T = 10³) one. In Nagendra et al. (2002), the radiation field is represented by the two Stokes parameters I and Q and there is no azimuthal Fourier decomposition of the angle-dependent PRD functions. For the thin slab, the model parameters are (T,a,ϵ,r,Γ_E/Γ_R) = (10,10^-3,10^-4,0,1) and for the thick slab they are (T,a,ϵ,r,Γ_E/Γ_R) = (10³,10^-3,10^-4,0,0). We have found that our new iterative methods yield emergent solutions that are in very good agreement with the results of the perturbation method used in Nagendra et al. (2002, differences in the ratio Q/I are at most 6%).

To study the convergence properties of the core-wing, FABFA, and scattering expansion methods, we introduce the quantities $\mathrm{(}{\mathit{c}}_{\mathit{Q}}^{\mathit{K}}{\mathrm{)}}^{\mathit{n}}\mathrm{=}{\max }_{{\mathit{\tau }}_{\mathrm{d}}\mathit{,x,\mu }}\begin{array}{}\mathrm{⎧}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎨}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎩}\end{array}\frac{\mathrm{\left|}\mathrm{\right(}{\mathrm{𝒮}}_{\mathit{Q}}^{\mathit{K}}{\mathrm{)}}_{\mathit{x\mu }}^{\mathit{n}\mathrm{+}\mathrm{1}}\mathrm{\left(}{\mathit{\tau }}_{\mathrm{d}}\mathrm{\right)}\mathrm{-}\mathrm{\left(}{\mathrm{𝒮}}_{\mathit{Q}}^{\mathit{K}}{\mathrm{)}}_{\mathit{x\mu }}^{\mathit{n}}\mathrm{\right(}{\mathit{\tau }}_{\mathrm{d}}\mathrm{\left)}\mathrm{\right|}}{\mathrm{\left(}\overline{\mathrm{𝒮}}\begin{array}{}\mathit{K}\\ \mathit{Q}\end{array}{\mathrm{)}}_{\mathit{x\mu }}^{\mathit{n}\mathrm{+}\mathrm{1}}\mathrm{\right(}{\mathit{\tau }}_{\mathrm{d}}\mathrm{)}}\begin{array}{}\mathrm{⎫}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎬}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎭}\end{array}\mathit{,}$(38)with n the iteration index, τ_d a depth-grid point, and $\mathrm{\left(}\overline{\mathrm{𝒮}}\begin{array}{}\mathit{K}\\ \mathit{Q}\end{array}{\mathrm{)}}_{\mathit{x\mu }}^{\mathit{n}\mathrm{+}\mathrm{1}}\mathrm{\right(}{\mathit{\tau }}_{\mathrm{d}}\mathrm{)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{\left[}\mathrm{\right|}\mathrm{\left(}{\mathrm{𝒮}}_{\mathit{Q}}^{\mathit{K}}{\mathrm{)}}_{\mathit{x\mu }}^{\mathit{n}\mathrm{+}\mathrm{1}}\mathrm{\right(}{\mathit{\tau }}_{\mathrm{d}}\mathrm{\left)}\mathrm{\right|}\mathrm{+}\mathrm{\left|}\mathrm{\right(}{\mathrm{𝒮}}_{\mathit{Q}}^{\mathit{K}}{\mathrm{)}}_{\mathit{x\mu }}^{\mathit{n}\mathrm{+}\mathrm{1}}\mathrm{\left(}{\mathit{\tau }}_{\mathit{d}\mathrm{+}\mathrm{1}}\mathrm{\right)}{\mathrm{\left|}}^{\mathrm{\right]}}\mathit{.}$(39)The $\mathrm{(}{\mathit{c}}_{\mathit{Q}}^{\mathit{K}}{\mathrm{)}}^{\mathit{n}}$ yield the maximum relative change from one iteration to the next for each component ${\mathrm{𝒮}}_{\mathit{Q}}^{\mathit{K}}$. We also introduce ${\mathit{c}}_{\mathrm{p}}^{\mathit{n}}\mathrm{=}{\max }_{\mathit{x,\mu }}\begin{array}{}\mathrm{⎧}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎨}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎩}\end{array}\frac{\mathrm{|}{\mathit{P}}_{\mathit{x\mu }}^{\mathit{n}\mathrm{+}\mathrm{1}}\mathrm{-}{\mathit{P}}_{\mathit{x\mu }}^{\mathit{n}}\mathrm{|}}{\mathrm{\left|}{\mathit{P}}_{\mathit{x\mu }}^{\mathit{n}\mathrm{+}\mathrm{1}}\mathrm{\right|}}\begin{array}{}\mathrm{⎫}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎬}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎭}\end{array}\mathit{,}$(40)where P = Q/I is the degree of linear polarization at the surface. The quantity ${\mathit{c}}_{\mathrm{p}}^{\mathit{n}}$ measures the progress of iteration on the surface polarization of the radiation field. For the ALI type methods, the convergence is tested using the criterion ${\mathit{c}}^{\mathit{n}}\mathrm{=}\max \mathrm{\left[}\mathrm{\right(}{\mathit{c}}_{\mathrm{0}}^{\mathrm{0}}{\mathrm{)}}^{\mathit{n}}\mathit{,}{{\mathit{c}}_{\mathrm{p}}^{\mathit{n}}}^{\mathrm{]}}\mathit{<}{\mathrm{10}}^{-8}\mathit{.}$(41)In the asymptotic regime of large n, the components ${\mathrm{𝒮}}_{\mathit{Q}}^{\mathit{K}}$ are slave modes of the component ${\mathrm{𝒮}}_{\mathrm{0}}^{\mathrm{0}}$, so all the $\mathrm{(}{\mathit{c}}_{\mathit{Q}}^{\mathit{K}}{\mathrm{)}}^{\mathit{n}}$ have more or less the same behavior as $\mathrm{(}{\mathit{c}}_{\mathrm{0}}^{\mathrm{0}}{\mathrm{)}}^{\mathit{n}}$ (see Nagendra et al. 1998). For the scattering expansion method, we perform two iteration tests, a first one on $\mathrm{(}{\mathit{c}}_{\mathrm{0}}^{\mathrm{0}}{\mathrm{)}}^{\mathit{n}}$ when calculating the component ${\mathrm{ℐ}}_{\mathrm{0}}^{\mathrm{0}}$ with an ALI method (first stage), and then a second one on ${\mathit{c}}_{\mathrm{p}}^{\mathit{n}}$, when calculating the components ${\mathrm{ℐ}}_{\mathit{Q}}^{\mathrm{2}}$ (second stage) with the iteration formula given in Eq. (36). The iterations are stopped when $\mathrm{(}{\mathit{c}}_{\mathrm{0}}^{\mathrm{0}}{\mathrm{)}}^{\mathit{n}}$ and ${\mathit{c}}_{\mathrm{p}}^{\mathit{n}}$ become smaller than 10^-8.

Figure 3 shows the variations in $\mathrm{(}{\mathit{c}}_{\mathrm{0}}^{\mathrm{0}}{\mathrm{)}}^{\mathit{n}}$ and ${\mathit{c}}_{\mathrm{p}}^{\mathit{n}}$ versus the iteration number. The model parameters used are (T,a,ϵ,r,Γ_E/Γ_R) = (2 × 10⁴,10^-3,10^-4,0,0). The panel (a) corresponds to the FABFA method, and the panel (b) to the core-wing method. The convergence behavior of these two new ALI methods is nearly similar to the standard Polarized ALI (PALI) methods (Nagendra et al. 1999). The approximation introduced for the core-wing method (see Eq. (22)) does not seem to affect the speed of convergence significantly. We also find that the speeds of convergence are about the same for angle-averaged and angle-dependent PRD. The main parameters that can affect the speed of convergence are the optical thickness of the line (there being faster convergence for optically thin than optically thick lines), and the frequency grid, coarser grids leading to faster convergence (Olson et al. 1986). The panel (c) of Fig. 3 shows the variation in $\mathrm{(}{\mathit{c}}_{\mathrm{0}}^{\mathrm{0}}{\mathrm{)}}^{\mathit{n}}$ and ${\mathit{c}}_{\mathrm{p}}^{\mathit{n}}$ for the scattering expansion method. The number of iterations in the second iteration stage is clearly small (~60) compared to the first iteration stage (~140).

To illustrate the convergence properties of the ALI methods, we show in Fig. 4 the convergence history of the four components of at x = 0 calculated with the core-wing method. To reduce the number of lines, every fourth iteration solutions are plotted. We note that ${\mathrm{𝒮}}_{\mathrm{0}}^{\mathrm{0}}$ for x = 0 satisfies the $\sqrt{\mathit{ϵ}}$ law at the surface. As the slab is not optically very thick (T = 2 × 10⁴), ${\mathrm{𝒮}}_{\mathrm{0}}^{\mathrm{0}}$ has not fully reached unity at the mid-slab. All the other components ${\mathrm{𝒮}}_{\mathit{Q}}^{\mathrm{2}}$ go to zero at the mid slab. In the wing frequencies, say x = 4, the rate of convergence for all the four components of is quite large (figure not shown here).

Figure 5 shows the convergence history of the components ${\mathrm{𝒮}}_{\mathit{Q}}^{\mathrm{2}}$ for x = 0 and μ = 0.025 calculated with the scattering expansion method, for the same atmospheric model as in Figs. 3 and 4. The dotted lines indicate the single scattered solution. For ${\mathrm{𝒮}}_{\mathrm{1}}^{\mathrm{2}}$ and ${\mathrm{𝒮}}_{\mathrm{2}}^{\mathrm{2}}$, the single scattered solution is very close to the converged solution, while for ${\mathrm{𝒮}}_{\mathrm{0}}^{\mathrm{2}}$ a few iterations are needed to reach the converged solution. In the last panel of Fig. 5, we show the convergence history for the ratio Q/I. We have not performed a systematic investigation of the variation in the speed of convergence with the slab thickness as was done for CRD in Frisch et al. (2009), but the calculations that we have performed suggest that the behavior observed with CRD will also hold for angle-dependent and angle-averaged PRD, namely, a very good approximation to the emergent polarization is provided by the single scattered solution for optically thin slabs (T ≪ 10) and optically very thick ones (T ≥ 2 × 10⁶), but is not reliable for intermediate optical thicknesses. As shown here, the problem can be cured by considering higher order terms in the scattering expansion. For example, for a very strong line such as Ca i 4227 Å  the Q/I profile can be fitted fairly well with a single scattered solution (Anusha et al. 2010), although the line core is somewhat underestimated. The fit could probably be improved by considering higher order terms in the scattering expansion.

In terms of computing resources, in particular computing time, the three iterative methods behave quite differently. In Table 1, we present the CPU time requirements for the three iterative methods in the case of angle-dependent PRD. The CPU times spent in specific parts of the computational process are listed, for each of the methods. The model used is the same as in Figs. 3–5. Since the computation of the azimuthal Fourier coefficients $\stackrel{}{{{\mathit{r}}_{\mathrm{X}}^{\mathrm{\left(}\mathit{Q}\mathrm{\right)}}}_{\mathrm{˜}}}$ are common to all the three iterative methods, they are excluded when estimating the CPU time requirements. The computing time for $\stackrel{}{{{\mathit{r}}_{\mathrm{X}}^{\mathrm{\left(}\mathit{Q}\mathrm{\right)}}}_{\mathrm{˜}}}$ can be huge (several hours) depending on the frequency and angle grids. For each method, the last but one column in Table 1 gives the number of iterations multiplied by the CPU time per iteration. This number of iterations is smaller than the results presented in Fig. 3 because we have applied the Ng acceleration to all the three iterative methods.

The scattering expansion method requires the construction of a scalar Λ^∗ operator needed to calculate ${\mathrm{ℐ}}_{\mathrm{0}}^{\mathrm{0}}$ in the first stage of the method. The computation of Λ^∗, ${\mathrm{ℐ}}_{\mathrm{0}}^{\mathrm{0}}$, and the single scattered solution $\mathrm{[}\begin{array}{}{}_{\mathrm{˜}}\\ {\mathrm{ℐ}}_{\mathit{Q}}^{\mathrm{2}}\end{array}{\mathrm{]}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}$ requires roughly 3/20th of the total computing time. Thus, most of the time is spent in calculating the higher order terms in the scattering expansion. For the core-wing method, almost the entire CPU time is spent in the iterative cycle itself. Each iterative cycle is about 20% shorter than the time per iteration of the scattering expansion method. In addition, a larger number of iterations are needed for convergence. It is the need to invert the large matrix A^-1 (see Eq. (19)) that considerably slows down the FABFA method. We also note that the time per iteration is about three times as long as with the core-wing method. The CPU times given in Table 1 correspond to the pure R_II case (Γ_E/Γ_R = 0). The inclusion of collisions does not change the CPU times for different stages of computational process. The only difference occurs in the number of iterations required for convergence. The larger the value of Γ_E/Γ_R, the smaller the total number of iterations.

For the core-wing and the scattering expansion methods, the computing time per iteration scales as N_d(N_x2N_μ), where N_d is the total number of depth points. For the FABFA method, the total computing time scales as N_d(N_x2N_μ)². For the FABFA and core-wing methods, the number of iterations can be reduced by using a Gauss-Seidel decomposition or SOR method instead of a Jacobi decomposition. In any case, the FABFA method will remain slow, by a large factor, compared to the core-wing method. We recall that for the perturbation method used in Nagendra et al. (2002), the computing time per iteration scales as N_d(N_x2N_μN_ϕ), with N_ϕ the number of azimuths needed to describe the azimuthal variation of the Stokes source vector. The number of iterations required for convergence is around 20.

6. Stokes profiles calculated with angle-dependent and angle-averaged PRD functions

We compare the Stokes parameters I, Q, and the ratio Q/I calculated with angle-averaged and angle-dependent PRD functions, for several atmospheric models. The Stokes parameters calculated with angle-dependent PRD functions are shown as solid lines and those calculated with angle-averaged PRD functions as dashed lines. The results are analyzed with the help of the ${\mathrm{ℐ}}_{\mathit{Q}}^{\mathit{K}}$ decomposition, which may be written as $\begin{array}{ccc}& & \mathit{I}\mathrm{\left(}\mathit{\tau ,x,\mu }\mathrm{\right)}\mathrm{=}{\mathrm{ℐ}}_{\mathrm{0}}^{\mathrm{0}}\\ & & \mathrm{+}\frac{\mathrm{1}}{\mathrm{2}\sqrt{\mathrm{2}}}\mathrm{(}\mathrm{3}{\mathit{\mu }}^{\mathrm{2}}\mathrm{-}\mathrm{1}\mathrm{)}{\mathrm{ℐ}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{-}\frac{\sqrt{\mathrm{3}}}{\mathrm{2}}\mathit{\mu }\sqrt{\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}}{\mathrm{ℐ}}_{\mathrm{1}}^{\mathrm{2}}\mathrm{+}\frac{\sqrt{\mathrm{3}}}{\mathrm{4}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}\mathrm{)}{\mathrm{ℐ}}_{\mathrm{2}}^{\mathrm{2}}\mathit{,}\\ & & \mathit{Q}\mathrm{\left(}\mathit{\tau ,x,\mu }\mathrm{\right)}\mathrm{=}\\ & & \mathrm{-}\frac{\mathrm{3}}{\mathrm{2}\sqrt{\mathrm{2}}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}\mathrm{)}{\mathrm{ℐ}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{-}\frac{\sqrt{\mathrm{3}}}{\mathrm{2}}\mathit{\mu }\sqrt{\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}}{\mathrm{ℐ}}_{\mathrm{1}}^{\mathrm{2}}\mathrm{-}\frac{\sqrt{\mathrm{3}}}{\mathrm{4}}\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{\mu }}^{\mathrm{2}}\mathrm{)}{\mathrm{ℐ}}_{\mathrm{2}}^{\mathrm{2}}\mathit{.}\end{array}$(42)We recall that the irreducible components ${\mathrm{ℐ}}_{\mathit{Q}}^{\mathit{K}}$ depend on τ, x, and μ. The contributions of the component ${\mathrm{ℐ}}_{\mathrm{1}}^{\mathrm{2}}$ go to zero at both the limb and the disk center, for both Stokes I and Stokes Q.

In Fig. 6, we show the emergent I, Q/I, and Q profiles at μ = 0.11 for angle-dependent and angle-averaged PRD functions. The angle-averaged solution is calculated with the PALI method of Fluri et al. (2003) and the angle-dependent solution with the core-wing method introduced here. In the left-hand side panels, T = 10 and the PRD used is an equal mixture of r_II and r_III. In the right-hand side panels, T = 2 × 10⁴ and the PRD used is pure r_II. Figure 6 clearly shows that the angle-averaged and angle-dependent solutions may differ significantly and that the angle-averaged emergent polarization may be smaller or larger than the angle-dependent one. We now try to explain the reasons for this, by considering the components ${\mathrm{ℐ}}_{\mathit{Q}}^{\mathit{K}}$.

In Fig. 7, we show the components ${\mathrm{ℐ}}_{\mathit{Q}}^{\mathit{K}}$ corresponding to the two atmospheric models in Fig. 6. The components with Q ≠ 0 are zero for the angle-averaged PRD functions. The components ${\mathrm{ℐ}}_{\mathrm{0}}^{\mathrm{0}}$ are essentially equal for the AD and AA cases when T = 10, while they slightly differ around x = 3 when T = 2 × 10⁴. Since Stokes I is dominated by ${\mathrm{ℐ}}_{\mathrm{0}}^{\mathrm{0}}$, it is also nearly independent of the choice of the PRD function. Stokes Q consists of the three components ${\mathrm{ℐ}}_{\mathrm{0}}^{\mathrm{2}}$, ${\mathrm{ℐ}}_{\mathrm{1}}^{\mathrm{2}}$, and ${\mathrm{ℐ}}_{\mathrm{2}}^{\mathrm{2}}$. We can observe that ${\mathrm{ℐ}}_{\mathrm{0}}^{\mathrm{2}}$ has very different values in the AD and AA cases, being smaller in absolute value in the AD case. This is a somewhat unexpected result. The reason for the difference between the AA and AD results is the use of $\begin{array}{}{}_{\mathrm{˜}}\\ {\mathrm{ℛ}}_{\mathrm{0}}^{\mathrm{2}}\end{array}\mathrm{\left(}\mathit{x,\mu ,}{\mathit{x}}^{\mathrm{\prime }}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{\right)}$ in the AD case and the angle-averaged redistribution function $\mathrm{(}{\mathrm{ℛ}}_{\mathit{x}{\mathit{x}}^{\mathrm{\prime }}}^{\mathit{AA}}{\mathrm{)}}_{\mathrm{0}}^{\mathrm{2}}$ in the AA case.

Equation (42) shows that for small values of μ, only ${\mathrm{ℐ}}_{\mathrm{0}}^{\mathrm{2}}$ and ${\mathrm{ℐ}}_{\mathrm{2}}^{\mathrm{2}}$ contribute to Stokes Q and that $\mathit{Q}\mathrm{\simeq }\mathrm{-}\mathrm{(}\sqrt{\mathrm{3}}\mathit{/}\mathrm{4}\mathrm{)}\mathrm{(}\sqrt{\mathrm{6}}{\mathrm{ℐ}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{+}{\mathrm{ℐ}}_{\mathrm{2}}^{\mathrm{2}}\mathrm{)}$. For T = 10, the two terms have the same sign and when added together give a value of Stokes Q that is slightly smaller in the AD case than in the AA case. For spectral lines with a small optical thickness, the role of PRD is essentially negligible. These lines can be modeled reasonably well with CRD (Sampoorna et al. 2010). Thus, one may not expect large differences between angle-averaged and angle-dependent polarization profiles. For T = 2 × 10⁴, ${\mathrm{ℐ}}_{\mathrm{0}}^{\mathrm{2}}$ and ${\mathrm{ℐ}}_{\mathrm{2}}^{\mathrm{2}}$ are mainly negative, with ${\mathrm{ℐ}}_{\mathrm{0}}^{\mathrm{2}}$ having a slightly larger absolute value than ${\mathrm{ℐ}}_{\mathrm{2}}^{\mathrm{2}}$. When added together, they give a negative value, but because of the minus sign, one ends up with a positive Stokes Q, which is however smaller than the corresponding AA Stokes Q.

We now examine how the differences between the AD and AA polarizations vary with the position on the disk, the optical thickness T of the slab, the thermalization parameter ϵ, the value of the continuum absorption, and the ratio Γ_E/Γ_R.

6.1. Center-to-limb variations

We know that the angle-dependent PRD functions become azimuthally symmetric at the disk center. Thus we can expect the AD and AA polarization profiles to become very close to each other when μ → 1. This is indeed what we observe in Fig. 8 where we show the emergent profiles of Stokes I and the ratio Q/I. In the bottom three panels, we show the absolute difference between the AA and AD Q/I profiles. As μ → 1, the absolute difference clearly goes to zero.

We show in Fig. 9 the dependence on μ of the components ${\mathrm{ℐ}}_{\mathit{Q}}^{\mathit{K}}$ and ${\mathrm{𝒮}}_{\mathit{Q}}^{\mathit{K}}$ at the surface. For ${\mathrm{𝒮}}_{\mathit{Q}}^{\mathit{K}}$, the dependence on μ comes from the AD redistribution functions $\begin{array}{}{}_{\mathrm{˜}}\\ {\mathrm{ℛ}}_{\mathit{Q}}^{\mathit{K}}\end{array}\mathrm{\left(}\mathit{x,\mu ,}{\mathit{x}}^{\mathrm{\prime }}\mathit{,}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{\right)}$ (see Eq. (13)). We recall that the source vector components are independent of μ for AA-PRD. For the components ${\mathrm{ℐ}}_{\mathit{Q}}^{\mathit{K}}$, the variations with μ are caused by the limb-darkening and the μ-variation in the source vector components. In the right panels of Fig. 9, we can observe that ${\mathrm{𝒮}}_{\mathrm{0}}^{\mathrm{0}}$ is almost independent of μ and we can verify that ${\mathrm{𝒮}}_{\mathit{Q}}^{\mathrm{2}}$, Q = 1,2, go to zero when μ is close to one. The variation in ${\mathrm{𝒮}}_{\mathrm{2}}^{\mathrm{2}}$ is rather monotonic, but that of ${\mathrm{𝒮}}_{\mathrm{1}}^{\mathrm{2}}$ is not. The component ${\mathrm{𝒮}}_{\mathrm{0}}^{\mathrm{2}}$ increases towards the disk center, as it is controlled by the magnitude of ${\mathrm{ℐ}}_{\mathrm{0}}^{\mathrm{0}}$ (see dotted lines in ${\mathrm{𝒮}}_{\mathrm{0}}^{\mathrm{2}}$ panel of Fig. 5). In the left panels of Fig. 9, one can verify that ${\mathrm{ℐ}}_{\mathrm{0}}^{\mathrm{0}}$ increases from the limb to the disk center, while the components ${\mathrm{ℐ}}_{\mathit{Q}}^{\mathrm{2}}$, Q = 1,2, go to zero at the disk center.

6.2. Effects of the optical thickness T and thermalization parameter ϵ

In Fig. 10, we compare the angle-averaged and the corresponding angle-dependent Stokes I, Q, and Q/I profiles for slabs with different values of the optical thickness T and the thermalization parameter ϵ. For the angle-dependent case, the calculations were performed with the core-wing method. For all the examples shown in Fig. 10, the optical thickness is equal to or larger than 2 × 10⁴. One can observe that the ratio Q/I is always smaller for the angle-dependent than for the angle-averaged PRD functions. This situation appears to be typical of optically thick lines in isothermal atmospheres (see also Nagendra et al. 2002).

In the left panels, ϵ = 10^-4. The three values chosen for T correspond to effectively thin (ϵT = 2), effectively thick (ϵT = 200), and semi-infinite like conditions (ϵT = 2 × 10⁴). The amplitudes of Q/I and Q in the near wing peaks decrease with increasing T. For Q/I, this is accompanied by a decrease in the differences between the angle-averaged and angle-dependent values. In the particular case of T = 2 × 10⁸, typical features of semi-infinite atmospheres can be observed. For example, the appearance of double peaks – one narrow peak in the near wings, and a second broader one in the far wings – is a typical behavior of Q/I in strong resonance lines such as that of the Ca ii K line (see e.g., Stenflo 1980; Holzreuter et al. 2006).

In the right panels of Fig. 10, the optical thickness is T = 2 × 10⁴. As ϵ increases from 10^-6 (effectively thin, ϵT = 2 × 10^-2) to 10^-2 (effectively thick, ϵT = 2 × 10²), the magnitudes of both I and Q increase because the number of photons that are emitted increases, as does the strength of the coupling between the radiation field and the thermal pool. However, the ratio Q/I decreases because the number of polarized photons is reduced by the thermalization of the radiation field (see e.g., Saliba 1986). For ϵ = 10^-6, the angle-averaged and angle-dependent profiles display significant differences in the Q/I profiles. This difference decreases when ϵ = 10^-2.

We also considered the case ϵ = 0 (conservative scattering). In this case, there are no internal sources of photons (ϵB = 0), but there is an incident field at the lower boundary I(τ = T,x,μ) = B. For this model, the primary source of photons G₀(τ) decreases towards the surface at τ = 0. For the emergent Stokes I and Q profiles, the differences between the angle-averaged and angle-dependent cases are insignificant. Some differences appear however in the ratio Q/I around the line center (0 < x < 4).

6.3. Effect of the continuum strength parameter r

All the results shown in the previous figures were obtained with a continuum absorption set to zero. In Fig. 11, we show the Q/I profile for values of the continuum strength parameter r set to 10^-8, 10^-6, and 10^-4. The atmospheric model is the same as in Fig. 6b except for the value of r. When r increases, one can observe a significant decrease in amplitude of the near wing peak, an effect already investigated in Faurobert (1988), and a gradual disappearance of the differences between the angle-averaged and angle-dependent Q/I profiles. For r = 10^-2 (result not shown here), the differences become essentially zero.

6.4. Effect of the elastic collisions Γ_E

A change in the value of Γ_E, modifies the coherency factor γ_coh defined in Sect. 5, hence the relative contributions of r_II and r_III. It is well established that the wings of Stokes I and the linear polarization are quite sensitive to the PRD mechanism (Nagendra et al. 2002). We show in Fig. 12 the Stokes parameters I and Q, and the ratio Q/I calculated with the angle-averaged and angle-dependent PRD functions with the core-wing method. We show two examples, Γ_E/Γ_R = 1, which corresponds to an even mixture of r_II and r_III, and Γ_E/Γ_R = 10, which corresponds to an uneven mixture of r_II and r_III with a dominant contribution from r_III. For D⁽²⁾, we assume the relation D⁽²⁾ = 0.5Γ_E. We can observe that the polarization rate decreases as Γ_E/Γ_R increases, a phenomenon discussed in detail in Nagendra et al. (2002). The differences between the AA and AD Stokes I and Q profiles decrease as Γ_E/Γ_R increases. For Γ_E/Γ_R > 100, these differences become insignificant.

7. Conclusions

We have investigated the Rayleigh scattering with angle-dependent partial redistribution (PRD) functions to evaluate the reliability of the usual angle-averaged approximation. The analysis has been carried out for a two-level atom with unpolarized ground level using the angle-dependent PRD functions established in Domke & Hubeny (1988); Bommier (1997a). The scattering medium has been assumed to be a plane-parallel, cylindrically symmetrical, isothermal slab. For this model, the polarized radiation field is cylindrically symmetrical (i.e. depends only on the inclination θ with respect to the direction perpendicular to the slab surface) and can be represented by the two Stokes parameters I and Q. The source terms in the transfer equations for I and Q have a complicated dependence on the angle θ, because it appears in the polarization phase matrix and also in the angle-dependent PRD functions. It was shown in HF10 that the polarized radiation field can be decomposed into four components that are cylindrically symmetrical and satisfy standard radiative transfer equations. The source terms still depend on θ but in a much simpler way. Thanks to this decomposition method, we have been able to construct three different iterative methods of solution for the calculation of the Stokes parameters. Two are of the accelerated lambda iteration (ALI) type; they are generalizations of the frequency-by-frequency (FBF) and core-wing methods originally formulated for scalar PRD transfer in Paletou & Auer (1995). The other is a scattering expansion (Neumann series) method developed for the Hanle effect with complete frequency redistribution in Frisch et al. (2009).

Crucial ingredients for the three methods are the azimuthal Fourier coefficients of order 0,1, and 2 of the angle-dependent Hummer’s (1962) PRD functions r_II and r_III. The calculations of these Fourier coefficients is quite time consuming, in particular for r_III, because very fine grids are needed to properly represent their sharp variations with the frequency and inclination of the incident and scattered beams. Once these coefficients have been calculated, and also an approximate solution of Stokes I (neglecting polarization) is obtained in the scattering expansion method, the core-wing and scattering expansion approaches both provide fast and accurate solutions. The scattering expansion method appears particularly interesting in problems of large dimensionality (very large number of frequency, angle, and depth grid points).

We have found that the angle-averaged PRD functions overestimate the emergent polarization rate Q/I between 10% and 30% for slabs with a large optical thickness (between 2 × 10⁴ and 2 × 10⁸), the largest differences occurring in the near wing peaks. For optically thin slabs, the differences are much smaller. We have not considered semi-infinite atmospheres. As the difference between the polarization rates obtained with angle-averaged and angle-dependent PRD functions remain fairly small, the sign of this difference may easily depend on the particular choice of the atomic and atmospheric model. It would certainly be interesting to consider realistic solar atmosphere models as has recently been done in Sampoorna et al. (2010). The numerical methods we have proposed here are able to do this. This type of analysis would certainly help us to understand the discrepancies between observed polarization rates and theoretical predictions (Anusha et al. 2010).

For the weak-field Hanle effect, the reliability of the angle-averaged approximation for PRD functions is clearly questionable, for Stokes Q and especially Stokes U (Nagendra et al. 2002). In the presence of a weak magnetic field, the polarized radiation field can be decomposed into a set of six irreducible components. For an angle-dependent PRD, these components have no cylindrical symmetry, but a Fourier azimuthal expansion allows one to construct an infinite set of integral equations similar to the equations given here (Frisch 2009). The core-wing or a scattering expansion method applied to a truncated set of these equations offer some hope in calculating the polarization with a reasonable amount of numerical work.

---

Acknowledgments

The authors are grateful to Dr. V. Bommier for providing a FORTRAN PROGRAM to compute angle-dependent type III redistribution function and to L. S. Anusha for helpful discussions.

References

All Tables

All Figures