Lithium

Quantum Mechanics
K. Konishi, G. Paffuti

Use the Mathematica function NVariationalBound to give an estimate of He and Li, S state energies.

ClearAll["Global`*"];
Off[General::spell, General::spell1,InterpolatingFunction::"dmval"];

The problem

General remarks

In the main text and in exercises we gave a simplified variational solution of helium ground state and excited s-states. The main hypothesis was a factorized form of wave functions, corrected possibly by antisymmetrization. In this notebook we give a mild generalization of the procedure, covering atoms constructed with s-states electrons in the ground state or in excited states for some particular cases. The limitations are due to the exclusion of higher angular momentum.

Let us call ψ[r,σ] a single particle wave function, r is the radial coordinate and σ its spin coordinate. A simple n-particles state built with such a kind of functions can be taken of the form

Ψ = (1/(Sqrt[n!]))𝒜 ((Subscript[ψ, 1])[(Subscript[r, 1]),(Subscript[σ, 1])]... (Subscript[ψ, n])[(Subscript[r, n]),(Subscript[σ, n])])

𝒜 stands for the antisymmetric sum on permutations. We will try to use these kind of states to give a simple estimate of energies for alkaline atoms.

The state () is just a determinant, properly normalized if, as we assume below, the functions (Subscript[ψ, i])are mutually orthogonal. These states are not general enough to describe general eigenstates of angular momentum and spin, as we know they must be combined in linear superpositions with Clebsch-Gordan coefficients to be eigenstates of L and S. One exception is the following: suppose that the ψ are single particle eigenstates of some effective one particle potential, we can "build" atoms by filling the single particle spectrum, with one particle per state to satisfy Pauli principle. This is in effect what we have done speaking of helium 1(s^2)state or 1s2s state. If the potential is a central potential its levels are grouped in degenerate shells, s-shells, p-shells etc., following angular momentum. We will consider only s-shell, thus at most atoms with 4 electrons, after the 2p state must be included. If the shell are completely filled, or only one electron is outside a filled shell, a single determinant is sufficient to describe the system, as there is not other determinants which can be constructed with these quantum numbers for single electrons.  Let us make some example.

Let us call α and β spinors describing respectively up and down eigenstates of (Subscript[s, z]). For ground state of the helium the only possible assignment of quantum number is described by the (Slater) determinant

det  (({{(Subscript[φ, 1s])[(Subscript[r, 1])](Subscript[α, 1]), (Subscript[φ, 1s])[(Subscript[r, 2])](Subscript[α, 2])}, {(Subscript[φ, 1s])[(Subscript[r, 1])](Subscript[β, 1]), (Subscript[φ, 1s])[(Subscript[r, 2])](Subscript[β, 2])}})) = (Subscript[φ, 1s])[(Subscript[r, 1])](Subscript[φ, 1s])[(Subscript[r, 2])]((Subscript[α, 1])(Subscript[β, 2]) - (Subscript[β, 1]) (Subscript[α, 2]))

This is not true for a generic excited 1s2s state. While (Subscript[S, z])= ± 1 have only one representative, like

det  (({{(Subscript[φ, 1s])[(Subscript[r, 1])](Subscript[α, 1]), (Subscript[φ, 1s])[(Subscript[r, 2])](Subscript[α, 2])}, {(Subscript[φ, 2s])[(Subscript[r, 1])](Subscript[α, 1]), (Subscript[φ, 2s])[(Subscript[r, 2])](Subscript[α, 2])}}))= (Subscript[α, 1])(Subscript[α, 2])((Subscript[φ, 1s])[(Subscript[r, 1])](Subscript[φ, 2s])[(Subscript[r, 2])]-(Subscript[φ, 2s])[(Subscript[r, 1])](Subscript[φ, 1s])[(Subscript[r, 2])])

The states with (Subscript[S, z])= 0 admit two possibilities

