Faddeyeva Function Overview

• The Faddeyeva function is a complex probability function that combines Gaussian error functions with Lorentzian profiles, playing a key role in spectroscopy and plasma diagnostics.
• Advanced computational methods like Algorithm 916 and modified Humlíček’s w4 ensure robust accuracy and efficiency across diverse complex domains.
• Optimized implementations in Fortran and similar languages support high-throughput simulations in radiative transfer and atmospheric science with controlled error bounds.

The Faddeyeva function, also known as the complex probability function, plays a central role in physical and mathematical modeling where the interplay between Gaussian and Lorentzian profiles is critical. Defined for complex arguments, it encapsulates both error-function behavior and complex-analytic properties, and underlies the Voigt profile's formulation used extensively in spectroscopy, atmospheric radiative transfer, plasma physics, and related fields. Efficient and robust computation of the Faddeyeva function is crucial because practical applications often require evaluation over large domains with strict accuracy constraints.

1. Mathematical Definition and Structural Properties

The Faddeyeva function $w(z)$ is succinctly defined by

$$w(z) = e^{-z^2}\ \mathrm{erfc}(-i z), \quad z = x + i y \in \mathbb{C}$$

where $\mathrm{erfc}$ denotes the complementary error function. Alternate formulations include the Cauchy principal value integral representation: $$w(z) = \frac{i}{\pi} \int_{-\infty}^{\infty} \frac{e^{-t^2}}{z - t}\ dt$$ which demonstrates analyticity with a branch cut solely on the real axis. Series and asymptotic expansions yield precise computational strategies for different $z$ regimes.

The Faddeyeva function is intimately related to Dawson’s integral via

$$w(z) = e^{-z^2} + \frac{2i}{\sqrt{\pi}}\ \mathrm{daw}(z)$$

and possesses representations in terms of the confluent hypergeometric function: $$w(z) = e^{-z^2} \left[ 1 + \frac{2i}{\sqrt{\pi}} z \cdot {}_1F_1\left(\frac12; \frac32; z^2\right) \right]$$ where ${}_1F_1$ is Kummer’s confluent hypergeometric function (Nijimbere, 2017).

Key functional properties include the differential recursion: $w'(z) = -2z\, w(z) + \frac{2i}{\sqrt{\pi}}$ and the symmetries

$$w(-z) = 2 e^{-z^2} - w(z), \quad \overline{w(z)} = w(-\overline{z})$$

which ensure analytic consistency across the complex plane.

2. Connections to the Voigt Profile and Applied Contexts

The real part of $$w(z) = e^{-z^2}\ \mathrm{erfc}(-i z), \quad z = x + i y \in \mathbb{C}$$0 directly yields the Voigt function,

$$w(z) = e^{-z^2}\ \mathrm{erfc}(-i z), \quad z = x + i y \in \mathbb{C}$$1

the canonical line-shape profile modeling Doppler and pressure broadening in spectroscopy and radiative transfer problems. The imaginary part, $$w(z) = e^{-z^2}\ \mathrm{erfc}(-i z), \quad z = x + i y \in \mathbb{C}$$2, is relevant in dispersion and resonance phenomena such as magnetic resonance and astrophysical plasma diagnostics. Applications routinely require $$w(z) = e^{-z^2}\ \mathrm{erfc}(-i z), \quad z = x + i y \in \mathbb{C}$$3–$$w(z) = e^{-z^2}\ \mathrm{erfc}(-i z), \quad z = x + i y \in \mathbb{C}$$4 evaluations per run in large-scale numerical simulations (Zaghloul, 2015). Robustness is mandated to avoid catastrophic loss of precision, especially near the real axis.

Further, the plasma dispersion (Fried–Conte) function and Jackson’s function are direct relatives: $$w(z) = e^{-z^2}\ \mathrm{erfc}(-i z), \quad z = x + i y \in \mathbb{C}$$5 with the structure of $$w(z) = e^{-z^2}\ \mathrm{erfc}(-i z), \quad z = x + i y \in \mathbb{C}$$6 underpinning these special functions in plasma wave theory and transport modeling (Nijimbere, 2017).

