# Numerical K-K — Practice *How the principal-value integrals actually get computed on measured, band-limited, noisy data — and what to do about the three unavoidable headaches: truncation, the singularity at ω′ = ω, and the need for subtractions when the response doesn't decay.* The convention choice affects signs only; algorithms below are written in the **physics convention** `e^{−iωt}` with UHP analyticity (Chapter 16). The K-K pair for a sufficiently decaying χ(ω) = χ′ + iχ″ is ``` χ′(ω) = (1/π) P ∫_{−∞}^{∞} χ″(ω′) / (ω′ − ω) dω′ χ″(ω) = −(1/π) P ∫_{−∞}^{∞} χ′(ω′) / (ω′ − ω) dω′ ``` ## The three problems 1. **Truncation.** Real measurement data covers a finite band [ω_min, ω_max], but the integral runs to infinity. 2. **Singularity.** The integrand has a 1/(ω′ − ω) pole at the evaluation point ω. 3. **Subtractions.** If χ(∞) ≠ 0 or χ diverges slowly at low frequency (DC conductivity, see Chapter 30), the unsubtracted form does not converge. The singly- or doubly-subtracted forms (Chapter 15) must be used. ## Dealing with truncation Three options, in order of increasing care: - **Zero extrapolation** — assume χ = 0 outside the measured band. Easy, but biases the reconstruction toward zero at the band edges. Usable only when the measured band genuinely contains the entire resonance structure. - **Match to known asymptotes.** Optical data must have n(ω → ∞) → 1 and ε(ω → ∞) → ε_∞ (Chapter 4). Dielectric data should match a Debye or Cole–Cole tail at low frequency. Fit the tails and use the analytic form as the integrand outside [ω_min, ω_max]. - **Lorentzian / Drude tails.** For broadband optical data, extend the measured χ″ with a Lorentzian oscillator fitted to the high-frequency tail and a Drude tail at low frequency for conductive media. The analytic K-K of each extension is known in closed form and can be added separately. The band-edge reconstruction is always the noisiest. A common diagnostic is to K-K-reconstruct χ′ from χ″, K-K again to recover χ″, and compare the round-trip residual. Large residuals at the band edges flag bad extrapolation. ## Dealing with the principal-value singularity Three standard numerical approaches: ### Maclaurin's method (skip the singular sample) Evaluate the integrand on a uniform grid ω′_k = ω_min + k Δω, and at the target frequency ω = ω_j the summand with k = j is **skipped** entirely: ``` χ′(ω_j) ≈ (2 Δω / π) · Σ_{k ≠ j, k+j odd} χ″(ω′_k) / (ω′_k − ω_j) ``` The `k + j odd` restriction is Maclaurin's trick: summing only over grid points whose index differs by an odd offset from j cancels the leading error of omitting the k = j term. The method is third-order accurate, trivial to code, and robust. ### Trapezoidal with midpoint shift Place the evaluation points ω_j at the midpoints of the integration grid ω′_k, so k = j is never reached. The singularity is avoided by construction, though the reconstructed χ′ is only on the shifted grid. ### FFT after sign-multiplier Interpolate measured data to a uniform ω grid, apply FFT, multiply by −j sgn(ω) (engineering) or +i sgn(ω) (physics), inverse FFT. The principal value is handled implicitly by the DFT periodisation. This is by far the fastest for large data and is the standard choice in ellipsometry and dielectric-spectroscopy software. The FFT method has one subtle pitfall: the data must be zero-padded enough that the implicit wrap-around does not alias the far tails into the signal band. A factor of 4× oversampling past the measured band is a safe default; more if the tails are heavy. ## Subtractions for bad asymptotic behaviour If χ(∞) = χ_∞ ≠ 0, use the singly-subtracted form (Chapter 15): ``` χ′(ω) − χ_∞ = (2/π) P ∫_{0}^{∞} [ω′ χ″(ω′) − ω χ″(ω)] / (ω′² − ω²) dω′ ``` This pulls an extra factor of ω into the denominator, improving high-frequency convergence. For χ with DC conductivity, separate the singular 1/ω piece first (see Chapter 30), K-K the remainder, and add the analytic contribution of the pole back in. ## Error estimation and iterative refinement The K-K map is a bounded operator on the appropriate Hardy space, so errors propagate linearly in a well-defined norm. Two practical estimators: - **Round-trip residual.** R(ω) = χ″_meas(ω) − K-K[K-K[χ″_meas]](ω). Plot and inspect near band edges. - **Cross-check with independent measurement.** Optical n and k measured independently by different techniques (transmission + reflection vs spectroscopic ellipsometry) should agree under K-K. Disagreement flags calibration issues. If the residual is large but follows a smooth pattern, iterative refinement of the tail fit often helps: re-fit the Lorentzian parameters to minimise the round-trip residual. ## Pseudocode: FFT–Hilbert K-K of a real sequence The simplest numerical K-K: given χ″ sampled on a uniform ω grid, reconstruct χ′ by computing the Hilbert transform of χ″ and negating. This assumes χ has been pre-subtracted so that χ(∞) = 0 and χ decays at both band edges. ```python import numpy as np def kk_imag_to_real(chi_imag, oversample=4): """ Reconstruct chi_real from chi_imag via FFT Hilbert transform. Physics convention: chi_real(w) = (1/pi) P int chi_imag(w')/(w'-w) dw'. chi_imag is assumed on a uniform grid covering the full support of chi. Oversampling factor pads the data before FFT to reduce wrap-around. """ N = len(chi_imag) M = oversample * N # zero-pad symmetrically pad = (M - N) // 2 y = np.zeros(M) y[pad : pad + N] = chi_imag # FFT, apply +i*sgn(w) multiplier (physics convention), # inverse FFT, take real part Y = np.fft.fft(y) k = np.fft.fftfreq(M) H = 1j * np.sign(k) chi_real_padded = np.real(np.fft.ifft(Y * H)) # strip the padding return chi_real_padded[pad : pad + N] def kk_real_to_imag(chi_real, oversample=4): """Inverse direction; note sign flip relative to imag→real.""" N = len(chi_real) M = oversample * N pad = (M - N) // 2 y = np.zeros(M) y[pad : pad + N] = chi_real Y = np.fft.fft(y) k = np.fft.fftfreq(M) H = -1j * np.sign(k) chi_imag_padded = np.real(np.fft.ifft(Y * H)) return chi_imag_padded[pad : pad + N] ``` For production code, add: - a tail-fit stage before the FFT, to extend chi_imag with Lorentzian (or other analytic) tails beyond the measured band; - a DC-subtraction stage to remove chi(∞) and any DC-conductivity pole before the Hilbert transform, re-added analytically afterward; - a windowing / apodisation stage near the band edges to suppress Gibbs ripple. ## When to use which method - **Fast, high-resolution data (thousands of points, broadband).** FFT method, with explicit tail matching and oversampling. - **Sparse, irregular data (tens of points, narrowband).** Maclaurin's method on the native grid, no interpolation. - **Data with strong DC pole.** Subtract the pole analytically, K-K the remainder by either method, add the pole back. - **Data with χ(∞) known non-zero.** Use the subtracted form (Chapter 15); the FFT method works unchanged on the subtracted quantity. ## See also - [15_subtracted_dispersion.md](15_subtracted_dispersion.md) - [28_discrete_hilbert_transform.md](28_discrete_hilbert_transform.md) - [30_pitfalls_limitations.md](30_pitfalls_limitations.md) ## References - K. Ohta and H. Ishida, "Comparison among several numerical integration methods for Kramers–Kronig transformation", *Applied Spectroscopy*, vol. 42, pp. 952–957, 1988. - F. W. King, *Hilbert Transforms*, vols 1–2, Cambridge University Press, 2009. - V. Lucarini, J. J. Saarinen, K.-E. Peiponen, E. M. Vartiainen, *Kramers–Kronig Relations in Optical Materials Research*, Springer, 2005 — Ch. 4 on numerics.