Real representations of the Lie algebra $\mathfrak{so}(3)$

Recently I wondered about the real (irreducible) representations of the Lie group SO(3) and its Lie algebra $\mathfrak{so}(3)$. It is well known, that the representations of SO(3) are odd-dimensional, there are no (irreducible) representations of SO(3) in even dimensions.

Usually, in textbooks the authors quickly move from SO(3) to the covering group SU(2) and derive the (complex) representations of the corresponding Lie algebra $\mathfrak{su}(2)$ using quantum mechanical angular momentum algebra and ladder operators.

There appears to be surprisingly little in the published literature on real-valued representations of $\mathfrak{so}(3)$. How do the real-valued representations look like? A google search lead me to

V. M. Gordienko
Matrix entries of real representations of the group O(3) and SO(3)
Siberian Mathematical Journal (2002), 43(1):36-46
DOI: 10.1023/A:1013816403253

who describes their construction; but see also

G. Itzkowitz, S. Rothman and H. Strassberg
A note on the real representation of SU(2,C)
Journal of Pure and Applied Algebra (1990), 69:285-294
DOI: 10.1016/0022-4049(91)90023-U

If I understand correctly, there is no ladder operator approach for real representations. Rather, Gordienko constructs real representations using homogeneous polynomials $P^{(N)}(x,y)$ in $x$ and $y$ with $x,y \in \mathbb{C}$. I.e.
$$\begin{align} P^{(1)}(x,y) =& a_1 x^2 + a_2 x y + a_3 y^2 \\ P^{(2)}(x,y) =& a_1 x^4 + a_2 x^3 y + a_3 x^2 y^2 + a_4 x y^3 + a_5 y^4 \notag\\ & \ldots \notag\\ P^{(N)}(x,y) =& a_1 x^{2N} + a_2 x^{2N-1} y + a_3 x^{2N-2} y^2 \notag\\ &+ \ldots + a_{N+1} x^{N} y^{N} + \ldots + a_{2N} x y^{2N-1} + a_{2N+1} y^{2N} \notag\\ \end{align}$$ Gordienko chooses the following basis in the vector space of polynomials $P^{(N)}(x,y)$
$$\begin{align} e^{(N)}_n(x,y) =& i^{N+1}\sqrt{\frac{(2N+1)!}{2(N-n)!(N+n)!}} \\ &\times(x^{N-n}y^{N+n} - (-1)^{-n}x^{N+n}y^{N-n}) \quad \mathrm{if}\; n \le -1 \notag\\ e^{(N)}_{n=0}(x,y) =&  i^{N} \frac{\sqrt{(2N+1)!}}{N!} x^N y^N \quad \mathrm{if}\; n = 0 \notag\\ e^{(N)}_n(x,y) =& -i^{N}\sqrt{\frac{(2N+1)!}{2(N+n)!(N-n)!}} \notag\\ & \times(x^{N+n}y^{N-n} + (-1)^{n}x^{N-n}y^{N+n}) \quad \mathrm{if}\; n \ge 1 \notag \end{align}$$ with $n = -N,\ldots,N$. Now let's take a generic rotation matrix from the group SU(2), parametrized by the three real numbers $\alpha$, $\beta$ and $\gamma$ and we write with $\delta\equiv\sqrt{1-\alpha^2-\beta^2-\gamma^2}$
$$M =  \begin{pmatrix} \delta - i\gamma & \beta  -i\alpha\\ -\beta - i\alpha &\delta + i\gamma \end{pmatrix} \equiv \begin{pmatrix} A^* & B^*\\ -B & A \end{pmatrix}$$ and apply it on the 2-dimensional complex vector (x,y), i.e.
$$(x',y') \equiv (x,y)  \begin{pmatrix} A^* & B^*\\ -B & A \end{pmatrix}=(A^*x-By, B^*x+Ay)$$ For the "rotated" basis $e^{(N)}_n(x',y')$ we obtain
$$\begin{align} e^{(N)}_n(x',y') =& e^{(N)}_n(A^*x-By, B^*x+Ay) \\ =& \sum_{m=-N}^{+N} T^{(N)}_{n,m}(A,B) \; e^{(N)}_m(x,y)\notag \end{align}$$ The $(2N+1)\times(2N+1)$ matrix $T^{(N)}_{m,n}(A,B)$ is the $(2N+1)$-dimensional representation matrix of SO(3) we're looking for.

