A rational approximation for efficient computation of the Voigt function in quantitative spectroscopy

We present a rational approximation for rapid and accurate computation of the Voigt function, obtained by residue calculus. The computational test reveals that with only $16$ summation terms this approximation provides average accuracy ${10^{- 14}}$ over a wide domain of practical interest $0<x<40,000$ and ${10^{- 4}}<y<{10^2}$ for applications using the HITRAN molecular spectroscopic database. The proposed rational approximation takes less than half the computation time of that required by Weideman's rational approximation. Algorithmic stability is achieved due to absence of the poles at $y \geqslant 0$ and $ - \infty<x<\infty $.


Introduction
The Voigt function is widely used and finds broad applications in many scientific disciplines [1,2,3,4,5]. It is commonly applied in Applied Mathematics, Physics, Chemistry and Astronomy as it describes the line profile behavior that occurs due to simultaneous Lorentz and Doppler broadening effects; the Lorentz broadening is observed as a result of the Heisenberg uncertainty principle and chaotic multiple collisions of the particles while the Doppler broadening appears due to velocity distribution of the particles.
The Voigt function can describe the spectral properties in the photon emission or absorption of atmospheric gases [6,7,8,9,10] and celestial bodies [11]. It is also widely used in crystallography [12] and can be utilized in many other spectroscopic applications, for example, to characterize the photo-luminescent properties of nanomaterials [13] or to determine the hyper structure of an isotope [14] and so on.
Mathematically, the Voigt function is a convolution integral of the Cauchy and Gaussian distributions [1,2,3,4,5] K (x, y) = y π ∞ −∞ e −t 2 y 2 + (x − t) 2 dt (1) and represents the real part of the complex probability function [2,3] where z = x + iy is a complex argument. The complex probability function can be expressed explicitly as a superposition of the real and imaginary parts W (x, y) = K (x, y) + iL (x, y), where its imaginary part is given by [2,3] L (x, y) Another closely related function is the complex error function, also known as the Faddeeva function [3,15,16,17,18,19] There is a relation between complex probability function (2) and complex error function (4) In order to describe spectral characteristics of a system with high resolution, intense computation is required. For example, in a line-by-line radiative transfer simulation to resolve some problems associated with inhomogeneity, the Earth's or other planetary atmosphere can be divided up to 1000 layers [6,7]. Taking into account that computation requires a nested loop procedure in order to adjust properly for the fitting parameters for each atmospheric layer that may contain many different molecular species, the total number of the computed points may exceed hundreds of millions. Since in a radiative transfer model the computation of spectral broadening profiles requires considerable amount of time, a rapid approximation of the Voigt function is very desirable [6,7]. Consequently, the rapid and accurate computation of the Voigt/complex error function still remains topical (see for example an optimized algorithm in the recent work [20]).
In this work we present a new rational approximation of the Voigt function for efficient computation. Due to absence of the poles at y 0 and −∞ < x < ∞ this rational approximation enables stability in algorithmic implementation.

Derivation of the rational approximation
The complex error function (4) can also be expressed in an alternative form as [21,22] w (x, y) where its real and imaginary parts are respectively. By changing sign of the variable x to negative in the last two equations above, we can see the symmetric properties of the complex error function Consequently, it follows that and It is worth noting that with these identities and equation (4) we can also obtain two interesting relations for the real and imaginary parts of the error function of complex argument as follows Consequently, we can write Re [w (x, y = 0)] = K (x, y = 0) ≡ exp (−x 2 ) and from the identity (6) it immediately follows that In our recent publication we have shown that a sampling methodology based on incomplete expansion of the sinc function leads to a new series approximation of the complex error function [23] where the coefficients are 75, h = 0.25, M = 5 and N = 23. As we can see, the integer on upper limit of the summation in this approximation is equal to 2 M −1 . However, this restriction can be omitted and application of the series approximation above can be generalized for any integer.
Consider the following limit for the sinc function [23] where we imply that the sinc function is defined as Change of the integer variable 2 M −1 → m max in this limit leads to This signifies that if the integer m max is large enough, it retains all properties required to approximate the sinc function that can be used for sampling (see sampling methodology in [23] for details). Consequently, we can generalize the approximation (8) of the complex error function for an arbitrary integer m max as follows where the corresponding coefficients are rewritten as Combining identity (7) and approximation (9) together at y = 0 yields an exponential function approximation Figure 1 shows the difference ε (t) between the original exponential function exp (−t 2 ) and its approximation As we can see from this figure, even with only m max = 16 summation terms the difference ε (t) is very small and remains within the narrow range ±5 × 10 −10 . This confirms a rapid convergence of the exponential function approximation that makes it suitable for numerical integration. Specifically, this series approximation can be further used to replace the original exponential function exp (−t 2 ) from the integrand in integral equation (1) as follows Consider the series approximation (10) of the Voigt function in more detail. The integrand in this integral is analytic everywhere over the entire complex plain except 2 + 4m max isolated points where singularities are observed. However, as we take a contour integral only on the upper complex plane, for example as a semicircle C ccw with infinite radius in counterclockwise (CCW) direction, the quantity of isolated points is reduced twice and becomes equal to 1 + 2m max .
Lastly, substituting the corresponding isolated points inside the domain enclosed by contour C ccw : into the Residue Theorem's formula that for our specific case is expressed in where f (t) is the integrand of integral (10), we find a new series approxima-tion of the Voigt function Since an algorithm involving complex numbers takes extra time, it would be very desirable to exclude them in computation. Thus, after some trivial rearrangements of the equation above, it can be represented in a simplified form as the series approximation consisting of the real variables and constants only where As the constants α m , β m and γ m are independent of the input parameters x and y, the obtained series (12) is a rational approximation.

