quantum mechanics - Why are rotation matrices always unitary operators?

Can someone explain why the rotation matrix is a unitary, specifically orthogonal, operator?

Answer

You can define and do the geometry several ways but I'd say the reasons are linearity, isometry and handedness (preservation of left/right handedness: this one is not needed to prove orthogonality so it's a bit more than what you asked for, but it is what sets rotations aside from other isometries). Handedness is sometimes rather loftily called chirality.

Intuitively, you need to think of a grid of $x$, $y$ and $z$ co-ordinates being ruled throughout the space taken up by the rotated object and also think of what happens to drawings and 3D sculptures in that space.

After the rotation, all the $x$, $y$ and $z$ gridlines are still orthogonal and not distorted. Distances between all mapped points are the same as what they were before the rotation, and so angles between vectors are left unchanged.

We know a rotation leaves at least one point in space fixed. So let's arbitrarily put our origin at such a point. Then the lack of global distortion in our grid shows that the transformation is linear. I say lack of "global distortion" because some nonlinear transformations (conformal ones) can also have zero local distortion - little drawings are undistorted - but beget distortion in big enough drawings and sculptures.

So, with the origin fixed, our transformation is linear and homogeneous. So our transformation $\mathcal{U}$ can be represented by a matrix $\mathbf{U}$ so that:

$$\mathcal{U}:\mathbb{R}^N\to \mathbb{R}^N:\;X\mapsto \mathbf{U}\,X$$

Now as discussed above, lengths of positions vectors stay the same, as do angles between position vectors. This means the inner product $\left = X^T\,Y = Y^T X$ between any pair of position vectors $X$ and $Y$ is unchanged. Therefore:

$$\left<\mathbf{U}\,X,\,\mathbf{U}\,Y\right> = (\mathbf{U}\,X)^T\,(\mathbf{U}\,Y) = X^T\,\mathbf{U}^T\,\mathbf{U}\,Y = \left = X^T\,Y$$

or rather:

$$X^T\,\left(\mathbf{U}^T\,\mathbf{U} - \mathbf{I}\right)\,Y = 0;\;\forall X,\,Y\in\mathbb{R}^N$$

it is now not hard to show, since we can put any pair of basis vectors $X$, $Y$ into the above equation, that we must have $\mathbf{U}^T\,\mathbf{U} = \mathbf{I}$ as an identity. Therefore the matrix must be orthogonal. Here naturally $\mathbf{I}$ is the identity matrix.

This essentially answers your question, but rotations are not the only orthogonal transformations. Reflexions are too; in $\mathbb{R}^3$ the orthogonal matrix $\operatorname{diag}(1,-1,1)$ reflects in the $x-y$ plane and it fulfills $\mathbf{U}^T\,\mathbf{U} = \mathbf{I}$. So let's go a little further. Any rotation of angle $\theta_0$ can be thought of as being joined to the identity transformation (rotation through angle of nought) through a continuous path of rotations, all about the same axis and with angles between $0$ and $\theta_0$. It belongs to the identity connected component of the group of all orthogonal transformations. Therefore:

$$\mathbf{U} = \exp(\theta\,\mathbf{H})$$

for some constant matrix $\mathbf{H}$. By imposing the orthogonality condition on the expression we get $\mathbf{U}$ orthogonal iff $\mathbf{H} = -\mathbf{H}^T$, i.e. $\mathbf{H}$ is skew-symmetric. This then is the general form of an $N$ dimensional rotation: it is a matrix of the form $\exp(\mathbf{H}_\theta)$ for some skew-symmetric $\mathbf{H}_\theta$. In three dimensions, the most general such matrix is:

$$\theta\,\mathbf{H} = \theta\,\left(\begin{array}{ccc}0& \gamma_z& -\gamma_y\\\gamma_z&0&\gamma_x\\\gamma_y&-\gamma_x&0\end{array}\right) $$

where $\gamma_x^2 + \gamma_y^2 +\gamma_z^2 = 1$, $(\gamma_x, \gamma_y,\gamma_z)$ is a unit vector defining the axis of rotation, as you can prove by finding the eigenvectors and values of $\mathbf{H}$ and showing that this vector is the eigenvector corresponding to an eigenvalue of 0 (therefore the exponential of $\exp(\theta\,\mathbf{H})$ has this vector as an eigenvector and its eigenvalue is $e^0 = 1$, i.e. it is an axis left invariant by the transformation). Also note that $\det \mathbf{U} = \exp(\operatorname{trace}(\theta\,\mathbf{H})) = 1$, as in the comments. This is the last ingredient, namely handedness I spoke of at the beginnning. A reflexion has a determinant of $-1$ and maps a right handed co-ordinate system into a left handed one and contrawise. You can find the wonted expressions for rotation operators using the Rodrigues formula grounded on the $\mathbf{H}$ matrix's characteristic equation: working through this reasoning in 3D: the three eigenvalues of $\mathbf{H}$ are $0,\, \pm i$, so by the Cayley-Hamilton theorem:

$$\mathbf{H}^3= -\mathbf{H}$$

which relationship is then used to simplify the exponential's Taylor series:

$$\exp(\theta\mathbf{H}) = \mathbf{I} + \theta \mathbf{H} + \frac{\theta^2}{2!}\mathbf{H}^2 + \cdots$$

leading to:

$$\mathbf{U} = \mathbf{I}+\sin\theta\,\mathbf{H} +(1-\cos\theta)\,\mathbf{H}^2$$

whence can be worked out the wonted formulas for a 3D rotation of angle $\theta$ about an axis defined by the unit vector $(\gamma_x, \gamma_y,\gamma_z)$.

In higher dimensions, a real valued skew-symmetric matrix $\mathbf{H}$ has the eigenvalue $0$ (possibly repeated) as well as imaginary eigenvalues in conjugate pairs $\pm i\,\theta_j$. It should be mentioned here that $\mathbf{U}$, being orthogonal is also normal (commutes with its adjoint - here equal to its transpose) and so it can always be diagonalised (has a strictly diagonal Jordan normal form) and its eigenvectors are all orthogonal. The rotation then has an invariant hyperspace given by the kernel (nullspace) of $\mathbf{H}$ - this is the generalization of the rotation axis in 3D and then one or more linearly independent 2D hyperplanes (indeed orthogonal hyperplanes, given normalness of $\mathbf{U}$), each spanned by the pair of eigenvectors corresponding to the eigenvalues $\pm i\,\theta_j$. So the idea of an "axis" is no longer really useful: in 3D it is useful because the nullspace of $\mathbf{H}$ must be precisely one dimensional. Sometimes authors require "rotations" to be transformations that leave all of $\mathbb{R}^N$ invariant aside from precisely one 2D hyperplane, but I don't think this is particularly useful because the composition of two such transformation is not then a rotation (unless the hyperplane is the same one for the two composed "rotations"): there is no group of rotations defined in this way. It's easier and more useful simply to talk of simply orthogonal transformations with unit determinant, i.e. members of the group $SO(N)$.

If, as Michael Brown's comment suggests, you are thinking of representations of rotations, then further discussions of the lifting of $SO(3)$ to its universal (in this case double) cover $SU(2)$ can be found in the second section "What the Lie Bracket does not "remember" about the group: Global Topology and the Fundamental Group" in my answer here and especially in the Stillwell references my answer gives.