Explanation: Instead of speaking of sphere parameterization, I would say that it is a variant of the Rodrigues formula for rotations, written with tensor notations.
Let us make the connection. If you apply your formula:
$$\begin{equation}
R_{ij}=\cos r\left(\delta_{ij}-\hat{\theta}_{i}\hat{\theta}_{j}\right)+\hat{\theta}_{i}\hat{\theta}_{j}+\sin r\sum\epsilon_{ijk}\hat{\theta}_{k},
\end{equation}$$
to a vector $V_i$, you get:
$$\begin{equation}
R_{ij}V_i=\cos r\left(\delta_{ij}-\hat{\theta}_{i}\hat{\theta}_{j}\right)V_i+\hat{\theta}_{i}\hat{\theta}_{j}V_i+\sin r\sum\epsilon_{ijk}\hat{\theta}_{k}V_i, \ \ \ (1)
\end{equation}$$
Have now a look at (https://en.wikipedia.org/wiki/Rodrigues'_rotation_formula). Do not stop at its main formula, but proceed a little further in the demontration where you will find:
$$V_{rot}=\cos\theta(V-V_{\parallel})+V_{\parallel}+\sin\theta \ k \times V \ \ \ (2) \ \ \text{(a version of Rodrigues' formula)}$$
(where $k$ is the unit vector giving the axis of rotation, and $V_{\parallel}=(k.V)k=k(k.V)=k(k^TV)$ is the projection of $V$ onto the axis directed by unit vector $k$).
With your notations, where $\hat{\theta}$ is also a unit vector, (2) can be written under the form:
$$V_{rot}=\cos r(IV-(\hat{\theta}\hat{\theta}^TV))+\hat{\theta}\hat{\theta}^TV+\sin r \ \hat{\theta} \times V$$
The identification between (1) and (2) is almost finished. It remains to understand the last term of (1) which is decrypted as a matrix applied to a vector, this matrix being the matrix associated with the cross product with $\hat{\theta}$:
$$\Lambda=\Lambda_{\hat{\theta}}=\begin{bmatrix}
0 & -\hat{\theta}_3&\hat{\theta}_2 \\
\hat{\theta}_3 & 0 & -\hat{\theta}_1 \\
-\hat{\theta}_2 & \hat{\theta}_1 & 0
\end{bmatrix} \ \ \text{which is such that } \ \ \Lambda V = \hat{\theta} \times V$$
One understands now that tensor $\epsilon_{ijk}$ is there for providing its $+1,-1$ or $0$ at the right places in matrix $\Lambda$ (see below the definition of this tensor).
What is asked to you is to prove that a rotation fulfills
$$RR^T=I$$
which is of course true (but also true as well for symmetries with respect to planes).
Of course, you can prove this formula using tensor notations; in my opinion, vector-matrix notations (i.e., non tensorial) are more easy to work with, at least for this issue.
Explanation about the Levi-Civita tensor: $\epsilon_{ijk}$ may take 3 values:
values $+1$ or $-1$ according to the parity of permutation $(i,j,k)$ in the case all $i,j,k$ are different,
value $0$ if two or more indices are identical.
