Two-electron states 2

Quantum Mechanics
                     K. Konishi, G. Paffuti

Study the problem of energy splitting in atoms with a configuration s-p (or s-d). Deduce radial integrals and LS coupling from the experimental data and check the consistency of the procedure.

Note: The first section of this notebook is identical to the first section of notebook [NB5]. The reader is referred to that notebook for a first introduction to the functions used here.

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 will 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 fixes the sign of the determinant. The quantum numbers are understood completely antisymmetric, in agreement with Pauli's 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 determine 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, then this construction obeys the Pauli principle. In the function coded below the Antisymmetry is assured 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 gives a list of possible states for Lz=0=Sz. These states are representative of the possible multipltes. Each element contains {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-particle operators (LS interaction, interaction 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-particle 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 couple, (i,j) = (r,t) (apart the order). Two kind of functions can appear. If i=r (then necessarily j = t as the couples 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 interaction

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 interaction

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 interaction

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];

s-p configurations

Quantum numbers and states

We can fix global variables for the two particles by fixParticles:

fixParticles["s","p"];

Electrostatic and LS interaction

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

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,eigV}=Eigensystem[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"]-> η};

indexEne = Ordering[eigx/.η-> 0.2];

{eigx,eig,eigV}={ eigx[[indexEne]],eig[[indexEne]], eigV[[indexEne]]};

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]))/.η->t/(1-t) ,{t,0,1}],{i,1,Length[eigx]}],PlotRange-> All,AxesLabel-> {"η/1+η","ΔE[i]/(Sqrt[((ΔE)^2)])"}]

Experimental Data

s-p configurations

Here we write the terms and the fine structure levels of some excited atoms in an s-p configuration. The energies are written with respect to the ground state and are expressed in 1/cm. For these data and similar ones the reader is referred to NIST, National Institute of Standards and Technology. On  line data are at

http : // physics.nist.gov/PhysRefData/ASD/index.html

The energies are written for states

(({)^3)(Subscript[P, 0])((,)^3)(Subscript[P, 1])((,)^3)(Subscript[P, 2])}  ,  (({)^1)(Subscript[P, 1])}

valExp = {
{"Mg",{21850.405,21870.464,21911.178},{35051.264}},
{"Ca",{15157.901,15210.063,15315.943},{23652.304}},
{"Sr",{14317.520,14504.351,14898.563},{21698.482}},
{"Ba",{12266.024,12636.623,13514.745},{18060.261}},
{"Zn",{32311.319,32501.399,32890.327},{46745.404}},
{"Cd",{30113.990,30656.087,31826.952},{43692.384}},
{"Hg",{37645.080,39412.300,44042.977},{54068.781}}};

If the approximation in the preceding is correct (configuration mixing are negligible etc.) these four energies are to be expressed through only 2 parameters, The radial integral G and the ratio η between Ls and electrostatic interaction.

deltavalExp[k_]:=Module[{med,molt,J,ene,deltaE},
ene = Flatten[valExp[[k,{2,3}]]];
J = {0,1,2,1};
molt = (2*J+1);
med = molt.ene/Tr[molt];
deltaE = ene - med]

getEtaFactor[deltaE_,eigIn_]:=Module[{eigNorm,δENorm,eta,etaM,eig,η,factorEM,factorEMMean},
eig = eigIn/.{(F^(x_))[a_,b_]-> 1,(G^(x_))[a_,b_]-> 1,ξ[c_]-> η};
eigNorm = eig/(Sqrt[eig.eig]);
δENorm = deltaE/(Sqrt[deltaE.deltaE]);
eta=Table[η/.FindRoot[eigNorm[[k]]== δENorm[[k]],{η,0,1}],{k,1,Length[eig]}];
etaM = Mean[eta];
eig=eig/.η-> etaM;
factorEM = deltaE/eig;
factorEMMean = Mean[factorEM];
{factorEM,factorEMMean,eta,etaM}
];

Here we make a loop on all atomic data and collect in the vector meanValues the computed quantities

meanValues = Table[getEtaFactor[deltavalExp[i],eig],{i,1,Length[valExp]}];

Let us order the names of the states in increasing energy order

string

stringaOrd = string[[{4,3,1,2}]]

These are predictions for G

tab1= Flatten/@meanValues[[All,{1,2}]];
TableForm[tab1,TableHeadings->{valExp[[All,1]],Flatten[{{stringaOrd,"G-Mean"}}]} ]

and these are predictions for η

tab1=Flatten/@meanValues[[All,{3,4}]];
TableForm[tab1,TableHeadings->{valExp[[All,1]],Flatten[{{stringaOrd,"η-Mean"}}]} ]

We see a good consistency for most elements but Hg, where mixing configuration are important.

Landè factors

The experimental gyromagnetic factor is extracted from the linear slope in B from energy levels, at least in principle. In perturbation theory this amount to compute

⟨ψ|μ∑ (Lz + 2 Sz) B|ψ⟩ ≡ gJ Jz (μ B)

on  eigenstates of H and Jz. In our computation it is sufficient to compute mean values of the sum of Lz+ 2 Sz for each electron.

If the LS coupling is relevant the eigenstates are no more eigenstates of L and S, so a deviation from Landè factors is expected, of order δELs /δE, where δE is the splitting between the terms (i.e. LS levels) in the Hamiltonian.

We computed in the previous section eigenstates of H with Jz = 0. To get eigenstates with Jz=1 and compute directly gJ it is sufficient to apply the rising operator Jp. We first compute the expression for vectors with Jz = 1 (basis vectors)

statesJz1 = normalize/@(Expand/@(opJp/@vec[[All,4]]));

The coefficients of these basis vectors to form an eigenstate of the Hamiltonian are those computed in the previous section, and are stored in the vector eigV. To get numerical values we substitute the values of G and η just obtained, stored in the matrix meanValues. In atomic units the coefficient entering in the LS interaction (ξ) is given by ξ = G η.

gJ=Table[0,{Length[valExp]}];
Do[
valueOgG =  meanValues[[i,2]];
xival =  meanValues[[i,2]]* meanValues[[i,4]];
eigVnum = eigV/.{(F^(x_))[a_,b_]-> 0,(G^(x_))[a_,b_]->valueOgG,ξ[c_]-> xival};
eigVnum = eigVnum/(Norm/@eigVnum);
gJ[[i]]=Table[Simplify[matOpB[Expand[eigVnum[[i]].statesJz1],Expand[eigVnum[[i]].statesJz1]]],{i,1,4}]/.Bmagn-> 1;
,
{i,1,Length[valExp]}];

TableForm[N[gJ],TableHeadings->{valExp[[All,1]],stringaOrd},TableSpacing-> {3,4} ]

Some experimental values, for Ba and Zn are given below :

gBa = {0, 1.45, 1.52, 1.02};  gZn = {0, 1.496, 1.505, unknown};