Compute and display (SubPlus[S]), (SubMinus[S]),  (Subscript[S, x]), (Subscript[S, y]), (Subscript[S, z]), (S^2) for a given spin (with hbar = 1).  Use this to display the matrix representations of these operators for each of: spin 1/2, spin 1, spin 3/2, spin 2, spin 5/2.  Use this to solve the (spin one) eigensystem in Sakurai pr. 4.12.

ClearAll[spinMatrices, showSpin, sopts] ;

(*s = 3/2 ;*)
(* compute the spin-s analogue of the Pauli matrices, and the corresponding raising and lowering operators *)
spinMatrices[s_] :=Module[{mValues, m, vec, raiseOrLower, raise, lower, splus, sminus, sz, szOp, sSquared},

dimension = 2 s + 1 ;

mValues := -# + s + 1 & /@ Range[dimension] ;

(*Map the values m, in|m> to the standard basis for (R^dimension)*)
(*For[i=1,i≤ dimension,i++, vec[mValues[[i]]] = UnitVector[dimension, i]
] ;*)
(vec[mValues[[#]]] = UnitVector[dimension, #]) &/@ Range[dimension] ;

(* Any invalid value is mapped to a zero vector (i.e. (SubPlus[S])|3/2>, or (SubMinus[S])|-3/2>*)
vec[n_] := UnitVector[dimension, 1] - UnitVector[dimension, 1] ;

(* For a basis vector |m> (not a linear combination) return the numeric value m *)
m[v_] := v.DiagonalMatrix[mValues] // Tr ;

(*s = 1 for raise*)
raiseOrLower[v_, sign_Integer] := Module[{mv},
mv = m[v] ;
Sqrt[(s -sign mv)(s + sign mv + 1)] vec[mv+ sign ]
] ;

raise[v_] := raiseOrLower[v,1] ;
lower[v_] := raiseOrLower[v,-1] ;
szOp[v_] :=  m[v] v ;

(* Form the matrix for an operator that acts on vec[m] values. *)
operatorToMatrix[o_] := (o[vec[#]] &/@mValues)// Transpose ;

(* Matrixes for (SubPlus[S]), (SubMinus[S]), and (Subscript[S, z]), (S^2) *)
splus = operatorToMatrix[raise] ;
sminus = operatorToMatrix[lower] ;
sz := operatorToMatrix[szOp] ;

sSquared := s(s+1)  IdentityMatrix[dimension] ;

(*return: (SubPlus[S]), (SubMinus[S]),  (Subscript[S, x]), (Subscript[S, y]), (Subscript[S, z]), (S^2)*)
{splus, sminus, (splus +  sminus)/(2),(splus -  sminus)/(2I),sz, sSquared}
]

sopts =  {14, "Text"}  ;
showSpin[s_] := Module[{sp},
sp =  spinMatrices[s] ;
Style[{#1," = ", #2  // MatrixForm  } &@@@( {{ "(SubPlus[σ])","(SubMinus[σ])","(Subscript[σ, x])","(Subscript[σ, y])","(Subscript[σ, z])","(σ^2)"  }, sp} // Transpose)
// Grid
, sopts
] 
]

showSpin[1/2]
showSpin[1]
showSpin[3/2]
showSpin[2]
showSpin[5/2]

(* Sakurai problem 4.12 *)
ClearAll[a,b, hbar, display, h, splus, sminus, sx, sy, sz,sq, ei, kets]
$Assumptions=(And@@(#>0&)/@{hbar, a,b});


{splus, sminus, sx, sy, sz,sq} =  spinMatrices[1] ;
h = (a sz.sz + b(sx.sx - sy.sz)) hbar^2 ;
display = {hbar -> ℏ, a -> "A", b->"B"} ;
h (*/hbar^2*)/. display // MatrixForm
{ei, kets} =Eigensystem[h] ;
ei /hbar^2/. display
((#(*/Norm[#]*) /. display // MatrixForm)&/@ kets) // Column