Clebsch.html

Clebsch-Gordan coefficients

The calculation of vector coupling coefficients is based on representing J_total -sqared = J ^2 in the basis of the single-particle angular momentum eigenstates and diagonalization. We make use of the matrix elements of the raising and lowering operators to assemble the matrix.

<j,m+1| J+ |j,m> = <j,m| J- | j,m+1> = -h_*sqrt((j-m)*(j+m+1))

We are using the decomposition J _tot^2 = J1 ^2 + J2 ^2 + 2 J1 _z * J2 _z + J1 _ * J2+ + J1 + * J2_

When assembling the matrix representation we keep in mind that the 5 terms on the RHS of the above equation separate into two groups:

the first three terms will contribute only on the diagonal, and the last two terms only one line away from the diagonal.

Suppose we want to couple J1 =3/2 and J2 =1. We work in units where h-bar h_ = 1 .

> restart; interface(warnlevel=0): with(linalg):

> j1:=3/2; j2:=1;

Normally, we are interested in picking an allowed valueof total angular momentum to which j1 and j2 couple (i.e., out of the range |j1-j2| .. (j1+j2) ).

However, the way things will work in the approach where J ^2 is diagonalized in the |j1 j2 m1 m2> basis, we will not use the specification of J , but rather find the CG coefficients for all possible J values. To leave the z -projection of total J as unrestricted as possible we determine the maximum allowed length:

> J:=j1+j2;

Our task is to represent J ^2 in the |j1 j2 m1 m2> basis with j1 and j2 specified through their quantum numbers. This means that the matrix representation of J ^2 can be labeled by m1 and m2 . However, m1 and m2 do not vary independently, but are restricted by the condition m1 + m2 = M .

Therefore, we now pick a projection out of the range M=-J..J . Let us pick anyone but the largest (or lowest) so that we will get some coupling. The reason for avoiding the largest projection is that only one combination of z -projections of J1 and J2 works in that case.

> M:=-J+j2+1;

The solution of the J ^2 eigenvalue/eigenvector problem will give us:

(i) through the eigenvalues the allowed quantum numbers of J ^2, i.e., h_^2 J (J+1) , with the possible choices of J that are compatible with the chosen M .

(ii) in anticipation of half-integer spins we pick the M value above such that it is compatible with J , i.e,. for integer J we pick integer M , and for half-integer J we pick half-integer M .

(iii) the eigenvector for a given J ^2-eigenvalue gives us the set of coupling coefficients that says how the uncoupled states combine to form an eigenvector of J ^2.

Now the question is: how do we find the correct range of (m1,m2) values by use of restrictions (such as M=m1+m2 ) ? This will determine the matrix size.

We determine by a counting method the matrix size in the variable ic , and the configurations are stored in cf :

> ic:=0: for m1 from -j1 to j1 do: for m2 from -j2 to j2 do: if m1+m2=M then ic:=ic+1: cf[ic]:=[m1,m2]: print(`combo `,ic,` involves [m1,m2]= `,cf[ic],` with m1+m2=M: `,M); fi; od: od: ic;

Now we understand that j2 (the smaller of the two) steps through all allowed values:

> Jsq:=matrix(ic,ic,0);

The matrix is a square matrix: we let the row and column indeces go through all allowed values of m1,m2 , i.e., each index goes from 1 to ic . We call them m1,m2 and m1p,m2p for row and column indices respectively.

We also need indeces to refer to the matrix entries (ranging from 1 to ic ) called im2 and im2p .

To understand the diagonal matrix entries is straightforward. They are simply the sum of the diagonal elements of the operators J1^2+J2^2+ 2*J1z*J2z , i.e., of the eigenvalues of these operators.

The off-diagonal elements require one to understand the following structure:

when m2=m2p+1 it follows from m1=M-m2=M-m2p-1 and m1p=M-m2p=m1+1 that m1=m1p-1 .

Vice versa: when m2=m2p-1 then m1=m2p+1 . This allows J1+J2- and J1-J2+ to have entries just above and below the diagonal.

> for irow from 1 to ic do: m1:=cf[irow][1]: m2:=cf[irow][2]:

> for icol from 1 to ic do: m1p:=cf[icol][1]: m2p:=cf[icol][2]:

> if m2=m2p then Jsq[irow,icol]:=j1*(j1+1)+j2*(j2+1)+2*m1*m2:

> elif m1=m1p-1 and m2=m2p+1 then Jsq[irow,icol]:=sqrt((j2-m2+1)*(j2+m2)*(j1+m1+1)*(j1-m1)): print(m1,m1p,m2,m2p,Jsq[irow,icol]);

