Inhomogeneous broadening in the time domain

2.1 Probability formalism for dispersion

This section aims to formulate the material dispersion through the concept of photon absorption probabilities (or broadening functions ^[3]) G _i (x), enabling generalization of classical dielectric laws from Lorentz broadening to arbitrary distributions. We start with a representation of complex relative permittivity in the time and frequency domains, connected via the Fourier transform^[4] (FT, F )

where, for generality, standard high-frequency permittivity (ɛ _∞) and conductivity terms (with DC electric conductivity σ _e) are assumed [53].

The dispersion terms χ ̂ i ( ω ) are defined as ideal unbroadened susceptibilities χ ̂ i 0 ( ω ) convolved (broadened by) absorption probabilities G _i (x), which must be valid Probability Density Functions (PDFs) [54], i.e., nonnengative with full probability support, G _i (x) ≥ 0 and ∫ − ∞ ∞ G i ( x ) d x = 1 ,

Each PDF G _i (x) is parameterized by the mean (μ _i ), variance σ i 2 , and/or other higher-order statistical moments and parameters. For now, we assume symmetric distributions, G _i (−x) = G _i (x), and zero means, μ _i = 0, so that the center frequency of susceptibility doesn’t change with broadening.

In the time domain, obtained via the inverse FT and applying the convolution theorem, Eq. (2) reads

where the symmetric broadening functions G _i (x) contribute through its characteristic functions (CF) φ _i (t) [54]. Standard CF properties include boundedness and zero-centered unity, |φ _i (t)| ≤ 1 and φ _i (0) = 1; moreover, if the PDF is symmetric, its CF is real-valued.

As a clear example, we reformulate the classical Lorentz oscillator using the proposed formalism

Here f, γ are oscillator’s strength and damping parameters, while Ω and ω 0 = Ω 2 − γ 2 are resonance and natural frequencies; G _L(x) and φ _L(x) are known PDF and CF of the Cauchy–Lorentz distribution [54], see Eq. (12); delta functions in χ ̂ 0 ( ω ) are omitted for simplicity⁴.

Substituting the general form of the unbroadened susceptibilities χ i 0 ( . ) from Eq. (6), derived later in Section 3.1, into Eqs. (2) and (3) we obtain the fundamental dispersion relation ^[5]

where G i ( x ) = G i ( x ) + ı H G i ( x ) is the complex PDF incorporating the Hilbert transform (HT, H ) of G _i (x) as the imaginary part, and represents broadening, while [a _i , ϕ _i , Ω _i ] are the amplitude, phase and resonant frequency parameters of the ideal unperturbed transition (see Figure 2).

Figure 2:

Real and imaginary parts of susceptibility χ ̂ ( ω ) (or conductivity σ ̂ ( ω ) ) for different broadening functions: zero broadening (ZB), Lorentzian homogeneous broadening (HB) and Gaussian inhomogeneous broadening (IB) with different types of dispersion: (ab) oscillator, (cd) relaxation, (ef) conductive media, according to the newly derived formulas in this work.

Equation (5) represents a powerful theoretical framework that generates physically consistent permittivity models for any probability distribution with known complex PDFs G i ( x ) and CFs φ _i (t).^[6] In Section 3, we show how to use the general formula (5) for common broadening functions – Lorentz, Gauss, and mixed Gauss–Lorentz (Voigt), and different dispersion types – oscillator, relaxation and conductive media. A comprehensive summary of all cases and formulas, highlighting new (derived in this work) and established known models, is provided in Table A (Appendix A).

2.2 Analytical constraints

Time-domain modeling requires the dielectric function to be physically consistent, ensuring analyticity in the upper half-plane, causality, time-reversal symmetry (T-symmetry), Kramers–Kronig (KK) consistency, passivity and proper decay at infinity to satisfy the sum rule.

Causality of the total permittivity (ɛ(t) = 0, ∀t < 0) in Eq. (1a) is ensured as long as the unbroadened functions χ i 0 ( t ) are causal; e.g., general form (5a) is causal.

T-symmetry and KK-consistency. The real and imaginary parts of each term in (1b) satisfy the time-reversal symmetry χ ̂ i ( − ω ) = χ ̂ i * ( ω ) and are related via the Hilbert transform (HT, H ), ensuring KK consistency,

