In mathematical jargon, atoms within 3D crystals exist in $\mathbb{R}^{3}$, or 3-space. In $\mathbb{R}^{3}$, the following $3 \times 3$ matrices $\textbf{R}(\phi)$ codify rotations by an angle $\phi$ about the $x$-, $y$-, and $z$-axes, respectively:
\begin{equation}
\textbf{R}_{x}(\phi)=\left[\begin{array}{ccc}
1 & 0 & 0 \\
0 & \cos(\phi) & -\sin(\phi) \\
0 & \sin(\phi) & \cos(\phi)
\end{array}\right] \\
\end{equation}
\begin{equation}
\textbf{R}_{y}(\phi)=\left[\begin{array}{ccc}
\cos(\phi) & 0 & \sin(\phi) \\
0 & 1 & 0 \\
-\sin(\phi) & 0 & \cos(\phi)
\end{array}\right] \\
\end{equation}
\begin{equation}
\textbf{R}_{z}(\phi)=\left[\begin{array}{ccc}
\cos(\phi) & -\sin(\phi) & 0 \\
\sin(\phi) & \cos(\phi) & 0 \\
0 & 0 & 1
\end{array}\right]
\end{equation}
In each matrix, the axis of rotation is specified by whichever row and column have been extracted from the identity matrix $\textbf{I}$. If we post-multiply any of these matrices with a column vector $\textbf{v}$ representing 3D coordinates in real space, we obtain a modified set of coordinates $\textbf{u}$ corresponding to the rotated point. Let's select $\textbf{R}_{z}(\phi)$ as an arbitrary example:
\begin{equation}
\begin{bmatrix}
\cos(\phi) & -\sin(\phi) & 0 \\
\sin(\phi) & \cos(\phi) & 0 \\
0 & 0 & 1
\end{bmatrix}
\underbrace{
\begin{bmatrix}
x \\
y \\
z
\end{bmatrix}}_{\textbf{v}}
=
\underbrace{
\begin{bmatrix}
x \cos(\phi) - y \sin(\phi) \\
x \sin(\phi) + y \cos(\phi) \\
z
\end{bmatrix}}_{\textbf{u}}.
\end{equation}
An underappreciated wealth of information is buried in the eigenvalues of these matrices, which we'll delve into here. To solve for $\lambda$, we need to find the roots of the characteristic polynomial $\operatorname{det}(\textbf{R}_{z}-\lambda \textbf{I})=0$. Since we're dealing with $3 \times 3$ matrices, a convenient strategy to use is Laplace expansion:\begin{equation}
\operatorname{det}(\textbf{R}_{z}-\lambda \textbf{I})=\left|\begin{array}{ccc}
\cos (\phi)-\lambda & -\sin (\phi) & 0 \\
\sin (\phi) & \cos (\phi)-\lambda & 0 \\
0 & 0 & 1-\lambda
\end{array}\right|=0,
\end{equation}
\begin{equation}
\operatorname{det}(\textbf{R}_{z}-\lambda \textbf{I})=(1-\lambda)\left|\begin{array}{cc}
\cos (\phi)-\lambda & -\sin (\phi) \\
\sin (\phi) & \cos (\phi)-\lambda
\end{array}\right|=0.
\end{equation}
At this stage, it's already clear that one of the eigenvalues is $\lambda = 1$. This signifies that whatever our rotation matrix is doing, it's leaving one particular dimension invariant. In this case, since we arbitrarily picked $\textbf{R}_{z}$, that corresponds to the $z$-axis. Let's refer to this specific eigenvalue as $\lambda_{3}$. To find the remaining two eigenvalues, let's evaluate the $2 \times 2$ determinant left over from the third-row Laplace expansion:
\begin{equation}
\operatorname{det}(\textbf{R}_{z}-\lambda \textbf{I})=(1-\lambda)\left[\cos ^{2} (\phi)-2 \lambda \cos (\phi)+\lambda^{2}+\sin ^{2} (\phi)\right]=0.
\end{equation}
Invoking the trigonometric identity $\cos ^{2} (\phi) + \sin ^{2} (\phi) = 1$, we find
\begin{equation}
\operatorname{det}(\textbf{R}_{z}-\lambda \textbf{I})=(1-\lambda)\left[\lambda^{2}-2 \lambda \cos (\phi) +1\right]=0.
\end{equation}
Applying the quadratic formula to the bracketed polynomial, we arrive at a weirdly profound result:
\begin{equation}
\lambda =\cos (\phi) \pm \sqrt{\cos^{2} (\phi)-1} =\cos (\phi) \pm \sqrt{-\sin ^{2} (\phi)} =\cos (\phi) \pm i \sin (\phi)= \mathrm{exp}(\pm i \phi).
\end{equation}
That's right, Euler's formula is lurking here too. In other words, $\textbf{R}_{z}$ has three eigenvalues: $\lambda_{1}=\mathrm{exp}(i\phi)$, $\lambda_{2}=\mathrm{exp}(-i\phi)$, and $\lambda_{3}=1$. Critically, notice how $\lambda_{1}$ and $\lambda_{2}$ represent complex conjugates:
\begin{equation}
\lambda_{1,2} = \mathrm{exp}(\pm i \phi).
\end{equation}
A useful fact about complex conjugates is that their sum is always real. We can demonstrate this by invoking Euler's formula:
\begin{equation}
\lambda_{1} + \lambda_{2} = \mathrm{exp}(i \phi)+\mathrm{exp}(-i \phi)=\mathrm{cos}(\phi) + i\,\mathrm{sin}(\phi) + \mathrm{cos}(-\phi) + i\,\mathrm{sin}(-\phi).
\end{equation}
Recall that cosine is an even function (i.e., in crystallographic parlance, centrosymmetric), so $\cos(\phi) = \cos(-\phi)$. On the other hand, sine is an odd function (i.e., antisymmetric), which means $-\sin(\phi) = \sin(-\phi)$. Substituting these identities, the imaginary sine components cancel, yielding
\begin{equation}
\lambda_{1} + \lambda_{2} = \mathrm{exp}(i \phi)+\mathrm{exp}(-i \phi)=2\,\mathrm{cos}(\phi).
\end{equation}
Piecing all these fragments together, we can derive an immensely powerful expression:
\begin{equation}
\sum_{j=1}^{3}\lambda_{j} = 1 + 2\,\mathrm{cos}(\phi).
\end{equation}
In other words, the sum of a rotation matrix's eigenvalues also encodes the angle of rotation $\phi$. Rearranging the previous equation to solve for $\phi$, we obtain
\begin{equation}
|\phi| = \cos^{-1}{\left(\frac{-1+\sum_{j=1}^{3}\lambda_{j}}{2}\right)}.
\end{equation}
Armed with this formula, we can easily deduce $\phi$ for any arbitrary 3D1 rotation matrix. Incidentally, you'll occasionally encounter this rule written as $\mathrm{tr}[\textbf{R}(\phi)]=1+2\,\mathrm{cos}(\phi)$. This holds true because of what we've just derived, since the trace of a matrix is equivalent to the sum of its eigenvalues (including degeneracies)!
Now let's plug in a few arguments for $\phi$ and analyze the resultant eigenvalues. In the trivial case where $\phi=0$, $\textbf{R}_{z}$ reduces to the identity matrix, and the eigenvalue $\lambda=1$ becomes triply degenerate:
\begin{equation}
\textbf{R}_{z}(0)=\left[\begin{array}{ccc}
\cos(0) & -\sin(0) & 0 \\
\sin(0) & \cos(0) & 0 \\
0 & 0 & 1
\end{array}\right]
=\left[\begin{array}{ccc}
1 & 0 & 0 \\
0 & 1 & 0 \\
0 & 0 & 1
\end{array}\right]=\textbf{I},
\end{equation}
\begin{equation}
\lambda_{3}=1, \lambda_{1}=\lambda_{2}=\mathrm{exp}(0)=1.
\end{equation}
These eigenvalues tell us that all three dimensions emerge unscathed from this transformation (obviously, since it's a rotation by 0 radians). Another special case arises when $\phi=\pi$. We obtain
\begin{equation}
\textbf{R}_{z}(\pi)=\left[\begin{array}{ccc}
\cos(\pi) & -\sin(\pi) & 0 \\
\sin(\pi) & \cos(\pi) & 0 \\
0 & 0 & 1
\end{array}\right]
=\left[\begin{array}{ccc}
-1 & 0 & 0 \\
0 & -1 & 0 \\
0 & 0 & 1
\end{array}\right],
\end{equation}
\begin{equation}
\lambda_{3}=1, \lambda_{1}=\mathrm{exp}(i\pi)=-1, \lambda_{2}=\mathrm{exp}(-i\pi)=-1.
\end{equation}
Here, the doubly degenerate eigenvalues $\lambda_{1} = \lambda_{2} = -1$ convey that the $x$ and $y$ coordinates become inverted after a rotation by $\pi$ radians about the $z$-axis. In both trivial cases, notice how $\lambda_{1}$ is equivalent to its complex conjugate, $\lambda_{2}$—a condition which is possible if and only if $\lambda_{1}$ and $\lambda_{2}$ are real. Now here's the really thought-provoking part: for any angle $\phi$ which isn't 0 or $\pi$, $\lambda_{1}$ and $\lambda_{2}$ become complex. As an illustrative example, let's consider $\phi = \tfrac{\pi}{2}$:
\begin{equation}
\textbf{R}_{z}(\tfrac{\pi}{2})=\left[\begin{array}{ccc}
\cos(\frac{\pi}{2}) & -\sin(\frac{\pi}{2}) & 0 \\
\sin(\frac{\pi}{2}) & \cos(\frac{\pi}{2}) & 0 \\
0 & 0 & 1
\end{array}\right]
=\left[\begin{array}{ccc}
0 & -1 & 0 \\
1 & 0 & 0 \\
0 & 0 & 1
\end{array}\right],
\end{equation}
\begin{equation}
\lambda_{3}=1, \lambda_{1}=\mathrm{exp}\left(\frac{\pi}{2}i\right)=i, \lambda_{2}=\mathrm{exp}\left(-\frac{\pi}{2}i\right)=-i.
\end{equation}
At first glance, these imaginary eigenvalues may seem forbiddingly obscure: not quite as amenable to a straightforward physical interpretation. However, a pinch of geometric intuition makes their meaning startlingly clear. Let's say we have some vector $\textbf{v}$ hanging out in the complex plane. What happens when we multiply $\textbf{v}$ by $i$? Algebraically, we just get the trivial result $i\textbf{v}$. Geometrically, however, multiplication of $\textbf{v}$ by $i$ is equivalent to rotating $\textbf{v}$ by $\tfrac{\pi}{2}$ radians in the complex plane (analogously, multiplication by $-i$ denotes a rotation of $-\tfrac{\pi}{2}$ radians). This is why the imaginary eigenvalues of $\textbf{R}_{z}(\tfrac{\pi}{2})$—i.e., the transformation matrix which encodes a rotation by $\tfrac{\pi}{2}$ radians—work out to $i$ and $-i$! More abstractly, for any nontrivial 3D rotation matrix, the two non-unity eigenvalues always divulge how that particular symmetry operation would manifest itself in the complex plane. An interesting ambiguity, though, is that the directionality of the rotation is irretrievably lost.
A scrollable table containing matrices and eigenvalues for several crystallographic symmetry operators is given below.
1. Notice how $\sum_{j=1}^{3} \lambda_{j}$ is explicitly specified here, because this equation isn't generalizable to rotations in either lower- or higher-dimensional spaces. For instance, in $\mathbb{R}^{2}$, we only need one $2 \times 2$ matrix to describe all possible rotations in 2-space, and that matrix has two eigenvalues: $\exp(\pm i \phi)$. Collapsing the third dimension gets rid of the dangling unity eigenvalue, so the trace formula in 2D reduces to $\sum_{j=1}^{2} \lambda_{j} = 2 \cos(\phi)$. If we expand to $\mathbb{R}^{4}$, we need six $4 \times 4$ matrices to describe all possible rotations in 4D. Complexity starts to ramp up pretty rapidly in 4-space, since now we also need two discrete angles, $\phi$ and $\psi$, to fully describe a proper 4D rotation. These $4 \times 4$ matrices have four eigenvalues: $\exp(\pm i \phi)$ and $\exp(\pm i \psi)$. If $\phi = \psi$, that's referred to as an isoclinic rotation, and in that particular case $\sum_{j=1}^{4} \lambda_{j} = 4 \cos(\phi)$. If $\phi \neq \psi$, $\sum_{j=1}^{4} \lambda_{j} = 2[\cos(\phi) + \cos(\psi)]$. Thankfully, we don't have to deal with 4D crystals... ↩