Implementation
Since the Voigt function is an even with respect to the parameter x and odd with respect to the parameter y: it is sufficient to consider the values x and y only from the I st and II nd quadrants in order to cover the entire complex plane. Consequently, in algorithmic implementation it is reasonable to take the second input parameter by modulus as |y| and compute the Voigt function according to the scheme Thus, if the parameter y is negative, we first take it by absolute value and, after computation, simply change the sign of the computed result to opposite. It should also be noted that taking the argument |y| is advantageous in implementation as it prevents computational overflow and enables an algorithmic stability (see Appendix A for details). The series approximation (12) alone covers the domain 0 < x < 40, 000 and 10 −4 < y < 10 2 , required in applications using the HITRAN molecular spectroscopic database [24]. In general, it provides accurate results while y 10 −6 . However, this approximation may be used only to cover a smaller domain 0 x 15 and 10 −6 y 15 that is considered most difficult for rapid and accurate computation.
In our recent publication we have shown that the following approximation (see equation (6) in [25]) can be effective for computation in the narrow domain 0 x 15 and 0 y < 10 −6 along x-axis: is the Dawson's integral. As argument x in the Dawson's integral is real, its implementation is not problematic and several efficient approximations can be found in literature [26,27,28].
When the input parameters x and y are large enough (say when the condition |x + iy| > 15 is satisfied), many rational approximations become effective for accurate and rapid computation. For example, the Gauss-Hermit quadrature or the Taylor expansion can be effectively implemented (see for example [4] for details).
A Matlab source code for computation of the Voigt/complex error function that covers the entire complex plane can be accessed through Matlab Central, file ID: #47801 [29]. This code has been developed by our research group and can be used for verification of the computed results. The domain divisions for computation of the Voigt function with complete coverage of the complex plane can be developed similarly.
In order to demonstrate the computational efficiency of the series approximation (12), the comparison with the Weideman's rational approximation has been made (see equation (38-I) and corresponding Matlab code in [19]). Such a choice is justified since the Weideman's approximation is one of the most rapid for computation of the Voigt/complex error function. The computational testing we performed by using a typical desktop computer shows that with same number of the summation terms m max = 16 (default integer in Matlab code in [19] is also 16), the algorithm based on series approximation (12) is faster in computation than that of based in the Weideman's rational approximation by factors about 2.2 and 2.7 for input arrays x and y consisting of 5 and 50 million elements, respectively (see the Matlab source code with implementation of the series approximation (12) in Appendix B). This is mainly because the Weideman's rational approximation computes simultaneously both the real K (x, y) and imaginary L (x, y) parts, while the rational approximation (12) computes only the real part K (x, y) of the complex error function w (x, y). It should be noted that in most practical applications the imaginary part L (x, y) (3) of the complex error function is not needed and simply ignored. Moreover, due to rapid convergence of the series approximation (12) we may decrease the number of the summation terms. In particular, at m max = 12 the computational acceleration of the Voigt function can be further gained by about 30%. Therefore, the application of the series approximation (12) may be advantageous especially for intense computations with extended input arrays.