For symmetric distributions G _i (x) = G _i (−x), convolution (2) holds these properties, provided the unbroadened functions χ ̂ i 0 ( ω ) satisfy them; e.g., this holds in the general form (5b) since HT commutes with convolution and anticommutes with reflection implying G i ( − x ) = G i * ( x ) .

Sum rules. In ultrafast TD modeling, physically accurate high-frequency asymptotic behavior is important. As ω → ∞, total permittivity in Eq. (1a) should approach a free electron gas behavior: (a) ε ̂ ′ ( ω ) ∼ 1 − ω p 2 / ω 2 , with (b) the imaginary part decaying faster than 1/ω, ω ε ̂ ″ ( ω ) → 0 , [51], [55], [56].

Condition (b) yields the sum rule ∑ i a i ⁡ sin ϕ i = σ e ε 0 − 1 , requiring that all contributions to 1/ω from non-zero phase (ϕ _i ≠ 0) terms (e.g., conductivity, Debye, or phase-relaxed Lorentz) cancel out. This sum rule is unaffected by broadening and holds as long as satisfied for the unbroadened permittivity. Moreover, the MiMOSA approximation (Section 3.4) also conserves the sum rule (b) exactly, since ∑ i ∑ j a i j ⁡ sin ⁡ ϕ i j = ∑ i a i ⁡ sin ϕ i follows from combining conjugate pole pairs in Eq. (25) and constraint ∑ _j B ^j = π ^−1/2.

Condition (a) in Voigt multi-term dispersion model (16) is satisfied asymptotically, in both exact and MiMOSA models, as ε ′ ( ω ) − 1 = O ( ω − 2 ) , assuming ɛ _∞ = 1 (often relaxed over a finite frequency ranges). The exact constant (ω _p) is defined by the sum rule ω p 2 = ∑ i a i ( Ω i ⁡ cos ϕ i + γ i ⁡ sin ϕ i ) which, in general, can depend on homogeneous broadening γ _i (if non-zero phases ϕ _i are involved) but remains unaffected by inhomogeneous broadening σ _i . For zero-phase systems (∀i ϕ _i = 0), the MiMOSA approximation (Section 3.4) also preserves this sum rule exactly, ∑ i ∑ j a i j Ω i j ⁡ cos ⁡ ϕ i j + γ i j ⁡ sin ⁡ ϕ i j = ∑ _i a _i Ω _i , see Eq. (25) for validation.

Passivity of the total permittivity (ɛ″(ω) ≥ 0, ∀ω ≥ 0) is easy to ensure in the general formulation (5b) by the passivity of individual terms, provided that all phases are zero (ϕ _i = 0) and the broadening functions G _i (x) are bell-shaped.^[7] When non-zero phases (ϕ _i ≠ 0) are present, individual terms may locally exhibit gain, compensated by other terms in the total sum. A representative class of examples are MiMOSA models in Section 3.4, where coupled oscillators with conjugate poles maintain overall passivity^[8] (see also Figure 5 in [50]).

3.1 Zero broadening (ZB)

ZB represents an idealized scenario with infinitely narrow spectral lines (G(x) = δ(x)) and infinite transition lifetimes. In the class of rational functions, the general form of a single-term unbroadened model is derived by taking the limit γ → 0⁺ in the HB case (11) resulting in

with [a, ϕ, Ω] being the amplitude, phase and oscillation frequency parameters.

The ZB formula (6) is consistent with the fundamental dispersion equation (5), where a delta function distribution is used as the PDF,

and represents zero scattering γ = 0⁺.^[10]

The phase parameter ϕ in (6) (also known as the loss angle) mixes the real and imaginary parts and allows the transition between two orthogonal cases: (ϕ = 0, Ω > 0) representing a classical oscillator and ( ϕ = − π 2 , Ω = 0 ) corresponding to a relaxation in the time domain. The two cases of the ZB formula (6), along with a case for conductive media, are addressed below.

Lossless Lorentz oscillator (ϕ = 0, Ω > 0), also called the Sellmeier model [57], has quadratically decaying real part and delta functions in absorption

Lossless Debye relaxation ( ϕ = − π 2 , Ω = 0 ) , is simply a DC conductivity term [58]

