Issue 
A&A
Volume 636, April 2020



Article Number  A24  
Number of page(s)  5  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/202037566  
Published online  08 April 2020 
Near optimal angular quadratures for polarised radiative transfer^{⋆}
^{1}
Astronomical Institute ASCR, v.v.i., Ondřejov, Czech Republic
email: jiri.stepan@asu.cas.cz
^{2}
Instituto de Astrofísica de Canarias, Vía Láctea s/n, 38205 La Laguna, Tenerife, Spain
email: jjaumeb@iac.es; jtb@iac.es
^{3}
Departamento de Astrofísica, Universidad de La Laguna (ULL), 38206 La Laguna, Tenerife, Spain
^{4}
Consejo Superior de Investigaciones Científicas, Spain
Received:
23
January
2020
Accepted:
27
February
2020
In threedimensional (3D) radiative transfer (RT) problems, the tensor product quadratures are generally not optimal in terms of the number of discrete ray directions needed for a given accuracy of the angular integration of the radiation field. In this paper, we derive a new set of angular quadrature rules that are more suitable for solving 3D RT problems with the short and longcharacteristics formal solvers. These quadratures are more suitable than the currently used ones for the numerical calculation of the radiation field tensors that are relevant in the problem of the generation and transfer of polarised radiation without assuming local thermodynamical equilibrium (nonLTE). We show that our new quadratures can save up to about 30% of computing time with respect to the Gaussiantrapezoidal product quadratures with the same accuracy.
Key words: methods: numerical / polarization / radiative transfer
The tables mentioned in Appendix A are only available at the CDS via anonymous ftp to cdsarc.ustrasbg.fr (130.79.128.5) or via http://cdsarc.ustrasbg.fr/vizbin/cat/J/A+A/636/A24
© ESO 2020
1. Introduction
The nonlinear radiative transfer problem out of local thermodynamic equilibrium (hereafter, nonLTE) 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 nonaxisymmetric and generally multidimensional 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 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 nonLTE solution.
In RT problems involving threedimensional (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 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 largescale problems in 3D models resulting from realistic radiationmagnetohydrodynamical 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 complexvalued function f(x) over the unit sphere 𝕊^{2} = {Ω ∈ ℝ^{3} : ∥Ω∥ = 1} can be approximated by a weighted sum of the function values sampled at a finite number of N discrete points on 𝕊^{2}. In case the quadrature is exact for the function f(Ω),
where are positive weights associated with the corresponding unit vectors Ω_{i}. The quadrature is said to be optimal if N is the smallest possible number, such that the quadrature integrates exactly a particular set of functions :
In one dimension, the subject of quadratures is very well understood (Press et al. 2007; Dahlquist & Björck 2008). The Npoint 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 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_{0}, I_{1}, I_{2}, I_{3}). It is important to realise that the amplitude of I_{0} is often much larger than that of the other parameters, but, in practice, we still need to be able to integrate these lowamplitude quantities to the accuracy of a few decimal places using the same quadrature.
In the case of nonaxisymmetric problems, and especially in the general 3D radiative transfer, a common choice of the quadrature over the unit sphere is a product of two onedimensional quadratures (or tensor product quadrature), namely the Gaussian quadrature in the cosine of inclinations and the equallyspaced 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
the key quantities in the polarised transfer case are the irreducible spherical tensors of the radiation field
where 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 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 nonLTE 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 mirrorsymmetric 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 into the basis of spherical harmonics,
where
are the realvalued spherical harmonics constructed from the complexvalued spherical harmonics
where 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 𝕊^{2}. While the functions to be integrated are complex valued, we use the real spherical harmonics for expansion of the realvalued Stokes parameters. The expansion coefficients (I_{j})_{ℓm} can be obtained from the inverse transformation
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
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 Stokesparameter coefficients
where the complexvalued factors
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 as
Since the righthand 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,
for any set of indices. Any quadrature, which satisfies the set of Eq. (13), integrates exactly the 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 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 nonlinear 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 doubleprecision computer arithmetics. Our condition for the quadrature to be considered exact is
that is, the maximum error of integration of any of the functions is smaller than 10^{−15}. We solved Eq. (13) as a leastsquares 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 nonconvex, 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 nonzero Stokes parameter is the specific intensity I_{0} and the only relevant geometrical tensor is . Therefore, the Eq. (13) take the following form:
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 nearoptimal 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 fullStokes case and the dashed lines show the unpolarised (Stokes I) quadratures, respectively. 

