Conventions

Defined exclusively in terms of quaternions, the Wigner $\mathfrak{D}$ matrices and spherical-harmonic functions (spin-weighted and scalar) are pretty simple, and easy to make internally consistent. However, it is important to establish which conventions are in use — especially in comparison to other source. Here, I carefully examine all the assumptions built in to the conventions for the spherical_functions package, and relate these choices to those made by other authors.

On this page, we simply introduce the basic conventions and notation for rotations in terms of the beautiful, elegant, efficient, and highly intuitive presentation with quaternions. We then compare to uglier, old-fashioned presentations in terms of spherical coordinates and the profoundly hideous Euler angles. On another page, we will see how to get directly from quaternions to highly efficient and accurate formulas for the Wigner $\mathfrak{D}$ matrices, with no use of Euler angles whatsoever. Similarly, we will be able to express spin-weighted spherical harmonics directly in terms of quaternions, though with a simple translation to and from standard spherical coordinates. This will allow us to derive simple rotation laws for the SWSHs and modes of a general decomposition in terms of SWSHs.

Quaternions, rotations, spherical coordinates

A unit quaternion (or “rotor”) $\rotor{R}$ can rotate a vector $\vec{v}$ into a new vector $\Rotated{\vec{v}}$ according to the formula \begin{equation*} \Rotated{\vec{v}} = \rotor{R}\, \vec{v}\, \rotor{R}^{-1}. \end{equation*} In principle, a unit quaternion obeys $\co{\rotor{R}} = \rotor{R}^{-1}$. In practice, however, there are cases where the system is (slightly slower, but) more stable numerically if the explicit inversion is used. And since the inversion is such a simple operation, we just use it.

The Wigner $\mathfrak{D}$ matrices can be derived directly in terms of components of the unit rotor corresponding to the rotation of which $\mathfrak{D}$ is a representation. That is, we don’t need to use Euler angles or rotation matrices at all. In fact, the resulting expressions are at least as simple as those using Euler angles, but are faster and more accurate to compute. The derivation is shown on this page.

When necessary to make contact with previous literature, we will use the standard spherical coordinates $(\vartheta, \varphi)$, in the physics convention. To explain this, we first establish our right-handed basis vectors: $(\basis{x}, \basis{y}, \basis{z})$, which remain fixed at all times. Thus, $\vartheta$ measures the angle from the positive $\basis{z}$ axis and $\varphi$ measures the angle of the projection into the $\basis{x}$-$\basis{y}$ plane from the positive $\basis{x}$ axis, as shown in the figure.

Now, if $\hat{n}$ is the unit vector in that direction, we construct a related rotor \begin{equation*} \rthetaphi = e^{\varphi \basis{z}/2}\, e^{\vartheta \basis{y}/2}. \end{equation*} This can be obtained as an quaternion module in python as

>>> import numpy as np, quaternion
>>> vartheta, varphi = 0.1, 0.2
>>> R_tp = quaternion.from_spherical_coords(vartheta, varphi)

Here, rotations are given assuming the right-hand screw rule, so that this corresponds to an initial rotation through the angle $\vartheta$ about $\basis{y}$, followed by a rotation through $\varphi$ about $\basis{z}$. (The factors of $1/2$ are present because $\rthetaphi$ is essentially the square-root of the rotation; the effect of this rotor produces the full rotations through $\vartheta$ and $\varphi$.) This rotor is particularly useful, because $\hat{n}$ and the two standard tangent vectors at that point are all given by rotations of the basis vectors $(\basis{x}, \basis{y}, \basis{z})$. Specifically, we have \begin{align*} \hat{\vartheta} &= \rthetaphi\,\basis{x}\,\rthetaphi^{-1}, \\\
\hat{\varphi} &= \rthetaphi\, \basis{y}\, \rthetaphi^{-1}, \\\
\hat{n} &= \rthetaphi\, \basis{z}\, \rthetaphi^{-1}. \end{align*} Note that we interpret the rotations in the active sense, where a point moves, while our coordinate system and basis vectors are left fixed. Hopefully, this won’t cause too much confusion if we express the original point with the vector $\basis{z}$, for example, and the rotated point is therefore given by some $\rotor{R}\, \basis{z}\, \rotor{R}$.

Euler angles

