Python/numba package for evaluating and transforming Wigner $\mathfrak{D}$ matrices and spin-weighted spherical harmonics directly in terms of quaternions, and in more standard forms

Project maintained by moble

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.

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 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.

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.