Notes: An assessment of some closed-form expressions for the Voigt function III: Combinations of the Lorentz and Gauss functions

TL;DR

This paper reviews various pseudo-Voigt functions, which are linear combinations of Lorentz and Gauss functions, assessing their accuracy as approximations for the Voigt function through comparison with reference values.

Contribution

It provides a comprehensive review and comparison of multiple closed-form pseudo-Voigt approximations for the Voigt function, highlighting their accuracy limitations.

Findings

Relative errors can reach several percent

Some approximations show significant deviations from reference values

The review clarifies the accuracy trade-offs of different formulas

Abstract

A variety of "pseudo-Voigt" functions, i.e. a linear combination of the Lorentz and Gauss function (occasionally augmented with a correction term), have been proposed as a closed-form approximation for the convolution of the Lorentz and Gauss function known as the Voigt function. First, a compact review of several approximations using a consistent notation is presented. The comparison with accurate reference values indicates relative errors as large as some percent.

Tables

Table 1. Table 1: Execution time (in μ ​ s 𝜇 s \rm\mu s ) measured by the %timeit function in the IPython interpreter. For all test 0 ≤ x ≤ 100 0 𝑥 100 0\leq x\leq 100 with n x = 10001 subscript 𝑛 𝑥 10001 n_{x}=10001 grid points. The tests have been performed on a desktop with an Intel x86_64 CPU “i7-4770” running at 3.4 GHz with 8192 KB cache size.

$y$	10.0	1.0	0.001
Kielkopf	459	424	393
Matveev	451	437	392
hum1wei24	252	326	332
weideman24	803	801	809
zpf16h	693	694	695

Equations

$$K(x,y)=\frac{y}{\pi }{\int }_{-\infty }^{\infty }\frac{{\mathrm{e}}^{-t^2}}{(x-t{)}^2+y^2}\mathrm{d}t$$

$$\displaystyle K(0,y)~{}=~{}\exp(y^{2})\bigl{(}1-\mathop{\rm erf}(y)\bigr{)}~{}=~{}\exp(y^{2})\mathop{\rm erfc}(y)~{}.$$

$$w(z)$$

$$=\exp (-z^2)\mathrm{erfc}(-\mathrm{i}z)withz=x+\mathrm{i}y.$$

$$K(x,y)~{}=~{}K(0,y)\left[\bigr{(}1-\eta(y)\bigl{)}G(x)~{}+~{}\eta(y)L(x)\right]\\$$

$$L(x)$$

$$G(x)$$

$$\eta$$

$$x_{\text{h}}$$

$$\displaystyle C(x,y)=0.016\,\eta(1-\eta)\,\Biggl{[}$$

$$\displaystyle-{10\over 10+(x/x_{\text{h}})^{(9/4)}}\Biggr{]}$$

$$K(x,y)=\frac{\sqrt{\ln 2}}{x_{\text{h}}}\left[(1-\eta )G(x)+\frac{\eta }{\sqrt{\pi \ln 2}}L(x)\right]$$

$$C(x,y)=$$

$$\displaystyle~{}~{}~{}\times\Biggl{[}0.066\,\exp{\left(-0.4\left({x\over x_{\text{h}}}\right)^{2}\right)}$$

$$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}-{1\over 40-5.5\,\left({x\over x_{\text{h}}}\right)^{2}+\left(\frac{x}{x_{\text{h}}}\right)^{4}}\Biggr{]}~{}.$$

$$x_{\text{h}}$$

$$+0.05y\left(1-\frac{2y}{y+\sqrt{y^2+4\ln 2}}\right)$$

$$=x_{\text{h}}^{\text{W}}+0.05y\left(1-\frac{y}{x_{\text{h}}^{\text{W}}}\right)$$

$$\eta$$

$$x_{\text{h}}$$

$$C(x,y)~{}=~{}\eta(1-\eta)\bigl{(}G(x)-L(x)\bigr{)}{k_{1}+k_{2}x^{2}\over 1+k_{3}x^{2}+k_{4}x^{4}}~{}.$$

$$K(x,y)=\frac{\sqrt{\ln 2}}{x_{\text{h}}}\left[(1-\eta )G(x)+\frac{\eta }{\sqrt{\pi \ln 2}}L(x)\right]$$

$$x_{\text{h}}$$

$$\eta$$

$$c_{\text{L}}$$

$$c_{\text{G}}$$

$$K(x,y)=\frac{c_{\text{L}}}{x_{\text{h}}\sqrt{\pi }}L(x)+c_{\text{G}}\frac{\sqrt{\ln 2}}{x_{\text{h}}}G(x)$$

$$x_{\text{h}}(y)$$

$$\beta$$

$$g_{\text{L}}(\nu -\hat{\nu },{\gamma }_{\text{L}})=$$