Why is $T^{(N)}_{m,n}(A,B)$ real-valued? Turns out, the basis vectors $e^{(N)}_n(x,y)$ are constructed in such a way that
$$e^{(N)}_m(x,y) = (e^{(N)}_m(-y,x))^*$$ holds, where $x^*$ denotes the complex conjugate of $x$. Gordienko shows that for basis vectors with this property the entries of the representation matrices $T^{(N)}_{n,m}(A,B)$ are indeed real. Specifically,
$$\begin{align} e^{(N)}_m(x',y') =& e^{(N)}_m(A^*x-By, B^*x+Ay) \\ =&  (e^{(N)}_m(-B^*x-Ay,A^*x-By))^* \notag\\ =&  (e^{(N)}_m(B^*(-x)+A(-y),A^*x-By))^* \notag\\ =&  (e^{(N)}_m(-y',x'))^* \notag\\ =&  (\sum_{n=-N}^N T^{(N)}_{m,n}(A,B)\;e^{(N)}_n(-y,x))^* \notag\\ =&  \sum_{n=-N}^N (T^{(N)}_{m,n}(A,B))^*\;(e^{(N)}_n(-y,x))^* \notag\\ =&  \sum_{n=-N}^N (T^{(N)}_{m,n}(A,B))^*\;e^{(N)}_n(x,y) \notag \end{align}$$ Thus, $T^{(N)}_{m,n}(A,B) = (T^{(N)}_{m,n}(A,B))^*$ since the $e^{(N)}_n(x,y)$ are linear independent; in other words, the matrix $T^{(N)}_{m,n}(A,B)$ is real.

representationofso3.m a is a rough-and-ready MATLAB programme (included in LieTools), which outputs numerical approximations of the n-dimensional representations of $\mathfrak{so}(3)$ (for n below about 90) using the method described above. For dimensions 3, 5 and 7 the output of representationofso3.m yields the following matrices. The 3-dimensional representation is familiar
$$\begin{equation} t_1 = \begin{pmatrix} 0 & 1 & 0 \\ -1 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}\\ t_2 = \begin{pmatrix} 0 & 0 & -1 \\ 0 & 0 & 0 \\ 1 & 0 & 0 \end{pmatrix}\\ t_3 = \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & 1 \\ 0 & -1 & 0 \end{pmatrix} \end{equation}$$ the 5-dimensional representation is $$\begin{equation} t_1 = \begin{pmatrix} 0 & 0 & 0 & -1 & 0 \\ 0 & 0 & \sqrt{3} & 0 & 1 \\ 0 & -\sqrt{3} & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 \\ 0 & -1 & 0 & 0 & 0 \end{pmatrix}\\ t_2 = \begin{pmatrix} 0 & 0 & 0 & 0 & -2 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & -1 & 0 & 0 & 0 \\ 2 & 0 & 0 & 0 & 0 \end{pmatrix}\\ t_3 = \begin{pmatrix} 0 & -1 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & -\sqrt{3} & 0 \\ 0 & 0 & \sqrt{3} & 0 & -1 \\ 0 & 0 & 0 & 1 & 0 \end{pmatrix} \end{equation}$$ and the 7-dimensional representation turns out to be
$$\begin{equation} t_1 = \begin{pmatrix} 0 & 0 & 0 & 0 & 0 & -\sqrt{3/2} & 0 \\ 0 & 0 & 0 & 0 & \sqrt{5/2} & 0 & \sqrt{3/2} \\ 0 & 0 & 0 & \sqrt{6} & 0 & -\sqrt{5/2} & 0 \\ 0 & 0 & -\sqrt{6} & 0 & 0 & 0 & 0 \\ 0 & -\sqrt{5/2} & 0 & 0 & 0 & 0 & 0 \\ \sqrt{3/2} & 0 & \sqrt{5/2} & 0 & 0 & 0 & 0 \\ 0 & -\sqrt{3/2} & 0 & 0 & 0 & 0 & 0 \end{pmatrix}\\ t_2 = \begin{pmatrix} 0 & 0 & 0 & 0 & 0 & 0 & -3 \\ 0 & 0 & 0 & 0 & 0 & 2 & 0 \\ 0 & 0 & 0 & 0 & -1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & -2 & 0 & 0 & 0 & 0 & 0 \\ 3 & 0 & 0 & 0 & 0 & 0 & 0 \end{pmatrix}\\ t_3 = \begin{pmatrix} 0 & -\sqrt{3/2} & 0 & 0 & 0 & 0 & 0 \\ \sqrt{3/2} & 0 & -\sqrt{5/2} & 0 & 0 & 0 & 0 \\ 0 & \sqrt{5/2} & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & \sqrt{6} & 0 & 0 \\ 0 & 0 & 0 & -\sqrt{6} & 0 & -\sqrt{5/2} & 0 \\ 0 & 0 & 0 & 0 & \sqrt{5/2} & 0 & -\sqrt{3/2} \\ 0 & 0 & 0 & 0 & 0 & \sqrt{3/2} & 0 \end{pmatrix} \end{equation}$$ The zero entries in the off-diagonals look somewhat strange, but the matrices are definitely real, skew-symmetric and obey the $\mathfrak{so}(3)$ commutation relations
$$[t_1,t_2] = t_3\\ [t_2,t_3] = t_1\\ [t_3,t_1] = t_2$$