Boyle (2016)

Summary

The Wigner $𝔇$ matrices of Boyle (2016) — which were the convention of this package before version 3.0 — are the complex conjugates of the ones used in the SphericalFunctions package, for integer and half-integer indices: $𝔇^{(ℓ)}_{m',m}(𝐑)|_{2016} = \overline{𝔇^{(ℓ)}_{m',m}(𝐑)}$. Because that paper defines the spin-weighted spherical harmonics as $(-1)^s \sqrt{(2ℓ+1)/4π}\, 𝔇^{(ℓ)}_{m,-s}(𝐑)|_{2016}$, its ${}_sY_{ℓ,m}$ agree with ours. The explicit expression for $𝔇$ in Eq. (35) of that paper contains two errors, which are described and corrected below; one of them affects only half-integer $ℓ$.

Boyle [1] argued that spin-weighted spherical functions should be defined as functions on the spin group $\mathrm{Spin}(3)$, represented by unit quaternions, rather than on coordinates of the 2-sphere. That paper is the origin of most of the conventions in this package, with one important exception: the convention for the $𝔇$ matrices themselves is the complex conjugate of the paper's — so that they agree with LALSuite, Wikipedia, Sakurai, and the other sources in our priority list — as explained on the conventions pages.

The paper defines the spin-weighted spherical harmonics as functions of a unit quaternion $𝐑$ by [Eq. (21)]

\[{}_sY_{ℓ,m}(𝐑) = (-1)^s \sqrt{\frac{2ℓ+1}{4π}}\; 𝔇^{(ℓ)}_{m,-s}(𝐑),\]