det  (({{(Subscript[φ, 1s])[(Subscript[r, 1])](Subscript[α, 1]), (Subscript[φ, 1s])[(Subscript[r, 2])](Subscript[α, 2])}, {(Subscript[φ, 2s])[(Subscript[r, 1])](Subscript[β, 1]), (Subscript[φ, 2s])[(Subscript[r, 2])](Subscript[β, 2])}}));    det  (({{(Subscript[φ, 1s])[(Subscript[r, 1])](Subscript[β, 1]), (Subscript[φ, 1s])[(Subscript[r, 2])](Subscript[β, 2])}, {(Subscript[φ, 2s])[(Subscript[r, 1])](Subscript[α, 1]), (Subscript[φ, 2s])[(Subscript[r, 2])](Subscript[α, 2])}})).

It is well known that the spin symmetric combination is the 0-component of the triplet S=1 while the antisymmetric spin combination is the singlet S=0.

We will consider, mainly,  only cases where only one Slater determinant is needed. This applies to Lithium states as 1(s^2)ns, or triplet Beryllium states etc.

Suppose that we have at our disposal a set of functions ψ which depend on some parameters. We can compute mean energy on these states and by minimizing with respect to the parameters obtain a variational upper bound on actual energy levels. Let us note that this is rigorous upper bound and must be satisfied by true energy levels, neglecting spin interactions and in general relativistic effects which we do not consider here.

The computation of mean energy can be pursued by using general formulas for operators between Slater determinants, as explained in the book (see chapter [*] ), but let us retrieve from scratch the formulas we need, in this way maybe the reader can grasp some intuition on general formulas.

We have to compute

ℳ =(1/(Sqrt[n!]))𝒜 ((Subscript[ψ, 1])[(Subscript[r, 1]),(Subscript[σ, 1])]... (Subscript[ψ, n])[(Subscript[r, n]),(Subscript[σ, n])]). H . (1/(Sqrt[n!]))𝒜 ((Subscript[ψ, 1])[(Subscript[r, 1]),(Subscript[σ, 1])]... (Subscript[ψ, n])[(Subscript[r, n]),(Subscript[σ, n])])

H is a symmetric operator, then all matrix elements for each component of the left hand side vector are the same, they are just a permutation of indices which is reabsorbed by changing the index names in the right hand side sum. A more formal way to say the same thing is the following. H is symmetric, thus it commutes with every permutation and in particular with the antisymmetrizer 𝒜. We can move 𝒜 on the right hand side and using the obvious property (𝒜^2) = n! 𝒜 we have

ℳ = (Subscript[ψ, 1])[(Subscript[r, 1]),(Subscript[σ, 1])]... (Subscript[ψ, n])[(Subscript[r, n]),(Subscript[σ, n])] . H .𝒜 ((Subscript[ψ, 1])[(Subscript[r, 1]),(Subscript[σ, 1])]... (Subscript[ψ, n])[(Subscript[r, n]),(Subscript[σ, n])])

which is what we said above. Our states have spin dependence factorized and as our H will not depend on spin, only the terms in the right hand side with the same spin part as the left hand side will survive in matrix elements. Let us make an example. If f and g are the wave functions for 1s and 2s states, a Lithium ground state function can be written as (the subindices are a short notation for variable, (Subscript[f, 1])means f[(Subscript[r, 1])] etc. )

(slater1={Table[(Subscript[f, i])(Subscript[α, i]),{i,1,3}],Table[(Subscript[f, i])(Subscript[β, i]),{i,1,3}],Table[(Subscript[g, i])(Subscript[α, i]),{i,1,3}]})//MatrixForm

The determinant contains many terms, but if the (non antisymmetrized) state is ( (Subscript[f, 1]) (Subscript[α, 1]) (Subscript[f, 2]) (Subscript[β, 2]) (Subscript[g, 3]) (Subscript[α, 3])) only those term with spin component  (Subscript[α, 1]) (Subscript[β, 2]) (Subscript[α, 3]) give a contribution to matrix elements, i.e.

Coefficient[Det[slater1],(Subscript[α, 1]) (Subscript[β, 2]) (Subscript[α, 3])]

In the following section we will write a few functions which extract the relevant orbital part from a given Slater determinant. The following task is to write matrix elements for H. This is simple for factorized states, it is sufficient to make integrals, and as general states are linear combination of these factorized wave functions our task will be relatively simple. Finally we will make some explicit model for our single particle wave functions and obtain bounds for energies.

Wave functions and interactions