Error analysis
In order to quantify accuracy of the series approximation (12), it is convenient to define the relative error as , y) is the reference. The highly accurate reference values can be generated, for example, by using the Algorithm 680 [18,30] or recently published Algorithm 916 [31]. Figures 2a and 2b show the logarithm log 10 ∆ of the relative error of the series approximation (12) at m max = 16. The domain required for coverage of the HITRAN molecular spectroscopic database is 0 < x < 40, 000 and 10 −4 < y < 10 2 [7,32] while the domain 0 x 15 and 10 −6 y 15 is the most difficult for accurate and rapid computation of the Voigt function. Therefore, we will consider the accuracy behavior within the HITRAN subdomain and narrow band domain 0 x 15 ∩ 10 −4 y 15 and 0 x 15 ∩ 10 −6 y 10 −4 separately as shown in Figs. 2a and 2b, respectively.
As we can see from Fig. 2a, within the HITRAN subdomain the accuracy of the series approximation is quite uniform and better than 10 −14 over most of this area. Although the accuracy deteriorates with decreasing y, it, nevertheless, remains high and better than 10 −9 . Another advantage is that the area where the accuracy deteriorates is relatively small. Particularly, the area where accuracy is worse than 10 −13 (yellow and red colors) is smaller than 2% of the domain's total area.
With randomly taken input parameters x and y, it is determined that the average accuracy over the domain of practical interest 0 < x < 40, 000 and 10 −4 < y < 10 2 is 10 −14 . Although the series approximation (12) can cover this domain accurately, it may be implemented only within domain 0 x 15 and 10 −6 y 15 that is the most difficult for accurate and rapid computation.
In the narrow band shown in the Fig. 2b, the accuracy deteriorates further with decreasing y. However, it still remains high and better than 10 −8 . In particular, the best and worst accuracies in the narrow band domain 0 x 15 ∩ 10 −6 y 10 −4 exceed 10 −10 (yellow color) and 10 −8 (dark red color), respectively.
In modern applications requiring the HITRAN molecular spectroscopic data-base, the accuracy of the Voigt function should be 10 −6 . Therefore, we may reduce the integer m max in the series approximation (12) from 16 to 12 in order to gain computational acceleration. The number of the summation terms, determined by the integer m max , is quite sensitive to the small parameter value h. We have found empirically that at m max = 12 the best accuracy can be achieved by taking h = 0.293. Figure 3a depicts the logarithm log 10 ∆ of the relative error of the series approximation (12) at m max = 12 in the HITRAN subdomain 0 x 15 and 10 −4 y 15. One can see that in the HITRAN subdomain the accuracy is better than 10 −8 . Figure 3b illustrates the logarithm log 10 ∆ of the relative error of the series approximation (12) at m max = 12 in the narrow band domain 0 x 15 and 10 −6 y 10 −4 . Despite only 12 summation terms involved in the series approximation (12), the accuracy within this domain is better than 10 −6 . For comparison, to achieve the same accuracy 10 −6 at y 10 −5 , the Weideman's approximation requires 32 summation terms (see Fig. 4 in [22] for details). Thus, we can see that the series approximation (12) may be useful and convenient in spectroscopic applications.

Conclusion
A rational approximation for rapid and accurate computation of the Voigt function is presented. With only 16 summation terms, the proposed rational approximation provides average accuracy 10 −14 in the domain of practical interest 0 < x < 40, 000 and 10 −4 < y < 10 2 that is needed for applications using the HITRAN molecular spectroscopic database. The computational test shows that the algorithm based on series approximation (12) is more rapid in computation than that of based on the Weideman's rational approximation by factor greater than 2. Algorithmic stability is achieved since the proposed series approximation (12) contains no poles at y 0 and −∞ < x < ∞.

Appendix A
According to definition of the κ-function (12) we can write the following identity and since where the right side of the identity (A.2) is the m th denominator of κ (x, y + ς/2), it is sufficient to show that the poles do not exist on the right side of the identity (A.1) when both variables x and y are real such that y 0 in order to prove that κ (x, y + ς/2) has no poles under same conditions. The proof is not difficult. Let us equate the left side of identity (A.2) to zero Since the constants C m , ς are real valued and since ς > 0, y 0, these solutions {x 1 , x 2 , x 3 , x 4 } must be always complex. However, the complex solutions {x 1 , x 2 , x 3 , x 4 } contradict our initial assumption that x is real. Similarly, four possible solutions of equation (A.3) with respect to the variable y are y 1 = −i (x − C m ) − ς/2, y 2 = i (x − C m ) − ς/2, y 3 = −i (x + C m ) − ς/2 and y 4 = i (x + C m ) − ς/2. Since the constants C m , ς are real and positive, the solutions {y 1 , y 2 , y 3 , y 4 } must be either complex at x = C m or negative and equal to −ς/2 at x = C m . However, the complex or negative solutions {y 1 , y 2 , y 3 , y 4 } contradict our initial assumption that y 0. Due to these contradictions we must conclude that there are no poles in identity (A.3) under the conditions {x, y} ∈ R such that y 0.
The absence of the poles signifies that while the arguments is taken by absolute value as |y|, the function κ (x, |y| + ς/2) will never encounter division to zero that leads to computational overflow. That is why taking the input parameter by absolute value as |y| is advantageous since this approach provides stable performance of the algorithm.