Analytical formula for numerical evaluations of the Wigner rotation matrices at high spins

The Wigner d function, which is the essential part of an irreducible representation of SU(2) and SO(3) parameterized with Euler angles, has been know to suffer from a serious numerical errors at high spins, if it is calculated by means of the Wigner formula as a polynomial of cos and sin of half of the second Euler angle. This paper shows a way to avoid this problem by expressing the d functions as the Fourier series of the half angle. A precise numerical table of the coefficients of the series is obtainable from a web site.

💡 Research Summary

The paper addresses a long‑standing numerical problem in the evaluation of the Wigner small‑d matrix elements, which are the essential building blocks of SU(2) and SO(3) rotation operators. The conventional expression, often called the Wigner formula, writes d⁽ʲ⁾ₘₖ(θ) as a finite polynomial in cos(θ/2) and sin(θ/2). While analytically exact, this representation suffers from catastrophic loss of significance when the angular‑momentum quantum number j becomes large (typically j ≳ 50). The individual terms in the polynomial scale as 2ʲ, so that for double‑precision arithmetic (53‑bit mantissa) the cancellation among huge terms destroys all meaningful digits already at j≈54, and even quadruple precision (113‑bit mantissa) fails beyond j≈114. This loss of precision is especially severe near θ = π/2 with m = k = 0, a regime that appears frequently in nuclear‑structure calculations, deformed‑mean‑field projections, and other high‑spin applications.

To overcome this, the author proposes a completely different representation: a Fourier series in the half‑angle. By exploiting the parity of m − k, the d‑function can be written as a sum over ν of either cos νθ or sin νθ, with ν ranging from a minimum value ν_min (determined by the parity of m − k and by j) up to j. The coefficients t^{ν}{j m k} are obtained by exploiting the orthogonality of the trigonometric basis: t^{ν}{j m k}= (1/2π)(1+δ_{ν0})∫₀^{4π} d^{j}{m k}(θ) f(νθ) dθ, where f is cos for even (m − k) and sin for odd (m − k). Substituting the Wigner polynomial into this integral yields a compact closed‑form expression (Eq. 16) that involves only factorials, binomial coefficients, and simple integrals I{λμ}=∫₀^{2π}cos λx sin μx dx. Crucially, the I_{λμ} vanish unless both λ and μ are even, which eliminates many terms and guarantees that each term in the series is of order unity or smaller. Consequently, the Fourier series does not suffer from the exponential growth that plagues the original polynomial.

The coefficients t^{ν}_{j m k} are still non‑trivial to compute for large j because intermediate quantities can be astronomically large. The author therefore evaluates them with arbitrary‑precision arithmetic using the computer algebra system MAXIMA, retaining enough digits to guarantee that the final 64‑bit floating‑point representation is accurate to the full mantissa. For j up to 100, the required precision reaches about 74 decimal digits. The resulting 64‑bit coefficients are then stored in plain‑text files; the total storage requirement grows roughly as j⁴, amounting to ~27 MiB for j_max = 50 and ~40 MiB for j_max = 100. The magnitude of the coefficients decreases with j, but even at j = 100 most coefficients are larger than 10⁻⁵, ensuring that they contribute meaningfully to the sum.

Evaluation of the trigonometric functions cos νθ and sin νθ is performed via a stable recurrence relation (Eq. 27) rather than direct calls to library functions. The recurrence keeps the arguments within