k
Using pure state occupation numbers of zero and one in the grand canoniÂ¬
cal density operator for the initial computation of G(E), then occupation
numbers determined from Eq. (2.42) on subsequent computations, a self
consistent set of occupation numbers can be sought.
2.4 Computational Considerations and Applications
The most time consuming step in the construction of the [2,1] Pade1
approximant to the electron propagator is the construction of the moment
matrices. The fourth moment matrix (given in Appendix 2) is particularly
difficult since it involves five unrestricted orbital summations plus
another two symmetry restricted, orbital summations of twoelectron inteÂ¬
grals. Using direct summation techniques, the time needed to construct
this matrix is roughly proportional to N^. This is a formidable computaÂ¬
tional problem, but one that must be accepted in favor of the more manageÂ¬
able matrix dimensions.
Fortunately, the problem is not as intimidating as it might seem
on first appearance. The "bruteforce" summation of twoelectron inteÂ¬
grals in the moment matrices resums certain partial sums which may appear
in more than one term or matrix element. These redundant summations can
be avoided with considerable savings in computer time by computing the
partial sums once and reusing them. Two specific partial sums we have
employed are
[ijkl] = I < s s1k1>
s, s'
(2.43)
44
{i j  kl} = l . (2.44)
s, s1
Since these partial sums contain a double summation which is performed
7 c
only once, the original N problem is effectively reduced to N . The
construction of the moment matrices is now comparable in difficulty to
the transformation of the twoelectron integrals from the primitive basis
to the computational (usually HartreeFock) basis which is also roughly
5
proportional to N .
When the number of twoelectron integrals is too large to be held
in core, their random access from peripheral storage becomes relatively
time consuming. The partial sums are much more efficiently constructed
from ordered lists of twoelectron integrals which can be read into
primary (core) storage when needed. For the partial sums defined above,
the twoelectron integrals must be sorted into ordered lists of the type
and where * indicates all orbital indices which yield a
nonzero integral for the corresponding i,jth distribution.
The integral sorts are performed using the Yoshimine sorting techÂ¬
nique (Yoshimine, 1973). Briefly summarized, this technique involves a
partition of available core into a number of buffers. Each buffer holds
integrals corresponding to a specific i,j distribution, e.g. .
(When the number of distributions is large, several may be held in each
buffer.) Reading through the twoelectron integral list, integrals are
then sorted into the appropriate buffers. As each buffer fills, it is
written to direct access, peripheral storage and assigned a record number.
All record numbers corresponding to integrals from the same buffer are
saved in a "chaining" array for that buffer. After the entire integral
list has been processed and all buffers have been dumped, it is then
45
possible to chain back through the direct access records, copying inteÂ¬
grals of the same distribution back into core. These integrals may then
be further sorted within distributions, e.g. k
finally saved sequentially on a peripheral storage device.
Diatomic nitrogen was the first molecule to be studied with the
[2,1] Fade' approximant. Owing to its abundance in the atmosphere, nitroÂ¬
gen has great chemical interest and has been extensively investigated
both experimentally and theoretically. It is an ideal test case for calÂ¬
culating ionization energies from a correlated, manyelectron formalism
such as propagator theory since both the HartreeFock and AE(SCF) approxÂ¬
imations incorrectly predict the order of the 3o^ and lir^ ionizations
(Cade et_ aj_., 1966). Only when correlation corrections are included is
the correct ordering obtained (Cederbaum and Domcke, 1977 and references
therein).
A double zeta, contracted basis of Gaussian type orbitals (GTO's)
was employed in this calculation. This basis consisted of Huzinaga's
9s,5p set of primitive orbitals (Huzinaga, 1965) which was contracted to
4s,2p (Dunning, 1970). This basis has been optimized by Dunning on the
nitrogen atom and is listed in Table 1. The corresponding one and two
electron integrals were calculated at the experimental internuclear
separation, R=2.068 a.u. (Herzberg, 1955), using the MOLECULE integral
program (Almlof, 1974 ).
The HartreeFock calculation and the twoelectron integral transÂ¬
formation were performed with the program GRNFNC (Purvis, 1973). The
HartreeFock total energy with this basis was E(HF)= 108.8782 a.u.
which is about 3 eV higher than the result of Cade et al. (1966).
There is also a discrepancy in the HartreeFock orbital energies. While
46
Table 1. Contracted Gaussian Basis for Nitrogen.
Nitrogen s orbitals
Contraction
Exponents
Coefficients
5909.4400
0.002001
887.4510
0.015310
204.7490
0.074293
59.8376
0.253364
19.9981
0.600576
2.6860
0.245111
7.1927
1.000000
0.7000
1.000000
0.2133
1.000000
Nitrogen
p orbitals
Contraction
Exponents
Coefficients
26.7860
0.018257
5.9564
0.116407
1.7074
0.390111
0.5314
0.637221
0.1654
1.000000
47
the calculation of Cade Â£t aj_. (incorrectly) predicted the In orbital
energy to be 0.53 eV below the 30^, this calculation predicts the ln^
energy to be 0.05 eV higher. The correct ordering of the 3o^ and Itt^
ionizations with this basis is merely fortuitous, since based on a total
energy criterion, the basis of Cade ejt aJL is more accurate.
The next step of the calculation involved the integral sorts, parÂ¬
tial summations, and the construction of the moment matrices. The poles
and spectral density were finally computed as outlined in the previous
section and are presented along with the [1,0] results in Table 2. The
ionization energies of both approximants seem to be upper bounds to the
experimental results of Siegbahn Â£t al_. (1969), but without exception,
the results of the [2,1] approximant are worse than the [1,0] approximant.
In an attempt to incorporate some ground state correlation into the grand
canonical density operator, new occupation numbers were computed from the
spectral density and the [2,1] approximant was recalculated. This calÂ¬
culation, however, yielded no significant improvements in the ionization
energies.
In order to ascertain whether the poor results from the [2,1] approxÂ¬
imant for nitrogen are representative of other calculations or just the
consequence of a pathological test case, the water molecule was chosen
for a second application. Similarly to the calculation for nitrogen, a
double zeta contracted basis of GTO's was also employed in this calculaÂ¬
tion. Huzinaga's 9s,5p primitive basis for oxygen and 4s primitive basis
for hydrogen were contracted with Dunning's coefficients to 4s,2p and
2s, respectively. The orbital exponents for the hydrogen atoms were
scaled by 1.14 to more realistically represent the effective nuclear
charge in the molecule,and the final basis appears in Table 3. Again,
48
Table 2. Principal Ionization Energies for the Nitrogen
Molecule Resulting from the [1,0] and [2,1]
Propagator Approximants.
Orbital
[l.oi
[2,11
Exp.3
%
427.7
472.4
409.9
2o
g
41.6
46.5
37.3
3og
17.0
30.4
15.5
liru
17.0
23.1
16.8
lau
427.6
478.8
409.9
2 a
u
21.0
30.9
18.6
E(HF) = 108.8782 H.
aSiegbahn et al. (1969).
Table 3.
Contracted Gaussian Basis for Water.
Hydrogen s sets
Exponents
Contraction
Coefficients
13.3615
2.0133
0.4538
0.1233
0.032828
0.231208
0.817238
1.000000
Oxygen s sets
Exponents
Contraction
Coefficients
7816.5400
1175.8200
273.1880
81.1696
27.1836
3.4136
9.5322
0.9398
0.2846
0.002031
0.015436
0.073771
0.247606
0.611832
0.241205
1.000000
1.000000
1.000000
Oxygen p sets
Exponents
Contraction
Coefficients
35.1832
7.9040
2.3051
0.7171
0.2137
0.019580
0.124189
0.394727
0.627375
1.000000
50
the integrals were computed with the MOLECULE program at the equilibrium
internuclear geometry, R(OH) = 1.809 a.u., jHOH = 104.5Â° (Benedict et al.,
1956). A total energy of E(HF)= 76.0082 a.u. was computed with the
llartreeFock portion of GRNFNC and was followed by the twoelectron inteÂ¬
gral transformation. Finally, the integral sorts and partial sums were
performed, the moment matrices constructed, and the poles and spectral
density obtained for the [2,1] approximant. The results for both the
[1,0] and [2,1] approximants are presented in Table 4 and appear to be
upper bounds to the experimental ionization energies. Once more, the
[2,1] results are consistently worse than the [1,0] results. A few
iterations on the occupation numbers yielded no significant improvements.
2.5 Eval_uation of the Moment Conserving Decoupling
Formally, the moment conserving decoupling is an attractive decouÂ¬
pling procedure. Being closely related to the Fade1 approximant method,
this decoupling allows the application of numerous results from the clasÂ¬
sical moment problem to propagator theory. In particular, it was proven
that the sequence of [N,N1] approximants converge uniformly to the
exact electron propagator in a given basis, and it was shown that these
approximants represent a variationally optimum choice of the inner proÂ¬
jection manifold. Why then are the results of the [2,1] approximant so
much worse than the results of the [1,0] approximant? To answer this
question, it is necessary to analyze the three approximations identified
in Section 1.3, namely: basis quality, density operator, and decoupling
procedure.
First of all, since computational economy and not high accuracy was
the criterion for the test calculations on nitrogen and water, polarization
51
Table 4. Principal Ionization Energies for Water Resulting
from the [1,0] and [2,1] Propagator Approximants.
Orbital
[1,01
[2,1]
Exp.a
lai
559.4
619.2
540.2
2al
37.0
44.7
32.2
3a ^
15.4
29.7
14.7
lbl
13.8
32.6
12.6
lb2
19.5
29.6
18.6
E(HF) = 76.0082 H.
aSiegbahn et al_. (1969).
52
functions were intentionally excluded from the basis sets. PolarÂ¬
ization functions are diffuse, virtual orbitals which can be very imporÂ¬
tant in describing electron relaxation and correlation (Purvis and Ohrn,
1974, Cederbaum and Domcke, 1977). It is reasonable to expect that the
addition of polarization functions will improve both the [1,0] and [2,1]
approximants to varying degrees; however, with the same basis and with
the same density operator, the larger inner projection manifold (if
judiciously chosen) should yield a more accurate decoupling. Since this
was not the situation in these test calculations, any improvements in
the basis sets did not seem worthwhile.
Second, it is possible that significant ground state correlation
may have been neglected with our choice of the grand canonical density
operator. With the spin orbital annihilation and creation operators
expanded in the HartreeFock basis and using pure state occupation
numbers of zero or one, this density operator yields the uncorrelated,
HartreeFock ground state average. Rather than explicitly correlating
the density operator (as e.g. through perturbation theory), an attempt
was made to estimate the effect of correlation through the selfconsisÂ¬
tent determination of the occupation numbers as described in Section 2.3.
This procedure was not pursued to true selfconsistency since each iterÂ¬
ation required a complete recalculation of the [2,1] approximant. It
was obvious, however, after the first few iterations that no significant
improvements had been obtained.
Based on the preceding implications, the third approximationthe
inner projection manifold truncationseems to be primarily responsible
for the poor numerical results. Owing to the complicated operator sums
in this manifold, an order analysis (as discussed in Section 1.6) is not
53
readily possible. Consequently, it is extremely difficult to identify
the problem with this decoupling procedure. It can only be concluded
that the number of moments conserved is not a useful criterion for decou
pling. This conclusion is consistent with the uniform convergence of
the [N.Nl] sequence since uniform convergence is not necessarily mono
tonic, but it suggests that more accurate decouplings require the incorÂ¬
poration of more information about the moment expansion than just the
moment matrices. The additional information needed is indeed available
and,in the next chapter, we will demonstrate how it may be extracted
using perturbation theory.
*Babu and Ratner (1972) reported the same conclusion which was based
on an application of their rational approximants to the Hubbard model.
CHAPTER 3
DIAGRAM CONSERVING DECOUPLINGS
3.1 The Diagrammatic Expansion Method
The superoperator formalism which is employed in the previous two
chapters is by no means the only formalism available to formulate decouÂ¬
pling approximations for the electron propagator. Two commonly used,
alternative methods are the functional differentiation method (see e.g.
Csanak et aj_., 1971) and the diagrammatic expansion method (see e.g.
Mattuck, 1967, Fetter and Walecka, 1971, or Cederbaum and Domcke, 1977).
Of these latter two methods, the diagrammatic expansion method has proven
to be particularly effective. This method avoids some of the algebraic
tedium involved in deriving propagator decoupling approximations by
establishing certain rules for constructing and manipulating diagrams
which represent the underlying algebraic structure.
The diagrammatic expansion of the electron propagator is usually
derived using timedependent perturbation theory. The Uelectron HamilÂ¬
tonian is partitioned into an unperturbed part plus a timedependent
perturbation
H = Hq + exp(c111 )V
(3.1)
where e is a small positive quantity. The unperturbed Hamiltonian, Hq,
is chosen to yield an exactly solvable, eigenvalue problem
(3.2)
54
55
and the time dependence of the unperturbed eigenstates is given by
14>0(t)> = exp(iHQt) $Q> . (3.3)
In order to simplify the remaining problem of finding the fully
perturbed eigenstates it is convenient to introduce the "interÂ¬
action representation" (Fetter and Walecka, 1971) by the transformation
Yj(t)> = exp(iHQt)Y(t)> . (3.4)
In this representation, the Schrodinger equation has the form
i Â§f 'fj(t)> = exp(e 11 j) V (t) jfj (t)> (3.5)
where
V(t) = exp(iHQt)Vexp(iHgt) . (3.6)
The time dependence of the interaction eigenstates can be expressed as
Â¡Yj(t)> = Ue(t,t0)Â¥j(t0)> (3.7)
where U (t,tg) is the timeevolution operator. Substituting Eq. (3.7)
into Eq. (3.5), the evolution operator is found to satisfy the differenÂ¬
tial equation
i f UE(t,t0) = exp(et)V(t)UÂ£(t,t0) (3.8)
with the initial condition
uE(t0,t0) = 1 â€¢ (3.9)
It is more convenient to solve for U (t,tg) by first transforming Eq.
(3.8) into an integral equation
56
t
UE(t,t0) = 1  i  dt1 exp(ct)V(t1)Ue(t1,t0) . (3.10)
This integral equation has the form of the Volterra equation of the
second kind (Lowdin, 1967) and is solved iteratively
t
UÂ£(t,t0) = 1  i j dt1 exp(et[)V(t1)
t0
t t^
+ (i)2 j dtj j dt? exp{e(t1+t2)}V(t1)V(t2)Ue(t2>t0) (3.11)
t0 t0
t t j tn_^
= L (i)n j dtj J dt2 . . .  dtn exp{e(t1+t2+ . . . tj)}
n_^ t t t
u0 0 0
x V(tj)V(t2) . . . V(tn) t>tj>t2> . . . > tn . (3.12)
Eq. (3.12) can be generalized slightly by modifying the limits of inteÂ¬
gration and introducing the time ordering operator, T,
t t t
Ue(tâ€™V7nH)n n! 1 ritl \ dt2 â€¢ â– ' 1 dtn
^ to t0
x exp{r.(t1 + t2[+ . . .  tn )} T[V(t1)V(tz) . . . V(tn) ] . (3.13)
The time ordering operator rearranges the product of perturbation operaÂ¬
tors such that the leftmost term is the latest in chronological order.
The perturbed eigenstates l'*'j(tg)> can now be expressed in terms of
the unperturbed eigenstates by noting that as tgÂ»+Â», 'fj(tg)>,iâ€™g>, and
as tg increases from Â°Â° to zero, the perturbation is "adiabatically
switched on
57
'1,I(0)> = Ue(o,Â»)$0> . (3.14)
According to a theorem of GellMann and Low (1951), if
UÂ£ (o, â€œ>) ] 4â€™q> _  'f j (0) >
llo <*0lue(oÂ»â€œ)lÂ®0> E <^0ITI(Â°)> (3'15)
exists, then it is an eigenstate of II
H11 j (0) > E  T j (0) >
Â«I.glâ€™fjÃOlT = <4.0Â¥j(0)> â€¢ (316)
These results can now be used to determine the electron propagator.
In Chapter 1, the propagator was defined as the ground state average of
a timeordered product of field operators in the Heisenberg representaÂ¬
tion
iGij(t)
(3.17)
Using Eq. (3.15) and the fact that 'l'^>=Â¥j(0)>, this average can be exÂ¬
pressed in the interaction representation as
<4>0Ue(Â«,t)T[ai(t)a!(0)]Ue(t,â€œ)
lGij(t) 
(3.18)
Using the expansion of the evolution operator (Eg. (3.13)) and taking
the limit ef0, it can be shown (Fetter and Walecka, 1971) that
oo , ,
1Gij(t) = nf0(â€˜Ã)n n! J dtl ' ' â€¢ J dtn
â€” OO â€” oo
<0T[V(t1) . . . V(tn)a.(t)a](0)]*0>
<'Ãâ€™gU(
X
(3.19)
58
The final step in the diagrammatic expansion method is to expand
the numerators of each term in Eq. (3.19) using Wick's theorem (Wick,
1950) and to represent them diagrammatically (e.g. Fetter and Walecka,
1971). The denominator of Eq. (3.19) must also be expanded and diaÂ¬
grammed, and when this is done, all disconnected diagrams
arising from the expansion of the numerator will cancel (Abrikosov
et al., 1965).
Formally, the diagrammatic expansion method and the superoperator
formalism appear strikingly dissimilar. The diagrammatic method is
formulated in the causal representation while the superoperator formalÂ¬
ism utilizes the energy representation. The diagrammatic method employs
a pictorial representation of the algebraic structure while the superÂ¬
operator formalism emphasizes the algebraic structure directly. Yet
the primary goal of each formalism is the same: an accurate prediction
of ionization energies and electron affinities. Therefore, the two
formalisms are inherently equivalent. It is our desire in this chapter
to explicitly demonstrate the equivalence between these two formalisms
and to reexamine the superoperator decoupling approximations in terms
of a diagrammatic analysis.
3.2 Perturbation Theory
The unifying feature of the diagrammatic expansion method and the
superoperator formalism is perturbation theory (Born and Ohrn, 1978).
Since the commutator product is distributive with respect to addition,
we can define a partitioning of the superoperator Hamiltonian into an
unperturbed part plus a perturbation,
(3.20)
59
One convenient partitioning, which will be shown to readily yield the
HartreeFock propagator as the unperturbed electron propagator, is the
MdllerPlesset partitioning (Muiller and Plesset, 1934). With this parÂ¬
titioning, H0 has the form
H = E Â£ a'a  h I
0 â€ž r r r , 11 r r1
i r j r
(3.21)
and the perturbation is expressed as
V= Z [l;a V ,a ,a 6 , ,aâ€ža ] + >s E
r,r',s,s' r r s s r s r r s r>r.
(3.22)
Of course, when the commutator product is formed for the superoperators,
the constant term in these definitions will cancel.
Other partitionings of the Hamiltonian may also be assumed and may
lead to superior convergence properties (Claverie et_ aj_., 1967). One
alternative partitioning which has been employed in the perturbation
calculation of correlation corrections to the total energy is the Epstein
Nesbet partitioning (Epstein, 1926, Nesbet, 1955a, 1955b). In propagator
applications, the work of Kurtz and Ohrn (1978) may be roughly interpreted
in terms of a partitioning where the unperturbed Hamiltonian incorporates
all relaxation contributions to the ionization energy. It is difficult
to define this unperturbed Hamiltonian explicitly, but it formally satisÂ¬
fies the eigenvalue equation
HÂ¿ak = AEk(SCF)ak (3.23)
in contrast to
Vk = ckak (3â€˜24'
for the MdllerPlesset partitioning. The method of Kurtz and Ohrn yields
60
excellent ionization energies and electron affinities with a simple
secondorder selfenergy, however it has not been formally analyzed in
detail.
Corresponding to the partitioning of the superoperator Hamiltonian,
we can introduce a partitioning of the operator space defined by the
projection superoperators 0 and P,
O = Iz ak)(ak = a)(a (3.25)
P = I  6 (3.26)
These superoperators operate on elements of the operator space through
the relations
0Xi = X ak)(a<Xi.) (3.27)
PX. = X  0Xi (3.28)
and are idempotent (0^ = O, P^ = P), selfadjoint (0â€˜ = 0, P1 = P), and
mutually exclusive (0P=P0=0). The superoperator 0 projects from an arÂ¬
bitrary operator product that part which lies in the model subspace,
i.e. that part which is spanned by the eigenelements of H^. The superÂ¬
operator P projects onto the orthogonal complement of the model subspace,
i.e. that part which we have no a priori knowledge about.
To obtain a perturbation expansion of the superoperator resolvent,
we consider its outer projection (Lowdin, 1965) onto the model subspace,
G(E) = 6(EIH)16
(3.29)
= 6(eih0v)1o
(3.30)
By iterating the identity
61
(AB)"1 = A'1 + A_1B(AB)_1
the inverse in Eq. (3.30) can be expanded as
(3.31)
G(E) = (EIH0)_1Ã“ + (EIH0)_10V(EIH0)_10
+ (EÃH0)â€œ10V(EÃH0)_1V(EÃH0)~16 + . . .
where the property
Nl 3
(3.32)
(H0,0]_ = 0
(3.33)
has been used. Now since 0 plus P form a resolution of the identity,
each resolvent of Hg occurring between perturbation superoperators, V,
can be rewritten as a sum of its projections on the model subspace and
the orthogonal complement,
(EIHg)'1 = (EIH0)_16 + (EIM0)_1P
Gq(E) + T0(E)
(3.34)
(3.35)
With this notation, Eq. (3.32) becomes
G(E) = G0(E) + G0(E)VG0(E) + GQ(E)V[GQ(E) + TQ(E)]VGQ(E)
+ Gq(E)V[G0(E) + T0(E)]V[G0(E) + TQ(E)]VGQ(E) + . . .
and can be resummed to yield
(3.36)
G(E) = G0(E) + Gq(E)[V + VT0(E)V + VTQ(E)VTQ(E)V + . . . ]G(E) (3.37)
Defining the reduced resolvent of the full superoperator Hamiltonian as
T(E) = P[aO + P(EIH)P]1P
= Tq(E) + Tq(E)VT(E) ,
(tfO)
(3.38)
(3.39)
62
Eq. (3.37) can be written in closed form
G(E) = Gq(E) + G0(E)[V + VT(E)V]G(E) . (3.40)
Alternatively, we can define wave and reaction superoperators through
the equations (cf. Lowdin, 1962, or Brandow, 1967)
W(E) = Ã + T(E)V (3.41)
t(E) = VW(E) . (3.42)
The reduced resolvent, wave, and reaction superoperators introduced
in this section are functions of the superoperators I, Hq, and V and as
a consequence, operate in a more complicated way. To apply a superoperÂ¬
ator function to an operator in the operator space, it must first be exÂ¬
panded in terms of the superoperators I, Hq, and V which are then sucÂ¬
cessively applied to the operator. For example,
W(E)X. = [Ã + T(E)V]X. (3.43)
= [I + Tq(E)V + T0(E)VTQ(E)V + . . . ]X. (3.44)
where
T0(E)VXi = [E~11 + E~2 Hq + E'3 H2 + . . . JPVX^ , (3.45)
etc.
3.3 Equivalence of the Superoperator Formalism and the Diagrammatic
Expansion Method
Eqs. (3.37) and (3.40) represent the superoperator form of the Dyson
equation (Dyson, 1949), and the reaction superoperator (Eq. (3.42)) can
be identified as the selfenergy. To demonstrate that Eq. (3.37)
63
corresponds term by term with the diagrammatic propagator expansion, we
must first form the operator average of G(E) to obtain the matrix Dyson
equation, next evaluate all necessary operator averages, and finally
diagram the resulting algebraic formulae. Owing to the complicated
operator averages that must be evaluated in third and higher orders of
the perturbation superoperator, the equivalence between these two forÂ¬
malisms has only been explicitly demonstrated through third order and
is assumed in all higher orders.
The matrix Dyson equation is obtained by forming the operator averÂ¬
age of G(E) with respect to the basis elements of our model subspace
G(E) = (aG(E)a) (3.46)
= Gq(E) + Gq(E)S(E)G(E) , (3.47)
where
KE) = (aVa) + (aVTQ(E)Va)
+ (aVT0(E)VT0(E)Va) + . . . . (3.48)
Since Hq was chosen to be the Fock superoperator, the appropriate denÂ¬
sity operator to employ in the evaluation of the operator averages is
the HartreeFock density operator. Realizing that the grand canonical
density operator (Eq. (1.44)) reduces to the HartreeFock density operÂ¬
ator when pure state occupation numbers of zero or one are chosen, we
shall employ this density operator.
Beginning with the evaluation of matrix elements for the unperturbed
propagator, G^(E), the HartreeFock propagator is easily obtained (cf.
Section 1.3)
64
VE>ij ' â€˜*jlÂ«EIHo> Â»1>
â– râ€™lajla,) t E_Z
, r\ * rt,.u * r\â€˜
VE>ij â– (Â£.,)%
(3.49)
(3.50)
(3.51)
(3.52)
The evaluation of each term in the selfenergy expansion requires the
initial evaluation of Va..,
Vai = U T. ^rr'I ssl>[ai,arar,as,as]_
r, r j s, s
5 [a . .aras]_
r,s,s'
l a a ,a  l a
i rs s i s s
y^jS^S s, s
(3.53)
(3.54)
With this result, the firstorder term (aj  Va ^) is obtained without much
additional effort
Â£Ã1J(E)ij = (â€œjlV^i) (3.55)
= !Ã¡ 7 Tr{p[aâ€ža ,a ,a .] }
r,s,s' J
l Tr{p[a ,at] } (3.56)
s, s1 J
= >: 6 ,  T. (3.57)
r,s s *
i(1)(E)ij = 0 . (3.58)
When the effective, singleparticle potential used in the unperturbed
problem is the HartreeFock selfconsistent field potential, all singleÂ¬
particle corrections vanish (Bartlett and Silver, 1975a).
65
The evaluation of the second and higherorder selfenergy matrices
requires the evaluation of T0(E)Vai and VT0(E)Va.. The first of these
quantities can be expanded as
T 0(E) V a. = (EIH0)"1PVai
(3.59)
= (EIH0)'1Vai  z (EÃH0)'1ak)(akVa.)
(3.60)
using Eq. (3.26). It follows from the previous result for (a^ Va^) that
the second term in Eq. (3.60) vanishes. The first term can now be evaluÂ¬
ated by expanding the resolvent of Hq and realizing that any operator
product is an eigenelement to Mq, i.e.
Vras'as = (WEs)aras'as '
(3.61)
Consequently, we obtain
T0(E>Vai =  z i(E+VcsEs,)'1 ajas,as
r, s s s
 Z (Ee )1a
s, s1
(3.62)
with the help of Eq. (3.45). The remaining application of V and the
i
average value evaluation is straightforward and yields
Z(2)(E)i:j = (aj\/T0(E)Va.)
(3.63)
= 7. >+1 >
r,s ,s â– (E+Ercses,) r s s' r sâ€™
(3.64)
for the matrix elements of the secondorder selfenergy.
The HartreeFock average is now obtained by choosing occupation
numbers of zero and one. An examination of the occupation number factor
in Eq. (3.64) reveals that with this restriction, it will be nonvanishing
66
only when the summation index r runs over occupied spin orbitals and s
and s' run over unoccupied spin orbitals or when r runs over unoccupied
spin orbitals and s and s' run over occupied spin orbitals. Denoting
a, b, c, . . . as summation indices over occupied spin orbitals; p, q,
r, . . . for unoccupied spin orbitals; and i, j, k, . . . for unspecified
spin orbitals, Eq. (3.64) can now be written as two terms which involve
restricted spin orbital summations
y(2)m = ,, r
Â¿ U ij 2 L (E+e e ~e )
a.p.q v a p V
+ 1. v
P.a.b
(3.65)
The conversion of Eq. (3.65) into diagrams is a straightforward proÂ¬
cedure for which we shall use the rules and diagram convention of Brandow
(1967) and Bartlett and Silver (1975b). This convention represents the
synthesis of the antisymmetrized vertices of the Uugenholtz (1957) or
Abrikosov (1965) diagrams with the extended interaction lines of the
Goldstone (1957) diagrams, and the rules for constructing these diagrams
are given in Table 5. The application of these rules to the terms in
Eq. (3.65) yields the following diagrams:
j. l
a,p,q ^VVV
1.0
y
(E+e E Â£,)
p a b
p, a , b
(3.66)
(3.67)
These diagrams are precisely the same as those obtained in the second
order diagrammatic expansion after a Fourier transformation into the
energy representation (Cederbaum and Domcke, 1977).
67
Table 5. Rules for Constructing SelfEnergy Diagrams.
1. Each antisymmetrized twoelectron integral factor in the
numerator is represented by an interaction line with a vertex
(dot) at both ends. The number of interaction lines denotes
the order of the term.
2. Using the Dirac braket notation, both indices in the bra are
represented by lines leaving a vertex while those of the ket
are represented by lines entering a vertex. There must be only
one outgoing and one incoming line per vertex, therefore, assign
the index of electron coordinate one to the left vertex and the
index of electron coordinate two to the right vertex of each
interaction line.
3. Summation indices running over hole states (occupied orbitals)
are directed downward, indices running over particle states
(unoccupied orbitals) are directed upward, and external indices
(not summed) are drawn horizontally.
To Check Diagrams:
4. The energy denominator of the diagrammed expression should be
obtained by first connecting the external lines and assigning
a factor of E to this directed segment. Second, imagine horiÂ¬
zontal lines drawn between each pair of interaction lines.
Each horizontal line corresponds to a multiplicative, denominator
factor obtained by summing the orbital energies of each hole
(downgoing) line that intersects it minus the sum of orbital
energies for particle (upgoing) lines that intersect it. Treat
the factor E of the connected external lines as an orbital energy.
5. Numerical factors should be obtained by assigning a factor of %
for each pair of equivalent internal lines. Equivalent internal
lines are two lines which begin on the same interaction line,
end on the same interaction line, and go in the same direction.
6. The overall sign factor should be obtained by assigning a factor
of minus one to each internal hole line segment and a minus one
to each closed loop.
68
The evaluation of the thirdorder selfenergy matrix is similar to
the secondorder matrix but much more tedious and the result is presented
in Appendix 3. As was done for the secondorder expression, the occupaÂ¬
tion numbers must again be restricted to zero and one to obtain the
HartreeFock average. When this restriction is made, the unrestricted
spin orbital summations in Appendix 3 will reduce to summations involving
occupied, unoccupied, and unspecified spin orbitals. Using the algebraic
identity
a a r i i
(Ea) (EbT ' Ta^bJ [ " TTO
it is possible to combine terms in such a way that expressions involving
only occupied and unoccupied spin orbital summations are obtained. These
expressions are presented in Appendix 1. The corresponding diagrams in
Appendix 1 again are precisely those occurring in the thirdorder, diaÂ¬
grammatic selfenergy expansion.
3.4 Diagram Conserving Decoupling
The wave and reaction superoperators identified with the help of
perturbation theory in Section 3.2 have special importance in the developÂ¬
ment of decoupling approximations for the electron propagator. As we
have already seen, the reaction superoperator generates the diagrammatic
selfenergy expansion. A truncation of this expansion offers one viable
decoupling scheme. The wave superoperator, on the other hand, has the
property of generating eigenelements to the full superoperator HamiltoÂ¬
nian from the eigenelements of the unperturbed superoperator Hamiltonian
(EIH)W(E)a = 0 .
(3.68)
(3.69)
69
This property is easily proven by first using Eq. (3.41) to expand W(E)
and then premultiplying both sides of Eq. (3.69) by P
P(EIH)W(E)a = P(EIH)a + P(EIH)T(E)Va . (3.71)
Using the identity
P(EIH)T(E) = P (3.72)
and the property Pa_ = 0, Eq. (3.71) simplifies to
P(EIH)W(E)a =  PVa + PVa = 0 (3.73)
which implies the validity of Eq. (3.69)
It is of interest at this point to show a connection between the
superoperator formalism and the Equations of Motion (EOM) method for
determining ionization energies (Simons and Smith, 1973). In this
method, one seeks solutions of the equation
[H,Q]_ = ojQ (3.74)
which is precisely Eq. (3.69). Here the operator Q is interpreted as a
correlated ionization operator that generates, in principle, the exact
(Nl)electron ion states from the exact Nelectron reference state.
One approach to solving Eq. (3.74) involves the application of Rayleigh
Schrodinger perturbation theory (Dalgaard and Simons, 1977). By partiÂ¬
tioning the Hamiltonian operator, expanding both the ionization operator,
Q, and the ionization energy, io, in terms of a perturbation parameter,
and collecting terms of the same order, a set of perturbation theory
equations are obtained. The solution of these equations yields an expanÂ¬
sion for Q which is analogous to the superoperator equation
70
h = W(E)a . (3.75)
The only difference is that E is replaced by which is a consequence
of using RayleighSchrodinger rather than BrilliounWigner perturbation
theory.
Returning now to the inner projection of the superoperator resolvent,
G(E) = (ah)(h(EIH)h)1{ha) (3.76)
we may view Eq. (3.75) as an alternative prescription for choosing the
inner projection operator manifold. Recalling from Section 1.6 that since
the density operator describing the unperturbed (model) problem does not
commute with the full Hamiltonian, the operator scalar product will not
in general exhibit Hermitian symmetry. Consequently, we define
(h=(aW+(E) (3.77)
and note that
(a_W+(E) f (W(E)aJ . (3.78)
Approximate electron propagator decouplings can now be obtained by
truncating the expansion of the wave superoperator,
W(E) = I + Tq(E)V + TQ(E)VT0(E)V + . . . . (3.79)
Truncation of this expansion, with only the superoperator identity, trivÂ¬
ially yields the HartreeFock propagator, therefore we next consider
W(E) = Ã + T0(E)V .
doting that the subspaces (a^J and (fjf^ = T^(E)Va^} are mutually
orthogonal, Eq. (3.76) can be readily solved for G~*(E)
(3.80)
/I
G."1 (E) = G^(E)  Â£(E) , (3.81)
where
7(E) = (aVT0(E)Va)(a.VT0(E)(EÃH)T0(E)Va)'1(aVT0(E)Va) (3.82)
Making the following identifications from Section 3.3:
(aJVT0(E)Va) = E(2)(E) , (3.83)
(a)VT0(E)(EÃH0)T0(E)Va) = Â£(2)(E) , (3.84)
and
(aVT0(E)VT0(E)Va) = Â£(3)(E) , (3.85)
Eq. (3.82) can be rewritten
Â£(E) = Â£^2^(E)[^2^(E)  I^3)(E)]'V2^(E) . (3.86)
Expanding the inverse of Eq. (3.86), we easily see that this selfÂ¬
energy approximant coincides with the diagrammatic expansion through
third order but additionally yields contributions to all higher orders.
If the exact selfenergy is rewritten as a moment expansion in terms of
a perturbation parameter, A,
X_1E(E) = z Ak(aV(T (E)V)ka) , (3.87)
k=0 U
we see that Eq. (3.86) represents the [1,1] Pade1 approximant to this
expansion. Owing to the close connection between Pade' approximants and
the inner projection technique as demonstrated in Chapter 2, this result
is not surprising. These Pade' approximants to the selfenergy, however,
will have entirely different convergence properties than those studied in
Chapter 2.
72
3.5 Approximations and App1ications
Computational applications of the [1,1] Pade' approximarit to the
selfenergy require the evaluation of the second and thirdorder selfÂ¬
energy matrices. The secondorder matrix is relatively easy to evaluate.
The thirdorder matrix, on the other hand, is exceedingly more difficult
and can presently be only approximately calculated without excessive
computational effort. An examination of the formulae in Appendix 3 reÂ¬
veals that unlike the fourth moment matrix in the moment conserving
decoupling, the thirdorder selfenergy matrix is energy dependent.
This additional complication makes the partial summation technique used
in the moment conserving decoupling ineffectual since the thirdorder
selfenergy matrix will generally need to be resummed with different
values of E hundreds of times in the search for poles of the propagator.
The first approximation that we will examine is the complete neglect
of the thirdorder selfenergy matrix. With this approximation, the
[1,1] approxiniant in Eq. (3.86) reduces to a secondorder truncation of
the diagrammatic selfenergy expansion,
Z.(E) = E^(E) . (3.88)
This secondorder selfenergy approximation is interesting not only beÂ¬
cause it contains the most important relaxation and correlation correcÂ¬
tions to Koopmans' theorem (in a perturbation theoretical sense), but
also because it exhibits the same analytic form as the exact selfenergy
(Hedin and Lundqvist, 1969, Cederbaum and Domcke, 1977). Furthermore,
since several secondorder, ionization energy calculations have been
reported in the literature, this approximation will afford both a conÂ¬
venient check of new computer code and the computational experience
necessary to implement more refined approximations.
73
The first computational application of this decoupling approximation
was to the water molecule using the same basis and internuclear geometry
as described in Section 2.4. The results of this calculation are preÂ¬
sented in Table 6 along with the Koopmans' theorem, AE(SCF), and experiÂ¬
mental values for the ionization energies. Two ionization energies have
been tabulated for the 2a^ ionization with their corresponding pole
strength (r^ of Eq. (1.70)) in parentheses. The occurence of two, relaÂ¬
tively strong propagator poles for this ionization represents a breakdown
in the quasiparticle description of inner valence ionizations (Cederbaum,
1977) and makes assignments of principal and shakeup ionizations ambig
â– k
uous. In general, the secondorder ionization energies are quite enÂ¬
couraging and represent significant improvements to each of the Koopmans'
values. Furthermore, these results are comparable in accuracy to the
AE(SCF) results but possess the convenience of being obtained in a single
calculation whereas the AE(SCF) results required six separate Hartree
Fock calculations.
The relatively poor agreement of the 3a^ and lb^ ionization energies
with the experimental values in Table 6 seems attributable to basis inÂ¬
completeness. Despite the lack of polarization functions, this suspicion
is supported by the facts that the 3a^ orbital is the highest occupied
orbital in that symmetry and that this basis contains only two contracted
Gaussian orbitals of symmetry. In order to study the basis dependence
of the secondorder selfenergy approximation, two additional calculations
*The ESCA spectrum of the water molecule (Siegbahn et al., 1969)
substantiates this phenomenon since the 2a, peak is quite broad and asymÂ¬
metric. Experimentally, it appears that the lower energy ionization
should have a larger pole strength (in contrast with the results of Table
5) since the peak is skewed to higher binding energies.
Table 6. Principal Ionization Energies of Water Computed
with the 14 CGTO Basis.
Orbital
Koopmans
AE(SCF)a
(2)
2(E)
Exp.b
lal
559.4
540.8
539.4
540.2
2a i
37.0
34.6
34.0 (.61)
32.6 (.28)
32.2
3d 2
15.4
13.0
12.9
14.7
ibi
13.8
11.0
10.8
12.6
lb2
19.5
17.8
18.1
18.6
aGoscinski ert al_. (1975).
bSiegbahn et al_. (1969).
75
were performed with larger basis sets. The first of these calculations
employed a 26 contracted orbital basis which augmented the original 14
orbitals (Table 3) with a set of porbitals on the hydrogen atoms and a
set of dorbitals on oxygenall with unit exponents. The HartreeFock
total energy obtained with this basis was E(HF)= 76.0459 H. The second
calculation employed a 38 contracted orbital basis which included all of
the orbitals in the 26 orbital basis plus an additional set of diffuse
porbitals on the hydrogen atoms (a = 0.25) and a set of diffuse dorbitÂ¬
als on oxygen (a = 0.40). This basis yielded a HartreeFock total energy
of E(HF)= 76.0507 H.
The most significant propagator poles calculated in the valence
region (0~40 eV) with each of the three water basis sets are presented in
Table 7 along with the secondorder results of Cederbaum (1973a). The
inclusion of polarization functions not only improves the 3a^ and lbj
ionization energies, it also reverses the relative pole strengths of the
two dominant 2a^ poles bringing the theoretical results into better agreeÂ¬
ment with experimental observations (see footnote on page 73).
Cederbaum's secondorder results were obtained with a basis comparable
in size and quality to the 26 orbital basis in Table 7. He deletes
several virtual orbitals from this basis before computing the ionization
energies, however. This approximation may account for the small discrepÂ¬
ancies between his results and those reported here.
The formaldehyde molecule was chosen for a second application of
the secondorder selfenergy approximation. Ionization energies were
calculated using two basis sets. The first consisted of Huzinaga's 9s,
5p primitive basis sets for oxygen and carbon (Huzinaga, 1965) contracted
to 4s and 2p functions with Dunning's contraction coefficients (Dunning,
76
Table 7. Basis Set Effects on the Ionization Energies of Water
Computed with a SecondOrder SelfEnergy Approximation.
Symmetry
14 CGTO's
26 CGTO's
38 CGTO's Ced.
al
36.5
(.005)
37.1
(.003)
34.0
(.607)
33.4
(.288)
33.2
(.231)
32.6
(.279)
32.1
(.592)
31.9
(.628)
32.9
12.9
(.913)
13.4
(.908)
13.5
(.903)
13.2
bl
34.9
(.005)
35.1
(.005)
10.8
(.909)
11.1
(.904)
11.2
(.900)
10.9
40.6
(.003)
40.8
(.004)
18.1
(.931)
18.0
(.922)
18.0
(.919)
17.7
E(HF)
76.0082
76.0459
76.0507
76.0419
aCederbaum (1973a).
77
1970). The orbital exponents of Huzinaga's 4s primitive basis for hydroÂ¬
gen were scaled by a factor of 1.2, and the resultant orbitals were conÂ¬
tracted to 2s functions as recommended by Dunning. The complete basis
appears in Table 8. The second basis augmented the first by the addition
of one set of porbitals on the hydrogen atoms and one set of dorbitals
on both the oxygen and carbon atoms. Unit exponents were chosen for the
porbitals on hydrogen while exponents of 0.8 were chosen for the dorbit
als. One and twoelectron integrals were computed with the MOLECULE
program (Almlof, 1974) at the experimental equilibrium geometry: R(C0)=
2.2825 a.u., R(CH)= 2.1090 a.u.,and )(HCH)= 116.52Â° (Oka, 1960, Takagi
and Oka, 1963), and the HartreeFock calculations and twoelectron inteÂ¬
gral transformations were performed with GRNFNC (Purvis, 1973). The
HartreeFock total energy for the smaller, 24 orbital basis (no polariÂ¬
zation) was E(HF)= 113.8257 H., and for the larger, 42 orbital basis
(with polarization) E(HF)= 113.8901 H. The HartreeFock orbital energies
and secondorder selfenergy results for both basis sets are presented
in Table 9 for the principal ionizations along with the secondorder
results of Cederhaum et al_. (1975) and the experimental values.
The results in Table 9 typify two general features of ionization
energy calculations. The first is that Koopmans1 theorem yields values
which are usually higher than experimental ionization energies. Second,
the inclusion of secondorder relaxation and correlation corrections
generally overcorrects the Koopmans' estimate and yields values which
are usually lower than experiment. For several ionizations in Table 9,
the secondorder deviations from experiment are as large as the Koopmans'
values only opposite in sign. Although it is possible that the larger,
polarized basis used in the second calculation may still lack adequate
Table 8. Contracted Gaussian Basis for Formaldehyde.
Carbon
s sets
Oxygen
s sets
Exponents
Contraction
Coefficients
Exponents
Contraction
Coefficients
4232.6100
0.002029
7816.5400
0.002031
634.8820
0.015535
1175.8200
0.015436
146.0970
0.075411
273.1880
0.073771
42.4974
0.257121
81.1696
0.247606
14.1892
0.596555
27.1836
0.611832
1.9666
0.242517
3.4136
0.241205
5.1477
1.000000
9.5322
1.000000
0.4962
1.000000
0.9398
1.000000
0.1533
1.000000
0.2846
1.000000
Carbon p sets
Oxygen p sets
Exponents
Contraction
Coefficients
Exponents
Contraction
Coefficients
18.1557
0.018534
35.1832
0.019580
3.9864
0.115442
7.9040
0.124189
1.1429
0.386206
2.3051
0.394727
0.3594
0.640089
0.7171
0.627375
0.1146
1.000000
0.2137
1.000000
Hydrogen
s sets
Exponents
Contraction
Coefficients
19.2406
0.032828
2.8992
0.231208
0.6535
0.817238
0.1776
1.000000
79
Table 9.
Principal
Ionization
Energies for
Formaldehyde.
24
42
(2)
(2)
Orbital
Koopmans
3(E)
Koopmans
2(E)
Ced.a
Exp.
lal
560.12
538.93
559.81
538.62

539.43b
2al
309.09
297.27
308.87
296.90

294.21b
3a 
38.94
33.66
38.18
32.56

34.2C
4a i
23.39
20.97
23.30
21.03

21.15c
5a ^
17.29
13.98
17.38
14.38
14.42
16.2d
lbl
14.56
13.83
14.45
13.72
13.50
14.5d
lb2
19.47
17.16
19.08
17.07
16.63
17.0d
2b2
12.06
9.04
11.93
9.30
9.25
10.9d
E(HF)
113.8257
113.8901
113.9012
aSecond order results of Cederbaum et al_. (1975).
bJolly and Schaaf (1976).
cHood et a]_. (1976).
^Estimated center of gravity (Cederbaum and Domcke, 1977) from spectrum
of Turner et al. (1970) .
80
polarization functions, the rather large discrepancies between the
secondorder results and experiment more probably indicate that third
(and higher) order selfenergy corrections are now sizable. The general
conclusion that a secondorder selfenergy approximation is inadequate
for an accurate calculation of ionization energies has been previously
concluded by Cederbaum (1973b) and necessitates a reexamination of the
approximation made in Eq. (3.88).
Rather than completely neglecting the thirdorder selfenergy matrix,
let us now consider an approximation that includes at least part of these
contributions. Which thirdorder selfenergy diagrams should be included?
There are two wellestablished results that are relevant to this quesÂ¬
tion: Studies of the electron gas model have shown that in the limit of
high electron density, the socalled ring diagrams dominate the selfÂ¬
energy expansion (Pines, 1961), while in the limit of low electron denÂ¬
sity, the socalled ladder diagrams dominate (Galitskii, 1958). In order
to determine whether atomic and molecular selfenergies can be approxiÂ¬
mated by specific thirdorder diagrams (e.g. rings or ladders), we need
to evaluate all thirdorder diagrams for some representative systems.
Cederbaum (1975) has done this for several simple systems and has found
that both ring an_d ladder diagrams dominate the thirdorder selfenergy.
This result implies that atoms and molecules lie somewhere between the
high and low density extremes. It is therefore essential to include
both ring and ladder diagrams in any thirdorder selfenergy approximaÂ¬
tion. These diagrams are
rings ladders
(3.89)
81
and correspond to the algebraic expressions labeled A, B, C, and D in
Appendix 1.
We include six additional diagrams in our thirdorder selfenergy
approximation because of the computational efficiency with which they
are evaluated. These diagrams are the energy independent diagrams
corresponding to expressions MR in Appendix 1. For these six diagrams,
it is feasible to employ the partial summation technique since they must
be evaluated only once.
Approximating the full thirdorder selfenergy matrix by only ring,
ladder, and constant energy diagrams, let us now consider the solution
of the Dyson equation with the [1,1] Pade1 approximant to the selfenergy
expansion. Owing to the fact that the inner projection manifold from
which the [1,1] approximant was derived is energy dependent (Eq. (3.80)),
the simple analytic form of the selfenergy eigenvalues, illustrated in
Fig. 2, is lost. Furthermore, the selfenergy poles are now given by
det (2{2^(E)  2^(E)) = 0 (3.91)
rather than by an eigenvalue problem and are consequently more difficult
to obtain. For these reasons, the pole search described in Chapter 1
and used with the secondorder selfenergy approximation is no longer
an efficient or reliable procedure. An alternative method of solution
used in the following applications was to use the HartreeFock orbital
energy as an initial guess to the propagator pole and to iterate Eq.
(1.67) to convergence. When convergent, this procedure invariably yields
82
a principal propagator pole and its corresponding pole strength. AlÂ¬
though the [1,1] selfenergy approximant does not quarantee a positive
pole strength, this was never a problem in any of the calculations reÂ¬
ported here.
The principal ionization energies for the water molecule were calÂ¬
culated using both the 14 and 26 CGTO basis sets in order to evaluate
the [1,1] Pade' approximant to the selfenergy expansion, and the results
appear in Table 10. The most significant feature of these results is that
each ionization energy has been shifted from its second order value to
higher energy and is now in better agreement with the experimental value.
It is further noticed that the valence ionization energies are still
smaller than the experimental values while the la^ (core) ionization
energy is now larger than experiment. Apparently, the diagrams included
in the thirdorder selfenergy matrix overestimate the actual relaxation
and correlation effects for this ionization.
Convergence difficulty was experienced for the 2a^ ionization
energy using the 14 orbital basis. A schematic plot of W9 (E) is pre
2a i
sented in Fig. 3 and reveals that there are no propagator poles in this
energy region. This anomaly is no doubt a consequence of some quirk in
the basis since the 26 orbital basis yields a very accurate 2a^ ionizaÂ¬
tion energy.
3.6 Evaluation of the Diagram Conserving Decoup 1ing
The algebraic structure of the superoperator formalism has been
successfully exploited in this chapter to yield several new insights into
the decoupling problem. The application of perturbation theory has
demonstrated that the electron propagator equation of motion can be
83
Table 10. Comparison of Principal Ionization Energies for Water
Obtained with the SecondOrder and the [1,1] SelfEnergies
Using the 14 and 26 CGTO Basis Sets.
14
26
Orbital
(2)
Â£(E)
[1,1]
(2)
S(E)
[1,1]
Exp.a
lal
539.4
541.6
539.2
540.9
540.2
2al
34.0 (.607)
32.6 (.279)
no convergence
(see text)
32.1
32.2
32.2
3a j
12.9
13.4
13.4
13.6
14.7
Ibi
10.8
11.1
11.1
11.3
12.6
lb2
18.1
18.4
18.0
no results
18.6
aSiegbahn et al. (1969).
Figure 3. A sketch of WÂ£ai in the energy region of the
2aj ionization ^obtained with the [1,1] self
energy approximant using the 14 CGTO basis.
85
1.215
 1.175
86
resummed to yield the equivalent of the diagrammatic expansion. This
resummation also allows the identification of wave and reaction superÂ¬
operators which have special importance in decoupling approximations.
We have shown that when the inner projection manifold of the superoperÂ¬
ator resolvent is chosen to consist of the firstorder truncation of the
wave superoperator, a [1,1] Pade' approximant to the selfenergy expanÂ¬
sion is obtained. This approximant is correct through third order and
contains a geometric approximation to all higher orders. In general,
the Nthorder truncation of the wave superoperator will yield an [N,N]
Pade1 approximant which is correct through the (2N+l)st order in the
selfenergy expansion. One final insight afforded by this decoupling
is the realization that electron correlation can be described exclusively
in the operator space. We argued in Chapter 1 that when the propagator
was defined as a singletime Green's function, the density operator was
arbitrary. We have now demonstrated in this chapter that any desired
order in the selfenergy expansion may be obtained using as a specific
choice,the uncorrelated, HartreeFock density operator.
Computational applications of the diagram conserving decoupling
have been encouraging. These applications have confirmed previous conÂ¬
clusions (Cederbaum, 1973b) that a secondorder selfenergy is generally
inadequate for obtaining accurate ionization energies. It is important
if not essential that thirdorder ring and ladder diagrams be included
in any selfenergy approximation although the errors arising from basis
incompleteness may be of equal magnitude and hence cannot be ignored.
The inclusion of the thirdorder ring, ladder, and constant energy diaÂ¬
grams in the [1,1] selfenergy approximant has succeeded in improving the
secondorder results but even these results are not consistently better
than the Koopmans1 theorem values.
8;
One important feature of the [1,1] selfenergy approximant is that
even though it is constructed from only the second and thirdorder selfÂ¬
energy matrices, it contains a geometric approximation of all higher
orders in the selfenergy expansion. Certainly, some fourth or higher
order terms may be just as important as thirdorder terms; therefore,
this approximation is highly desirable. The fourth and higherorder
terms arising from the [1,1] selfenergy approximant, however, are not
readily analyzed diagrammatically. In fact, being a purely algebraic
approximation, the [1,1] approximant may not yield any valid fourth or
higherorder diagrams. Given the fact that ring and ladder diagrams
dominate the thirdorder selfenergy matrix, one can argue that they may
also dominate the higher orders of the selfenergy expansion. An approÂ¬
priate modification of this decoupling scheme might then allow the sumÂ¬
mation of these specific diagrams in all orders. Approximations of this
type are referred to as renormalized decouplings and are examined within
the superoperator formalism in the next chapter.
CHAPTER 4
RENORMALIZED DECOUPLINGS
4.1 Renormalization Theory
In Chapter 3, we tacitly assumed that the application of perturbaÂ¬
tion theory to the calculation of ionization energies and electron
affinities was valid and that the resulting selfenergy expansion was
convergent. Historically however, it was discovered that in both the
nuclear manybody problem and the electron gas model, the simple selfÂ¬
energy expansions are divergent. In order to remove these divergencies,
it is necessary to sum certain appropriate classes of diagrams to all
orders. This method of partial summations is known as renormalization
theory (see e.g. Kumar, 1962 or Mattuck, 1967) and may be viewed as an
analytic continuation of the perturbation expansion. Although a variety
of renormalization procedures exist, such as propagator renormalizations,
interaction renormalizations, and vertex renormalizations, the distincÂ¬
tions mainly depend on the types of diagrams included in the partial
summation and are not particularly important for our consideration.
One renormalization that we are already familiar with is the [1,1]
selfenergy approximant derived in the preceding chapter. In fact, any
rational selfenergy approximant may be regarded as a renormalization
since its geometric expansion will approximate all orders of the perturÂ¬
bation expansion. One problem encountered with the [1,1] approximant
and that occurs in general for rational approximants derived via purely
algebraic considerations is that their geometric expansions may contain
88
89
no readily identifiable diagrams (at least beyond the lowest orders).
Since specific diagrams often dominate the selfenergy expansion (such
as ring and ladder diagrams for atoms and molecules) it is valuable to
investigate whether the superoperator formalism can be adapted to yield
renormalized selfenergy expressions that sum specific diagrams. The
solution as we shall see is rather simple and involves a restriction in
the types of operator products allowed to span the orthogonal complement
of the model subspace. As a specific example, the two particleone hole
TammDancoff approximation (2ph TDA), (Schuck et aj_., 1973, Schirmer
and Cederbaum, 1978), is derived from an effective interaction which is
logically obtained by a projection of the perturbation superoperator onto
the subspace spanned by 2ph type operators (Born and Ohrn, 1979).
Finally, the diagonal approximation to the full 2ph TDA selfenergy
previously derived and applied to the calculation of ionization energies
is shown to neglect terms which, in fact, are diagonal and are necessary
to prevent an overcounting of all diagrams containing diagonal ladder
parts.
4.2 Derivation of the 2ph TDA and Diagonal 2ph TDA Equations
Recalling some of the results of the previous chapter, we had obÂ¬
tained the matrix Dyson equation
G(E) = Gq(E) + EgfEWEME) , (4.1)
where the selfenergy matrix, T_(E), had the following expansion
1(E) = (aV+VT0(E)V+VT0(E)VT0(E)V+ . . . a)
(4.2)
90
Introducing the reduced resolvent of the full superoperator Hamiltonian,
T(E), which is just a projection of the superoperator resolvent on the
orthogonal complement
T(E) = P[aO+(EÃH0)PPVP]_1P , (4.3)
the selfenergy expansion was written in closed form
2(E) = (aJVa) + (aJVT(E)Va) . (4.4)
It was further shown that when the grand canonical density operator is
used to evaluate the operator averages, the firstorder term vanishes.
When P is the exact projector of the orthogonal complement,
P = Ã  6 , (4.5)
the term PVP in Eq. (4.3) is responsible for generating the operator
products that span this subspace. The expansion of this term from the
inverse and its repeated application in each order of the perturbation
expansion yields larger and larger operator products which are only
limited by the number of electrons in the reference state. If instead
of allowing all possible operator products, we restrict them to some
simple types which occur in each order, it may be possible to identify
and sum specific diagrams in all orders of the perturbation expansion.
The restriction of the operator products in the orthogonal compleÂ¬
ment is achieved by approximating the orthogonal projector as
P = f)(f (4.6)
where the manifold {f} contains the desired operator products. The proÂ¬
jector P now has the effect of projecting from the perturbation expansion
91
only those operator products which lie in the subspace spanned by {f}.
The approximation to P in Eq. (4.6) must of course preserve the properÂ¬
ties of the exact projector and should be idempotent, selfadjoint, and
orthogonal to 0.
Our previous experience with the operator product decouplings
suggests that the set of triple operator products {a^a^a^} be chosen
as a first approximation to P. There is a stronger motivation for using
this operator product, however. If the thirdorder ring and ladder diaÂ¬
grams in Eq. (3.70) are examined, it can be seen that between any two
interaction lines there occurs only two particle lines (upgoing) and one
hole line (downgoing) or vice versa. This implies that the intermediate
or virtual states that are represented by these diagrams consist of only
2ph or 2hp excitations of the reference state. Both of these excitaÂ¬
tions are described with the tripleoperator products.
The set of triple products {a^a^a^} is not orthogonal to the simple
operators of the model subspace, hence these two subspaces must be or
thogonalized. Using the GramSchmidt orthogonalization procedure (see
e.g. Pilar, 1963), we define
f i, = N. ,2 [a, a, a E(a la. a .a )a ]
klm klm1 k 1 m ' n1 k 1 nr nJ
l +
1. j [a, a,a +6, a,6, ,a ]
klm1 k 1 m km k 1 kl km
(4.7)
(4.8)
where
N. , = +
klm k k 1 k ni 1 m
(4.9)
The projector
P = E 
k,l ,m
fklm^fkliJ
l
(4.10)
92
is now idempotent, selfadjoint, and orthogonal to 0. The projection of
the perturbation superoperator on the 2ph subspace, Pl/P, which occurs in
Eq. (4.3) can now be regarded as an effective interaction. The expansion
of Eq. (4.3) with P defined as in Eq. (4.10) should yield all diagrams
containing 2ph and 2hp excitations of the reference state.
The necessary operator averages needed to evaluate PVP are:
(4ID
and
(ak'al'Vl5akalaJ = Nk'l'm'{6kk'(1)
)<5^, ()
+<5m'1()+6nil,()
+6kl+5krn}
+Nklm{6m1kâ€™+51,k,} (4.12)
Substituting these expressions and performing some cancelation yields:
E
c',1
1
N'!f N
, k 1 in
' 2
k' 1
'm
.6kk
i(l
> 
m
11 'm' ^klm
1 '<
m
E
â– ' >1
1
tcf m
, klm
k'l
'm
,
1mkm'>6^ ^
â– (

m>)lfk'l
'm'^kltJ
1 '<
mâ€˜
?:
;M
1 ,m'
N''5 N
, klm
k'l
'm
,
11 i Ik1'>6
1 1 mm
,(

l>)lfk'l
'm' ^kklnJ
1 '<
m*
E
M'
3m'
N,'5 W
klm
1,
'2
k'l
'nr
,
'mllkl'.dim
â– (

fk111
'm' ^klrJ
1 '<
ni'
E
M '
sm*
n:â€™5 N
klm
2
k'l
1 m'
1 ! km'>6ml
(

l>)lfkâ€˜T
m' ^kkliJ
1 < m 1 1 < m1
(4.13)
93
Additional simplification can be achieved at this point by anticommuting
the operators aj, and am, in I f^, j ,mâ– ) Â°f the last two terms in Eq. (4.13)
with the appropriate change of sign,
If
k' 11 m
If
k' m11
(4.14)
After interchanging dummy indices 1'+>â– mâ€™, the fourth and fifth terms beÂ¬
come equal to the second and third terms, respectively, and since the
diagonal terms l'= m1 vanish, the summations with 11< m' and 1'> m1 can
be combined as unrestricted summations over 1' and m'. Since the first
term in Eq. (4.13) is obviously symmetric in 1' and m1, the restriction
l'< m' in that term may be removed by multiplying the sum by a factor
of The remaining restriction, 1
ing another factor of ^ since
PVP = E E NkÃmNk,l'm,Ã,5,5kk,(1"")
K j I jlil K ^ I )III
1 < m
6,,, ()6 , ()} I f, ,, , ,)(f. , I
11 11 k m ' 11 mm 1 k 1 ' 1 k 1 m k1m1
(4.15)
is symmetric in these indices as can be verified by interchanging l++m
and relabeling dummy indices 1 '<+ m1. Expanding the ketbra superoperator
(4.16)
'S ' f I f V 11 >m1) (f
k,l,m k1,11,m1 K 1 m
klnr
out of the inverse, evaluating the remaininq operator averages, and
resumming the expansion yields the 2ph TDA selfenergy
2ph TDA
>:(E)iJ 
k,1,m k1,11 ,m
v, ,NklrnNkTm*<1kl llm>
x {(f j (EIH0)f)(f Vf)}j[)(k,rm,
(4.17)
where
94
))} (4.18)
Although different in appearance, this selfenergy expression is formally
the same as that obtained by Purvis and Ohrn (1975a) using the operator
product decoupling and by Cederbaum (1975) and Schirmer and Cederbauin
(1978) using the diagrammatic method. The present derivation clearly
illuminates the parallelism between the two formalisms.
Owing to the large dimension of the {f, , } operator subspace and
the associated difficulty in diagonalizing (fVf), computational appliÂ¬
cations of the 2ph TDA have usually involved additional approximations.
One approximation which has facilitated computational applications is
known alternatively as the shifted Born collision (SBC) approximation
(Purvis and Ohrn, 1974, 1975a) or the diagonal 2ph TDA (Cederbaum,
1974, 1975 and Cederbaum and Domcke, 1977). This purportedly "diagonal"
approximation restricts the spinorbital summation indices in Eq. (4.13)
to kâ€™ = k, l' = l, and m'=m thereby neglecting the last two summations
and yielding the following selfenergy expression
(4.19)
where
A=(l)()(). (4.20)
mi Km k i
By neglecting the last two summations in Eq. (4.13), however, this approxÂ¬
imation actually neglects some diagonal contributions to (fkâ– t *mâ– IVfkim)â€¢
As we have explicitly demonstrated in the derivation of Eq. (4.17) and
95
(4.18), the 2ph TDA selfenergy sums are symmetric in both 1,m and
1',m'; consequently, the last two summations in Eq. (4.13) contain preÂ¬
cisely the same contributions as the second and third summations, reÂ¬
spectively. If the diagonal approximation (k'=k, l'=l, m'=m) is made in
Eqs. (4.17) and (4.18), this symmetry is properly accounted for and the
resulting selfenergy expression is
(4.21)
where
A=!2^ml  ml>(l)()(). (4.22)
ni i k m k i
Eq. (4.22) differs from Eq. (4.20) by the factor of % in the first term.
The inclusion of this factor in the diagonal 2ph TDA is necessary to
prevent an overcounting of diagrams with diagonal ladder parts (see next
section) and typically shifts ionization energies 0.30.4 eV higher in
energy (Born and Ohrn, 1979).
4.3 Diagrammatic Analysis
In order to determine precisely which selfenergy diagrams are inÂ¬
cluded in the 2ph TDA selfenergy, it is necessary to expand the effecÂ¬
tive interaction matrix, (Â£V_f), from the inverse in Eq. (4.17) and diaÂ¬
gram the resulting algebraic expressions in each order. The expansion
of Eq. (4.17) yields the following terms in lowest orders
2ph TDA
7(E),. = 8 >: >: N
1J U 1 m L ' 1 ' m 1
k, 1, m k1,1 ', m'
96
+ Z
K,X,P
%Y?ltVv:?IVrv
+
â€¢ ]â– <1 'm1   jk1 > .
(4.23)
As was done in Section 3.3, the terms in Eq. (4.23) can be simplified by
first restricting the summations over all spin orbitals to summations
over occupied or unoccupied spin orbitals such that the occupation
1. U
number factor NklmNk'rm' is nonvanishing. Doing this in the first term
of Eq. (4.23) and then summing over the delta functions yields
is 7
a,p,q
(E+e e E )
a p q'
(4.24)
+H 7
p,a,b
(E+e e e, )
p a b;
(4.25)
which are the same secondorder selfenergy diagrams as obtained in Eqs.
(3.66) and (3.67).
Restricting and spin orbital summations in the second term of Eq.
(4.23) yields the following expressions
h z
a,p,q
7
b,r,s
(fbrslÃfaPq)
(4.26)
{!VpV
lEVfSy
+!S Z
a,p,q
7
r,b,c
^rbcl^apq)
(4.27)
+H Z
p,a,b
7
d,q,r
(fcqr!ÃfPab)
^Vv^b7
(E+e c e )
c q ry
(4.28)
+â€˜s 7
p,a,b
7
q,c,d
Vd^W
(4.29)
(E+VVebJ
TF+VVcdT
Now substituting Eq. (4.18) for the effective interaction matrices, we
find that the delta functions in Eq. (4.18) further restrict the spin
97
orbital summations in such a way that only expressions (4.26) and (4.29)
are nonvanishing. After some simplification, the nonvanishing contribuÂ¬
tions are found to be
h 2
a
h 2
a ,b
!S 2
a ,b
!i 2
a,b,c,d
+*2 2
a,b,c
+!2 2
a,b,c
PÂ»q,r,s
(E+WV
he+vvs)
2
P.q>r
arxpr[  jb>
TeTVVV
,^E+EbVer)
2
P>q>r
larxrql jb>
v a p q;
E+c.e c )
b r q'
2
P
1 dcxcd  ,ip>
(E+epecÂ£d)
2
p>q
2
p,q
< i p 1 abxqa
1 pcxcb   jq>
(E+e e,e )
q b c
(4.30)
(4.31)
(4.32)
(4.33)
(4.34)
(4.35)
Expressions (4.30) and (4.33) are the only ones which represent valid
thirdorder diagrams as written, however, by interchanging dummy indices
p +* q in Eq. (4.32) and a Â« b in Eq. (4.35), expressions (4.31) and
(4.32) can be combined to yield
2
a, b
7
P,q,r
(E+Vrp'ViE+Cb:VÂ£rT
(4.36)
and expressions (4.34) and (4.35) combined to yield
2
a,b,c
2
p>q
(4.37)
98
which now correspond to the two thirdorder ring diagrams as indicated.
These results verify that one of our original objectives, which was to
include all thirdorder ring and ladder diagrams, has been achieved.
The diagrammatic analysis of the third term in Eq. (4.23) proceeds
in the same way as that of the first two terms, but since the effective
interaction matrix appears twice, it involves considerably more algebra.
For this reason, we simply display the resulting diagrams in Fig. 4.
It is significant to realize that the fourthorder diagrams in Fig. 4
include not only ring and ladder diagrams but also mixed diagrams which
consist of both ring and ladder parts. In the thirdorder analysis, the
first term in Eq. (4.18) was responsible for yielding the ladder diagrams
while the second and third terms yielded the ring diagrams. If we thereÂ¬
fore denote the first term as a ladder part and the second and third
terms as ring parts, the mixed diagrams in fourth order are found to
arise from the product of a ladder part and a ring part. Inducing the
results of the fourthorder analysis to higher orders, we conclude that
our second objective, which was to sum all ring and ladder diagrams in
all orders of the selfenergy expansion, has been exceeded: not only
are all ring and ladder diagrams included in the 2ph TDA selfenergy,
but also the mixed diagrams which exhibit both ring and ladder parts.
The diagonal 2ph TDA selfenergy may also be analyzed diagrammat
ically. This analysis is even simpler than for the full 2ph TDA since
the denominator shifts are now scalars rather than matrices. A comparÂ¬
ison of diagrams obtained with the denominator shift in Eq. (4.30)
versus that in Eq. (4.22) will reveal the significance of the factor of
â€™s. Considering the approximation in Eqs. (4.21) and (4.22) first, we
obtain the following thirdorder expressions
Figure 4
Fourthorder selfenergy diagrams arising from
the 2ph TDA.
100
k Â£
a,p,q
(4.38)
(E+Â£aÂ£peq)2
i, j,
4a,b,p (EWeb)2
u 9
a,p,q (E+e e e )
a p q
j. < ip I abxpb pbxab jp>
a,b,p (E+epEaeb)2
(4.39)
(4.40)
(4.41)
The only difference between these diagrams and the thirdorder diagrams
of the full, 2ph TDA is that the incoming lines on the middle interacÂ¬
tion line have the same labels as the outgoing lines. Now analyzing the
approximation in Eqs. (4.19) and (4.20), we obtain the following third
order expressions
(4.42)
' 2 L
a,p,q
(E+eaâ€™ep~Eq)2
â€” '2 1
a ,b,p
(E+Ep'Â£aÂ£b)2
(4.43)
T,
a,p,q
(4.44)
(E+Â£aEpEq)2
+ T,
, 7
(4.45)
102
Expressions (4.44) and (4.45) are identical to (4.40) and (4.41) respecÂ¬
tively; however, expressions (4.42) and (4.43) both differ by a factor
of from (4.38) and (4.39). This discrepancy is a direct consequence
of the missing factor in the denominator shift, Eq. (4.20), and leads
to an overcounting of these thirdorder diagonal ladder diagrams since,
by rule 5 in Table 5, there should be a factor of % for each pair of
equivalent lines. Similarly, it is rather easy to show that this approxÂ¬
imation overcounts all higher order diagrams containing this diagonal
ladder part.
4.4 Computational Applications and Evaluation of the Diagonal 2ph TDA
SelfEnergy
The main attraction of the diagonal 2ph TDA selfenergy for the
calculation of ionization energies and electron affinities is its pseudo
secondorder structure. Computational experience with the diagram conÂ¬
serving decouplings has taught us that the secondorder selfenergy
approximant is both easily constructed and evaluated. The diagonal 2ph
TDA requires only the additional evaluation of Eq. (4.22) which merely
shifts the secondorder selfenergy poles. As a consequence the diagonal
2ph TDA selfenergy mimics the exact selfenergy by possessing only
simple poles. Another consequence of the pseudo secondorder structure
is that the energy dependence will have the simple analytic form illusÂ¬
trated in Fig. 2. This property, which was absent in the [1,1] self
energy approximant, simplifies the pole search for the electron propagaÂ¬
tor. In spite of the simplicity of its pseudo secondorder structure,
however, the diagonal 2ph TDA selfenergy incorporates diagonal ring,
ladder, and mixed diagrams in all orders of the selfenergy expansion
103
as was shown in the previous section. This property encourages speculaÂ¬
tion that significantly more relaxation and correlation will be accounted
for in this decoupling than with the secondorder decoupling, but the
accuracy of the corresponding ionization energies or electron affinities
can only be evaluated via actual computational applications. As was
done with the diagram conserving decouplings, ionization energies for
the water and formaldehyde molecules were computed using the diagonal
2ph TDA selfenergy.
For the water molecule, calculations were performed with both the
14 and 26 contracted Gaussian orbital basis sets described in Chapters
2 and 3. The principal ionization energies computed with the diagonal
2ph TDA and the diagonal 2ph TDA plus thirdorder constant energy
(2) (2) (3)
diagrams (denoted 2(E)SHIFT and 2(E)SHjpT + Â£(Â»)) are presented in Table
11. Comparing these results with those in Table 9 for the secondorder
and [1,1] selfenergy approximants reveals that each ionization has been
shifted to higher energy. This shift has led to a significant improveÂ¬
ment in the valence ionization energies (3ap, lb^, and lb,,) which are
now within approximately 0.5 eV of the experimental results. For the
inner valence (2a^) and core (la^) ionizations, however, this energy
shift leads to worse agreement. In addition, the diagonal 2ph TDA
results for the la^ ionization exhibit an enormous basis dependence.
The addition of polarization functions in the 26 orbital basis has
yielded nearly 13 eV in additional relaxation. The most probable explanaÂ¬
tion for this basis dependence is that the 2ph TDA selfenergy poles
are determined by HartreeFock orbital energies and twoelectron inteÂ¬
grals rather than by orbital energies alone as with the secondorder
selfenergy. The orbital energies are rather insentitive to basis
changes, whereas the twoelectron integrals are not.
Table 11. Water Results Obtained with the Diagonal 2ph TDA and Diagonal 2ph TDA plus Constant
ThirdOrder SelfEnergies.
14 CGTO's 26 CGTO's
Koopmans
(2)
z(e)shift
(2)
2(E)
(3)
SHIFT +
Koopmans
(2)
z(e)shift
(2:
E(E]
1 (3)
'shift +
Exp.a
lal
559.45
554.91
556.04
559.40
542.78
543.51
540.2
2al
37.04
34.62
35.00
36.52
33.69
33.93
32.2
3a j
15.43
13.61
13.95
15.66
14,12
14.31
14.7
lbl
13.78
11.64
11.97
13.67
11.83
12.00
12.6
lb2
19.52
18.51
18.80
19.34
18.46
18.61
18.6
aSiegbahn al_. (1969).
105
Calculations for the formaldehyde molecule were performed with the
24 and 42 orbital basis sets described in Chapter 3. Ionization energies
were computed using the diagonal 2ph TDA and are presented in Table 12.
The thirdorder constant energy diagrams were not evaluated for this
molecule. Similar to the water results, the diagonal 2ph TDA results
for formaldehyde are also consistently higher in energy than the second
order results (Table 9). This shift considerably improves the valence
ionization energies; however, the average deviation from the experimental
results remains approximately 0.8 eV. The core ionizations within the
diagonal 2ph TDA suffer a small deterioration in accuracy but do not
exhibit the extreme basis dependence which was observed in the water
calculations. Part of the discrepancies between the diagonal 2ph TDA
and the experimental results can certainly be eliminated by further basis
saturation; however, in the next chapter, we will propose that even ioniÂ¬
zation energies of 1.0 eV accuracy are usually sufficient to unambiguously
interpret photoelectron spectra if combined with a calculation of relaÂ¬
tive photoionization intensities.
Cederbaum and coworkers have recently developed computer programs
which implement the full, nondiagonal 2ph TDA to the selfenergy and
have reported several molecular applications (Cederbaum ert aj_., 1977 ,
Schirmer ert aj_., 1977, Cederbaum Â£t aj_., 1978, and Schirmer et al.,
1978). In these calculations, they claim only a 1.0 eV accuracy and rely
heavily on vibrational analyses to assist with the interpretation of
photoelectron spectra. The offdiagonal matrix elements seem to have
little importance in the valence region because the propagator poles are
relatively wellseparated. In the inner valence and core regions, howÂ¬
ever, where principal ionization poles and shakeup poles overlap and
106
Table 12. Formaldehyde Results Obtained with the Diagonal 2ph
TDA SelfEnergy.
24CGT0's 42 CGTO's
Koopmans
(2)
Z^SHIFT
Koopmans
(2)
^^SHIFT
Exp.
lal
560.12
542.86
559.81
542.09
539.43'
2a i
309.09
299.82
308.87
299.31
294.21'
3a j
38.94
35.60
38.18
34.31
34.2b
4a i
23.39
21.75
23.30
21.79
21.15*
5a j
17.29
14.99
17.38
15.33
16.2C
lbl
14.56
14.11
14.45
14.03
14.5C
lb2
19.47
17.91
19.08
17.77
^4
O
n
2b2
12.06
9.92
11.93
10.11
10.9C
aJo 11y and Schaaf (1976).
bHood et al_. (1976).
cEstimated center of gravity (Cederbaum and Domcke, 1977) from
spectrum of Turner et al_. (1970).
d14.38(8) VIP (Turner etal., 1970).
107
interact, level shifts and intensity changes are observed. In these
regions, even the nondiagonal 2ph TDA is not fully satisfactory since
ionstate relaxation and holehole interactions, neither of which are
described by this selfenergy approximation, may also be important
(Wendin, 1979).
CHAPTER 5
PHOTOIONIZATION INTENSITIES
5.1 Introduction
The evaluation of each decoupling approximation in the preceding
chapters was based on the comparison of propagator poles to experimental
ionization energies. This criterion represents a particular bias since
it does not reflect the quality of the FeynmanDyson amplitudes (defined
in Section 1.1). The FeynmanDyson amplitudes determine the spectral
density function (Linderberg and Ohrn, 1973)
A(x,x 1 ;E) =
f r fk(x)f (x
1)6(EEk)
f gk(x)g (x')6(EEk)
E > p
E <
(5.1)
which contains a plethora of useful information. This is evidenced by
the relation of the spectral density to the singleparticle, reduced
density matrix (Linderberg and Ohrn, 1973)
y(x,x') = [ A(x,x';E)dE .
JC
(5.2)
It is important, therefore, to choose a decoupling approximation which
not only yields accurate ionization energies but also an accurate spectral
density. The quality of the spectral density is somewhat more difficult
to evaluate since there are no experimental data with which it can be
directly compared. One evaluation procedure, however, might involve the
calculation of averages of specific oneelectron operators from the
reduced density matrix. The averages could then be compared with
108
109
experimental results. Another procedure which is more closely related
to our goal of interpreting photoelectron spectra is the calculation of
photoionization intensities or crosssections. A theoretical prediction
of relative photoionization intensities can simplify orbital assignments
when the ionization energies are not accurately given. A theoretically
predicted variation in relative intensities with photon energy is parÂ¬
ticularly useful if photoelectron spectra are available with different
ionization sources (Katrib et aj_., 1973). In the following sections,
we derive computational expressions which relate the FeynmanDyson ampliÂ¬
tudes to the total photoionization crosssection, discuss the major
approximations assumed, and then present several applications.
5.2 Derivation of Computational Formulae for the Total Photoionization
CrossSection
The differential crosssection for photoionization derived from
firstorder, timedependent perturbation theory using a semiclassical
model for the interaction of radiation and matter is (Bethe and Salpeter,
1957, Kaplan and Markin, 1968, Smith, 1971)
da 4/L
dfi*
c IA01 w
(5.3)
In this equation, is the momentum of the ejected photoelectron and
it = S it. is the vector potential. For a closedshell system, the initial
k K
state N> can be represented by an antisymmetrized Nelectron wavefunction
[N>  0q(Xj,X2,   â€¢ x^j) (5.4)
and the final state N,s> is represented by
110
;N,s> = (N/2)2 0AS [v(kf,r)a(c)Ã's0(x1,x2, . . . x^)
 v^.i^eteH^Ãx^x^ . . . xN_j)] . (5.5)
Here $ and $ are the two, doublet spin components of the (Nl)
SCt S p
electron ion, v(k^,r) denotes the wavefunction of the photoelectron, and
is an antisymmetrizer
,1 r
11
Â°AS N PkNJ '
(5.6)
The form of Eq. (5.5) guarantees that the singlet spin symmetry of the
system is preserved.
Evaluating the matrix element in Eq. (5.3) yields (Purvis and Ohrn,
1975a)
= SZ J ?s gs(r)dr
+ /? J v (Ãf,r)ps(r)dr (5.7)
where
N'^CrMi) =  (Xj, . . . â€™tN_1)Â«â€™0(x1, . . . xN1 ;r,x(c))
x dXj . . . dxN_j (5.8)
and
N'1;ps(r)x(c) = (Nl)  * (Xj, . . . xN_1)^1 ^ Â®0(xr â€¢ â€¢ â€¢ xN_j;
r,x(c))dXj . . . dxN_1 . (5.9)
The first term in Eq. (5.7) relates the FeynmanDysori amplitude g (r) to
the photoionization crosssection and the second term arises from the
nonorthogonality of v(lt^,r) and i>g. When v(k^,r) is strongly orthogonal
to $q, this term vanishes. Even when v(Â£f,r) is not strongly orthogonal
to $q, the first term in Eq. (5.7) will dominate if the kinetic energy
Ill
of the photoelectron is much greater than its binding energy (Rabalais
et al., 1974), therefore we shall neglect this term:
8^kf)
cAq2Ã¼)
v (kf,r)Ã$ gs(r)dr 2
(5.10)
Further simplification can be obtained by neglecting retardation of
the photoelectron momentum by the photon momentum (also called the dipole
approximation in analogy with photon induced transitions in the discrete
spectrum, Bethe and Salpeter, 1957, Steinfeld, 1974). The vector potenÂ¬
tial t has the following plane wave decomposition in terms of the inciÂ¬
dent photon momentum f and polarization n
Ãs = \tQ\ exp (i^ â– f$)n . (5.11)
If the wavelength of the incident photon is large compared with molecular
dimensions, the exponential may be approximated by the first term in its
Taylor series expansion (unity)
With this approximation, the differential photoionization crosssection
becomes
8tt2 I Ãtf 
CO)
n â€¢ ?2 .
where
(5.13)
?   v (k>f,r)Vgs(r)dr .
(5.14)
112
Owing to the random orientation of molecules in a gaseous sample,
the experimentally observed photoionization intensities in the molecular
reference frame represent an average over all incident photon directions.
Furthermore, if the incident photon beam is unpolarized, we must also
average over photon polarizations. Making the appropriate averages in
Eq. (5.13) (Smith, 1971), we obtain
^r~J Â¿ I {l"iâ€¢ V + â€¢ (515>
Since the two polarization directions, and n^, and the incident photon
direction, kykj, form a righthanded system of axes, we can write
nr = iv V + V V + V V / iV
and Eq. (5.15) becomes
(5.16)
V
dfi
8tt^ I ]
V
C(i)
8"2*fl
CM
8n2  (
Â¿I
If!
iV .
,2 Â±
8tt
Â£ â– Â£\Â¿ / 1< ridn
0) 1 1 to1 w
{1  cos e }dtt
OJ 0)
(5.17)
(5.18)
(5.19)
In order to evaluate ?2, some form for the photoelectron wave
function v(k^.,r) must now be chosen. In principle, the photoelectron
wavefunction could be obtained by the solution of the BetheSal peter
equation for the polarization propagator where the superoperator resolÂ¬
vent has been modified to include the timedependent interaction of the
radiation and matter fields (see e.g. Csanak et al., 1971). A solution
113
of this type would require the use of continuum functions as well as
discrete functions in the molecular basis and is not yet feasible owing
to several formal and practical difficulties. Alternatively, we seek a
simple but accurate, analytic representation of the photoelectron wave
function. For photoionization of atoms or molecules with high (tetraÂ¬
hedral or octahedral) symmetry, the electronic potential of the ion is
nearly spherically symmetric, and v(k^,r) may be asymptotically repreÂ¬
sented by a plane wave plus incoming Coulomb waves (see e.g. Smith, 1971).
For molecules with lower symmetry, distortions of the electronic potenÂ¬
tial enormously complicate the nature of the incoming waves, in this
derivation, the incoming waves are neglected, and the photoelectron is
simply represented by the plane wave part
v(j
(5.20)
The applicability and implications of this approximation will be disÂ¬
cussed in the following section. With this choice, Eq. (5.14) can be
integrated by parts and yields
(5.21)
(5.22)
In our computational scheme, the FeynmanDyson amplitudes in Eq.
(5.22)are represented by a linear combination of HartreeFock orbitals.
The HartreeFock orbitals can be decomposed into contracted Cartesian
Gaussian functions on each atomic center, and each contracted Gaussian
function can be further decomposed into a sum of primitive Gaussian
functions. Ultimately therefore, the FeynmanDyson amplitudes can be
114
represented as some linear combination of primitive Gaussian functions
on each atomic center
z
a, k
o I
ak
^ak
(?  ta) .
(5.23)
Here a represents the sum over atomic centers and k the sum over primitive
Gaussian functions on each center. Transforming the photoelectron posiÂ¬
tion vector to the coordinate frames of each atomic center by the subÂ¬
stitution r' = rft , Eq. (5.22) can be reexpressed in terms of the primÂ¬
itive computational basis
iÂ£f(2TT)'3/2 exp(iicf â€¢ ft )4  exp(ilcf â– ?' )4>ak(1"' )dr' .
(5.24)
The integral over r' in Eq. (5.24) represents a Fourier transform of
each primitive Gaussian function, and formulae for its evaluation are
derived by Kaijser and Smith (1977).
The final step in this derivation is to integrate the differential
photoionization crosssection (Eq. (5.19)) over all photoelectron direcÂ¬
tions in the solid angle di2f
Since Â£ involves a sum over atomic centers, Eq. (5.25) will contain both
one and twocenter contributions (Schweig and Thiel, 1974)
as =
3ao
Â£
a,k
Â£
6,1
'ak
C61 Qak,Bl
(5.26)
The onecenter terms (a=B) have the form
115
where O^fi^) represents the Fourier transform
*ak(fy = (2tt)_3/2 J exp(iiak(r')dr'
and the twocenter terms (a^S) are given by
Qak,81 = I *ak(ltf)*Blil?f)exPfil'f ' ^aB)dnf
(5.27)
(5.28)
(5.29)
where = Sa  Owing to the orthogonality of the spherical harÂ¬
monics which describe the angular dependence of the Fourier transforms,
the onecenter terms are easily evaluated. The twocenter terms are
complicated slightly by the exponential factor, but these too can be
evaluated analytically. For this purpose, it is convenient to define
the integral
1(1, 12m2^
(0f,<Â¡>f)y1^(0f,iÃf)exp(iÂ£f â€¢ ftâ€žR)dilf, (5.30)
where yq m and ^mg are real spherical harmonics (Harris, 1973) repre
the angular
Using the expansion
senting the angular dependence of the transforms and respectively.
exp(Â±itf â€¢ Â¡taS)
m 1
4tt Â£ Â£ (Â±i)
1=0 m=l
JVkfRaf
)klm(ef,f)ylm(0,lf), (5.31)
Eq. (5.30) becomes
m 1 , m.mâ€žm
1(11m1 ,l2m2) = 4ir^ ^Mi) JV kfRap^Cl ^ 2 1 ^In/0â€™^â€™
(5.32)
In these equations, (6^,<Â¡i^) represents the photoelectron direction in
the atomic reference frame, (0,<>) represents the direction of R^gin the
molecular frame, ji(k) is a spherical Bessel function of 1th order,
mjn^m
and Cq ^ q g q is a ClebshGordan coefficient defined by (Harris, 1973)
116
rmln,2m
L1 1 1
12
{ ^Ijm^f ,4>f2m2 (efâ€™'Ãâ€™fÃyimief^dSif
(5.33)
When the distance vector 1$ â€ž coincides with the molecular zaxis, the
Â«6
expansion in Eq. (5.31) simplifies to
exp(Â±iÂ£f â€¢ iÃa6) = (4rr)'s r, (Â±i)1(21+l)*sj1(kfRaB)yio(0f,Â«f) (5.34)
and Eq. (5.30) becomes
IW
i j rn i m o 0
inim1,l?m?) = (4Tr)â€™s6_ 2 (i)1 (21+1) % V,
1 1 Â¿ Â¿ mlm2 1 = 111g 11121
(5.35)
For arbitrary directions of ft . which do not coincide with the molecular
ap
zaxis, a new coordinate frame can be defined in which S . does coincide
Cip
with this axis, and the spherical harmonics can be transformed to this
coordinate frame using the familiar rotation matrices (Schweig and Thiel,
1974). In this way, Eq. (5.35) can be used to evaluate all twoccnter
integrals arising from any molecular geometry.
5.3 Discussion of Approximations
Relatively little work has been devoted to the theoretical calculaÂ¬
tion of photoionization crosssections for molecules compared with that
for atoms (see e.g. Marr, 1967, Steward, 1967, Kelly, 1976). The major
impediments until recently have been a lack of sufficiently accurate
molecular wavefunctions and an absence of accurate analytic representaÂ¬
tions for the photoelectron wavefunction (Kaplan and Markin, 1967,
Schweig and Thiel, 1974). With the development of efficient, molecular
integral and HartreeFock programs, HartreeFock wavefunctions are now
readily available for a large number of molecules. This availability
in turn has stimulated several theoretical calculations of molecular
117
photoionization crosssections using the frozen orbital approximation
(Rabalais et aK, 1974, Dewar et aj_., 1975, Allison and Cavell, 1978,
and Cavell and Allison, 1978).
In the frozen orbital approximation, the Nelectron reference state
is assumed to be the HartreeFock ground state, and the ion states are
constructed by removing the orbital corresponding to the ionized elecÂ¬
tron. This approximation neglects both correlation corrections in the
ground and ion states and ion state relaxation (see Introduction).
Following Purvis and Ohrn (1975a), we have replaced the frozen orbital
approximation by a manyelectron treatment which incorporates both relaxÂ¬
ation and correlation corrections. This treatment derives from the use
of the FeynmanDyson amplitudes obtained from the electron propagator to
compute the photoionization crosssection.
The form of the photoelectron wavefunction is still a major problem
in the calculation of molecular photoionization crosssections and repreÂ¬
sents the most critical step in our derivation. The plane wave approxiÂ¬
mation was chosen for its simplicity, however, it has several serious
limitations. The most serious limitation is its failure to correctly
predict experimentally observed angular distributions (Bethe and Salpeter,
1957, Schweig and Thiel, 1974). This deficiency is not readily apparent
when the differential crosssection is averaged over all photoelectron
directions, and in our computational applications, we compute only spherÂ¬
ically averaged, total crosssections. Lohr (1972) has shown that by
retaining the second term in Eq. (5.7), qualitatively correct angular
distributions may be obtained. The retention of this term is equivalent
to orthogonalizing the plane wave to every oneelectron function in the
Nelectron ground state and is known as the orthogonalized plane wave
(OPW) approximation.
118
A second limitation of the plane wave approximation is the implicit
neglect of electrostatic interactions between the photoelectron and the
molecular ion. For ionization processes near threshold where the photoÂ¬
electron leaves with low kinetic energy, these interactions are espeÂ¬
cially important, and the wavefunction of the photoelectron exhibits a
decrease in wavelength as r 5. The OPW approximation again partially
corrects this deficiency (Lohr, 1972) but exhibits a rather abrupt change
in wavelength which is more characteristic of a shortrange potential
rather than the longrange Coulomb potential. Other representations of
low energy photoelectrons, e.g. multicentered Coulomb wave expansions,
have also been proposed (Tuckwell, 1970). As the kinetic energy of the
photoelectron tends to higher energies, the effect of electrostatic
interactions on the total crosssection becomes negligible, and the plane
wave approximation becomes sufficiently accurate (Bethe and Salpeter,
1957).
Unfortunately, as the plane wave approximation becomes more suitable
at high photoelectron energies, the dipole approximation (Eq. (5.12))
deteriorates and must be reexamined. Bethe and Salpeter (1957) have
shown that when retardation effects are included in the calculation of
differential photoionization crosssections, the lowest order correction
is proportional to (v/c), where v is the photoelectron velocity. This
correction, however, only effects the angular distribution and vanishes
when the differential crosssection is averaged over all incident photon
directions. Higherorder corrections for retardation are proportional
to (v/c) as are the relativistic corrections, but these are usually
negligible even when the photon wavelength is comparable to the molecular
dimensions.
119
5.4 Computational Applications
In this section, the relative photoionization intensities for the
water and acetylene molecules are computed using the FeynmanDyson ampliÂ¬
tudes obtained in the secondorder, diagram conserving and the diagonal
2ph TDA, renormalized decoupling approximations. Comparative calculaÂ¬
tions corresponding to a Mg photon source (tim = 1253.6 eV), for which
the plane wave approximation is expected to be good, and a He (II) photon
source (tioo = 40.81 eV), for which the plane wave approximation is expected
to be less accurate, are presented. These results are further compared
with intensities obtained using the frozen orbital approximation in order
to assess the magnitude of manyelectron correlation and relaxation corÂ¬
rections. For acetylene, the dependence of intensity on photon energy
for the valence orbitals is plotted, and orbital and density difference
plots are presented and discussed.
The relative photoionization intensities for water corresponding to
a Mg photon source and a He (II) photon source are presented in Tables
13 and 14, respectively. These results were obtained using the 26 conÂ¬
tracted Gaussian orbital basis at the experimentally determined equilibÂ¬
rium geometry as described in detail in Section 3.5 and are scaled relaÂ¬
tive to the 3aj intensity. As expected, the relative intensities comÂ¬
puted for the Mg source in Table 13 compare reasonably well with those
obtained experimentally (Rabalais et al_., 1974). The orthogonal i zed
plane wave results of Rabalais et al_. (1975) are presented for compariÂ¬
son and also correctly order the relative intensities. It is interesting
to note the larger discrepancy between the lb^ and lbg intensities in
our calculations compared to Rabalais et al_. than in the 2a^ and 3a^
intensities. Since for Xray photons the OPW corrections are not expected
Table 13. Relative Photoionization Intensities for Water
Excited by Mg
Ka (1253.6 eV).
(2)
(2) (3)
Orbital
F0a
2(E)
SHIFT +
0PWD
Exp.*3
2a i
4.69
1.32
4.15
4.50
3.84
2.57
3a i
1.00
1.00
1.00
1.00
1.00
Ibl
0.62
0.52
0.55
0.18
0.36
lb2
0.40
0.35
0.37
0.09
0.31
aFrozen orbital approximation.
^Rabalais et al. (1974).
121
Table 14. Relative Photoionization Intensities for Water
Excited by He(II) (40.81 eV).
Orbital
F0a
(2)
2(E)
(2) (3)
Z^SHIFT +
0PWb
Exp.C
2a i
0.55
0.28
0.66
0.84
0.26
not observed
3a j
1.00
1.00
1.00
1.00
1.00
lbl
1.16
1.15
1.15
0.86
0.96
lb2
0.83
0.84
0.83
1.10
0.80
aFrozen orbital
approximation.
bRabalais
et a]_
. (1975)
cRabalais
et al
. (1974)
122
to be large, this discrepency is most likely a basis set effect. Our
basis contained dtype polarization functions on oxygen as well as ptype
polarization functions on the hydrogen atoms, whereas Rabalais et al.
have used only a minimal basis.
The relative intensities presented in Table 14 do not correctly
match the experimentally observed intensities. This is not particularly
surprising owing to the inadequacy of the plane wave approximation at
low photon energies. It is more surprising that the OPW results also
yield an incorrect ordering of the intensities. The failure of the OPW
approximation indicates the necessity for a more accurate photoelectron
wavefunction.
Manyelectron correlation and relaxation corrections appear negliÂ¬
gible in the intensities of both Tables 13 and 14. The deviations beÂ¬
tween the frozen orbital approximation (FO), the secondorder selfenergy
(2)
approximation (2(E)), and the diagonal 2ph TDA with thirdorder constant
(2) (3)
energy diagrams (2(E)j^jpj + 2(Â°Â°))> are small and may be attributed to
differences in the photoelectron momenta. Different ionization energies
obtained with different decoupling approximations were used to compute
the photoelectron momenta from the conservation of total energy
lcf = [2(m  I.P.JP (in atomic units) , (5.36)
consequently, deviations in the I.P.'s result in deviations in K,.
The relative intensities computed with the diagonal 2ph TDA plus
thirdorder constant energy diagrams for the Mg K photon source are
plotted in Fig. 5, In order to compare the theoretical spectrum to the
experimental ESCA (Electron Spectroscopy for Chemical Analysis) spectrum
(Siegbahn et al_., 1969) which is sketched as an insert, the ionization
Figure 5. A plot of the theoretical ESCA spectrum for the valence
ionizations of water. The experimental spectrum of
Siegbahn eT (1969) is sketched in the insert.
[T'l; >j(i
il'i,: â€˜NVii1
a1. Mi'ii
m
125
lines were given a Lorentzian width estimated from the experimental
spectrum. The major features of the experimental spectrum are well
reproduced with the exception of the weak peak at approximately 23 eV.
This peak has been attributed to the ionization of 2a^ electrons by the
Mg K(.Â£3 4 satellite in the photon source (Siegbahn et Â¡TL, 1969) and
therefore should not appear in the theoretical spectrum since it was
calculated for a purely monochromatic source.
A recent experimental study of relative photoionization intensities
of acetylene with various photon sources (Cavell and Allison, 1978)
motivated our examination of this molecule. As preparatory steps to the
calculation of relative intensities, HartreeFock and electron propagator
calculations were performed at the experimental equilibrium geometry,
R(CC) = 2.279 a.u. and R(CH) = 2.005 a.u. (Buenker et a_l_. , 1967), with
two different basis sets. The first basis consisted of Huzinaga's 9s,5p
primitive basis for carbon and 4s basis (unsealed) for hydrogen (Huzinaga,
1965) contracted with Dunning's coefficients (Dunning, 1970) to 4s,2p
on carbon and 2s on hydrogen (see Table 8). The complete molecular basis
consisted of 24 contracted Gaussian orbitals and yielded a HartreeFock
total energy of E(HF) = 76.7948 H. The second basis augmented the first
by the addition of a set of d functions on each carbon atom (oj = 0.60)
and a set of p functions on each hydrogen atom (oip = 0.75). The addition
of these diffuse, polarization functions brought the size of the basis
to 42 contracted Gaussian orbitals and yielded a HartreeFock total
energy of E(HF) = 76.8267 H.
Valence ionization energies were computed with the secondorder
selfenergy approximation and the diagonal 2ph TDA for each basis. The
results for the 24 orbital basis are presented in Table 15 and the
126
Table 15.
Valence Ionization
(24 CGTO's).
Energies for Acetylene
Orbital
(2)
Koopmans Z(E)
(2)
^^SHIFT
Ced.a
r 1
Exp.
2ag

38.29

28.9
27.6
Shakeup
2ag
28.42
25.47
26.48
23.9
23.5
3og
18.55
16.39
17.08
16.4
16.8
11.32
11.20
11.29
10.8
11.4
2a
u
20.81
18.21
19.12
18.0
18.7
aCederbaum
et al.
(1978).
^Cavell and
Al 1ison
(1978).
127
results for the 42 orbital basis appear in Table 16. In each table the
full, nondiagonal 2ph TDA results of Cederbaum et aj_. (1978) and the
experimental results of Cavell and Allison (1978) are included for comÂ¬
parison. In contrast to previous calculations reported here, the ioniÂ¬
zation energies obtained for acetylene are larger (with the exception of
the lr ionization) than the experimental values. A breakdown in the
quasiparticle picture for the 2a^ ionization which is evidenced by the
shakeup pole may account for the particularly poor results for these
ionizations. The offdiagonal 2ph TDA contributions considerably imÂ¬
prove these two ionization enerqies (Cederbaum et al_., 1978), however,
the shakeup energy still disagrees with the experimental value by more
than an electron volt.
Relative photoionization intensities for acetylene are presented in
Table 17 for a Mg photon source and in Table 18 for a He (II) photon
source. In both tables the intensities were scaled relative to the In
u
intensity. As for water, the intensities corresponding to the Mg
source are in reasonable agreement with experimental intensities and show
only a slight diminution when correlation and relaxation effects are
included. The OPW result of Cavell and Allison (1978) for the 2a
9
ionization in Table 17 seems unexplicably high, and the 3a ^ ionization
is not observed experimentally (Cavell and Allison, 1978).
The relative intensities for acetylene computed with the lie (II)
source exhibit better agreement with the experimental results than did
the water results. Here, only the relative intensities of the two
weakest ionizations, the 2a^ and 2au, were reversed. The OPW results
predict the correct ordering but attribute a much weaker intensity to
the 2
128
Table lfi. valence Ionization Energies for Acetylene
(42 CGTO's).
Orbital
Koopmans
(2)
1(E)
2 o

38.20
Shakeup
2ag
28.07
24.56
3o
18.46
16.56
3
lTru
11.16
11.13
2a
u
20.90
no results
(2)
(E)
1,1 ;SHIFT
Ced.a
Exp.b
34.45
28.9
27.6
25.73
23.9
23.5
17.24
16.4
16.8
11.24
10.8
11.4
19.47
18.0
18.7
aCederbaum et. aj_. (1978).
bCavell and Allison (1978).
129
Table 17. Relative Photoionization Intensities for Acetylene
Excited by Mg Kq (1253.6 eV).
Orbital
F0a
(2)
E(E)
(2)
z<^e^shift
0PWb
r b
Exp.
2ag
18.90
17.75
17.79
26.86
13.3
3ag
0.51
0.48
0.48
0.65
not observed
1.00
1.00
1.00
1.00
1.00
2a
u
7.74

7.50
11.18
9.5
aFrozen orbital approximation.
bCavell and Allison (1978).
130
Table 18. Relative Photoionization Intensities for Acetylene
Excited by He(II) (40.81 eV).
Orbita 1
F0a
(2)
2(E)
(2)
r^SHIFT
0PWb
c b
Exp.
2a
9
0.45
0.36
0.72
0.09
0.26
3ag
0.45
0.44
0.87
0.77
0.68
llTu
1.00
1.00
1.00
1.00
1.00
2au
0.29

0.50
0.41
0.29
aFrozen orbital approximation.
bCavell and Allison (1978).
131
The relative intensities in Table 17 computed with the diagonal
2ph TDA and Mg source have been plotted in Fig. 6. As in Fig. 5,
the lines have been given Lorentzian widths and the experimental specÂ¬
trum appears as an insert. The most striking difference between these
two spectra is the absence of the intense shakeup in the theoretical
spectrum, hot only was this peak predicted to lie 7 eV higher than the
experimental peak in energy, it also yielded a relative intensity several
orders of magnitude smaller than the 2a^ intensity (at 23.4 eV). The
weak experimental peak at about 14.5 eV arises from the ionization of
2
absent in the theoretical spectrum. The 3a^ peak which should occur at
16.8 eV is not observed experimentally but can be identified as a shoulder
on the 2ou peak in the theoretical spectrum.
Figure 7 shows the dependence of the photoionization crosssection
on photon energy from 020(1 eV. Over this energy range, ionization from
the 1tiu orbital is predicted to be most probable, and this is verified
by the experimental data of Cavell and Allison (1978) obtained at 21.2 eV
(He ( I )), 40.8 eV (He (II)), and 151.4 eV (Zr M.,). Another interesting
feature of this figure is the slight minimum exhibited by the 2a^ curve
around 125 eV. Minima of this type in the photoionization crosssection
were first predicted by Cooper (1962) and are referred to as Cooper
minima. Cooper minima have been observed in the crosssections for a
number of atoms by means of Xray and ultraviolet absorption spectroscopy
(see e.g. Codling, 1976); however, at present there are no photoelectron
spectroscopic data of this type owing to the limitation of photon sources
which are available. The application of synchrotron radiation to photoÂ¬
electron spectroscopy should soon offer a means of studying the energy
Figure 6. A plot of the theoretical ESCA spectrum for the valence
ionizations of acetylene. The experimental spectrum of
Cavell and Allison (1978) is sketched in the insert.
133
Figure 7.
A plot of the
photon energy
in the region
photoionization crosssections versus
for the valence orbitals of acetylene
0200 eV.
136
dependence of the orbital crosssections in a continuous energy interval
from 0 to about 200 eV (Codling, 1973).
Orbital plots are presented for the 2o^, 30^, and Itt^ FeynmanDyson
amplitudes in Fig. 8. Since the dominant component in each of these
amplitudes is the corresponding HartreeFock orbital (i.e. 2g , 3c?g, and
1ttu, respectively) the correlation and relaxation corrections are not
readily observable. In order to examine these manyelectron contributions
more readily, density difference plots between the amplitudes of Fig. 8
and their corresponding HartreeFock orbital were made and are presented
in Figs. 911. Since the FeynmanDyson amplitudes are not normalized,
it is necessary to normalize them before computing the density difference.
In all of these plots, the HartreeFock density is negative which means
any positive distortions imply density enhancement in the propagator amÂ¬
plitude while negative distortions imply density diminution. Figure 9
shows a density diminution in both the CC and CH bonding regions with
a density enhancement in the antibonding region. Figure 10 shows an
enhancement in the CC pi bonding, and Fig. 11 shows an enhancement in
the CH bonding with a slight diminution of antibonding character.
It is apparent from the numerical results presented in this section
that the calculated relative photoionization intensities, at least within
the plane wave approximation, are not very sensitive to improvements in
the FeynmanDyson amplitudes. It is likely that more accurate photoelecÂ¬
tron wavefunctions will improve the sensitivity of these quantities but
not at the orthogonalized plane wave level. Although the OPW approximaÂ¬
tion is more .justifiable formally, the results of Rabalais et_ aj_. ( 1974)
and Cavell and Allison (1978) do not exhibit any marked improvements to
the plane wave results in the two molecules studied here. Qualitative
Figure 8. Orbital plots for the 2au, 3c>g, and liru FeynmanDyson
amplitudes of acetylene. Approximate scales on the
plots are Â±0.72 a.u. for 2Ã¼u, Â±0.28 a.u. for 3aa, and
Â±0.29 a.u. for liru.
138
Figure 9. A density difference plot between the 3oq FeynmanDyson
amplitude and the 3og HartreeFock orbital of acetylene.
Each contour represents a density increment of about
1.5 x 104 atomic units.
140
m
i.i
\  *
t'k
Ikik
H&a
â€™x'V'&iif
â– :U*'
/Â§v
c'iAlr:.' ' â€™
VMv ip^si P
'tÂ» "if
Figure 10. A density difference plot between the liru FeynmanDyson
amplitude and the 1ttu HartreeFock orbital of acetylene.
Each contour represents a density increment of about
3.6 x 10~3 atomic units.
1
zn
Figure 11. A density difference plot between the 2au FeynmanDyson
amplitude and the 2au HartreeFock orbital of acetylene.
Each contour represents a density increment of about
4.4 x 10~4 atomic units.
HO D H
145
changes in the propagator amplitudes may be observed via density
difference plots; however, more quantitative comparisons must await
other computational applications.
CONCLUSIONS AND EXTENSIONS
The primary objective of this dissertation has been the investigaÂ¬
tion of alternative decouplings which may allow more accurate calculaÂ¬
tions of the electron propagator with less computational expenditure
than the operator product decouplings. Contrary to previous contentions
that decouplings of the propagator equations of motion and choice of
reference state averages are interrelated (Oddershede and Jorgensen,
1977), we have demonstrated that they are in fact independent approximaÂ¬
tions, hence justifying the systematic improvement of one or the other.
We have chosen to examine the decoupling approximation while maintaining
a HartreeFock reference state since these reference states have a simple
theoretical representation and are efficiently generated by a number of
available computer programs. Since HartreeFock reference states ignore
electron correlation (by definition), the correlation and relaxation
corrections to the electron propagator must be exclusively incorporated
through the decoupling approximation.
The superoperator formalism which was originally introduced as a
notational convenience has a rich algebraic structure that has been
largely unexploited. As we demonstrate in Chapter 3, the superoperator
Hamiltonian may be partitioned and a perturbation theory may be developed
in the operator space. The perturbation expansion of the superoperator
resolvent readily yields the Dyson equation and allows an identification
of wave and reaction superoperators. The definition of these auxiliary
superoperators elucidates the parallelism between the superoperator
146
147
formalism and the diagrammatic expansion method and permits a unified
conceptual framework.
When the superoperator resolvent is approximated by an inner proÂ¬
jection, various choices of the inner projection manifold which correÂ¬
spond to different decouplings may be viewed as different approximations
to the wave superoperator in the following equation (Born and Ohrn, 1977)
[h) = W(E) aj .
In particular, the operator product decoupling corresponds to the choice
W = I + {a^a^}
W = Ã + {a^} + ia+aja^}
* * * )
the moment conserving decoupling corresponds to
W = Ã + H
W = Ã + H + H2
... 5
and the diagram conserving decoupling of the Pade' type corresponds to
W(E) = 1 + T0(E)V
W(E) = I + T0(E)VT0(E)V
Alternatively, the superoperator resolvent may be approximated by
an outer projection in which case decoupling approximations correspond
to various approximations to the reaction (or selfenergy) superoperator.
Two decouplings of this type which have been presented are the simple
148
diagram conserving decoupling which corresponds to
t(E) = V + VT(E) V
wi th
T(E) = T0(E)
T(E) = T0(E) + T0(E)VT0(E)
and the renormalized decoupling which corresponds to
t(E) = V + VT(E)V
with
T(E) = Tq(E) + T0(E)MT(E)
and where M is an effective potential obtained by projecting the perturÂ¬
bation superoperator onto a subspace of the full operator space. When
this subspace is chosen to consist of simple operator products, this
decoupling is formally equivalent to the operator product decoupling,
however, other approximations may be envisioned.
Each of the abovementioned decoupling approximations has been
studied in detail computationally and has been evaluated on the basis
of a comparison of the propagator poles with experimental ionization
energies. The moment conserving decouplings discussed in Chapter 2 were
found to yield unacceptably poor numerical results. It is plausible
that since the moment matrices involve averages of various powers of
the full superoperator Hamiltonian, the moment expansion may not be conÂ¬
vergent anywhere in the complex plane. If this is the case, analytic
continuation cannot be performed, and the Pade1 approximant method will
not improve convergence. In any event, it may be generally concluded
149
that the number of conserved moments is no indication of the accuracy
of this decoupling.
A more successful approach is the diagram conserving decouplings
presented in Chapter 3. A simple secondorder selfenergy truncation
was found to overcorrect the HartreeFock orbital energy estimate of
the ionization energy and yielded results that were generally too small.
The inclusion of thirdorder diagrams, particularly rings and ladders,is
necessary to obtain reasonable agreement with experiment. Energy shifts
resulting from the addition of diffuse, polarization functions to the
computational basis were observed to be nearly the same magnitude as the
shifts obtained from thirdorder ring and ladder diagrams, hence care
must be taken to insure basis saturation.
The renormalized 2ph TDA decoupling was found to be the most satisÂ¬
factory decoupling studied. Although only a diagonal approximation was
implemented, principal ionization energies accurate to within an electron
volt were generally obtained. Shakeup energies were less accurately
described owing to the neglect of ion state relaxation and holehole
correlation in the 2ph TDA.
The numerical results of Chapter 5 indicate that the calculation of
relative photoionization intensities using a plane wave approximation
for the photoelectron is not a sensitive reflection of the quality of
the FeynmanDyson amplitudes. The severity of the plane wave approximaÂ¬
tion most likely obliterates the manyelectron relaxation and correlation
corrections. More accurate photoelectron representations such as orthog
onalized plane waves or Coulomb waves should afford a more accurate evalÂ¬
uation of the FeynmanDyson amplitudes, or alternatively, oneelectron
properties may be computed using the singleparticle, reduced density
matrix obtained from the spectral density.
150
There are several other possible extensions of this research that
should be mentioned. Most obvious is an application of these decoupling
ideas to the polarization propagator. Currently, decouplings of this
propagator are based on operator products. The odd or Fermionlike
operator products in the electron propagator theory must be replaced by
even or Bosonlike operator products to account for the particleconservÂ¬
ing nature of the polarization propagator theory, but with this minor
modification, the application of the diagram conserving and renormalized
decouplings also seems feasible.
More accurate decouplings are still needed for the electron propaÂ¬
gator in order to better describe shakeup energies. With this goal in
mind, Cederbaum and coworkers have recently implemented the full, nonÂ¬
diagonal 2ph TDA by exploiting spin and orbital symmetry to reduce the
dimension of the 2ph operator subspace. Preliminary results indicate
that this approximation provides some improvement but is still inadequate.
The need for higher operator products, particularly quintuple products,
has been demonstrated by Bagus and Viinikka (1977) in Cl studies and may
necessitate the implementation of a 3p2h renormalization.
Another attractive avenue for extension is a combination of the
2ph TDA renormalization with the Hamiltonian partitioning of Kurtz and
Ohrn. Some work along this line is in progress by Ortiz and Ohrn (1979)
and preliminary results are very encouraging. One criticism of this
approximation is that if one interprets it as an alternative partitioning
of the superoperator Hamiltonian in which the eigenvalues of the unperÂ¬
turbed Hamiltonian correspond to AE(SCF) energies (as in Eq. (3.23)),
then the denominators of the selfenergy which are obtained from
(EIHq)"* (some operator product)
151
should also involve AE(SCF) energies. For example, the secondorder
selfenerqy denominator would be
[E + AEjJSCF)  AE] (SCF)  AE^SCF)]'1 .
Since it is probably not feasible to perform AE(SCF) calculations for
all doublet (Nl) and (N+l)electron states, simple approximations to
the AE(SCF) energies such as
ak>
or
might be considered.
Finally, the problem of open shell reference states has not been
considered here. The retention of a restricted HartreeFock reference
state in this regard will be a tremendous simplification, and the various
spin and orbital symmetry couplings between the reference state and
operator manifold may be treated with the theory of tensor operators.
2SI
(M11j>I\ A)
3  4 ? ~ * * + "â€™)(** â€œ * i ~ â€œ * * '/) râ€”T *y A A '~7
(hi  .(Â»)
wn \i>
â€” if) ryryAAry ^
~ >  q 3 t "a)(q3  *^3 ~ J > + ;./) râ€”7 ry A A ^7 b
<(v"l \â€˜1â€˜) A AAA A [
d > 3 ? t ** >) ( >  *A+.7) ry A A A A? b
<Â«i j<*<7iT".'> A.xAAA i
S qifA)( 3 ~  J3 ~* //) f7 r~7 A A A b
aaaaa Â¡
â€˜ * 3  > + .7H* * â€œ * + 7) A A ry A A
{â€™If 11{
>  B3~tA+.7)(Â°3a?ti3Ay) A A y7 A A
ATT^/Tk^111 >jA kAAA A
A'iii'iX/''i i'iÂ«kâ€˜i"''nA) kAAAA [
(Â°5<5'>+3)(â€˜3i5"5+,V)
_("7 I 1 ar.i >"<s:.rI j Â¿./)(ZÃœ?"' u I) '
XI0N3ddV
OÃ)
(O)
un
(O)
(N)
(S 53)(_2 ^ 'Ã*+_*) r~^ r~7 r~7 r~J r~~r Z
(,,J\  Ãlâ€œ)(flU I  *â– / ><*/?"  ><> ^ I
3 q> + P>) ,~~7r~7rY^7^V ^
<<7w I biiypil I  r/.Ã^ (/'  Â«/) AAAiA i
(* * ~ C3)( 5 râ€”7 rj !â€”y r~7 Z
AA.Ã.ÃA [
(SMf'ÃS'jtS) e
(B 3  J J  J5 + a>)(Â° ^ v 3  fl s â– * L>) A^r7r7r7'r'~7 ^
Ã±^Ãw"ÃÃ''.o AAAAA [
UM)
+ ">)(*2_ Ã3*5â€™_"3)
(1/1/1 \b'l)(l, 11 Ã/Â«)(i/l â€¢ I:
33333: f
(1)
(I/ 
â€¢ .Y)
33333
Â£91
APPENDIX 2
â– 45Ãj
+ ;; (e. + e, + e + e ,  e )
r,s,s' J
x [,Ã + !5]
+ X <1111  jr>
r,s,s1,1,1 '
x [h  !j + J5  *5
+   h
r s 1 1 1 s s s
+ !s]
+ X <1 's   jl>
r,s ,s 1,1,11
x [2 +  
1 s rls rl rl
+ +2 
rl s r 1 1 ' Is r
  2 
+ ]
+ X <1111 rs'>
r,s,s1,1,11
x [h  