Lossless Drude model ( ϕ = − π 2 , Ω = 0 , a = ε 0 ω p 2 , χ ( . ) → σ ( . ) ) is handled as a Debye case (9) but with a switch from susceptibility χ(.) to conductivity σ(.),^{[11] , [12]}

Here ω _p is a plasma frequency – a characteristic point where the lossless Drude permittivity ε ( ω ) = 1 − ω p 2 ω 2 switches from metallic to dielectric behavior [59].

The delta function terms in (6b), often omitted in the literature, represent degenerate distributions of zero width and play a key role in the convolution formalism. When the ZB model (6b) is convolved with a PDF G(x), the absorption of an oscillator (ϕ = 0) is directly linked to the function G(x) as

This is why, for example, a Gaussian distribution produces a Gaussian lineshape in the absorption. In the case of relaxation/conduction, the lineshape (of χ ̂ ( ω ) or σ ̂ ( ω ) ) is zero-centered and is rotated by ϕ = − π 2 from the imaginary part to the real part, e.g., χ ̂ ′ ( ω ) = a π G ( ω ) , as shown in Figure 2.

3.2 Homogeneous broadening (HB)

HB represents the natural linewidth broadening that affects all atoms or molecules equally, due to the finite lifetime τ = γ ⁻¹ of excited states (uncertainty principle [2]). The general form of single-term HB dispersion, also known as the critical point model [60], represents an arbitrary rational function^[13]

The HB case (11) can be derived by either convolving (“blurring”) the ideal unbroadened susceptibility χ ̂ 0 ( ω ) in (6) with the Cauchy–Lorentz PDF G _L(x), or substituting complex function G L ( x ) into the general formulation (5), where

As expected, the parameter substitutions (outlined in parenthesis) reduce the general HB formula (11) to the classical Lorentz [61], Debye [62], and Drude [59] dispersion models, as shown below.

Lorentz oscillator ϕ = 0 , Ω = ω 0 2 − Γ 2 4 , a = f Ω , γ = Γ 2 is conventionally formulated with doubled broadening Γ = 2γ, the natural frequency ω 0 = Ω 2 + γ 2 instead of resonance frequency Ω, and oscillator strength f = aΩ instead of amplitude a,

Debye relaxation ( ϕ = − π 2 , Ω = 0 , a = Δ ε τ , γ = 1 τ ) conventionally uses parameters of relaxation time τ = γ ⁻¹ and permittivity jump Δɛ = aτ,

Drude model ( ϕ = − π 2 , Ω = 0 , a = ε 0 ω p 2 , χ ( . ) → σ ( . ) ) is classically parameterized by the plasma frequency ω _p and broadening γ and is obtained as the Debye case of the conductivity function¹¹

In the Drude case, convolution with unbroadened susceptibility (10b) is unphysical but valid for its unbroadened conductivity function (10a), with χ(.) restored from σ(.) afterward.¹²

3.3 Inhomogeneous broadening (IB)

IB arises from statistical distribution of microscopic resonant frequencies Ω affected by local environmental variations and the Doppler shift.^[14] As a result, the observed spectral broadening deviates from the ideal Lorentzian lineshape to a mix of both – natural lifetime-based (HB) defined by γ and statistical (e.g., Gaussian) broadening defined by variance σ ² (Figure 1), leading to the general Gauss–Lorentz model

where w(z) is the Faddeeva (Kramp) function [63].

The IB formula (16) is obtained as a convolution of the unbroadened response χ ̂ 0 ( ω ) with both the Cauchy PDF (G _L, see Eq. (12)) and the Gaussian PDF (G _G) defined by

corresponding to the time-dependent scattering γ(t) = σ ² t/2. In the presence of multiple broadening mechanisms, the probability theory for the sum of random variables dictates that the PDFs are convolved, while their CFs are multiplied,^[15] and so the IB formula (16) can also be obtained from the general formula (5) using the Gauss–Lorentz (Voigt) PDF/CF

where the scattering function has both – the constant and the linear correction terms, γ(t) = γ + σ ² t/2.^[16]

The new Voigt formula is consistent with all limiting cases: σ → 0⁺ gives classical Lorentzian models (11), γ → 0⁺ gives pure Gaussian lineshape typical for strong disorder, while σ, γ → 0⁺ gives the ZB case (6). These transitions are easy to see through the general formula (5) and distributions equations (7), (12), (17) and (18).^[17]

