Twoーelectron states 1

Quantum Mechanics
                     K. Konishi, G. Paffuti

Use the Clebsch-Gordan coefficients to write directly the states with definite angular momentum and spin for a system of two electrons. Write the form of one- and two- electron operators. In particular consider the electrostatic interaction, the LS coupling and the interaction with an external magnetic field.

Warning: This notebook has several global variables. It is better to restart from this point if some parameters are changed.

Global variables and initialization

Unprotect[Lorb,Lpart,partEquiv,p$p];
ClearAll["Global`*"];
Off[General::spell, General::spell1];
Off[ClebschGordan::"phy",ClebschGordan::"tri"];
<<Combinatorica`
Protect[Lorb,Lpart,partEquiv,p$p];

Two electron states

The notation used in this notebook is the same as in notebooks [NB3, NB4].

Kinematics

We shall use the notation

v[{(Subscript[m, 1]),(Subscript[S, 1]),(Subscript[m, 2]),(Subscript[S, 2])}]

a state with fixed projections along z for L and S of each electron, i.e. a Slater determinant. (Subscript[m, i]),(Subscript[S, i])are the quantum numbers of the i-th electron. The list in v[] is ordered and this fix the sign of the determinant. The quantum numbers are understood completely antisymmetric, in agreement with Pauli principle. The basis vectors are assumed with unit norm.

Simple operations

In this notebook we will consider both equivalent and non equivalent electrons and we need some tools to distinguish between the two cases and some simple linear operators acting on them. Here we list some functions coded below and briefly describe their use

Assigning orbitals

The notebook is intended to cover different cases, so a function is provided which transforms a conventional input for electrons in quantum numbers:

fixParticles[orb1,orb2]

This function assigns angular momenta and determines if the two electrons are or not equivalent. It is useful to set some global variables Lorb , partEquiv,  Lpart, p$p which refer respectively ...  . The function works both with a single particle input or two particles. The input can be ["s","p"], or ["1s","2p"] etc.

The construction of L, S eigenstates

The eigenstates of L and S are simply constructed with the Clebsch Gordan coefficients as

|L,Lz;S,Sz⟩ = (∑_(sz,lz))C[(1/2),sz,(1/2),(Subscript[S, z])-sz,S,(Subscript[S, z])]C[l1,l1z,l2,(Subscript[L, z])-l1z,L,(Subscript[L, z])]v[{l1,sz,(Subscript[L, z])-l1z,(Subscript[S, z])-sz}]

C[] are Clebsch Gordan coefficients.The addition rules for third components of L and S has been taken into account. By definition the Slater determinants v[] are antisymmetric. This construction obeys Pauli's principle. In the function coded below the Antisymmetry is ensured by applying orderVector. The same construction works for J eigenvectors, combining eigenvectors of L and S. We supply to functions

For two particles states there are two possibilities

The function allowedStates below give a list of possible states for Lz=0=Sz. These states are representative of the possible multipltes. Each element contain {L, S, vector}. The function uses arguments stored in the global variables, so the function fixParticles must be called before this one.

The similar function allowedJStates computes the admissible states in the subspace Jz = 0. This subspace, for two electrons, is representative for all possible fine structure states.

Code

Electron configuration

orderVector[x_]:= orderVector[Expand[x]];
orderVector[x_. v[a_] + y_.]:= x   orderVector[v[a]]+orderVector[y];
orderVector[v[k_]]:=Module[{vett,k1,signk},
If[partEquiv== True,
k1 = Partition[k,2];
signk = Signature[k1];
vett=signk*v[Flatten[Sort[k1]]];,
vett = v[k]];
vett];
orderVector[0]=0;

normalize[x_]:= Module[{cf,zzz},
 cf = Apply[List,Expand[x]+zzz]/.{v[a_]-> 1,zzz-> 0};
x/(Sqrt[cf.cf])]

mom["s"]=0;mom["p"]=1;mom["d"]=2;mom["f"]=3;

fixParticles[particelle_,seconda_:{}]:=Module[{ll,ord,cond,pp,bool},
(* we order angular momenta *)
pp =Flatten[{particelle,seconda}];
ll = mom/@(StringTake[#,-1]&/@pp);
ord = Ordering[ll];
ll = ll[[ord]];
pp = pp[[ord]];
If[Length[pp] <2,bool = False;,bool=(pp[[1]]== pp[[2]]);];
Unprotect[Lorb,Lpart,partEquiv,p$p];
partEquiv = bool;
Lorb = Table[0,{2*Length[pp]}];
Table[Lorb[[2i-1]]=ll[[i]];,{i,1,Length[pp]}];
Lpart = ll;
p$p = pp;
Protect[Lorb,Lpart,partEquiv,p$p];
];

allowedStates[]:=Module[{vett,totstati = {}},
Do[
If[(partEquiv== False)||EvenQ[L+S],
vett = eigVLS[L,0,S,0];
totstati = Append[totstati,{L,S,vett}];]
,
{S,1,0,-1},{L,Lpart[[2]]+Lpart[[1]],Lpart[[2]]-Lpart[[1]],-1}];
totstati]

allowedJStates[JzIn_:0]:=
Module[{Jz,vett,totstati = {}},
Jz = JzIn;
Do[
Do[
If[(partEquiv== False )|| EvenQ[L+S],
vett = eigVJ[J,Jz,L,S];
totstati = Append[totstati,{J,L,S,vett}];];
,
{J,L+S,Abs[L-S],-1}];
,
{S,1,0,-1},
{L,Lpart[[2]]+Lpart[[1]],Lpart[[2]]-Lpart[[1]],-1}];
totstati];

Angular momentum

oplz[x_.  v[k_]]:= x* Sum[k[[i]]*v[k],{i,1,Length[k]-1,2}]
oplz[0]=0; oplp[0]=0;oplm[0]=0;
oplm[x_.  v[k_]]:= 
x*Sum[(Sqrt[Lorb[[i]]*(Lorb[[i]]+1) - k[[i]](k[[i]]-1)])*
orderVector[v[ReplacePart[k,k[[i]]-1,i]]],{i,1,Length[k]-1,2}]
oplm[0]=0;
oplp[x_.  v[k_]]:= 
x*Sum[(Sqrt[Lorb[[i]]*(Lorb[[i]]+1) - k[[i]](k[[i]]+1)])*
orderVector[v[ReplacePart[k,k[[i]]+1,i]]],{i,1,Length[k]-1,2}]
oplp[0,x_]=0;
opLp[a_]:= Expand[Distribute[oplp[a]]]
opLm[a_]:= Expand[Distribute[oplm[a]]]
opLz[a_]:= Expand[Distribute[oplz[Expand[a]]]]
opL2[a_]:= opLp[opLm[Expand[a]]] + opLz[opLz[Expand[a]]] - opLz[Expand[a]]

opsz[x_.  v[k_]]:= x* Sum[k[[i]]*v[k],{i,2,Length[k],2}]
opsz[0]=0;
opsm[x_.  v[k_]]:= x*Sum[(Sqrt[(1/2)*((1/2)+1) - k[[i]](k[[i]]-1)])*
orderVector[v[ReplacePart[k,k[[i]]-1,i]]],{i,2,Length[k],2}]
opsm[0]:=0;
opsp[x_.  v[k_]]:= x*Sum[(Sqrt[(1/2)*((1/2)+1) - k[[i]](k[[i]]+1)])*
orderVector[v[ReplacePart[k,k[[i]]+1,i]]],{i,2,Length[k],2}]
opsp[0]=0;
opSp[a_]:= Expand[Distribute[opsp[a]]]
opSm[a_]:= Expand[Distribute[opsm[a]]]
opSz[a_]:= Expand[Distribute[opsz[Expand[a]]]]
opS2[a_]:= opSp[opSm[Expand[a]]] + opSz[opSz[Expand[a]]] - opSz[Expand[a]]

opJp[a_]:=opLp[a] + opSp[a];
opJm[a_]:=opLm[a] + opSm[a];
opJz[a_]:=opLz[a] + opSz[a];
opJ2[a_]:= opJp[opJm[Expand[a]]] + opJz[opJz[Expand[a]]] - opJz[Expand[a]];

Eigenstates of L, S and of J

eigVLS[L_,Lz_,S_,Sz_]:=Module[{l1,l2},
l1 = Lpart[[1]];l2 = Lpart[[2]];
Expand[normalize[Simplify[Sum[ClebschGordan[ {(1/2),sz},{(1/2),Sz-sz},{S,Sz}]* Sum[ClebschGordan[{l1,m},{l2,Lz-m},{L,Lz}]orderVector[v[{m,sz,Lz-m,Sz-sz}]],{m,-l1,l1}],{sz,-(1/2),(1/2)}]]]]
];

eigVJ[J_,Jz_,L_,S_]:=Expand[normalize[Simplify[Sum[ClebschGordan[{L,Lz},{S,Jz-Lz},{J,Jz}]eigVLS[L,Lz,S,Jz-Lz],{Lz,-L,L}]]]]

Dynamics

We are interested in two kind of operators : twoーparticle operators (electrostatic energy) and oneーarticle operators (LS interactions, interactions with a magnetic field). Below we describe briefly how these operators are constructed in the notebook.

Two particles operator

The matrix elements of a two particles operator between two Slater determinants is given by

(v[i,j]| (Subscript[H, I]) | v[r t]) = (ij| (Subscript[H, I]) |r t ) - (ij | (Subscript[H, I]) | t r)

This is a particular (and trivial) case of a general formula given in the text. We are interested in particular to the case of electrostatic interaction.

Using multipole expansion it is possible to write for this case (see text):

(ij| (Subscript[H, I]) |r t ) = (Subscript[δ, si sr])(Subscript[δ, sj st])(∑_k)(Subscript[c, k])[Li mi;Lr mr](Subscript[c, k])[Lt mt;Lj mj](R^k)[ij;rt]

where coefficients ck are a combination of Clebsch Gordan coefficients (written explicitly in the text and in the code below). Conservation of angular momentum and parity bound the sum over k to |(Subscript[L, a])-(Subscript[L, b])|≤ k ≤ (Subscript[L, a])+(Subscript[L, b]).

(R^k) are integrals of radial wave functions

(R^k)[ij;rt]= ∫(Subscript[P, i])[r1](Subscript[P, r])[r1](1/(Min[r1,r2]))(((((Min[r1,r2])/(Max[r1,r2]))))^k)(Subscript[P, j])[r2](Subscript[P, t])[r2]dr1dr2 .

We are considering matrix elements within a fixed electronic configuration, i.e. fixed radial wave functions. This means that, as a pair, (i,j) = (r,t) (apart the order). Two kind of functions can appear. If i=r (then necessarily j = t as the pairs are equal) the integral is named F, is i=t (then j=r) the integral is named G. If the two particles are equivalent (i.e. the same orbital) then clearly only F appear. As an instance R[s,p;s,p] = F[s,p]; while R[s,p; p, s] = G[s,p]. Usually G is the exchange integral.

Below we write a code for a function opEE which gives the formulas ( , ) and extend it to a linear combination of determinants. The functions radialIntegrals and ck define the corresponding coefficients (the first one only as symbols).

One-particle operators : LS interactions

Single particle operators between the same configuration have non vanishing operators only between the same state or between states which differ at most by one orbital. This is shown in general in the text,but it is trivial to verify for a two particle states, due to the orthogonality of the orbitals.

LS interactions

This is true in particular for LS interaction. With

(Subscript[H, LS]) = V((Subscript[r, 1])) (Subscript[L, 1])(Subscript[S, 1]) + V((Subscript[r, 2])) (Subscript[L, 2])(Subscript[S, 2])

the only non zero matrix elements are of the form

(v[m1 s1,m2 s2] | (Subscript[H, LS]) | v[m1 s1,m2 s2] ) = (Subscript[ξ, 1])(m1 s1|L S|m1 s1) + (Subscript[ξ, 2])(m2 s2|L S|m2 s2);
(v[m1 s1,m2 s2] | (Subscript[H, LS]) | v[m1 s1,k2 σ2] ) =  (Subscript[ξ, 2])(m2 s2|L S|k2 σ2);

ξ are radial integrals while the LS matrix elements can be easily computed with the decomposition

L S  = (1/2)(SubPlus[L])(SubMinus[S]) + (1/2)(SubMinus[L])(SubPlus[S]) + (Subscript[L, z])(Subscript[S, z])

The formula () is implemented below by the function matOpLS.

Magnetic interactions

The interaction with a magnetic field B has the form (B is the product of the magnetic field and of the Bohr magneton in this formula):

(Subscript[H, B]) = (∑_electrons)((Subscript[L, z]) + 2(Subscript[S, z]))B

The only non vanishing matrix elements are those diagonal on Slated determinants (eigenstates of single particle Lz and Sz) and for two particle one has immediately

(v[ij]|(Subscript[H, B]) | v[ij])= B((mi + 2si)+ (mj + 2 sj))

The formula () is implemented in the function matOpB below.

Code

funzRad[k_,x1_,x2_,x3_,x4_]:=
If[x1 == x3 || x1== x2,ris = (F^(ToString[k]))[x1,x2],(G^(ToString[k]))[x1,x2]];

c[k_,l1_,m1_,l2_,m2_]:= 
Quiet[Check[(((-1))^m2)(Sqrt[(((2l1+1)(2l2+1))/(((2k+1))^2))])*
ClebschGordan[{l1,-m1},{l2,m2},{k,m2-m1}]*
ClebschGordan[{l1,0},{l2,0},{k,0}],0]]

The operator is extended bilinearly on superpositions of v[].

opEE[x_. v[a_] + y_.,z_. v[b_] + t_.]:= x  * z opEE[v[a],v[b]]+z*opEE[y,v[b]] +
x opEE[v[a],t] + opEE[y,t];
opEE[0,x_]:= 0; opEE[x_,0]:= 0;

Let us remember that angular momenta are ordered as  L1 ≤ L2.

opEE[v[a_],v[b_]]:= Module[{L1,L2,mi,mj,mr,mt,si,sj,sr,st,elem1,elem2,fatt},
mi = a[[1]];mj = a[[3]];mr=b[[1]];mt=b[[3]];
si=a[[2]];sj=a[[4]];sr=b[[2]];st=b[[4]];
L1 = Lpart[[1]];L2=Lpart[[2]];
fatt = KroneckerDelta[si,sr]*KroneckerDelta[sj,st];
If[fatt== 0,elem1=0;,elem1=fatt*
Sum[c[k,L1,mi,L1,mr]c[k,L2,mt,L2,mj]funzRad[k,p$p[[1]],p$p[[2]],p$p[[1]],p$p[[2]]],
{k,0,2 L1,2}];];
fatt = KroneckerDelta[si,st]*KroneckerDelta[sj,sr];
If[fatt== 0,elem2=0;,
elem2=fatt*
Sum[c[k,L1,mi,L2,mt]*c[k,L1,mr,L2,mj]*
funzRad[k,p$p[[1]],p$p[[2]],p$p[[2]],p$p[[1]]],
{k,L2-L1, L1+L2,2}];];
elem1-elem2];

singmatOpJp[J_,i_,j_]:= If[i== j+1,(Sqrt[J*(J+1) - j(j+1)]),0];
singmatOpJm[J_,i_,j_]:= If[i== j-1,(Sqrt[J*(J+1) - j(j-1)]),0];
singmatOpLS[L_,n1_,n2_]:= (1/2)singmatOpJp[L,n1[[1]],n2[[1]]]*singmatOpJm[(1/2),n1[[2]],n2[[2]]]+
 (1/2)singmatOpJm[L,n1[[1]],n2[[1]]]*singmatOpJp[(1/2),n1[[2]],n2[[2]]];

matOpLS[0,x_]:= 0; matOpLS[x_,0]:= 0;
matOpLS[v[k1_],v[k2_]]:=Module[{diff,diffnum,ris},
diff = k1 - k2;
diff = Partition[diff,2];
diffnum = Length[Select[diff,#≠ {0,0}&]];
Which[
diffnum > 1, ris = 0;,
diffnum == 0, ris = Sum[ξ[p$p[[j]]]k1[[2j-1]]*k1[[2j]],{j,1,Length[k1]/2}];,
diffnum == 1,ris = Sum[ξ[p$p[[j]]]singmatOpLS[Lpart[[j]],{k1[[2j-1]],k1[[2j]]},{k2[[2j-1]],k2[[2j]]}],
{j,1,Length[k1]/2}];
];
Return[ris];
];
matOpLS[x_. v[a_] + y_.,z_. v[b_] + t_.]:= x  * z matOpLS[v[a],v[b]]+
z*matOpLS[y,v[b]] +x matOpLS[v[a],t] + matOpLS[y,t];

prodScalare[x_. v[a_] + y_.,z_. v[b_] + t_.]:= x  * z prodScalare[v[a],v[b]]+
z*prodScalare[y,v[b]] +x prodScalare[v[a],t] + prodScalare[y,t];
prodScalare[0,x_]:= 0; prodScalare[x_,0]:= 0;
prodScalare[v[a_],v[b_]]:= If[a== b,1,0];

matOpB[0,x_]:= 0; matOpB[x_,0]:= 0;
matOpB[v[k1_],v[k2_]]:=Module[{ris},
If[k1 ≠ k2, ris = 0;,
ris = Sum[k1[[i]],{i,1,Length[k1],2}]+2*
Sum[k1[[i]],{i,2,Length[k1],2}];];
Bmagn*ris];
matOpB[x_. v[a_] + y_.,z_. v[b_] + t_.]:= x  * z matOpB[v[a],v[b]]+
z*matOpB[y,v[b]] +x matOpB[v[a],t] + matOpB[y,t];

Example 1 : Two p - p equivalent electrons

Quantum numbers and states

Let us consider the case of two equivalent p electrons. From Pauli we know tat in this case only terms 3P((,)^1)S((,)^1)D are allowed.

We can fix global variables for the two particles by fixParticles (the principal quantum number is necessary only for two non equivalent electrons):

fixParticles["p","p"];

All terms (i.e. LS levels eigenspaces) contains the states Lz=0  and Sz=0 (as the total number of electrons is even). Here we list them using eigVLS

Do[Print[{L,S,eigVLS[L,0,S,0]}];,{L,0,2},{S,0,1}]

We can check that the states are in effect eigenstates of (L^2) and (S^2):

{Sum[Abs[Simplify[opL2[eigVLS[L,0,S,0]] - L(L+1)eigVLS[L,0,S,0]]],{L,0,2},{S,0,1}],
Sum[Abs[Simplify[opS2[eigVLS[L,0,S,0]] - S(S+1)eigVLS[L,0,S,0]]],{L,0,2},{S,0,1}]}

The non zero vectors are directly given by the function allowedStates[].

allowedStates[]

Electrostatic interactions

namesMultiplets = {"S","P","D","F","G","H"};

These are the states in the subspace (Lz=0 , Sz=0):

vec = allowedStates[];

Using the operator opEE we can immediately write the matrix for electrostatic interaction on these states:

matHI = Table[0,{Length[vec]},{Length[vec]}];
Do[
matHI[[i,j]]=Simplify[opEE[vec[[i,3]],vec[[j,3]]]];
matHI[[j,i]]=matHI[[i,j]];,{i,1,Length[vec]},{j,1,i}]

MatrixForm[matHI]

The reader can check that this matrix is identical to the one in notebooks [NB3, NB4].

LS interaction

To write the matrix elements of LS interaction let us first generate the possible vectors (remember that the format is {J, L, S, vett}:

vec = allowedJStates[];

The interaction matrix is obtained directly from matOpLS

matLS = Table[0,{Length[vec]},{Length[vec]}];
Do[
matLS[[i,j]]=Simplify[matOpLS[vec[[i,4]],vec[[j,4]]]];
matLS[[j,i]]=matLS[[i,j]];
,{i,1,Length[vec]},{j,1,i}]

Let us compute the multiplicity of the terms, (2J+1), and order all terms in decreasing order in J

vmolt=Table[2*vec[[i,1]]+1,{i,1,Length[vec]}];
indexJ = Ordering[vmolt,Length[vmolt],Greater];

{vec,matLS,vmolt} = {vec[[indexJ]],matLS[[indexJ,indexJ]],vmolt[[indexJ]]};

We collect name for J eigenstates

string = Table[""^ToString[2 vec[[i,3]]+1](Subscript[namesMultiplets[[vec[[i,2]]+1]], vec[[i,1]]]),{i,1,Length[vec]}];

and write the formal form of the LS matrix interaction

TableForm[matLS,TableHeadings->{string,string},TableSpacing-> {3,3}]

As expected the LS interaction do not commute with L and S but is diagonal in J.

We can compute on these states the electrostatic energy with opEE

electrostaticEnergy=Simplify[Table[opEE[vec[[i,4]],vec[[i,4]]],{i,1,Length[vec]}]];

As usual is simpler to consider the energies by subtracting the (weighted) mean

meanEnergy = Simplify[(vmolt.(Tr[matLS,List] + electrostaticEnergy))/Tr[vmolt]]

matH = matLS+DiagonalMatrix[electrostaticEnergy - meanEnergy];

The eigenvalues are found with the function Eigenvalues.

eig=Eigenvalues[matH]

To have an idea of their behavior let us consider them as function of x= ξ/F (i.e. LS interaction over electrostatic interaction). This can be done by the following substitution:

eigx = eig/.{(G^("1"))["s","p"]-> 1,(F^("2"))["p","p"]-> 1,ξ["p"]-> x};

This is the behavior as a function  of x. On the lhs of the graph we have a pure Russel Sunders coupling, while for large x the LS interaction prevails and the coupling is jj.

Show[Table[Plot[(eigx[[i]]/(Sqrt[eigx.eigx]))/.x->t/(1-t) ,{t,0,1}],{i,1,Length[eigx]}],PlotRange-> {{0,1},{-1,1}}]

Magnetic interaction

The Hamiltonian matrix on eigenstates of J. Jz is immediately written with the operator matOpB. We select with allowedJStates the eigenstates with Jz=1

vec=allowedJStates[1];

and compute the matrix

matB = Table[0,{Length[vec]},{Length[vec]}];
Do[
matB[[i,j]]=Simplify[matOpB[vec[[i,4]],vec[[j,4]]]];
matB[[j,i]]=matB[[i,j]];
,{i,1,Length[vec]},{j,1,i}]

As in previous example we order the states in Jz and print the matrix

vmolt=Table[2*vec[[i,1]]+1,{i,1,Length[vec]}];
indexJ = Ordering[vmolt,Length[vmolt],Greater];

{vec,matB,vmolt} = {vec[[indexJ]],matB[[indexJ,indexJ]],vmolt[[indexJ]]};

TableForm[matB,TableHeadings->{string,string},TableSpacing-> {3,3}]

It is easily seen that the diagonal elements are the one expected form Landé gyromagnetic factor.

The position of vectors with non zero J is

nonzero = Drop[Range[Length[vec]],Flatten[Position[vec[[All,1]],0]]];

{valJ,valL,valS}=vec[[nonzero,#1]]&/@{1,2,3};

gFactorLande = 1 +( valJ(valJ+1)-valL(valL+1) + valS(valS+1))/(2valJ(valJ+1))

Exercise

The reader can repeat the computation for other pairs, such as for non equivalent p-p electrons, or s-p pairs etc.