Open with DEXTER 
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
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 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_{0})_{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 (I_{j})_{ℓ} defined in Eq. (16) of the spherical harmonic components of Stokes I (black) and Stokes Q (red) at the unit optical depth for disccentre 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). 

Open with DEXTER 
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. 

Open with DEXTER 
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 online 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 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 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 nonlinearity and nonconvexity 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 nonLTE radiative transfer models of the solar atmosphere.
We note that according to the results of Sect. 3, any quadrature that is capable of exact integration of the spherical harmonics up to the order L can be used for the angular integration of the intensity.
Acknowledgments
J.Š. acknowledges the financial support of the grant 1920632S 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
 Abramowitz, M., & Stegun, I. 2014, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables (Martino Fine Books) [Google Scholar]
 Ahrens, C., & Beylkin, G. 2009, Proc. R. Soc. A, 465, 3103 [NASA ADS] [CrossRef] [Google Scholar]
 Branch, M. A., Coleman, T. F., & Li, Y. 1999, SIAM J. Sci. Comput., 21, 1 [CrossRef] [Google Scholar]
 Beentjes, C. H. L. 2015, Quadrature on a Spherical Surface. Technical Report (Oxford: University of Oxford) [Google Scholar]
 Bommier, V. 1997, A&A, 328, 726 [NASA ADS] [Google Scholar]
 Brauchart, J. S., & Grabner, P. J. 2015, J. Complex., 31, 293 [CrossRef] [Google Scholar]
 Carlson, B. G. 1963, in Methods in Computational Physics, eds. B. Alder, & S. Fernbach, 1, 1 [Google Scholar]
 Carlsson, M., Hansteen, V. H., Gudiksen, B. V., et al. 2016, A&A, 585, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dahlquist, G., & Björck, Å. 2008, Numerical Methods in Scientific Computing (Philadelphia: Society for Industrial and Applied Mathematics), 1 [CrossRef] [Google Scholar]
 Landi Degl’Innocenti, E., & Landolfi, M. 2004, Polarization in Spectral Lines (Dordrecht: Kluwer) [CrossRef] [Google Scholar]
 Lebedev, V. I. 1976, Zh. Vychisl. Mat. Mat. Fiz., 16, 293 [Google Scholar]
 Leenaarts, J., & Carlsson, M. 2009, ASP Conf. Ser., 415, 87 [Google Scholar]
 Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 2007, Numerical Recipes (New York: Cambridge University Press) [Google Scholar]
 Štěpán, J., & Trujillo Bueno, J. 2013, A&A, 557, A143 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Štěpán, J., & Trujillo Bueno, J. 2016, ApJ, 826, L10 [NASA ADS] [CrossRef] [Google Scholar]
 Uitenbroek, H. 1998, ApJ, 498, 427 [NASA ADS] [CrossRef] [Google Scholar]
 Xiao, H., & Gimbutas, Z. 2010, Comput. Math. Appl., 59, 663 [CrossRef] [Google Scholar]
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 tensors (see Table A.1). We also show one quadrature for just the specific intensity (Table A.2)^{1}. The tables of the other quadratures shown in Fig. 1 are available at the CDS.
Fullstokes quadratures.
StokesI quadrature for L = 9, n = 6.
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 firstoctant 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 fullStokes case and the dashed lines show the unpolarised (Stokes I) quadratures, respectively. 

Open with DEXTER  
In the text 
Fig. 2. Statistical distribution of the amplitudes (I_{j})_{ℓ} defined in Eq. (16) of the spherical harmonic components of Stokes I (black) and Stokes Q (red) at the unit optical depth for disccentre 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). 

Open with DEXTER  
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. 

Open with DEXTER  
In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.