NOTE: I've run this notebook with the correction to sympy's trig simplification routines found in this pull request, which has not yet made it into a released version. I just ran python -c "import sympy; print(sympy.__file__)" to find where on my system the actual files are, then edited .../simplify/fu.py as given in the PR.
In [1]:
import csv
import sympy
from sympy import sin, cos
from sympy.parsing.mathematica import mathematica
from sympy.physics.quantum.spin import Rotation
from sympy.abc import _clash
import numpy as np
import quaternion
import spherical_functions as sf
The Mathematica documentation for WignerD says
WignerD[{j, m_1, m_2}], psi, theta, phi] gives the Wigner D-function $D^j_{m_1, m_2}(\psi, \theta, \phi)$.The Wigner D-function $D^j_{m_1, m_2}$ gives the matrix element of a rotation operator parametrized by Euler angles in a $2 j+1$-dimensional unitary representation of a rotation group when parameters $j, m_1, m_2$ are physical, i.e. all integers or half-integers such that $-j \leq m_1,m_2 \leq j$.
The Wolfram Language uses phase conventions where $D^j_{m_1, m_2}(\psi, \theta, \phi) = e^{i m_1 \psi + i m_2 \phi} D^j_{m_1, m_2}(0, \theta, 0)$.
WignerD[{j, m1, m2}, psi, theta, phi] == (-1)^(m1 - m2) Conjugate[WignerD[{j, -m1, -m2}, psi, theta, phi]
WignerD[{j, m1, m2}, psi, theta, phi] == (-1)^(m1 - m2) WignerD[{j, m2, m1}, phi, theta, psi]
There are no more specifics about what the Euler angles mean in this function's documentation, but the documentation for EulerMatrix[{alpha, beta, gamma}] says that it "gives the Euler 3D rotation matrix formed by rotating by $\alpha$ around the current $z$ axis, then by $\beta$ around the current $y$ axis, and then by $\gamma$ around the current $z$ axis." This is ambiguous, but
we later find that EulerMatrix[{alpha, beta, gamma}, {a, b, c}] is equivalent to $R_{\alpha, a} R_{\beta, b} R_{\gamma, c}$. Evidently, "current" refers to the rotating body axes, and so EulerMatrix[{alpha, beta, gamma}] is what I would write in quaternion form as
\begin{equation}
e^{\alpha \hat{z}/2} e^{\beta \hat{y}/2} e^{\gamma \hat{z}/2} = e^{\gamma \hat{z}''/2} e^{\beta \hat{y}'/2} e^{\alpha \hat{z}/2}
\end{equation}
I've created a CSV file with the analytic expressions for j from 0 through 5, using this code:
SetDirectory[NotebookDirectory[]];
Export[
"conventions_mathematica.csv",
Flatten[
Table[{j, m1, m2, ToString[WignerD[{j, m1, m2}, psi, theta, phi], InputForm]},
{j, 0, 5}, {m1, -j, j}, {m2, -j, j} ], 2],
TableHeadings -> {"j", "m1", "m2", "WignerD[{j, m1, m2}, psi, theta, phi]"}
];
I'll be comparing these expressions to SymPy's, and then evaluating them to compare to the results from spherical_functions.
The SymPy documentation is unclear and a little self-contradictory. The main docstring for sympy.physics.quantum.spin.Rotation says that it
Defines the rotation operator in terms of the Euler angles defined by the z-y-z convention for a passive transformation. That is the coordinate axes are rotated first about the z-axis, giving the new x'-y'-z' axes. Then this new coordinate system is rotated about the new y'-axis, giving new x''-y''-z'' axes. Then this new coordinate system is rotated about the z''-axis. Conventions follow those laid out in Varshalovich.
alpha: First Euler Anglebeta: Second Euler anglegamma: Third Euler angle
The docstring for sympy.physics.quantum.spin.WignerD says
The Wigner D-function gives the matrix elements of the rotation operator in the jm-representation. For the Euler angles $\alpha$, $\beta$, $\gamma$, the D-function is defined such that: \begin{equation} \left\langle j,m| \mathcal{R}(\alpha, \beta, \gamma ) |j',m' \right \rangle = \delta_{jj'} D(j, m, m', \alpha, \beta, \gamma) \end{equation} Where the rotation operator is as defined by the Rotation class.
The Wigner D-function defined in this way gives: \begin{equation} D(j, m, m', \alpha, \beta, \gamma) = e^{-i m \alpha} d(j, m, m', \beta) e^{-i m' \gamma} \end{equation} Where
dis the Wigner small-d function, which is given byRotation.d.The Wigner small-d function gives the component of the Wigner D-function that is determined by the second Euler angle. That is the Wigner D-function is: \begin{equation} D(j, m, m', \alpha, \beta, \gamma) = e^{-i m \alpha} d(j, m, m', \beta) e^{-i m' \gamma} \end{equation} Where
dis the small-d function. The Wigner D-function is given byRotation.D.
j: Total angular momentumm: Eigenvalue of angular momentum along axis after rotationmp: Eigenvalue of angular momentum along rotated axis
Again, this is all pretty ambiguous, regarding exactly which angle is supposed to go with which rotation, but my best guess is that it looks like this: \begin{equation} e^{\gamma \hat{z}''/2} e^{\beta \hat{y}'/2} e^{\alpha \hat{z}/2} = e^{\alpha \hat{z}/2} e^{\beta \hat{y}/2} e^{\gamma \hat{z}/2}, \end{equation} which is precisely the same as my interpretation of Mathematica's convention.
Here's the notation SymPy uses for the D matrices:
In [10]:
Rotation.D(j, m, mp, alpha, beta, gamma)
Out[10]:
Although the Mathematica documentation doesn't explicitly relate its WignerD and EulerMatrix functions, I think enough of Mathematica to guess that they at least use consistent conventions. And spherical_functions explicitly takes a quaternion object, so to the extent that I use Euler angles at all, we can stay consistent in this way.
So first, we check the rotation matrix that comes out of quaternion via Euler angles:
In [49]:
class QuaternionFromEuler(object):
def __init__(self, alpha, beta, gamma):
# This is essentially copied from the quaternion code
self.w = cos(beta/2)*cos((alpha+gamma)/2)
self.x = -sin(beta/2)*sin((alpha-gamma)/2)
self.y = sin(beta/2)*cos((alpha-gamma)/2)
self.z = cos(beta/2)*sin((alpha+gamma)/2)
q = QuaternionFromEuler(psi, theta, phi)
sympy.Matrix([
[simplify(1 - 2*(q.y**2 + q.z**2)), simplify(2*(q.x*q.y - q.z*q.w)), simplify(2*(q.x*q.z + q.y*q.w))],
[simplify(2*(q.x*q.y + q.z*q.w)), simplify(1 - 2*(q.x**2 + q.z**2)), simplify(2*(q.y*q.z - q.x*q.w))],
[simplify(2*(q.x*q.z - q.y*q.w)), simplify(2*(q.y*q.z + q.x*q.w)), simplify(1 - 2*(q.x**2 + q.y**2))]
])
Out[49]:
This is precisely the same matrix as Mathematica returns from EulerMatrix[{ψ, θ, ϕ}], which would suggest to me that my Euler conventions are the same as Mathematica's.
In [2]:
j, m1, m2 = sympy.symbols('j, m1, m2', integer=True)
#psi, theta, phi = sympy.symbols('psi, theta, phi', real=True)
In [3]:
with open('conventions_mathematica.csv', 'r') as csvfile:
reader = csv.reader(csvfile)
header = next(reader, None)
print('Mathematica header:', ', '.join(header))
mathematica_wignerD = {
tuple(int(s) for s in row[:3]): sympy.sympify(mathematica(row[3]))
for row in reader
}
Unfortunately, I can't get sympify to correctly use locals, so I have to just grab all the symbols that it created in the previous cell, as follows:
In [4]:
free_symbols = list(set(symbol for jm1m2 in mathematica_wignerD for symbol in mathematica_wignerD[jm1m2].free_symbols))
In [7]:
sorted(free_symbols, key=lambda s: str(s))
Out[7]:
In [8]:
phi, psi, theta = sorted(free_symbols, key=lambda s: str(s))
In [9]:
j, m, mp = sympy.symbols('j, m, mprime', integer=True)
alpha, beta, gamma = sympy.symbols('alpha, beta, gamma', real=True)
As always, SymPy isn't good enough at simplifying trig functions, so I have to jump through some extra hoops:
In [11]:
half_angle_replacements = ([]
+[(sin(theta/2)**n, ((1-cos(theta))/2)**(n//2)) for n in [2, 4, 6, 8]]
+[(cos(theta/2)**n, ((1+cos(theta))/2)**(n//2)) for n in [2, 4, 6, 8]]
+[(sin(theta/2)*cos(theta/2), sin(theta)/2)]
#+[(sin(3*theta/2), 3*sin(theta/2)-4*sin(theta/2)**3)]
#+[(sin(5*theta/2), 5*cos(theta/2)**4*sin(theta/2)-10*cos(theta/2)**2*sin(theta/2)**3+sin(theta/2)**5)]
)
def simplify(difference):
from sympy import trigsimp, expand
difference = trigsimp(expand(sympy.simplify(difference), trig=True))
difference = sympy.simplify(difference.subs(half_angle_replacements, simultaneous=True))
return difference
Finally, we can go through and check each and every expression (even though we could have assumed some symmetries to skip certain combinations) to ensure that the Mathematica expression returned by
WignerD[{j, m1, m2}, psi, theta, phi]is identical to the SymPy expression returned by
Rotation.D(j, m1, m2, -psi, -theta, -phi)
In [12]:
for j, m1, m2 in mathematica_wignerD:
mathematica_value = sympy.expand(sympy.trigsimp(sympy.simplify(mathematica_wignerD[(j, m1, m2)])), trig=True)
sympy_value = sympy.expand(sympy.trigsimp(sympy.simplify(Rotation.D(j, m1, m2, -psi, -theta, -phi).doit())), trig=True)
ratio = sympy.simplify(mathematica_value/sympy_value)
mathematica_value, sympy_value = sympy.fraction(ratio)
mathematica_value = sympy.simplify(sympy.simplify(mathematica_value.subs(half_angle_replacements)).subs(half_angle_replacements))
sympy_value = sympy.simplify(sympy.simplify(sympy_value.subs(half_angle_replacements)).subs(half_angle_replacements))
difference = simplify(mathematica_value - sympy_value)
print('Checking (j, m1, m2) = ({0}, {1}, {2})'.format(j, m1, m2))
if difference:
display(mathematica_value, sympy_value, difference)
print()
All cases show agreement.
I'm not quite sure how to interpret this weird sign difference. Flipping the signs is one of the things you do when inverting a rotation, but you also flip the order of the angles. This could essentially be done here as well if we also flip the order of m1 and m2 — except that we need an additional factor of $(-1)^{m_1+m_2}$. So we could think of this as saying that one of these provides the D matrix for the inverse rotation of the other, and they swap the order of the m arguments, and there's a (Condon-Shortley) phase difference.
I'll just quickly verify that Rotation.D actually satisfies this symmetry:
In [13]:
for j in range(6):
for m1 in range(-j, j+1):
for m2 in range(-j, j+1):
#print(j, m1, m2)
difference = sympy.simplify(Rotation.D(j, m1, m2, -psi, -theta, -phi).doit()
- (-1)**(m1+m2)*Rotation.D(j, m2, m1, -phi, -theta, -psi).doit())
if difference:
display(difference)
So another way to say this is that SymPy takes the inverse rotation, and returns the transpose with that weird phase.
I find both Mathematica's and SymPy's descriptions to be ambiguous, but it looks like Mathematica's is closer to my thinking — except that other places in the documentation make me think that their Condon-Shortley phases are weird, so I'll play around with that until I get some agreement. I'll check by simply evaluating on random numbers.
In [23]:
np.random.seed(1234)
for _ in range(100): # Test for 100 sets of random Euler angles
ψ, θ, ϕ = np.random.rand(3) * np.array([2*np.pi, np.pi, 2*np.pi])
for j, m1, m2 in mathematica_wignerD:
mathematica_value = mathematica_wignerD[(j, m1, m2)].subs({psi: ψ, theta: θ, phi: ϕ}).evalf()
spherical_functions_value = (-1)**(m1+m2) * sf.Wigner_D_element(quaternion.from_euler_angles(ψ, θ, ϕ), j, m1, m2)
diff = abs(mathematica_value - spherical_functions_value)
if diff > 3e-13:
print(j, m1, m2, ψ, θ, ϕ, diff)
As expected, the only difference is in the overall phase. In short, we have numerical equality between Mathematica's expression
WignerD[{j, m1, m2}, psi, theta, phi]and spherical_function's expression
(-1)**(m1+m2) * sf.Wigner_D_element(quaternion.from_euler_angles(psi, theta, phi), j, m1, m2)So, the main conclusion is this:
The Wigner $\mathfrak{D}$ functions contained in Mathematica and
spherical_functionsagree except for an overall phase factor of $(-1)^{m_1 + m_2}$.