Atomic Model

A simple model for the structure of few-electron atoms is presented. A model potential for the helium atom is introduced (it represents a good approximation to the so-called Hartee-Fock potential for the helium atom). The eigenenergies and eigenfunctions of this potential are determined by

a) the variational method for the lowest-lying eigenstate for a given angular momentum symmetry;

b) the numerical solution of the radial Schroedinger equation;

c) a matrix diagonalization using a Slater-type orbital basis.

A simple model potential for the ground state of helium can be obtained as follows: assume that in the ground state the wavefunction corresponds to an 1s^2 configuration (one spin-up, one spin-down electron). What potential should the electron experience?

The nuclear attraction -2/r and a repulsive potential due to the presence of the other electron.

What is the repulsive potential due to the other electron (which has the same 1s-wavefunction, it has just opposite spin projection)?

We have an chicken-and-egg problem here: given a potential, we can determine the |1s> state, and calculate the repulsive part of the potential correctly. This is the objective of a self-consistent field calculation. Then we would need to repeat the calculation of the new |1s> state using the new potential. Repeating this procedure one can come up with a |1s> state and corresponding potential that provide the lowest energy.

We break the closed circuit by stating an approximate answer for the potential in helium:

> V:=r->-2/r+1/r*(1-exp(-3.36*r)*(1+1.665*r));

This is the potential experienced by an 1s-electron in the helium atom: at short distances the electron just feels the full nuclear attraction (-2/r) as the second expression approaches a constant as r goes to zero. At large distances the electronic repulsion term goes like 1/ r and screens the nucleus by one unit so that the overall potential goes like -1/ r .

We can check for the consistency by a simple variational calculation:

We start with an unnormalized 1s state that depends on a 'charge' parameter (cf. the hydrogen-like wavefunction)

> chi:=r*exp(-beta*r);

We are using the radial wavefunction (i.e. r*R_nl(r)) so that the radial kinetic energy is just proportional to the second derivative.

> assume(beta>0);

> E1sT:=Int(chi*(-1/2*Diff(chi,r$2)+V(r)*chi),r=0..infinity)/Int(chi^2,r=0..infinity);

> E1sT:=value(E1sT);

> plot(E1sT,beta=1.55..1.8);

> beta0:=fsolve(diff(E1sT,beta),beta=1..2);

> E0:=subs(beta=beta0,E1sT);

Note that the beta-value is less than 2. The value of beta=2 would be obtained if there was only one 1s-electron (the hydrogen-like solution for the ground state of the He+ ion). The fact that the electron wants to have a slightly more diffuse wavefunction reflects the so-called inner screening: the electrostatic repulsion when combined with the simple model of two identical electrons (apart from the spin projection) leads to two electrons which repel each other on average (independent electron model). Both electrons experience the same common central potential and are bound by the same eigenenergy.

Now we would like to check two things:

1) how accurate is this variational solution?

2) how close is the used potential to the potential produced by this 1s-wavefunction?

Part 1 can be answered by using dsolve[numeric]. For part 2 we need to solve the Poisson equation.

Let us start with the first question: To solve the SE we start the integration not at zero, but at some small value of r=eta to avoid the singularity.

> eta:=10^(-8);

> IC:=phi(eta)=eta,D(phi)(eta)=1;

> Et:=-0.9042;

> SE:=-1/2*diff(phi(r),r$2)+(V(r)-Et)*phi(r);

> sol:=dsolve({SE,IC},phi(r),numeric,output=listprocedure):

> phir:=subs(sol,phi(r)):

> plot(phir,0..10);

> E1s:=Et;

The numerically exact solution for the radial SE yields a slightly lower eigenvalue of E_1s = -0.9042 a.u. (1 a.u. = 27.12 eV) compared to the variational result of E_1s_v =-0.897 a.u.. This means that the true eigenfunction is somewhat different from the hydrogenic form.

The numerical eigenfunction is not normalized properly. To compare the graphs we simply change the normalization of the variational answer. In fact, the variational state was normalized such that the derivative of the function equals unity at r =0.

> eval(subs(r=0,diff(subs(beta=beta0,chi),r)));

> with(plots):

Warning, the name changecoords has been redefined

> P1:=plot(phir,0..6,color=red):

> P2:=plot(subs(beta=beta0,chi),r=0..6,color=blue):

> display([P1,P2]);

Note that the agreement for the energy was at the level of 7/900 , i.e. in the 1 % range. The deviation between the wavefunctions is in the 10 % range.

Nevertheless, we can state that the simple hydrogenic wavefunction catches the main feature of the numerical solution, namely the most likely location in r for the 1s electron.

The result for the binding energy for one of the two helium electrons is the first quantity that can be compared with experiment:

The ionization potential of helium is measured to be 24.481 eV (cf. R. Liboff, table 12.2). Our answer is:

> -Et*27.12*_eV;