We will denote a non antisymmetrized state by a list of doublets, the first component for spin the second for orbital part. 1 will mean first function, 2 the second function and so on. For the spin we can use the symbols we want, as this variable will cancel out, for simplicity we will use 0, 1, to denote spin down and up. For example, the lithium state above will be written as

stateLi = {{0,1},{1,1},{0,2}};

The third electron has the wave function number 2, the first two electrons wave function number 1.

The function makeWaveFunction below enumerates all permutations of the initial states and compares their spin content (the first components of the lists) with original spin components, if they agree the orbital part is added to the original state, with the relative sign of the permutation. We assign a tag to the products of wave functions, 𝓀, for ket. κ[1,1,2] will means the product of first function in the variables 1 and 2 and of second function for the variable number 3.

makeWaveFunction[state_]:=
Module[{spin0,ψ0,ψ,spin,sign1,vec,allstates},
spin0 = state[[All,1]];
ψ0 = state[[All,2]];
sign1=Signature[state];
allstates = Permutations[state];
ψ = 𝓀@@ψ0;
Do[vec = allstates[[k]];
spin = vec[[All,1]];
If[spin == spin0,
ψ=ψ+sign1*Signature[vec]*𝓀@@vec[[All,2]];];,
{k,2,Length[allstates]}];
ψ]

As an example

makeWaveFunction[stateLi]

Which coincides with the above result : (Subscript[f, 1]) (Subscript[f, 2]) (Subscript[g, 3]) - (Subscript[g, 1])(Subscript[f, 2]) (Subscript[f, 3]) .

Now we have to write matrix elements. The Hamiltonian is made of two kind of operators: one particle operators, like kinetic energy and interaction with nuclear charge:

(Subscript[H, 1])= (∑_(i=1)^n)(  (1/2) (Subscript[p, i]^2)  - (Z/(Subscript[r, i])))

and two particle operators, i.e. electrostatic interactions between electrons (we use atomic units)

(Subscript[H, 2]) = (∑_(i<j))(1/(|(Subscript[x, i])-(Subscript[x, j])|))

The wave functions of variables not included in integrals of terms in (Subscript[H, 1])or (Subscript[H, 2])must coincide between bra and ket states, as we assume our wave functions are orthogonal. Then to compute matrix elements we have just to do the following (for factorized wave functions).

For (Subscript[H, 1])

Pick a variable, drop this variable form the ket 𝓀, if the rest is equal between bra an ket, compute the integral and add the result to the energy, otherwise skip the step.

For (Subscript[H, 2])

Generate all possible couples. Pick a couple, drop this couple form the ket 𝓀, if the rest is equal between bra an ket, compute the integral and add the result to the energy, otherwise skip the step.

We do not really compute integrals, just give a symbol to them, we will call 𝒯 the single particle integrals and J the two particle integrals. As an instance J[1,1,1,2] will mean

∫f[(Subscript[x, 1])]f[(Subscript[x, 2])] (1/(|(Subscript[x, 1])-(Subscript[x, 2])|))f[(Subscript[x, 1])]g[(Subscript[x, 2])]

Here are the functions :

OneBody[a_𝓀,-b_𝓀]:= -OneBody[a,b];

OneBody[a_𝓀,b_𝓀]:= Module[{rest1,rest2,couples,ris},
ris = 0;
Do[
rest1 =Drop[a,{i1}];
rest2 = Drop[b,{i1}];
If[rest1== rest2,ris = ris + 𝒯[a[[i1]],a[[i1]]];];,
{i1,1,Length[a]}];
ris]

TwoBody[a_𝓀,-b_𝓀]:= -TwoBody[a,b];

TwoBody[a_𝓀,b_𝓀]:= Module[{rest1,rest2,couples,ris},
couples= Flatten[Table[{i,j},{i,1,Length[a]-1},{j,i+1,Length[b]}],1];
ris = 0;
Do[
{i1,i2} = couples[[tc]];
rest1 = Drop[Drop[a,{i1}],{i2-1}];
rest2 = Drop[Drop[b,{i1}],{i2-1}];
If[rest1== rest2,ris = ris + J[a[[i1]],a[[i2]],b[[i1]],b[[i2]]];];,
{tc,1,Length[couples]}];
ris]