> elif m1=m1p+1 and m2=m2p-1 then Jsq[irow,icol]:=sqrt((j1-m1+1)*(j1+m1)*(j2+m2+1)*(j2-m2)): print(m1,m1p,m2,m2p,Jsq[irow,icol]);

> fi: od: od:

> print(Jsq);

To find out the allowed total angular momentum for which the chosen M sublevel can be realized we calculate the eigenvalues and keep in mind that they follow h_^2*J*(J+1) .

> eigenvalues(Jsq);

Among the three possible eigenvalues of J ^2 we do find the previously selected maximum configuration J=j1+j2 , as well as the minimum J=j1-j2 , and one in-between.

> J*(J+1);

> abs(j1-j2)*(abs(j1-j2)+1);

We obtain the mixing coefficients that tell us the following:

the uncoupled basis states labeled by [j1 m1] and [j2 m2] can be coupled to [J M] values specified by fixed M and by calculated J (from the diagonalization of J ^2). The meaning of the eigenvector entries is that the coeffcients tell how much of the product of [j1 m1] and [j2 m2] states is needed with m1,m2 given by the specification in table cf . Thus the first eigenvector component is the admixing coefficient for the (m1, m2) values given in cf[1] , the second eigenvector entry is the one for the product of [j1 m1] with [j2 m2] with m1 and m2 given in cf[2] , etc.

> vec:=eigenvects(Jsq);

$vec := [15/4, 1, {vector([6^(1/2), 1, -2*2^(1/2)])}], [3/4, 1, {vector([-1/2*6^(1/2), 1, -1/2*2^(1/2)])}], [35/4, 1, {vector([1/6*6^(1/2), 1, 1/2*2^(1/2)])}]$

The eigenvector corresponding to the chosen value of J when normalized gives the coupling coefficients.

> nops([vec]);

> for i from 1 to nops([vec]) do: if vec[i][1]=J*(J+1) then myvec:=op(vec[i][3]); print(`unnormalized eigenvector: `,myvec); fi: od:

> CG:=map(combine,normalize(myvec));

The normalization is important to maintain unitarity, and thus a probabilistic interpretation. If we prepare the system in a [J, M, j1, j2] eigenstate, and then measure j1_z and j2_z , then the square of the CG coefficient will provide us the probabilities to measure a particular set of allowed (m1, m2) values. First we check:

> add(CG[i]^2,i=1..vectdim(CG));

> print(cf);

> CG[1]^2,CG[2]^2,CG[3]^2;

The statement is that when we are in the |J=5/2, M=-1/2, j1=3/2, j2=1 > eigenstate, then the probabilities are:

10% to be in the (m1,m2) = (-3/2, 1) combination

60% to be in the (m1,m2) = (-1/2, 0) combination

30% to be in the (m1,m2) = (1/2, -1) combination.

Exercise1 :

Go through the calculation for the other allowed M values and the same J value, i.e., J=j1+j2 . Calculate the probabilities to find the system in a particular (m1,m2) configuration and display graphically using arrows for the (j1,m1) and (j2,m2) angular momentum vectors how the (J M) state comes about.

Exercise2 :

Look at the other J values for which the set of (j1 m1) and (j2 M-m1) configurations couples to (J M) and construct the probabilities as well as the graphical representations discussed in exercise 1.

Exercise 3:

Pick another combination of j1 and j2 , and repeat the calculations for exercises 1 and 2.

A program for Clebsch-Gordan coefficients

One can arrive at a general formula given, e.g., in William J. Thompson: Angular Momentum , Wiley 1994, eq. 7.51. One introduces first the related (3j) coefficient for coupling three angular momenta to zero. It involves a sum over k (which is carried out as a sum over 2*k to allow half-integer angular momentum projection).

> ThreeJ:=proc(j1,m1,j2,m2,j3,m3) local k,tk,tkmin,tkmax,res,phas,n1,n2,d1,d2,d3,term;

> if m1+m2+m3 <> 0 then RETURN(0) fi;

> tkmin:=2*max(0,j2-j1-m3);

> tkmax:=2*min(j3-m3,j3-j1+j2); res:=0;

> phas:=(-1)^(tkmin/2+j2+m2);

> for tk from tkmin to tkmax by 2 do: k:=tk/2;

> n1:=(j2+j3+m1-k)!; #print("n1= ",n1);

> n2:=(j1-m1+k)!; #print("n2= ",n2);