For any atom the eigenenergy of the highest occupied orbital should equal the negative of the ionization potential. Our result is indeed quite close.

There is another quantity that can be measured, namely the total energy of the atom (equal to the sum of ionization energies for both electrons). This is not simply twice the eigenenergy, since after ionizing on of the two He-electrons, the other is left in a hydrogen-like state (energy of -2 a.u. due to Z=2, and E_1s=-Z^2/2).

Combining this answer with the calculated 1s binding energy we have for the total binding energy:

> (Et-2)*27.12*_eV;

The electrostatic repulsion:

Now we look at the question as to what potential is associated with the approximate wavefunction. For this we need a solution to the Poisson equation for a spherically symmetric charge distribution based on the properly normalized variational wavefunction.

> A1s:=1/sqrt(int(chi^2,r=0..infinity));

> rho:=subs(beta=beta,(A1s*chi)^2);

> int(rho,r=0..infinity);

Note that the 4Pi from the integration over theta and phi are cancelled by the square of the spherical harmonic Y00!

The solution to Poisson's equation in multipole expansion leads to the monopole expression:

> V0:=unapply(simplify(int(rho,r=0..R)/R+int(rho/r,r=R..infinity)),R);

> plot([V(r)+2/r,subs(beta=beta0,V0(r))],r=0..5,color=[red,blue],view=[0..5,0..2]);

We recognize that the electronic repulsion in our helium atom of independent electrons is modelled after the simple hydrogen-like solution to the problem. The above potential shows how asymptotically the charge distribution of electron 1 screens one of the protons for electron 2 (located at a large r -value).

A sophisticated central-field or Hartree-Fock calculation take the electrostatic repulsion due to the numerically obtained charge density and re-calculates the eigenenergy/eigenfunction until convergence is achieved.

Electronic excitations

We proceed to calculate the energy levels for the 2s and 2p states. We simply assume that we can use the potential obtained for the ground state, and calculate the energy spectrum for this potential. For the 2s-state we will not carry out a variational calculation as the energy will not be guranteed to be above the exact eigenenergy for the given potential. We repeat our trial-and-error procedure to find a radial function with one node at r >0:

> IC:=phi(eta)=eta,D(phi)(eta)=1;

> Et:=-0.15768;

> SE:=-1/2*diff(phi(r),r$2)+(V(r)-Et)*phi(r);

> sol:=dsolve({SE,IC},phi(r),numeric,output=listprocedure):

> phir:=subs(sol,phi(r)):

> plot(phir,0..20,-0.5..0.5);

> E2s:=Et;

Now let us carry out the calculation for the 2p state. First we carry out the hydrogen-like variational calculation.

> l:=1;

The trial function is hydrogen-like: the radial function picks up a factor of r^l, and we add the centrifugal potential to the radial Hamiltonian:

> chi:=r^(1+l)*exp(-beta*r);

> E2pT:=Int(chi*(-1/2*Diff(chi,r$2)+(V(r)+1/2*l*(l+1)/r^2)*chi),r=0..infinity)/Int(chi^2,r=0..infinity);

> E2pT:=value(E2pT);

> plot(E2pT,beta=0.45..0.6);

> beta0:=fsolve(diff(E2pT,beta),beta=0..2);

> E2pV:=subs(beta=beta0,E2pT);

Now we verify the variational calculation by a numerical solution:

> IC:=phi(eta)=eta^(l+1),D(phi)(eta)=(l+1)*eta^l;

> Et:=-0.12699;

> SE:=-1/2*diff(phi(r),r$2)+(V(r)+1/2*l*(l+1)/r^2-Et)*phi(r);

> sol:=dsolve({SE,IC},phi(r),numeric,output=listprocedure):

> phir:=subs(sol,phi(r)):

> plot(phir,0..25);

> E2p:=Et;

The answer is just slightly below the value obtained by the variational calculation with a hydrogenic wavefunction.

> (E2p-E1s)*27.12*_eV;

We learn that the model predicts an excitation energy of slightly more than 21 eV. What does this have to do with reality?

Unfortunately the connection is not very direct: inspection of a spectroscopic table for the He atom reveals that it is complicated, and that many levels correspond to the primitive (1s 2p) configuration: this is the result of complications arising from coupling orbital and spin angular momenta for the two electrons.

Nevertheless, if we are generous, we can compare the calculated excitation energy with the average result which in some sense ignores the spin-orbit interactions (these are also called magnetic or fine structure interactions). The model does provide a realistic assessment for the average transition energy (which is in the UV regime). For an understanding of the visible spectrum one cannot ignore the fine structure.

One important piece of information that does emerge from this calculation is the lifting of the ( nl )-degeneracy observed in the hydrogen atom spectrum: the 2s and 2p states are no longer degenerate, and thus, there is the possibility of observing photons that correspond to 2p-2s de-excitations.

