© ESO 2020

1. Introduction

The non-linear radiative transfer problem out of local thermodynamic equilibrium (hereafter, non-LTE) must be solved by an iterative solution of the equations of statistical equilibrium (ESE) and the radiative transfer equation (RTE). In most applications in astrophysics, the numerical solution of the RTE is the most computationally demanding part of the problem. This is especially true for the non-axisymmetric and generally multi-dimensional problems.

In the radiative transfer codes based on the concept of short or long characteristics on Cartesian grids, the discrete set of rays along which the formal solution of the RTE is calculated are the angular quadrature points on a unit sphere. The values of Stokes parameters at the quadrature points are used, with appropriate weights, to obtain the scalar mean intensity J or the radiation field tensors $J_Q^K$ in the general polarisation transfer case (see e.g. Landi Degl’Innocenti & Landolfi 2004). Since these quantities enter the ESE, their accuracy is crucial for the accuracy of the atomic density matrix elements and, consequently, for the accuracy of the whole non-LTE solution.

In RT problems involving three-dimensional (3D) model atmospheres and scattering polarisation, little is known about the suitability of different quadrature rules. In general, the higher the number of rays is in the quadrature, the better the accuracy of the $J_Q^K$ tensor integration is. On the other hand, the overall computing time (CPU time) increases approximately proportionally to the number of ray directions. This becomes especially important for solving large-scale problems in 3D models resulting from realistic radiation-magnetohydrodynamical simulations where the CPU time for one formal solution can be of the order of tens of thousands of hours (e.g. Štěpán & Trujillo Bueno 2016).

Our goal in this paper is to find an angular quadrature for the solution of the transfer of polarised radiation that guarantees any given accuracy with the minimum number of points, that is, an optimal quadrature. While recently, significant progress has been made in the mathematical research on optimal quadratures in general, information was still lacking as to the specific field of polarised radiative transfer. This paper tries to fill in this gap. Since the quadratures we derived are optimal or near optimal for the problems arising in the theory of polarised radiative transfer, they are generally more efficient than the currently used quadratures. In Sect. 2 we formulate quantitatively the problem of optimal quadrature for spectropolarimetry. In Sect. 3 we derive the key system of equations and we discuss their numerical solution. In Sect. 4 we present our results for various levels of accuracy. Finally, in Sect. 5 we summarise our findings.

2. Formulation of the problem

The integral of a real or complex-valued function f(x) over the unit sphere 𝕊² = {Ω ∈ ℝ³ : ∥Ω∥ = 1} can be approximated by a weighted sum of the function values sampled at a finite number of N discrete points ${\left\{{\mathbf{\Omega }}_i\right\}}_{i=1}^N$ on 𝕊². In case the quadrature is exact for the function f(Ω),

$$\begin{array}{c}\hfill {\displaystyle {\int }_{{\mathbb{S}}^2}f(\mathbf{\Omega })\mathrm{d}\mathbf{\Omega }=\sum _{i=1}^N{w}_{i}f({\mathbf{\Omega }}_i),}\end{array}$$(1)

where ${\left\{w_i\right\}}_{i=1}^N$ are positive weights associated with the corresponding unit vectors Ω_i. The quadrature ${\{w_i,{\mathbf{\Omega }}_i\}}_{i=1}^N$ is said to be optimal if N is the smallest possible number, such that the quadrature integrates exactly a particular set of functions ${\left\{{\phi }_k\right\}}_{k=1}^M$:

$$\begin{array}{c}\hfill {\displaystyle {\int }_{{\mathbb{S}}^2}{\phi }_k(\mathbf{\Omega })\mathrm{d}\mathbf{\Omega }=\sum _{i=1}^N{w}_i{\phi }_k({\mathbf{\Omega }}_i),1\le k\le M.}\end{array}$$(2)