The states will be in general  superposition of factorized states, but the whole result is simply obtained with the distributive property. As an example for our lithium state

stateLi = {{0,1},{1,1},{0,2}};

ψ = makeWaveFunction[stateLi]

The bra vector is represented by

ψ0=𝓀@@stateLi[[All,2]]

and we can write, using as an example infix notation:

singleParticle=Distribute[ψ0~OneBody~ψ]

twoParticles=Distribute[ψ0~TwoBody~ψ]

Now we have to assign an actual procedure to compute matrix elements. We assume to work with reduced radial wave functions, in this case

f𝒯[x_ a_𝒯]:= x f𝒯[a];
f𝒯[a_𝒯]:= Integrate[(1/2)((((∂_r)f[a[[1]],r]))^2) - (Z/r) ((f[a[[1]],r])^2),{r,0,∞},GenerateConditions-> False];

fJ[x_ a_J]:= x fJ[a];
fJ[a_J]:= (Integrate[Integrate[f[a[[1]],r1]f[a[[3]],r1](1/r1)
f[a[[2]],r2]f[a[[4]],r2],{r2,0,r1},GenerateConditions-> False]
,{r1,0,∞},GenerateConditions-> False]+
Integrate[Integrate[f[a[[1]],r1]f[a[[3]],r1](1/r2)
f[a[[2]],r2]f[a[[4]],r2],{r1,0,r2},GenerateConditions-> False]
,{r2,0,∞},GenerateConditions-> False]);

We assign the properties

f𝒯[x_ a_𝒯] := x f𝒯[a];    fJ[x_ a_J] := x fJ[a];

to express linearity with respect multiplication for a number.

The form of the integral for two particles interactions come from multipole expansion, see text and notebook [*],

(1/(|(Subscript[r, 1])-(Subscript[r, 2])|))= (∑_(ℓ=0)^∞)(1/(Subscript[r, >]))(((((Subscript[r, <])/(Subscript[r, >]))))^ℓ)(Subscript[P, ℓ])[Cos(γ)];

with  (Subscript[P, ℓ]) Legendre polynomials and

(Subscript[r, >])=Max[(Subscript[r, 1]),(Subscript[r, 2])];  (Subscript[r, <])=Min[(Subscript[r, 1]),(Subscript[r, 2])];   (Subscript[r, 1])·(Subscript[r, 2]) = (Subscript[r, 1])(Subscript[r, 2])Cos[γ].

For s - states and for reduced wave functions the formula () imply

J[i,j,k,l]= ∫(Subscript[dr, 1])(Subscript[dr, 2]) (Subscript[f, i])[(Subscript[r, 1])](Subscript[f, j])[(Subscript[r, 2])] (1/(Subscript[r, >]))(Subscript[f, k])[(Subscript[r, 1])](Subscript[f, l])[(Subscript[r, 2])]

as the reader can easily check, and this is what is written in the definition of fJ.

Radial wave functions

To perform actually the integrals we have to choose a set of wave functions. We have used orthogonality, then we have to choose an orthogonal set. For 1s, 2s, 3s reduced wave functions we can take (up to a normalization factor)

f[1,r] = r Exp[-Z1 r];
f[2,r] = r Exp[-Z2 r](1 + b r);
f[3,r] = r Exp[-Z3 r](1 + c1 r + c2 (r^2))

The coefficients b, c1, c2 must be determined by imposing orthogonality

b1=xvar/.First[Solve[Integrate[(r^2)Exp[- Z1 r] Exp[-Z2 r](1+ xvar r),{r,0,∞},GenerateConditions-> False]== 0,{xvar}]];

{c1,c2}= {xvar,yvar}/.First[Solve[{Integrate[(r^2)Exp[- Z2 r] (1+b1 r)Exp[-Z3 r](1+ xvar r + yvar (r^2)),{r,0,∞},GenerateConditions-> False]== 0,
Integrate[(r^2)Exp[- Z1 r] Exp[-Z3 r](1+ xvar r + yvar (r^2)),{r,0,∞},GenerateConditions-> False]== 0},{xvar,yvar}]];

Now we normalize the f' s

norm1= Integrate[(r^2) Exp[- Z1 r]^2,{r,0,∞},GenerateConditions-> False];