3. Computational Methodologies: Algorithms and Approximations

Numerical evaluation of $$w(z) = e^{-z^2}\ \mathrm{erfc}(-i z), \quad z = x + i y \in \mathbb{C}$$7 must be accurate, stable, and efficient across a vast complex domain. Several major algorithmic frameworks have been established:

Algorithm 916 (Zaghloul & Ali, 2011) and Successors

Algorithm 916 employs:

• Region partitioning into small-$$w(z) = e^{-z^2}\ \mathrm{erfc}(-i z), \quad z = x + i y \in \mathbb{C}$$8 (power series), moderate-$$w(z) = e^{-z^2}\ \mathrm{erfc}(-i z), \quad z = x + i y \in \mathbb{C}$$9 (hybrid power series and continued fraction), and large-$\mathrm{erfc}$0 (asymptotic Laplace continued fraction) zones.
• Salzer’s expansion for $\mathrm{erfc}$1 is used to derive convergent, stable series representations for both $\mathrm{erfc}$2 and $\mathrm{erfc}$3.
• Dynamically monitored convergence, symmetry exploitation, and efficient recurrence for exponentials.
• Accuracy vs efficiency trade-off, using a user parameter (“tiny”), controlling the ultimate precision and computational cost (Zaghloul et al., 2011, Zaghloul, 2015).

Efficiency improvements (Algorithm 916 v2) introduced fixed-cycle summation (double precision achievable with $\mathrm{erfc}$4 terms), precomputed exponentials, vectorized grid handling, and aggressively tuned region boundaries, yielding up to $\mathrm{erfc}$5 speedup at highest accuracy, and nearly $\mathrm{erfc}$6 for low-accuracy modes (Zaghloul, 2015).

Humlíček’s w4 Algorithm and Region-Reform

Humlíček’s “w4” partitions the $\mathrm{erfc}$7 plane into four regions (asymptotic, exp-times-rational, Maclaurin expansion, and rational approximations) using simple border functions, selecting analytic formulas tailored for each (Zaghloul, 2015). However, cancellation in Region II (exp-times-rational) causes severe loss of precision for $\mathrm{erfc}$8, violating the nominal $\mathrm{erfc}$9 relative error claim. The minimal reform is to reroute small-$$w(z) = \frac{i}{\pi} \int_{-\infty}^{\infty} \frac{e^{-t^2}}{z - t}\ dt$$0 inputs to the more stable Region IV approximation:

• Original border: $$w(z) = \frac{i}{\pi} \int_{-\infty}^{\infty} \frac{e^{-t^2}}{z - t}\ dt$$1 (Region II).
• Revised: enforce $$w(z) = \frac{i}{\pi} \int_{-\infty}^{\infty} \frac{e^{-t^2}}{z - t}\ dt$$2 for Region II; else, revert to Region IV formula. This change eliminates the "failure tongue" along $$w(z) = \frac{i}{\pi} \int_{-\infty}^{\infty} \frac{e^{-t^2}}{z - t}\ dt$$3 with negligible performance penalty, restoring sub-$$w(z) = \frac{i}{\pi} \int_{-\infty}^{\infty} \frac{e^{-t^2}}{z - t}\ dt$$4 accuracy across the domain (Zaghloul, 2015, Zaghloul, 2015).

Algorithm 985 (Zaghloul, 2017): Piecewise Rational and Asymptotic

Algorithm 985 synthesizes and blends classical rational approximations (Hui, Humlíček) with Laplace continued-fraction asymptotics. The $$w(z) = \frac{i}{\pi} \int_{-\infty}^{\infty} \frac{e^{-t^2}}{z - t}\ dt$$5-plane is partitioned into regions by $$w(z) = \frac{i}{\pi} \int_{-\infty}^{\infty} \frac{e^{-t^2}}{z - t}\ dt$$6 and $$w(z) = \frac{i}{\pi} \int_{-\infty}^{\infty} \frac{e^{-t^2}}{z - t}\ dt$$7 cuts, applying:

• Asymptotic convergents (Laplace continued fraction) for $$w(z) = \frac{i}{\pi} \int_{-\infty}^{\infty} \frac{e^{-t^2}}{z - t}\ dt$$8,
• Humlíček’s Region IV rational for small $$w(z) = \frac{i}{\pi} \int_{-\infty}^{\infty} \frac{e^{-t^2}}{z - t}\ dt$$9 and $z$0,
• Hui’s rational $z$1 elsewhere. Piecewise selection gives uniform accuracy $z$2 in both $z$3 and $z$4 with competitive per-evaluation cost (Zaghloul, 2017).

Poppe & Wijers’s Algorithm 680 and Corrections

Algorithm 680 uses a hybrid Taylor and Laplace continued fraction method with boundaries tuned for maximum throughput and typically $z$5–$z$6 digit precision. However, original region boundaries force continued fraction use for $z$7 and $z$8, where catastrophic cancellation destroys accuracy. The “upward Taylor–Dawson” patch corrects this; for $z$9 and $$w(z) = e^{-z^2} + \frac{2i}{\sqrt{\pi}}\ \mathrm{daw}(z)$$0, a $$w(z) = e^{-z^2} + \frac{2i}{\sqrt{\pi}}\ \mathrm{daw}(z)$$1-term Taylor expansion about the real axis, with derivatives of $$w(z) = e^{-z^2} + \frac{2i}{\sqrt{\pi}}\ \mathrm{daw}(z)$$2 calculated from recurrence relations involving Dawson’s integral, restores full precision (Zaghloul, 2019).

4. Accuracy, Efficiency Trade-Offs, and Benchmarks

Optimization of $$w(z) = e^{-z^2} + \frac{2i}{\sqrt{\pi}}\ \mathrm{daw}(z)$$3 evaluation requires clear trade-offs among accuracy, execution time, and robustness, particularly for high-volume scenarios:

• Algorithm 916 and its efficiency-improved variant allow explicit user control of significant digits ($$w(z) = e^{-z^2} + \frac{2i}{\sqrt{\pi}}\ \mathrm{daw}(z)$$4). Empirical tuning matches the number of series cycles and the asymptotic boundary ($$w(z) = e^{-z^2} + \frac{2i}{\sqrt{\pi}}\ \mathrm{daw}(z)$$5) to the desired accuracy (see Table below).
• Fortran implementations are provided with precision-kind selection, vectorization, handling of scalar or array arguments, and explicit status codes (Zaghloul, 2015).

Mode	Z_border	Ncycles	Run time (s)	Relative Error
sdgts=13	$$w(z) = e^{-z^2} + \frac{2i}{\sqrt{\pi}}\ \mathrm{daw}(z)$$6	13	1.55	$$w(z) = e^{-z^2} + \frac{2i}{\sqrt{\pi}}\ \mathrm{daw}(z)$$7
sdgts=8	$$w(z) = e^{-z^2} + \frac{2i}{\sqrt{\pi}}\ \mathrm{daw}(z)$$8	9	0.70	$$w(z) = e^{-z^2} + \frac{2i}{\sqrt{\pi}}\ \mathrm{daw}(z)$$9
sdgts=4	$$w(z) = e^{-z^2} \left[ 1 + \frac{2i}{\sqrt{\pi}} z \cdot {}_1F_1\left(\frac12; \frac32; z^2\right) \right]$$0	7	0.49	$$w(z) = e^{-z^2} \left[ 1 + \frac{2i}{\sqrt{\pi}} z \cdot {}_1F_1\left(\frac12; \frac32; z^2\right) \right]$$1

Benchmarks:

• For $$w(z) = e^{-z^2} \left[ 1 + \frac{2i}{\sqrt{\pi}} z \cdot {}_1F_1\left(\frac12; \frac32; z^2\right) \right]$$2 grid points, Algorithm 916 (13 digits) achieves $$w(z) = e^{-z^2} \left[ 1 + \frac{2i}{\sqrt{\pi}} z \cdot {}_1F_1\left(\frac12; \frac32; z^2\right) \right]$$3 s with accuracy better than $$w(z) = e^{-z^2} \left[ 1 + \frac{2i}{\sqrt{\pi}} z \cdot {}_1F_1\left(\frac12; \frac32; z^2\right) \right]$$4; the reformed w4 (4 digits) attains $$w(z) = e^{-z^2} \left[ 1 + \frac{2i}{\sqrt{\pi}} z \cdot {}_1F_1\left(\frac12; \frac32; z^2\right) \right]$$5 s with maximal error below $$w(z) = e^{-z^2} \left[ 1 + \frac{2i}{\sqrt{\pi}} z \cdot {}_1F_1\left(\frac12; \frac32; z^2\right) \right]$$6.
• In double precision, Fortran implementations of the improved algorithm consistently outperform legacy codes in both speed and accuracy, matching or exceeding all competitors except in the lowest-accuracy regime, where Humlíček’s w4 is slightly faster ($$w(z) = e^{-z^2} \left[ 1 + \frac{2i}{\sqrt{\pi}} z \cdot {}_1F_1\left(\frac12; \frac32; z^2\right) \right]$$7 μs per call) but no longer outperforms the reformed variant in robustness (Zaghloul, 2015).

5. Implementation Considerations and Best Practices

Efficient and correct usage involves the following:

• Always set the accuracy-control parameter explicitly to avoid local defaults inappropriate for derivative-sensitive or high-precision tasks.
• For real-world large-scale applications (e.g., atmospheric line-by-line codes, plasma transport solvers), Fortran bindings are recommended for maximal throughput; Matlab or Scilab implementations suit prototyping and moderate grid sizes.
• For low-accuracy settings suitable when data-model error dominates (e.g., Voigt profile modeling limited by molecular parameters), sdgts of $$w(z) = e^{-z^2} \left[ 1 + \frac{2i}{\sqrt{\pi}} z \cdot {}_1F_1\left(\frac12; \frac32; z^2\right) \right]$$8–$$w(z) = e^{-z^2} \left[ 1 + \frac{2i}{\sqrt{\pi}} z \cdot {}_1F_1\left(\frac12; \frac32; z^2\right) \right]$$9 provides large speed gains with controlled, guaranteed error below ${}_1F_1$0. For double-precision spectral modeling or derivative computation, prefer ${}_1F_1$1–${}_1F_1$2 digits.
• Verification against reference codes such as Algorithm 916 with tight tolerances is recommended after integration or adaptation of any rational-approximation or continued-fraction-based code.
• Whenever the input domain is dominated by large ${}_1F_1$3, further speed-ups can be realized by adjusting asymptotic zone boundaries and region testing following local performance evaluation (Zaghloul, 2015).

6. Summary of Current Best Practices and Limitations

Current state-of-the-art guarantees robust, uniform accuracy for the Faddeyeva function throughout the complex plane by explicit domain partitioning and careful algorithm selection:

• High-accuracy needs are best met by Algorithm 916 and its successors—no known catastrophic breakdowns across any ${}_1F_1$4 regime.
• Reformed regional approaches to Humlíček’s w4 and Algorithm 680 remove previously intractable loss-of-precision zones near ${}_1F_1$5 without performance penalty.
• Algorithm 985 offers a streamlined alternative for moderate- to high-accuracy needs that blends optimal rational and continued-fraction formulas (Zaghloul, 2017).

Persisting limitations include residual slowdowns for extreme accuracy settings (where series must be pushed to full machine precision and beyond) and the requirement for careful branch and boundary management to preclude switching artifacts in transition zones. No single rational approximation is globally sufficient; domain-adapted hybrid schemes remain mandatory for modern high-throughput scientific computing.

7. References to Key Computational Schemes

The following references detail the foundational and specialized algorithms:

• Algorithm 916 (Zaghloul & Ali, 2011) (Zaghloul et al., 2011)
• Efficiency improvements and Fortran translation (Zaghloul, 2015)
• Algorithm 985: robust rational–asymptotic blend (Zaghloul, 2017)
• Humlíček’s w4 algorithm and region-reform (Zaghloul, 2015)
• Algorithm 680 correction: upward Taylor patch (Zaghloul, 2019)
• Analytical and hypergeometric representations (Nijimbere, 2017)

These works provide explicit pseudocode, parameter tables, benchmark survey results, and full implementation guidance suitable for immediate research- and production-level deployment.