Unprotect[Lp]
Needs["Utilities`CleanSlate`"]
CleanSlate[]
Remove["Global`*"]

Angular momentum

Cell 1: Define the function cp[j,m] which gives 

	(SubPlus[C])(j,m) = ℏ(Sqrt[j(j+1)-m(m+1)])
	
and  the function cn[j,m] which gives

	(SubMinus[C])(j,m) =(Sqrt[j(j+1)-m(m-1)])

Clear[cp, cm]
cp[j_, m_] := ℏ*Sqrt[j*(j + 1) - m*(m + 1)]
cm[j_, m_] := ℏ*Sqrt[j*(j + 1) - m*(m - 1)]

Test that they work properly. The TableForm function makes the  display look nicer.

{cp[j, m],cm[j,m]}//TableForm[#,TableSpacing->{1,5}]&

Check that (SubPlus[C])(j,j) = (SubMinus[C])(j,-j) =0

{cp[j, j], cm[j, -j] // FullSimplify }

Note that if m > j (which is not possible in physics), the function cp becomes imaginary. In the definition of the function there is no safeguard to prevent the user from using m > j. If you do this by mistake you will get an imaginary output that should tip you off.

cp[2,3]

Cell 2:  Define the operator (SubPlus[(OverHat[L])]) acting on |j,m⟩. Use the equation:

	(SubPlus[(OverHat[L])]) |j,m⟩ = (SubPlus[C])(j,m)|j, m+1⟩
	
with 

	(SubPlus[C])(j,m) = ℏ(Sqrt[j(j+1)-m(m+1)])
	
Note that  (SubPlus[(OverHat[L])]) |j, j⟩ = (SubPlus[C])(j,j)|j, m=j+1⟩ = 0 because (SubPlus[C])(j,j) = 0


For |j, m⟩ I use the symbol ket[j,m].  

(SubPlus[(OverHat[L])])|j, m⟩ is given by Lp[ket[j,m]]

For |j,m⟩ use the symbol ket[j,m].  
The function Lp[ket[j,m]] produces (OverHat[L])|j,m⟩ and should be used to act only on an argument of the form 
	a ket[j,m] + b ket[j1,m1]+...
where a, b can be any expressions.  In other words, it is designed to act correctly on a coherent superposition of kets  (denoted ket[a,b] and representing |a,b⟩ )

Below I give rules that define how Lp (which is Mathematica for (SubPlus[(OverHat[L])]) ) acts on a function of the form ket[j,m] (which is Mathematica for |j,m⟩)

Unprotect[Lp];
Clear[Lp];
 Lp::badArg=" Abs[m]=`1` is larger than j= `2` ";

Lp[ket[j_,m_]]:=Module[{cp,out },
If[ m>j||m<-j, Message[Lp::badArg,m,j], ℏ*Sqrt[j*(j + 1) - m*(m + 1)]*ket[j,m+1] ,ℏ*Sqrt[j*(j + 1) - m*(m + 1)]*ket[j,m+1]]
      ] (*if m and j are symbols the first argument in If does not evaluate to True or False and the fourth argument  in If is evalauated*)

(* define Lp to factor scalars that multiply  the ket*)
Lp[ a_*ket[j_, m_]] := a Lp[ket[j, m]]/;FreeQ[a,ket]

(*Make Lp distributive when the argument is a sum*)

Lp[a_+ b_] :=Lp[a]+Lp[b]

Mathematica explanation. I wrote the function so that it gives me an error message if I use an erroneous argument in ket[j,m] such as m < - j or m > j.

To make the output look like the equations in the lecture, use this replacement rule:

kR=ket[j_,m_]:>StringForm["|``,``⟩",j,m]

Test that this function for (OverHat[(SubPlus[L])]) works properly

Lp[ket[j,m]]/.kR

Lp[a ket[2,1] + b ket[3,2]]/.kR

Lp[ket[2,3]]

Lp[a ket[j,m] + b ket[1,ms]]/.kR

Lp[ket[j,m] + b ket[1,ms]]/.kR

Lp[a ket[j,m] +  ket[1,ms]]/.kR

Note that if ms > 1, the ket |1, 1+ms⟩ should be zero. We have not provided a rule that implements this automatically. It is not possible to implement such a rule for a symbolic argument ms, since there is no way to tell whether ms is negative or positive. If j and m have numerical values one can implement such rules and tie them to the function ket.

Check that (Subscript[(OverHat[L]), + ])(expression |j,m⟩) = expression (Subscript[(OverHat[L]), + ])|j,m⟩ for a few expressions

Lp[(a+Sqrt[b])ket[j,m]]/.kR

{Lp[(a+Sqrt[b])ket[j,m]]/.kR,Lp[(x^2+3) ket[j,m]]/.kR  }//TableForm[#,TableSpacing->{5,0}]&

Check that it gives correct results for the electron spin:

Lp[ket[1/2,-1/2]]/.kR

Check that  (Subscript[(OverHat[L]), + ])| j, j⟩ = 0. This is automatically ensured because (SubPlus[C])(j,j) = 0

Lp[(x^2+3) ket[j,j]]/.kR

Act on a sum of three kets:

Lp[a ket[j,m]+(x+y)^2 ket[j1,m1]+ d ket[j,m]] /.kR

Cell 3: Define the operator (SubMinus[(OverHat[L])]) acting on |j,m⟩.  

 	(SubMinus[(OverHat[L])]) |j, m⟩ = (SubMinus[C])(j,m)|j, m-1⟩
with
	(SubMinus[C])(j,m) =(Sqrt[j(j+1)-m(m-1)])

Note that 
	(SubMinus[(OverHat[L])]) |j, -j⟩ = 0

because (SubMinus[C])(j,-j)=0 

For |j,m⟩ I use the symbol ket[j,m].  

(SubMinus[(OverHat[L])])|j,m⟩ is given by Lm[ket[j,m]]

The Mathematica syntax for Lm is the same as that for Lp

Unprotect[Lm];
Clear[Lm];
 Lm::badArg=" Abs[m]=`1` is larger than j= `2` ";
Lm[ket[j_,m_]]:=Module[{cp,out },
If[ m>j||m<-j, Message[Lm::badArg,m,j],ℏ*Sqrt[j*(j + 1) - m*(m - 1)]*ket[j,m-1] ,ℏ*Sqrt[j*(j + 1) - m*(m - 1)]*ket[j,m-1]]
      ]

(* define Lm to act on a ket and factor scalars that multiply  the ket*)
Lm[a_*ket[j_, m_]] := a*Lm[ket[j, m]]/;FreeQ[a,ket]

(*Make Lp distributive over the sum*)

Lm[a_+b_] := Lm[a]+Lm[b]

Perform the type of tests you used for (SubPlus[(OverHat[L])])

Lm[ket[j,m]]/.kR

Lm[a ket[j,m]]

Lm[(a+Sqrt[x]) ket[y,z]]

{Lm[ a ket[j,-j]] , Lm[ a ket[j,-j]] //FullSimplify/.kR}//TableForm[#,TableSpacing->{5,0}]&

Note that if you do not use FullSimplify to simplify the expression, you might think that it is not equal to zero

Lm[ a ket[j,-j]] /.kR

Lm[Sqrt[x] ket[j,m1]]/.kR

Clear[x,a]
Lm[a ket[j,-j]+Sqrt[x] ket[j,j]]/.kR

The result should be simplified because the term multiplying a is zero but the square root has not been evaluated

Lm[a ket[j,-j]+Sqrt[x] ket[j,j]]/.kR//FullSimplify

Cell 4: Define the operators (Subscript[(OverHat[L]), x]), (Subscript[(OverHat[L]), y]),(Subscript[(OverHat[ L]), z]) acting on |j,m⟩.  Use

	(Subscript[(OverHat[L]), z]) |j,m⟩ = m ℏ |j,m⟩

 	
(Subscript[(OverHat[L]), x]) =(((SubPlus[(OverHat[L])])+(SubMinus[(OverHat[ L])]))/2)
 and 
 	
(Subscript[(OverHat[L]), y]) =(((SubPlus[(OverHat[L])])-(SubMinus[(OverHat[L])]))/(2i))
 
 with (Subscript[(OverHat[L]), + ]) and (SubMinus[(OverHat[L])])  defined in Cells 2 and 3.

Unprotect[Lz];
Clear[Lz];
Lz::badArg=" Abs[m]=`1` is larger than j= `2` ";
Lz[ket[j_,m_]]:=  Module[{}, If[ m>j||m<-j, Message[Lp::badArg,m,j],ℏ*m*ket[j,m] ,ℏ*m*ket[j,m]]  ]

(*define Lz to act on ket and factor any expression that multiplies the ket*)
Lz[a_*ket[j_, m_]] := a*Lz[ket[j, m]]

(*Make Lm distributive over the sum*)

Lz[a_+b_] :=  (Lz[a]+Lz[b])//Expand 
Protect[Lz];

Test that Lz  works correctly

Lz[ket[j,m]]/.kR

Lz[a b c  ket[j,m]]/.kR

Lz[a ket[j,m]+ b c ket[j,m]]

Lz[a  b ket[1,-1]]

Define Lx (the operator  (Subscript[(OverHat[L]), x])) and Ly (the operator  (OverHat[(Subscript[L, y])]))

Unprotect[Lx];
Clear[Lx]

Lx[x_]:=(1/2)*(Lp[x]+Lm[x])//Expand
Protect[Lx];

Unprotect[Ly];
Clear[Ly]
Ly[x_]:=(1/(2 I))*(Lp[x]-Lm[x])//Expand
Protect[Ly];

Tests:

Lx[ket[j,m]]/.kR

Ly[ket[j,m]]/.kR

Do not use the syntax below because ket is not defined to give zero when m is out of range

Lx[ket[j,m]]/.kR/.{j->1/2,m->1/2}

That worked without problems but look at the result obtained below

Lx[ket[j,m]]/.kR/.{j->1/2,m->3/2}

The ket in the output should be zero but it is not because we did not provide rules for that.

Lx[ket[j,m]]/.kR

Do some more tests

Lx[ket[1/2,1/2]]

Ly[ket[1/2,-1/2]]/.kR

Lx[ket[1,-1]]/.kR//FullSimplify

Ly[ket[1,-1]]/.kR//FullSimplify

Lx[ket[j,j]]//FullSimplify

Lx[ket[1,0]]/.kR

Ly[ket[1,0]]/.kR

Lx[a ket[j,m]+b ket[j1,m1]]/.kR

Lx[ket[1,1]]/.kR

Cell 5: Define the scalar product involving angular momentum kets

Define the function
	
	scalar[ket[j,m],ket[j',m']]

that calculates 

	⟨j,m|j',m'⟩

Clear[scalar]
(*a rule which makes the scalar product zero if the quantum numbers are different and 1 if they are the same: *)
scalar[ket[j_,m_],ket[jp_,mp_]]:=Which[j===jp&&m===mp, 1,True,0];
(*a rule that factors a scalar from the first argument of the scalar product and complex conjugates the scalar:  *)
scalar[a_*ket[j_,m_],  ket[jp_,mp_]]:= Conjugate[a]* scalar[ket[jp,mp],ket[j,m]]/;FreeQ[a,ket];  
(*a rule that factors a scalar from the second argument of the scalar product*)
scalar[ ket[j_,m_], b_*ket[jp_,mp_]]:= b*scalar[ket[jp,mp],ket[j,m]]/;FreeQ[b,ket];  ;  
(* a rule that factors scalar from both arguments*)
scalar[ a_* ket[j_,m_], b_*ket[jp_,mp_]]:= Conjugate[a]*b*scalar[ket[j,m],ket[jp,mp]]/;FreeQ[a,ket]&&FreeQ[b,ket];  
(* rules that factors expressions of scalars*)
 scalar[a_*ket[j,m]+b_*ket[j1,m1],c_*ket[j2,m2]]:=Conjugate[a]*c*scalar[ket[j,m],ket[j2,m2]]+ Conjugate[b]*c* scalar [ket[j1,m1],ket[j2,m2]];
scalar[ket[j_,m_],0]:=0

scalar[0,ket[j_,m_] ]:=0

scalar[a_+b_,c_]:= scalar[a,c]+scalar[b,c];

scalar[a_,b_+c_]:= scalar[a,b]+scalar[a,c];

Run tests

scalar[ket[j,m],ket[j,m]]

scalar[ket[j,m],ket[j,m1]]

scalar[ket[j,m],ket[j1,m]]

scalar[a ket[j,m],b ket[j,m]]

scalar[b ket[j,m]+a ket[j,m],c ket[j,m]]

scalar[d ket[j,m]+ e ket[j,m],a ket[j,m]+c ket[j,m]]

scalar[a ket[j1,m1],b ket[j1,m1] ]

scalar[d ket[j,m],b ket[j,m]+c ket[j1,m1]]

scalar[ket[j,m],bx ket[j,m]]

scalar[(ax  +I 3)ket[j,m],  ket[j,m]]

scalar[(ax  +I 3)ket[j,m],  b ket[j,m]]

scalar[Sqrt[a+I b] ket[j,m],ket[j,m]]

scalar[a ket[j,m],b ket[j,m]]

scalar[ket[j,m],Sqrt[a+ I b] ket[j,m1]]

The function scalar should be usable for calculating the matrix elements 

	⟨j, m |  (OverHat[O]) | (Subscript[j, 1]),(Subscript[m, 1])⟩
	
where (OverHat[O]) is any of the operators (Subscript[L, x]) , (Subscript[L, y]), (Subscript[L, z]), (L^2), (SubPlus[L]), or  (SubMinus[L]).

Let us find some matrix elements.

Find the matrix element ⟨j,m|(SubPlus[L])| j, m-1⟩. I calculate both (SubPlus[L])| j, m-1⟩ and ⟨j,m|(SubPlus[L])| j, m-1⟩ to compare them and make sure that they both work properly

{Lp[ ket[j,m-1]]/.kR,scalar[ket[j,m],Lp[ ket[j,m-1]]]/.kR}

See below more examples

Lp[ket[j,m]]
scalar[ket[j,m],Lp[ ket[j,m]]]

In the second term below we have a more complex situation, which scalar handles well.

{Lp[ket[j,m-1]],
scalar[(2+I 3) Sqrt[g] ket[j,m],c Sqrt[a+b] Lp[ket[j,m-1]]]}

Lz[ket[j,m]]/.kR

Lp[ket[j,m-1]]/.kR

{Lz[ket[j,m]]/.kR,
Lp[ket[j,m-1]]/.kR,
scalar[(2+I 3) Sqrt[g]  Lz[ket[j,m]],c Sqrt[v] Lp[ket[j,m-1]]]/.kR}//MatrixForm

scalar[(2+I 3) *ket[j,m],c *ket[j,m]]

Cell 6: Calculate the matrices of (SubPlus[(OverHat[L])]) , (SubMinus[(OverHat[L])]), (Subscript[(OverHat[L]), x ])and (Subscript[(OverHat[L]), y]) for a given value of j using

	⟨j, m| (OverHat[O])|j,m'⟩
	
for the matrix elements of the operator (OverHat[O]), where (OverHat[O]) is one of (SubPlus[(OverHat[L])]) , or (SubMinus[(OverHat[L])]),or  (Subscript[(OverHat[L]), x ]),or(Subscript[(OverHat[L]), y]) or (Subscript[(OverHat[L]), z]) 

This means that we are calculating the matrix elements of these operators  in the basis set   |j,m⟩ , m=-j, -j+1,....j-1, j.  To avoid confusion we need to be consistent when we define vectors and matrices and use the same basis set for both. The order of the kets in the basis set matters and should be the same in all definitions.

In the matrix  ⟨j,m|O|j’,m’⟩ ≡ ⟨m|O|m’⟩  (we dropped j since it is not changing in the calculations below (see lecture notes)). The values of m label rows and m’ labels columns. 

Calculate the eigenvalues and eigenvectors of (Subscript[(OverHat[L]), x])

For j = 1, m = -1, 0, 1  the matrix of an operator is labeled as below

(({{ma[-1,-1], ma[-1,0], ma[-1,1]}, {ma[0,-1], ma[0,0], ma[0,1]}, {ma[1,-1], ma[1,0], ma[1,1]}}))
The components of a vector in this space are

	|ψ⟩ = (Subscript[ψ, 1])|-1⟩ + (Subscript[ψ, 2])|0⟩ + (Subscript[ψ, 3])|1⟩

Check that the line of code given below labels the matrix correctly

Table[ma[m,mp],{m,-1,1,1},{mp,-1,1,1}]//MatrixForm

1. Definitions of matrices for (SubPlus[(OverHat[L])]) , (SubMinus[(OverHat[L])]), (Subscript[(OverHat[L]), x ])and (Subscript[(OverHat[L]), y])

The matrix of (SubPlus[(OverHat[L])]) for arbitrary j

Clear[LpMatrix]
LpMatrix[j_]:=Table[scalar[ket[j,m], Lp[ket[j,mp]]],{m,-j,j,1},{mp,-j,j,1}]

LpMatrix[3]//MatrixForm

Define the matrix of (SubMinus[(OverHat[L])]) for arbitrary j

Clear[LmMatrix]
LmMatrix[j_]:=Table[scalar[ket[j,m], Lm[ket[j,mp]]],{m,-j,j,1},{mp,-j,j,1}]

Take a look at the matrix LmMatrix (this is (Subscript[L, x]))

LmMatrix[2]//MatrixForm

Define the matrices of (Subscript[(OverHat[L]), x]) and (Subscript[(OverHat[L]), y]) for arbitrary j

Clear[LxMatrix,LyMatrix]
LxMatrix[j_]:=(1/2)*(LpMatrix[j]+LmMatrix[j])
LyMatrix[j_]:=(1/(2 I))(LpMatrix[j]-LmMatrix[j])
LzMatrix[j_]:=Table[scalar[ket[j,m], Lz[ket[j,mp]]],{m,-j,j,1},{mp,-j,j,1}]

Take a look at them. Note that they are very sparse matrices.  Try to explain why.

LxMatrix[1]//MatrixForm

LxMatrix[2]//MatrixForm

LyMatrix[2]//MatrixForm

I calculate below the matrix of (Subscript[(OverHat[L]), z]). This is diagonal since the basis set consists of the eigenstates of (Subscript[(OverHat[L]), z])

LzMatrix[2]//MatrixForm

2. The eigenvalues and the eigenvectors of (Subscript[(OverHat[L]), x]).

{{lx[1],lx[2],lx[3]},{vx[1],vx[2],vx[3]}}=Eigensystem[LxMatrix[1]]

The eigenvalues are -ℏ, 0, ℏ.  They are the same as the eigenvalues of (Subscript[L, z]), as they should be.

These eigenstates are not normalized. Normalize them. Note that because there are no complex coefficients involved, I can use the scalar product (Dot) provided by Mathematica.

ma=Table[vx[i]=vx[i] /Sqrt[vx[i].vx[i ]],{i,1,3}]

I could have used the function Normalize.

The general theory says that the eigenvectors form an unitary matrix. Check that this is true.
The absolute value of the determinant of the matrix should be equal to 1.

Det[ma]//Abs

The inverse of the matrix should be equal to the transpose (true for any unitary matrix).

Inverse[ma]==Transpose[ma]

In the ket space these eigenstates are  given by the usual formula. If {(Subscript[x, 1]), (Subscript[x, 2]), (Subscript[x, 3 ])} is one of the eigenvectors then the eigenket is 
	|ψ⟩ = (Subscript[x, 1]) |1,-1⟩ + (Subscript[x, 2]) |1,0⟩ + (Subscript[x, 3]) |1,1⟩
Display below in a table the eigenvalues and the corresponding eigenkets.

kets={ket[1,-1],ket[1,0],ket[1,1]};
TableForm[Table[{lx[i],k[i]=vx[i].kets/.kR},{i,1,3}], TableHeadings->{None,{"eigenvalue Lx","     eigenvector"}},TableSpacing->{3,10}]

Cell 7:  The eigenvalues and the eigenvectors of the Hamiltonian
	
	(OverHat[H])= (Subscript[(OverHat[H]), 0]) + (OverHat[V]) 
with	
	(OverHat[V]) = -(e/2 (Subscript[m, e])) B (Subscript[(OverHat[L]), x]) 
	
in the basis  
	
	{|(Subscript[L, z]);2, 0,0⟩, |(Subscript[L, z]);2, 1,-1⟩, |(Subscript[L, z]);2, 1,0⟩, |(Subscript[L, z]);2, 1,1⟩}
		
Here e is the charge of the proton and B is the magnetic field in the direction OX.  In the symbol  |(Subscript[L, z]);n, j,m⟩, n indicates that this is an eigenstate of the Hamiltonian (Subscript[H, 0]) of the hydrogen atom for a state of the energy  (-Ry/(n^2)); |(Subscript[L, z]);n, j,m⟩ is also an eigenstate of the angular momentum squared with eigenvalue (ℏ^2)j (j+1) ; |(Subscript[L, z]);n, j,m⟩ is also an eigenstate of (Subscript[(OverHat[L]), z ])with the eigenstate mℏ.  

The presence of (Subscript[L, z ])in |(Subscript[L, z]);n, j,m⟩  reminds  you that the ket is an eigenstate of (Subscript[(OverHat[L]), z]).   Since we work in the subspace generated by n=2 and the eigenstates of (Subscript[(OverHat[L]), z]) we don't really need to keep (Subscript[L, z ])and 2 in |(Subscript[L, z]);2, j,m⟩ and write only |j,m⟩.  This will be represented in Mathematica by ket[j,m]

The matrix (Subscript[H, 0]) is diagonal in this basis set so  we only need the matrix elements of the perturbation proportional to (Subscript[L, x]).

The matrix representing (OverHat[V])=-((e (Subscript[B, x]))/(2 (Subscript[m, e])))(OverHat[(Subscript[L, x])]) in the basis set  { | (Subscript[L, z]);1,-1⟩,|(Subscript[L, z]):1,0⟩, |(Subscript[L, z]);1,1⟩}

Clear[ket, mat]
LxMa=Table[scalar[ket[1,m],Lx[ket[1,mp]] ], {m,-1,1,1},{mp,-1,1,1}]

LxMa//MatrixForm

The matrix of (Subscript[H, 0]) is diagonal and has the values - Ry/4 as diagonal elements. Generate it

H0= IdentityMatrix[3]*(- Ry/4)

H0-(e B/(2(Subscript[m, e])))LxMa//MatrixForm

Find the eigenvalues and the eigenvectors of (OverHat[H])

Eigenvalues[ H0-(e B/(2(Subscript[m, e])))   LxMa]//FullSimplify//MatrixForm

{eigenVal,eigenVec}=Eigensystem[ H0-(e B/(2(Subscript[m, e])))   LxMa]//FullSimplify

kets={ket[1,1],ket[1,0],ket[1,-1]};
TableForm[Table[{eigenVal[[i]],k[i]=eigenVec[[i]].kets/.kR},{i,1,3}], TableHeadings->{None,{"eigenvalue Lx","     eigenvector"}},TableSpacing->{3,10}]

I did not bother to check if the eigenvectors are normalized.  

The point here is that the eigenvalues do not depend on the fact that the field is oriented along the OX direction or along the OZ direction.  However, since the basis set consisted of the eigenstates of (Subscript[L, z]) the calculation would have been much simpler if we took the Z-axis along the field.

The addition of two angular momenta 
	could be the sum of two spins
or 
	could be the sum of two orbital angular momenta
or 
	could be a spin plus an angular momentum.