norm2= Integrate[(((r Exp[-Z2 r](1+b1 r)))^2),{r,0,∞},GenerateConditions-> False];

norm3= Integrate[(((r Exp[-Z3 r](1+c1 r + c2 (r^2))))^2),{r,0,∞},GenerateConditions-> False];

f[1,r_]:= (r/(Sqrt[norm1]))Exp[-Z1 r];
f[2,r_]:= (r/(Sqrt[norm2]))Exp[-Z2 r](1+b1 r);
f[3,r_]:= (r/(Sqrt[norm2]))Exp[-Z3 r](1+c1 r + c2 (r^2));

We have all the elements to perform some computation.

Helium

Let us first check the functions on the known case of Helium. This is the ground state:

groundStateHe = {{0,1},{1,1}};

Bra and ket will coincide but let us define them independently:

ψ0=𝓀@@groundStateHe[[All,2]];
ψ = makeWaveFunction[groundStateHe];

singleParticle=Distribute[ψ0~OneBody~ψ]

twoParticles=Distribute[ψ0~TwoBody~ψ]

We do not have any exchange integral for ground state, as expected.

This is the one-particle contribution to energy :

ene1=Distribute[f𝒯[singleParticle]]

and this is the two-particle contribution

ene2=Distribute[fJ[twoParticles]]

The reader will recognize the same expressions given in the text.

The variational estimation for ground state is obtained by minimizing with respect our free parameter Z1:

Simplify[Minimize[ene1+ ene2,{Z1}]]

This is the well known result.

Lithium

Let us consider now the lithium ground state.

groundstateLi = {{0,1},{1,1},{0,2}};

We can proceed as before :

ψ0=𝓀@@groundstateLi[[All,2]];
ψ = makeWaveFunction[groundstateLi];

singleParticle=Distribute[ψ0~OneBody~ψ]

twoParticles=Distribute[ψ0~TwoBody~ψ]

The last term is an exchange integral between the 1s electron and the 2s electron. Energies now are more complicated

ene1=Distribute[f𝒯[singleParticle]]

and this is the two particle contribution

ene2=Distribute[fJ[twoParticles]]

(ene=(ene1+ene2)/.Z-> 3);

min2sLi=FindMinimum[ene,{Z1,Z2}]

A more accurate Hartree Fock calculation, to be done in the chapter on atomic levels, will give for the energy, in atomic units

E = -7.4327008 a.u.

in agreement with the experimental data. The present result is in agreement with an error of about 2/700.  From this result and form previous calculation on helium we can obtain ionization potential, as lithium without an electron is a two electrons atom with a nuclear charge of 3,  J =  E[He, Z=3] - E[Li] :

ionizationLithium =- min2sLi[[1]] - (((Z-5/16))^2)/.Z-> 3

The experimental result, using 1 a.u. = 2 Ry is :

Ry =1/2;
experimental= 0.39628407Ry

We see that even with our very elementary model we have been able to reproduce ionization potential to few percent level.

Excited state

Let us try a more difficult case, the excited 3 s state.

excitedstateLi = {{0,1},{1,1},{0,3}};

We can proceed as before :

ψ0=𝓀@@excitedstateLi[[All,2]];
ψ = makeWaveFunction[excitedstateLi];

singleParticle=Distribute[ψ0~OneBody~ψ]

twoParticles=Distribute[ψ0~TwoBody~ψ]

The last term is an exchange integral between the 1s electron and the 2s electron. Energies now are more complicated

ene1=Distribute[f𝒯[singleParticle ]];

and this is the two particle contribution

ene2=Distribute[fJ[twoParticles]];

We substitute our values for Z1 and Z2 and find :

We substitute values for Z and for Z2, as previously computed :

(ene=(ene1+ene2)/.{Z-> 3,min2sLi[[2,2]]});

min3s=FindMinimum[{ene, Z3>0 && Z1>0 },{Z1,Z3},AccuracyGoal-> 4]

The Hartree Fock computation gives

For the energy gap between 2 s and 3 s states we estimate

min3s[[1]] - min2sLi[[1]]

The experimental value is

0.2479204 Ry

We see that for excited state the order of magnitude for energy is correct but the value is inaccurate.