and gives an explicit expression for $𝔇$ directly in terms of the components of the quaternion [Eq. (35)]. Writing $𝐑 = R_s + R_a$ with $R_s = W + Z𝐤$ and $R_a = Y𝐣 + X𝐢$ (the parts that commute and anticommute with $𝐤$), and treating these as complex numbers $R_s = r_s e^{iϕ_s}$ and $R_a = r_a e^{iϕ_a}$, the expression is a sum over $ρ$ of terms proportional to $r_s^{2ℓ-(m-m')-2ρ}\, r_a^{(m-m')+2ρ}\, e^{i[(m+m')ϕ_s + (m-m')ϕ_a]}$, with an alternative form when $r_a > r_s$; Appendix A of the paper describes how to evaluate the sum stably, and that algorithm is transcribed below as WignerDElement. Since $e^{iϕ_s}$ and $e^{iϕ_a}$ are the half-angle phases $e^{i(α+γ)/2}$ and $e^{i(γ-α)/2}$, the $e^{+i(m'α + mγ)}$ dependence is the complex conjugate of ours.

Errata in Eq. (35) of Boyle (2016)

Eq. (35) of the paper, as printed, contains two errors. The code below implements the corrected expression, which agrees with the naive sum of Eq. (33), from which Eq. (35) was derived, and with every other source on these pages.

In the two branches for $r_s < ϵ$ and $r_s < r_a$, the sign $(-1)^{ℓ+m}$ should be $(-1)^{ℓ-m}$. The two are equal for integer $ℓ$, but opposite for half-integer $ℓ$ — for which the paper states that the expression is valid. The correct sign is that of the only term of Eq. (33) that survives as $r_s → 0$, namely $ρ = ℓ - m$.
In the two general branches, the powers of $r_s$ and $r_a$ in front of each sum already include the first term of the sum, at $ρ = ρ_1$ or $ρ = ρ_3$, so that the summands should be $(-1)^ρ (r_a^2 / r_s^2)^{ρ-ρ_1}$ and $(-1)^ρ (r_s^2 / r_a^2)^{ρ-ρ_3}$, rather than $(-r_a^2 / r_s^2)^ρ$ and $(-r_s^2 / r_a^2)^ρ$. As printed, the expression is wrong — for integer $ℓ$ as well — whenever $m' > m$ in the branch for $r_a ≤ r_s$, or $m' < -m$ in the branch for $r_s < r_a$.

The tests below confirm the corrected expression by comparing WignerDElement with this package for integer and half-integer indices, and show explicitly that the printed sign fails for half-integer $ℓ$.

The same paper defines left and right operators $L$ and $K$ [Eqs. (42)–(43)] with $K_z\, {}_sY_{ℓ,m} = -s\, {}_sY_{ℓ,m}$ [Eq. (47)], and identifies $\eth = -K_-$ and $\bar{\eth} = K_+$ [Eq. (46)]; in the present conventions $K = -R$, so this is the statement $R_z\, {}_sY_{ℓ,m} = s\, {}_sY_{ℓ,m}$, $\eth = R_+$, and $\bar{\eth} = -R_-$ of the summary page.

Implementing formulas

Because the WignerDElement function is also needed on the Varshalovich page — where it serves as one of two independent references for half-integer $ℓ$ — we define it in a test module that both pages can use.

using Quaternionic

Compute a single Wigner-D matrix element for half-integer or integer $(ℓ, m', m)$, following Eq. (35) and Appendix A of Boyle (2016), with the two corrections to Eq. (35) described above: the sign $(-1)^{ℓ-m}$ in the branches for $r_s < ϵ$ and $r_s < r_a$, and summands measured from the first term of each sum. R is a Rotor, and ℓ, m′, and m are the indices of the Wigner-D matrix element. The indices must all be integers or all be Rational with denominators of 2.

function WignerDElement(R::Rotor{T}, ℓ::I, m′::I, m::I) where {T, I}    # If `I` is Rational, check that the denominators are 2    if I <: Rational        if (denominator(ℓ) != 2) || (denominator(m′) != 2) || (denominator(m) != 2)            error("The indices ℓ, m′, and m must all be integers or all be half-integers")        end    end    # Convert to twice the input values for half-integer support    L  = Int(2ℓ)    M′ = Int(2m′)    M  = Int(2m)    if L > 16        error(            "The maximum supported ℓ for this function is 8; " *            "larger numbers become numerically unstable.\n" *            "Consider using `SphericalFunctions.D` instead."        )    end    # Simple helper for 0+0im    zeroCT = zero(Complex{T})    if abs(M′) > L || abs(M) > L        return zeroCT    end    let π = T(π)        # Split input `R` into its two complex components and extract magnitude and phase        Rₛ = Complex(R[1], R[4])        Rₐ = Complex(R[3], R[2])        rₛ = abs(Rₛ)        rₐ = abs(Rₐ)        ϕₛ = angle(Rₛ)        ϕₐ = angle(Rₐ)        # Check simple limiting cases        if rₐ ≤ 4eps(rₛ)            if M′ != M                return zeroCT            else                return cis(M * ϕₛ)            end        elseif rₛ ≤ 4eps(rₐ)            if -M′ != M                return zeroCT            else                return cis(M * ϕₐ) * (((L - M) % 4 == 0) ? 1 : -1)            end        elseif rₐ ≤ rₛ            λ = -(rₐ/rₛ)^2            ρₘᵢₙ = max(0, (M′ - M)÷2)            κ = √(                    T(factorial((L + M)÷2) * factorial((L - M)÷2))                    / T(factorial((L + M′)÷2) * factorial((L - M′)÷2))                ) *                binomial((L + M′)÷2, ρₘᵢₙ) * binomial((L - M′)÷2, (L - M)÷2 - ρₘᵢₙ)            if (ρₘᵢₙ % 2) != 0                κ = -κ            end            ρₘₐₓ = min((L + M′)÷2, (L - M)÷2)            N₁ = L + M′ + 2            N₂ = L - M + 2            N₃  = M - M′            total = one(T)            for P in reverse(2ρₘᵢₙ+2:2:2ρₘₐₓ)                total *= λ * ((N₁ - P)*(N₂ - P)) / (P*(N₃ + P))                total += one(T)            end            return κ *                (rₛ ^ (L - (M - M′)÷2 - 2ρₘᵢₙ)) *                (rₐ ^ ((M - M′)÷2 + 2ρₘᵢₙ)) *                cis((M + M′)÷2 * ϕₛ + (M - M′)÷2 * ϕₐ) *                total        else # rₛ < rₐ            λ = -(rₛ/rₐ)^2            ρₘᵢₙ = max(0, -(M′ + M)÷2)            κ = √(                    T(factorial((L + M)÷2) * factorial((L - M)÷2))                    / T(factorial((L + M′)÷2) * factorial((L - M′)÷2))                ) *                binomial((L + M′)÷2, (L - M)÷2 - ρₘᵢₙ) * binomial((L - M′)÷2, ρₘᵢₙ)            if (((L - M)÷2 - ρₘᵢₙ) % 2) != 0                κ = -κ            end            ρₘₐₓ = min((L - M′)÷2, (L - M)÷2)            N₁ = L - M′ + 2            N₂ = L - M + 2            N₃  = M + M′            total = one(T)            for P in reverse(2ρₘᵢₙ+2:2:2ρₘₐₓ)                total *= λ * ((N₁ - P)*(N₂ - P)) / (P*(N₃ + P))                total += one(T)            end            return κ *                (rₐ ^ (L - (M + M′)÷2 - 2ρₘᵢₙ)) *                (rₛ ^ ((M + M′)÷2 + 2ρₘᵢₙ)) *                cis((M - M′)÷2 * ϕₐ + (M + M′)÷2 * ϕₛ) *                total        end    endend

Equation (21), the spin-weighted spherical harmonics as functions of a rotor:

function ₛYₗₘ(s, ℓ, m, R::Rotor{T}) where {T}    (-1)^s * √((2ℓ+1) / (4T(π))) * WignerDElement(R, ℓ, m, -s)end

Tests

We can now test the functions against the equivalent functions from the SphericalFunctions package.

using Quaternionic

We will need to test approximate floating-point equality, so we set absolute and relative tolerances (respectively) in terms of the machine epsilon:

ϵₐ = 100eps()ϵᵣ = 1000eps()

The algorithm is accurate up to

ℓₘₐₓ = 8

so we test up to that point, on a set of rotors that includes the identity, the basis rotations, rotors near the special cases of the algorithm, and random rotors. The special cases apply when $r_a$ or $r_s$ is below $4ϵ$, so the rotors near them are offset by $3ϵ$, inside those thresholds, and by $6ϵ$, which puts most of them just outside, in the general branches:

Rs = [    Rotor{Float64}(1);    [Rotor{Float64}(𝐯) for 𝐯 ∈ (imx, imy, imz)];    [exp(3eps() * 𝐯) for 𝐯 ∈ (imx, imy, imz)];    [Rotor{Float64}(𝐮) * exp(3eps() * 𝐯) for 𝐮 ∈ (imx, imy, imz) for 𝐯 ∈ (imx, imy, imz)];    [exp(6eps() * 𝐯) for 𝐯 ∈ (imx, imy, imz)];    [Rotor{Float64}(𝐮) * exp(6eps() * 𝐯) for 𝐮 ∈ (imx, imy, imz) for 𝐯 ∈ (imx, imy, imz)];    randn(rng, Rotor{Float64}, 20)]

For integer indices, the 2016 matrices are the complex conjugates of ours:

for R ∈ Rs    for (ℓ, 𝔇ˡ) ∈ SphericalFunctions.DCalculator(R, ℓₘₐₓ), m′ ∈ -ℓ:ℓ, m ∈ -ℓ:ℓ        @test Boyle2016.WignerDElement(R, ℓ, m′, m) ≈ conj(𝔇ˡ[m′, m]) atol=ϵₐ rtol=ϵᵣ    endend

Because Eq. (21) uses $𝔇_{m,-s}$ without a conjugate — where our definition uses $\overline{𝔇_{m,-s}}$ — the spin-weighted spherical harmonics agree:

import Quaternionic: from_spherical_coordinatesfor (θ, ϕ) ∈ θϕrange(rng, Float64, 7)    R = from_spherical_coordinates(θ, ϕ)    for (ℓ, Yˡ) ∈ SphericalFunctions.sYlmCalculator(θ, ϕ, 4, -2:2)        for s ∈ -min(ℓ, 2):min(ℓ, 2), m ∈ -ℓ:ℓ            @test Boyle2016.ₛYₗₘ(s, ℓ, m, R) ≈ Yˡ[s, m] atol=ϵₐ rtol=ϵᵣ        end    endend

For half-integer indices, the 2016 matrices are also the complex conjugates of ours. Our calculator labels half-integer blocks with HalfOddIntegers, which serve as indices but refuse arithmetic with other kinds of number; the 2016 formulas do that arithmetic, so the label is converted to a Rational first.

for R ∈ Rs    for (J, 𝔇ᴶ) ∈ SphericalFunctions.DCalculator(R, 15//2)        J = Rational(J)        for M′ ∈ -J:J, M ∈ -J:J            @test Boyle2016.WignerDElement(R, J, M′, M) ≈ conj(𝔇ᴶ[M′, M]) atol=ϵₐ rtol=ϵᵣ        end    endend

The sign printed in the $r_s < ϵ$ branch of Eq. (35), $(-1)^{ℓ+m}$, is wrong for half-integer indices, as described in the errata above. On rotors with $r_s = 0$, the nonzero elements are $±e^{2imϕ_a}$, and the printed sign gives the negative of the correct value for every one of them:

for R ∈ (Rotor{Float64}(imx), Rotor{Float64}(imy), Rotor(0.0, 0.6, 0.8, 0.0))    ϕₐ = angle(Complex(R[3], R[2]))    for (J, 𝔇ᴶ) ∈ SphericalFunctions.DCalculator(R, 15//2)        J = Rational(J)        denominator(J) == 2 || continue        for M ∈ -J:J            printed = (-1)^Int(J + M) * cis(2M * ϕₐ)            @test printed ≈ -conj(𝔇ᴶ[-M, M]) atol=ϵₐ rtol=ϵᵣ        end    endend

They also satisfy two properties that any representation of $\mathrm{Spin}(3)$ must satisfy — the representation property $𝔇(𝐑_1 𝐑_2) = 𝔇(𝐑_1)\, 𝔇(𝐑_2)$ and the sign change $𝔇(-𝐑) = -𝔇(𝐑)$ for half-integer $ℓ$. The comparison against an independent closed form is on the Varshalovich page.

for R₁ ∈ Rs[1:8], R₂ ∈ Rs[end-4:end]    for J ∈ (1//2, 3//2, 5//2, 7//2)        𝔇₁ = [Boyle2016.WignerDElement(R₁, J, M′, M) for M′ ∈ -J:J, M ∈ -J:J]        𝔇₂ = [Boyle2016.WignerDElement(R₂, J, M′, M) for M′ ∈ -J:J, M ∈ -J:J]        𝔇₁₂ = [Boyle2016.WignerDElement(R₁ * R₂, J, M′, M) for M′ ∈ -J:J, M ∈ -J:J]        @test 𝔇₁₂ ≈ 𝔇₁ * 𝔇₂ atol=ϵₐ rtol=ϵᵣ        𝔇₋ = [Boyle2016.WignerDElement(-R₁, J, M′, M) for M′ ∈ -J:J, M ∈ -J:J]        @test 𝔇₋ ≈ -𝔇₁ atol=ϵₐ rtol=ϵᵣ    endend

These successful tests show that the $𝔇$ matrices of Boyle (2016) are the complex conjugates of those defined by the SphericalFunctions package, for integer and half-integer indices, and that the spin-weighted spherical harmonics agree.

This page was generated using Literate.jl.