In one dimension, the subject of quadratures is very well understood (Press et al. 2007; Dahlquist & Björck 2008). The N-point Gaussian quadratures are perhaps the most prominent example of the quadratures that integrate the polynomials up to the order 2N − 1 exactly. This quadrature has also been successfully used in numerical RT. On the other hand, integration on the spherical surface, or angular integration, is still a subject of active mathematical research. The most common approach to the problem is to choose ${\left\{{\phi }_k\right\}}_{k=1}^M$ to be a set of spherical harmonics up to a certain order L (e.g. Lebedev 1976; Brauchart & Grabner 2015; Beentjes 2015), but alternative approaches exist (e.g. Ahrens & Beylkin 2009). For more general methods of numerical integration of functions in ℝ^d, see Xiao & Gimbutas (2010).

In polarisation transfer problems, we deal with angular dependencies of the four Stokes parameters, namely the specific intensity I, two components Q and U of the linear polarisation, and the circular polarisation V. These parameters can be grouped into a formal vector and indexed as (I, Q, U, V) ≡ (I₀, I₁, I₂, I₃). It is important to realise that the amplitude of I₀ is often much larger than that of the other parameters, but, in practice, we still need to be able to integrate these low-amplitude quantities to the accuracy of a few decimal places using the same quadrature.

In the case of non-axisymmetric problems, and especially in the general 3D radiative transfer, a common choice of the quadrature over the unit sphere is a product of two one-dimensional quadratures (or tensor product quadrature), namely the Gaussian quadrature in the cosine of inclinations and the equally-spaced trapezoidal quadrature in azimuths (GT quadrature in the following text). While this quadrature is easy to construct, it is not optimal due to unnecessary accumulation of points near the “poles” (see Fig. 1 of Beentjes 2015). Another popular set of quadrature rules in numerical radiative transfer, which is not of the tensor product type, is the one of Carlson (1963). The RT codes allowing 3D radiative transfer, such as RHSC3D (Uitenbroek 1998), Multi3D (Leenaarts & Carlsson 2009), and PORTA (Štěpán & Trujillo Bueno 2013), typically implement one of these quadratures. In this paper, we do not discuss the Carlson quadratures because our calculations indicate that they fail to correctly integrate even weakly anisotropic polarised radiation fields. Instead, we compare our newly derived quadratures with the GT ones.

While in the unpolarised RT transfer theory, the important quantity is the mean radiation field

$$\begin{array}{c}\hfill J={\displaystyle {\int }_{{\mathbb{S}}^2}I(\mathbf{\Omega })\mathrm{d}\mathbf{\Omega },}\end{array}$$(3)

the key quantities in the polarised transfer case are the irreducible spherical tensors of the radiation field

$$\begin{array}{c}\hfill J_Q^K={\displaystyle \sum _{j=0}^3{\int }_{{\mathbb{S}}^2}{I}_j(\mathbf{\Omega }){\mathcal{T}}_Q^K(j,\mathbf{\Omega })\mathrm{d}\mathbf{\Omega },}\end{array}$$(4)

where ${\mathcal{T}}_Q^K(j,\mathbf{\Omega })$ is the geometrical tensor defined in Table 5.6 of Landi Degl’Innocenti & Landolfi (2004) or Table 1 of Bommier (1997).

It is natural to expand the angular dependence of the Stokes parameters into the basis of spherical harmonics and to restrict the basis to a certain maximum order L. We then define the optimal quadrature of the order L as a quadrature with a minimum number of points N, which integrates exactly all of the $J_Q^K$ tensors for the Stokes parameters expansion up to the order L. The suitable value of L depends on the problem to be solved. Our experience with the 3D non-LTE calculations in the solar atmosphere shows that L ≈ 10 should be sufficient for such kind of applications. We leave the discussion of this topic for Sect. 4 and for the following paper in this series.

