Faddeeva function
The Faddeeva function, denoted w(z), is a special function of complex analysis defined as w(z) = e^{-z^2} \erfc(-iz), where \erfc(z) = \frac{2}{\sqrt{\pi}} \int_z^\infty e^{-t^2} \, dt is the complementary error function.[1] This scaling by e^{-z^2} provides numerical stability, particularly for large |z|, and distinguishes it from the unscaled complex error function.[1] An equivalent integral representation is w(z) = e^{-z^2} \left( 1 + \frac{2i}{\sqrt{\pi}} \int_0^z e^{t^2} \, dt \right).[1] Named after Russian mathematician Vera Nikolaevna Faddeeva, the function originated from tables computed by Faddeeva and Nikolai N. Terent'ev in 1954 to evaluate the probability integral for complex arguments, essential for early computational methods in numerical analysis.[2] It later gained prominence in the 1960s through its appearance in standard references on special functions and its adoption in plasma physics.[2] The function is also known as the Kramp function or complex probability function in some contexts.[3] As an entire function analytic everywhere in the complex plane, the Faddeeva function satisfies functional equations such as w(\bar{z}) = \overline{w(z)} for real parts and exhibits asymptotic behavior w(z) \sim \frac{i}{\sqrt{\pi} z} as |z| \to \infty in |\arg z| < \frac{3\pi}{4}.[1] It connects to other special functions, including Dawson's integral via F(z) = \frac{i \sqrt{\pi}}{2} \left( e^{-z^2} - w(z) \right).[4] and the Voigt profile, expressed using the real part of w(z) for a complex argument z related to the broadening parameters.[1] Numerical algorithms for its computation, such as those using continued fractions or series expansions, ensure high precision across the complex plane.[3] In applications, the Faddeeva function models dispersion in hot plasmas, as in the plasma dispersion function Z(\zeta) = i \sqrt{\pi} \, w(\zeta), crucial for analyzing linearized waves and instabilities.[2] It also arises in atomic and astrophysical spectroscopy through the Voigt profile for line broadening due to Doppler and pressure effects, aiding diagnostics in neutron diffraction, laser physics, and radiative transfer.[2] Further uses include electromagnetic wave propagation in ionized media and deconvolution of motional broadening in scattering experiments.[2]Definition
Mathematical definition
The Faddeeva function, denoted as w(z), is defined for a complex argument z by the equation w(z) = e^{-z^2} \operatorname{erfc}(-i z), where \operatorname{erfc} denotes the complementary error function.[3] This definition holds directly for \operatorname{Im}(z) \geq 0, with the function extended to the entire complex plane via analytic continuation, for example using the functional equation w(-z) = 2 e^{-z^2} - w(z).[3] The Faddeeva function thus serves as a scaled version of the complex complementary error function, incorporating the Gaussian factor e^{-z^2} to facilitate computations in various applied contexts.[3] As an entire function, w(z) is analytic everywhere in \mathbb{C}.[1]Relation to the error function
The Faddeeva function w(z) provides a natural extension of the real error function to the complex plane, unifying various forms of the error function under a single analytic expression valid for complex arguments z \in \mathbb{C}. It is defined in direct relation to the complementary error function as w(z) = e^{-z^2} \operatorname{erfc}(-i z), where \operatorname{erfc}(z) = 1 - \operatorname{erf}(z) and \operatorname{erf}(z) = \frac{2}{\sqrt{\pi}} \int_0^z e^{-t^2} \, dt. This identity allows computation of the complex complementary error function via \operatorname{erfc}(z) = e^{-z^2} w(i z), thereby extending the classical error function, originally defined for real arguments in probability and heat conduction problems, to broader applications in physics requiring complex-domain evaluations.[1][5] An equivalent integral representation links w(z) to the imaginary error function \operatorname{erfi}(z), defined as \operatorname{erfi}(z) = -i \operatorname{erf}(i z) = \frac{2}{\sqrt{\pi}} \int_0^z e^{t^2} \, dt: w(z) = e^{-z^2} \left( 1 + \frac{2i}{\sqrt{\pi}} \int_0^z e^{t^2} \, dt \right) = e^{-z^2} \left( 1 + i \operatorname{erfi}(z) \right). This form highlights how w(z) incorporates the oscillatory behavior of \operatorname{erfi}(z) for imaginary arguments, essential for modeling phenomena like wave propagation where real error functions alone are insufficient. The relation also connects to the standard error function via \operatorname{erf}(z) = 1 - e^{-z^2} w(i z), ensuring consistency across real and complex cases.[1][5] Furthermore, w(z) relates to the scaled complementary error function \operatorname{erfcx}(z) = e^{z^2} \operatorname{erfc}(z), which arises in asymptotic analyses and avoids underflow in numerical computations for large |z|: \operatorname{erfcx}(z) = w(i z). This equivalence facilitates efficient evaluation in scientific computing libraries, as \operatorname{erfcx}(z) maintains numerical stability for real positive arguments while w(z) handles the full complex plane. Through these identities, the Faddeeva function serves as a foundational tool for extending error function theory to complex analysis, with applications in spectroscopy and quantum mechanics relying on its analytic continuation properties.[5]Properties
Real and imaginary parts
The Faddeeva function w(z) for complex argument z = x + i y with y > 0 decomposes into real and imaginary parts as w(z) = U(x, y) + i V(x, y), where U(x, y) and V(x, y) are the Voigt functions. The real part U(x, y) corresponds to the Voigt profile, which physically represents the convolution of a Gaussian distribution (modeling Doppler broadening) with a Lorentzian distribution (modeling natural or pressure broadening), a fundamental form in spectral line shapes.[5] The imaginary part V(x, y) relates to the associated Hilbert transform or dispersion component of this convolution. These functions exhibit even and odd symmetries in the real variable: U(x, y) = U(-x, y) and V(x, y) = -V(-x, y).Sign inversion and scaling relations
The Faddeeva function satisfies a sign inversion relation given by w(-z) = 2 e^{-z^2} - w(z) for any complex argument z. This functional equation arises from the complementary property of the complementary error function, \operatorname{erfc}(u) + \operatorname{erfc}(-u) = 2, combined with the defining relation w(z) = e^{-z^2} \operatorname{erfc}(-i z). The relation is particularly useful in numerical computations, as it allows evaluation in the lower half of the complex plane using values from the upper half, where algorithms often perform more stably.[6] The Faddeeva function is also connected to the scaled (or repeated) complementary error function, defined as \operatorname{erfcx}(z) = e^{z^2} \operatorname{erfc}(z), through the identity w(z) = \operatorname{erfcx}(-i z). This equivalence follows directly from substituting the definition of \operatorname{erfcx} into the expression for w(z), noting that (-i z)^2 = i^2 z^2 = -z^2. For real positive scaling factors a > 0, the function transforms as w(a z) = e^{-a^2 z^2} \operatorname{erfc}(-i a z), or equivalently w(a z) = \operatorname{erfcx}(-i a z); no further algebraic simplification exists, but this form preserves the scaling in the argument while adjusting the exponential prefactor. A special case occurs for purely imaginary arguments z = i y with y > 0 real, where w(i y) = \operatorname{erfcx}(y), yielding a real-valued result that connects to applications in spectroscopy and diffusion problems.[7] Although the Faddeeva function w(z) is an entire function—analytic everywhere in the finite complex plane with no branch cuts—its integral representation requires careful consideration of the argument: w(z) = \frac{1}{i \pi} \int_{-\infty}^{\infty} \frac{e^{-t^2}}{t - z} \, dt, \quad \Im(z) > 0. This contour integral converges for arguments in the upper half-plane. For real z, the representation holds in the sense of the Cauchy principal value, ensuring continuity across the real axis. The principal value interpretation avoids singularities on the integration path and maintains consistency with the entire function property.[8]Derivatives
The first derivative of the Faddeeva function w(z) is given by \frac{\mathrm{d} w}{\mathrm{d} z}(z) = \frac{2 i}{\sqrt{\pi}} - 2 z \, w(z). This relation follows directly from the definition of w(z) as e^{-z^2} \operatorname{erfc}(-i z) and the known derivative of the complementary error function, and it forms the basis for the function's role in solving first-order linear differential equations of the form w' + 2 z w = \frac{2 i}{\sqrt{\pi}}.[9] Higher-order derivatives satisfy the recurrence relation \frac{\mathrm{d}^{n+1} w}{\mathrm{d} z^{n+1}}(z) = -2 z \frac{\mathrm{d}^n w}{\mathrm{d} z^n}(z) + 2 n \frac{\mathrm{d}^{n-1} w}{\mathrm{d} z^{n-1}}(z), for n = 1, 2, \dots, with initial conditions from the zeroth and first derivatives (w^{(0)}(z) = w(z) and the first derivative above). This three-term recurrence arises by repeated differentiation of the first-order differential equation satisfied by w(z) and enables efficient computation of higher derivatives without explicit series expansions.[9] In the context of series expansions, the higher derivatives of w(z) connect to Hermite polynomials through the relation to the error function, where the nth derivative involves terms like H_{n-1}(z) e^{-z^2}, providing a polynomial structure for analytic manipulations. The recurrence facilitates solving higher-order linear differential equations in applications such as boundary value problems in diffusion or quantum mechanics, where repeated differentiation of the governing equation yields systems closed under this relation.[9] For symmetry checks, the sign inversion property w(-z) = 2 e^{-z^2} - w(z) implies that odd-order derivatives satisfy w^{(2k+1)}(-z) = -w^{(2k+1)}(z), consistent with the recurrence.Representations
Integral representation
The Faddeeva function possesses a fundamental integral representation given by w(z) = \frac{i}{\pi} \int_{-\infty}^{\infty} \frac{e^{-t^{2}}}{z - t} \, dt for \Im(z) > 0, where the path of integration is along the real axis; for \Im(z) > 0, the pole at t = z lies off the path in the upper half-plane.[8] This form, originally derived in the context of higher transcendental functions, allows for analytic continuation to the entire complex plane. This representation holds as an entire function. An equivalent representation, symmetric about the origin, is w(z) = \frac{2z}{\pi i} \int_{0}^{\infty} \frac{e^{-t^{2}}}{t^{2} - z^{2}} \, dt, \quad \Im(z) > 0. This alternative integral emphasizes the even nature of the integrand for real t > 0 and facilitates certain computational approaches while maintaining the same convergence properties. When z = x is real, the integral requires the Cauchy principal value to handle the singularity at t = x: w(x) = \frac{i}{\pi} \, \mathrm{P.V.} \int_{-\infty}^{\infty} \frac{e^{-t^{2}}}{x - t} \, dt, \quad x \in \mathbb{R}. This principal value ensures the integral converges, as the contributions from symmetric limits around the pole cancel in a specific manner due to the odd symmetry of the kernel. The integral representation reveals an interpretation of the Faddeeva function in terms of the Hilbert transform of the Gaussian e^{-t^2}. For real x, w(x) = i \, H(e^{-t^2})(x), where H denotes the Hilbert transform H(f)(x) = \frac{1}{\pi} \mathrm{P.V.} \int_{-\infty}^{\infty} \frac{f(t)}{x - t} \, dt. More broadly, the form embodies a Gaussian convolution with a simple pole kernel, providing the analytic continuation of the Voigt profile in the complex plane and underscoring its role in dispersive phenomena.[8]Series expansions
The Faddeeva function admits a power series expansion around z = 0 that converges for all finite complex z, reflecting its status as an entire function. This expansion is particularly effective for numerical evaluation when |z| is small, typically |z| \lesssim 5, beyond which other representations are preferred for efficiency. The series is derived from the relation to the complementary error function and can be expressed in terms of a confluent hypergeometric function or an explicit sum. One compact form involves the Kummer confluent hypergeometric function {}_1F_1: w(z) = e^{-z^2} \left[ 1 + \frac{2 i z}{\sqrt{\pi}} \, {}_1F_1\left( \frac{1}{2}; \frac{3}{2}; z^2 \right) \right], where {}_1F_1(a; b; \zeta) = \sum_{n=0}^{\infty} \frac{(a)_n}{(b)_n} \frac{\zeta^n}{n!} with Pochhammer symbols (a)_n = a(a+1) \cdots (a+n-1) (and (a)_0 = 1), a = 1/2, b = 3/2, and \zeta = z^2. The radius of convergence is infinite due to the entire nature of {}_1F_1. The Pochhammer symbols admit closed forms using double factorials: (1/2)_n = (2n-1)!! / 2^n and (3/2)_n = (2n+1)!! / 2^n, yielding (1/2)_n / (3/2)_n = 1/(2n+1). Thus, the series simplifies to {}_1F_1\left( \frac{1}{2}; \frac{3}{2}; z^2 \right) = \sum_{n=0}^{\infty} \frac{z^{2n}}{(2n+1) n!}. Substituting gives the equivalent explicit expansion w(z) = e^{-z^2} \left[ 1 + \frac{2 i z}{\sqrt{\pi}} \sum_{n=0}^{\infty} \frac{z^{2n}}{(2n+1) n!} \right]. This form is obtained by substituting the power series for \erf(iz) into the definition w(z) = e^{-z^2} \erfc(-iz) = e^{-z^2} [1 + \erf(iz)].[10] The terms in the sum can be generated efficiently via recurrence relations, avoiding direct computation of factorials for large n. Let t_0 = 1 and, for n \geq 1, t_n = t_{n-1} \cdot \frac{(2n-1) z^2}{n (2n+1)}, so the sum is \sum_{n=0}^{\infty} t_n. These recurrences stem from the standard three-term relation for confluent hypergeometric series terms and ensure numerical stability for successive approximations. In practice, the number of terms required decreases with smaller |z|; for example, around 20 terms suffice for |z| \approx 3 to achieve double-precision accuracy. The derivatives of w(z), which generate higher-order Taylor coefficients via w(z) = \sum_{k=0}^{\infty} \frac{w^{(k)}(0)}{k!} z^k, are related to Hermite polynomials through the differential equation satisfied by w(z), allowing recursive computation of coefficients in the full power series expansion.Asymptotic approximations
For large values of |z| in the sector |arg z| < 3π/4, the Faddeeva function w(z) admits the following asymptotic expansion: w(z) \sim \frac{i}{z \sqrt{\pi}} \sum_{m=0}^{\infty} (-1)^m \frac{(1/2)_m}{z^{2m}}, where (1/2)_m = \frac{(2m-1)!!}{2^m} is the Pochhammer symbol with the double factorial defined as (2m-1)!! = 1 \cdot 3 \cdot 5 \cdots (2m-1) for m \geq 1 and (-1)!! = 1.[11] This expansion is derived directly from the asymptotic behavior of the complementary error function via the relation w(z) = e^{-z^2} \operatorname{erfc}(-i z), where the leading term arises from substituting the argument -i z into the known asymptotic series for \operatorname{erfc}(\zeta) with \zeta = -i z.[11] The series is divergent but asymptotic, meaning that truncating at the optimal number of terms (before the terms begin to increase) provides an approximation whose error decreases as |z| increases. Specifically, if the series is truncated after the nth term, the remainder satisfies |\operatorname{erfc}(z) - s_n(z)| \leq (1 + \chi(n+1)) times the (n+1)th term in magnitude, where \chi(k) is a monotonically increasing function bounded by \chi(k) \approx 2^{2k-1}/\sqrt{\pi k} for large k; this bound extends to w(z) through the relation to \operatorname{erfc}. For |ph z| ≤ π/4, the remainder has the same sign as the first neglected term when ph z = 0, and for π/4 ≤ |ph z| < π/2, it is bounded by \csc(2 |\ph z|) times the first neglected term.[11][12] The sector of validity |arg z| < 3π/4 is delimited by Stokes lines at arg z = ±3π/4, across which the dominant exponential behavior changes due to the Stokes phenomenon; beyond these lines, the expansion requires re-expansion to account for exponentially small subdominant contributions that become relevant.[11][12] An equivalent representation in continued fraction form, which converges to the same asymptotic series, is w(z) = \frac{i /(z \sqrt{\pi})}{1 + \cfrac{1/(2 z^2)}{1 + \cfrac{3/(2 z^2)}{1 + \cfrac{5/(2 z^2)}{1 + \cfrac{7/(2 z^2)}{1 + \ddots}}}}} . Here, the partial quotients are a_k = (2k-1)/(2 z^2) for k = 1, 2, 3, \dots, and this form is particularly useful for numerical evaluation in the outer region of the complex plane where |z| is large (typically |z| > 5–10, depending on precision requirements).[13] The convergents of this continued fraction provide successive approximations whose errors decrease monotonically, with the error after n stages being O(1/z^{2n+1}) as |z| → ∞ in the sector -π/4 < arg z < 5π/4, though optimal truncation aligns with the asymptotic series behavior.[13] This representation originates from the Laplace continued fraction associated with the integral form of the error function and has been adapted for efficient computation of w(z).[13]Applications
In plasma physics and spectroscopy
In plasma physics, the Faddeeva function is essential for deriving the dielectric permittivity ε(ω, k) in kinetic theory descriptions of warm, Maxwellian plasmas, where it appears through the plasma dispersion function in the dispersion relations for longitudinal waves and instabilities.[14] The longitudinal dielectric function takes the form ε(ω, k) = 1 - \frac{ω_p^2}{2 k^2 v_th^2} Z'(ζ), with ζ = ω / (k v_th √2) and the derivative Z' of the plasma dispersion function, enabling the analysis of phenomena like Landau damping.[15] The plasma dispersion function Z(ζ) is related to the Faddeeva function w(z) by the equation Z(ζ) = i √π w(ζ), allowing efficient numerical evaluation of dispersion relations in simulations of plasma waves.[16] In spectroscopy, the Faddeeva function underpins the Voigt profile, which models the absorption or emission lineshape as a convolution of Gaussian (Doppler) and Lorentzian (pressure or natural) broadening mechanisms, critical for interpreting spectral data from atomic and molecular transitions.[3] The Voigt function V(x; σ, γ) is proportional to the real part of the Faddeeva function evaluated at a complex argument: Re[w(x + i y)], where x = (ν - ν_0)/(√2 σ) is the normalized frequency detuning, y = γ/(√2 σ) incorporates the Lorentzian half-width γ relative to the Gaussian standard deviation σ, and ν_0 is the central frequency.[17] This representation facilitates accurate fitting of observed spectra in astrophysical and laboratory plasmas, capturing the asymmetric broadening effects without explicit convolution integrals.[18] Applications extend to neutron transport simulations, where the Faddeeva function enables direct, on-the-fly Doppler broadening of resonance cross-sections during Monte Carlo particle tracking, improving efficiency over pre-tabulated libraries by accounting for temperature-dependent thermal motion in nuclear reactors.[19] In atomic physics simulations, particularly for high-pressure plasmas, it models complex line profiles by incorporating the Faddeeva function to resolve Stark and van der Waals broadening alongside Doppler effects, aiding diagnostics in fusion and lighting device research.[18]In electromagnetics and other fields
In electromagnetics, the Faddeeva function is integral to modeling groundwave propagation for amplitude modulation (AM) radio signals over lossy terrain. The classical Sommerfeld problem of a vertical Hertzian dipole radiating above a planar earth with finite conductivity and permittivity leads to an attenuation function expressed in terms of the Faddeeva function, accounting for the complex integration over the Sommerfeld integral. This formulation captures the lateral wave and surface wave contributions essential for predicting signal strength in medium-wave broadcasting, where groundwaves provide reliable coverage up to hundreds of kilometers.[20][21] In materials science, the Brendel-Bormann model utilizes the Faddeeva function to parameterize the complex dielectric function of amorphous oxides and solids in the infrared regime. By representing the dielectric response as a superposition of damped oscillators with a Gaussian distribution of resonance frequencies, the model incorporates the Faddeeva function to yield absorption profiles with Gaussian shapes while ensuring compliance with Kramers-Kronig relations. This approach improves the accuracy of optical constant extraction for thin films and bulk materials compared to traditional Lorentzian oscillators.[22] In nuclear physics, the Faddeeva function enables efficient computation of temperature-dependent neutron cross-sections through Doppler broadening via the Voigt profile. In multipole representations, such as the windowed multipole method, it evaluates the convolution of Gaussian thermal motion with Lorentzian resonances, requiring multiple invocations per energy point for accurate on-the-fly broadening in Monte Carlo simulations. For nuclides with numerous poles, this can involve thousands of Faddeeva evaluations to resolve cross-section tails.[23][24] Beyond these areas, the Faddeeva function's integral representation supports convolution analyses in signal processing, interpreting Gaussian-Lorentzian blends for profile fitting in filtered signals. In heat transfer through complex media, it appears in analytical solutions for transient conduction near encapsulated heat sources, such as nuclear fuel pellets, where the dimensionless temperature field is expressed using the function to handle diffusive boundaries.[25]History
Early development and tabulation
The Faddeeva function emerged in the early 1950s as part of efforts to extend the real-valued probability integral, or error function, to complex arguments for solving practical problems in mathematical physics, including heat conduction in media with variable properties. Vera N. Faddeeva, a prominent Soviet mathematician specializing in computational methods, recognized the need for accurate numerical values of this function during her work on integral equations. Motivated by its integral representation, which generalizes the Gaussian integral to complex domains, Faddeeva initiated the computation and tabulation of the function, known as the W-function, to enable its use in engineering and physical applications.[26] In collaboration with Nikolai M. Terentyev, Faddeeva produced comprehensive numerical tables covering a broad range of complex arguments, computed with precision up to five decimal places using available mechanical calculators and early analytical techniques. Their joint effort addressed the challenges of evaluating the function for arguments where real and imaginary parts varied independently, providing essential data for researchers lacking direct computational tools. This work marked a foundational step in making the Faddeeva function accessible beyond theoretical derivations.[27] The resulting publication, Tables of Values of W-Functions and of the Probability Integral for Complex Argument (1954), compiled by Faddeeva and Terentyev under the editorship of V. A. Fok, was released by the State Publishing House for Physical-Mathematical Literature in Moscow. Spanning detailed grids for the function's real and imaginary parts, the tables facilitated immediate adoption in Soviet scientific computations and became a standard reference for complex error function evaluations.[27] In the 1950s, as Soviet plasma physics advanced amid thermonuclear research, these tables gained early prominence in the literature for computing the plasma dispersion function, crucial for analyzing wave stability in collisionless plasmas via the Vlasov equation. Works in kinetic theory from this era, including those exploring Landau damping, relied on Faddeeva and Terentyev's data to numerically resolve dispersion relations with complex frequencies.[28][14]Naming and algorithmic advancements
The Faddeeva function received its formal name in 1990 through the work of Gerard P. M. Poppe and Cornelis M. J. Wijers, who introduced the term in their algorithmic implementation published in the ACM Transactions on Mathematical Software, honoring the contributions of V. N. Faddeeva to the study of the complex error function in her 1961 monograph with N. M. Terent'ev. Prior to this, the function was commonly referred to as the complex error function or the plasma dispersion function, reflecting its roles in error function extensions and plasma physics applications. In some contexts, it is also known as the Kramp function, named after the mathematician Christian Kramp for his early work on related special functions, though this nomenclature remains less prevalent.[29] Algorithmic advancements for computing the Faddeeva function began in the late 1960s with Walter Gautschi's development of a method using power series expansions and asymptotic approximations for the complex error function, implemented as ACM Algorithm 363, which provided efficient evaluation up to machine precision for arguments in the first quadrant of the complex plane. In the 1970s, J. Humlíček advanced this further by proposing rational approximations tailored for the Voigt function, a real-valued special case of the Faddeeva function, segmenting the complex plane into regions for optimized polynomial and rational fits, achieving rapid computation with relative errors below 10^{-5} in spectroscopic applications.90152-5) These methods laid the groundwork for handling the function's oscillatory behavior and branch cuts. The 1990 algorithm by Poppe and Wijers marked a significant improvement, combining continued fractions, power series, and asymptotic expansions to cover the entire complex plane with high accuracy (up to 18 decimal digits) and reduced computational cost compared to Gautschi's approach, as detailed in ACM Algorithm 680. Building on this, J. A. C. Weideman introduced in 1994 a compact rational approximation valid across the whole complex plane, adjustable in degree for balancing speed and precision, which simplified implementations while maintaining errors under 10^{-10} for moderate arguments. Recent developments include a 2020 method using modified trapezoidal rules for integral representations, offering superior accuracy for large imaginary parts with convergence rates independent of the argument magnitude, and a 2024 recurrence-based algorithm that enhances stability and efficiency for scaled variants, outperforming prior codes in benchmarks for high-precision needs.Numerical implementations
Key algorithms
The evaluation of the Faddeeva function w(z) for complex arguments z requires numerical algorithms that balance accuracy, computational efficiency, and robustness across the complex plane. Key methods leverage its integral representation, series expansions, and asymptotic behaviors to handle different regions of z, particularly distinguishing between small and large |z|, while using domain symmetries.[30] The Poppe-Weideman algorithm employs a hybrid approach, utilizing power series expansions for small |z| (typically |z| < 1) and continued fraction approximations for larger |z|, with a tunable parameter to control truncation and ensure convergence. This method achieves relative errors below $10^{-13} across the first quadrant (\operatorname{Re}(z) \geq 0, \operatorname{Im}(z) > 0), extending to the full plane via reflection symmetries, and performs reliably even for extreme arguments like |\operatorname{Re}(z)| > 10^4 or small imaginary parts down to $10^{-20}. Its strengths include high precision suitable for scientific computing and avoidance of overflow issues in asymptotic regimes, though it requires more operations for very large |z|.[31] Humlicek's approximation approximates the real and imaginary parts of w(z) separately using rational functions, often in the form of partial fraction expansions with fixed coefficients, optimized for the Voigt profile (the real part of w(z)) in spectroscopy applications. The 12-term variant (cpf12) delivers 5–6 significant digits of accuracy (relative error < $5 \times 10^{-6}) over most of the complex plane, with modifications for small \operatorname{Im}(z) and asymptotic handling for |\operatorname{Re}(z)| + \operatorname{Im}(z) > 15. It excels in speed for repeated evaluations in Voigt-related computations, such as line broadening models, due to its low polynomial degree and Horner scheme evaluation, but accuracy degrades near the origin for very small \operatorname{Im}(z).[32][33] Trapezoidal rule methods compute w(z) by directly approximating its integral representation \int_0^\infty \exp(-t^2) / (z + it) \, dt (or equivalent forms) using modified trapezoidal quadrature on a truncated real-line contour in the complex plane, incorporating Cauchy's residue theorem to account for poles. With 11–12 quadrature points and step size h = \sqrt{\pi}/(N+1), these achieve absolute and relative errors below $2 \times 10^{-15} in the upper half-plane, with exponential convergence improving accuracy per additional point by a factor of about 23. They are particularly robust for all z \in \mathbb{C}, using symmetries, and offer high precision without series divergence issues, though they may be slower for bulk computations due to the fixed number of integral evaluations.[30] Comparisons among these algorithms highlight trade-offs: Poppe-Weideman provides superior accuracy (up to 13 digits) and reliability compared to Humlicek's 5–6 digits, but Humlicek is faster (often 2–5 times) for Voigt-specific tasks with moderate precision needs. Trapezoidal rules match or exceed Poppe-Weideman in accuracy (15 digits) while being competitive in speed for single evaluations, but require more setup for vectorized use; all handle different domains effectively through quadrant-specific formulas, with trapezoidal methods offering the most uniform performance across domains.[33][30]Software libraries
The Faddeeva function is implemented in various scientific computing libraries, often as part of broader special functions modules, leveraging algorithms such as continued-fraction expansions and series approximations for complex arguments. A prominent open-source C++ implementation, known as the Faddeeva Package, provides high-accuracy computation (typically 13 or more significant digits) of the Faddeeva function w(z) = e^{-z^2} \mathrm{erfc}(-iz) and related error functions like \mathrm{erf}(z), \mathrm{erfc}(z), \mathrm{erfcx}(z), \mathrm{erfi}(z), and Dawson's integral, using a hybrid approach that combines a continued-fraction method for large |z| with a power series method for small |z| or values near the real axis.[31] This package serves as the core for several downstream libraries and wrappers. In Python, the SciPy library includes the functionscipy.special.wofz(z), which directly computes w(z) for complex z, incorporating the Faddeeva Package starting from version 0.12.0 to ensure efficient and accurate evaluation across the complex plane.[34] For R users, the RcppFaddeeva package (archived on CRAN since 2022 but available from source or GitHub) wraps the Faddeeva Package's C++ code, exposing functions like Faddeeva_w(z) to compute w(z) and associated error functions with relative error control, enabling seamless integration into R workflows for statistical and physical simulations.[35]
Julia's SpecialFunctions.jl module provides the Faddeeva function via faddeeva(x) = erfcx(-im * x), where erfcx is the scaled complementary error function, built on the Faddeeva Package for complex arguments and supporting high-precision computations in scientific applications.[36] In the Wolfram Language (used in Mathematica), the Faddeeva function is accessible through the built-in Erfc[z] for complex z, as w(z) = \exp(-z^2) \cdot \mathrm{Erfc}[-I z], with the system handling analytic continuation and numerical stability natively. Additionally, a dedicated FaddeevaW function is available in the Wolfram Function Repository for explicit evaluation.[37]
For C and C++ developers, the standalone libcerf library offers an efficient implementation of the Faddeeva function alongside Dawson, Voigt, and complex error functions, based on the Faddeeva Package and optimized for self-contained use in numerical simulations without external dependencies.[38] MATLAB and GNU Octave users can employ wrappers from the Faddeeva Package, which compile to MEX files or Octave plugins for direct computation of w(z) with configurable relative error tolerances.[31] These implementations prioritize accuracy and performance, often achieving machine precision while avoiding overflow in the exponentially scaled terms.