As before, we simplify the general IB formula (16) for the oscillator, relaxation, and conductivity cases, specifying the corresponding parameter substitutions.

Gauss-Lorentz oscillator (ϕ = 0, Ω = ω 0 2 − Γ 2 4 , a = f Ω , γ = Γ 2 ) with strength f, natural frequency ω ₀ and two (HB and IB) broadening parameters Γ, σ reads

The first causal Gaussian oscillator model (Γ = 0⁺) was derived in 2006 [64] using a causality tip from [65]. The first attempts to formulate the Gauss–Lorentz model date back to the late 1970s [66], with a later reproduction [67] usually referred to as the Brendel-Bormann (BB)-model. Unfortunately, the non-causal BB model, incompatible with TD, remains broadly adopted by experimentalists to fit material responses to infra-red light, e.g., [68], [69], [70]. Causal corrections with logarithmic terms, rational approximations, and ongoing discussions of the physical validity of the BB model can be found in Refs. [71], [72], [73].

Our new Gauss–Lorentz (GL) model (19) fixes all the issues with the previous BB formulation [67] (see the details in Appendix B). When σ → 0⁺, the GL formula gives the classical Lorentz model with 1/(1 + x ²) absorption lineshape, (13). The case of Γ → 0⁺ gives a causal Gaussian oscillator with exp[−x ² ln2] absorption lineshape¹⁷ [64], while the mixed case (Γ, σ > 0) yields the Voigt profile – a new causal formulation with a time-dependent scattering function γ(t) = γ + σ ² t/2, first mentioned by Kim et al. [24], [25], and consistent with [74].

Gauss–Debye relaxation ( ϕ = − π 2 , Ω = 0, a = Δ ε τ , γ = 1 τ ) with the IB parameter σ additionally to the standard relaxation time τ and permittivity jump Δɛ reads

The limits σ → 0⁺ and γ → 0⁺ give the classical Debye (14) and a new Gauss relaxation model, respectively.¹⁷

The Gauss–Debye model (20) has not been shown in the literature. Known generalizations to the Debye relaxation – the Cole–Cole, Cole–Davidson, and Havriliak–Negami models [75], [76] – remain inaccessible to efficient TD simulations, and will be addressed in future work.

Gauss-Drude model ( ϕ = − π 2 , Ω = 0, a = ε 0 ω p 2 , χ(.) → σ(.)) with IB broadening parameter σ additionally to the plasma frequency ω _p and HB γ reads

The limits σ → 0⁺ and γ → 0⁺ give the classical Drude model (15b) and a new Gauss conductive model, respectively.¹⁷ The derivation is based on the broadening formalism (Section 2) for the conductivity function. Comparing classical Drude model (15b) to the new Gauss conductivity model (21b), we observe lineshape change and find that disordered analogue of classic conductivity ε 0 ω p 2 / γ , is ε 0 ω p 2 π / ( σ 2 ) . In the TD, linear argument of the exponential decay in σ(t) becomes quadratic, and in χ(t) changes to a complementary error function (erfc).

Corrections to the classical Drude model have been widely studied, including empirical frequency-domain formulations with fractional derivatives, effective mass parameter, and modified scattering functions [77]. However, the causal Gauss–Drude model introduced here has never been presented.

Figure 2 illustrates the final Voigt formula (16) for three cases of broadening – ZB, HB (σ = 0⁺) and IB (γ = 0⁺) for three types of dispersion: oscillator, relaxation, and conductive media. HB and IB curves are matched at the peak maximum and full-width-half-maximum (FWHM) of the real (for relaxation) or imaginary (for oscillator) parts, demonstrating the deviation of the “heavy-tail” Lorentzian lineshape 1 1 + x 2 from the Gaussian lineshape e − x 2 ln 2 . For conductive media, the DC conductivity ( σ ̂ (0)) and a plasma crossover point ω = ω p are matched.