There are additional constraints on the quadratures, which we require for practical reasons. Firstly, none of the rays are allowed to be parallel to any of the XYZ axes of the reference frame. The formal solution of the RTE and the calculation of the ray intersections with the mesh planes is more efficient if all the rays are not parallel to any axes. Secondly, we impose the natural assumption that the quadrature is rotationally invariant with respect to 90° rotations around the Z axis and that it is mirror-symmetric with respect to the Z = 0 plane. This way, we can restrict the search for the quadrature points to just n = N/8 rays in the first octant and to obtain the rest of the rays by simple rotations and/or reflection. We note that this choice of the quadrature symmetries is not the only possible one and more investigations should be done in the future regarding other constraints.

3. Calculation of the quadratures

We start our derivation by expanding the Stokes parameters ${\{I_j(\mathbf{\Omega })\}}_{j=0}^3$ into the basis of spherical harmonics,

$$\begin{array}{c}\hfill I_j(\mathbf{\Omega })={\displaystyle \sum _{\ell =0}^{\infty }\sum _{m=-\ell }^{\ell }{(I_j)}_{\ell m}{Y}_{\ell m}(\mathbf{\Omega }),}\end{array}$$(5)

where

$$\begin{array}{c}\hfill Y_{\ell m}(\mathbf{\Omega })=\{\begin{array}{cc}i\sqrt{2\pi }[Y_{\ell }^m(\mathbf{\Omega })-{(-1)}^m{Y}_{\ell }^{-m}(\mathbf{\Omega })]\hfill & \mathrm{if}m<0\hfill \\ 2\sqrt{\pi }{Y}_{\ell }^0\hfill & \mathrm{if}m=0\hfill \\ \sqrt{2\pi }[Y_{\ell }^{-m}(\mathbf{\Omega })+{(-1)}^m{Y}_{\ell }^m(\mathbf{\Omega })]\hfill & \mathrm{if}m>0\hfill \end{array}\end{array}$$(6)

are the real-valued spherical harmonics constructed from the complex-valued spherical harmonics

$$\begin{array}{c}\hfill Y_{\ell }^m(\chi ,\theta )=\sqrt{\frac{2\ell +1}{4\pi }\frac{(\ell -m)!}{(\ell +m)!}}{P}_{\ell }^m(cos\theta ){\mathrm{e}}^{im\chi },\end{array}$$(7)

where $P_{\ell }^m(\cos \theta )$ are the associated Legendre polynomials (Abramowitz & Stegun 2014). In this paper, we adopt the common definition in which the inclination angle θ is measured from the positive axis Z and the azimuth χ is measured from the positive axis X. The spherical harmonics form an orthonormal basis for functions on 𝕊². While the functions to be integrated are complex valued, we use the real spherical harmonics for expansion of the real-valued Stokes parameters. The expansion coefficients (I_j)_ℓm can be obtained from the inverse transformation

$$\begin{array}{c}\hfill {(I_j)}_{\ell m}={\displaystyle {\int }_{{\mathbb{S}}^2}{I}_j(\mathbf{\Omega })Y_{\ell m}(\mathbf{\Omega })\mathrm{d}\mathbf{\Omega }.}\end{array}$$(8)

In this work, we assume that the angular dependence of I_j(Ω) is such that decomposition into spherical harmonics can be truncated at a finite order L. With this assumption, Eq. (5) becomes

$$\begin{array}{c}\hfill I_j(\mathbf{\Omega })={\displaystyle \sum _{\ell =0}^L\sum _{m=-\ell }^{\ell }{(I_j)}_{\ell m}{Y}_{\ell m}(\mathbf{\Omega }).}\end{array}$$(9)

This is usually a reasonable assumption in stellar photospheres and chromospheres, such as the solar ones, in which the angular dependence is not particularly peaked at some angles. In other cases of astrophysical interest, it may be desirable to derive different quadrature rules. The value of L depends on the particular problem, the desired accuracy, and the constraints on the CPU time. In general, the larger L is, the higher the accuracy is of the numerical solution.

Using the expansion of I_j(Ω) from Eq. (9) in Eq. (4), we obtain the expression for the radiation field tensor in terms of the Stokes-parameter coefficients

$$\begin{array}{c}\hfill J_Q^K={\displaystyle \sum _{j=0}^3\sum _{\ell =0}^L\sum _{m=-\ell }^{\ell }{(I_j)}_{\ell m}y(j,\ell ,m,K,Q),}\end{array}$$(10)

where the complex-valued factors

$$\begin{array}{c}\hfill y(j,\ell ,m,K,Q)={\displaystyle {\int }_{{\mathbb{S}}^2}{Y}_{\ell m}(\mathbf{\Omega }){\mathcal{T}}_Q^K(j,\mathbf{\Omega })\mathrm{d}\mathbf{\Omega }}\end{array}$$(11)

can be analytically calculated for any combination of the indices. From the definition of the quadrature in Eq. (1) and by applying the expansion of the Stokes parameters, we can express $J_Q^K$ as

$$\begin{array}{cc}\hfill J_Q^K& ={\displaystyle \sum _{j^{\prime }=0}^3\sum _{i=1}^N{w}_i{I}_{j^{\prime }}({\mathbf{\Omega }}_i)}{\mathcal{T}}_Q^K(j^{\prime },{\mathbf{\Omega }}_i)\hfill \\ \hfill & ={\displaystyle \sum _{j^{\prime }=0}^3\sum _{{\ell }^{\prime }=0}^L\sum _{m^{\prime }=-{\ell }^{\prime }}^{{\ell }^{\prime }}{(I_{j^{\prime }})}_{{\ell }^{\prime }{m}^{\prime }}\sum _{i=1}^N{w}_i{Y}_{{\ell }^{\prime }{m}^{\prime }}({\mathbf{\Omega }}_i){\mathcal{T}}_Q^K(j^{\prime },{\mathbf{\Omega }}_i).}\hfill \end{array}$$(12)

Since the right-hand sides of Eqs. (10) and (12) must be equal for any realisation of the radiation field, that is, for any values of (I_j)_ℓm and (I_j′)_ℓ′m′, it follows that the coefficients multiplying these radiation field components must be the same in both equations, that is,

$$\begin{array}{c}\hfill y(j,\ell ,m,K,Q)={\displaystyle \sum _{i=1}^N{w}_i{Y}_{\ell m}({\mathbf{\Omega }}_i){\mathcal{T}}_Q^K(j,{\mathbf{\Omega }}_i)}\end{array}$$(13)

for any set of indices. Any quadrature, which satisfies the set of Eq. (13), integrates exactly the $J_Q^K$ tensors for any realisation of the radiation field. The optimal quadratures up to the order L are the ones with the minimum value of N. The functions $Y_{\ell m}(\mathbf{\Omega }){\mathcal{T}}_Q^K(j,\mathbf{\Omega })$ play the role of the functions φ_k in Eq. (2).

The set of Eq. (13) for the unknown values of w_i and directions Ω_i is non-linear and can only be solved by numerical methods. Moreover, given L, the value of N(L) is also unknown. In practice, knowing the value of n(L − 1) = N(L − 1)/8, we start searching for n(L) = n(L − 1) and, in case no solution is found, we proceed by increasing n(L) by 1 until the quadrature is found. Apart from n unknown weights w_i, the unknown directions are represented by the inclination θ_i ∈ (0° ,90° ) and azimuths χ_i ∈ (0° ,90° ) in the first octant defined by the positive axes X, Y, and Z. Given the necessarily numerical solution of the problem, the exactness of the quadrature is limited by the double-precision computer arithmetics. Our condition for the quadrature to be considered exact is

$$\begin{array}{c}\hfill \underset{j,\ell ,m,K,Q}{max}\parallel y(j,\ell ,m,K,Q)-{\displaystyle \sum _{i=1}^N}{w}_i{Y}_{\ell m}({\mathbf{\Omega }}_i){\mathcal{T}}_Q^K(j,{\mathbf{\Omega }}_i)\parallel <{10}^{-15},\end{array}$$(14)

that is, the maximum error of integration of any of the $Y_{\ell m}(\mathbf{\Omega }){\mathcal{T}}_Q^K(j,\mathbf{\Omega })$ functions is smaller than 10⁻¹⁵. We solved Eq. (13) as a least-squares minimisation problem using the trust region reflective algorithm (TRF, Branch et al. 1999). The initial guess of the quadrature is chosen randomly, but since the problem is non-convex, that is, it suffers from the existence of local minima, it is useful to start with a guess that is as close to the global minimum as possible. Therefore, the initial guess used as a starting point for the TRF method is chosen as the most precise one out of the number of randomly generated quadratures. It can be easily shown that thanks to the symmetry of the quadrature with respect to rotation by integer multiples of 90° around the Z axis, a significant number of Eq. (13) are identically satisfied, which simplifies the numerical solution.

In the case of unpolarised radiative transfer, the problem is significantly simplified. The only non-zero Stokes parameter is the specific intensity I₀ and the only relevant geometrical tensor is ${\mathcal{T}}_0^0(0,\mathbf{\Omega })=1$. Therefore, the Eq. (13) take the following form:

$$\begin{array}{c}\hfill y(0,\ell ,m,0,0)={\displaystyle \sum _{i=1}^N{w}_i{Y}_{\ell m}({\mathbf{\Omega }}_i).}\end{array}$$(15)

The optimal quadrature for unpolarised RT is then simply the one that is optimal for the integration of spherical harmonics.

4. Results

We derived the quadrature rules with an increasing accuracy up to L = 15. The number of these near-optimal quadrature points per octant as a function of L is shown in Fig. 1. For comparison purposes, we also show the number of points per octant in the GT quadrature needed for exact integration up to the same order.

Fig. 1. Comparison of the number of quadrature points per octant, n(L) = N(L)/8, of the GT quadratures (black lines) and the new quadratures (red lines) as a function of the order L. The solid lines correspond to the full-Stokes case and the dashed lines show the unpolarised (Stokes I) quadratures, respectively.

It follows from the figure that with the exception of very small orders, the new quadratures perform better than the GT ones. For some orders, the number of rays is about 30% smaller at the same level of accuracy. A similar behaviour is found for the unpolarised case, which is shown with the dashed lines.

In order to provide some insight as to how to choose L in practical calculations, we define the quantity

$$\begin{array}{c}\hfill {(I_j)}_{\ell }=\frac{1}{2\ell +1}\sqrt{\frac{\sum _{m=-\ell }^{\ell }{(I_j)}_{\ell m}^2}{{(I_0)}_{00}^2}},\end{array}$$(16)

which characterises, at a given point in the 3D model, the expansion coefficients of the given order ℓ for each Stokes parameter. In Fig. 2 we show box plots of (I_j)_ℓ in a particular 3D model of the solar atmosphere. The coefficients decrease approximately exponentially with ℓ; in other words, the role of the coefficients rapidly decreases with their order. An important fact to realise is that the anisotropy of the specific intensity often dominates some of the $J_Q^K$ components, in part, because the intensity is usually significantly higher than the other Stokes parameters. In order to accurately account for the contribution of all of the Stokes parameters, it is important to choose L such that the amplitudes (I₀)_L of the intensity are smaller than at least the first amplitudes (I_j)_ℓ with j >  0 of the other Stokes parameters. In the case of the Ca I 4227 Å line, this means L ⪆ 6.

Fig. 2. Statistical distribution of the amplitudes (Ij)ℓ defined in Eq. (16) of the spherical harmonic components of Stokes I (black) and Stokes Q (red) at the unit optical depth for disc-centre line of sight of the Ca I 4227 Å line in a 3D snapshot model of the solar atmosphere (Carlsson et al. 2016) calculated using a very fine GT quadrature (20 × 20 points per octant).

In Fig. 3 we show how the rays are distributed in the L = 9 and L = 15 quadratures. The rays are more evenly distributed over the unit sphere in contrast to the GT quadrature, hence we avoided oversampling in some regions.

Fig. 3. Distribution of the quadrature points over the unit sphere in the first octant for the case of the L = 9 (left) and L = 15 (right) quadratures. The red areas of the points are proportional to their weights.

In Appendix A, some of the new quadratures with different levels of accuracy are tabulated. All the new quadratures indicated in Fig. 1 can be found in the on-line material.

5. Discussion and conclusions

We derived a set of new quadratures for the angular integration of polarised radiation fields. As a criterion for optimality, we used the condition of exact integration of the $J_Q^K$ radiation field tensors up to a given order of expansion of the Stokes parameters in the basis of spherical harmonics. For the maximum orders up to L = 15 that we studied, our quadratures can save up to 30% of ray directions, hence approximately 30% of computing time compared to the GT quadratures of the same accuracy.

Our derivation of near optimal angular quadratures was based on the $J_Q^K$ tensors that appear in the complete frequency redistribution (CRD) theory of spectral line formation (see Landi Degl’Innocenti & Landolfi 2004). In the general case of partial redistribution (PRD), other quadratures may prove to be more suitable. We also note that in some situations, such as the case of plasma structures embedded in outer stellar envelopes, it may be advantageous to consider different quadratures to properly quantify the impact of the highly anisotropic illumination of the structure from the stellar surface. This type of investigation is out of the scope of the present paper.

Due to the non-linearity and non-convexity of the problem, there is not a unique solution to the problem. Moreover, we cannot prove analytically that the resulting quadratures are really optimal in terms of the minimal number of quadrature points, hence we can only claim that our empirical approach provides nearly optimal quadratures. Our numerical experiments indicate that quadratures with different symmetry properties with respect to rotations and reflections of the coordinate system may be of future interest in some applications. In the following paper, we will demonstrate a practical comparison of the new quadratures with the GT quadratures in realistic 3D non-LTE radiative transfer models of the solar atmosphere.

Acknowledgments

J.Š. acknowledges the financial support of the grant 19-20632S of the Czech Grant Foundation (GAČR) and from project RVO: 67985815 of the Astronomical Institute of the Czech Academy of Sciences. J.T.B. acknowledges the funding received from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (ERC Advanced Grant agreement No. 742265).

References

Appendix A: Tables of quadratures

As an example of the quadratures derived in this paper, here we show three quadratures of increasing accuracy for the full set of Stokes parameters and the $J_Q^K$ tensors (see Table A.1). We also show one quadrature for just the specific intensity (Table A.2)¹. The tables of the other quadratures shown in Fig. 1 are available at the CDS.

The n = N/8 quadrature points in the tables correspond to the points in the first octant. The points in the remaining Z >  0 octants were generated by rotating the first-octant points by 90, 180, and 270°, while preserving their weights. The points in the Z <  0 octants were generated by reflection of all the Z >  0 points with respect to the Z = 0 plane and by preserving their weights.

All Tables

All Figures

	Fig. 1. Comparison of the number of quadrature points per octant, n(L) = N(L)/8, of the GT quadratures (black lines) and the new quadratures (red lines) as a function of the order L. The solid lines correspond to the full-Stokes case and the dashed lines show the unpolarised (Stokes I) quadratures, respectively.
In the text

	Fig. 2. Statistical distribution of the amplitudes (Ij)ℓ defined in Eq. (16) of the spherical harmonic components of Stokes I (black) and Stokes Q (red) at the unit optical depth for disc-centre line of sight of the Ca I 4227 Å line in a 3D snapshot model of the solar atmosphere (Carlsson et al. 2016) calculated using a very fine GT quadrature (20 × 20 points per octant).
In the text

	Fig. 3. Distribution of the quadrature points over the unit sphere in the first octant for the case of the L = 9 (left) and L = 15 (right) quadratures. The red areas of the points are proportional to their weights.
In the text