$$g_{\text{G}}(\nu -\hat{\nu },{\gamma }_{\text{G}})=$$

$$g_{\text{V}}(\nu -\hat{\nu },{\gamma }_{\text{L}},{\gamma }_{\text{G}})=$$

$$x=\sqrt{\ln 2}\frac{\nu -\hat{\nu }}{{\gamma }_{\text{G}}}andy=\sqrt{\ln 2}\frac{{\gamma }_{\text{L}}}{{\gamma }_{\text{G}}}.$$

$$\frac{{\gamma }_{\text{L}}}{{\gamma }_{\text{V}}}$$

$$\frac{{\gamma }_{\text{G}}}{{\gamma }_{\text{V}}}$$

$$k_1$$

$$k_3$$

$$t_0$$

Full text

02.60.-x Numerical approximation and analysis

02.70.-c Computational techniques

††journal: J. Quant. Spectros. & Rad. Transfer, received December 2018, accepted January 2019; doi: 10.1016/j.jqsrt.2019.01.017

1 Introduction

Rapid yet accurate computation of the Voigt function [1] is a challenge in many fields of the physical sciences. Because the convolution integral of a Lorentz and Gauss function does not have an analytical solution, approximations have been discussed in numerous papers. Whereas many modern “state-of-the art” algorithms evaluate the closely related complex error function (also known as complex probability function, plasma dispersion function, or Fadde(y)eva function, cf. e.g., 2, 3, 4) utilizing sophisticated numerical techniques, “simple” closed-form expressions still appear to be attractive.

Linear combinations of the Lorentz and/or Gauss functions have been suggested by several authors. Flores-Llamas et al. [5] proposed a sum of the Lorentzian and its derivatives and Melcher and Gerth [6] fitted the Voigt functions by “generalized Lorentz functions”. McLean et al. [7] have further developed an approximation originally suggested by Martin and Puerta [8, 9], Puerta and Martin [10, 11] and proposed a superposition of four Lorentzians [for a recent assessment see 12]. Linear combinations of a Lorentz and Gauss function (sometimes called “pseudo-Voigt” function and occasionally including a correction term) have been suggested or used by Whiting [13], Matveev [14], Peyre and Principi [15], Kielkopf [16], Wertheim et al. [17], Thompson et al. [18], Teodorescu et al. [19], Ida et al. [20] and Liu et al. [21].

In this note we present an assessment of closed-form expressions for the Voigt function using a combination of Lorentz and Gauss functions. After a short review of the basic definitions in the next subsection, we describe several combinations using a consistent notation (in chronological order). The results of our tests using an accurate Voigt function code as reference are presented in Section 3. The codes have been implemented in Python and a Scientific Python (scipy.org) implementation has been used as a reference. In the final Section 4 we provide a summary and some conclusions.

2 Theory

2.1 The Voigt function

The Voigt function (normalized to $\sqrt{\pi}$) is defined by

[TABLE]

where $x$ is a measure of the distance to the center peak, and $y$ is essentially the ratio of the Lorentzian and Gaussian width, $y=\sqrt{\ln 2}\gamma_{\text{L}}/\gamma_{\text{G}}$. At the line center $x=0$ the Voigt function can be expressed as the exponentially scaled complementary error function

[TABLE]

The Voigt function is symmetric, i.e. $K(-x,y)=K(x,y)$ and essentially reduces to the Lorentz function for large $|x+\mathrm{i}y|$. The Voigt function is identical to the real part of the complex function

[TABLE]

2.2 The Whiting [13] approximation

To our knowledge the first approximation of the Voigt function using a combination of Lorentz and Gauss functions

[TABLE]

with

[TABLE]

is due to Whiting [13]. The weight factor is given by the ratio of the widths of the Lorentz and Voigt profiles or

[TABLE]

An improved approximation is obtained by adding a correction term (i.e. $K(x,y)\longrightarrow K(x,y)+C(x,y)$) that is given by “kind of Lorentzians and Gaussians”,

[TABLE]

By definition, the approximation (5) is exact in the center at $x=0$. Furthermore, both approximations are exact for the limiting cases of pure Lorentz and Gauss functions. According to the author, “this approximation matches the Voigt profile within 5 per cent at worst and is generally within 3 per cent or less.”

2.3 The Matveev [14] approximation

The approximation is given by

[TABLE]

with a correction term

[TABLE]

The weight $\eta$ is defined as in (8) and the half width is given by a refinement of Whiting’s approximation (9)111Note that a prefactor $\gamma_{\text{L}}$ for the correction term is missing in Matveev’s Eq. (5) and is correctly inserted in Titov and Haus [22, after Eq. (12)]. Furthermore, the factor $\eta$ in (12) is not given in Titov and Haus [22, Eq. (12c)].

[TABLE]