Matrix representation

We proceed with the calculation of approximate spectra in an angular momentum symmetry sector. For this purpose we first define ourselves a suitable basis. We defined so-called Slater-type orbitals for a given l -symmetry, and use a Gram-Schmidt procedure to orthogonalize them.

Our radial functions are real-valued, and therefore we have the simple inner product and normalization:

> IP:=(f,g)->int(f*g,r=0..infinity):

> NO:=xi->xi/sqrt(IP(xi,xi));

> STO:=(n,l,beta)->r^(n+l)*exp(-beta*r);

> L:=0;

> N:=8;

> B1:=[seq(STO(n,L,17/10),n=1..N)];

These states are linearly independent, but not orthogonal!

> IP(B1[1],B1[2]);

The Gram-Schmidt procedure takes a list of functions (state vectors) and orthonormalizes them using the two procedures.

> GS:=proc(vecs) local i,n,j,res,xi;

> n:=nops(vecs); res:=[NO(vecs[1])];

> for i from 2 to n do:

> xi:=vecs[i]-add(IP(vecs[i],res[j])*res[j],j=1..i-1);

> res:=[op(res),NO(xi)]; od: end:

> B1ON:=GS(B1):

> IP(B1ON[5],B1ON[5]);

> IP(B1ON[4],B1ON[5]);

> B1ON[2];

> plot([seq(B1ON[i],i=1..4)],r=0..10,color=[red,blue,green,black]);

Now we recycle the code from MatrixRep.mws to generate the Hamiltonian matrix: we define procedures to calculate the kinetic and potential energy matrix elements:

> Tkin:=(phi,psi)->-1/2*int(expand(phi*diff(psi,r$2)),r=0..infinity);

> Vpot:=(phi,psi)->int(expand(phi*psi*(V(r)+L*(L+1)/(2*r^2))),r=0..infinity);

> with(LinearAlgebra): Digits:=15:

> HM:=Matrix(N):

We make use of the symmetry of the hamiltonian matrix: it allows to save almost a factor of 2 in computation time.

> for i from 1 to N do: for j from 1 to i do: HM[i,j]:=Tkin(B1ON[i],B1ON[j])+Vpot(B1ON[i],B1ON[j]); od: od:

> for i from 1 to N do: for j from i+1 to N do: HM[i,j]:=HM[j,i]: od: od:

> HMf:=map(evalf,HM):

> evals:=Eigenvalues(HMf, output='list'):

> ev_s:=sort(map(Re,evals));

We see that the chosen basis set can reproduce the 1s and 2s states, but that it does not have a prediction for a bound 3s eigenstate.

Exercise 1:

Increase the matrix size and observe the stability of the 1s and 2s eigenvalues. Do you find an acceptable 3s eigenvalue?

Exercise 2:

Change the value of the parameter that controls the Slater type orbital (STO) basis from the chosen value of 17/10. Find the best basis for the 3s eigenvalue.

Now the eigenfunctions. We need a sorting procedure to arrange the result from the eigenvector calculation in proper order.

> VE:=Eigenvectors(HMf,output='list'):

> Vp:=[seq([Re(VE[i][1]),VE[i][2],map(Re,VE[i][3])],i=1..nops(VE))]:

> Min:=proc(x,y); if type(x,numeric) and type(y,numeric) then if x<=y then RETURN(true): else RETURN(false): fi; elif type(x,list) and type(y,list) and type(x[1],numeric) and type(y[1],numeric) then if x[1]<=y[1] then RETURN(true): else RETURN(false): fi; elif convert(x,string)<=convert(y,string) then RETURN(true): else RETURN(false): fi: end:

> VEs:=sort(Vp,Min):

Suppose that we would like to see the eigenfunctions corresponding to the four lowest-lying eigenvalues. The eigenvector for a given eigenvalue contains the expansion coefficients for the expansion of the eigenstate in terms of the chosen basis states.

> for i from 1 to 4 do:

> psi0:=add(VEs[i][3][j]*B1ON[j],j=1..N):

> No:=1/sqrt(int(expand(psi0^2),r=0..infinity));

> phi_a[i]:=add(No*VEs[i][3][j]*B1ON[j],j=1..N): od:

> plot([seq(phi_a[i],i=1..4)],r=0..15,color=[red,blue,green,black]);

The basis functions do not have a good span at large r where the higher states wish to reside. A smaller value of beta should improve matters.

For this reason one mixes usually STO's with different values of beta. The Gram-Schmidt process takes care of the orthonormalization.

Exercise 3:

Calculate the eigenfunctions for a few low-lying eigenstates in the L=2 and L=3 symmetry sectors. Observe how the eigenenergies approach the hydrogenic result in this case of E_n=-1/(2n^2). Can you explain this behaviour? Hint: graph the effective potential, and compare with that of a hydrogen atom.

>