Magnetic Ion in Crystalline
Electric Field

(a)

1. Define (Subscript[Y, lm]) and \(TraditionalForm\`\(Y\_lm\%*\)\)≡Yc

Y[l_,m_,u_,v_]:=SphericalHarmonicY[l,m,u,v];
Yc[l_,m_,u_,v_]:=(-1)^m*Y[l,-m,u,v];

2. The position vectors Ri (i=1,...,6)  in spherical coordinates are given by Ri=(d,θi,ϕi). The angles can be given using a table structure as follows

ang={{Pi/2,0},{Pi/2,Pi/2},{Pi/2,Pi},{Pi/2,3*Pi/2},
     {0,0},{Pi,0}}

ang (from `angle') is a 6×2 matrix; θi=ang[[i,1]] and  ϕi=ang[[i,2]].

3. You will have to compute the sum Σi\(TraditionalForm\`\(Y\_lm\%*\)\)(θi,ϕi) , from i=1 to 6. You can do this by making use of the Sum `Mathematica' command. You can learn about this command by typing ?Sum. The function s[l,m] defined below returns the value of the above sum for a specified l and m.

s[l_,m_]:=Sum[Yc[l,m,ang[[i,1]],ang[[i,2]]],{i,1,6}]

Note: You may discard the constant term (which contains  (Subscript[Y, 00])) in V since this term causes only an overall energy shift. In fact, by doing so is equivalent to redefine the origin of the energy scale and, what is the most important, which makes the forthcoming computations easier.

(b)

1. In the spherical tensor operators from the expansion of V you will have to replace x with Lx  etc. by following the symmetrizing procedure indicated in the text of the problem. For this you can use the transformation rules feature of `Mathematica'. Suppose that F[x,y,z] contains terms like x*z2 , x*y*z , (x+i*y)4 ,etc. Then you will use the replacement rules

F[lx_,ly_,lz_]=Expand[F[x,y,z]]/.{x^2*z^2->(lx.lx.lz.lz+lx.lz.lx.lz+
               lx.lz.lz.lx+lz.lz.lx.lx+lz.lx.lz.lx+
               lz.lx.lx.lz)/6,
               x*y*z->(lx.ly.lz+lx.lz.ly+ly.lx.lz+
               ly.lz.lx+lz.lx.ly+lz.ly.lx)/6,
               (x+I y)->lp,...}

where lp=L+:=Lx+i*Ly as you will see below.

Warning: Do not use the Simplify feature at this time (why?). Also be aware that you have to deal with more than one F function; therefore label them properly, e.g. use F40 to denote the spherical tensor operator F corresponding to l=4 and m=0.

2. Now you will need some knowledge how to do matrix algebra on `Mathematica'. Here they are
		(i) Define the function d[x] which displays a matrix x in matrix form

d[x_]:=MatrixForm[x]

(ii) Define the 5×5 zero matrix O (you may find useful to type in ?Table)  and display it

o=Table[0,{i,5},{j,5}];d[o]

(iii) Define the lp=L+ and lm=L- operators and display them (you may find it useful to type in       
?Transpose)

lp=o;lp[[1,2]]=2;lp[[2,3]]=Sqrt[6];
lp[[3,4]]=Sqrt[6];lp[[4,5]]=2;d[lp]

lm=Transpose[lp];d[lm]

(iv) Define lx=Lx and ly=Ly

lx=(lp+lm)/2;d[lx]

ly=(lp-lm)/2/I;d[ly]

(v) Define lz=Lz

lz=DiagonalMatrix[2,1,0,-1,-2];d[lz]

3. In order to evaluate C[2,l] for given l just use directly eq.(7) from the text of the problem for the particular case m=0 and M=M'=0. For evaluating the matrix element on the left hand side in eq.(7) you may use the command

Integrate[Expand[Yc[2,0,u,v]*Y[l,0,u,v]*Y[2,0,u,v]*Sin[u]],
{u,0,Pi},{v,0,2*Pi}]

(c)

1. Next step is to obtaine the matrix form of V. You can do this simply by puting things together. Remember that now lx etc. are already defined and therefore by typing e.g. F[lx,ly,lz] the output will be a matrix. At this stage the Simplify command might be useful. Also, it is not a bad idea to divide the matrix V with a = (e2<r4>)/R5, i.e. to measure the energy in units of a.

2. All what you have left to do is to diagonalize V by using the command Eigenvalues[V]. You can also determine the corresponding eigenvectors by typing Eigenfunctions[V].