An adaptive Gaussian quadrature for the Voigt function

We evaluate an adaptive Gaussian quadrature integration scheme suitable for the numerical evaluation of generalized redistribution in frequency functions. The latter are indispensable ingredients for “full non-LTE” radiation transfer computations, assuming potential deviations of the velocity distribution of massive particles from the usual Maxwell–Boltzmann distribution. A ﬁrst validation is made with computations of the usual Voigt proﬁle.


Introduction
Radiation transfer is, by essence, a difficult problem (e.g., Rutily & Chevallier 2006), as well as a question of very large relevance in astrophysics. It relies indeed on complex nonlinear light-matter interactions (see e.g., Hubeny & Mihalas 2014;Rutten 2003).
At the very heart of the problem lies the issue of how photons may scatter on these moving massive particles constituting the atmosphere being studied. The usual literature classifies these processes as either, "complete" redistribution in frequency (CRD), or "partial" redistribution in frequency (PRD). Complete redistribution in frequency means complete decorrelation of absorbed and emitted photons, while PRD considers a fraction of coherent scattering, up to the limit of full coherence (see, e.g., Sect. 10 of Hubeny & Mihalas 2014). The vast majority of astrophysical problems are solved, still, within the frame of CRD, for which further simplifications are the equality of emission and absorption profiles -the latter also usually known a priori -which also leads to the independence of the so-called source function vs. frequency.
Besides, and more generally, while nonequilibrium distributions of photons, meaning potential departures from the Planck law, have been routinely considered since the late 60's, a very limited number of studies tried to further push the description of the physical problem, by questioning to what extent the most often assumed Maxwell-Boltzmann velocity distribution of the massive particles onto which photons scatter may remain valid? (see e.g., Oxenius & Simonneau 1994, and references therein).
Non-Maxwellian velocity distribution functions (hereafter, VDF) have been studied (e.g., Scudder 1992, and further citations) or evidenced in natural plasma (see e.g., Jeffrey et al. 2017, for a recent study about solar flares). Such departures from Maxwell-Boltzmann VDFs have also been considered in the radiative modeling of spectral lines formed in neutral planetary exospheres (e.g., Chaufray & Leblanc 2013), where these authors introduced so-called κ VDFs into their photon scattering physical model. However, such non-Maxwellian VDFs are still known ab initio before solving the radiation transfer problem. The more general issue of computing self-consistent nonequilibrium distributions for both photons and massive particles -for whose associated problem we coined the phrase "full non-LTE radiation transfer" -remains quite an open question in astrophysics, although a few studies have already been conducted in the past (see e.g., Borsenberger et al. 1986Borsenberger et al. , 1987. Hereafter, we provide a first numerical tool that allows us to go further in this direction, enabling further computations of generalized redistribution functions. Moreover, the numerical scheme we evaluated may also be of more general interest, for other topics of numerical (astro)physics.

Redistribution in frequency
As an illustrative but important example, we focus here on the case of coherent scattering in the atomic frame of reference, for a spectral line of central wavelength ν 0 . We also assume that only "natural" broadening is at play for the upper energy level of, typically, a resonance line with an infinitely sharp lower level. Therefore, we consider an elementary frequency redistribution function r(ξ , ξ) such that: where ξ and ξ are, respectively, the incoming and the outgoing frequencies of a photon, and δ is the usual Dirac distribution, together with: The latter is a Lorentzian profile, with damping parameter Γ, resulting from the "natural width" of the upper atomic state of the transition at ν 0 . If we assume that the angular redistribution associated with the scattering event is isotropic, such a case of radiation damping and coherence in the atom's frame refers to the standard case "II-A" in the nomenclature of Hummer (1962; see also Hubeny & Mihalas 2014).
Once the elementary scaterring processes have been defined in the atomic frame of reference, we have to consider, for a further practical implementation into a radiative transfer problem, the collective effects induced by the agitation of a pool of massive particles populating the atmosphere. This is precisely in this "jump" to the observer's frame of reference, because of Doppler shifts such as: where ν is the observed frequency, n may be either the incoming or the outgoing direction of a photon, and u the velocity of the massive particle onto which the scattering takes place, that some assumption has to be made about the VDF of the massive atoms (or molecules) present in an atmosphere, under given physical conditions. Detailed derivations of R II-A can be found in the classical literature about redistribution functions, from Henyey (1940) 1 to Hummer (1962). Standard redistribution functions have been first derived assuming that the VDF of the atoms scattering light is a Maxwell-Boltzmann distribution. Then, but more generally, any macroscopic redistribution function in the observer's frame suitable for implementation into the numerical radiative transfer problem would result from the further integration along each velocity components u i (hereafter normalized to the most probable velocity v th. = √ 2kT/m), characterizing the movement of the scattering atoms, and therefore considering these changes of frequencies due to the associated Doppler shifts as expressed by Eq. (3). The latter phenomenon is usually refered to as Doppler, or thermal, broadening.

The numerical problem
We aim to generalize computations of redistribution functions in order to be able to compute VDFs self-consistently with the radiation field. Therefore, we need a robust numerical approach to repeatedly perform numerical integrations like: but we would need to consider f (u 1 ) to be non-analytic, and, at first, (slightly) departing from the Maxwellian standard VDF. Indeed, physical conditions leading to small departures from a Gaussian VDF were already identified and discussed by Oxenius (1986), and they would correspond to a non-LTE gas of moderate optical thickness. It should also be noted that, for a preliminary study, we assume a self-consistent VDF solution of the problem that may still be decomposed as f 1 (u 1 ) f 2 (u 2 ) f 3 (u 3 ). However, before exploring potential departures from Gaussianity, we need to adopt a robust enough numerical strategy in order to numerically evaluate integrals such as Eq. (4), a task which is notoriously difficult even with Maxwell-Boltzmann VDFs. It is very easy to verify that, for instance, a standard Gauss-Hermite (GH) quadrature, even at high rank k, fails at properly computing a somewhat simpler expression like the Voigt 3 function given in Eq. (13). We display, in Fig. 1, the comparison between a GH integration and the new numerical scheme that is presented hereafter.

Adaptive Gaussian quadrature
We start by following the scheme proposed by Liu & Pierce (1994), which is based on the classical GH quadrature. The latter is indeed suitable for integrations such as: Then, the GH quadrature is: where the nodes y i are the zeros of the kth order Hermite polynomial, and w i the corresponding weights. Tabulated values of both nodes and weights can be found very easily, and they are also available for various programming language. We use numpy's (Oliphant 2006) function polynomial.hermite.hermgauss, and a GH of order k = 70 for all results presented hereafter. The main drawback of such a standard quadrature is that function f is scanned at the very nodes y i irrespective of the range where it may have its most significant variations.
However, Liu & Pierce (1994) proposed that, should a function g be integrated, one may define:   where N is the usual Gaussian function: so that one can write: and, finally: This adaptive Gaussian quadrature scheme (AGQ) makes it possible to use the original nodes and weights of the GH quadrature, but somewhat zooms in on these domains where function g has its most significant variations.
We come back to this choice in the following section, and show that a somewhat largerσ value is more suitable for the special case of the Voigt profile.

AGQ tests with the Voigt function
We now consider the normalized Voigt function hereafter defined as: and which satisfies: We highlight that several authors use H for the Voigt profile normalized to √ π, but U instead of our H normalized to unity (see e.g., Hubeny & Mihalas 2014, their Sect. 8). We also use h(y; a, u) for the integrand of Eq. (13). For this numerical integration, three main regimes should be considered, depending on the values of u, according to the respective amplitudes of the Gaussian and the Lorentzian components of the integrand. For the line core region such as |u| < 2, we use a slightly modified AGQ, for which we use a value of σ larger than the one suggested in the original article of Liu & Pierce (1994), given that, for the Voigt function, j = 2(1 + 1/a 2 ). We display, in Fig. 2, the new quadrature nodes, marked with crosses, centered at the Lorentzian peak located at y ≈ 1.7 and using 3σ instead of the value suggested in the original prescription of Liu & Pierce (1994). The nodes of the standard GH quadrature (at the same order) are displayed as dots. They extend too far away and clearly "miss" the large amplitude Lorentzian peak, and therefore the dominant contribution to the integral.
Secondly, we perform a double AGQ scheme for the near wing regions such as 2 < |u| < 4, and for which two discernable peaks of comparable amplitudes result from, respectively, the Lorentzian and the Gaussian components of the convolution (we hereafter refer to u 2 and u 4 for these two boundary values). In such a case, we use both the centering and integration range controls provided by the original AGQ for evaluating the contribution from each component of the integrand separately. For the Lorentzian component, we therefore do the same as when |u| < 2, but we add to this part of the integral the contribution of the nearby Gaussian peak using another AGQ, centered at 0, and of specificσ G adapted to the width of the known Gaussian component of the integrand. More generally, it should be evaluated using criteria allowing a significant covering of the VDF component of the integrand; a relevant ∆y may be such that, at this value, the VDF contribution is already down to a few percent of its maximum amplitude. This also depends, according to Eq. (11), on the maximum value of the nodes coming along with the GH quadrature, at a given order. Finally, overlap with the nearby Lorentzian component should be avoided. For our example, for which the larger node y i 11 (at order k = 70), we usedσ G = 1/8. With such a value, the Gaussian component of the VDF was reasonably well covered, without overlap with the nearby Lorentzian contribution. In Fig. 3, the two distinct sets of nodes, based on the same original GH quadrature nodes, at same order, are displayed by different colored crosses. Finally, for the far wing, where |u| > u 4 , and when the Lorentzian peaks fade out, the usual Gauss-Hermite quadrature is satisfactory.
Results obtained using our double AGQ quadrature scheme are displayed in Fig. 4, for different values of a ranging from 0.001 to 10 −6 , that is more likely regimes expected for our next computations. Maximum relative error computations using the Faddeeva function method as a benchmark, and the scipy.special.wofz Python function, are at most of a few percent, as displayed in Fig. 5. It should also be noted that the latter was obtained using a u 4 value of 4.25, instead of the fiducial value of four, also indicating in which direction a further fine-tuning could be worked out, if necessary, by considering u 2 and u 4 as slowly varying functions of a. Sometimes, we can still notice small discontinuities at the changes of regimes, at u 2 and u 4 . We believe, however, that should our procedure be used for Voigt profile computations and radiative modeling, such small and very local discontinuities will not impair further computations of these scattering integrals entering the equations of the statistical equilibrium.
This new numerical scheme is particularly efficient for small values of a, typically lower than 0.01, where other schemes may fail (see for instance the discussion in Schreier 2018 about the implementation of Voigt1D in the astropy package in Python, using the McLean et al. 1994 method). But, first of all, it is certainly suitable for our next applications of such a numerical integration scheme, and for physical conditions leading to very sharp Lorentzian peaks. We could also test the sensitivity of our scheme to the order of the initial Gauss-Hermite quadrature. For instance, for a < 0.01, we could go down to orders 40-50 without any significant loss of accuracy.
For larger values of a, typically more than 0.1, we noticed that no intermediate scheme between the original Liu & Pierce (1994) at line core, and the Gauss-Hermite in the wings appears necessary. However, the transition value between the two regimes should be adapted to the value of a, in a two to four range.

Conclusion
We have tested a suitable numerical strategy for our first step towards "fully non-LTE" radiative transfer calculations, and the computation of generalized frequency redistribution functions. We modified the original strategy of Liu & Pierce (1994), but also applied it to a non-unimodal distribution.
Our numerical scheme does not pretend to compete with these numerical methods implemented for the very accurate computation of the Voigt function (see, e.g., Schreier 2018, and references therein) since our aim lies elsewhere, such as in exploring departures from Gaussian VDFs. It is, however, providing very good results as compared to reference computations, such as the one using the Faddeeva complex function. Relative Fig. 5. Relative error between our computations with the double AGQ method, for a = 10 −4 , and a reference computation using the Faddeeva complex function. errors down to a few percent are systematically reported in the near wings region, and we believe that further fine tuning could be achieved to reach an even better accuracy. This is, however, not the scope of our study, which aims to compute generalized redistribution functions, after selfconsistent computations of respective distributions of both massive particles and photons under various physical conditions. In that respect, our main concern is for a proper "capture" of the expected very sharp, and therefore very large-amplitude Lorentzian peaks. And we believe that the principle of our numerical integration scheme should remain valid for the more easy to track contribution from the velocity distribution function, even for computed perturbations from a Gaussian shape.
As a final remark, we are also aware that computations with non-Gaussian functions convolved with a Lorentzian may also be doable, using a Fourier transform-based method (e.g., Mendenhall 2007).