Matrices of finite order are dense in SO(n)

Idea. An orthogonal matrix is a direct sum of plane rotations and identity blocks in a suitable orthonormal basis. The angles are free continuous parameters, so we may replace them by rational multiples of $2\pi$.

Step 1: the normal form. A real orthogonal matrix $U$ is normal ($UU^T=U^TU=I$), hence unitarily diagonalisable over $\mathbb{C}$ with eigenvalues of modulus $1$. The real spectral theorem for normal matrices (quoted) gives an orthogonal $Q$ with
$$U=Q\,\mathrm{diag}\bigl(R(\theta_1),\dots,R(\theta_k),\,I_s,\,-I_t\bigr)\,Q^T,\qquad 2k+s+t=n,$$
where the blocks $R(\theta_j)$ correspond to the pairs of non-real eigenvalues $e^{\pm i\theta_j}$ and $I_s$, $-I_t$ to the eigenvalues $1$ and $-1$.

Since $\det U=1$, the number $t$ of eigenvalues $-1$ is even. Group them in pairs: a $2\times2$ block $-I_2$ equals $R(\pi)$. So, after absorbing them into the rotation blocks, we may assume $t=0$ (this changes $k$ but not the form), and
$$U=Q\bigl(R(\theta_1)\oplus\dots\oplus R(\theta_k)\oplus I_s\bigr)Q^T.$$
Note that conjugating by any $Q\in O(n)$, even with $\det Q=-1$, keeps us inside $SO(n)$, because $\det(QXQ^T)=\det X$.

Step 2: approximating. Choose rationals with $\theta_j^{(m)}\in2\pi\mathbb{Q}$ and $\theta_j^{(m)}\to\theta_j$ as $m\to\infty$, and put
$$U_m=Q\bigl(R(\theta_1^{(m)})\oplus\dots\oplus R(\theta_k^{(m)})\oplus I_s\bigr)Q^T\in SO(n).$$
The eigenvalues of $U_m$ are $e^{\pm i\theta_j^{(m)}}$ and $1$, all roots of unity. Writing $\theta_j^{(m)}=2\pi p_j/q_j$, we get $U_m^N=I$ with $N=\mathrm{lcm}(q_1,\dots,q_k)$, so $U_m$ has finite order.

Step 3: convergence. The operator norm is invariant under orthogonal conjugation, and the norm of a block diagonal matrix is the maximum of the norms of the blocks. Hence
$$\|U_m-U\|=\max_j\|R(\theta_j^{(m)})-R(\theta_j)\|=\max_j2\Bigl|\sin\frac{\theta_j^{(m)}-\theta_j}2\Bigr|\longrightarrow0.$$

Check (numerical). For a random $U\in SO(5)$ computed with the real Schur form, rounding the two angles to multiples of $2\pi/N$ gave $\|U_N-U\|\approx0.28,\ 0.027,\ 0.002$ for $N=10,100,1000$ and $\|U_N^N-I\|\approx10^{-12}$.

Remark. The same argument works in $O(n)\setminus SO(n)$: there the normal form has an odd number of eigenvalues $-1$, which is kept fixed while the rotation angles are approximated. So the elements of finite order are dense in all of $O(n)$.

Lesson. A normal form with free continuous parameters turns a density question about a group into a density question about $\mathbb{R}$ (or the torus).