Issue 
A&A
Volume 633, January 2020



Article Number  A111  
Number of page(s)  4  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/201937116  
Published online  17 January 2020 
An adaptive Gaussian quadrature for the Voigt function
^{1}
Université de Toulouse, UPSObservatoire MidiPyrénées, CNRS, CNES, IRAP, Toulouse, France
email: frederic.paletou@univtlse3.fr
^{2}
CNRS, Institut de Recherche en Astrophysique et Planétologie, 14 av. E. Belin, 31400 Toulouse, France
^{3}
Université de Toulouse, Faculté des Sciences et de l’Ingénierie, 31062 Toulouse Cedex 9, France
^{4}
CNRS, Cesbio, UPS–CNES–IRD, 18 av. E. Belin, 31401 Toulouse, France
Received:
15
November
2019
Accepted:
10
December
2019
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 nonLTE” radiation transfer computations, assuming potential deviations of the velocity distribution of massive particles from the usual Maxwell–Boltzmann distribution. A first validation is made with computations of the usual Voigt profile.
Key words: methods: numerical / radiation mechanisms: general / line: profiles
© F. Paletou et al. 2020
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1. 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 lightmatter 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 socalled 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).
NonMaxwellian 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 socalled κ VDFs into their photon scattering physical model.
However, such nonMaxwellian VDFs are still known ab initio before solving the radiation transfer problem. The more general issue of computing selfconsistent nonequilibrium distributions for both photons and massive particles – for whose associated problem we coined the phrase “full nonLTE 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. 1986, 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.
2. 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 “IIA” 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 v 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_{IIA} 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 ), 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.
3. The numerical problem
We aim to generalize computations of redistribution functions in order to be able to compute VDFs selfconsistently with the radiation field. Therefore, we need a robust numerical approach to repeatedly perform numerical integrations like:
where x′ and x are the usual incoming and outgoing reduced^{2} frequencies in the observer’s frame, Δν_{D} the Doppler width defined as (ν_{0}/c)v_{th.}, and γ the diffusion angle between incoming and outgoing directions in the plane defined by u_{1} and u_{2}. For the Maxwell–Boltzmann case, we should indeed use:
but we would need to consider f(u_{1}) to be nonanalytic, 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 nonLTE gas of moderate optical thickness. It should also be noted that, for a preliminary study, we assume a selfconsistent 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.
Fig. 1. Failure of a standard Gauss–Hermite quadrature of order k = 70 (green), as compared to the almost superimposed results from, respectively, the method using the Faddeeva complex function (dark), and our alternative double adaptive Gaussian quadrature scheme, for a normalized Voigt profile with a = 10^{−2}. 
4. 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 𝒩 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.
The choice of and is of importance. Liu & Pierce suggested adopting to be the mode of g, and , where:
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.
5. 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 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.
Fig. 2. Example of distribution of nodes for an initial Gauss–Hermite quadrature of order k = 70 (dots), and for our adaptive Gaussian quadrature centered at the Lorentzian peak (crosses). 
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 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 . 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.
Fig. 3. Respective distributions of nodes, marked by crosses of different colors, from an initial Gauss–Hermite quadrature of order k = 70 when the Gaussian and Lorentzian peaks are of comparable amplitude, here for u around three. 
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 finetuning 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.
Fig. 4. Voigt profiles H(a, u) computed with our double AGQ scheme for, respectively, a = 0.001, 10^{−4}, 10^{−5} and 10^{−6} (and decreasing wing values). Small discontinuities are still noticeable at the transition values about two and four. This should not, however, impair any standard scattering integral computation. 
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. 
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.
6. Conclusion
We have tested a suitable numerical strategy for our first step towards “fully nonLTE” 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 nonunimodal 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 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 largeamplitude 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 nonGaussian functions convolved with a Lorentzian may also be doable, using a Fourier transformbased method (e.g., Mendenhall 2007).
It should be noted that there is a “typo” or mistake in this article, more specifically in its Eq. (3). A cos^{2}(α/2) term appears, instead of the correct csc^{2}(α/2), where his α is our γ.
The original article of W. Voigt (1912) published Sitz. Ber. Bayer. Akad. München (in German) can be found online at http://publikationen.badw.de/de/003395768
Acknowledgments
Frédéric Paletou is grateful to his radiative transfer sensei, Dr. L. H. “Larry” Auer, with whom we started discussing about these issues long time ago. We also thank an anonymous referee for his valuable remarks on an earlier version of this article.
References
 Borsenberger, J., Oxenius, J., & Simonneau, E. 1986, J. Quant. Spectr. Rad. Transf., 35, 303 [NASA ADS] [CrossRef] [Google Scholar]
 Borsenberger, J., Oxenius, J., & Simonneau, E. 1987, J. Quant. Spectr. Rad. Transf., 37, 331 [NASA ADS] [CrossRef] [Google Scholar]
 Chaufray, J.Y., & Leblanc, F. 2013, Icarus, 223, 975 [NASA ADS] [CrossRef] [Google Scholar]
 Henyey, L. G. 1940, Proc. Natl. Acad. Sci., 26, 50 [NASA ADS] [CrossRef] [Google Scholar]
 Hubeny, I., & Mihalas, D. 2014, Theory of Stellar Atmospheres (Princeton: Princeton University Press) [Google Scholar]
 Hummer, D. G. 1962, MNRAS, 125, 21 [NASA ADS] [CrossRef] [Google Scholar]
 Jeffrey, N. L. S., Fletcher, L., & Labrosse, N. 2017, ApJ, 836, 35 [NASA ADS] [CrossRef] [Google Scholar]
 Liu, Q., & Pierce, D. A. 1994, Biometrika, 81, 624 [Google Scholar]
 McLean, A. B., Mitchell, C. E. J., & Swanston, D. M. 1994, J. Quant. Spectr. Rad. Transf., 69, 125 [Google Scholar]
 Mendenhall, M. H. 2007, J. Quant. Spectr. Rad. Transf., 105, 519 [NASA ADS] [CrossRef] [Google Scholar]
 Oliphant, T. E. 2006, A Guide to NumPy (USA: Trelgol Publishing) [Google Scholar]
 Oxenius, J. 1986, Kinetic Theory of Particles and Photons  Theoretical Foundations of nonLTE Plasma Spectroscopy (Berlin: Springer) [CrossRef] [Google Scholar]
 Oxenius, J., & Simonneau, E. 1994, Ann. Phys., 234, 60 [NASA ADS] [CrossRef] [Google Scholar]
 Rutily, B., & Chevallier, L. 2006, EAS Pub. Ser., 18, 1 [CrossRef] [Google Scholar]
 Rutten, R. J. 2003, Radiative Transfer in Stellar Atmospheres, Lecture Notes, 8th edn., Utrecht University [Google Scholar]
 Schreier, F. 2018, J. Quant. Spectr. Rad. Transf., 213, 13 [NASA ADS] [CrossRef] [Google Scholar]
 Scudder, J. D. 1992, ApJ, 398, 299 [NASA ADS] [CrossRef] [Google Scholar]
All Figures
Fig. 1. Failure of a standard Gauss–Hermite quadrature of order k = 70 (green), as compared to the almost superimposed results from, respectively, the method using the Faddeeva complex function (dark), and our alternative double adaptive Gaussian quadrature scheme, for a normalized Voigt profile with a = 10^{−2}. 

In the text 
Fig. 2. Example of distribution of nodes for an initial Gauss–Hermite quadrature of order k = 70 (dots), and for our adaptive Gaussian quadrature centered at the Lorentzian peak (crosses). 

In the text 
Fig. 3. Respective distributions of nodes, marked by crosses of different colors, from an initial Gauss–Hermite quadrature of order k = 70 when the Gaussian and Lorentzian peaks are of comparable amplitude, here for u around three. 

In the text 
Fig. 4. Voigt profiles H(a, u) computed with our double AGQ scheme for, respectively, a = 0.001, 10^{−4}, 10^{−5} and 10^{−6} (and decreasing wing values). Small discontinuities are still noticeable at the transition values about two and four. This should not, however, impair any standard scattering integral computation. 

In the text 
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. 

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.