(* This Mathematica notebook contains the definitions of several tensor equivalent operators (Subscript[O, 2k]^q) (1≤k≤3 and -2k≤q≤2k) in momentum space, assuming J = L + S is a good quantum number. Its contents are provided freely and are intended for educational use only. The author accepts no liability for the quality of the information provided or for it being correct, complete or up to date. *)

(* The operators defined below are the easiest to program and are those with the closest relationship to the (Subscript[C, q]^((k))) crystal-field tensor operators which appear in the parent Hamiltonian. Their definition and properties can be found in D. Smith and J.H.M. Thornley, Proc. Phys. Soc. 89,779 (1966) *)

osmith[k_,q_,j_,m_,mm_]:=(((-1))^(j-m))ThreeJSymbol[{j,-m},{k,q},{j,mm}]((1/(2^k))(Sqrt[(((2j+k+1)!)/((2j-k)!))]))

(* Buckmaster operators' definition differs by that of Smith and Thornley by a numerical factor only: for details see H.A. Buckmaster, Can. J. Phys. 40, 1670 (1962) *)

conv={{((Sqrt[6])/4),((Sqrt[6])/2),(1/2),-((Sqrt[6])/2),((Sqrt[6])/4)},{((Sqrt[70])/16),((Sqrt[35])/4),((Sqrt[10])/8),((Sqrt[5])/4),(1/8),-((Sqrt[5])/4),((Sqrt[10])/8),-((Sqrt[35])/4),((Sqrt[70])/16)},{((Sqrt[231])/32),((3(Sqrt[77]))/16),((3(Sqrt[14]))/32),((Sqrt[105])/16),((Sqrt[105])/32),((Sqrt[42])/16),(1/16),-((Sqrt[42])/16),((Sqrt[105])/32),-((Sqrt[105])/16),((3(Sqrt[14]))/32),-((3(Sqrt[77]))/16),((Sqrt[231])/32)}};obuck[k_,q_,j_,m_,mm_]:=((osmith[k,q,j,m,mm])/(conv[[k/2,k+q+1]]))

(* The widely-used Stevens' operator equivalents, originally introduced in K.W.H. Stevens, Proc. Phys. Soc. 65, 209 (1952) can be expressed in a simple way in terms of the above Buckmaster operators. *)

ostevens[k_,q_,j_,m_,mm_]:=If[q<0,(1/(2I))(obuck[k,-q,j,m,mm]-obuck[k,q,j,m,mm]),(1/2)(obuck[k,q,j,m,mm]+obuck[k,-q,j,m,mm])]

(* Example: the matrix representing the fourth-order rank crystal-field potential for a system with one single f-electron (e.g. (Ce^(3+)), ground multiplet (^2)(Subscript[F, 5/2])) is calculated below. *)

J=5/2;O40=Table[ostevens[4,0,J,m,mm],{m,-J,J},{mm,-J,J}];O44=Table[ostevens[4,4,J,m,mm],{m,-J,J},{mm,-J,J}];MatrixForm[O40+5O44]