Euler angles are a terrible set of coordinates for the rotation group. Compared to the other three standard presentations of rotations (rotation matrices, axis-angle form, and the closely related unit quaternions), Euler angles present no advantages and many severe disadvantages. Composition of rotations is complicated, numerically slow and inaccurate, and essentially requires transformation to a different presentation anyway. Interpolation of Euler angles is meaningless and prone to severe distortions. Most damningly of all are the coordinate singularities (gimbal lock). To summarize, Euler angles are absolutely — and by a wide margin — the worst way to deal with rotations.

We can work entirely without Euler angles. (My own work simply never uses them; the $\mathfrak{D}$ matrices are written directly in terms of rotors.) Nonetheless, they are ubiquitous throughout the literature. And while we can work entirely without Euler angles, it can sometimes be useful to compare to other results. So, to make contact with that literature, we will need to choose a convention for constructing a rotation from a triple of angles $(\alpha, \beta, \gamma)$. We therefore define the rotor \begin{equation*} R_{(\alpha, \beta, \gamma)} = e^{\alpha\, \basis{z}/2}\, e^{\beta\, \basis{y}/2}\, e^{\gamma\, \basis{z}/2}. \end{equation*} This can be obtained as a quaternion object in python as

>>> import numpy as np, quaternion
>>> alpha, beta, gamma = 0.1, 0.2, 0.3
>>> R_euler = quaternion.from_euler_angles(alpha, beta, gamma)

Note that the rotations are always taken about the fixed axes $\basis{z}$ and $\basis{y}$. Also, this is in operator form, so it must be read from right to left: The rotation is given by an initial rotation through $\gamma$ about the $\basis{z}$ axis, followed by a rotation through $\beta$ about the $\basis{y}$ axis, followed by a final rotation through $\alpha$ about the $\basis{z}$ axis. This may seem slightly backwards, but it is a common convention — in particular, it is the one adopted by Wikipedia in its Wigner-D article.

It is worth noting that the standard right-handed basis vectors $(\basis{x}, \basis{y}, \basis{z})$ can be identified with generators of rotations usually seen in quantum mechanics (or generally just special-function theory) according to the rule \begin{align*} \frac{\basis{x}}{2} &\mapsto -i\, J_x, \\\
\frac{\basis{y}}{2} &\mapsto -i\, J_y, \\\
\frac{\basis{z}}{2} &\mapsto -i\, J_z. \end{align*} This is important when relating quaternion expressions to expressions more commonly seen in the literature. In particular, with this identification, we have the usual commutation relations \begin{align*} \left[\frac{\basis{x}}{2}, \frac{\basis{y}}{2}\right] = \frac{\basis{z}}{2} &\mapsto [J_x, J_y] = i\, J_z, \\\
\left[\frac{\basis{y}}{2}, \frac{\basis{z}}{2}\right] = \frac{\basis{x}}{2} &\mapsto [J_y, J_z] = i\, J_x, \\\
\left[\frac{\basis{z}}{2}, \frac{\basis{x}}{2}\right] = \frac{\basis{y}}{2} &\mapsto [J_z, J_x] = i\, J_y. \end{align*} And in any case, this certainly clarifies what to do with expressions like the following from Wikipedia: \begin{equation*} \mathcal{R}(\alpha,\beta,\gamma) = e^{-i\alpha\, J_z}\, e^{-i\beta\, J_y} e^{-i\gamma\, J_z}, \end{equation*} which shows that my interpretation of Euler angles is the same as Wikipedia’s.

Comparison with Mathematica

I find it nearly impossible to decipher the meaning behind Mathematica’s documentation (or anyone else’s for that matter), so the only comparison I’m willing to make is between actual results when using different codes. Specifically, the Mathematica expression

WignerD[{j, m1, m2}, psi, theta, phi]

results in a quantity that is identical to the output from my code with this expression:

(-1)**(m1+m2) * spherical_function.Wigner_D_element(quaternion.from_euler_angles(psi, theta, phi), j, m1, m2)

That is, if you chose actual values for j, m1, m2, psi, theta, phi, the numbers output by these two expressions would be identical (within numerical precision).

Note that besides the different syntax, the only real difference is the factor of $(-1)^{m_1+m_2}$. This is presumably related to the choice of Condon-Shortley phase.