Spin matrix and (Subscript[S, z]) matrix relative to a specific normal vector n.  This is based on PauliMatrixSpinOperators.nb from ps8 p4, hyperfine interaction, but switched up for the Fermion/Boson interaction final exam question p6.

ClearAll[n, hbar, s, sz,display,nplus, nminus, absket, theta, phi, braket, szSpinBasis, matrixRepresentation, normalizedEigenVectors, szKets, szPlus, szMinus, u, uInverse, spinMatrixInSzBasis, nPlusSzBasis, nMinusSzBasis,h, h0, v, a, b, fullfieldEnergies, z2]
On[Assert] ;

(*hack: had to omit theta = 0, π points to get the S eigenkets to simplify*)
$Assumptions=(And@@(#>0&)/@{hbar, a, b})&& theta < Pi && phi ≤ (2Pi) && theta >0 && phi ≥ 0;

n = {Sin[theta]Cos[phi],Sin[theta]Sin[phi],Cos[theta]} ;

display := {hbar -> ℏ, theta -> θ, phi -> ϕ, a -> "J", b-> "B", E -> "e", Complex[0,-1] -> -i, I -> i} ;
s = (hbar/2)n[[#]] PauliMatrix[#] & /@Range[3] // Total ;
sz = (hbar/2) PauliMatrix[3] ;

(* grab the two eigenvectors ordered low to high *)
normalizedEigenVectors[m_] := Module[{energies,eigenvectors, norms},
{energies,eigenvectors} = (m /. hbar -> 1)// Eigensystem ; 
Assert[(energies // First) < (energies // Last)] ;
norms = Norm[eigenvectors[[#]]] & /@ Range[2] ;

(eigenvectors[[#]]/norms[[#]] // FullSimplify)&/@ Range[2] 
]  ;

{nminus, nplus} = normalizedEigenVectors[s] ;

braket[f_, g_] := Conjugate[f].g ;
matrixRepresentation[m_, one_, two_] := {{braket[one, m.one],braket[one, m.two]},
{braket[two, m.one],braket[two, m.two]}} ;

(((hbar/2)PauliMatrix[3] == matrixRepresentation[s, nplus, nminus])// FullSimplify) // Assert
szSpinBasis = matrixRepresentation[sz, nplus, nminus] // FullSimplify ;

{szMinus,szPlus} = normalizedEigenVectors[szSpinBasis] ;
(((hbar/2)PauliMatrix[3] == matrixRepresentation[szSpinBasis, szPlus, szMinus])// FullSimplify) // Assert

u = {szPlus, szMinus}  // Transpose ;
uInverse = u // Inverse ;

((u . (PauliMatrix[3].uInverse) // FullSimplify) == (szSpinBasis/(hbar/2))) // Assert
((uInverse. (szSpinBasis/(hbar/2)).  u // FullSimplify) == PauliMatrix[3]) // Assert

spinMatrixInSzBasis =(uInverse. s.  u // FullSimplify) ;
nPlusSzBasis = (u// Inverse) . nplus ;
nMinusSzBasis = (u// Inverse) . nminus ;

(((hbar/2)PauliMatrix[3] == matrixRepresentation[spinMatrixInSzBasis, nPlusSzBasis, nMinusSzBasis])// FullSimplify) // Assert

z2 = DiagonalMatrix[{0,0}] ;
h0 = (a hbar^2/2) ArrayFlatten[{{PauliMatrix[3],z2, z2},{z2,z2,z2},{z2,z2,-PauliMatrix[3]}}];

v = b ArrayFlatten[{{szSpinBasis,z2, z2},{z2, szSpinBasis, z2},{z2, z2,szSpinBasis}}];
h = h0 + v ;

fullfieldEnergies =Eigenvalues[ h ] ;

Style[ {#1," = ",(#2 // FullSimplify)/. display  // TraditionalForm} &
@@@{
{"(OverHat[n])", n // MatrixForm},
(*{"(OverVector[S]) · (OverHat[n])", Row[{hbar/2,(s 2/hbar)// MatrixForm }]},*)
{"|(OverVector[S]) · (OverHat[n]);+>",nplus// MatrixForm},
{"|(OverVector[S]) · (OverHat[n]);->",nminus// MatrixForm},
(*{"(OverVector[S]) · (OverHat[z])", Row[{hbar/2,(sz 2/hbar)// MatrixForm }]},*)
{"(S basis) : (OverVector[S]) · (OverHat[z])", Row[{hbar/2,(szSpinBasis 2/hbar)// MatrixForm }]},
{"(S basis) : |(OverVector[S]) · (OverHat[z]);+>", szPlus // MatrixForm},
{"(S basis) : |(OverVector[S]) · (OverHat[z]);->", szMinus // MatrixForm},
(*{"U = [ |(OverVector[S]) · (OverHat[z]);+> |(OverVector[S]) · (OverHat[z]);-> ]", u // MatrixForm},*)
{"(Subscript[H, 0])", Row[{a hbar^2/4,(h0 /(a hbar^2/4))// MatrixForm }]},
{"V", Row[{b hbar/2,(v /(b hbar/2))// MatrixForm }]},
{"(full field) E", Column[fullfieldEnergies] }
(*,
{"((Subscript[S, z]) basis) : S' = (SuperDagger[U]) S U", Row[{hbar/2,(spinMatrixInSzBasis 2/hbar)// MatrixForm }]},
{"((Subscript[S, z]) basis) : |n;+>",nPlusSzBasis// MatrixForm},
{"((Subscript[S, z]) basis) : |n;->",nMinusSzBasis// MatrixForm}*)

} // Grid,{ 14, "Text"}] 

(*Column[fullfieldEnergies] /. display // TraditionalForm*)