EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 633 (January 2020) A&A, 633 (2020) A111 Full HTML
Open Access
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/0004-6361/201937116
Published online 17 January 2020

© F. Paletou et al. 2020

Licence Creative CommonsOpen 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 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. 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:

r ( ξ ′ , ξ ) = φ ( ξ ′ ) δ ( ξ ′ − ξ ) , $$ \begin{aligned} r(\xi \prime ,\xi ) = \varphi (\xi \prime ) \delta (\xi \prime -\xi ), \end{aligned} $$(1)

where ξ′ and ξ are, respectively, the incoming and the outgoing frequencies of a photon, and δ is the usual Dirac distribution, together with:

φ ( ξ ′ ) = ( Γ π ) 1 ( ξ ′ − ν 0 ) 2 + Γ 2 · $$ \begin{aligned} \varphi (\xi \prime ) = \left( \frac{\Gamma }{\pi } \right) {\frac{1}{(\xi \prime - \nu _0)^2 + {\Gamma }^{2}}}\cdot \end{aligned} $$(2)

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:

ν = ξ + ν 0 c n . v , $$ \begin{aligned} \nu = \xi + \frac{\nu _0}{c} {\boldsymbol{n} . \boldsymbol{v}}, \end{aligned} $$(3)

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 RII-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 ui (hereafter normalized to the most probable velocity v th . = 2 k T / m $ \mathit{v}_{\mathrm{th.}}=\sqrt{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.

3. 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:

H 1 ( x ′ , x , γ ) = ∫ − ∞ + ∞ f ( u 1 ) d u 1 [ ( x + x ′ 2 ) sec ( γ / 2 ) − u 1 ] 2 + [ a Δ ν D sec ( γ / 2 ) ] 2 , $$ \begin{aligned} H_{1}(x\prime ,x,\gamma ) = \displaystyle \int _{-\infty }^{+\infty } { { {f(u_1)\mathrm{d}u_1} \over { \left[ (\frac{x+x\prime }{2}) \sec (\gamma /2) - u_1 \right]^2 + \left[ \frac{a}{{\Delta \nu _{\rm D}}} \sec (\gamma /2)\right]^2 } } }, \end{aligned} $$(4)

where x′ and x are the usual incoming and outgoing reduced2 frequencies in the observer’s frame, ΔνD the Doppler width defined as (ν0/c)vth., and γ the diffusion angle between incoming and outgoing directions in the plane defined by u1 and u2. For the Maxwell–Boltzmann case, we should indeed use:

f ( u 1 ) = 1 π e − u 1 2 , $$ \begin{aligned} f(u_1)=\frac{1}{\sqrt{\pi }} \mathrm{e}^{-u_1^2}, \end{aligned} $$(5)

but we would need to consider f(u1) 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 f1(u1)f2(u2)f3(u3).

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 Voigt3 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.

thumbnail 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:

I = ∫ − ∞ + ∞ f ( y ) e − y 2 d y . $$ \begin{aligned} I=\displaystyle \int _{-\infty }^{+\infty } { f({ y}){ {\mathrm{e}^{-{ y}^2}\mathrm{d}{ y}} } } . \end{aligned} $$(6)

Then, the GH quadrature is:

∫ − ∞ + ∞ f ( y ) e − y 2 d y ≃ ∑ i = 1 k w i f ( y i ) , $$ \begin{aligned} \displaystyle \int _{-\infty }^{+\infty } { f({ y}){ {\mathrm{e}^{-{ y}^2}\mathrm{d}{ y}} } } \simeq \sum _{i=1}^{k} {{ w}_i f({ y}_i),} \end{aligned} $$(7)

where the nodes yi are the zeros of the kth order Hermite polynomial, and wi 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 yi 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:

h ( y ) = g ( y ) N ( y ; μ ̂ , σ ̂ ) , $$ \begin{aligned} h({ y})= { {{g}({ y})} \over {\mathcal{N}({ y};\hat{\mu },\hat{\sigma })} } , \end{aligned} $$(8)

where 𝒩 is the usual Gaussian function:

N ( y ; μ ̂ , σ ̂ ) = 1 σ ̂ 2 π e − 1 2 ( y − μ ̂ σ ̂ ) 2 , $$ \begin{aligned} \mathcal{N}({ y};\hat{\mu },\hat{\sigma }) = \frac{1}{\hat{\sigma } \sqrt{2\pi }} \mathrm{e}^{-\frac{1}{2} (\frac{{ y}-\hat{\mu }}{\hat{\sigma }})^2 } , \end{aligned} $$(9)

so that one can write:

∫ − ∞ + ∞ g ( y ) d y = ∫ − ∞ + ∞ h ( y ) N ( y ; μ ̂ , σ ̂ ) d y , $$ \begin{aligned} \displaystyle \int _{-\infty }^{+\infty } { {g}({ y}) \mathrm{d}{ y} } = \displaystyle \int _{-\infty }^{+\infty } { h({ y}){\mathcal{N}({ y};\hat{\mu },\hat{\sigma })} \mathrm{d}{ y} } , \end{aligned} $$(10)

and, finally:

∫ − ∞ + ∞ g ( y ) d y ≃ ∑ i = 1 k w i π h ( μ ̂ + 2 σ ̂ y i ) . $$ \begin{aligned} \displaystyle \int _{-\infty }^{+\infty } { g({ y}) \mathrm{d}{ y} } \simeq \sum _{i=1}^{k} {{{{ w}_i} \over {\sqrt{\pi }}} h(\hat{\mu }+ \sqrt{2} \hat{\sigma } { y}_i)} . \end{aligned} $$(11)

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 μ ̂ $ \hat{\mu} $ and σ ̂ $ \hat{\sigma} $ is of importance. Liu & Pierce suggested adopting μ ̂ $ \hat{\mu} $ to be the mode of g, and σ ̂ = 1 / j $ \hat{\sigma}=1/\sqrt{j} $, where:

j = − ∂ 2 ∂ y 2 log g ( y = μ ̂ ) . $$ \begin{aligned} j = - \frac{\partial ^2}{\partial { y}^2} \log g({ y}=\hat{\mu }) . \end{aligned} $$(12)

We come back to this choice in the following section, and show that a somewhat larger σ ̂ $ \hat{\sigma} $ 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:

H ( a , u ) = 1 π ( a π ) ∫ − ∞ + ∞ e − y 2 d y ( u − y ) 2 + a 2 , $$ \begin{aligned} H(a,u) = \frac{1}{\sqrt{\pi }} \left( \frac{a}{\pi } \right) \int _{-\infty }^{+\infty } { \mathrm{e}^{-{ y}^2} \mathrm{d}{ y} \over {(u-{ y})^2 + a^2} } , \end{aligned} $$(13)

and which satisfies:

∫ − ∞ + ∞ H ( a , u ) d u = 1 . $$ \begin{aligned} \int _{-\infty }^{+\infty } {H(a,u) \mathrm{d}u} = 1 . \end{aligned} $$(14)

We highlight that several authors use H for the Voigt profile normalized to π $ \sqrt{\pi} $, 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 σ ̂ $ \hat{\sigma} $ larger than the one suggested in the original article of Liu & Pierce (1994), given that, for the Voigt function, j = 2(1 + 1/a2). 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 σ ̂ $ 3\hat{\sigma} $ 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.

thumbnail 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 u2 and u4 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 $ \hat{\sigma}_{\mathrm{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 yi ≃ 11 (at order k = 70), we used σ ̂ G = 1 / 8 $ \hat{\sigma}_{\mathrm{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| > u4, and when the Lorentzian peaks fade out, the usual Gauss–Hermite quadrature is satisfactory.

thumbnail 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 u4 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 u2 and u4 as slowly varying functions of a. Sometimes, we can still notice small discontinuities at the changes of regimes, at u2 and u4. 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.

thumbnail 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.

thumbnail 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 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 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 self-consistent 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).

1

It should be noted that there is a “typo” or mistake in this article, more specifically in its Eq. (3). A cos2(α/2) term appears, instead of the correct csc2(α/2), where his α is our γ.

2

The usual reduced frequency x is the difference between frequency ν in the observer’s frame and the central frequency ν0 of a spectral line, divided by the Doppler width ΔνD.

3

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

  1. Borsenberger, J., Oxenius, J., & Simonneau, E. 1986, J. Quant. Spectr. Rad. Transf., 35, 303 [NASA ADS] [CrossRef] [Google Scholar]
  2. Borsenberger, J., Oxenius, J., & Simonneau, E. 1987, J. Quant. Spectr. Rad. Transf., 37, 331 [NASA ADS] [CrossRef] [Google Scholar]
  3. Chaufray, J.-Y., & Leblanc, F. 2013, Icarus, 223, 975 [NASA ADS] [CrossRef] [Google Scholar]
  4. Henyey, L. G. 1940, Proc. Natl. Acad. Sci., 26, 50 [NASA ADS] [CrossRef] [Google Scholar]
  5. Hubeny, I., & Mihalas, D. 2014, Theory of Stellar Atmospheres (Princeton: Princeton University Press) [Google Scholar]
  6. Hummer, D. G. 1962, MNRAS, 125, 21 [NASA ADS] [CrossRef] [Google Scholar]
  7. Jeffrey, N. L. S., Fletcher, L., & Labrosse, N. 2017, ApJ, 836, 35 [NASA ADS] [CrossRef] [Google Scholar]
  8. Liu, Q., & Pierce, D. A. 1994, Biometrika, 81, 624 [Google Scholar]
  9. McLean, A. B., Mitchell, C. E. J., & Swanston, D. M. 1994, J. Quant. Spectr. Rad. Transf., 69, 125 [Google Scholar]
  10. Mendenhall, M. H. 2007, J. Quant. Spectr. Rad. Transf., 105, 519 [NASA ADS] [CrossRef] [Google Scholar]
  11. Oliphant, T. E. 2006, A Guide to NumPy (USA: Trelgol Publishing) [Google Scholar]
  12. Oxenius, J. 1986, Kinetic Theory of Particles and Photons - Theoretical Foundations of non-LTE Plasma Spectroscopy (Berlin: Springer) [CrossRef] [Google Scholar]
  13. Oxenius, J., & Simonneau, E. 1994, Ann. Phys., 234, 60 [NASA ADS] [CrossRef] [Google Scholar]
  14. Rutily, B., & Chevallier, L. 2006, EAS Pub. Ser., 18, 1 [CrossRef] [Google Scholar]
  15. Rutten, R. J. 2003, Radiative Transfer in Stellar Atmospheres, Lecture Notes, 8th edn., Utrecht University [Google Scholar]
  16. Schreier, F. 2018, J. Quant. Spectr. Rad. Transf., 213, 13 [NASA ADS] [CrossRef] [Google Scholar]
  17. Scudder, J. D. 1992, ApJ, 398, 299 [NASA ADS] [CrossRef] [Google Scholar]

All Figures

thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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 (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.