> d1:=k!*(j3-j1+j2-k)!; #print("d1= ",d1);

> d2:=(j3-m3-k)!; #print("d2= ",d2,k,j1,j2,m3);

> d3:=(k+j1-j2+m3)!; #print("d3= ",d3);

> term:=phas*n1*n2/(d1*d2*d3);

> phas:=-phas;

> res:=res+term; od;

> simplify(res)*(-1)^(j1-j2-m3)*sqrt((j3+j1-j2)!*(j3-j1+j2)!*(j1+j2-j3)!*(j3-m3)!*(j3+m3)!/((j1+j2+j3+1)!*(j1-m1)!*(j1+m1)!*(j2-m2)!*(j2+m2)!)); end:

> CGC:=(j1,m1,j2,m2,J)->combine(ThreeJ(j1,m1,j2,m2,J,-m1-m2)*sqrt(2*J+1)*(-1)^(j1-j2+m1+m2));

For our example we have j1 j2 and the range of m2 and m1 :

> j1,j2,[-1,0,1],[M-(-1),M,M-1]; #the matrix was done by stepping with m2 [-1,0,1].

>

> CG[1];

> CGC(j1,M+1,j2,-1,J);

> %-%%;

> CG[2];

> CGC(j1,M,j2,0,J);

> simplify(%-%%);

> CGC(j1,M-1,j2,1,J);

> CG[3];

> %-%%;

In a nutshell the meaning of the calculation is:

|JM j1 j2> = |5/2, -1/2, 3/2, 1> = add( CG[M-m2,m2]* |j1,M-m2> |j2 m2> ,m2=-1..1)

> J,M,j1,j2;

There are sum rules for CG-coefficients given in Thompson's book (or other books on group theory/angular momentum algebra), which can be checked by explicit calculations.

Application

One useful application of CG coefficients or (3j) symbols is in the evaluation of angular integrals over three spherical harmonics. W.J. Thompson lists it as eq. (7.102) and shows how it follows as a special case of an integral over three Wigner rotation matrices. We start with a definition of the spherical harmonics.

> with(orthopoly);

> Plm:=proc(theta,l::nonnegint,m::integer) local x,y,f;

> x:=cos(theta);

> if m>0 then f:=subs(y=x,diff(P(l,y),y$m));

> else f:=subs(y=x,P(l,y)); fi;

> (-1)^m*sin(theta)^m*f; end:

> Plm(theta,3,2);

For the spherical harmonics we don't need the Plm's with negative argument.

> Y:=proc(theta,phi,l::nonnegint,m::integer) local m1;

> m1:=abs(m); if m1>l then RETURN("|m} has to be <= l for Y_lm"); fi;

> exp(I*m*phi)*Plm(theta,l,m1)*(-1)^m*sqrt((2*l+1)*(l-m1)!/(4*Pi*(l+m1)!)); end:

> Y(theta,phi,3,0);

> Y(theta,phi,3,1);

> No:=(l,m)->int(int(evalc(Y(theta,phi,l,m)*conjugate(Y(theta,phi,l,m))),phi=0..2*Pi)*sin(theta),theta=0..Pi);

> No(1,0);

Now suppose we are interested in the following integral:

> AI:=(J,M,j1,m1,j2,m2)->int(int(evalc(Y(theta,phi,j1,m1)*Y(theta,phi,j2,m2)*conjugate(Y(theta,phi,J,M))),phi=0..2*Pi)*sin(theta),theta=0..Pi);

There are restrictions on the angular momentum parameters: (i) triangle rule on (j1,j2,J) ; (ii) j1+j2+J = even integer; (iii) m1+m2=M . For example:

> AI(2,1,1,0,1,1);

In terms of (3j) symbols the integral can be calculated as:

> AI3j:=(J,M,j1,m1,j2,m2)->simplify((-1)^M*sqrt((2*j1+1)*(2*j2+1)*(2*J+1)/(4*Pi))*ThreeJ(j1,m1,j2,m2,J,-M)*ThreeJ(j1,0,j2,0,J,0));

> AI3j(2,1,1,0,1,1);

Exercise 4:

Verify the short-cut formula by comparing it against the symbolic computation of the integral. Go to large values of the angular momentum parameters and check its computational advantage.

Integrals of this type arise in atomic physics in the computation of multipole expansions of the potential of a charge distribution, and in the computation of dipole matrix elements (e.g., in the radiative decay) in which the Cartesian components of the position vector are represented through r*Y(1,m) , e.g., x+I*y=r*Y(1,1) .

>