What are possible ways to construct J-matrices (higher order Pauli matrices)? I'm looking for possible ways to construct $J$-matrices. $J$-matrices are the higher-order version of Pauli matrices. Pauli matrices are suited for spin-1/2 systems, while J-matrices can be for any spin system, 1 or 3/2 or 2 or anything else.
My favorite way:
My favorite way for constructing J-matrices is using Wigner D-matrices. Starting from the eigenvalues of the $J_z$ operator (assuming $z$ is the quantization axis), which are very trivial to systematically derive. For example, for a spin-1 system:
$$J_z=\left(\begin{matrix}1 & 0 & 0\\
0 & 0 & 0\\
0\text{} & 0 & -1
\end{matrix}\right)$$
Now we rotate this operator with Wigner D-matrices to get a hypothetical vector that points along $z$ to get once along $x$ and once along $y$. This will create $J_x$ and $J_y$.
However, this assumes that I know the Wigner D-matrices, which I don't (because my problem is computational, and I don't want to calculate the Wigner D-matrices).
So my question is: Is there a simpler way to derive the $J$ matrices? The simplest way I would love is if there's a way to derive them from the Pauli-matrices.
 A: As a matter of fact, there is a much easier way to derive those matrices. What you are really after is the matrix representation of the $\mathfrak{su}(2)$ Lie algebra generators $J_x$, $J_y$ and $J_z$ in the irreducible $\mathrm{SU}(2)$ representation with spin $j$, given by
$$
J_z|j,m\rangle = m |j,m\rangle,\quad J_\pm|j,m\rangle = C_\pm(j,m) |j,m\pm 1\rangle,
$$
with
$$
J_\pm=J_x \pm i J_y
$$
and
$$
C_\pm(j,m)=\sqrt{(j\mp m)(j\pm m+1)}.
$$
The matrix representations of $J_z$ and $J_\pm$ in this basis are
$$
J_z=\begin{pmatrix}
j\\
& j-1\\
&& \ddots\\
&&& -j+1\\
&&&& -j
\end{pmatrix}
$$
$$
J_+=
\begin{pmatrix}
0 & C_+(j,j-1)\\
& 0 & C_+(j,j-2)\\
&& \ddots\\
&&& 0 & C_+(j,-j)\\
&&&& 0
\end{pmatrix}
$$
$$
J_-=
\begin{pmatrix}
0 \\
C_-(j,j) & 0\\
&& \ddots\\
&& C_-(j,-j+2) &0\\
&&&  C_-(j,-j+1) &0
\end{pmatrix},
$$
and you can recover $J_x$ and $J_y$ as
$$
J_x=\frac{J_+ + J_-}{2},\quad J_y=\frac{J_+ - J_-}{2i}.
$$
For example, with $j=1$ you would get
$$
J_z=\begin{pmatrix}1 & 0 & 0\\
0 & 0 & 0\\
0& 0 & -1
\end{pmatrix},\quad
J_x=\frac{\sqrt{2}}{2}\begin{pmatrix}0 & 1 & 0\\
1 & 0 & 1\\
0& 1 & 0
\end{pmatrix},\quad
J_y=\frac{\sqrt{2}}{2}\begin{pmatrix}0 & -i & 0\\
i & 0 & -i\\
0& i & 0
\end{pmatrix}.
$$