Without correction term (12) Matveev [14] reports a “greatest error of $~{}\approx 25\%$ at $\eta=0.1$ and $x\approx 3$”. With correction the maximum error at the line center does not exceed $0.6\%$, and for $x/x_{\text{h}}>6$ the error lies within the limits of $1\%$” for any $y$. In the intermediate frequency regime “the error nowhere exceeds $3\%$.”

2.4 The Kielkopf [16] approximation

Without correction term this approximation is identical to (5), but with the weight and half width defined as

[TABLE]

The correction contains the difference of a Gaussian and Lorentzian multiplied with a rational function of $x$ (see Online B for numerical values)

[TABLE]

This approximation “is accurate to the order of 0.0001 of the peak value of the function”.

2.5 The Thompson et al. [18] approximation

In contrast to the three approximations discussed so far the weight factors of the Thompson et al. [18] (and Liu et al. [21], next subsection) approaches are defined by power series of $y$. The pseudo-Voigt function is written as

[TABLE]

with the Voigt half width and weight (see Online B)

[TABLE]

According to Ida et al. [20] the maximum deviation of about $1.2\%$ is found at $y/(y+\sqrt{\ln 2})\approx 0.5$.

2.6 The Liu et al. [21] approximation

Introducing a dimensionless parameter $d=(y-\sqrt{\ln 2})/(y+\sqrt{\ln 2})$ and approximating the weights (see Eq. (34) in the Online Appendix) as

[TABLE]

the pseudo-Voigt function is written as

[TABLE]

with the half width given by the Olivero and Longbothum [23] approximation

[TABLE]

According to the abstract [21] “the maximum errors of width, area, and peak … are 0.01%, 0.2%, and 0.55%, respectively.”

3 Results

Ignoring the correction terms Eqs. (10), (12), and (17), all algorithms use a weighted sum of the Lorentz and Gauss function, where the Lorentz weight tends to one for $y\gg 1$ and zero for $y\ll 1$. Fig. 1 compares the weights of the Lorentz function for all five approximations. The weights of Whiting and Matveev are identical by construction, Eq. (8), and the Thompson et al. and Liu et al. weights appear to be very similar (with differences for small and large $y$ only). For $y\approx 1$ (i.e. equal width of the Lorentz and Gauss function) all weights are approximately 0.7. Note that except for Liu et al. the Gaussian weight is simply given by $1-\eta$.