These plots effectively illustrate the physical interpretation of the complex PDF G ( x ) = G ( x ) + ı H { G ( x ) } . For an oscillator with zero phase (ϕ = 0), the real part of G ( x ) defines the absorption spectrum ( χ ̂ ″ ( ω ) ) , producing symmetric peaks G(ω ± Ω), while the imaginary part generates KK-consistent terms H { G } ( ω ± Ω ) in χ ̂ ′ ( ω ) . When the phase ϕ is nonzero, the real and imaginary parts of the susceptibility become mixed. In the case of relaxation (ϕ = −π/2), χ ̂ ′ and χ ̂ ″ are effectively swapped, and the two resonant peaks coalesce into a single peak. The conductive case is identical to the relaxation case, with the susceptibility function χ ̂ ( ω ) replaced by the conductivity function σ ̂ ( ω ) .

3.4 Minimax approximation (MiMOSA)

When the broadening function G(x) is non-Lorentzian, the general dispersion formula (5) falls outside of the class of the rational functions of argument s = −ıω, and cannot be immediately translated into auxiliary differential equations, making it challenging to construct short discretization stencils and coupling to time-domain solvers, such as FDTD. The solution for efficient TD implementation of non-Lorentzian dispersion was first developed for a pure Gaussian oscillator [50], and employs minimax optimization to generate the shortest possible time stencil for a given error (MiMOSA).

Derivation of MiMOSA (Mini-max optimized semianalytical approximation) models for a general dispersion formula (5) starts with a minimax rational approximation^[18] of the complex PDF G x = G x + ı H G x ,

For example, for the Voigt distribution (18), the approximation coefficients [B ^j , C ^j ] are calculated for the Faddeeva function w(z), with n being the number of approximation poles, Figure 3(c),^[19]

Substituting the approximation (22) into the general susceptibility formula (5) gives a set of FDTD-compatible analytically derived dispersion terms

where parameters [a ^j , ϕ ^j , Ω ^j , γ ^j ] are directly connected to the parameters of the exact single-term model [a, ϕ, Ω, γ, σ] (5) and the approximation constants [B _j , C _j ] in (22) as

The MiMOSA model (24)–(25) retains the single-oscillator form [78], with only its envelope modified by approximation (compare to the exact susceptibility χ(t)),

This identity arises by substituting coefficients (25) into the time-domain expression (24) and combining the conjugate pole pairs. It ensures that the model remains physically consistent, without introducing nonphysical oscillations.

Due to the equioscillation theorem [79], the minimax solution provides the shortest rational polynomial approximation (corresponding to most compact numerical stencil), with the approximation error spread evenly across the entire frequency domain. In the Voigt case σ 2 π G ( z σ 2 ; γ = 0 ) = w ( z ) , the error converges exponentially with number of poles (n) throughout the upper half-plane [50]. As a result, even two poles (n = 2) already give a few percent error, which is sufficient for many applications, such as initial optimization or ellipsometry characterization. Using three poles (n = 3) drives the FD relative error below 1 %, making the approximation indiscernible from the experimental data (see Appendix C for details).

The MiMOSA models (24) fit the class of rational functions that can be coupled efficiently to the TD Maxwell’s solvers. Detailed ADE and RC numerical schemes for this class of dispersion can be found in Refs. [50], [80], [81] for second-order accurate TD solvers, and in Refs. [82], [83], [84] for higher-order schemes. We recommend using the universal compact scheme, which minimizes computational cost per dispersion term and enables easy switching between different second-order accurate ADE and RC formulations. An FDTD code implementing six such schemes is available in Ref. [50].

Compared to the approximations derived in the 1950s by reincarnating the minimax methods for rational polynomials and the more recent literature on the rational approximations to Faddeeva/Kramp/plasma dispersion function (or their real/imaginary parts) [85], [86], [87], [88], [89], [90], [91], [92], [93], [94], [95], [96], [97], [98], [99], [100], [101], [102], our MiMOSA method achieves impressive < 1 % error with just 2–3 terms (and thus the minimal number of additional equations in the numerical model), while preserving the necessary analytical properties of the dielectric function including causality and the sum rules.

A related computational approach has been recently proposed in Ref. [103], where an ab initio integral dispersion formulation is presented and subsequently transformed into a rational function through a quadrature approximation, thereby preserving the physical meaning of the main model parameters. However, no alternative approximation technique achieves the same minimal number of additional equations as MiMOSA for a given maximal error across the entire upper half-space.