As discussed in Schreier [24] and in our previous assessments of simple closed-form approximations of the Voigt function [25, 26, 12], the range of $y$ values encountered in molecular spectroscopy and atmospheric and astrophysical applications spans many orders of magnitude. In Fig. 2 (left) we compare the pseudo-Voigt approximations with reference values for $w(z)$, Eq. (3), obtained with the wofz code (algorithm originally based on Poppe and Wijers [27, 28] and later refined with ideas from Zaghloul and Ali [29]; Scientific Python (http://scipy.org) implementation scipy.special.wofz with at least 13 significant digits according to the documentation).

The function values shown on the left appear to be in reasonably good agreement with the reference. However, significant problems show up in the relative errors $|K-K_{\text{wofz}}|/K_{\text{wofz}}$ (right side) and for all approximations the maximum error is larger than 1 percent.

In the line center ($x=0$) the Whiting and Kielkopf approximations are exact by definition (assuming that the exponentially scaled complementary error function $K(0,y)=\mathop{\rm erfce}(y)$, Eq. (2), is evaluated exactly). The other approximations have errors of some percent for $y=1$, but smaller errors for smaller and larger $y$. Note that in the line center the Whiting and Kielkopf corrections vanish, i.e. $C(0,y)=0$, in contrast to the Matveev correction.

In the line wings only Matveev’s approximation (with and without correction) has relative errors decreasing with increasing $|x|$. Evaluating the Voigt function for very large $x$ ($x\leq 2000$ for $y=10$ and $x\leq 100$ for $y=1$) indicates that for all other approximations the errors become constant for large $x$. For Kielkopf this asymptotic error is about $10^{-4}$ for $y=10$ and $0.02$ for $y=0.1$ and $0.01$.

For small $y\leq 0.1$ Matveev’s approximation appears to be superior (with errors less than a few percent), however, for $y>1$ the correction term worsens the accuracy, and the results shown in the top of Fig. 2 have been obtained without the correction (12). For Kielkopf and Matveev, the largest errors occur for intermediate values of $x$, whereas the other approximations always fail for large $x$.

The contour plots of relative errors shown in Fig. 3 essentially confirm these findings. Except for the Liu et al. [21] Lorentz-Gauss combination all approximations have at least a small region where the relative accuracy is better than $10^{-4}$ (see [24] for a discussion of the $10^{-4}$ criterion). For Whiting and Kielkopf a relative high accuracy is achieved near the origin for small and large $y$, for Matveev and Thompson only for small $y$.

The maximum relative error indicated in the title of all subplots identifies Kielkopf’s code (with correction) as the most accurate approximation. However, this maximum error is relatively large for modern standards (about six percent in the wings), and for small $y$ Matveev’s code appears to be better.

Despite the significant accuracy problems of all pseudo-Voigt approximations it is nevertheless instructive to test the numerical performance of the codes. Simple tests within the IPython interpreter [30] indicate that the Kielkopf and Matveev approximations are somewhat slower than the optimized combination of the Humlíček [31] and Weideman [32] rational approximations [24]. However, evaluation of HNO3 cross sections in the microwave for a series of pressures and temperatures corresponding to Earth’s atmosphere in the $0\,\text{--}\,120\rm\,km$ altitude range is about a factor two slower with these two pseudo-Voigt approximations compared to the Humlíček–Weideman combination. For details see the online appendix D.

4 Summary and Conclusions

Closed-form expressions for the Voigt function based on combinations of Lorentz- and Gauss-type functions show significant accuracy problems, with relative errors in the percent range. Note that in this study we have only considered “pseudo-Voigt” approximations based on an analytical, closed-form expression for the weight of the Lorentz and Gauss components (as a function of $y$, the ratio of the Lorentz and Gauss width). In several studies pseudo-Voigt approximations have been used for analysis of experimental data, where the weight has been estimated by least squares fitting [e.g. 17, 33].

Our conclusions now are therefore similar to those given in Schreier [25, 26]: Closed-form expressions as presented here might be desirable for certain applications, but their quality is limited. In general approximations based on modern state-of-the-art numerical methods, e.g. rational approximations as discussed in Humlíček [34, 31], Weideman [32] and Schreier [24, 35] are recommended.

Acknowledgments

Financial support by the Deutsche Forschungsgemeinschaft — DFG (project SCHR 1125/3-1) is gratefully acknowledged.

Appendix A Translations

The Lorenz, Gauss, and Voigt profiles are defined as

[TABLE]

with normalization $\int g(\nu,\dots)\,\mathrm{d}\nu=1$ and half width at half maximum (HWHM) $\gamma$. The dimensionless variables of the Voigt function $K$ are defined as ratios

[TABLE]

In the definition of the pseudo-Voigt functions the following ratios of the widths are used frequently

[TABLE]

Appendix B Numerical Constants

The constant used in the Kielkopf [16] half width approximation is $k_{\text{e}}=0.0990$ and the coefficients of the correction term are

[TABLE]

The coefficients of the half width expansion (19) of Thompson et al. [18] are

[TABLE]

and the weight expansion (20) is defined with

[TABLE]

The coefficients of the Lorentz and Gauss weights used by Liu et al. [21] are

[TABLE]

Appendix C Matveev

As indicated in the footnote of subsection 2.3, some differences show up in the original work by Matveev [14] and in Titov and Haus [22]. Whereas the $y\propto\gamma_{\text{L}}$ factor in the correction for the width is required for dimension reasons, our tests depicted in Fig. 4 indicate that the weight factor $\eta$ (or $\zeta$ in the original work) has been forgotten by Titov and Haus [22].

Appendix D Computational efficiency

For a first, preliminary test of the speed of pseudo-Voigt approximations we have used the IPython [30] builtin “magic” function %timeit:

In [1]: from pseudoVoigt import *

In [2]: x=numpy.linspace(0.,100.,10001); y=1.0

In [3]: %timeit kielkopf(x,y)

424 $\mu$s $\pm$ 273 ns per loop (mean $\pm$ std. dev. of 7 runs, 1000 loops each)

For the three cases $y=10.0$, $y=1.0$, and $y=0.001$ (see Table 1) the Humlíček–Weideman combination suggested in Schreier [24] is somewhat faster than the Kielkopf and Matveev approximations. However, evaluation of one or two exponential(s) (Gaussians) and two fractions for a single $x$ by the Kielkopf and Matveev codes appears to be faster than the numerous multiplications required for the “brute-force” Weideman [32] 24-term rational approximation or the Humlíček [34] rational approximation “zpf16” generalized to 16 terms [35].

For a more realistic assessment, molecular cross sections are computed as required for high resolution atmospheric radiative transfer modeling, i.e. HNO3 cross sections in the $16\,\text{--}\,17\rm\,cm^{-1}$ interval resulting from the superposition of 2376 lines in $6\,\text{--}\,27\rm\,cm^{-1}$ are computed for a series of pressure and temperature pairs corresponding to altitudes $0\,\text{--}\,120\rm\,km$ in Earth’s atmosphere (see Schreier [24, 35] for more details). In Python the total time required with the Kielkopf and Matveev approximations is more than a factor two larger than with the Humlíček–Weideman combination.