Electron-initiated chemistry in polyatomic
molecules
by
S. Tonzani
Laurea in Chemistry, Universita’ di Perugia, Italy, 2001

A thesis submitted to the
Faculty of the Graduate School of the
University of Colorado in partial fulfillmentof the requirements for the degree ofDoctor of Philosophy
Chemical Physics
2006

This thesis entitled:
Electron-initiated chemistry in polyatomic molecules
written by S. Tonzani
has been approved for the program in Chemical Physics

Chris H. Greene

Robert Parson

Date

The final copy of this thesis has been examined by the signatories, and we find thatboth the content and the form meet acceptable presentation standards of scholarly
work in the above mentioned discipline.

iii
Tonzani, S. (Ph.D., Chemical Physics )Electron-initiated chemistry in polyatomic molecules
Thesis directed by Prof. Chris H. Greene

We develop techniques to describe theoretically the process of electron collision
with a polyatomic molecule and the interplay of many degrees of freedom in the evolution of the complex. In particular, we characterize resonant processes involving large
molecules, believed to be important in radiation damage to nucleic acids. We also
examine the electronic resonances in the DNA and RNA bases, and in the backbone
constituents. Dissociative recombination reactions, in which a polyatomic ion recombines with an electron and quickly dissociates, are also considered. We elaborate aframework to describe the reaction mechanism and the role of the Rydberg states in
the recombination, show results for HCO+ and compare them with experimental datain order to resolve some of the existing controversies for this reaction.

iv

Acknowledgements

The research work presented here is the result of a few years of my time, but alsoof the collaboration and efforts of many other people. I would like to acknowledge firstof all my advisor Chris Greene, whom besides being a great scientist is a great person,and has been able to deal with me also at the beginning when I produced graphs very
sparingly. He is also one of the few physicists I know that is not repelled by moleculesof more than two atoms. Thanks also to Steve Leone for pointing me, when I arrived
at JILA and I was TAing for him, towards Chris’ group.
Throughout these years I have collaborated successfully with many people. I
would like to thank in particular Robin Santra, from whom I learned a great deal about
quantum chemistry and a lot of other stuff during his stay at JILA. Also my collaboratorsin the dissociative recombination project: Åsa Larson, Viatcheslav Kokoouline andRoman Čurik. I would also like to thank Åsa for my visit to Stockholm and Slava formy stay in Orlando. I would like to thank Jonathan Tennyson and his group in London
for my visit there to learn the R-matrix code, and in general for useful discussions. Also
I acknowledge the collaboration and discussions with Leon Sanche and Laurent Caronat the University of Sherbrooke and discussions with Paul Burrow. An advantage ofworking in the Greene group is the variety of science that I came in contact with, and I
would like to thank all the other group members that have been here during my stay. A
special thanks goes to some of my best friends during these years: Daniele Bortolotti and
Josh Dunn, that besides being great for partying with, have also helped me innumerable

v
times with my work and to understand science beyond my area of expertise. Thanks
also to Alessandro Volpi who showed me where JILA and Boulder were on the worldmap. Thanks also to all the people at JILA and in particular Peter Ruprecht without
whom the computational work would have been a lot slower.
All of this would not have been possible without my family, it is hard to be half aworld away from them. And thanks to Jennie who, among all the very good stuff theseyears in Colorado have brought me, is the very best.
The work here described was performed with financial support from Departmentof Energy, Office of Science, National Science Foundation, the Keck Foundation and
NERSC through allocation of computing resources.

vi

Contents

Chapter
1 Introduction

1

2 Methods: Electron scattering

6

2.1

Scattering theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

2.2

Approximations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

8

2.2.1

R-matrix method . . . . . . . . . . . . . . . . . . . . . . . . . . .

11

2.2.2

Finite elements . . . . . . . . . . . . . . . . . . . . . . . . . . . .

13

2.2.3

Potential and DFT . . . . . . . . . . . . . . . . . . . . . . . . . .

15

2.3

Photoionization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

21

2.4

Other methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

25

2.4.1

Ab initio methods . . . . . . . . . . . . . . . . . . . . . . . . . .

25

2.4.2

Rmatrix-UK . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

26

QDT and long range interactions . . . . . . . . . . . . . . . . . . . . . .

28

2.5.1

External potentials . . . . . . . . . . . . . . . . . . . . . . . . . .

31

2.5.2

Coulomb potential . . . . . . . . . . . . . . . . . . . . . . . . . .

33

2.5.3

Dipole physics . . . . . . . . . . . . . . . . . . . . . . . . . . . .

34

2.5.4

Perturbation theory . . . . . . . . . . . . . . . . . . . . . . . . .

36

2.5

3 Methods: vibrations, rotations and dissociative recombination
3.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

39
39

vii
3.2

Ab initio surfaces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

42

3.3

Frame transformation and adiabatic methods . . . . . . . . . . . . . . .

47

3.4

Siegert states . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

49

3.5

Model Hamiltonian including the Renner-Teller coupling . . . . . . . . .

52

3.6

Reaction matrix and O’Malley formalism for the DR cross section . . .

54

4 Results: electron scattering, small molecules

63

4.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

63

4.2

Computational Details . . . . . . . . . . . . . . . . . . . . . . . . . . . .

63

4.3

Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

65

4.3.1

Quantum defect calculations . . . . . . . . . . . . . . . . . . . .

69

4.3.2

XeF6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

75

Photoionization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

78

4.4.1

CO2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

78

4.4.2

SF6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

81

HCO+ Rmatrix-UK calculations . . . . . . . . . . . . . . . . . . . . . .

83

4.5.1

86

4.4

4.5

Perturbation results . . . . . . . . . . . . . . . . . . . . . . . . .

5 Results: electron scattering, DNA bases
5.1

90

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

90

5.1.1

Elastic cross sections . . . . . . . . . . . . . . . . . . . . . . . . .

93

5.1.2

Positions and widths of resonances . . . . . . . . . . . . . . . . .

96

5.1.3

Resonance molecular structures . . . . . . . . . . . . . . . . . . . 103

5.1.4

Uracil - Rmatrix UK . . . . . . . . . . . . . . . . . . . . . . . . . 112

5.2

DNA backbone . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115

5.3

Results: THF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119

5.4

Results: H3 PO4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129

5.5

Outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130

viii
6 Results: dissociative recombination

132

6.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132

6.2

Model 1: no Renner-Teller couplings . . . . . . . . . . . . . . . . . . . . 134

6.3

6.2.1

Ab initio surfaces . . . . . . . . . . . . . . . . . . . . . . . . . . . 135

6.2.2

Vibrational dynamics . . . . . . . . . . . . . . . . . . . . . . . . 137

6.2.3

Autoionization widths . . . . . . . . . . . . . . . . . . . . . . . . 139

6.2.4

Electron scattering and cross section . . . . . . . . . . . . . . . . 140

Model 2: inclusion of Renner-Teller couplings . . . . . . . . . . . . . . . 141
6.3.1

Estimation of the effect of autoionization and CO vibration on
the cross section . . . . . . . . . . . . . . . . . . . . . . . . . . . 147

6.3.2
6.4

Comparison with our previous theoretical study . . . . . . . . . . 152

Outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155

7 Summary

156

Bibliography

158

Appendix
A Finite element matrices

168

B Related publications

171

ix

TablesTable
3.1

Comparison of adiabatic vibrational energies to 3D calculations for Model
2 of HCO+ . From Ref. [130]. . . . . . . . . . . . . . . . . . . . . . . . .

5.1

59

Energies, widths and largest partial waves (and percentages) of the resonances discussed in the text for the DNA and RNA bases. From Ref.
[194]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5.2

98

Energies, widths and largest partial waves of the resonances discussed in
the text for the DNA backbone components. From Ref. [195]. . . . . . . 129

6.1

Comparison of vibrational frequencies for HCO+ using an adiabatic approximation and 3D calculations for Model 1. . . . . . . . . . . . . . . . 139

6.2

Comparison of calculated and experimental autoionization widths for
HCO+ . From Ref. [104]. . . . . . . . . . . . . . . . . . . . . . . . . . . . 140

6.3

Landau-Zener parameters and probabilities for typical avoided crossingsfor HCO+ . From Ref. [130]. . . . . . . . . . . . . . . . . . . . . . . . . . 149

xFiguresFigure2.1

Finite element grid for CO2 . From Ref. [192]. . . . . . . . . . . . . . . .

16

2.2

Electron-molecule interaction potential for N2 . From Ref. [192]. . . . . .

22

3.1

Schematic representation of direct and indirect DR processes. . . . . . .

43

3.2

HCO+ potential energy surfaces as functions of CH stretching and bending. From Ref. [104]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

45

3.3

Jacobi coordinates and adiabatic curves for HCO+ vibrational states. . .

50

3.4

Eigenquantum defects of the electronic states of HCO+ we include in ourstudy. From Ref. [130]. . . . . . . . . . . . . . . . . . . . . . . . . . . .

55

3.5

Adiabatic potential curves for Model 2 of HCO+ . From Ref. [130]. . . .

57

4.1

Structure of the finite element matrix Γ. From Ref. [192] . . . . . . . .

66

4.2

e-N2 integrated elastic cross section. From Ref. [192]. . . . . . . . . . .

70

4.3

e-CO2 integrated elastic cross section. From Ref. [192].

71

4.4

e-CO2 integrated elastic cross section using a potential with no adjustable. . . . . . . . .

parameters. From Ref. [191]. . . . . . . . . . . . . . . . . . . . . . . . .

72

4.5

e-CO2 differential elastic cross section. From Ref. [191]. . . . . . . . . .

73

4.6

e-C2 H4 integrated elastic cross section. From Ref. [192]. . . . . . . . . .

74

4.7

Eigenquantum defects for HeH at 0.1 eV. From Ref. [192]. . . . . . . . .

76

xi
4.8

The integrated elastic cross section for e-XeF6 shown as a function ofenergy.

4.9. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

79

Resonance wavefunction map for e-XeF6 , for the resonance at 2.5eV.
Plotted is the largest eigenvalue of the time-delay matrix. At the center
it is visible the Xe atom, whereas the four fluorines in the plane arelocated at 90 degrees with respect to each other, between the yellow andblue areas, since the wavefunction has a node at the nuclei. . . . . . . .

80

4.10 CO2 photoionization cross section for the three highest nondegenerate
initial orbitals. From Ref. [191]. . . . . . . . . . . . . . . . . . . . . . . .

82

4.11 e-SF6 integrated elastic cross section. . . . . . . . . . . . . . . . . . . . .

84

4.12 SF6 photoionization cross section for the 7 highest nondegenerate initialorbitals. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

85

4.13 R-matrix UK calculations of eigenquantum defects of HCO at 0.01eV. .

87

4.14 Comparison of perturbation theory and R-matrix propagation for the
eigenquantum defects of HCO with multipole potentials outside the Rmatrix box.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

89

5.1

Schematics of the DNA damage from low energy electrons.

. . . . . . .

92

5.2

Ground state equilibrium structure of the DNA bases. From Ref. [194].

94

5.3

Resonant wavefunction and grid for uracil at 2.4 eV. From Ref. [194]. .

95

5.4

Partial elastic cross section for the 5 DNA and RNA bases. From Ref.
[194]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5.5

Total time-delay for the DNA bases described in the text. From Ref.
[194].

5.6

99. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100

Comparison of calculated time-delay with experimental data for uracil.
From Ref. [194]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104

xii
5.7

Comparison of calculated time-delay with experimental data for adenine.
From Ref. [194]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105

5.8

Comparison of calculated time-delay with experimental data for guanine.
From Ref. [194]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106

5.9

Dominant time-delay eigenvector at 0.2 eV for uracil and first resonancewavefunction at 1.3 eV. From Ref. [194]. . . . . . . . . . . . . . . . . . . 110

5.10 Uracil virtual orbital associated to the first resonance. From Ref. [194].

111

5.11 First and second resonance wavefunctions for adenine. From Ref. [194].

113

5.12 Third resonance wavefunction for adenine. From Ref. [194]. . . . . . . . 114
5.13 Eigenphase sums for e-uracil scattering calculated with the R-matrix UK
code. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116
5.14 Energies of the ground state of the uracil molecule and its first few resonant anionic states. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
5.15 Structures of tetrahydrofuran, phosphoric acid, cyclopentane and deoxyribose. From Ref. [195]. . . . . . . . . . . . . . . . . . . . . . . . . . 120
5.16 Integrated elastic cross section for e-THF and time-delay analysis, comparison to experimental data. From Ref. [195]. . . . . . . . . . . . . . . 123
5.17 Fixed nuclei partial integrated elastic cross sections for e-THF and ecyclopentane. From Ref. [195]. . . . . . . . . . . . . . . . . . . . . . . . 126
5.18 Resonance wavefunctions for e-THF and comparison to virtual orbitalsof THF. From Ref. [195]. . . . . . . . . . . . . . . . . . . . . . . . . . . 127
5.19 e-H3 PO4 partial cross section and time-delay analysis. From Ref. [195].
6.1

131

Potential energy curves and quantum defects for HCO+ from our ab initio
study. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138

6.2

Dissociative recombination cross section of HCO+ obtained for Model 1.
From Ref. [104]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142

xiii
6.3

Dissociative recombination cross section obtained for Model 2. From Ref.
[130]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143

6.4

Reaction rate for DR of HCO+ obtained with Model 2, compared toexperiment. From Ref. [130]. . . . . . . . . . . . . . . . . . . . . . . . . 146

6.5

Resonance potential curves for HCO and vibrational thresholds. From
Ref. [130]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148

6.6

Detail of the spectrum of Fig. 6.5. From Ref. [130]. . . . . . . . . . . . 151

A.1 Top: finite element basis functions from Eqs. A.1-A.4, the four polynomials are plotted in the same order as in the text. Bottom: nodal structure
for one cubic sector in the rescaled variables. . . . . . . . . . . . . . . . 169

Chapter 1

Introduction

It is an exciting time to work in the field of electronic collisions. This area of
research has undergone a major change over the last ten years, mainly due to experimental breakthroughs that have opened the field to the exploration of complex chemical
processes, both in the gas phase and in condensed phases. Electron-driven processes
are in fact fundamental in many areas, from radiation damage to living tissue [24],to plasma etching and processing of semiconductors, lighting applications (and plasmaTVs), scanning tunneling microscope (STM) [171] induced chemistry, chemistry of theatmosphere and interstellar clouds [125].
As it is possible to see, many of these processes are quite far from an electronic
collision with a simple molecule in the gas phase, since they happen in complex environments like plasmas or condensed phases or biological cell environments. Nevertheless,the electron-molecule interaction drives all of them and its understanding is absolutely
fundamental to developing a correct description of the process.
Traditionally the electron collision community has focused its efforts on an extremely narrow slice of this pie, that is to say gas phase processes involving small
molecules (many diatomics) or atoms. Now, these are still interesting processes, but
the wealth of possibilities that are rapidly opening up call for a broader spectrum of2
investigation, and we will try to undertake this effort from a theory standpoint for a
few of the areas cited above.
For the understanding of complex processes it is, again, fundamental to have thesupport of a theoretical description that can do what experiment often cannot, namely
separate different contributions and find out which is the most important, in other words
figuring out the mechanism. Unfortunately, this is a very hard task, in the sense that
performing calculations involving continua interacting with large molecules is very muchat the cutting-edge both of current theoretical approaches and of computational power.
Moreover, the numerous advances in the very closely related field of quantum chemistrythat deals with bound states, are very slowly being implemented in this field. In otherwords, there is not yet a “continuum Gaussian”, referring to the famous suite [150] of
ab initio quantum chemistry programs for bound states. Some groups, however, areexpending significant effort in this direction, see Refs. [126, 136, 149, 206]. Since there
is no “silver bullet”, in most cases targeted approaches have to be applied to get themain physics of the process, even if they are somewhat coarse grained.
An entirely different set of problems is presented by the interaction of the metastableelectron-molecule complex with its surrounding environment, in particular solvents (seeRef. [55] for a review on electron-initiated chemistry in water) and matrices in condensed
environments and solids. In cluster environments (see for example studies of electron
attachment to CO2 clusters recently carried out [46, 177]) the continuum electron interaction with the solvent could also lead to interesting cage effects for the electrons as
it happens for molecules [197]. This can lead to a vast array of effects (negative ion
stabilization, resonant interactions, fast energy transfers) whose exploration has juststarted.

3
The experimental advances that have made possible the exploration of these new
exciting areas of research are deserving of some comment. In the text we will talkfrequently about storage ring experiments, see Ref. [105] for example, which allow the
study of collisions between electrons and charged molecules in a controlled environment,
and at very high energy resolution. These will be referred to in the dissociative recombination part of this thesis. Also thin film techniques [24] have allowed the study ofprocesses that occur in condensed phase, often coupled with the precision and control
that is possible using an STM (scanning tunneling microscopy) [171] to study processes
at the atomic scale. In dissociative processes, COLTRIMS (cold target recoil ion momentum spectroscopy) and other momentum imaging techniques [103] have provided
a more detailed picture of the Coulomb breakup of a system. The new and rapidly
evolving field of high harmonic generation (HHG) from molecules and more generally
the topic of ultrafast intense laser pulses that interact with molecules, have led to newtechniques (like “orbital imaging” [86]) and the intriguing possibility of following theelectronic dynamics (say, in a chemical reaction) as the nuclei move in a pump-probe
experiment [13]. Molecular electronics is another fascinating subject of topical interest,
where techniques used in electron scattering theory can be applied [38].
Many processes are possible, both during and after the period when an electroninteracts with a molecule, and it is not possible to offer a comprehensive treatment of allof them. We will concentrate mainly on gas phase processes, even though many of ourresults will have relations to condensed phase experimental results. In the gas phase anelectronic collision can produce elastic and inelastic scattering, vibrational or rotational
excitation, and then dissociation (processes of dissociative attachment or recombination,
depending on whether the target is a neutral or an ion), electron-impact ionization (so-

4called e-2e processes), and energy transfers between different degrees of freedom which
are extremely effective because unencumbered by selection rules as in photoabsorption.
For example it is possible through electronic collisions to activate species that are notreactive as singlets but will be highly reactive in a triplet state. It is also possible
to reach many dissociative states, giving these collisions a great energy efficiency intransferring energy to other degrees of freedom of the molecule, such as vibrations, also
considering that the colliding electron is indistinguishable from the target electrons andthus will strongly couple to them.
We will focus here on some of these processes; our final goal is to treat electron
scattering from a general polyatomic molecule and to describe processes that so far
theory has not been able to tackle. Ch. 2 starts with a general introduction to electron
scattering theory and illustrates the approximations contained in our methods [191, 192]
and the other techniques by which we treat electronic dynamics. Then we will apply
this description to many systems in Ch. 4 to test its capabilities and we will show its
predictive power [192]. In particular we will then concentrate on the role of electronicresonances in radiation damage to nucleic acids in Ch. 5. This problem is important inview of its medical applications and it has raised a great deal of interest among boththeorists [17] and experimentalists [2, 66, 152, 164]. Many groups around the world areworking on explaining how low energy electrons below the DNA ionization threshold
can lead to one or more strand breaks in the double helix structure of the nucleic
acid, with obvious nefarious consequences such as mutations. We will try to give some
contribution to this effort. Some of the material in Chapters 2 and 4 has been published
in Ref. [192], whereas the majority of the material in Ch. 5 has been published in Refs.
[194, 195] Another part of this thesis deals with dissociative recombination processes,

5
important in the description of the chemistry of interstellar clouds and in the elaborationof star formation models (see Ref. [125] and references therein). We will formulate ourtheoretical description in Ch. 3 and discuss corresponding results in Ch. 6 for the HCO +ion, which we have chosen for its abundance in interstellar clouds. Recombination ofpolyatomic ions is a very complicated process to treat since it involves a coupling ofmany electronic and vibrational degrees of freedom, and, since ions are involved as target
states, the treatment of an infinite number of neutral molecule Rydberg states. Few suchcalculations have been performed to date, and we will show how our results [104, 130]
compare with experiments in merged beams [144] and storage rings [56]. Rydberg states
are highly excited states of atoms or molecules, in which the external electron roams far
from the core and it looks a lot like a continuum electron, we will examine techniques
to deal with the dynamics of these states such as quantum defect theory (QDT). We
will show how Rydberg states can be important and even determinant in some DR
reactions. We will also talk about symmetry-distortion effects such as Jahn-Teller [9]and Renner-Teller [155], vibronic couplings and conical intersections, and how they canprovide a preferential pathway to fast dissociation. Some of the material in Chapters 3
and 6 has been published in Refs. [104, 130].
We will also consider, whenever possible, the outlook on new processes, the possible continuation of this work and the application of the techniques we developed here tolook at a variety of processes in clusters, strong laser field, and more complex biologicalenvironments.
Atomic units (see Ref. [181]) will be used throughout this thesis, except whereotherwise specified.

Chapter 2

Methods: Electron scattering

This chapter describes all the methods we use to tackle the electronic dynamics of a molecular system plus an external interacting electron. After a brief general
introduction to electron-molecule scattering theory, Sec. 2.2 describes an electron scattering approach that will be applied to large molecular systems like the DNA bases in
Ch. 5 and all the approximations that are contained in this framework. The following
sections describe the methods we have employed for dealing with small molecules, also
taking into account their vibrational degrees of freedom, like ab initio quantum chemistry calculations (Sec. 2.4.1) for Rydberg states and the R-matrix UK code of Sec.
2.4.2. Finally a brief account of quantum defect theory (QDT) ideas that have beenused to treat the long-range dynamics of the electron-molecule interaction problem will
be shown in Sec. 2.5, and the treatment of the electrostatic multipolar interaction atlong range by various methods will be described.

7

2.1

Scattering theory
The interaction of an electron with a polyatomic molecule is modeled by the full

Hamiltonian of the compound system:
Ĥ = −

~2
X
X
X
Zα Zβ
1
1X∇
Zα
1X~2
Rα
+
+
∇r i −−. (2.1)~~~2 i
2 α Mα|~r−ri − Rα | α>β | Rα − Rβ | j>l j ~rl |i,α | ~

This operator contains both the nuclear and electronic degrees of freedom, indicatedrespectively with Greek and Latin indices. We will indicate in the following as “target”
electrons the N bound electrons present in the molecule or molecular ion, while thescattered electron will be added to form an N + 1 electrons compound system. Here wedeal with the electronic problem alone, within the Born-Oppenheimer (BO) approximation, namely freezing the nuclei in some definite configuration (usually the equilibrium
configuration) while solving for the electronic wavefunction. The treatment of vibrations can be carried out by repeating the electronic calculations for different values of
the nuclear positions, followed by vibrational averaging or a vibrational frame transformation description. [47, 77]. Vibrational dynamics will be considered in more detail in
Chapters 3 and 6 where we will show also how to include non-BO effects [95, 104].
In principle, even considering just the electronic degrees of freedom the problem isstill a very complicated many-body interaction. The scattered electron is indistinguishable from the target electrons inside the molecule and it can in principle electronically
excite the molecule. For this reason different approximations will be required to makethe description of different processes feasible.
The final wavefunction of the system will have to obey scattering boundary conditions for the escaping electron, therefore at large distance from the molecule it shouldlook like an incident plane wave plus a spherical scattered wave, in order to be able

8
to apply the usual formulas to calculate cross sections and scattering amplitudes, as
described in many quantum mechanics textbooks, see for example Ref. [165]. This canbe done by matching to a long-range analytical form, as we will see in Sec. 2.5.
At this point a comparison of the electron collision problem and the more traditional bound state problem is opportune. The continuum problem is more complicatedbecause of the scattering boundary conditions that have to be imposed on the wavefunctions, but in principle the two are very similar many body problems. The quantum
chemistry community has made tremendous advances in the past 20 years (coupled cluster theory, DFT, linear scaling methods and so on, see Ref. [181] for an introduction
to the subject), and some of these are only now starting to be applied to the collision
problem. It will be auspicious for the future to foster a closer collaboration between
these two areas of research, especially as the interest in continuum processes involving
large biomolecules and nanostructures continues to strengthen.

2.2

Approximations
The starting point of the approximation procedure is to develop an ansatz for

the wavefunction of the total system. One of the simplest and most used methods is a
close coupling expansion [10], in which the N + 1 electron wavefunction is expanded ineigenstates of the target multiplied by an unknown function for the scattering electron.
Since usually the electronic excitations of the target are limited to the first few excitedstates, these expansions can be kept reasonably small, usually limited to a numberof electronic target states of 10 or fewer. The unknown function that multiplies thetarget wavefunctions is a continuum orbital for the scattered electron. These terms
in the expansion only represent channels that are open (this is not always true, e. g.

9in quantum defect theory (MQDT) as we will see in Sec. 2.5), allowing one electron
to escape to infinity. To describe closed channels another sum, this time over square
integrable functions (L2 functions) will have to be added to the wavefunction ansatz.
The close-coupling concept leads to a trial variational wavefunction that depends
parametrically on the nuclear coordinates, and is an antisymmetrized product of thetarget and scattering electron wavefunctions:
Ψγ = AX
γ0

Φγ 0 (x1 ...xN , R)φ0,γ 0 (x) +Xχj (x, x1 ...xN , R)

(2.2)

j

where x represents the spatial (~r) and spin (σ) coordinates of the scattering electron,
x1 ...xN are the coordinates of the target electrons, and R the nuclear coordinates. In Eq.
2.2, γ represents the set of quantum numbers that fully describe the state of the system,and the sum over γ 0 allows for different configurations of the compound system (target+ scattered electron) to contribute. The second sum represents the antisymmetrized
square integrable terms, which are also called polarization and correlation functions,because they allow in some sense a relaxation of the target electronic density under theaction of the incoming electron, whereas the first sum would keep the target electrons
strictly confined to their bound state basis set. In principle, the configuration statefunctions (CSF) that correspond to the Φ and χ in Eq. 2.2 contain spin degrees of
freedom, and the total wavefunction has to be antisymmetric. Since we will be mostly
dealing with closed-shell molecules the total spin state will be a singlet and we only
write explicitly the spatial degrees of freedom. In the following, integrals imply a traceover the spin degrees of freedom, where appropriate. In the remainder of this thesis thexi spatial and spin coordinates in Eq. 2.2 will be replaced by the spatial coordinates r~i .
If only the ground state configuration γ 0 in this sum is retained, the approximation

10
made is called static exchange. We can obtain the static exchange equation for the
scattered electron multiplying the Schrödinger equation from the left by the target
wavefunction and integrating over the variables of the bound electrons:Z

dr~1 ...dr~N Φ(r~1 ...r~N )(Ĥ − E)AΦ(r~1 ...r~N )φ0 (~r)(2.3)

This is valid for any form of the target wavefunction. In particular we consider aconfiguration-interaction (CI) type form:
Φ(r~1 ...r~N ) =X

ci χi (r~1 ...r~N )

(2.4)

i

where ci are the CI coefficients and the χi is a Slater determinant. For closed shelltarget molecules, we can simplify this equation considerably since the scattering orbital
will be forced to be orthogonal to the target orbitals for the Pauli exclusion principle.
We then obtain

N
X
X
1~2
φki (~r)(− ∇+ Vs − E)φ0 (~r) =
ci cj
2
i,j
k=1

Z

dr~0

φ∗kj (r~0 )φ0 (r~0 )| ~r − r~0 |(2.5)

where the φj (j ≥ 1) are the target molecular orbitals of a closed-shell molecule. The
electrostatic potential Vs is the averaged Coulomb interaction of the scattered electronwith all the other electrons and the nuclei
Vs (~r) =

X
i,j

ci cj

N Z
X

k=1

dr~0

φ∗ki (r~0 )φkj (r~0 ) X
Zα
.
−
0
~α |
~
r−R
| ~r − r |
α |~

(2.6)

The term on the right hand side of Eq. 2.5 is referred to as the exchange potential.
So far, we have presented the equations for static exchange calculations, whichmeans complete neglect of correlations. We will see that, when we use a polarization
potential, that term will contain a short range cutoff that can be interpreted as mimicking a correlation term. We deem it more appropriate to present here uncorrelated

11
equations, since this simplifies the derivation of the approximations that have been implemented here, and also because we never perform a CI calculation on the target +scattered electron system. For this reason it would not be pertinent to present fully
correlated equations. Furthermore, it is possible to see that the static potential (andalso the exchange potential in LDA approximation) depends only on the electron density, which allows us to treat the target molecule at a correlated level without changing
the formalism, since this only amounts to changing the Hartree-Fock electron density
to a correlated one. We will discuss the effect of this change in section 2.2.3. The methods we develop here have been published in Refs. [191, 192, 194] and the associatedcomputer codes are available through the CPC library [1].

2.2.1

R-matrix method

The R-matrix method is a well-established tool for problems where the continuum
portion of the spectrum of a Hamiltonian must be treated. In its usual implementation,
it involves diagonalization of the (Bloch-modified) Hamiltonian operator in a box subject
to some fixed boundary condition obeyed by the basis orbitals. The R-matrix box
partitions the space into two regions. There is an internal reaction zone, within which
all the short-range interactions are confined, and an external zone, where instead either
no potential is present or else there is a long range Coulomb or dipole potential (orboth), and the behavior of the solutions of the Schrödinger equation is very simple. In
some studies, other long-range multipole potentials are included in the external zone.
[12, 45] We use the R-matrix method in the eigenchannel form [74]. In this case we seek
those stationary states for which the logarithmic derivative of the wavefunction at thesurface of the R-matrix box is constant at every point. Refs. [48, 72, 109] derive a new12variational principle,b≡−∂log (rΨβ )Ψ∗ (E − Ĥ − L̂)ΨdV,= 2 VR∗∂r
V Ψ δ(r − r0 )ΨdV
R

(2.7)

for the logarithmic derivative of the wavefunction, starting from the Rayleigh-Ritz variational principle for the energy. If Ψ is discretized in some basis set inside a sphericalbox, within which all the short range dynamics is localized, this results in a generalized
eigenvalue problem for b:~ = 2(EO − H − L)C~~ = ΛCbΓC(2.8)

where Λ is the overlap of the basis functions calculated on the surface of the R-matrix
box and L̂ is the Bloch operator,defined as
1∂L̂ = δ(r − r0 ) r,
2
∂r

(2.9)

O is the volume overlap of the basis set and r0 is the radius of the box. The eigenvector~ represents the expansion coefficients of the basis set used. Both Γ and Λ are definedCin appendix A for the finite element basis set used in this work. It is possible to partition
the basis functions in two subspaces, closed and open, depending on whether their value
at the surface of the box is zero or nonzero [74]. This allows us to reduce the burden of
the solution of Eq. 2.8 to the easier task of solving a much smaller eigenvalue problemof type~ o = (Γoo − Γoc Γ−1~~ΩC
cc Γco )Co = Λoo Co b

(2.10)

in the open functions subspace, in addition to the large auxiliary system of equations:~ c = −Γco C
~o
Γcc C

(2.11)

13
where the subscripts indicate the matrix blocks. At the boundary of the R-matrix box
Ψ is matched to an external solution depending on the long range tail of the potential(Bessel functions for neutral molecules, Coulomb functions for molecular ions). This
allows us to calculate the reaction matrix K, from which the scattering matrix is derived
asS=

1 + iK
1 − iK(2.12)

Scattering cross sections can then be calculated in the standard manner. One difference
from the usual treatment performed in quantum mechanics textbooks [165] for spherically symmetric potentials, will be in this case the fact that we want to average overthe molecular orientations. From the scattering matrix S, the total (elastic) integrated
scattering cross section can be calculated as [41]:σ=π X| δll0 δmm0 − Slm,l0 m0 |22E lm,l0 m0(2.13)

where E is the energy of the incident electron and the S matrix indices are angularchannels, since here we do not include any excited target states [41]. A similar treatment
is possible in the case one wants to derive differential cross sections, and we can calculatethem using this approach employing formulas 46-47 of Ref. [41], one example is given
in Fig. 4.5.

2.2.2

Finite elements

The essence of the finite element method is the use of a basis set that is definedover small local regions. By this we mean that each basis function is nonzero only
within a small region, and it has a simple polynomial form. By using many “sectors”
or “elements” (the volume over which the local basis function is defined) though, it is14
possible to reproduce very complex features of the solutions to the differential equationof interest. An introduction to these basis sets and their practical implementation canbe seen in Refs. [19, 21]. We discretize Ψ using finite element polynomials in all three
dimensions. The basis set is a direct product of 4 cubic Hermite polynomials defined
locally in each sector for each dimension. The use of a spherical coordinate grid, inr, θ, φ makes the boundaries of the sectors simple and the three-dimensional integrals(the main bottleneck of these calculations) faster to calculate. The grid is represented
in Fig. 2.1.
The finite element basis set is composed of piecewise polynomials, which provides
advantages over a global basis representation. In particular one can treat potentials andwavefunctions of complicated form by simply reducing the size of the elements in which
the polynomials are defined, in those areas where fine features arise. In our case thebasis functions are third order Hermite polynomials, which allow us to achieve function
and derivative continuity, while still permitting a simpler implementation comparedto higher order polynomials. Each polynomial is defined in a hexahedral sector (acube in the rescaled variables used for the evaluation of the integrals), and since thewavefunction is discretized in terms of finite elements in all three dimensions, the basis
set is a direct product of 4 polynomials in each dimension per sector, which means 64
basis functions are defined in each sector.
In finite element analysis (FEA) the polynomials are matched with the ones inneighboring sectors to ensure functional and derivative continuity (and mixed derivative
continuity also, in multidimensional FEA). Each sector has 8 physical nodes (at theedges of the cube) and the basis functions defined in the sector have coefficients (to be
determined by the solution of the Schrödinger equation) that represent the value of the

15
wavefunction, or its derivatives, at the nodal point. In the language of finite element
analysis, a node is the vertex of one of the sectors into which the three-dimensional spaceis divided. The matching at the boundary of each sector is imposed when assembling
the global Hamiltonian matrix from the local ones. The global index of functions that
correspond to the same node and quantity (e. g. derivative) in neighboring sectors has
to be the same. Hence their matrix elements have to be summed together. Details of
the procedure are given in Appendix A.

2.2.3

Potential and DFT

Using an approach derived from Refs. [42, 139], we approximate the exchange
integral (that is nonlocal), by a local form using free electron gas (FEG) orbitals [93],i.e. plane waves, for the target molecule and using the first order Born approximation~

φ0 = N eik·~r(2.14)

for the scattered electron. The arbitrary normalization constant N is unimportant andit disappears as soon as we express the exchange functional as a product of a local
exchange potential times the scattered wave. After these substitutions are made, it ispossible to evaluate the integral on the right hand side of Eq. 2.5 analytically, obtaining
a local potential of the form
2
Vex (~r) = − kF F (η),π

(2.15)

whereas the Fermi momentum kF (the momentum of the electron that is at the top of
the Fermi sea in a free electron gas) is:kF (~r) = [3π 2 ρ(~r)]1/3 .

(2.16)

16

3

2.5

2

1.5

θ

1

0.5

0
0

1

2

3

4

5

6

7

8

r (a.u.)Figure 2.1: From this two dimensional cut in the radius r and the polar angle θ of thefinite element grid (for a CO2 target), the finer mesh is apparent near the oxygen nuclei
localized at r = 2.19 a.u. and θ = 0 and π respectively, while the carbon is located at
the center of the grid (r = 0). A drawback of this type of separable grid in sphericalcoordinates is that it is more fine than necessary at 1.5 < r < 3 and 0.3 < θ < 2.8
radians. From Ref. [192].

17
The other functions present in Eq. 2.15 are
F (η) =

¯1 + η ¯
1 1 − η2¯+
log ¯¯24η
1 − η¯
¯

η=

¯

k
,
kF

(2.17)(2.18)

where k is the modulus of the momentum of the scattered electron. It should be noticed
that the exchange potential in Eq. 2.15 is energy dependent.
Many functionals of this form exist [139], with minor differences in the expression for k, the scattering electron wavenumber. The functional we have used mostsuccessfully is the Hara exchange [81] wherek=q2(E + I) + kF2(2.19)and I is the ionization energy of the molecule while E is the energy of the incident
electron. This emerges from the assumption that the scattered electron and the electronin the highest energy bound state (the Fermi electron, which has momentum kF ) movein the same potential field; Vex then depends only on ~r, through the electron density
ρ(~r), as a local potential, and on the energy, through the functional dependence of the
momentum k as approximated in Eq. 2.19.
We have also experimented with other functional forms of the exchange interaction(still based on a FEG approximation). One in particular is the Slater exchange [173],
derived by averaging the function F (η) over the momenta of all the electrons up to the
Fermi level, which has often been used to calculate bound states in atoms and molecules.
However the results using Slater exchange are unsatisfactory for continuum scattering
states, presumably owing to the neglect of the energy dependence in this model.
Since our main goal is to treat low energy scattering processes (0-10 eV) we
linearize the energy dependence of the functional in Eq. 2.17, in order to calculate the18exchange potential matrix elements at all energies at once. For a molecule like CO 2 , thematrix element calculation requires around 2 hours on an Alpha 500 Mhz workstation.
The next step is the solution of the linear system and the determination of the scatteringobservables, which requires approximately 15 minutes per energy desired, for a basis set
size of 33000. This step is trivially parallelizable, of course.
We have recently added to our computer code the capability to use a parameterfree correlation-polarization potential [35, 62], based on density functional theory (DFT)
ideas. As shown in Ref. [102] the polarization-correlation contribution is physicallyrelated to the distortion-relaxation effect on the molecule generated by the incoming
electron. This is extremely important for an accurate description of the scattering
process. The long range part of this potential is a simple multipole expansion, of which
we retain only the induced dipole polarization terms:
Vpol = −1[α0 + α2 P2 (cos θ)]
2r4(2.20)

where α0 and α2 are the totally symmetric and nontotally symmetric components ofthe polarizability tensor, and are calculated ab initio using electronic structure codes.
While the quadrupole term in the polarizability tensor has 5 components, we neverresort to its full form, since we have verified numerically that the second term in thepolarization potential (α2 = 3αzz −T r(α)) is already not very important, and we omit italtogether when calculating scattering from larger molecules. To extract the molecularelectron density, gradient, Laplacian, and the electrostatic potential needed to construct
the scattering potentials, we use GAUSSIAN 98 [150], but MOLPRO [203] has thesame capabilities and exploits the same file formats, as does PC-GAMESS [71]. The
latter is also freely available.

19
In the volume where the electronic density of the target is not negligible, thispotential is nonlocal. The interaction can be approximated again as a local potential,different forms of which have been suggested in the literature. The one we use is based on
DFT (specifically on the LYP potential of Ref. [107]) and it has yielded reliable results inthe work of Gianturco and coworkers [117]. This form makes use of the electron density,
its gradient and Laplacian, which have to be calculated for each target molecule. The
short and long range potentials are matched unambiguously at the innermost crossing
point, whose radius is dependent on the angles.
A derivation of this short-range correlation-polarization potential can be sketchedalong the lines of Refs. [35, 62, 107]. We start by choosing a particular form for thecorrelated wavefunction for an N-electron system. Taking xi to denote both the spin(σi ) and spatial (~ri ) coordinates of the electron, this reads:
Φ(x1 ...xN ) = ΦHF (x1 ...xN )

N
Y

i<j

[1 − φ(~ri , r~j )](2.21)

where
2 2

φ(~
ri , r~j ) = e−β r [1 − Θ(R)(1 + r/2)](2.22)

is a correlation function that obeys the cusp condition at the coalescence point. HereR = 21 | ~ri + ~rj |, and r =| ~ri − ~rj |. β = qρ1/3 can be related to the excluded volume(see Eq. 5 of Ref. [62]) and q has to be determined. The form chosen for the correlation
function determines
P2 (~r1 , ~r2 , r~0 1 , r~0 2 ) = P2HF (~r1 , ~r2 , r~0 1 , r~0 2 )[1 − φ(~r1 , ~r2 ) − φ(r~0 1 , r~0 2 ) − φ(~r1 , ~r2 )φ(r~0 1 , r~0 2 )](2.23)
for the second-order (or two-electron) density matrix traced over spin, where P 2HF is
the second-order Hartree-Fock density matrix. Making now the approximation that

20
P1 (~r1 , r~0 1 ) = P1HF (~r1 , r~0 1 ) for the first order (or one-electron) density matrix, the correlation energy is derived as
Ec =

1
2

Z

P2HF (~r1 , ~r2 )[φ(~r1 , ~r2 )2 − 2φ(~r1 , ~r2 )]1
dr1 dr2
r12(2.24)and Θ becomes a function only of β (see Ref.[35]). The correlation energy is finally
recovered by expanding P2HF (~r1 , ~r2 ) to second order in the interparticle coordinate. At
this point though, the correlation energy involves a gradient of P2HF (~r1 , ~r2 ). Throughsome further work, Lee Yang and Parr [107] recovered a final expression that is only
a function of the one-electron density, and therefore more akin to the spirit of DFT.
Essentially the steps required involve first expressing the second-order Hartree-Fockdensity matrix in terms of the first order one (see Eq. 4 of Ref. [107]) and secondly
expressing the gradient in terms of the Hartree-Fock local kinetic energy. This finallygives an expression like Eq. 10 in Ref. [107]. At this point, to obtain the correlation
potentials we use, it is necessary to perform a functional derivative of the correlationenergy with respect to the electron density:
Vc =

8
∂
Ec = −a(F 0 ρ + F ) − abCF ρ5/3 (G0 ρ + G)
∂ρ
3

ab 00
[G ρ | ∇ρ |2 +G0 (3 | ∇ρ |2 +2ρ∇2 ρ) + 4G∇2 ρ]4ab− [3G00 ρ | ∇ρ |2 +G0 (5 | ∇ρ |2 +6ρ∇2 ρ) + 4G∇2 ρ]
72−(2.25)F = (1 + dρ−1/3 )−1
G = F −5/3 ecρ
CF =

−1/3

(2.26)

3(3π 2 )2/3
10

for closed shell molecules, these expressions are identical to the ones in Ref. [62] and areshown here for completeness. The parameters a, b, c, d are chosen by a fitting procedure

21
to the helium atom, and were shown to be effective in predicting energy correlations
for atoms and molecules [35, 107]. Note that it is possible to include rigorously thepolarization-correlation physics in an effective potential approach as in Ref. [128].
All the information needed to construct the local potential described in Eqs. 2.252.27 can be extracted from standard ab initio quantum chemistry codes; in this work
we have used GAUSSIAN 98. The electrostatic potential and the electron density(needed to construct the exchange functional) for the target molecule are calculated ona uniform cubic grid at a CI (singles and doubles) level for the molecules presented here.
The difference in using an electrostatic potential and density calculated at the RHF levelor at the CI level for CO2 at its equilibrium geometry amounts roughly to a differenceof 10% in the calculated eigenphase sum and overall magnitude of the elastic crosssections. These calculations usually require a minimal amount of time, of the orderof ten minutes per nuclear geometry for CO2 on the aforementioned computational
platform. The potentials are then interpolated on the three-dimensional quadrature
grid using fifth order splines. We avoid interpolating over the nuclear singularities by
premultiplying the potential at each grid point by the product of the distances of thepoint from all the atoms, then performing the interpolation over this quantity that does
not possess Coulomb singularities any longer.

2.3

Photoionization
The relationship between electron scattering and photoionization is a close one[41]. The scattering wavefunction relative to a photoionization process is essentially thetime-reverse of the one for the scattering of an electron with a positive ion [74]. Thus
it is useful to examine the resonant structures of the electron with the ionic state, since220

Potentials (a.u.)

-2

-4

-6

-8

-10
-2.5

VCoulomb
Vexchange
Vpolarization
-2

-1.5

-1

-0.5

0

0.5

1

1.5

2

2.5

R (a.u.)Figure 2.2: The three terms of the potential for a N2 molecule. The exchange potential
is large only at the nuclei (at r = −1.094 and r = 1.094 a.u. in the equilibrium
configuration of the molecule) where the static potential is singular, so V ex is always
much smaller than Vs . On the other hand the polarization potential becomes importantin the outer zone, where the electron density of the molecule goes to zero. From Ref.
[192].

23
photoionization allows to distinguish between different initial wavefunction symmetries.
In the independent electron approximation, the photon extracts the electron from justa single Hartree-Fock orbital, the initial state. The relevant matrix element for thephoton absorption from a molecule is:~ | ψ0 >,M0,f =< ψf | eik·r ²̂ · ∇(2.27)

where ψf is the final state. In our case the final state will have only partial wave
indices, since our scattering calculation is done at the static exchange level. For opticaltransitions, the dipole approximation can be used, leading to dipole matrix elements:
(−)

(−)

d0,l0 m0 =< ψl0 m0 | ²̂ · ~r | ψ0 >,(2.28)

using the length form of the dipole operator, where ²̂ is the polarization of the incident
photon, ψ0 the initial wavefunction (the orbital from which the electron is ionized), and(−)ψl0 m0 the final state of the continuum electron. This is formed from R-matrix eigenstatesψl0 m0 usingψ (−) = ψ(I − iJ)−1 ,(2.29)

where I and J are matrices generated from the Wronskians of the solutions inside andoutside the R-matrix box, as described in Eqs. 2-31 and 2.41 of Ref. [74], generating a
wavefunction with incoming wave boundary conditions:(−)ψ l 0 m0 =1
1†
Ylm (θ, φ)( √ fl+ δll0 δmm0 − √ fl− Slm,l0 m0 ),
i 2
i 2
lm

X

(2.30)

where f ± are incoming and outgoing wave solutions (exponential functions with unit
coefficient), respectively, and S † is the hermitian conjugate of the scattering matrix.
(−)

The total cross section can be easily calculated from the d0l0 m0 matrix elements by24
performing a rotational average. This simplifies after expressing the dipole operator as²̂ · ~r =

X 4π
mγ

3

∗
rY1mγ (θ, φ)Y1m(θ0 , φ0 ).
γ(2.31)

where the first spherical harmonic refers to the vector ~r, while the second refers to thepolarization vector in the body frame, and mγ gives the projection of the polarization
vector on the molecular z axis. After angular integration, the final cross section is:
σ=N

4πω X (−) (−)∗
d 0 0d 0 0
3 l0 m0 0,l m 0,l m(2.32)

where N is the number of degenerate electrons that can be knocked out by the photon.
For example, N = 4 for a fully occupied π orbital in CO2 . The cross section output
from the code has to be multiplied by the orbital occupation number. As an alternative
to Eq. 2.28, the velocity form of the dipole matrix element,(−)

d0,l0 m0 =1(−)
~ | ψ0 >,< ψl0 m0 | ²̂ · ∇ω

(2.33)

can be used in the code, although this will be slower because of the need to differentiate
many continuum wavefunctions. In fact, in the case of many initial wavefunctions, it
is too expensive to store all of their gradients, so the final wavefunctions have to bedifferentiated at each energy point.
Notice that in photoionization of a closed-shell molecule, the construction of our
exchange potential will not be strictly correct, since we only derived our exchange
functional for closed shell molecules and the resulting ion is a doublet. The final error
made in this case will however be quite small, due to the dominance of the electrostatic
ionic potential everywhere in space.
The calculation of dipole matrix elements for a returning electron recolliding witha molecule, therefore with boundary conditions opposite to photoionization, is useful also

25in the three step model [113] calculation of high harmonic generation (HHG) spectra [36],that are generated by femtosecond and sub-femtosecond pulses. This area of research is
growing at a very rapid pace and interesting applications to map the electronic dynamics
of molecules undergoing chemical reactions have already been published [13]. This isripe for a future application of the methods we are developing here.

2.4

Other methods
In the work included in this thesis, there are many different kinds of processes

that we are trying to model, and they will require a number of different approaches. The
accuracy needed in calculations for small molecules is much higher than for the large
DNA bases, where the approximate treatment we have just discussed is good enough,since more accurate calculations are at present too large to perform. In the following
sections a description of the other methods that have been used in this work to calculatescattering quantities will be given and the different schemes will be compared to assesswhen they should be used and what their weaknesses are. To develop methods that
can accurately solve the electron-molecule scattering problem is a serious job in itself,
and the approaches and programs we describe have been developed by large teams ofscientists over the past twenty or so years.

2.4.1

Ab initio methods

A possible alternative route to calculate quantum defect parameters, for scatteringof an electron from an ion, is to use bound state quantum chemistry programs, whichallow very large configuration interaction (CI) expansions to be used, and therefore highaccuracy for the first few low-lying excited states of the N + 1 electron system [132].

26
We have used this approach to calculate Rydberg state Born-Oppenheimer surfaces
and K-matrices relevant for HCO+ dissociative recombination calculations [104], i.e. by
calculating the first few Rydberg states (up to n = 4) and then neglecting the energy
dependence of the higher state eigenquantum defects.
The positive side of this approach is the high precision it is possible to reachusing very large configuration spaces (of the order of a million). On the other hand, thisprocedure is limited by two factors: the higher lying excited states will be less and lessconverged (and therefore the Rydberg states above, say, n = 4 will be incorrect) andthe scattering quantities that are of interest here (K-matrices through equation 3.1) will
have to be obtained through a fitting procedure, which can be both cumbersome andnonunique. Clearly, above the ionization threshold this approach cannot be used, since
there would be a variational collapse of the states onto the ionic surface energies, andan incorrect description of the resonances, due to the bound state boundary conditions
imposed that cannot represent the continuum correctly. This approach will be used in
the dissociative recombination calculations of Chapters 3 and 6.

2.4.2

Rmatrix-UK

Small molecules often need more accuracy than we are able to provide with themethods just outlined. For this purpose we have used a different approach that involves
the R-matrix UK polyatomic codes developed by J. Tennyson and coworkers [135, 136].
The main difference with the approach we have described so far is that the solutionof the inner region is a truly many-particle problem, developed starting from quantum
chemistry configuration interaction (CI) software. Within the limit of the basis sets
used, these codes can give an exact solution to the electron interaction with the many

27
electrons present in the molecule, at least in principle.
The calculations are usually done in this case starting from a Hartree-Fock mean
field description of the target molecule, followed by a CI step to get better target
wavefunctions. Many target states might be needed if electronic excitations are considered. The next step involves adding an electron to the molecule and diagonalizing
the full Hamiltonian, including the interaction of the scattering electron with the target
molecule. A wavefunction of the form:
ψk =

X
ij

φi (x1 ...xN )uij (xN +1 )aijk +X

χi (x1 ...xN )bik

(2.34)

i

is assumed, where φi (x1 ...xN ) are target CI configuration state functions, uij (xN +1 ) arecontinuum wavefunctions, while χi (x1 ...xN ) are L2 functions, a and b are the variationalcoefficients. This is essentially a close coupling expansion with polarization/correlation
configurations that allow us to describe the added electron in more compact configurations inside the target virtual orbital space, in a sense allowing the target to relax
under the influence of the scattered electron, much as our model polarization potential
2.2.3 does. This wavefunction form goes beyond the static exchange approximation inthat it allows electronic excitation of the target. We use this approach to calculatepotential energy surfaces and K-matrices for dissociative recombination calculations.
The limitation of this approach is one of size, in that when the number of electrons
in the scattering problem becomes large the scaling of CI calculations (even with inclusion of just single and double excitations) goes like N 6 , where N is the number ofsingle-electron basis functions, and is proportional to the number of electrons. Alreadya simple HF calculation scales like N 4 . This unfavorable scaling prevents this approachfrom being extended to carry out accurate calculations for molecules with more than28
about 10 atoms, especially if what is of interest is the effect of bond breaking on thecontinuum states. Homolytic bond breaking (in which each of the fragments retains
a bonding electron) inherently requires more than one reference configuration, sincethe Hartree-Fock wavefunction dissociates incorrectly in this case, and in general many
configurations will be needed for an accurate description of the process.
Two other methods that are based on the same philosophy of using a full manyelectron description of the electron-molecule interaction, and which are competitive with
the Rmatrix-UK code. These are the Schwinger variational approach of McKoy and coworkers [205] and the Kohn variational code of McCurdy-Rescigno and coworkers [156],
which share roughly the same advantages and limitations as the code we just described.
For the future, promising methods that use a DFT-based description of the interaction, using more accurate functionals than we have at present [148], which importantly
eliminate the local density approximation, are being developed and applied to scattering
problems [198]. The scaling of these methods is much better than wavefunction methods, and they might be the only way to treat interactions of a continuum electron withvery large targets in the future, to sufficient accuracy. Much as in bound-state quantum chemistry calculations, the success of DFT calculations could prove critical for
expanding the range of systems that can be treated, especially for new and extremelyinteresting applications to biomolecules [24].

2.5

QDT and long range interactions
The nature of the electron scattering problem is such that we can devise twozones in configuration space that have distinct properties. In the short range zone, theproblem is essentially many body in nature, whereas the long range zone effectively

29
reduces to a one-body problem of scattering from a potential. When the outer electron
is outside the area where the target electron density is confined, there is no nonlocal
exchange interaction. The first zone is the most complicated to treat, and we haveseen how R-matrix techniques for example can be used for that purpose. The outerzone instead has a feature that is very important: at sufficiently large distances, one
can solve for scattering of a particle in the long range potential and get an analyticalform of the wavefunction. The computational cost of solving the Schrödinger equationout to very large distances is thus reduced, and equally importantly, it factors out
analytically from the S-matrix the long-range part of the phase shift. This key idea
allows us to obtain quantities like S-matrices, phase shifts, time-delays, etc., which
depend only on the short-range potential, and therefore are much less energy dependent
than the ones that also include the long-range phase accumulation. This can in turn savemuch computation time when dealing with processes for which the scale of the energy
dependence of the short-range part of the S-matrix is much weaker than for the longrange part. In the treatment of electron collisions with rotating/vibrating molecules
this concept is useful, since the energy scales of the different degrees of freedom are so
different, as we will see in Ch.3. This is one of the fundamental ideas of quantum defecttheory (QDT), which will be important in the following, especially when we will deal
with dissociative recombination reactions.
QDT is a tool that allows to tackle a class of problems that are essentially intractable with any other method, for example the coupling of many degrees of freedom
with the infinite number of Born-Oppenheimer potential energy surfaces that arise in
scattering from a positive molecular ion. It is conceptually very simple, but in the literature [74, 75, 77, 169] it has been usually obscured by complexities associated with a30large number of parameters, basis pairs and so on.
The key ideas are essentially three. The first one we mentioned already is thedivision of the scattering electron coordinate space into two areas, short range and longrange. The long range part can be treated analytically and effectively removed from theS-matrix to make it a smooth quantity. The second idea is that the boundary conditions
on closed channels should be imposed as late as possible in the course of a calculations,again in order to obtain quantities that are smooth in energy, since the sharp energy
dependence at ionization thresholds has been eliminated. The third idea is that there
is a correspondence between the phase shift at positive energy and the parameter that
identifies the position of a bound state (the effective quantum number), at negative
energy. For Coulomb systems that automatically possess high lying Rydberg boundstates this is easy to see in that the bound states will look more and more like the ionicstate to which they converge plus a loosely bound electron. In the case of single channel
problems, their energies will be given by the Rydberg formula~ = U + (Q)~ −Un,i,Λ (Q)i1~ 2
2(n − µi,Λ (Q))

(2.35)~ represents all the vibrational coordinates. In the case of multichannel problems
where Qthe quantization condition Eq. 4.5 has to be applied to determine the position of thebound states. The interaction of the electron with the ionic core in these very extended
and high energy weakly bound states will be more and more similar to that of a free,
positive energy electron. In Eq. 2.35, n is the principal quantum number, i is theindex of the ionic electronic state, while Λ is the projection of the electronic angular~ is the quantum defect for thatmomentum on the molecular symmetry axis and µiΛ (Q)geometry and molecular symmetry. Λ is a good quantum number for linear molecules,

31
and we will use it in HCO+ to label states in bent geometries as well, even though theclassification will be approximate. A derivation that shows that the quantum defect informula 2.35 can be seen as a low-energy electron scattering phase shift can be found in
Sec. IIB of Ref. [74]. An important consequence of treating the bound-state physics asa scattering process is that we can use the S-matrix to determine bound-state positions,when the problem is multichannel in nature:†(Scc− e2iβ )B = 0(2.36)and the roots of the determinantal equation will give the Rydberg energies. In this† is the part of the S-matrix relative to the closed channels (in which fragequation Scc

mentation all the way to the asymptotic region is energetically forbidden), β is a diagonal
matrix of long range phase parameters [74, 75]. Note also that all channels are closed
when the energy level is truly bound. In an energy range where open channels are alsopresent, see Eq. 3.10. It is important to restate that QDT has to do with the long-range
part of the scattering process, and cannot give by itself the short-range quantities which
are needed to complete the modeling of the interaction. These must be calculated byeither ab initio bound state quantum chemistry calculations or R-matrix (or similar
approaches) for scattering processes. Both these avenues will be explored in Ch. 3.

2.5.1

External potentials

Many external potentials have been treated from a QDT standpoint: Coulomb,
dipole [75], Coulomb plus dipole [180], polarization [199] and dispersion forces in atomatom scattering [27]. Here we will be mostly concerned with Coulomb and dipole potentials. The long range form of the wave function will be always some linear combination

32
of a base pair consisting of a regular (at the origin) and an irregular solution of the
Schrödinger equation in the pure long range potential. Since a linear second order differential equation has two linearly-independent solutions, any solution can be expressed
as linear combinations of those two. We will call these functions f and g respectively,
without specifying their precise functional form, which will vary with the external potential. In the case of zero long-range potential, the external solutions are the familiar
spherical Bessel functions.
Given these two functions, matching at the boundary of the internal region with
the internal R-matrix of Sec. 2.2.1 generates the K-matrix of Eq. 2.12 in the followingway:
f − f 0RK=g − g0K(2.37)

here f , g, f 0 , g 0 are diagonal matrices of the matching functions and their derivatives
at the matching point, as described in Sec. IID of Ref.[74]. In the following we will
always assume that only one electron at a time can be in the external region, since theenergy will always lie below the threshold for double escape and we will not consider
such processes further. In electron scattering, the threshold for double escape is the firstionization threshold, which can be as low as 10 eV in some of the neutral molecules weconsider in this thesis. We neglect the double escape process nonetheless throughout
this work.

33
2.5.2

Coulomb potential

In the case of a Coulomb potential, the wavefunction at long range for the external
electron will be a combination of Coulomb functions [169]ψ l 0 m0 =1X
Ylm (θ, φ)[fl (r)δll0 δmm0 − glm (r)Klm,l0 m0 ]r lm

(2.38)

this form of the wavefunction is said to obey K-matrix boundary conditions. The
functions fl and gl are energy normalized regular and irregular Coulomb functions, thatbehave like sin(kr + lnkr + η) and − cos(kr + lnkr + η) for r → ∞, where η is the longrange Coulomb phase shift (see Eq. 2.39 in Ref.[74]). It is important to notice that thefull wavefunction for the system will be a Slater determinant of the wavefunctions in Eq.
2.38. Alternatively, the wavefunction 2.38 can be reorganized in terms of the S-matrix
solutions, using a complex linear combination of fl and gl , or in the eigenchannel form
involving eigenstates of the K-matrix, showing more effectively the short range phaseshift µα that is not incorporated in the Coulomb functions:ψα =1X
Ylm (θ, φ)Ulm,α [fl (r) sin πµα − gl (r) cos πµα ]r lm(2.39)

This phase shift µα is termed the eigenquantum defect.
At negative energy, in the bound state region, the form of the wavefunction is effectively unchanged, what varies is the asymptotic form of the Coulomb wavefunctions,which now behave as decaying and growing exponentials. Of course the exponentiallygrowing part is unphysical, but the growth is killed when the correct boundary conditions are imposed at the end of the calculation. All of this makes this long range formof the wavefunction very convenient and flexible. Another base pair that appears in thefollowing is (f 0 , g 0 ), which are analytic functions of the energy and will be important

34
at very low negative energies since, for ν < l, (f, g) are complex. The transformation
between these two pairs of solutions is linear and it is defined in Eqs. 2.15a and b ofRef. [74], the long range parameters of the transformation can be found in Ref. [170]where (f, g) are called (s, −c) and (f 0 , g 0 ) are called (f, −g).
2.5.3

Dipole physics

Since many of the molecules we consider have large electric dipole moments, there
is a need to consider the long range effect of the dipole field on the scattered electron.
Two possible options might be considered, either extending the boundary of the Rmatrix box far out to a region where the dipole potential is very small, which would be
extremely time-consuming for our calculations, or matching to outer region functions
adapted to the dipole interaction. We choose this second route and, since in the region
where only the dipole potential is dominant the Schrödinger equation is
−

1 d2
1
− E + (l2 − 2D cos θ)ΩN = N (N + 1)ΩN
2 dr2
2

(2.40)

following the example of Clark [33], we define the matrix of the operator:(l2 − 2D cos θ)ΩN = N (N + 1)ΩN

(2.41)

where l is the angular momentum operator, θ is the angle between the incoming electron
and the dipole direction, D is the dipole moment, N (N + 1) and ΩN are eigenvalues
and eigenfunctions. We expand ΩN in a basis of spherical harmonics to diagonalize thesystem in Eq. 2.41. The order of the spherical Bessel functions that are matched in theouter region will be now N (not an integer, in general) rather than the usual orbital
angular momentum quantum number l. Since the dipole moments of the molecules

35
in question are very large, the dipole plus centrifugal potential can become attractive
for the first few channels which may thus have a complex N . In such cases we define
N = −1/2 + iβ and the matching functions will become [33]:j̄N (kr) =rn̄N (kr) = −π
1
Im(Jiβ (kr))
2r sinh 12 πβ(2.42)π1
Re(Jiβ (kr))
2r cosh 21 πβ(2.43)

r

where Jiβ is a cylindrical Bessel function. This allows us to keep the functions in Eq.
2.42 always real, which therefore gives real K-matrices.
It should be mentioned, however, that at extremely low energies these functions
oscillate rapidly in energy as sin(β ln kr), giving rise to K-matrices that are not smoothfunctions of energy near threshold. Defining the base pair as in Ref. [75] solves the
problem, but since we are not interested in energies below about 0.5 eV usually, thefunctions in Eqs. 2.42-2.43 will be sufficient.
After the matching is carried out, an analytic propagation procedure (see Eq.
8-16 of Ref. [33]) is used to get the scattering wavefunctions at a large radius. This isimportant since then the scattering observables, phase shifts and cross sections, must be
expressed in terms of spherical Bessel functions of integer order, not the dipole-adapted
ones in Eqs. 2.42-2.43.
Other multipole contributions outside the R-matrix box must be described in adifferent way, specifically with a propagation or perturbation method. We have tested
an R-matrix propagation approach [15, 127] to include the quadrupole contribution.
For the DNA components this contribution has been seen as very small, and therefore
it has been neglected in the final results, while it has been included in some of the
UK-R-matrix results we have obtained for HCO+ .

36
The dipole plus centrifugal potential is attractive if the value of the dipole moment
is larger than a critical value (Dc =1.625 Debye for a nonrotating dipole). In this casethe dipole interaction can bind the electron all by itself. In general, when rotation is
included, the critical value of the dipole moment to have a bound state predominantly inthe dipole field is around [88] 2-2.5 D and the number of dipole-bound states is finite. In
the case of uracil [43, 164, 176], such a dipole-bound state exists at roughly 0.1 eV below
the neutral ground state energy, at the equilibrium geometry of the target molecule.

2.5.4

Perturbation theory

The multipolar potential outside the R-matrix box can also be seen as a perturbation, since it is small enough usually that this approximation will be accurate. To
this extent, an approach based on a Green function has been developed [73, 175]. We
consider the partitioning of the Hamiltonian into a main part and a perturbationĤ = Ĥ0 + V(2.44)

and we use a solution obeying K-matrix boundary conditions to Ĥ0 of the type:
M (r) = f0 (r) − g0 (r)K0(2.45)

where all the quantities are matrices and in particular f0 and g0 are diagonal matrices ofunperturbed long-range functions and K0 is the unperturbed K-matrix. We can define
the Green function for the system to be:(Ĥ0 − EI)G(r, r 0 ) = −δ(r − r 0 )I

(2.46)

where I is the identity matrix, then the solution of the full Hamiltonian will involve an
integral over G and V .

37
To calculate the Green function, we can start from a form of the typeG(r, r0 ) = P (r)QT (r0 )A f or r < r 0(2.47)G(r, r0 ) = Q(r)P T (r0 )B f or r > r 0(2.48)

where P and Q are two sets of linearly-independent solutions of the equation. Then
we have to impose continuity at r = r 0 and derivative discontinuity of −2µI at thesame point, where µ is the reduced mass, which will be unity in our case. Once theserequirements are met, we end up with:G− (r, r0 ) = −2µP (r)[W (Q, P )]−1 QT (r0 )(2.49)

where W (Q, P ) is the Wronskian of the solutions, the advanced Green function will look
like the transpose of the previous expression. Now, for the unperturbed problem

P (r) = (f0 (r) − g0 (r)K0 )(I + K02 )−1/2(2.50)

Q(r) = (f0 (r)K0 + g0 (r))(I + K02 )−1/2(2.51)

the Wronskian will be, after some algebra:
W (P, Q) =2µIπ(2.52)and the solution will formally be:
M (r) = P (r) +Z ∞G(r, r0 )V (r0 )M (r0 )(2.53)Z ∞G(r, r0 )V (r0 )P (r 0 ).

(2.54)

0and to first order
M (r) = P (r) +0

38
Now, we substitute in 2.54 the Green function 2.48, and P, Q from Eq. 2.51. Upon
introducing the integrals
Z ∞If f =
If g =Igg =f0 (r0 )V (r0 )f0 (r0 )dr0

(2.55)

f0 (r0 )V (r0 )g0 (r0 )dr0

(2.56)

g0 (r0 )V (r0 )g0 (r0 )dr0 ,(2.57)

K pert = AB −1(2.58)

r0
Z ∞

r0
Z ∞
r0

the perturbed K-matrix will be:

where
A = K0 − π(If f − K0 Igf − If g K0 + K0 Igg K0 )B = I + πK0 (If f − K0 Igf − If g K0 + K0 Igg K0(2.59)

when the perturbing potential vanishes we recover the unperturbed form. The advantageof this form of the K-matrix over, for example, the treatment in Ref. [12] is that thepoles are shifted with respect to the unperturbed poles. We will show in Ch. 4 how wellthis method works, compared to an R-matrix propagation, for positive and negative
energies, for a complicated system like an electron in the field of HCO+ .

Chapter 3

Methods: vibrations, rotations and dissociative recombination

This chapter describes the methods we employ to treat vibrational motion andnuclear dynamics in general. Particular attention is paid to dissociative recombination
reactions (in this chapter and Ch. 6). Sec. 3.2 shows how we obtain the electronicparameters adopted in the calculations, and it describes the different representations(adiabatic/diabatic) of the electronic states and their use. Symmetry-distortion effects,
such as Jahn-Teller and Renner-Teller, are considered, as well as the way they can accelerate the dissociative recombination (DR) rate. A brief survey of frame transformation
techniques is given in Sec. 3.3, while our treatment of the molecular dissociation isshown in Sec. 3.4. Our formulation can be adapted in principle to any kind of molecule,but some specific choices (especially the vibrational coordinates and the form of themodel K-matrix) are tailored for HCO+ , which is the system we have studied in detail.

3.1

Introduction
The dissociative recombination (DR) of small molecular ions that collide withelectrons plays an important role in interstellar diffuse and dense clouds, see Ref. [53]for a recent review on theory and experiment in this field. It is well known that theseclouds constitute building material for new stars. The importance of DR cross sections

40
as parameters in star formation models, for such astrophysically relevant ions as H +3and HCO+ , is one of the reasons why DR has been extensively studied in laboratory
experiments [105]. For diatomic molecules the process is well understood and describedtheoretically [65, 89]. Until recently [94, 95, 143] though, theory was unable to model
the DR process in triatomic ions, except in cases where a neutral dissociative statecrossed the ionic Born-Oppenheimer surface in the Franck-Condon region, generatinga rapid rate [143]. In triatomics, a key complication is that multiple vibrational androtational degrees of freedom must generally be taken into account. Beyond the greater
computational burden of treating more dimensions quantum mechanically, the additionof these new degrees of freedom can also lead to new conceptual issues, related to thedegeneracy of vibrational or rovibrational levels in certain triatomic ions (e.g., in both+H+
3 and HCO ). This can cause an intrinsic instability of the corresponding neutral

molecules that forces them to distort away from the symmetric geometry.
Consider an incident electron that interacts with a closed-shell triatomic ion having a degenerate vibrational mode. If the symmetry group Γ of the resulting neutral
complex has at least one degenerate irreducible representation, the electronic partial
wave components with angular momentum l > 0 typically contain at least one such representation, whereby the corresponding electronic states of the neutral molecule are degenerate in the clamped-nuclei approximation. Due to the Jahn-Teller theorem [9, 115],
this degenerate electronic state strongly interacts with a degenerate vibrational state ofthe Γ group. If the molecule is nonlinear, the interaction leads to the Jahn-Teller effect,a conical intersection, in which the two states have an off-diagonal coupling which islinear in the distortion coordinate. If the molecule is linear, the off-diagonal couplingof lowest order is quadratic, giving rise to a Renner-Teller [126] (glancing) intersection.

41
It has been well established in the literature that these distortion effects are a generalcharacteristic of potential energy surfaces for polyatomic molecules, and taking theminto account has explained many spectra and fast reaction rates. Conical intersections[20, 208] can generate a “funnel” effect that leads to rapid passage of a wavepacket
from the upper surface to the lower surface, thus providing a preferred pathway for thereaction to occur.
In electron-induced reactions, electron capture can be followed by dissociation or
autoionization of the recombined system. Since the degenerate electronic state causes
the instability of the symmetric configuration, the recombined molecule quickly distorts
away from the symmetry point to remove the degeneracy, after which the autoionization channel typically becomes energetically closed, and the molecule eventually finds
a pathway along which it can dissociate (if there is an electronic state of the molecule
which is open for dissociation). In a system for which no direct crossing of a neutral
dissociative state with the ground state of the ion is present, this mechanism can become dominant, and it can still produce large DR cross sections. This is expected to beparticularly true for molecules containing hydrogen. A schematic illustration of thesetwo cases is shown in Fig. 3.1.
−
Recently Kokoouline and Greene [95] have demonstrated this occurs in H +3 +ecollisions. For HCO+ , a closed-shell linear ion in its ground state, the picture is similar.
The lowest doubly-degenerate vibrational E states are coupled to the electronic states
E1 (npπ) of the neutral system through the Renner-Teller interaction, resulting in a large
probability for recombination.
Although the Renner-Teller coupling is well known, electron-molecule collisions
where the coupling is present have not yet been studied until recently. This is possibly

42
due to the complexity of the problem because the molecule must be at least triatomic sothat one has to treat three or four degrees of freedom simultaneously. In a recent study[126] McCurdy et al. have made an important step in understanding the role of RennerTeller physics in collisions between the electron and the CO2 molecule. The authors have
employed a time-dependent framework to describe the nuclear motion of the negative2
ion CO−
2 formed during the collision. The nuclei move in the two-component A1 and
2B

1 electronic complex potential corresponding to the doubly-degenerate2Πu molecularstate at linear configuration. The two components are coupled by the Renner-Tellerinteraction.

3.2

Ab initio surfaces
To treat the electronic degrees of freedom and the Rydberg state dynamics wework in the framework of multichannel quantum defect theory (MQDT). This allows
us to treat Rydberg states of arbitrarily high principal quantum number, provided we
can extract a few geometry-dependent quantum defect parameters from accurate abinitio quantum chemistry or scattering calculations. In our model we use a further
simplification, imposing the energy independence of the quantum defects, a statementthat we will show to be accurate, except possibly for the very lowest (n = 2 for HCO + )states, which have considerable valence character.
In MQDT the information regarding the potential energy surfaces of Rydberg
states is expressed in a compact and convenient fashion using an effective K-matrix.
This can be derived from the Lippmann-Schwinger equation [141]111− K=V − VKπ
π E − H0

(3.1)

43

DirectIndirect

v+=20

E
n=20

v+=0

v+=0

n=3

n=2

valence state

valence state

R

Figure 3.1: Schematic representation of direct and indirect DR processes. For directprocesses a crossing of the ion (black curve) ground state and the neutral dissociative
state (in red) leads to quick dissociation. For an indirect mechanism, in an adiabatic
picture, the neutral dissociative state lies below the ionic potential and cannot lead todirect dissociation, since there is no classically-allowed Condon point for capture. The
electron can be captured into a Rydberg state (dashed green curve) converging to an
excited vibrational threshold (full green curve) and, through subsequent decay throughthe Rydberg manifold and additional vibrational excitation, it can find a dissociative
pathway and give a high DR rate.

44
where V is a matrix that represents the electronic interaction between the electron andthe ion, H 0 is the ionic Hamiltonian. We can approach the solution of this equation
perturbatively, and in the case in which K can be considered energy independent, as in
the present case, it is sufficient to stop the perturbative expansion at first order, leadingto1− Ki,j =< ψi | V | ψj > +O(V 2 )
π

(3.2)

where i, j are two different electronic states. This matrix is usually multiplied by -π inorder to agree with the usual reaction matrix of MQDT [95]. Therefore the K-matrix
corresponds to the matrix of diabatic (potential) couplings. It is possible to extract such
quantities by fitting the adiabatic (Born-Oppenheimer) potential energy surfaces of the
Rydberg states to a model K-matrix with off-diagonal couplings. It is important tohave a diabatic description of the electronic states because when a conical (or RennerTeller) intersection is present the different potential energy curves have a true crossing.
In an adiabatic representation the couplings would be infinite, since in a Landau-Zener
picture [101] the couplings are inversely proportional to the difference of the potentials
between the two states. In a diabatic representation the couplings remain finite at every
point.
The electronic dynamics is tackled calculating the Rydberg potential energy surfaces with quantum chemistry codes, and then fitting them to a model K-matrix tobe used in a QDT calculation. Details of the procedure can be found in Sec. 6.2.1.
Some new results for HCO Rydberg states potential energy surfaces, achieved using a
Green function method, can be found in Ref. [49], they compare very well with our
calculations.

45

Figure 3.2: HCO+ ground state ionic and HCO ab initio potential energy curves as afunction of (a) CH stretch at the equilibrium rCO = 2.0877 a.u. and θ = 180 degrees and(b) bending at the equilibrium distances for CO and CH (rCH = 2.0 a.u.). The filled
symbols indicate MRCI potential energy curves of HCO states of A0 symmetry, while
open symbols are states of A” symmetry. The corresponding potential energy curves
obtained using Koopmans’ theorem are shown as solid and dashed lines respectively.
From Ref. [104].

46
We have performed two calculations of DR in HCO+ , which we will call Model 1
and 2, that are described in detail in this chapter. The results are shown in Ch. 6. The
main difference between the two models is the treatment of the Renner-Teller couplings,
which are not included in Model 1, whereas they are present in Model 2.
We assume that the incident electron is initially captured into one of the neutral
Rydberg states, which can be approximately characterized by the principal quantum
number n, the orbital angular momentum l, and its projection on the molecular axis
λ. In our model we have included only the npπ −1 ,npπ +1 , npσ, nsσ and ndσ states.
The symbols ±1 imply two opposite-sense electronic angular momenta associated withthe different signs of λ. Linear combinations of the nπ ±1 states give states npπ 0 and
npπ 00 , symmetric and antisymmetric with respect to reflection in the plane containing
the molecular axis. The other three states are symmetric. These states were established
in Ref. [104] to be the most important for DR, since they exhibit the largest dependence
on the vibrational coordinates.
In the first model we have fitted our ab initio electronic surfaces to a K-matrix
having the following structure: Ksσ,sσ Kpσ,sσK(Q) =  Kdσ,sσ




0




0

Ksσ,pσKsσ,dσ00Kpσ,pσ Kpσ,dσ

00

Kdσ,pσ Kdσ,dσ

0

0

0

0

Kπ−1 ,π−1

0

0

0

0

Kπ+1 ,π+1

where Q represents all the vibrational coordinates.











,









(3.3)

The important characteristic of this model K-matrix is the absence of RennerTeller coupling, which implies that the π-block is diagonal, while the sigma block is fully

47
coupled. Renner-Teller coupling was initially omitted because of the observed nature of
the Rydberg states surfaces as a function of the bending angle, shown in Fig. 3.2. In
particular it appeared to suggest that the Renner-Teller splitting was large for n = 2
states, but very small for all other Rydberg states. The higher Rydberg states are
responsible for the capture of the electron in the DR process, therefore Renner-Teller
coupling was thought to be uninfluential for the reaction rate. The fitting procedure is
described in Sec. 6.2.1.

3.3

Frame transformation and adiabatic methods
The frame transformation technique, which has similarities to the adiabatic nu-

clei approximation (also called Chase [31] or infinite-order sudden approximation [147])
assumes that the timescales for electronic motion and vibrational/rotational motion arevery different, as in the Born-Oppenheimer approximation. This would seem at first
not to hold for very high Rydberg states, where the Rydberg electron is so slow that its
motion has a period similar to a molecular vibration or longer. In fact, this approximation holds up very well for Coulomb systems no matter how high the principal quantum
number n is. This is because the Rydberg state has an elliptic “orbit”, therefore evenfor very high n the electron spends a large amount of time outside the ion core, andvery little time inside. Effectively the speed with which the electron passes throughthe core is almost energy independent, due to the large Coulomb acceleration that theelectron receives close to the ionic core, which is always much larger than the Rydberg
state energy. For neutral systems this does not happen, and indeed close to thresholdthe frame transformation treatment is less accurate. See, however, Ref. [76] for e-N 2
scattering.

48
We treat the nuclear dynamics in a two-step adiabatic framework, using Jacobi
coordinates. This two step process is dictated by the need to calculate many vibrational
states for the ion, which a full three dimensional calculation (which we have carried out
anyway to check our results) would not have allowed. The coordinates are chosen to
mimic as much as possible the normal modes of the molecule, since this improves theadiabaticity of the curves, probably because the couplings between the modes are minimal. We had initially tried to represent this ion using hyperspherical coordinates, but
the simple adiabatic picture failed completely. In H+
3 [95] the hyperspherical coordinates are similar to the normal modes and the adiabatic approximation holds well. The
variable R represents the CH distance, while R is the distance between the oxygen atomand the center of mass of the other two, which is very close to the CO distance. theangle θ is the acute angle between these two vectors:
R = rCHr=

mC rCO + mH rOH.
mC + m H(3.4)

The Hamiltonian in Jacobi coordinates for J 6= 0 reads [188]
HV = −∂h̄ ∂ 2 ∂
h̄ 1
1
1 ∂
∂
h̄ ∂
(R2
)−
(r
)− (
+
)
(sin θ ) (3.5)
2222µ1 R ∂R∂Rµ2 r ∂r∂r
2 µ1 Rµ2 r sin θ ∂θ∂θ

1 1
cot2 θ csc2 θ 2 cot θ2
2
HV R = [
(J
+
J
)
+
(
+
)J +
(Jx Jz + Jz Jx )xy
2 µ1 R 2µ1 R 2µ2 r 2 x µ1 R 2
h̄ 1
cot θ
∂
+
+
)Jy ]
(
2
i µ1 R ∂θ
2(3.6)

where Jx , Jy , Jz are the Cartesian components of the total angular momentum operator.
The need to include J 6= 0 states, even in the absence of rotational dynamics, is dueto the fact that for a linear molecule the odd quanta of the bending vibrations are49
present only for odd J. This means that considering only states arising from a J = 0
calculation would lead to neglecting the lowest vibrational excited state (01 1 0), thequantum numbers are illustrated in Tab. 3.1. In the calculations we include only J = 0
and J = 1 and neglect the rotational spacing.
We use a spline basis [11] to carry out the target vibrational calculations in all
three dimensions, first solving the two dimensional r − θ problem and then carryingout the R-dependent calculation on the adiabatic curves obtained in the first stage.
A schematic illustration of the molecule with the Jacobi coordinates we use and theadiabatic vibrational curves as a function of rCH is given in Fig. 3.3.
The vibronic K-matrix is obtained through a two-step frame transformation. The
first step consists in taking the overlap of the electronic K-matrix with the adiabatic
CO-bending wavefunctions
DE

0 ,ν 0.
KΛνr νϑ ,Λ0 νr0 νϑ0 (RCH ) = φνr ,νϑ |KΛ,Λ0 (Q)|φν∇
ϑ

(3.7)

In practice this integral is done on the quantum defect matrix, µ = π −1 tan−1 K, since
the K-matrix has poles.

3.4

Siegert states
In this first study, we have used Siegert (outgoing-wave) boundary conditions[80, 189, 190] in the C-H degree of freedom. This allows the description of outgoing
dissociative flux, i.e. predissociation. If outgoing wave boundary conditions are imposed at a finite distance, which is possible for finite range potentials (or for long-range
potentials as suggested in Ref. [162]), one can use an L2 (square-integrable) basis set torepresent the Siegert states. Using these boundary conditions, a quadratic eigenvalue

50

O
H
φ

r

R
C
Adiabatic potential energy curves for HCO+
0.18
0.16
0.14

Energy (a.u.)

0.12
0.1
0.08
0.06
0.04
0.02
0
1.4

1.6

1.8

2

2.2

2.4
rCH (a.u.)

2.6

2.8

3

3.2

3.4

Figure 3.3: Top: Jacobi coordinates used to model HCO+ vibrational dynamics, R isthe adiabatic coordinate and it is coincident with rCH . Bottom: adiabatic vibrational
curves as functions of the rCH adiabatic coordinate. The figure represents the 40 lowest
eigenvalues of the adiabatic Hamiltonian.

51
problem in k, the wave vector(H − ikL − k 2 O)~c = 0(3.8)

is recovered, where L is the matrix of the Bloch operator, and O is the overlap of thebasis functions. Eq. 3.8 can be solved iteratively. It can also be dealt with in one stepby solving an eigenvalue problem of dimension twice the dimension of the basis space[190]. The second step in the frame transformation is then to use the Siegert states(φνR ) to calculate the full S-matrix
E

DSΛνr νϑ νR ,Λ0 νr0 νϑ0 νR0 = φνR |SΛνr νϑ ,Λ0 νr0 νϑ0 |φνR0 .

(3.9)

The channel elimination procedure [74], in which physical (exponentially decaying)boundary conditions are imposed in the closed channels is then applied. If we arrange
the full S-matrix in open and closed channel partitions, we have [74]:
S phys = Soo − Soc (Scc − e−2iβ )−1 Sco(3.10)

and we obtain the physical S-matrix (S phys ) of dimension No × No where No is thenumber of open channels [74]. The matrix β in Eq. 3.10 is diagonal with elementsqβj = π/ 2(E − Ej ) and it represents the Coulomb long-range phase parameter; E isthe total energy, and Ej is the threshold of the j−th channel. When the Siegert states
are used, the physical S-matrix will not be unitary and the difference from unitarity will
be related to the DR cross sectionπσ(Eel ) =
2E

Ã

1−

No
X
i=1

†phys
phys
(E)
Si,i
0 (E)Si0 ,i!

,(3.11)

where Eel = E − Ei . This nonunitarity stems from the fact that the outgoing flux,
which will be lost, is assumed to go exclusively into dissociative channels, which are52
not included explicitly in the S-matrix. Notice that, following the definition of theSiegert states, their orthogonality relation involves a surface term, and therefore theframe transformation performed in Eq. 3.9 will have a similar surface term as in Eq. 6
of Ref. [80]. The results for this model are shown in Sec. 6.2.

3.5

Model Hamiltonian including the Renner-Teller coupling
Since the results obtained with the previous model K-matrix have not been sat-

isfactory, we have built a new treatment that has as a focal point the inclusion of theRenner-Teller coupling. We also freeze the CO degree of freedom, and we include many
vibrational angular momenta, instead of the two included in our previous treatment.
The starting point of the present theoretical approach is the molecular Hamiltonian H
of HCO, which we represent as H = Hion + Te + Hint , where Hion is the ionic Hamiltonian, Te is the kinetic energy of the incident electron, and Hint describes the electron-ion
interactions. We assume that the ion is in its ground electronic state. Hion and Hint
depend on the four internuclear coordinates Q = {RCH , RCO , θ, ϕ}, where RCH and
RCO are the C-H and C-O internuclear distances, θ is the bending angle, which is zerofor linear configurations. The electronic energies are independent of the angle ϕ, which
represents the azimuthal orientation of the bending.
As in Model 1 we consider 5 electronic states. Therefore, Hint + Te assumes
a block-diagonal form with an infinite number of 5 × 5 blocks corresponding to n =2, 3, · · · , ∞, in addition to the continuum. In the derivation of the model we assume
that there are no couplings between manifolds with different principal quantum numbers;
later, the MQDT treatment will relax this constraint.
Since the only off-diagonal couplings included in our analysis are those among

53
the three np states, for simplicity we specify only these states in the formulas below.
The np sub-block of Hint + Te in the basis of the npπ −1 , npσ and npπ +1 states
has the form [98]


δeiϕ

γe2iϕ

Eσ

δeiϕ

γe−2iϕ δe−iϕ

Eπ

 Eπ



Hint (Q) =  δe−iϕ










,




(3.12)

where Eσ and Eπ are the electronic energies of the npσ and npπ ±1 states at the linear
ionic configuration; δ and γ are the real, non-Born-Oppenheimer coupling elements. We
denote both the (nπ − nπ) and (nπ − nσ) couplings as Renner-Teller (RT) couplings,
whereas in some previous studies, only the former is denoted by this term. The couplings δ and γ depend on RCH , RCO , and θ. The diagonalization of [Te + Hint ](Q) is
accomplished by the unitary transformation matrix U [98]:




1 

U=√ 
2



eiϕ

eiϕ w−

0√2w+−e−iϕ e−iϕ w−

with the abbreviations

w± =

q

w =

p



eiϕ w+ √,− 2w− e−iϕ w+(3.13)



(1 ± ∆/w)/2; ∆ = (Eσ − Eπ − γ)/2;∆2 + 2δ 2 .

(3.14)

When diagonalized, the Hamiltonian becomesU † [Hint + Te ]U = diag{Vπ00 , Vσ , Vπ0 } ,(3.15)

where
Vπ00= Eπ − γ ,Vπ 0= (Eσ + Eπ + γ)/2 − w.

Vσ = (Eσ + Eπ + γ)/2 + w ,(3.16)

54
The adiabatic potential energy surfaces Vπ0 ,π00 ,σ (Q) are known from ab initio calculations(see Ref. [104] for a detailed description). From Eq. (3.16) we obtain
γ = Eπ − Vπ00 ,w = (Vσ − Vπ0 )/2 .

(3.17)

Therefore, the matrices U in Eq. (3.13) and Hint in Eq. (3.12) are obtained from
Vπ0 ,π00 ,σ (Q).

3.6

Reaction matrix and O’Malley formalism for the DR crosssectionAgain here we employ multichannel quantum defect theory, and the first step is toconstruct the reaction matrix K. First we introduce the diagonal quantum defect matrixµ, whose nonvanishing elements µi are related to the diagonal element of the potential
surface by Vi = −1/[2(n − µi )2 ], where i is the electronic state index. Therefore, thediagonal form of the reaction matrix is directly obtained from the ab initio calculation.
The coupling δ in Eq. 3.12 is linear in θ for small θ. Correspondingly, thesplitting between adiabatic electronic energies Vπ00 and Vπ0 is quadratic with respect toθ, and, therefore, the splitting between quantum defects µπ00 and µπ0 is also quadratic.
This can be seen in Fig. 3.4, where µπ00 is almost constant because it has a different
symmetry and is unchanged due to the Renner-Teller effect. The quantum defect µ π0
is approximately quadratic due to the Renner-Teller effect until θ ≈ 50o , beyond which
the linear approximation for δ(θ) is not adequate. In a preliminary calculation we usedsuch a linear dependence of δ(θ), because we have found that we can use Eqs. 3.13, 3.14,
3.15, 3.16, and 3.17 and fit the model Hamiltonian accurately to the ab initio energies
without the additional approximation of linearity of δ(θ). The result of the preliminary

55

0.6

pσ
dσ
sσ −
pπ
+
pπ

Quantum defects

0.4

0.2

0

-0.2

-0.4

0

20

40

o

60

Bending angle θ ( )80

Figure 3.4: The eigenquantum defects of the electronic states included in the present
study as functions of the bending angle θ for fixed RCH = 2.0 a.u. and RCO = 2.0877
a.u. Labeling the states with electronic momentum projection λ is an approximationand for large values of θ it is not appropriate. However, the four sσ, pσ, dσ, and pπ 0states are always uncoupled from pπ 00 states for nonlinear geometries. In the absenceof the non-diagonal matrix elements in Eq. 3.12 the pπ 0 and pπ 00 states would have thesame quantum defects. The Renner-Teller parameters are completely determined by
the geometry dependence of pπ and pσ defects (Eq. 3.17). From Ref. [130].

56
calculation is consistent with the present calculation using the accurate fit to the abinitio energies.
The next step is to bring the reaction matrix into the representation of the Hamiltonian of Eq. 3.12. Since the transformation matrix U is known from Eq. 3.15, the K
matrix in that representation reads
K = U tan(πµ) U † .

(3.18)The matrices µ̂, Hint + Te , K, U are diagonal for the nsσ and ndσ states. We used
quantum defects from Ref. [104] with n = 4 for nsσ and n = 3 for the other states.
Figure 3.4 gives the quantum defects as functions of θ, with RCH and RCO fixed at theionic equilibrium values.
Once the reaction matrix Ki,i0 is obtained, the DR treatment is along the lines
of Refs. [94] with the following differences. We choose here the RCH distance as thedissociative adiabatic coordinate. This assumption is justified because experimentally
the H+CO channel is largely dominant at low energies. It has been long known that the
CH+O channel is endothermic by 0.17 eV and could have an energy barrier, while the
OH product is not observed (see Ref. [144] and references therein). Recently, branchingratios have been measured [57] for the DCO+ dissociative recombination, confirmingthat the D+CO channel is by far the dominant, with a dissociation branching ratio of
0.88. We keep the RCO distance frozen at its ionic equilibrium value. This assumption
was made to simplify the treatment but it could result in an underestimation of the DR
cross section. In order to account for the CO vibration one might use the hypersphericalcoordinates similar to Ref. [104], or else normal mode coordinates as in McCurdy et al.
[126]. Although the CO bond is frozen, the Renner-Teller coupling physics is included

57

Ion energies (eV)

0.6

0.5

0.4

0.3

mϕ=0
mϕ=±1
mϕ=±2
mϕ=±3

0.2

0.1

1.8

2

2.2

RCH (a.u.)

2.4

Figure 3.5: Several adiabatic potential curves for the HCO+ ion are shown as functionsof the RCH distance for different projections mϕ of the vibrational angular momentum.
The curves of the same color correspond to the vibrational states with the same quantum
number mϕ but different L. The energies of the lowest vibrational states with mϕ = 0
and mϕ = ±1 are shown by two horizontal lines. From Ref. [130].

58
in our model; we show below that it does increase the DR cross section significantly.
We will briefly address below the possible influence of the CO vibration.
Consequently, for every RCH distance we calculate matrix elements of the “large”
reaction matrix K(RCH )Kj,j 0 (RCH ) = hΦmϕ ,L |Ki,i0 (Q)|Φm0ϕ ,L0 iθ,ϕ ,(3.19)

where Φmϕ ,L (RCH ; θ, ϕ) are vibrational wavefunctions of HCO+ parametrically dependent on RCH , while RCO is frozen. They are eigenfunctions of Hion with eigenvalues+
Um(RCH ). Each index j in K can be represented as {i, mϕ , L}, where mϕ specifies
ϕ ,L

the projection of the vibrational angular momentum on the molecular axis, while theindex L distinguishes different vibrational states with the same mϕ . The rotation of the+whole molecule is not considered. Fig. 3.5 shows the vibrational eigenenergies U m
asϕ ,La function of RCH .
It is informative to compare the ionic vibrational energies obtained in our simplified approach with the exact calculation accounting for all four vibrational coordinates[154]. We cannot compare the absolute values of the energies because the CO stretch
is frozen, however, we can compare the energies of a few vibrational levels {v 1 v2l , v3 }with respect to the energy E = 0 of the ground vibrational level {000 0}. Table 3.1
compares the vibrational energies obtained in this study with accurate results from Ref.
[154]. The overall errors of the present calculated vibrational levels are about 12 %or less. Since the vibrational wave function error can thus be deduced to be of order√

0.12 ≈ 0.35, whereby a conservative estimate of the error in our calculations of the

final cross section would be of order 70%.
The matrix K(RCH ) is used to obtain the potential curves Ua (RCH ) of HCO

59{v1 v2l , v3 }
100 0
011 0
020 0
031 0
040 0

Present calculation
2933
740
1460
2200
2925

Puzzarini et al.[154]3090.40830.7
1641.142458.93256.94

Table 3.1: The table compares the accuracy of the vibrational energies obtained in thepresent study that neglects CO vibrations with the “exact” calculation from Ref. [154].
v1 indicates CH stretching, v2 is the bending quantum number, v3 is CO stretching andl the is vibrational angular momentum. Since the CO bond is frozen, we don’t provide
energies for excited v3 modes. The overall energy error is about 12 %, which translates
into an error of about 35 % for the vibrational wave functions. The energies are givenin cm−1 . From Ref. [130].

60
as described in Ref. [94]. The curves in general have non-zero autoionization widths
Γa (RCH ).
The potential curves Ua (RCH ) and their autoionization widths Γa (RCH ) are thenused to calculate the DR cross section. The cross section is calculated in a manner
similar to the procedure described in Refs. [94], with some adaptations. Specifically,
we start from Eq. (5.19) of Ref. [142] for the cross section of the process of dissociative
attachment of the electron to a neutral molecule, which applies equally to the dissociative recombination process considered here. The resulting approximation adapted tothis situation reads (see Eq.6 of Ref.[94]):
σDA =2π 2 Γa (RCH ) +|χ (RCH )|2 e−ρ(E) .
ko2 |U 0 (RCH )| o

(3.20)

Here e−ρ(E) ≤ 1 describes the survival probability, which could be less than unity if
the system has a substantial probability to autoionize before it dissociates; k o is theasymptotic wave number of the electron incident on vibrational level o of the target
molecule. The distance RCH is a dissociation coordinate understood to be evaluatedat the Condon point, which depends on the total energy as well as the initial target
vibrational level o; χ+o (RCH ) is the initial vibrational wave function of the ion. The
above formula is appropriate for the capture into a resonant potential curve that isenergetically open for direct adiabatic dissociation. If the corresponding resonant stateis closed, i.e. bound with respect to dissociation, it requires modification. To this end,
we adapt Eq. (4.2) of Bardsley [16] to our present situation involving indirect DR. The
indirect process proceeds via capture into a bound (typically Rydberg) state, which
eventually predissociates. For a case involving a single incident electron partial wave
(e.g. pπ 0 , for definiteness) and a single ionic target state, the fixed-nuclei autoionization

61
width of the resonance potential curve will be denoted Γa (RCH ). Once the vibrationalmotion in the resonance potential is quantized into a vibrational resonance level, withradial wavefunctions χres (RCH ), it should be remembered that only a subset o of thetarget vibrational levels v will be energetically open, at any given total energy E. Each
resulting quantized resonance acquires a partial autoionization width Γ a,o0 into an open
vibrational channel o0 , which is given approximately by
Γa,o0 = |

Z ∞
0

q

res
2
χ+
o0 (RCH ) Γa (RCH )χ (RCH )dRCH | .

(3.21)The sum of these is then the total resonance autoionization width (within this approximation, neglecting Rydberg level perturbations of the MQDT type, which could
sometimes produce complex, non-isolated resonances),Γa,tot =X

Γa,o0(3.22)

o0

while the total linewidth Γ of this quantized resonance also includes its predissociation
partial linewidth, Γd : Γ = Γd + Γa,tot . In this notation, the contribution to the total DR
cross section will be the following, if only the ground vibrational level o of the target ispopulated:
σ res (E) =Γa,o Γd
2π 2 1.
ko2 2π (E − Eres )2 + 41 Γ2(3.23)

Our model does not account for vibrational motion along the CO bond and,
therefore, positions of resonances cannot be compared directly with the experiment.
Moreover, the resolution in the available experiments is usually not sufficient to resolve
individual resonances. This suggests that we should average the cross section overenergy. Each energetically-closed ionic bound vibrational state c generates a Rydberg
series of resonances ²nc , numbered with the effective principal quantum number nc . The

62cross section averaged over the energy interval between two resonances with energies
E = ²n and ²n+1 is given by
hσ(E)i =1∆n + ∆n+1

Z E+∆n+1

σ(E 0 )dE 0 ,E−∆n∆n+1 =²n+1 − E.
2(3.24)

Extending the limits of the integral to infinity, one can easily obtain the following:hσi ≈2π 2 Γa,o Γd
(
)
ko2
Γ∆

(3.25)

where the ∆ = ∆n + ∆n+1 ≈ 1/n3c . Note that the quantity in parentheses here approaches a constant at sufficiently high effective principal quantum numbers nc in therelevant closed channel c, because each partial and total width in the parentheses of thisformula should become proportional to ∆ in this limit. The total cross section is thencalculated by summing up the average contribution from all Rydberg states (a sum overclosed ionic channels). In the limit where Γd >> Γa,tot , this gives
hσi =q
2π 2 X res
2 3
|hχ
(R
)|
Γa (RCH )|χ+CHo (RCH )i| nc .
ko2 c(3.26)The projection M = mϕ + λ of the total angular momentum on the molecularaxis O–C is a conserved quantity in our model. We calculate all the resonances and thecross section for all allowed values of M . Since we are considering only sσ, pσ, dσ, pπ 0 , π 00
quantum defects and the initial ion is in its ground state mϕ = 0, the only possible
values of M for the final state are 0 and ±1. Resonance energies and widths are ofcourse dependent only on |M |. Finally, we obtain for the cross section:
hσ total i = hσ M =0 i + 2hσ M =1 i .
The results from this model are shown in Sec. 6.3.

(3.27)

Chapter 4

Results: electron scattering, small molecules

4.1

Introduction
The approach developed in sections 2.2-2.3 is now tested in calculations of elec-tron scattering by N2 and CO2 , classic benchmarks in this field [120, 139, 156], because
their elastic cross sections exhibit striking features that can be challenging to reproduce. These first results were published in Ref. [192]. We have moved to larger target
molecules, specifically SF6 , XeF6 and C60 fullerene (which is not presented here), tounderstand how well our approach scales for these large systems. These results are in
sections 4.3.2-4.4.2. Photoionization calculations are described in Sec. 4.4.

4.2

Computational Details
Three-dimensional integrals are expensive to calculate in general so it is highly

desirable to minimize the time spent in their calculation. For the sectors that do notcontain a nucleus it is possible to use just 4 Gauss-Legendre points of integration in eachcoordinate, since doubling the number of points changes the calculated eigenphase shifts
by only about 10−6 radians, while increasing the computation time by approximately anorder of magnitude. Particular caution has to be observed when integrating over sectors

64that contain a nucleus. We have found it important in general to have a finite element
vertex on each Coulomb singularity, in order to obtain correct results, and to use more
integration points. In these sectors we use 20 integration points in each dimension since
we found that the convergence of the phase shifts in this case is, as in the previous case,about 10−6 . It is important also to provide sector boundaries near the K-shell radius foreach atom in the molecule, since close to the nucleus the wavefunction varies rapidly. A
fine grid is needed to represent the wavefunction oscillations near the nuclei, otherwise
the phase shifts for the lowest, penetrating partial waves will be unconverged.
The sparse structure of the finite element matrices (see Fig. 4.1) can be exploitedwith great advantage from the beginning. No matter how the grid is defined, each basis
function has matrix elements with at most 216 functions. This allows us to know thesparsity pattern of the matrix Γ in Eq. 2.8 (which we show again here for convenience)~ = 2(EO − H − L)C~~ = ΛCbΓC(4.1)

in advance and store just the nonzero elements, with a reduction in memory cost of
approximately two orders of magnitude. This economy is crucial to allow us to perform
three dimensional calculations in the first place.
The dimension N of the eigensystem in Eq. 2.8 is, for CO2 , of the order of40000, whereas for the open subspace it is only 100 or less. N increases rapidly withthe complexity and spatial extension of the molecular potential, but the sparsity of thematrices is high (about 0.5% full for N ∼ 40000 ), and it increases with the dimension of
the system. To solve the linear system in Eq. 2.11 we use direct sparse LU factorization
solvers (SuperLU [87] at the beginning, now more extensively the Pardiso [168] solver).
We have tried in the past to use iterative methods, such as biconjugate gradient [68]

65
, with different preconditioners to speed up the solution of the linear system. The
one we have found to work best has been an incomplete Choleski factorization, which
reduces drastically the number of iterations with respect to a diagonal preconditioner,the Γ matrix in Eq. 2.11 is not, in fact, diagonally dominant. Iterative methods are
slower than direct factorization: in the tests we have performed normally the direct
method is faster by an order of magnitude, but they are convenient if there are memory
limitations. Since the factorization of a sparse matrix does not preserve the sparsity
pattern, the factorized matrices present storage problems, since a fill-in factor of around
10-15 is common for these systems. We have decided to not use iterative methods since
our linear system is often ill-conditioned and in these cases reaching convergence is very
difficult, while direct methods always achieve the correct solution. This code therefore
needs a large amount of memory for large target calculations, for which we have used
Opteron 64-bit workstations with great success.

4.3

Results
In the first few calculations we show here, also published in Ref. [192], we haveused an adjustable parameter in the polarization potential of Eq. 2.20, which here isused both at short and long distances:
Vpol = −r 6
1[α0 + α2 P2 (cos θ)][1 − e−( rc ) ].
4
2r

(4.2)

It is appropriate to carry out the expansion up to these first two terms, and neglect
the rest of the dipole polarizability tensor, since the target molecules for which thispotential is used are either linear or approximately so. Here rc is a distance parametercomparable to the range of the target charge distribution. When high accuracy is

66

0
100
200
300
400
500
600
700
800
900
0

100

200

300

400

500

600

700

800

900

Figure 4.1: Structure of the finite element matrix Γ for a small test case of dimension
900. Notice the great sparsity of the matrix, which increases with the dimension of thematrix. From Ref. [192].

67
needed for resonance positions in some applications, rc can be determined empirically[139] to reproduce the energies of one or more resonances of interest. In the calculations
for N2 , CO2 and HeH we used this simple adjustable potential instead of the moresophisticated and ab initio DFT potential described in section 2.2.3. This allows atuning of the resonant features, if needed for particular purposes, but it is detrimentalof the predictive power of the approach. In the calculations of Sec. 4.3.2 and following,
we use the DFT polarization potential.
The strong and narrow Πg resonance at 2.4 eV in N2 is reproduced in our calculations at the right energy, provided we use a physically reasonable cutoff radius r c = 2.8
a.u. for the polarization potential. In Fig. 4.2 are shown the results. Other theoreticalstudies of this resonance, definitely one of the most studied in electron-molecule scattering, can be found in Refs. [85, 137, 202], where exact treatments of exchange are
performed. The resonance is reproduced also at the static exchange level (without using
a long range polarization potential), but at an energy of 5.4 eV.
For CO2 the main feature in the total elastic cross section is a Πu resonance at 3.8
eV. To reproduce it at the correct energy we have to tune the polarization cutoff radius
to 2.4 a.u.. This feature is present also at the static exchange level, at 8 eV. It should benoticed that in calculations made using an exact form of static exchange this resonancelies at 6 eV [63, 118, 134, 156]. The polarization dependence becomes much more pronounced than in N2 , as one expects from the larger spatial extension, the larger numberof electrons and the greater asymmetry of this molecule. The scattering cross section
for this system is shown in Fig. 4.3 . Fig. 2.2 shows how the polarization potential of
N2 is appreciable just outside the region where the main part of the electronic density is
located. Therefore the value of the cutoff radius for the polarizability potential, which

68
is the only adjustable parameter in the model, is reasonable. The results are always ingood agreement with previous theory, as shown in the figures.
The present calculations have been performed for the molecular targets only attheir equilibrium distances. Vibrational effects tend to broaden these resonances inexperimental elastic scattering cross sections, and they also give rise to more structured
resonance peaks, which are not considered in this work.
The values of the polarizabilities used in these calculations are [139] α0 = 11.89a30
and α2 = 4.19a30 for N2 and [102] α0 = 17.9a30 , α2 = 9.19a30 for CO2 . It should
be pointed out that accurate static polarizability coefficients α0 and α2 in Eq. 2.20
can also be extracted from ab initio calculations, although to get the correct values,
extremely diffuse basis sets have to be used. Generally, the low-lying shape resonancespresent in these small molecules are spatially highly localized, which allows the radiusof the R-matrix box to be kept small, around 8 to 14 a.u. for the present calculations.
In the case of the third neutral molecule that we present here, ethylene, thesituation is more complicated. Since the target is now nonlinear it is more difficult to
describe in a discrete basis set and it is more expensive computationally to calculate
the scattering cross section. Nevertheless we are able to reproduce the features of theelastic cross section for this molecule. We find good agreement with the energies of theresonances and with the overall cross section magnitude, compared with previous theory
and also experimental data, although the vibrational effects again tend to broaden theresonance peak.
It should be noticed parenthetically that if we neglect exchange altogether incalculations for all of the molecules presented here, the cross sections are qualitatively
wrong, with resonances far lower in energy than the experimental ones and in the wrong

69
symmetry channels. This is due to the fact that some of the target electrons are notbound anymore in the resulting potential without exchange, because the static potential
is not attractive enough. For example, in electron scattering from N 2 the Πg resonanceat 2.4 eV is recovered at 5.4 eV for a static exchange calculation in our model potential.
Using only the static potential, there will be a Σg resonance around 2 eV, which is dueto the σu electrons of the target being above the ionization threshold. The exchange
potential is basically an attractive potential (local in our model), so adding it to thestatic potential produces the correct number of bound states for the target; consequently
the scattering resonances are generated by capture of the electron in truly unoccupied
molecular orbitals of the target. A more systematic study of the behavior of the crosssections, when different parts of the potential are neglected altogether, can be found in
Ref. [102].

4.3.1

Quantum defect calculations

It has been shown [186] that use of a local density approximation can often be
effective in calculating molecular quantum defects, for bound or scattering states, forsmall closed-shell target molecules. It is possible to calculate quantum defects from ascattering calculation carried out near zero energy. The key step is to diagonalize theK-matrix
Kii0 =XT
Uiα tan πµα Uαi
0

(4.3)

α

and then utilize the relationship between the quantum defect and the scattering phaseshift [74, 169],
δα = πµα .

(4.4)

70

50
Present results
Morrison and Collins
experimental

45

40

σ (Å 2)

35

30

25

20

15

10

5
0

1

2

3

4

5

6

Energy (eV)Figure 4.2: The integrated elastic cross section versus energy for electron-N 2 scattering, is compared to the theoretical results of Morrison and Collins [139] and to theexperimental cross section of Kennerly [90]. From Ref. [192].

71

40
Present results
Rescigno and McCurdy
Morrison and Collins
experimental

35

30

2

σ(Å )

25

20

15

10

5

0
0

1

2

3

4

5

6

7

8

9

10

Energy (eV)Figure 4.3: The integrated elastic cross section for scattering of electrons from CO 2 isshown as a function of the incident electron energy. The present results are comparedwith previous theory from Rescigno et al. [156] and Morrison and Lane [102], whereas
the experimental results are those of Szmytkowski et al. [182]. From Ref. [192].

72

Elastic cross section for e−CO2 (Å2)

35

present calculation
Rescigno SEP
Morrison SEP
Rescigno SE

30
25
20
15
10
5
0
0

5

10

15

20

25

Energy (eV)Figure 4.4: The e-CO2 elastic cross section at low energies using our completely abinitio potential. The prominent structure is the well-known πu resonance which, inexperiment, has been found at 3.8 eV. Often the resonances we calculate with ourapproach are roughly 2-2.5 eV higher if compared to Fig. 4 in Ref. [192]. Our results
are compared with those of Rescigno et al. [156] (exact static exchange and exchange
plus polarization) and Morrison and Collins [139] (model exchange plus semiempiricalpolarization). The resonance width is comparable to the SEP results of Rescigno et al..
From Ref. [191].

73

5.5

5.4 eV
10.9eV21.8 eV
32.6 eV

Differential cross section for e−CO2 (Å)2

5
4.5
4
3.5
3
2.5
2
1.5
1
0.5
0
0

20

40

60

80
100
angle (degrees)

120

140

160

180

Figure 4.5: The e-CO2 differential cross section is shown as a function of scattering
angle at four low electron impact energies. The results can be compared with those ofRefs. [63, 64, 138, 156]. From Ref. [191]. We calculate this differential cross section,
averaged over molecular orientations, using the formalism derived in Ref. [41].

74

35

30

25

2

σ(Å )

20

15

10

Present results
Schneider
Sueoka and Mori
Winstead and McKoy
Panajotovic

5

0
0

2

4

6

8

10

12

14

16

Energy (eV)Figure 4.6: The integrated elastic cross section for electron-C 2 H4 scattering shown as afunction of energy, is compared to previous theoretical results of Winstead et. al. [205]and of Schneider et. al.. [167] The experimental results are the ones of Panajotovic
et. al. [145] and of Sueoka and Mori [179]. From Ref. [192].

75
Accordingly, quantum defects can be extracted from electron-scattering calculations at
positive or negative energies if energy-dependence of the parameters is investigated.
These quantum defects can then be used to determine the Born-Oppenheimer potential
energy curves of the Rydberg states converging to the various ionization thresholds
through the Rydberg formula, Eq. 2.35. These potential energy curves can then be
exploited through MQDT techniques, to extract dynamical information to model, forexample, dissociative recombination [95], as we show in Ch. 3. Here we show an exampleof how well this approach works for a simple diatomic molecule.
We compare our results to the work of Sarpal and Tennyson [163] which made no
approximation about the nature of the electron-molecule potential, and the agreement
is generally very good. The quantum defects represented in Fig. 4.7 are the mostimportant ones, partial waves (l > 2) having very small phase shifts at the low energies
considered here. In electron scattering from an ionic target we must account for thefact that heteronuclear target molecules like HeH+ have a dipole moment, so we must
transform from the center of mass frame to a new frame defined by the center of charge(the proton in this case), which allows us to match to simple Coulomb functions at theboundary of the box. Otherwise multipole potentials have to be included in the external
region.

4.3.2

XeF6

The xenon fluorides are among the few stable molecules formed by this noble gas.
XeF6 has a structure that can be seen as a distorted octahedron [28]. The distortionis commonly ascribed to the lone pair of electrons that remains on the xenon atom,
after it has formed 6 bonds with the fluorine atoms. Since this system includes a heavy

76

0.2
0.1
0

µ (modulo 1)

-0.1
-0.2
-0.3
-0.4
-0.5

p-Σ
p-Π
s-Σ
p-Σ Sarpal and Tennysonp-Π Sarpal and Tennysons-Σ Sarpal and Tennyson

-0.6
-0.7
1

1.5

2

2.5

3

3.5

4

R (a.u.)Figure 4.7: Our calculated eigenquantum defects at an energy of 0.1 eV for the HeH
molecule calculated with our method are compared with those of Sarpal and Tennyson[163] as a function of the internuclear distance R. From Ref. [192].

77
atom, it is hard to treat with more sophisticated methods because of the large number ofelectrons and consequently the large basis set requirements. Our approach instead scales
roughly the same way as for SF6 . Notice that relativistic effects could be significant in
Xe, see for example Ref. [28], but we neglect them altogether.
For convenience, and because there are no experimental data to compare to forthis system, we have distorted the molecule back to octahedral geometry. The elastic
cross section results are shown in Fig. 4.8. The resonant structures are similar to SF 6 ,since the two molecules have similar electronic structure and, in this case at least, thesame geometry.
The distortion of the octahedral cage could lead to interesting effects on theresonances that might be worth exploring in the future. Many highly coordinatedcompounds show “cage” effects [39] in their electronic structures and interactions with
continuum electrons. Molecules like C60 , which have an internal cavity and therefore apotential energy barrier in the motion of an electron between the internal and the external regions, can trap electrons in this kind of resonance. A detailed study could uncoverinteresting properties. In cluster environments, caging effects have been measured forheavy particles, for example in photoionization of negative ions in a CO2 cluster [146].
Likewise these effects could be pronounced for an electron as well, leading to a sensible
modification of single-molecule shape resonances. Through the trapping of an electron
inside the cage, it might be possible to extend the lifetime of the diffuse anionic states,
while the continuum electron scatters off the cluster cage. These phenomena will bereserved for future studies using the same approaches adopted here for single molecules.
A spatial map of the wavefunction, calculated as described in Sec. 5.1.3, for thelowest energy resonance is shown in Fig. 4.9, as a cut through the x − z plane of the

78molecule, which shows clearly the positions of the Xe nucleus at the center and the four
fluorine atoms in the plane around it, with the areas of greater electron density shownin blue and yellow.

4.4

Photoionization
In photoionization studies, it is possible to detect resonances of the ion that

couple to different initial states by ionizing different electrons in the molecule, somethingthat cannot be done in scattering experiments. We have studied two molecules inparticular, CO2 and SF6 , described in detail below. All the cross sections are multipliedby the number of degenerate electrons of that symmetry, as described in Sec. 2.3.
Here, as in the previous electron scattering studies, we assume a fixed geometry for themolecule, namely the equilibrium geometry, and neglect vibrational effects altogether.
This leads to resonances that are sharper than in experimental data, where the effectof vibrational averaging is important. The cross sections are summed over the different
possible molecular orientations, as shown in Ref. [41]. The purpose of these calculations
is to show their feasibility and predictive power within our approach, detailed analysisof the spectra and their resonances have been discussed previously in the literature andwill not be examined in depth here.

4.4.1

CO2

We have tested our photoionization approach first for CO2 , and our results are
in Fig. 4.10. The ground electronic state of CO2 is 1 Σg and therefore when an electron
is extracted, in the single electron approximation, from a molecular orbital, the symmetry of the residual molecular ion will be the same as that of the orbital the ionized

79

70

2

Elastic cross section e−XeF6 ((Å) )

65

60

55

50

45

40

35

30
0

2

4

6

8

10

12

14

Energy (eV)Figure 4.8: The integrated elastic cross section for e-XeF6 shown as a function of energy.

80

Wavefunction for XeF6
15 ψ
2.5
2
1.5
1
0.5
0
-0.5
-1
-1.5
-2
-2.5

10

z

5

0

-5

-10

-15
-15

-10

-5

0

5

10

15

x

Figure 4.9: Resonance wavefunction map for e-XeF6 , for the resonance at 2.5eV. Plottedis the largest eigenvalue of the time-delay matrix. At the center it is visible the Xe atom,
whereas the four fluorines in the plane are located at 90 degrees with respect to eachother, between the yellow and blue areas, since the wavefunction has a node at thenuclei.

81
electron occupied in the target molecule. The curves in Fig. 4.10 are labeled accordingly. Comparing with the results of Ref. [119] we see quite good agreement, although
the resonance position is shifted about 0.8eV lower, due to the different description of
exchange in the two models, and also the fact that our calculations include polarization.
It is possible to see that the σu resonance at a kinetic energy of 21.8 eV couples to thegerade initial orbitals but not to the πu orbital, since the dipole matrix element is zeroin that case by symmetry. A detailed analysis of CO2 photoionization can be found in
Ref. [119].

4.4.2

SF6

In the case of SF6 we have also done both a scattering and photoionization study.
The scattering cross section for e+SF6 is shown in Fig. 4.11 and it compares favorablywith a more sophisticated theoretical model [59] and experimental data [91], both forthe overall magnitude and for the resonance position, which are as usual somewhat
shifted to higher energies with respect to experiment. The photoionization study is, atthe moment, preliminary, but it nevertheless reproduces correctly many of the features
present in experiment and other theoretical studies [207]. Our results are shown in Fig.
4.12. The main differences between our results and Ref. [207] are the sharp high energy(20 eV kinetic energy of the electron) resonances in the 5t1u and 4t1u channels, which
probably will be smeared out to a larger final width by vibrations. These resonances areprominently featured in the calculations of Ref. [112] as well, where a model potential
similar to ours is used, suggesting that this is a model-dependent feature. SF6 has 70
electrons and 35 filled molecular orbitals, therefore here we show data for the 7 orbitalgroups highest in energy. The ionization energies of these orbitals are reported in Ref.

82

14

2
Π
2 u
Σ
2 g
Πg

Photoionization cross section (Mb)

12
10
8
6
4
2
0
0

5

10

15
20
25
Electron kinetic energy (eV)30

35

Figure 4.10: The CO2 photoionization cross sections are shown versus the kinetic energyof the photoelectron, for the three highest valence orbitals (πu , πg , σg ), which generate
the doublet final states in the figure. Agreement with the results of Lucchese et al. [119]is very good. In our calculation the σu resonance occurs roughly 1 eV higher than in
Ref. [119]. From Ref. [191].

83[207]. Probably the most evident feature is the low-energy t2g resonance in the ungerade
channels, which reproduces other literature results [207, 112]. A more detailed analysisof this spectrum can be found in Ref. [207].

4.5

HCO+ Rmatrix-UK calculations
As a first example of results for the electronic degrees of freedom that we canextract from the R-matrix UK package to use subsequently in DR dynamics calculations,
we have performed a scattering calculation with HCO+ as a target. The results that weshow in Fig. 4.13 represent the quantum defect parameters at positive energy, obtainednear threshold (0.01 eV), as a function of the CH stretching and bending coordinates
respectively. The calculations were performed using a basis set from Ref. [185] and 4
target states: two singlets and two triplets, the lowest in energy.
Around a bending angle of 30 degrees a state crosses the threshold and interfereswith the π states, this is a Rydberg state relative to an excited electronic threshold.
Notice also the Renner-Teller splitting of the two π states, which at this energy is quite
small, as opposed to the large splittings present in the n = 2 Rydberg states shown in
Fig. 3.2. The π states are not exactly degenerate at linear geometry because our orbital
basis for the CI does not include all the orbitals obtained in a HF calculation of thetarget. If the CI included all the orbitals, the degeneracy would be exact, as we havetested, but the configuration space would become very large, and many of those highenergy orbitals are not needed to gain an accurate description of the scattering process
at the low energies we are interested in here. Our quantum defect results compare
favorably with those taken from bound state calculations of Ref. [49] using a Greenfunction approach and the ones we calculated with MRCI methods, shown in Sec. 6.2.1

8440
present results
Gianturco et al.
experimental

38

2

Elastic cross section e−SF6 ((Å) )

36
34
32
30
28
26
24
22
20
0

5

10

15

20

25

Energy (eV)Figure 4.11: The e-SF6 elastic scattering cross section shown as a function of energy, iscompared to the calculations of Ref. [59] and experiment by Kennerly et al [91].

85

Photoionization cross section for SF6
1000

t1g
5t1u
t2u
eg
t2g
4t1u
5a1g

Cross section (Mb)

100

10

1

0.1
0

5

10
15
Electron kinetic energy

20

25

Figure 4.12: SF6 photoionization cross sections versus photoelectron energy for theseven highest groups of valence orbitals indicated in the figure, which generate doublet
final states. Comparison with the results of Ågren et al. [207] is good. In our calculation
the high energy resonances in the t1u and t2u channels are much sharper than in Ref.
[207], presumably due to the differences in the model potentials used.

86
and Ref. [104].

4.5.1

Perturbation results

As a test of the perturbative approach developed in Sec. 2.5.4 we perform calculations on scattering of an electron from HCO+ , at negative energies, to map out thequantum defect parameters. When the appropriate quantization condition is applieddet[K + tan(πβ)] = 0,

(4.5)

with β as in Eq. 3.10, the geometry-dependent quantum defects will determine thepositions of the Rydberg states. HCO+ has an electric dipole moment and also a strong
quadrupole moment, therefore it is necessary to consider these contributions outside theR-matrix radius to obtain correct observables.
The R-matrix calculations we use to test this approach have been performed
with R0 =10a0 and only 1 target state, for the rest they are similar to what describedin the previous section, the geometry of the molecule is taken to be the ground state
equilibrium linear configuration. The highest partial wave included is a g-wave. In Fig.
4.14 a comparison is shown between different methods of dealing with the multipoles:
by perturbation theory as described in Sec. 2.5.4 (blue dots) with integrals going up to
Rf =100a0 , by R-matrix propagation [15, 127] (green dots) up to Rf =65a0 and neglecting
them altogether (red dots). The agreement between the simple perturbative approachand the propagation is very good down to -0.015 a.u. Below this energy, the effect ofstrongly closed channels is visible in the abrupt oscillations of the quantum defects,
which are smooth throughout the energy range shown in the figure when the multipoles
are omitted. The energy at which the effect of the closed channels starts to be strong is

87
1
0.9

Quantum defect (modulo 1)

0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
1

1.2

1.4

0

10

20

1.6

1.8

2
2.2
rCH (a.u.)

2.4

2.6

2.8

3

70

80

90

1
0.9

Quantum defect (modulo 1)

0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
30

40

50

60

Bending angle (degrees)

Figure 4.13: Top: R-matrix calculation of eigenquantum defects at threshold (0.01 eV)
for HCO+ plus an electron, as a function of the CH stretching coordinate, with the CO
coordinate and bending angle kept at their ionic ground state equilibrium values. The
squares represent the A0 states, while the asterisks represent the A00 states. Bottom:
R-matrix calculation of eigenquantum defects at threshold (0.01 eV) for HCO + plus
an electron, as a function of the bending angle, zero degrees represents the linear ionic
equilibrium geometry.

88
a function of Rf , in general the higher this value the higher the energy at which theseeffects are visible. We would like to point out that using the perturbative approachdeveloped in Ref. [12] the results obtained are very similar, except that strong oscillatory
behaviors result at very negative energies.
In general it is not possible, in a system as complex as the one we have considered,to obtain quantum defect parameters that are smooth in an energy range that extends
from threshold down to very negative energies with a single method. Both perturbationand propagation will at some point be affected by the strongly closed channels. Even if
the closed channels are treated as such by channel elimination, Eq. 3.10, they still affectthe remaining quantum defects as sharp avoided crossings at the closed channel energy.
It is probably a better solution, especially in the case of polar molecules, to separate
the very negative energy part of the spectrum and treat it in a different fashion.

89

unperturbed
propagation up to 65a0
perturbation

0.4

Phase shift

0.2

0

-0.2

-0.4
-0.03

-0.025-0.02-0.015
Energy (au)-0.01

-0.0050

Figure 4.14: Comparison between unperturbed results (red dots) at R 0 =10 a0 , perturbation theory (blue dots) with integration up to Rf =100a0 , and propagation from
R0 up to 65 a0 (green) for quantum defects of HCO+ plus an electron at negative energy. The perturbed values are very close to the correct ones obtained with R-matrix
propagation, at low energies the effect of strongly closed channels is important. The
R-matrix calculation includes up to l = 4, which will generate closed channels below
E = −1/(2n2 ) = −0.0156 a.u. for n = 5.

Chapter 5

Results: electron scattering, DNA bases

5.1

Introduction
In recent years, increasing importance has been associated with electron-induced

chemical processes in biological environments, especially in relation to radiation damage
to nucleic acids (DNA and RNA). These processes are initiated by the interaction ofionizing radiation with living tissue, generating possibly mutagenic and carcinogenic
byproducts. A wide variety of subsequent ionization, excitation and energy transfer
processes can generate possibly mutagenic and carcinogenic byproducts, which can affectmany molecular species in the complex cell environment.
The important work of Sanche and coworkers [24, 111, 159] has shown that damage to nucleic acids from ionizing radiation [140] (single and double strand breaks inparticular) can be generated through a mechanism involving low energy electron attachment to the nucleic acid and subsequent bond breaking due to energy transfer toa vibrational mode of the temporary anion formed in the electron capture step. The
ionizing radiation produces high energy primary electrons (with energies of the orderof a keV), which in turn, through electron-impact ionization, produce many low energy
secondary electrons. In the electron-impact ionization process, the scattered electron

91
loses part of its kinetic energy, while another electron is ejected, with energy much
lower than the first one. It has been shown that the secondary electrons have an energy
distribution concentrated between 0 and 20 eV [34]. A schematic representation of theradiation damage process is shown in Fig. 5.1.
If the electron energy is higher than the ionization threshold for DNA (between
7.85 and 9.4 eV, as measured for the DNA bases [114]), then the nucleic acid can be
ionized and damage produced through the subsequent rearrangement and dissociationof the cation [54]. If the electron energy is lower, damage can still be generated througha negative ion-mediated mechanism, that starts with the capture of the electron in amolecular resonance, followed by a transfer of energy and electron density into a weakbond that subsequently ruptures. Beyond this generally accepted framework, there aremany controversial issues, such as the locus of the initial capture site, [24] the evolutionof the anion produced and which bond actually ruptures.[17, 210]
In the past few years many studies have been devoted to understanding the mechanisms triggered by the low-energy electrons, including especially their capability tocause strand breaks [5, 17, 78, 121, 164]. A first general feature on which there is wide
agreement is that the electron capture is mainly due to the DNA and RNA bases. Thesemolecules have extended aromatic systems, whereby there are numerous low-lying unoccupied π ∗ orbitals into which an electron can be captured. The capture often occurs
into a shape resonance, a temporary anion, in the range of energies between 0 and 15 eV,
where experiments have found signatures of electron-induced damage to nucleic acids.
The simplest of these bases are thymine, cytosine, uracil (pyrimidines, monocyclic) and adenine and guanine (purines, bicyclic and generally larger than pyrimidines). Their structures are shown in Fig. 5.2. In this chapter we will present theoret-

92

Figure 5.1: Schematics of the DNA damage from low energy electrons: highly energetic
ionizing radiation (in red) impacts a molecule in the cell (water in this case) extracting
keV primary electrons (in blue), which in turn extract low energy secondary electrons(in black) through electron-impact ionization. These impact DNA and are captured by
the bases (shown in the inset is a uridine nucleotide, it will be replaced by thymidine inDNA), then an energy transfer through the sugar (in green) breaks a bond and generates
a strand break. The B-DNA structure is taken from Ref. [83].

93
ical predictions of cross sections and of even greater relevance, collisional time delays,
for elastic electron scattering from these large molecules. Determination of the location, width, and electronic structure of resonances for a single target molecule is an
important step towards understanding and possibly modeling the complex dynamics ofDNA, which consists of multiple components. Specifically, besides the bases, there arealso the sugar backbone, the phosphates, and the solvation water [22], that probably
plays a major role in stabilizing the temporary anions [70, 176]. No previous theoreticalor experimental study of low-energy electron scattering from all of the DNA and RNA
bases is available for comparison (although a study at intermediate energy has been
carried out recently [140]), but our method has proven its reliability in the study ofsmall molecules [192] and more extended systems like C60 , SF6 and XeF6 , as we haveshown in Sec. 4.3.2 and 4.4.2.
Much experimental work also has been carried out on dissociative electron attachment from the DNA bases [2, 66, 121], to understand what fragments are generated.
We will discuss the possible connections with measured dissociation branching ratiosthat we can infer from the examination of the spatial shape and nodal surfaces of theresonant wavefunctions.

5.1.1

Elastic cross sections

To our knowledge there are no available experimental data or calculations of low
energy electron scattering from the complete set of DNA bases. A study of electron
attachment has been presented in Ref. [5] and the resonance positions are clearlymarked. Compared to these results, our calculations show resonances shifted typicallyby about 2 eV higher in energy, but the energy spacing of the resonances is comparable

94

Figure 5.2: Ground state equilibrium structures of the molecules considered in thissection The red atoms are oxygens, the orange circles represent carbons, the blue atoms
are nitrogens while the small circles are hydrogens. From Ref. [194].

95

Figure 5.3: Blowup of the resonant wavefunction surface for uracil at 2.4 eV (Fig. 5) nearan oxygen atom. This surface is a cut through the plane containing the nuclei, shownin the region of the oxygen atom attached to C2 . The dots are the actual grid points ofour calculation, while the coarsely spaced grid is a linear interpolation of the resonant
wavefunction (eigenstate of Eq. 5.1) and the finely spaced grid that shows a cusp is theK-shell orbital of oxygen, here represented to show how our grid is designed around thenuclei. The colorbar on the right side refers only to the resonant wavefunction, while
the K-shell orbital is colored just to make it more visible. The K-shell orbital is shifted
upward by 0.04 to make it more clearly visible, therefore its base is not a nodal surface.
From Ref. [194].

96
to what is observed in the experiment. Moreover the relative values of the widths ofsuccessive resonances resemble the measured widths. There is also a theoretical study atintermediate energies [140], and calculations of electron scattering from uracil [60, 70];in the following we compare these results to ours.
We have already mentioned that the heterocyclic DNA bases have many low-lying
unoccupied orbitals, so it is not surprising that their elastic cross sections for electron
scattering exhibit many shape resonances. These can be viewed as a capture of thescattering electron into one of these antibonding orbitals to form a short lived negative
ion state [40, 116].
Since all these molecules have, in their equilibrium configuration, only one symmetry element - reflection through the molecular plane - we will characterize the resonances
as being of σ type (no node in the plane) or π type (when they have instead a node in
the plane) rather than using the A0 and A00 labels as is customary for the Cs group.

5.1.2

Positions and widths of resonances

A general comparison of partial elastic cross sections for all five of these molecules
is shown in Fig. 5.4, while in the following we give a more detailed description andcompare with information available in the literature. Also, a plot of total time-delays(see also Sec. 5.1.3 for details) is provided in Fig. 5.5 to show the resonances in moredetail.
Since we are dealing with polar molecules, use of the fixed-nuclei approximation
as it stands makes the partial wave expansion of the forward scattering amplitude divergent. Due to the long-range nature of the dipole interaction, in fact, all partial waves
would contribute to the scattering process, causing an infinite differential cross section

97
in the forward direction and therefore infinite integral cross sections. There is a method,
extensively discussed in the literature [32, 61], to deal with this problem by means ofa Born closure formula, which yields a finite integral cross section once molecular rotations are included. We will not pursue this further, since existing experiments are
not likely to deal with such detailed rotational structures. Therefore our cross sections
and time-delays include only up to lmax = 10 and omit all higher partial waves. The
correction would be proportional to the dipole moment and inversely proportional to
the smallest rotational spacing. For the DNA bases the dipole moment is large, while
the rotational spacing is small. Therefore the correction can be quite large especially atvery low energy. The correction would thus tend to mask the resonant structures, which
are the most interesting observables and which have been measured in experiments. Allof the calculated cross sections grow rapidly when the incident electron energy decreases
below 1 eV, which is a signature of the role played by the dipole field in pulling in theelectron and which is very common in electron scattering from polar molecules [5, 69].
A comparison of our resonance patterns with the electron transmission spectroscopy (ETS) data of Burrow and coworkers [5, 164] can be found in Figs. 5.6-5.8,
for a comparison of energies, widths and spacing. All of the resonances obtained in our
calculations are listed in Table 5.1.

5.1.2.1

UracilIn the cross section of uracil we find 3 resonances, at 1.3 eV (of width 0.2 eV),
at 3.5 (0.6 eV wide) and a very broad resonance at 6.3 eV. The resonance at 1.3 eV
is dominated by the l = 3 partial wave (50%) and has contributions from l = 1 (35%)
and l = 2 (11%). at 3.5 eV the main partial wave is d (66%), at 6.3 eV f-wave is the

98

Molecule
UracilThymineCytosine

Adenine

GuanineEnergy (eV)
1.3
3.5
6.3
1.5
3.8
6.8
1.6
4.0
8.1
2.4
3.2
4.4
9.0
2.4
3.8
4.8
8.9
12.0

Width (eV)
0.2
0.6
0.9
0.3
0.6
1.0
0.3
0.7
0.8
0.2
0.2
0.3
0.5
0.2
0.25
0.35
0.6
1.0

Partial wave
3 (50%)
2 (66%)
3 (64%)
3 (53%)
2 (62%)
3 (61%)
3 (51%)
2 (68%)
3 (63%)
2 (65%)
2 (62%)
3 (51%)
5 (53%)
2 (46%)
2 (44%)
4 (38%)
5 (33%)
4,5 (both 23%)

Table 5.1: Energies, widths and largest partial waves (and percentages) of the resonances
discussed in the text for the DNA and RNA bases. From Ref. [194].

99

Elastic cross section e−DNA bases (Å2)

200

UracilThymineCytosine
Adenine
Guanine

180
160
140
120
100
80
60
0

2

4

6

8

10

12

14

Energy (eV)Figure 5.4: Partial elastic cross section for the 5 DNA and RNA bases described in the
text. Calculations involve partial waves up to l = 10 and the dipole physics outside theR-matrix box is included exactly, for those partial waves. From Ref. [194].

Total time delay

100

30

Guanine030

Adenine

030

Cytosine

030

Thymine030

Uracil

0
0

2

4

6
8
Energy (eV)1012Figure 5.5: Total time-delay for the DNA bases described in the text. From Ref. [194].

101
dominant contribution (64%).
In the work of Gianturco et al. (see Ref. [60]) three π ∗ resonances are found atenergies of 2.2, 3.5 and 6.5 eV. The second and third π ∗ resonances from that work
fall at lower energies than ours, a somewhat surprising discrepancy since the theoreticalmodels are very similar.
The contribution of the dipole field at distances larger than 12 Bohr is neglected
in Ref. [60], but we have noticed that this influences only the overall magnitude of thecross sections (roughly an increase of 20% at very low energy, that is reduced to about
5% around 10eV), the dipole physics only weakly affects the resonance positions andwidths.
Resonances are measured in Ref. [5] to occur at 0.3, 1.5 and 3.8 eV. They are allassigned as π resonances [164], so our results should be shifted by about 1.1 eV down,
whereas the spacing between the resonances is larger than experiment. The relative
resonance widths are similar to Ref. [164], in that the first resonance is very narrow,the second broader and the third very broad, a comparison is shown in Fig. 5.6, where
an integration of the experimental data has been performed to show more clearly theresonance positions and widths.

5.1.2.2

Cytosine

For cytosine we find 3 main resonances, a very sharp one at 1.7 eV (width 0.3 eV),then at 4.0 eV (width 0.7 eV) and a third at 8.1 eV (width 0.8 eV). The dominant angular momentum character of the resonances is the same as for the three corresponding
resonances of uracil. Comparing the resonance positions with the experimental data ofRef. [5] we see the same general trend already observed with uracil, namely an overall

102
shift higher than experiment of all resonances by about 1.4 eV. Interestingly, the firsttwo resonances are measured to occur at an energy lower than in uracil, a trend that is
not verified in our calculations. However the third resonance is at a much larger energy
than in uracil and thymine, which also happens in the experiments.

5.1.2.3

Thymine

The scattering cross section for thymine is closely similar to uracil, which is not
surprising in view of their close structural similarities, this applies to both the magnitude and the position of the resonances, which are slightly shifted to higher energies.
Specifically, we find resonances at 1.5 eV (width 0.3 eV) at 3.8 eV (width 0.6 eV) andat 6.8 eV (width 1 eV).

5.1.2.4

AdenineThe electron scattering spectrum for adenine presents many resonances, as expected due to the complexity of the target structure. Also very interesting is the factthat the cross section drops sharply at energies below 2 eV, a behavior opposite to thatfound for the other molecules, if we do not consider the dipole physics outside the Rmatrix box, whereas a zero-energy peak appears in the full calculation, a possible signof a dipole bound state right below threshold.
The first resonance occurs at 2.4 eV (width 0.2 eV), the second at 3.2 eV (sharp,
width 0.2 eV), then another centered at 4.4 eV (0.3 eV wide), while at 9 eV we have abroader resonance of width 0.5 eV. The dominant partial wave of the first two resonancesis l = 2 (65% and 62% respectively). The third resonance is l = 3 at 51% and l = 4
at 33% The resonance at 9 eV is dominantly l = 5 (53%) with an l = 3 contribution

103(22%).
Compared to experiment we have a shift of all resonances roughly 1.5 eV higher,
as in guanine, in this case the spacings are correct (about 1 eV between the first threeresonances, while the fourth falls too high in energy and it is not measured in experiment). Also the experimental widths of the first three resonances are very similar, as
in our data. A comparison with the data of Ref. [5] is shown in Fig. 5.7.

5.1.2.5

GuanineFor guanine we find 4 resonances: at 2.4 eV (width 0.2 eV), at 3.8 eV (width
0.25eV), a third at 4.8 eV (width 0.35 eV), then at 8.9 eV (width 0.6 eV) and a broadresonance around 12 eV. Each of the first three resonances has strong contributions
from d, f, and g-waves. At 2.4 eV the contributions are 46% for l = 2 and 37% for l = 3,
for the second resonance l = 2 is 44% while l = 4 is 32%, the third is 38% of l = 4
character and 35% of l = 3. The resonance at 8.9 eV is 33% h-wave, 28% f-wave and20% g-wave. At 12 eV the composition is: l = 4 and l = 5 equally at 23%, while l = 6
contributes a further 13%.
Comparison to experimental data (see Fig. 5.8) shows again a shift of 1.5 eV
overall, while the resonance spacing is well reproduced, and the second and third resonances fall at higher energies with respect to adenine, as in our calculations. Also thewidths seem close to experiment.

5.1.3

Resonance molecular structures

From the shapes and nodal structures of the resonant states it is possible toattempt a discussion of the dissociation patterns observed experimentally, if we consider

104

5

integrated data from Burrow and co-workers
scaled time delay, shifted down by 1.1 eV
fπ

4

dπ

arbitrary units

fπ
3

2

1

0
0

1

2

3

4

5

6

7

8

9

E(eV)

Figure 5.6: Comparison of our calculated time-delay with the experimental scatteringcross section data of Burrow and co-workers [164] for uracil. The arrows indicate resonance positions from the present work, while labels show the dominant partial wave ofresonance. The time-delay curve is shifted downward by 1.1 eV to have the position of
the first resonance coincide with experimental data. From Ref. [194].

105

14
dπ integrated data from Burrow and co-workers
scaled time delay, shifted down by 1.5 eV
dπ

12

fπ

arbitrary units

10

8

6

4

2

0
0

0.5

1

1.5

2

2.5

3

3.5

4

E(eV)

Figure 5.7: Comparison with experimental data of Burrow and co-workers [5] for adenine. The arrows indicate resonance positions from the present work, while labels show
the dominant partial wave of resonance. The time-delay curve is shifted downward by
1.5 eV to have the position of the first resonance coincide with experimental data. From
Ref. [194].

106

14
integrated data from Burrow and co-workers
scaled time delay, shifted down by 1.5 eV
dπ

12

dπ
gπ

arbitrary units

10

8

6

4

2

0
0

0.5

1

1.5

2

2.5

3

3.5

4

E(eV)

Figure 5.8: Comparison with experimental data of Burrow and co-workers [5] for guanine. The arrows indicate resonance positions from the present work, while labels show
the dominant partial wave of resonance. The time-delay curve is shifted downward by
1.5 eV to have the position of the first resonance coincide with experimental data. From
Ref. [194].

107
the resonant states as being precursors for dissociative states. Caution must be used,though, in drawing conclusions from this analysis, because this involves a certain degreeof speculation. In fact, to establish once and for all the dissociation patterns of thesecomplicated systems, scattering calculations at many different geometries would haveto be carried out, and the nuclear dynamics should be included, At present this iscomputationally too expensive to contemplate. The first two resonances observed for
uracil are shown in Fig. 5.9. The quantity plotted is a projection on the molecular planeof the eigenvector corresponding to the maximum eigenvalue of the time-delay matrix[74]
Q = iS

dS†
,
dE

(5.1)

where S is the scattering matrix. At the energy where the time delay of the resonanceis a maximum, this eigenvector constitutes the dominant contribution to the resonant
structure, since it corresponds to the partial wave that experiences the maximum timedelay in the scattering process. For sufficiently narrow resonances one eigenvalue isalways dominant, making the resonance analysis much easier. The eigenvectors of thetime-delay matrix are complex, so we adopt a phase factor such that the highest peak of
the wavefunction is a purely real number. We then find that the resonance wavefunction
is real everywhere, to a good approximation (the imaginary part is about 10 −6 smaller
than the value of the real part), and we plot only the real part. We analyze in detail
only the cases of uracil and adenine, the other pyrimidines being very similar to theformer and guanine to the latter.
The nodal patterns for uracil are very similar to the ones showed in Ref. [60],which is not surprising since the approximations made in that work are similar to ours,

108
as already discussed, therefore we show only the first two resonances. Incidentally, we
notice the close resemblance of these resonant wavefunctions to the first virtual orbitalsof uracil from a HF calculation performed with a small basis set (6-31G*), see forexample Fig. 5.10. In Ref. [60] also a very low σ ∗ resonance is found at 0.012 eV. Ourcross section grows substantially at low energy. If we plot the eigenstate corresponding
to the largest eigenvalue of the time-delay matrix (as described in Sec. 5.1.3), as in Fig.
5.9, corresponding to this low energy range, it looks similar to Fig. 5 in Ref. [60], withthe main differences being that in our wavefunction the N3 -H bond has a node, there
is a large excess charge on N3 and on the oxygen attached to C4 , while another nodal
surface cuts diagonally from C2 to C5 . This state anyway does not appear to be so
relevant in the experimental data [164], where mainly the π ∗ resonances are detected.
We can also see that there is accumulation of electronic density (the peaks of
the wavefunction) on the ring structure, and that many of the ring bonds have nodal
surfaces cutting through them, so capture in these resonant states can be reasonably
thought as leading to a fragmentation of the molecule in which the aromatic ring is broken. Experimental dissociation patterns are illustrated for Br-uracil in Ref. [2], where
evidence for breaking of the ring structure lies in the peaks at 1.6 and 3.5 eV producedby (OCN)− and other fragments. These fragments can be generated by capture into
shape resonances, appearing in our calculations at 1.3, 3.5 and 6.3 eV respectively. In
particular there is a nodal surface in the 3.5 eV resonance that encloses the C 4 -N3 bond,which could generate a CN− fragment.
Since our calculations do not take into account core-excited states or vibrations,
our results do not include any Feshbach resonant structures. These appear to cause atleast some of the patterns observed in experiment, as in the case of uracil, [164] and they

109
will presumably constitute the dominant trapping pathways for energies higher than 7
eV, where the number of electronic Feshbach resonances starts to become very large, as
shown in Fig. 5.13.
In the case of uracil, we looked carefully for a σ ∗ resonance that might be similar
to the state shown in Fig. 3 of Ref. [164] around 3 eV at equilibrium geometry,a dissociative state most likely responsible for N1 -H bond cleavage. The Ref. [164]calculation was performed by scaling Hartree-Fock continuum orbital energies, so no
information about the width was provided. Such a state, taking into account an expected
shift of 1-2 eV upward in our calculations, should have appeared at around 5-6 eV, and it
was not found. This is probably due to the fact that this resonance is extremely broad,since it is also not seen even in experiment [164]. Moreover in calculations carriedout using complex absorbing potentials, in connection with Green’s function methods[160, 51], for similar systems (like benzene [50]), analogous σ resonances were extremely
hard to detect. They became narrower (around 1eV width at equilibrium) only when
the relevant hydrogen was substituted with a heavy atom like chlorine; this was also
demonstrated experimentally in the case of Cl-uracil in Ref. [164].
For adenine there is less experimental information available to compare. In Ref.
[66] it is stated that the dominant breakup channel for low energy (0-4eV) electron
attachment leads to hydrogen atom loss, and very prominent resonant structures arepresent in the range 1 to 3 eV. Looking at the resonance wavefunction maps in Fig.
5.11-5.12 we can see that there is no significant buildup of electronic density on any of
the hydrogens, consistently with the fact that the negative charge stays on the molecular
frame, and therefore there is no H− formation.
The first few unoccupied molecular orbitals that can be obtained from a Hartree-

6

110
0.08
0.06

4

H
0.04

y (Bohr)

2

N3

O

0.02

O

0

C4

C2

−0.02

0

−0.04

N1

−2

C

−0.06

5

H

C6

−0.08

H

−0.1

−4

−0.12

H
−6
−6

−4

−2

0

2

4

6

x (Bohr)
5
0.1
4

H

3

y (Bohr)

2

N3

O

1

O

0.05

C

C

4

2

0
0
−1

N1

−2

C5

H

−3

H

C6

−0.05

−4
−5

H
−5

−4

−3

−2

−1

0

1

2

3

4

5

−0.1

x (Bohr)Figure 5.9: Top: Eigenvector associated with the dominant eigenvalue of the time-delay
matrix in uracil at a low but nonresonant energy, as described in the text, for scattering
at 0.2 eV. Bottom: first resonance wave function in uracil. Eigenvector associated withthe dominant eigenvalue of the time-delay matrix, as described in the text, for the
resonance at 1.3 eV. The green squares correspond to the nuclei positions. From Ref.
[194].

111

Figure 5.10: Uracil: structure of the virtual orbital associated to the resonance at 1.3
eV. The energy of this virtual orbital is 3.42 eV when using a 6-31G** basis set. The
orientation of the molecule is the same as in the previous plots. The black and white
lobes correspond to opposite signs. It is possible to notice the node in the molecularplane that makes this a π ∗ orbital. From Ref. [194].

112
Fock calculation for adenine, as in the case of uracil, are extremely similar in their nodal
structures to our resonant wavefunctions, so these shape resonances can be viewed quite
reasonably as the trapping of the scattered electron in a virtual orbital.
Most of the C-C and C-N bonds have nodal surfaces passing through them, Thismight suggest that other channels that involve the breakup of C-C and C-N bonds could
also be available at these energies, although probably they are less important than thehydrogen loss products.

5.1.4

Uracil - Rmatrix UK

We have studied uracil in particular with the more accurate R-matrix UK code,
described in Sec. 2.4.2. The aim of this study was to determine potential energy curves
as a function of the N1 H bond stretching coordinate, since in a study by Burrow et al.
[164] it has been suggested that this coordinate can lead to dissociation of the hydrogenin a simple one-dimensional picture, without the need to consider couplings with outof-plane motion or other normal modes.
We first have performed a simple calculation for the equilibrium geometry, shownin Fig. 5.13 using the 4 lowest states of the target (two singlets and two triplets) inthe close coupling expansion and a 6-31G** basis set. The CI used here is forcibly very
limited, with only 8 electrons allowed to move, therefore we know that the dissociation
limit will not be well described. The polarization that we allow for the target is also
somewhat limited, therefore these results will be relatively similar to static exchange
calculations. The results we obtain for the resonance positions are similar to our previous
static exchange plus polarization calculations, see also the discussion in Sec. 5.3. In
this case we can separate the different symmetry contributions in A0 (symmetric with

113
0.025

6

0.02

H

4

H

N

0.015

7

C

8

0.01

N

y (Bohr)

2

H

C

5

C

N

6

9

0

0.005
0

C

H

4

−0.005

N

−2

1

N

−0.01

C2

3

−4

−0.015
−0.02

H

−6
−6

−4

−2

0

−0.0252

4

6

x (Bohr)
0.06
6
0.05

H

4

H

N7

C

0.04

8

y (Bohr)

0.03

N

2

H

C

5

6

9

0

0.02

C

N

0.01

C4

H

N1

−2

N3

0

−0.01

C2

−4

−0.02

H

−6

−0.03
−6

−4

−2

0

2

4

6

x (Bohr)Figure 5.11: Top: 1st resonance wavefunction for adenine. Eigenvector associated withthe dominant eigenvalue of the time-delay matrix, as described in the text, at 2.4 eV.
Bottom: 2nd resonance wavefunction for adenine. Eigenvector associated with thedominant eigenvalue of the time-delay matrix, for the resonance at 3.2 eV. From Ref.
[194].

114

0.6

6

H

4

7

C8

N

2

y (Bohr)

0.4

H

N

C5

C

N

6

9

0

0.2

H
0

C

H

4

N

−2

−0.2

1

N

C

3

−4

2

−6

H

−0.4

−0.6
−6

−4

−2

0

2

4

6

x (Bohr)Figure 5.12: Adenine: 3nd resonance wavefunction. Eigenvector associated with thedominant eigenvalue of the time-delay matrix, as described in the text, for the resonanceat 4.4 eV. From Ref. [192].

115
respect to the molecular plane) and A00 (antisymmetric), therefore it is possible to show
easily how the resonant structures are all in the A00 (or π) channel. It is also possible toevince from the data that the first three π ∗ resonances are shape resonances, while theothers are Feshbach type resonances. According to the results of a simple Hartree-Fockcalculation, in Ref. [164] it was stated that a dissociative σ ∗ resonant curve should
be present at low energy. Performing scattering calculations we found no indicationof this resonance, which means that its width is extremely broad (>1eV), as statedpreviously. The results as a function of the N1 H stretch, shown in Fig. 5.14, are similar
to the curves in Ref. [164], in that the π ∗ resonances are parallel to the neutral groundstate. The position of the σ ∗ resonance, though, cannot be determined. A resonant
state becomes bound roughly at the same position predicted by the HF calculations,and this is probably the σ ∗ state predicted in Ref. [164]. The behavior of the boundpart of this curve is incorrect, and a larger configuration space (we used roughly 10000
configurations in the scattering calculation) will be needed to get a quantitative result
for these potential curves. Qualitatively though, our results confirm the predictions ofRef. [164], with the exception that a direct capture in the resonant state is not possible,
so the dissociation has to happen either through a coupling with the dipole boundstate (of the same symmetry) through vibrational Feshbach resonances, or through an
out-of-plane motion that couples it to the π ∗ resonances.

5.2

DNA backbone
In this section (see also Ref. [195]), we analyze the behavior of the other subunits

in the macromolecule, with the aim of investigating the electron interaction with themoieties that constitute the sugar-phosphate backbone, and verify whether they could

116

Eigenphase sums for electron - Uracil
4

A’ symmetry
A’’ symmetry

Eigenphase sum (radiants)

3
2
1
0
-1
-2
-3
-4
-5
0

2

4

6
8
Energy (eV)

10

12

14

Figure 5.13: Eigenphase sums for the two symmetries (A0 and A00 ) of the uracil moleculeat the ground state equilibrium geometry as a function of energy. It is possible to noticethat the A0 channel has no low energy resonances, while the A00 has 3 shape resonancesand many Feshbach resonances above 8eV. Note that the downward discontinuities are
simply branch changes, rather than resonances.

117
4
3.5
3

Energy (eV)

2.5
2
1.5
1
0.5
neutral ground state
dipole bound state
π*
σ*

0
−0.5
0.6

0.8

1

1.2

1.4

1.6

1.8

2

2.2

2.4

2.6

N1H distance (Å)2

12

Neutral ground stateR−matrix pole and bound state
First A’’ resonance
Second A’’ resonance
A’ resonance

10Energy (eV)

8

6

4

2

0
0.8

1

1.2

1.4
1.6
1.8
N1−H distance (Å)

2

2.2

2.4

Figure 5.14: Top: potential energy curves for neutral uracil and the first few anionicstates calculated at the Hartree-Fock level, based on calculations from from Ref. [164].
Bottom: energies of the ground state of the uracil molecule and its first few resonant
anionic states. The π resonances are parallel to the ground state, while the σ statecrosses it. Although the unbound part of the σ state is too wide to be detected as aresonance, at the N1 H distance where this state is just above the ionization threshold,a small resonance is seen, which is represented as a square in the figure. The longdashed curve represents the first R-matrix pole when it is above the ionization threshold,whereas after it becomes bound it is the energy of the bound state.

118
play a role in the electron capture stage of the radiation damage process. In any case, it
has been shown [153] that low energy electrons can damage the DNA sugar, deoxyribose(2-deoxy-D-ribose), and this interaction could be important in the cell environment. To
model the sugar, we use tetrahydrofuran (THF, whose formula is C4 H8 O), which issimilar to deoxyribose, except that the latter has -OH groups attached to the C 1 andC3 and a CH2 OH side chain attached to the C4 (the latter two substituents are linked
to a phosphate group in the DNA molecule, while the former is linked to a base), they
are substituted here by hydrogens. We use THF instead of deoxyribose because in DNA
the OH groups are fundamentally modified by the phosphate and base that are attached
to it (in fact the one attached to the C1 disappears altogether), so we think that thesolution would be either to consider a whole nucleotide or use THF, which mimics
only the ring structure and is less affected by the rest. Another reason for this choiceis that a recent study on dissociative electron attachment (DEA) from all the DNA
components [6] has measured the average gas-phase DEA cross section from the bases
as being very similar to the cross section per base for supercoiled DNA measured in Ref.
[25]. This implies that the DEA cross sections for the sugar and phosphate groups aremuch smaller than those of the bases, which is what happens in the case of THF, but
not for a hydroxy-substituted THF, where the DEA cross section is larger than for thebases. This suggests that THF could be a better molecule to model the sugar moiety inDNA. The structures of these compounds are represented schematically in Fig. 5.15. Inpractice, the DNA backbone can be thought of approximately as constituted by THF
molecules to which the bases are attached, linked together by phosphate groups.
For the phosphate there are no previous results available for comparison, but
for THF recently there have been both new experimental [131, 209] and theoretical [23]

119
results, with which we will compare our calculations. We will also attempt a comparisonof the resonant wavefunctions to virtual orbitals from a Hartree-Fock calculation, inorder to gain some more insights into the capture process and its possible consequenceson the anion evolution. The conclusions suggested by this analysis lead us to draw a
possible link between our calculations and experimental data on dissociative electron
attachment.

5.3

Results: THF
In Fig. 5.16 our results are compared to the low-resolution experimental dataof Zecca et al. [209] and the theoretical results of Bouchiha et al., [23] obtained likeours without performing a Born closure to consider the effect of the dipole field on thehigher partial waves. We also plot the data in the form of a time-delay, [194] to makethe resonance position and width more evident. The total time-delay is the sum of theeigenvalues of the hermitian matrix Q of Eq. 5.1. In Fig. 5.16 we also plot the resonant
channels, particular eigenvalues of the time-delay matrix that are larger than the others
and show a Lorentzian behavior that sets them apart from the rest. We rescaled thetotal time-delay by dividing it by a factor, to show more easily everything in the samegraph.
The order of magnitude of our cross section is not different from the results of Ref.
[23], also calculated without Born closure. Since our calculations do not include excitedstates of the target, we will not be able to detect core-excited resonances (these were
included in the calculations of Ref. [23]). The experimental total cross section contains
also rotational and vibrational excitations and electronically inelastic channels, therefore
it is to be expected to be larger than the elastic cross section. Our calculated elastic cross

120

H

OH

H
C

O

OH
O

OH
(I)

(II)

(III)

Figure 5.15: Top: Three dimensional structures of phosphoric acid and tetrahydrofuran.
The small circles are hydrogen atoms. Bottom: Schematic structures of cyclopentane(1), THF (2), and the DNA sugar deoxyribose (3), that show the similarities between
these compounds. The hydrogen atoms that fill the carbon valences are not shown.
From Ref. [195].

121
sections lie instead above the experimental total cross sections of Ref. [209], and this ispartly due to the approximate nature of our model. The comparison with experiment isdifficult due to the effect of higher partial waves (with l > 10). The Born correction has
been calculated in Ref. [23], and the cross section becomes higher than experiment, so
we would expect a similar effect for our data. Due to the low resolution, the experimental
results below 5 eV are (in the words of the authors) to be considered only indicative.
The resonant structure that appears around 8.6 eV in our results looks similar to theanalogous experimental feature, also taking into account that our model usually predicts
resonances about 1-2 eV higher than their experimental position, as verified for the DNA
bases [194] and also for molecules like carbon dioxide, benzene and SF6 [193], becauseof the approximations adopted. This is due mostly to the adoption of a local exchange
model, since for example, in a system like CO2 [63, 102] this approximation causes a
resonance shift of around 2 eV or so for model static-exchange, compared to all-electron
static-exchange calculations, that is not possible to compensate with the polarizationcorrelation potential, which indeed shifts the resonance by a similar amount for thismolecule (2.7 eV compared to 2 eV) with respect to correct polarization. Also themodel polarization is approximate, but the model exchange constitutes the largest error
in our calculations. Indeed, if we perform a purely static-exchange calculation (see Fig.
5.16) the position of the resonance (maximum of the time-delay) is at 11.2eV, 2.6 eVhigher than the result including polarization, and it is composed of three overlapping
resonances. For what we said above about the shifts generated by the model potential
with respect to exact exchange, this resonance should be around 9 eV in an all-electron
calculation; hence it should probably have been visible in the calculations of Ref. [23]that go up to 10 eV in energy. When polarization-correlation is added, this resonance122
energy should fall below 10eV, and thus be clearly visible in the data of Ref. [23].
From the fact that our resonances are shifted too high in energy one might havethought that the widths would turn out to be larger than the experimental ones. We
have verified in the past for many systems [193, 194] that our calculated widths arecomparable to their correct values.
In fact no resonance was found in this energy range (0-10 eV) in the other published theoretical study [23]. Through a Lorentzian fit, we predict a width of about 2.3
eV for this resonance. At the static-exchange level, the resonance has a larger width,of about 3.0 eV. A higher energy resonance is found at 14.1 eV and it is even broader(about 2.5 eV). The partial wave contributions are mainly l = 4 (with a contribution of70%) for the lower energy resonance, while for the higher energy resonance the largest
are l = 5 (50%) and l = 3 (30%, as can be seen in Tab. 5.1). The partial wave decomposition will depend somewhat on the choice of the expansion center. We always
perform the decomposition around the center of mass of the molecule, which seems areasonable choice. The exception is the calculation of HeH+ quantum defects in Sec.
4.3.1, where the center of charge was instead used.
To better understand the reason for the difference between our results and theones in Ref. [23] and to make sure that our results are consistent with other known
results, we have first of all performed a time-delay analysis [74, 174] of our scattering
data, shown on the bottom of Fig. 5.16. This has allowed us to establish the presenceof three superimposed resonances, two of which peak at 8.1 eV (with one significantly
more intense than the other) while the third peaks around 9eV, and its intensity falls
between the other two. Due to the large width of these resonances, the weakest of thechannels in Fig. 5.16 has a time-delay eigenvalue which is not much larger than the

123

Elastic cross section (Å)2

70

Static exchange + polarization
Static exchange
Zecca et al. (exper.)
Bouchiha et al. (theor.)

60
50
40
30
20
10
0
0

5

10

20

Energy (eV)total time delay
1
2
3
4

30
25
Time delay

15

20
15
10
5
2

4

6

8

10
12
14
Energy (eV)

16

18

20

Figure 5.16: Top: partial elastic electron scattering cross sections from THF (solid line),
the dot-dashed line represents the static-exchange results. Calculations involve partial
waves up to l=10, for which the dipole physics outside the R-matrix box is included
exactly. Comparison with theoretical (dotted) [23] and experimental (dashed) [209]results, the open squares are experimental data points. Bottom: time-delay plot tohighlight the presence of resonances. The full curve is the total time delay, rescaled
by a factor of 6 to show all curves on the same graph more easily, while the numbered
curves correspond to the few highest eigenvalues that exhibit resonant behavior. From
Ref. [195].

124
nonresonant eigenvalues.
An existing calculation for cyclopropane [37] (C3 H6 ) shows that this cycloalkane
has a cross section with a similar resonant structure at 6 eV. It is superimposed on abroader resonance at higher energy and other nonresonant contributions that give a very
wide plateau, although in this case the first feature is not a composite resonance. The
symmetry of the cyclopropane resonance in Ref. [37] is even with respect to the planeof the three carbon atoms, and the same occurs in THF, in the sense that although thiscompound is not planar, our resonance structures do not change sign above and below
the ring bonds. It is plausible then that due to its larger frame and density of virtual
orbitals, THF is able of supporting more quasi-degenerate resonant states at a similarenergy.
Since there are no other existing experimental data on THF, we have performed
calculations for cyclopentane (C5 H10 ), which is very similar to THF except that here
the oxygen atom is substituted with a -CH2 - group. For this molecule, experimental
data [92] show a broad shape resonance at 8eV. The similarity of the two moleculessuggests that many features of the cross section will be analogous. The results arecompared to THF in Fig. 5.17. The cross section closely resembles that of THF, with atriple resonance superimposed on a broad higher energy resonance centered around 14
eV, which produces a wide plateau. The electron scattering cross sections for alkanes
and cycloalkanes of different sizes have very similar behaviors, [92, 178] namely thesame resonant structures and a magnitude dependent on the number of carbon atoms.
Our evidence thus suggests that THF, being very similar to a cycloalkane, shares thesecommon features, as suggested also in Ref. [209]. Also the cross section for cyclopentane
is, in our calculations, larger than in THF, a trend that seems to be confirmed by the125experimental data [178, 209]. The cyclopentane cross section peak is roughly 48 Å2which is close to the experimental value [92] for the total cross section. Recent theoretical
calculations [196], performed on a THF molecule distorted to achieve a planar geometry
for the furanose ring, show results similar to ours, with a double shape resonance in asimilar energy range.
As in Ref. [194] we attempt to correlate the spatial shape of the resonant wavefunction with the Hartree-Fock virtual orbitals obtained with a small (6-31G**) basis
set. In this case the analysis is complicated by the fact that the molecules are not planar,
so a projection on a two-dimensional surface would not pass through all the nuclei and,since the resonance is wide, different contributions mix and overlap. However, we haveestablished that the resonant wavefunction at its peak resembles closely a virtual orbitalwith energy of 7 eV (orbital 23). At higher energy (around 9.5 eV, when the second
resonant contribution becomes dominant), the resonance spatial structure appears similar to an orbital of 11 eV of energy (orbital 31). We show this in Fig. 5.18. Since there
is a great degree of resonance overlap, the relationship of the resonant wavefunctions
to the virtual orbitals we show becomes even clearer when observing the other virtual
orbitals, because they are actually very different from our wavefunctions.
Recently, dissociative electron attachment (DEA) cross sections have been measured [6] for all the DNA subunits. These cross sections display two prominent resonances at 6 and 8 eV, of width around 1 eV each for THF, which could correlate withthe calculated shape resonance, although clearly much complicated vibrational dynamics happens between electron capture and dissociation that we do not take into accounthere, and in the results of Ref. [23] there are electronically inelastic resonances in thisenergy region which could lead to similar outcomes in DEA.

126

Elastic cross section (Å)

255

cyclopentane
THF

50
45
40
35
30
25
0

2

4

6

8

10

12

14

16

18

20

22

Energy (eV)Figure 5.17: Partial fixed nuclei elastic cross sections for electron scattering by THF
and cyclopentane. The two molecular structures are similar (THF has an oxygen atom
instead of a -CH2 - group) and so are their electron scattering cross sections. Cyclopentane has essentially zero dipole moment, and therefore the low energy part of its crosssection does not rise the way it does in THF. From Ref. [195].

127

M
OOOOLLLLDDDDEEEENNNN
M
M
M
M
M
MOOOLLLDDDEEENNN

O

C

H
H

defaults used
Edge = 10.12 Contour Euclid Psi = 23

CH

4

H

1

CH

C

3

2

H

H

H

CONTOUR VALUE
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
............

0.00000
0.02500
0.05000
0.07500
0.10000
0.12500
0.15000
0.17500
0.20000
0.22500
0.25000
0.27500
0.30000
0.32500
0.35000-0.02500-0.05000-0.07500-0.10000-0.12500-0.15000-0.17500-0.20000-0.22500-0.25000MOOOOLLLLDDDDEEEENNNN
M
M
M
M
M
MOOOLLLDDDEEENNN
H

O

H C

defaults used
Edge = 10.12 Contour Euclid Psi = 31

C HH

4

CH 3
H

1

CONTOUR VALUE

C2

H

H

Figure 5.18: (Left: Time delay eigenfunctions for THF on resonance at 8eV (top, where
the widely-spaced dotted curve (labeled “2”) in Fig. 5.16 is dominant) and 9.5eV(bottom, where instead the dashed curve (labeled “1”) in Fig. 5.16 dominates). A sliceof the three dimensional eigenfunctions is shown, on the plane that contains OC 1 C4 ,while C3 is above the plane and C2 below. The red and yellow contours identify positive
areas, while the blue are negative areas of the real part of the wavefunction. Right:
virtual orbitals for energies around 7 (top) and 11 eV (bottom) from HF at the 6-31G**
level. The full lines identify positive areas and the broken lines negative areas. The
correspondence between the two top plots with each other and between two bottom
plots is very pronounced, allowing us to identify the two main contributions to theresonance. Although the molecule is not planar this projection is much easier to readthan the three dimensional structures. From Ref. [195].

128
Several electron-energy-loss (EEL) studies of THF [26, 108] and cyclopentane [7]have been performed, and some attempts have been made at assigning the resonances,also in relation to the position of the excited states from absorption spectroscopy. In Ref.
[26] a gas-phase EEL peak for THF at 6.6 eV is assigned as a core-excited resonance witha parent ion of Rydberg character, centered on the oxygen atom. A similar experimental
peak in cyclopentane in Ref. [7] is assigned as a shape resonance, on the basis ofa molecular orbital argument, with the electron wavefunction distributed around thewhole ring. Our calculations seem to agree with this second assignment, for two reasons:first, the resonance is very similar in character in THF and cyclopentane (see Fig. 5.17),
and therefore it seems implausible that it derives from excitations on the oxygen; second,
the resonant wavefunctions in Fig. 5.18 are found to be delocalized over the entire ringas had been speculated in Ref. [7] for cyclopentane. At higher energies core-excitedcontributions are most likely predominant, as was suggested in Ref. [108].
It is difficult to attempt a prediction of what the evolution of the resonances
found in elastic scattering could be, when the nuclear dynamics is included. At present
this incorporation of nuclear motion would be computationally prohibitive, but a first
analysis can be attempted by inspecting our calculated results. From the spatial structures of the resonant wavefunctions it seems likely that the resonances could initiate abreak-up of the ring, because of the presence of many nodal surfaces cutting throughthe ring bonds. For example, the experimental DEA cross sections for deoxyribose inRef. [153] exhibit a peak at around 6 eV for production of the anion C3 H5 O−
3 , whichcould be generated by a resonant structure that shows a node cutting through the ringfrom C3 to C1 , similar to the top left plot in Fig. 5.18, given that deoxyribose hasside groups attached to C4 and C3 that could lead to production of the aforementioned

Molecule
THF
H3 PO4

Energy (eV)
8.6
14.1
7.7
12.5

Width (eV)
2.3
2.5
2.0
1.5

Partial wave
4 (70%)
5 (50%)
4 (42%)
4 (41%)

129

Table 5.2: Energies, widths and largest partial waves of the resonances discussed in the
text for the DNA backbone components. From Ref. [195].

anion.

5.4

Results: H3 PO4
In the case of phosphoric acid there are no available theoretical or experimental

data to compare with, and consequently our analysis of the results will be somewhat
more limited. The cross sections show two prominent resonances, one at 7.7 eV andthe other at 12.5 eV, with respective widths 2.0 eV and 1.5 eV. The main partial wave
contributions are l = 4 and l = 3 at 40% and 30% respectively for both resonances,
as stated in Tab. 5.2. At low energies there is a sharp rise in cross section due tothe dipole moment. Cross section and total time delay are illustrated in Fig. 5.19.
It is difficult to say why the second resonance is narrower than the first, since this is
somewhat counterintuitive. We will limit ourselves to pointing out that in a complicatedpolyatomic molecule, shape resonance widths depend not only on the total energy, butalso the size of the probability density near the escape regions in the multidimensional
potential. It is possible, therefore, that a higher-energy resonance could be narrower
than one at lower energy. Resonances broader at lower energy than the higher energy
ones can be seen, for example, in C60 fullerene [117, 206].
In Ref. [6] the DEA cross section has been measured for trimethylphosphate(PO4 (CH3 )3 ), and it displays a very wide and prominent resonant structure at 7.5 eV,

130
to which our 7.7 eV resonance could correlate. No branching ratio has been measured forthis compound. The structure of phosphoric acid presents three identical -OH groups
and therefore the resonances will show structures with a similar probability for the
scattering electron to end up in each of these groups, making for complicated spatial
profiles, which in the end have not proven to be very illuminating. For these reasons wewill not attempt here an analysis of the spatial structure of the resonances.
As in THF, the shape resonances are quite broad, so the anions will be relatively
short lived, and probably they will autoionize back (possibly with vibrational excitation)
before having the possibility to trigger a breakup or energy transfer to another DNA
subunit. However, the presence of water, a structural component in biological DNA,can act to stabilize the anions, and probably will significantly influence the lifetimesof these resonances in the cell nucleus since the backbone is in closer contact with thesurrounding environment, compared to the bases that lie inside the double helix.

5.5

OutlookThe next aspect of the DNA radiation problem we want to explore is the coherent

interaction of DNA subunits and their effect on the electron scattering/capture process.
This can be accomplished in a multiple scattering framework (see e.g. Refs. [29, 41])
using the scattering data we gathered for the single components of the macromolecule,and it is currently under way in collaboration with Profs. L. Sanche and L. Caron at
the University of Sherbrooke.

131

14
12
10
Time delay

8
6
4
2
0
−2
−4
−6
0

5

10

15

20

25

30

20

25

30

Energy (eV)90

Elastic cross section (Å)

2

80
70
60
50
40
30
20
10
0
0

5

10

15
Energy (eV)Figure 5.19: e-H3 PO4 partial cross section (top) and time-delay analysis (bottom).
Calculations again involve partial waves up to l=10 and the dipole physics outside theR-matrix box is included exactly for those partial waves. Two broad resonances arepresent, at 7.7 and 12.5 eV, the cross section is smaller than in the DNA bases andcomparable to THF. From Ref. [195].

Chapter 6

Results: dissociative recombination

6.1

Introduction
In dissociative recombination, we have studied extensively the HCO+ molecule,which is interesting for many reasons. In particular, it is an important astrochemical
species, the ion is linear and therefore the neutral exhibits Renner-Teller coupling, itpossesses conical intersections, the reaction mechanism is controversial and previous
theoretical calculations of the rate coefficient are orders of magnitude lower than theexperimental values. We illustrate these points more extensively below.
The HCO+ molecular ion is one of the key ions in the chemistry of interstellar
clouds. The major destruction channel of this ion is believed to be dissociative recombination (DR) by low-energy electrons. Therefore it is crucial to study electron collisions
with HCO+ ions. Numerous experimental studies have been carried out on this reaction using different techniques, such as stationary afterglows [110], flowing afterglows[3, 4, 106], a merged beam study [144] and, very recently, a storage ring experiment [56].
The thermal rate coefficient at 300 K is found in several studies to be about 2.0·10 −7
cm3 s−1 [106]. Its temperature dependence is stronger than the T −0.5 dependence usuallyfound when DR is driven by the direct mechanism where an excited neutral diabatic

133
state is crossing the ionic ground state close to its minimum. Therefore it has been
proposed that there is no direct mechanism for DR of HCO+ and instead the process
is driven by electron capture into ro-vibrationally excited Rydberg states [110]. The
temperature dependence of the DR rate has been measured in Ref. [144] as T −1.2 ,which should be viewed as schematic, since it is not consistent with the known analytical threshold behavior. The DR cross section has been measured in the two ion-beamexperiments [56, 144]. The cross section measured by Le Padellec et al. [144] is up to
an order of magnitude smaller than that measured by Geppert et al. [56]. Absorption
measurements [4], and recently the storage ring experiment [56], have shown a complete
dominance of dissociation into H+CO. It has been long known that the CH+O channel
is endothermic by 0.17 eV and could have an energy barrier, while the OH product is notobserved (see Ref. [144] and references therein). Recently, branching ratios have beenmeasured [57] for the DCO+ dissociative recombination, confirming that the D+CO
channel is by far the dominant one, with a ratio of 0.88.
Relevant in this context is also the series of experimental publications on autoionization of HCO Rydberg states carried out by Grant et al. [122, 123, 124, 157, 211],with which we will compare some of our results.
The DR mechanism of HCO+ is controversial. Even though the temperaturedependence of the DR rate suggests an indirect mechanism, the two existing theoreticalstudies predict different pathways. Kraemer and Hazi [99, 100] carried out complete
active space SCF (CASSCF) calculations at linear geometry on the ground state of
HCO+ and some of the excited states of HCO relevant for DR. No indications were foundof a neutral state crossing the ionic ground state near the minimum. However, in thestudies by Talbi et al. [184] carried out at linear geometry, a direct mechanism for DR of

134
HCO+ was proposed on the basis of multi-reference configuration interaction (MRCI)
calculations with localized orbitals. Using a diabatization procedure, the adiabatic
potential curves generated a repulsive neutral state crossing the ionic ground statenear the minimum. This diabatization procedure seems questionable, as discussed forexample in Ref. [18] where instead a multi-step mechanism has been proposed, in which
the electron capture and dissociation proceed through a series of ∆v = 1 transitions
between Rydberg states of HCO.
One of the main goals of this study is to establish the reaction mechanism. To
do this, accurate quantum chemistry calculations of both valence and Rydberg statesof HCO have been performed. The Renner-Teller effect is strong in the ground ( X̃ 2A0 )and first excited (Ã 2A00 ) states of HCO [67], so we need to assess its importance forthe electron capture and predissociation in DR, as the similar Jahn-Teller effect was
found to be crucial in DR of H+
3 [94, 95]. Also this system shows a conical intersection[185], which in many molecules can provide a fast route towards the reaction products.
Its influence on the dissociative recombination process needs to be understood as well.
Conical intersections are understood now to be frequently occurring [208] rather than
exceptions.

6.2

Model 1: no Renner-Teller couplings
This section illustrates the results obtained for this reaction using Model 1 as

described in Sec. 3.2-3.4. The description of the ab initio results for the electronic
potential energy surfaces is identical for the two models.

135
6.2.1

Ab initio surfaces

To calculate the potential energy surfaces of HCO+ and HCO we used the MOLPRO program package [203]. In order to calculate all the Rydberg states desired, a basis
set with enough diffuse basis functions is necessary. We started with the aug-cc-VPTZ
basis set [44] to which we added diffuse s-, p- and d- functions manually. For eachnew diffuse function, the exponent was decreased by a factor of 2.5 for the s- and pfunctions, and by a factor of 3.0 for the d- functions. From Koopmans’ theorem [181]we know that the virtual molecular orbitals of the HCO+ ion with negative energies
will approximately describe the Rydberg orbitals of the neutral HCO molecule. We
therefore added diffuse functions to the basis set until the desired Rydberg states wereconverged. In the final calculation, we added up to 9s functions, 9p functions and 5dfunctions. From the orbital energies, we concluded that up to 12 Rydberg states of A 0
symmetry and 6 Rydberg states of A00 symmetry were converged. Then we calculated
the potential energy surface of the ion, the 8 lowest states of A0 symmetry and the 4lowest states of A00 symmetry by starting with a HF calculation of the ion followed by
a MRCI calculation. The advantage of starting with a HF calculation of the ion evenwhen the neutral is studied, is that then we know, according to Koopmans’ theorem,that the MOs will be suitable for describing the Rydberg states. The calculation may
not be ideal for the valence states. The CI wavefunction is represented in a space of 4000
internal configurations and more than 1.2 million external configurations consisting of
single and double excitations of the reference wavefunctions. We used the same active
space for the calculations of the ion and neutral and therefore we believe that the relative energy between the charged and neutral states will be correct and the calculations

136
will be balanced. We need to calculate the potential energy surfaces as functions of all
three internal coordinates (rCH , rCO and θ, the bending angle) and this would be very
time-consuming using the MRCI method. We therefore once again employed Koopmans’ theorem and approximated the energies of the Rydberg states using the ionicenergy calculated with the MRCI method plus the orbital energy of the corresponding
Rydberg state.
Since the orbital energies are negative, the Rydberg potential energy surfaces will
fall at a lower energy than the ion potential. In this approach we neglect the correlation
between the Rydberg electron and the ion core. The correlations among the electrons
sitting on the ion are included since we use the correlated MRCI potential for the ion.
Comparing the potential energy surfaces calculated using the MRCI method against
those obtained from Koopmans’ theorem, the MRCI results were lower in energy because
they incorporate more electron correlation. Fig. 6.1 shows the calculated potential
energy surfaces of HCO+ and HCO at linear (θ=180 degrees) and bent (θ=150 degrees)
geometries and rCO = 2.0877 a0 . The squares are the energies calculated with the MRCI
method and the lines are from the approximation using Koopmans’ theorem.
The ion is the curve with black squares. The filled squares (solid lines) are statesof A0 symmetry, while open squares (dotted lines) are states of A00 symmetry. Noticethat there is no curve-crossing of a repulsive neutral state with the ionic ground stateclose to its minimum. The first excited neutral repulsive state crosses the ionic groundstate roughly 2 eV above threshold, at an rCH distance of roughly 4 a0 , well outside theregion important for DR at the low energies that are of interest here, therefore we haveneglected it in the following. The repulsive state is of Σ symmetry at linear geometry
and lies below the ionic curve at all rCH distances. We therefore conclude that DR of

137
HCO+ is driven by an indirect mechanism whose initial step is electron capture into theRydberg states. The effective quantum numbers are then calculated from the Rydberg
formula, Eq. 2.35.
The agreement between the quantum defects calculated using the neutral MRCI
potentials and Koopmans’ theorem is remarkably good. We see that the states show
Rydberg character (the effective quantum numbers in one series differ by 1) at least at
small CH distances. Notice that some effective quantum numbers are close to 3.0. These
are assumed to be d-states with almost zero quantum defects. The effective quantum
numbers for linear and bent geometries are also shown in Fig. 6.1.
The next step was to fit the ion potential. The rCH and rCO dependencies were
fit to Morse functions, and linear functions were used to fit the rCH dependence of theMorse potential parameters. The ion-potential was fitted between 1.0 a 0 ≤ rCH ≤ 10.0
a0 , 90 degrees ≤ θ ≤ 180 degrees and for five rCO distances (1.7877, 1.9, 2.0877, 2.3877
and 2.55 a0 ). Then the quantum defects were fit to a model in which the σ levels haveoff-diagonal couplings, while the π part of the matrix is diagonal ( Sec. 3.2). For the
diagonal elements a fourth-order Taylor expansion was used, whereas the couplings wereassumed to have a Gaussian form.

6.2.2

Vibrational dynamics

To calculate the vibrational frequencies of the HCO+ ion, we have used the potential energy surface calculated by Palmieri et al. [154], since this was better suited todescribe the ground state of the ion than the one we calculated. This is because ourab initio calculations are optimized for the Rydberg states energies, not just for the
ground state surface. The difference between the two calculations is not very large, asE(au)

ν138

r CH(au)r CH(au)Figure 6.1: Potential energy curves (bottom) and quantum defects (top) for HCO. The
dots are MRCI calculations, while the lines are calculated using Koopmans’ theorem
for the Rydberg states. The black line in the potential plots is the ground state of theion. The solid lines and filled dots represent A’ states, whereas the dashed lines andempty dots are A” states. The left side of the figure has been calculated at a fixed angle
θ = 180 degrees (equilibrium value for the ion), while the right side has θ = 150 degrees.
In both cases rCO is fixed at its equilibrium value of 2.0877 a0

139
Table 6.1: Comparison of vibrational frequencies (in cm−1 ) for vibrational transitions in
HCO+ , calculated using the Palmieri ionic potential surface [154]. The second column
reports values from Ref. [154], while the third is an adiabatic calculation and the fourth
is a full three-dimensional diagonalization. The last column shows the frequencies obtained in an adiabatic framework using the potential we calculated, which we have used
only for the Rydberg states in this work and not to obtain ionic vibrational frequencies,
since the Palmieri potential was better suited for this purpose.
(νCH , νθ , νCO )
020
001
100
040
021
002
120

νP almieri
1641
2182
3090
3256
3822
4345
4679

νAdiabatic
1630
2195
2930
3188
4087
4514

ν3D
1629
2189
3082
3238
3815
4363
4658

νAdiab. Larson potential
1739
2299
3150
3460
4044
4678

is shown in Tab. 6.1 where also a comparison between the adiabatic model and full
three-dimensional calculations is shown. Since the coordinates mimic quite closely thenormal modes of the molecules, the nonadiabatic couplings are limited and the adiabatic
frequencies are largely correct for the lowest vibrational states.
As was mentioned in Ch. 3 we use Jacobi coordinates for the vibrational dynamics, and we only consider one fragmentation channel, due to the H+CO channel being
largely dominant, as stated in Sec. 6.1.

6.2.3

Autoionization widths

To test the wave functions and quantum defects used in the present study wecompare calculated autoionization widths with the measurements of Grant et al. [122,123, 124, 157, 211] using high resolution double-resonance spectroscopy, as a function of
the vibrational normal modes CO stretching and bending [123]. The vibrational wave
functions for the C-H stretch have either zero boundary conditions on the left (box

140
Table 6.2: Calculated energies and widths of the autoionizing HCO resonances. Here nis the approximate principal quantum number and (νCH , νθ , νCO ) is the ionic vibrationalmode of the resonance. Γbox is the width calculated using box-states, while ΓSiegert iscalculated using Siegert states. The widths measured by Grant et al. [123] are displayedas Γexp . From Ref. [104].
Energy (H)8.65·10−5
3.13·10−3
2.12·10−4
2.30·10−3
6.50·10−3
2.30·10−4
7.53·10−4

n
6
7
7
8
12
12
13

(νCH , νθ , νCO )(1,0,0)(1,0,0)(0,0,1)(0,0,1)(1,0,0)(0,1,0)(0,1,0)Γbox (cm−1 )
1
3
0.08
0.026
0.018
0.013
0.012

ΓSiegert (cm−1 )
3
6
2
1
0.7
0.6
0.45

Γexp (cm−1 )
10
8
2.5
1.5
1.25

states, which describe only autoionization processes) or Siegert boundary conditions(the Siegert states described above, which allow for escaping flux due to predissociation).
The resulting widths are displayed in Table 6.2. The measured widths are up to an orderof magnitude larger than the calculated widths. The widths calculated using Siegert
states are also an order of magnitude larger than those calculated with box states.
This seems to indicate that predissociation is a determining factor in the autoionization
dynamics. Notice also that the largest autoionization widths we calculate are for the
CH stretch, for which no experimental data are available. For the other two normal
modes the widths are comparable between each other (larger for CO stretching thanbending, similarly to experiment), but they are very sensitive to predissociation, while
this is true to a lesser extent for the CH stretching mode.

6.2.4

Electron scattering and cross section

In this first calculation, which includes no Renner-Teller coupling, we have used
the Siegert states method (Sec. 3.4) for calculating the DR cross section. The results

141
are converged when we include 40 adiabatic states and 9 C-H Siegert states for each
adiabatic potential. To converge these states, we need about 40 grid points in r and ϑ
and 100 points for R. The resulting DR cross section is shown in Fig. 6.2.
It should be noted that the calculated cross section is overall smaller than the
cross section measured by Le Padellec et al. [144], especially for E ∼ 0.01 eV where thetheoretical cross section drops several orders of magnitude. We have tried to estimate theeffect of nonadiabatic couplings on the DR rate, using the complex absorbing potential[161] and complex scaling [133] methods. In both cases the effects have turned out tobe quite small, therefore the adiabatic approximation is quite good for this system.

6.3

Model 2: inclusion of Renner-Teller couplings
In this section we show results for Model 2, described in Sec. 3.5-3.6. The

main differences with respect to the previous model are the inclusion of Renner-Teller
couplings, freezing of CO motion, and a different method to calculate the DR crosssection.
We have calculated DR cross sections for two different initial vibrational states0χ+v (RCH ) of the ion: the ground state v = 0 (which can be viewed as {00 0}) and thefirst excited state v = 1 (which can be viewed as {011 0}) whose energies are represented
in Fig. 3.5 with horizontal dotted lines. The energy difference between the two states
is about 0.1 eV. Figure 6.3 shows the two cross sections with dot-dashed (v = 0) anddashed (v = 1) lines. Fig. 6.3 shows that the DR cross section depends strongly on
the initial vibrational state. A similarly strong dependence has also been observed inexperiments with other molecular ions. Assuming that the initial experimental vibrational state distribution is in thermodynamic equilibrium, it is necessary to average the142

Figure 6.2: Dissociative recombination cross section obtained with the Siegert states
method, using the model with no Renner-Teller couplings. Calculated (full line) andmeasured (circles) [144] cross sections for DR of HCO+ . To compare with the measured
cross section, the calculated cross section has been convoluted with a 1 meV (dottedline) and 3 meV (dashed line) Gaussian energy distribution, this smoothes a large partof the Rydberg resonances that are not seen in experiment, since the resolution is toolow. Here we do not consider rotational structure, which would increase the number ofresonances. From Ref. [104].

143

experiment by Le Padellec et al.
calculation by Larson et al.
ground mϕ=0 vibr. state

-13

10

2

DR cross section (cm )excited |mϕ|=1 vibr. st
-14

10

-15

10

-16

10

-17

10

0.001

0.01

Energy (eV)0.1

Figure 6.3: The figure shows the present theoretical DR cross section (dashed and dotdashed lines) for HCO+ as a function of the incident energy E. The initial vibrational
state is the ionic ground state for the dot-dashed curve and the first excited state forthe dashed curve. The experimental [144] (cross symbols) and previous theoretical [104](thin solid line) cross sections are also shown for comparison. The theoretical curves
include the averaging over the electron energy distribution according to the proceduredescribed by Eq. (2) of Ref. [96] with ∆E⊥ = ∆E|| = 3 meV. From Ref. [130].

144
cross section over the thermal distribution to compare with the experiment. The firstexcited vibrational state has vibrational symmetry different from the ground state. The
vibrational angular momenta for the states are mϕ = 1 and 0. Thus, the deexcitation
process may be too slow to reach equilibrium, resulting in a vibrational temperaturethat is higher than the electron temperature (300K in the experiment of Ref. [144] ).
The vibrational temperature in the experiment of Le Padellec et al. [144] is not known.
If it was in equilibrium with the electron temperature at 300K, the contribution fromthe excited vibrational states would be small. But for a larger vibrational temperature,for example 1000 K, the averaged cross section would be about a factor of 1.5 larger thanthe cross section for the ground vibrational state in Fig. 6.3. The fact that in several
different experiments the measured DR rate ranges over values from 0.65 − 3 × 10 −7
cm3 /s at 300 K (see Ref. [3, 56, 106, 110, 158] and Fig. 6.4 below) might conceivably
derive from differences in the initial vibrational populations. Determination of the actual experimental vibrational distribution and/or controlling it in HCO+ could be an
important step in understanding DR in small polyatomic ions.
The present theoretical DR cross section is approximately a factor of 2 smaller
than the experimental data. Assuming that this reflects a limitation of the present theoretical description, this might derive from our approximation that freezes the CO bondlength. However, the expected error of order 70% caused by our adiabatic approximation in the RCH coordinate should also be kept in mind when assessing the implicationsof this discrepancy. Still, if the CO bond is allowed to vibrate, which makes the ion more
floppy, the probabilities of capturing and predissociation will presumably be increased.
Quantitatively, releasing the CO bond increases the density of HCO resonance statesthat can be populated in the electron-ion collision and, correspondingly, the sum in Eq.

145
3.26 is expected to increase.
In afterglow plasma experiments, the measured observable is the DR rate constantthermally averaged over the kinetic energy distribution of colliding electrons and ions,α(kT ) =< vσ >, where v is the velocity of the electron, is obtained from the DR
cross section shown in Fig. 6.3 according to Eq. (7) of Ref. [94]. The resulting
theoretical rates (thermally averaged over the electron energies, for each of the two initialvibrational states) are shown in Fig. (6.4) and compared with available experimental
measurements.
Although the theoretical DR cross section is still smaller by a factor of 2 thanthe experimental results of Le Padellec et al. [144], we find that inclusion of RennerTeller coupling significantly increases the theoretical DR cross section. The previous
theoretical study by Larson et al. [104] that omitted Renner-Teller coupling, but which
did account for the CO bond vibrations, gave a DR cross section significantly smaller
than the present result. This is evidence for the important role of Renner-Teller couplingin the HCO+ DR process.
Since the Born-Oppenheimer potential surfaces are the same for HCO (HCO + )and DCO (DCO+ ) we did a similar calculation for DR in DCO+ . The resulting DR
cross sections for the ground and first excited vibrational level are very similar to theones presented in Fig. 6.3 but smaller by a factor of 1.5. A similar dependence of the
DR rate on the isotopologue masses was found in theory and experiment for the H +3,
H2 D+ , D2 H+ , and D+
3 ions [97]. Notice that the lifetime of the first excited vibrational
state of DCO+ is supposed to be very long, of the order of some ms [204], therefore itis even more important in this case to control the initial vibrational population.

-7

10

3

thermal rate (cm /s)

146

-8

10

10

theory, ground vibr. state
theory, excited vibr. state
Le Padellec et al. 97
Poterya et al. 2005
Geppert et al. 2005
Adams et al. 84
Amano et al. 90
Ganguli et al. 90
Gougousi et al. 97
Laube et al. 98
Leu et al. 73
Rowe et al. 92

100

Te (K)

1000

Figure 6.4: Theoretical and experimental DR thermal rates for HCO + . The dot-dashedand dashed lines are the theoretical rates obtained for the ground and first excited
vibrational states of the ion. The experimental rate (solid line) is obtained from theexperimental cross section, Ref. [144]. The results from a number of other experimentswith HCO+ [3, 56, 106, 110, 158] are also shown. From Ref. [130].

147
6.3.1

Estimation of the effect of autoionization and CO vibration on the
cross section

We should mention some approximations made in the present study that canaffect the theoretical DR cross section. First, our treatment has not accounted forthe possible autoionization after the electron is captured by the ion, before the neutral
molecule has time to predissociate: in Eq. (3.20), the survival probability e−ρ(E) was
set to be 1. The inclusion of autoionization decreases the calculated DR cross section.
The effect of autoionization on the DR cross section can be estimated as follows. Figure
6.5 shows the calculated resonance curves Ua (RCH ) as described above and in Ref. [94].
Above the lowest ionic curve each of these potential curves has in general a non-zero
width Γa (RCH ), which represents the adiabatic autoionization rate. If there is morethan one open ionic channel |ji for a given Ua (RCH ), the neutral molecule can decay(i)into the i-th channel with a partial width Γa (RCH ), whereP(i)i Γa (RCH ) = Γa (RCH ).

Once the electron is captured, the system can evolve into two competing pathways:
resulting in either autoionization or dissociation. The relative probabilities per unit oftime for the two processes provide an estimate for the survival probability e −ρ(E) .
The largest partial autoionization widths are in general the widths of the resonance curves Ua (RCH ) belonging to Rydberg series associated with a few nearest closed
ionic channels. Such resonance curves are very similar in shape to their parent ion
potential curves, when the principal quantum number n is high. Consequently, in theenergy range of interest, these curves Ua (RCH ) are closed with respect to adiabatic dissociation. The dissociation can occur only through coupling to true dissociative states(i.e. the process denoted as predissociation in molecular spectroscopy). The predisso-

148

Figure 6.5: Resonance potential curves Ua (RCH ) of the neutral molecule states having
magnetic quantum number M = 1. The potential curves were calculated in the energy
range where at least one electron-ion channel Ui+ (RCH ) (for the given RCH ) is open.
Thus, the lower bound of the calculated curves corresponds to the ground potential curve+
U0,0(RCH ) of the ion. The potential curves of the ion are also shown by dashed thicklines (see also Fig. 3.5). The Rydberg series are crowded just below the ionic potential
curves. For clarity, visible gaps between the Rydberg series and the ionic curves havebeen intentionally introduced artificially, to show the behavior of the resonance curves
belonging to higher ionic channels. The horizontal line represents the energy of theground vibrational level. From Ref. [130].

149
RCH a.u.
2.056
2.11
1.93
2.02

Energy (eV)
0.112
0.144
0.34
0.283

α (a.u.)1.18 × 10−5
2.9 × 10−5
4.0 × 10−5
1.0 × 10−5

∆F (a.u.)
0.00410.01110.00360.0022∆E(a.u.)0.00590.0047
1-37×10−4
1-37×10−4

Pa0 ,a0.9 × 10−4
2 × 10−4
14 − 83 × 10−4
1.2 − 7.2 × 10−4

Table 6.3: The table demonstrates Landau-Zener parameters and probabilities for typical avoided crossings. The first two examples correspond to avoided crossings situateddeeply below the energy of the ground vibrational level. The third and fourth examples
correspond to energies around and above the ground vibrational level. Thus, ∆E could
be very different for the last two examples giving different probabilities. In general, for
any total energy of the system there are always some avoided crossings with small ∆E
at the corresponding turning points. From Ref. [130].

ciation probability can be estimated using the Landau-Zener model. When the neutral
molecule is vibrating along the Ua (RCH ) curve, every time it passes through an avoided
crossing with another curve Ua0 (RCH ), it can jump to the corresponding state |a0 i via
an adiabatic transition. We view such an adiabatic transition as the pathway to predissociation. The probability for adiabatic passage through an avoided crossing is givenby the Landau-Zener formula [101]Ã

2πα2Pa0 ,a = 1 − exp −
∆F

r

m
2∆E

!

≈

2πα2
∆F

r

m
2∆E

(6.1)

where α is the non-diagonal coupling element between diabatic states, i.e. the states thatwould cross if the coupling was absent. Numerically, α is equal to half of the adiabatic
potential curve splitting at the avoided crossing. ∆F is the absolute difference in slopes(net classical force) between crossing diabatic potential curves, ∆E is the classical kinetic
energy at the avoided crossing, m is the reduced mass of the system. The resonancecurves shown in Fig. 6.5 demonstrate numerous avoided crossings.
Table 6.3 indicates the order of magnitude of the parameters in Eq. 6.1. This
table shows two examples of avoided crossings that lie well below the energy at which

150
the electron can be captured. At such large values of ∆E, the Landau-Zener probability
is very small ( 10−4 ) and depends only weakly on the electron energy. The other two
examples are taken from the region around the left turning point where the velocity
of motion can be small or large depending on the energy of the electron. In cases like
this, we calculated the probability Pa0 ,a for a range of values of the kinetic energy ∆E.
There are many such avoided crossings near left turning points and comparatively few
at right turning points.
During half (τ1/2 ) of an oscillation period, the system goes through nc avoided
crossings. We estimate that nc is about 10 for a typical resonance curve. Therefore, thetotal predissociation probability Pd is somewhere in the range 0.001-0.01. In the above
estimation of the predissociation probability we have assumed that only two statesinteract at each avoided crossing. However, in reality many of the avoided crossing
cannot be considered as strictly two-state ones: three or more states |ai can participate.
The autoionization probability Pion during time τ1/2 is Pion ∼ τ1/2 Γ, where Γ isthe total autoionization width. The largest Γ is about 10−5 a.u., τ1/2 ∼ π/ω, whereω/2 is the frequency of oscillations in the ground vibrational level, h̄ω/2 = 0.007 a.u.
Thus, Pion ∼ πΓ/ω ≈ 0.002. This is presumably an upper bound on the autoionization
probability.
As one can see, the autoionization probability could be competitive with predissociation in our model if we take the smallest ratio Pd /Pion ∼ 0.5, but it is more
likely that the ratio is of the order of 2-3 and therefore, more favorable to predissociation. In the previous theoretical study of DR in H+
3 [94], the ratio was significantly
larger. An important difference with the H+
3 study is the two degrees of freedom taken
into account when adiabatic potentials Ui+ were calculated. As a result, the density of

151

0.13
0.128

Energy (eV)

0.126
0.124
0.122
0.12
0.118
0.116

1.9

1.92

1.94

1.96

1.98

RCH (a.u.)

2

2.02

2.04

2.06

Figure 6.6: A part of the spectrum of Fig. 6.5 showing the interval of RCH corresponding
to the left turning point of the vibrational motion. For convenience, the energy of the+ground ionic curve U0,0
(RCH ) was subtracted from all the data. In the figure, curves
with positive and small negative slopes correspond to Rydberg states converging to thenearest ionic thresholds. The probability to capture the electron tends to be high forsuch states, diminishing however at high principal quantum numbers as n−3
c . From Ref.
[130].

152
resonance curves was much higher (factor 20-40). The main contribution to the highdensity of resonance curves was from the states with small principal quantum number
corresponding to highly excited ionic channels. This increased significantly the probability of predissociation. In fact, in the case of HCO, the vibrating CO bond should also
increase the density of steep resonance curves in Figs. 6.5 and 6.6. This is because if
the vibration along CO is quantized, each of the curves shown in Fig. 3.5 will produce
a series of additional ionic curves with different CO quanta. Therefore, the number ofionic potentials will be larger and, correspondingly, the density of resonance curves andavoided crossings in Figs. 6.5 and 6.6 will be increased too. This should increase theDR cross section (more states to be captured to) and the ratio P d /Pion (more avoided
crossings).
Another approximation is that the n = 2 Rydberg states have dissociative character and are poorly described by this model. However, we believe that the dominant
DR pathways are triggered by an initial capture into higher Rydberg states, rather thanbeing directly captured into n = 2 states.

6.3.2

Comparison with our previous theoretical study

In our previous theoretical study of DR in HCO+ , published in Ref. [104], theRenner-Teller coupling was not taken into account because it appeared to have a small
effect on the potential surfaces of excited molecular states of the HCO molecule. However, Renner-Teller coupling involves degenerate π states and therefore should play animportant role in electron-ion scattering especially when pπ and pσ molecular potentials approach closely to each other. Although Renner-Teller coupling was not includedin the previous study, some nonadiabatic effects were considered there. The quantum

153
defects µ(Q) for sσ, pσ, dσ states obtained from ab initio calculation demonstrate sharp
variations at some configurations Q of the three HCO nuclei. Therefore, instead of using the three adiabatic defects µ(Q), a numerical diabatization procedure was applied
to give a Q-dependent 3x3 matrix µ(d) (Q) of quantum defects. The Q dependence ofthe matrix is weaker and the matrix has more information about nonadiabatic effects
than the diagonal matrix µ(Q). On the other hand the nonadiabatic effects represented
by non-diagonal elements of µ(Q) might have a smaller effect on the dynamics of theelectron compared to the nonadiabatic couplings involving degenerate π states. Indeed,
the DR cross section obtained in the previous study [104] is much smaller than the crosssection of the present treatment. To verify that the increase in the DR cross section isindeed due to the Renner-Teller coupling, we artificially set the non-diagonal elementsin the Hamiltonian of Eq. 3.12 to zero and calculated again the cross section for M = 0.
The resulting DR cross section is about a factor of three smaller than the cross section
shown in Fig. 6.3 and on average it is close to the result of the previous study. It is
worth mentioning that the calculation of Ref. [104] accounts for all three vibrational
degrees of freedom but the very high vibrational levels that generate the steep Rydberg
n = 2 curves in Fig. 6.5 (which have large widths and therefore lead to an increasein the capture probability) were omitted; these two effects would tend to balance eachother, and might give roughly the same cross section as in the present study when offdiagonal couplings are switched off. We speculate that if the CO bond vibration was
included in the present treatment the effect of the Renner-Teller coupling would likely
be even more pronounced. To summarize our discussion of this point, we conclude that
the Renner-Teller effect appears to be the most important nonadiabatic coupling in the
dissociative recombination of HCO+ .

154
The theoretical cross section obtained in this study is inversely-proportional tothe incident energy below 0.1 eV (see Fig. 6.3). Above 0.1 eV the cross section falls
more rapidly with increasing energy. This is caused by the fact that at 0.1 eV an
additional ionic channel becomes open. That ionic channel is responsible for the large
probability to capture the electron. Above that ionization threshold, the density of
Rydberg states (in which capture can occur) is much lower, as can be seen from Fig.
6.5, and therefore the capture probability drops. Moreover, resonance states above thationization threshold have another competing decay channel for autoionization, which
steals flux away from the DR observable. For the ground vibrational level {00 0 v3 } (dotdashed line in Fig. 6.3 ), that additional open channel is {011 v3 }, for the first excitedlevel there are two nearby levels: {020 v3 } and {022 v3 }. A similar behavior of the crosssection was observed in DR of H+
3 [94, 95]: the cross section drops down once a newionic channel becomes open. The cross section from the previous theoretical study of
HCO+ (the solid line in Fig. 6.3) does not show that behavior. It behaves smoothly
as the energy crosses the new ionization threshold. This is because the Renner-Teller
effect was not accounted for: hence transitions of the type {000 v3 } ←→ {011 v3 } areforbidden by the different symmetries of those vibrational states.
Another difference with the previous study is the manner in which the final crosssection is calculated. In Ref. [104], the scattering matrix Si,i0 for electron-ion collisions
was explicitly calculated, but only in the electron-ion subset of channel space with nochannel indices explicitly referring to dissociation. All of the channels |ii are vibronic
states. However, some of the channels |ii are open for dissociation, and as a result, theSi,i0 matrix is not completely unitary. The defect from unitarity was used to calculate
the dissociative flux and the corresponding DR cross section. The resulting cross section

155
from this method has a rich structure with many Rydberg series of resonances. The
cross section was then averaged over an appropriate thermal distribution of electronenergies. The averaged cross section still has a number of resonances that are not seenin the comparatively low resolution experiment. In the present paper we have employeda different averaging procedure, as was discussed above. This procedure gives a verysmooth curve for the cross section.

6.4

Outlook
The framework we have developed for dissociative recombination reactions forsmall molecules is quite general and applicable to a wide variety of targets. It is definitely
possible to improve upon these calculations, especially in the calculation of electronicparameters like Rydberg state potential energy surfaces and couplings. To this end weare currently working at calculating these quantities from R-matrix calculations, which
will also allow a better understanding of the problem and avoid the ambiguities related
to the choice of the model for the fit, which have plagued us at the beginning of the
HCO+ study.

Chapter 7

SummaryThis thesis has dealt with the theoretical description of various chemical processes
triggered by electrons in polyatomic molecules. As we have pointed out, the different
reaction environments require different approaches. We have developed, in Ch. 2, aframework that, though approximate, enabled us to make inroads in the problem of
scattering of an electron from a general polyatomic molecule in Ch. 4 and to give someinsight on the problem of radiation damage to DNA (see Ch. 5), of interest to thebiomedical community. In Chapters 4 and 5 we have shown how our approach is able
to determine scattering observables for a wide variety of large polyatomic molecules,
including CO2 , SF6 , the DNA bases, and tetrahydrofuran, many of which have notbeen tackled before by theory because they are so large.
In the part of this thesis that deals with dissociative recombination reactions,
we have first developed a theoretical framework that is quite robust and able to deal
with general small molecules in Ch. 3. Through the use of ab initio calculations, Rmatrix methods, and then quantum defect theory (QDT) we deal with the vibronicinteractions that arise in triatomics subject to conical intersections or Renner-Teller
effect. This approach is then applied to DR of the molecular ion HCO+ in Ch. 6
and the results are encouragingly close to available experimental data, although due to157
the complication of including many degrees of freedom some have been frozen and theresults can certainly be further improved.
Whenever possible we have given an outlook on the future developments that
could arise from this work, for example on the description of high harmonic generationfrom molecules, which is well under way.

Bibliography[1] Computer Physics Communications library. http://www.cpc.qub.ac.uk.
[2] H. Abdoul-Carime, M. A. Huels, F. Bruning, E. Illenberger, and L. Sanche. J.
Chem. Phys., 113:2517, 2000.
[3] N. G. Adams and L. M. Babcock. Astrophys. J., 434:184, 1984.
[4] N. G. Adams, C. R. Herd, M. Geoghegan, D. Smith, A. Canosa, J. C. Gomet,
B. R. Rowe, J. L. Queffelec, and M. Morlais. J. Chem. Phys., 94:4852, 1991.
[5] K. Aflatooni, G. A. Gallup, and P. D. Burrow. J. Phys. Chem. A, 102:6205, 1998.
[6] K. Aflatooni, A. Scheer, and P. D. Burrow. Phys. Rev. Lett., submitted, 2006.
[7] M. Allan and L. Andric. J. Chem. Phys., 105:3559, 1996.
[8] M. Allan and T. Skalicky. J. Phys. B, 36:3397, 2003.
[9] H. A. Jahn anad E. Teller. Proc. Roy. Soc. A, 161:220, 1937.
[10] A. Arthurs and A. Dalgarno. Proc. Roy. Soc. (London Ser.) A, 256:540, 1960.
[11] H. Bachau, E. Cormier, P. Decleva, J. E. Hansen, and F. Martin. Rep. Prog.
Phys., 64:1815, 2001.
[12] N. R. Badnell and M. J. Seaton. J. Phys. B, 32:3955, 1999.
[13] S. Baker, J. S. Robinson, C. A. Haworth, H. Teng, R. A. Smith, C. C. Chirila,
M. Lein, J. W. G. Tisch, and J. P. Marangos. Science, 312:424, 2006.
[14] R. Balog and E. Illenberger. Phys. Rev. Lett., 91:213201, 2003.
[15] K. L. Baluja, P. G. Burke, and L. A. Morgan. Comp. Phys. Comm., 27:299, 1982.
[16] J. N. Bardsley. J. Phys. B, 1:365, 1966.
[17] R. Barrios, P. Skurski, and J. Simons. J. Phys. Chem. B, 106:7991, 2002.
[18] D. R. Bates. J. Phys. B, 25:5479, 1992.
[19] K. J. Bathe. Finite Element Procedures. Prentice-Hall, Englewood Cliffs, N. J.,
1996.

159
[20] I. B. Bersuker. Chem. Rev., 101:1067, 2001.
[21] Bioengineering
Institute,
The
University
ofAuckland,New
Zealand.
Fem/bemnotes,
February
1997.
http://www.bioeng.auckland.ac.nz/cmiss/fembemnotes/fembemnotes.pdf.
[22] G. M. Blackburn and M. J. Gait. Nucleic Acids In Chemistry And Biology. Oxford
University Press, Oxford, 1996.
[23] D. Bouchiha, J. D. Gorfinkiel, L. G. Caron, and L. Sanche. J. Phys. B, 39:975,
2006.
[24] B. Boudaiffa, P. Cloutier, D. Hunting, M. A. Huels, and L. Sanche. Science,
287:1658, 2000.
[25] B. Boudaiffa, P. Cloutier, D. Hunting, M. A. Huels, and L. Sanche. Radiat. Res.,
157:227, 2002.
[26] L. J. Bremner, M. G. Curtis, and l. C. Walker. J. Chem. Soc. Faraday Trans.,
87:1049, 1991.
[27] J. P. Burke, Jr., C. H. Greene, and J. L. Bohn. Phys. Rev. Lett, 81:3355, 1998.
[28] Christian Buth, Robin Santra, and Lorenz S. Cederbaum.
119:7763, 2003.

J. Chem. Phys.,[29] L. G. Caron and L. Sanche. Phys. Rev. Lett., 91:113201, 2003.
[30] CERN. Cernlib. http://wwwasd.web.cern.ch/wwwasd/cernlib/.
[31] D. M. Chase. Phys. Rev., 104:838, 1956.
[32] C. W. Clark. Phys. Rev. A, 16:1419, 1977.
[33] C. W. Clark. Phys. Rev. A, 20:1875, 1979.
[34] V. Cobut, Y. Fongillo, J. P. Patau, T. Goulet, M.-J. Fraser, and J.-P. Jay-Gerin.
Radiat. Phys. Chem., 51:229243, 1998.
[35] R. Colle and A. Salvetti. J. Chem. Phys., 19:1404, 1983.
[36] P. B. Corkum. Phys. Rev. Lett., 71:1994, 1994.
[37] R. Curik and F. A. Gianturco. J. Phys, B, 35:717, 2002.
[38] M. Cek, M. Thoss, and W. Domcke. Phys. Rev. A, 70:125406, 2004.
[39] J. L. Dehmer. J. Chem. Phys., 54:4496, 1972.
[40] J. L. Dehmer. ACS Symposium Series, 263:139, 1984.
[41] D. Dill and J. L. Dehmer. J. Chem. Phys., 61:692, 1974.
[42] D. Dill and J. L. Dehmer. Phys. Rev. A, 21:85, 1980.

160[43] O. Dolgounitcheva, V.G. Zakrzewski, and J.V. Ortiz. Chem. Phys. Lett., 307:220,
1999.
[44] T. H. Dunning. J. Chem. Phys., 90:1007, 1989.
[45] B. D. Esry, C. H. Greene, and J. P. Burke, Jr. Phys. Rev. Lett., 83:1751, 1999.
[46] I. I. Fabrikant and H. Hotop. Phys. Rev. Lett., 94:063201, 2005.
[47] U. Fano. J. Opt. Soc. Am., 65:979, 1975.
[48] U. Fano and C. M. Lee. Phys. Rev. Lett., 31:1573, 1973.
[49] S. Feuerbacher and R. Santra. J. Chem. Phys., 123:194310, 2005.
[50] S. Feuerbacher and R. Santra. private communication, 2005.
[51] S. Feuerbacher, T. Sommerfeld, R. Santra, and L. S. Cederbaum. J. Chem. Phys.,
118:6188, 2003.
[52] G. A. Fiete and E. J. Heller. Rev. Mod. Phys., 75:933, 2003.
[53] A. I. Florescu-Mitchell and J. B. A. Mitchell. Phys. Rep., 430:277, 2006.
[54] A. F. Fuciarelli and J. D. Zimbrick, editors. Radiation Damage in DNA:
Structure/Function Relationships at Early Times. Battelle, Columbus, OH, 1995.
[55] B. C. Garrett and D. A. Dixon. Chem. Rev., 105:355, 2005.
[56] W. D. Geppert. unpublished, 2005.
[57] W. D. Geppert, R. Thomas, A. Ehlerding, J. Semaniak, F. Osterdahl, M. af Ugglas, N. Djuric, A. Paald, and M. Larsson. Faraday Disc., 127:425, 2004.
[58] F. A. Gianturco and R. R. Lucchese. J. Chem. Phys., 108:6144, 1998.
[59] F. A. Gianturco and R. R. Lucchese. J. Chem. Phys., 114:3429, 2001.
[60] F. A. Gianturco and R. R. Lucchese. J. Chem. Phys., 120:7446, 2004.
[61] F. A. Gianturco, S. Meloni, P. Paioletti, R. R. Lucchese, and N. Sanna. J. Chem.
Phys., 108:4002, 1998.
[62] F. A. Gianturco and J. A. Rodriguez-Ruiz. Phys. Rev. A, 47:1075, 1993.
[63] F. A. Gianturco and T. Stoecklin. J. Phys. B, 29:3933, 1996.
[64] J. C. Gibson, M. A. Green, K. W. Trantham, S. J. Buckman, P. J. O. Teubner,and M J Brunger. J. Phys. B, 32:213, 1999.
[65] A. Giusti. J. Phys. B, 13:3867, 1980.
[66] S. Gohlke, H. Abdoul-Carime, and E. Illenberger. Chem. Phys. Lett., 380:595,
2003.

161
[67] E. M. Goldfield, S. K. Gray, and L. B. Harding. J. Chem. Phys., 99:5812, 1993.
[68] G. H. Golub and C. F. Van Loan. Matrix Computation. Johns Hopkins University
Press, Baltimore, MD, 1996.
[69] J. D. Gorfinkiel, L. A. Morgan, and J. Tennyson. J. Phys. B, 35:543, 2002.
[70] A. Grandi, F. A. Gianturco, and N. Sanna. Phys. Rev. Lett., 93:048103, 2004.
[71] A.
A.
Granovsky.
PC
GAMESS
http://classic.chem.msu.su/gran/gamess/index.html.

version

7.0.

[72] C. H. Greene. In J.S. Briggs, H. Kleinpoppen, and H.O. Lutz, editors,
Fundamental Processes of Atomic Dynamics. Plenum, New York, 1988.
[73] C. H. Greene. private communication, 2005.
[74] C. H. Greene, M. Aymar, and E. Luc-Koenig. Rev. Mod. Phys., 68:1015, 1996.
[75] C. H. Greene, U. Fano, and G. Strinati. Phys. Rev. A, 19:1485, 1979.
[76] C. H. Greene and C. Jungen. Phys. Rev. Lett., 55:1066, 1985.
[77] C.H. Greene and C. Jungen. Adv. Atom. Mol. Phys., 21:51, 1985.
[78] J. Gu, Y. Xie, and H. F. Schaefer III. J. Am. Chem. Soc., 127:1053, 2005.
[79] S. L. Guberman and A. Giusti-Suzor. J. Chem. Phys., 95:2602, 1991.
[80] E. L. Hamilton and C. H. Greene. Phys. Rev. Lett., 89:263003, 2002.
[81] S. Hara. J. Phys. Soc. Jpn., 27:1009, 1969.
[82] N. A. Harris and C. Jungen. Phys. Rev. Lett., 70:2549, 1993.
[83] H.R.Drew, R.M.Wing, T.Takano, C.Broka, S.Tanaka,
R.E.Dickerson. Proc.Natl.Acad.Sci.U.S.A., 78:2179, 1981.

K.Itakura,and[84] W. M. Huo and D. Brown. Phys. Rev. A, 60:295, 1999.
[85] W.M. Huo, M. A. P. Lima, T. L. Gibson, and V. McKoy. Phys. Rev. A, 36:1642,
1987.
[86] J. Itatani1, J. Levesque, D. Zeidler, H. Niikura, H. Pepin, J. C. Kieffer, P. B.
Corkum, and D. M. Villeneuve. Nature, 432:867, 2004.
[87] J. Gilbert J. W. Demmel and X. S. Li. Superlu users’ guide. Technical Report
UCB/CSD-97-944, EECS Department, University of California, Berkeley, 1997.
[88] K. D. Jordan and F. Wang. Ann. Rev. Phys. Chem., 54:367, 2003.
[89] C. Jungen and S. C. Ross. Phys. Rev. A, 55:R2503, 1997.
[90] R. E. Kennerly. Phys. Rev. A, 21:1876, 1980.

162[91] R. E. Kennerly, R. A. Bonham, and M. Mc Millan. J. Chem. Phys., 70:2039, 1979.
[92] M. Kimura, O. Sueoka, A. Hamada, and Y. Itikawa. Adv. Chem. Phys., 3:537,
2000.
[93] W. Kohn and L. J. Sham. Phys. Rev., 140:A1133, 1965.
[94] V. Kokoouline, B. D. Esry, and C. H. Greene. Nature, 412:891, 2001.
[95] V. Kokoouline and C. H. Greene. Phys. Rev. A, 68:012703, 2003.
[96] V. Kokoouline and C. H. Greene. J. Phys.: Confer. Ser., 4:74, 2005.
[97] V. Kokoouline and C. H. Greene. Phys. Rev. A, 72:022712, 2005.
[98] H. Koppel, W. Domcke, and L. S. Cederbaum. J. Chem. Phys., 74:2945, 1981.
[99] W. P. Kraemer and A. U. Hazi. In Molecular Astrophysics. State of Art andFuture Directions. Proceedings of the NATO Advanced Research Workshop, page575. 1985.
[100] W. P. Kraemer and A. U. Hazi. In Dissociative Recombination: Theory,
Experiment and Applications, page 62. World Scientific, Singapore, 1988.
[101] L. D. Landau and E. M. Lifshitz. Quantum mechanics: non-relativistic theory.
Butterworth Heinemann, Burlingto, MA, 2003.
[102] N. F. Lane. Rev. Mod. Phys., 52:29, 1980.
[103] C. M. Laperle, J. E. Mann, T. G. Clements, and R. E. Continetti. Phys. Rev.
Lett.,, 93:153202, 2004.
[104] A. Larson, S. Tonzani, R. Santra, and C. H. Greene. J. Phys. Conf. Ser., 4:148,
2005.
[105] M. Larsson. In C. Y. Ng, editor, Photoionization and photodetachment, page 693.
World Scientific, Singapore, 2000.
[106] S. Laube, A. Le Padellec, O. Sidko, C. Rebrion-Rowe, J. B. A. Mitchell, and B. R.
Rowe. J. Phys. B, 31:2111, 1998.
[107] C. Lee, W. Yang, and R. G. Parr. Phys. Rev. B, 37:785, 1988.
[108] M. Lepage, S. Letarte, M. Michaud, F. Motte-Tollet, M.-J. Hubin-Franskin,
D. Roy, and L. Sanche. J. Chem. Phys., 109:5980, 1998.
[109] H. LeRouzo and G. Raseev. Phys. Rev. A, 29:1214, 1984.
[110] M. T. Leu, M. A. Biondi, and R. Johnsen. Phys. Rev., 8:420, 1973.
[111] P.L. Levesque, M. Michaud, and L. Sanche. Nucl. Instr. and Meth. in Phys. Res.
B, 225:230, 2003.
[112] H. J. Levinson, T. Gustafsson, and P. Soven. Phys. Rev. A, 19:1089, 1979.

163[113] M. Lewenstein, Ph. Balcou, M. Yu. Ivanov, A. L’Huillier, and P. B. Corkum. Phys.
Rev. A, 49:2117, 1994.
[114] S. G. Lias, R. D. Levin, and S. A. Kafafi. Ion energetics data. In P.J. Linstrom
and W.G. Mallard, editors, NIST Chemistry WebBook, NIST Standard Reference
Database Number 69. National Institute of Standards and Technology, Gaithersburg, MD, 2005. http://webbook.nist.gov.
[115] H. C. Longuet-Higgins. Advances in Spectroscopy, 2:429, 1961.
[116] D. Loomba, S. Wallace, D. Dill, and J. L. Dehmer. J. Chem. Phys., 75:4546, 1981.
[117] R. R. Lucchese, F. A. Gianturco, and N. Sanna. Chem. Phys. Lett., 305:413, 1999.
[118] R. R. Lucchese and V. McKoy. Phys. Rev. A, 25:1963, 1982.
[119] R. R. Lucchese and V. McKoy. Phys. Rev. A, 26:1406, 1982.
[120] M. G. Lynch, D. Dill, J. Siegel, and J. L. Dehmer. J. Chem. Phys., 71:4249, 1979.
[121] F. Martin, P. D. Burrow, Z. Cai, P. Cloutier, D. Hunting, and L. Sanche. Phys.
Rev. Lett., 92:068101, 2004.
[122] E. E. Mayer and E. R. Grant. J. Chem. Phys., 103:10513, 1995.
[123] E. E. Mayer, H. G. Hedderich, and E. R. Grant. J. Chem. Phys., 108:1886, 1998.
[124] E. E. Mayer, H. G. Hedderich, and E. R. Grant. J. Chem. Phys., 108:8429, 1998.
[125] B. J. McCall and T. Oka. Science, 287:1941, 2000.
[126] C. W. McCurdy, W. A. Isaacs, H. D. Meyer, and T. N. Rescigno. Phys. Rev. A,
67:042708, 2003.
[127] N. P. Mehta, B. D. Esry, and C. H. Greene. unpublished, 2006.
[128] H. D. Meyer, S. Pal, and U. V. Riss. Phys. Rev. A, 46:186, 1992.
[129] K. W. Meyer, C. H. Greene, and B. D. Esry. Phys. Rev. Lett., 78:4902, 1997.
[130] I. A. Mikhailov, V. Kokoouline, A. Larson, S. Tonzani, and C. H. Greene. Phys.
Rev. A, submitted, 2006.
[131] A. R. Milosavljevic, A. Giuliani, D. Sevic, M. J. Hubin-Franskin, and
B. Marinkovic. Eur. Phys. J. D, 35:411, 2005.
[132] I. Mistrik, R. Reichle, U. Muller, H. Helm, M. Jungen, and J. A. Stephens. Phys.
Rev. A, 61:033410, 2000.
[133] N. Moiseyev. Phys. Rep., 302:211, 1998.
[134] L. Morgan. Phys. Rev. Lett., 80:1873, 1998.
[135] L. A. Morgan, J. Tennyson, and C. J. Gillan. J. Phys. B, 30:4087, 1997.

164[136] L. A. Morgan, J. Tennyson, and C. J. Gillan. Com. Phys. Comm., 114:120, 1998.
[137] M. Morrison, B.C. Saha, and T. L. Gibson. Phys. Rev. A, 36:3682, 1987.
[138] M. A. Morrison, N. F. Lane, and L. A. Collins. Phys. Rev. A, 15:2186, 1977.
[139] M.A. Morrison and L. A. Collins. Phys. Rev. A, 17:918, 1978.
[140] P. Mozejko and L. Sanche. Radiat. Environ. Bioohys., 42:201, 2003.
[141] V. Ngassam, A. Florescu, L. Pichl, I.F. Schneider, O. Motapon, and A. SuzorWeiner. Eur. Phys. J. D, 26:165, 2003.
[142] T. F. O’Malley. Phys. Rev., 150:14, 1966.
[143] A. E. Orel and K. C. Kulander. Phys. Rev. Lett., 71:4315, 1993.
[144] A. Le Padellec, C. Sheehan, D. Talbi, and J. B. A. Mitchell. J. Phys. B, 30:319,
1997.
[145] R. Panajotovic, M. Kitajima, H. Tanaka, M. Jelisavcic, J. Lower, L. Campbell,
M. J. Brunger, and S. J. Buckman. J. Phys. B, 36:1615, 2003.
[146] J. M. Papanikolas, V. Vorsa, M. E. Nadal, P. J. Campagnola, J. R. Gord, andW. C. Lineberger. J. Chem. Phys., 97:7002, 1992.
[147] G. A. Parker and R. T. Pack. J. Chem. Phys., 68:1585, 1978.
[148] R. G. Parr and W. Yang. Density Functional Theory of Atoms and Molecules.
Oxford University Press, Oxford, 1989.
[149] K. Pfingst, H. T. Thummel, and S. D. Peyerimoff. J. Phys. B, 25:2107, 1992.
[150] J. A. Pople et al. GAUSSIAN98. Pittsburgh, 1998. http://www.gaussian.com.
[151] G. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling. Numerical
Recipes in Fortran. Cambridge University Press, Cambridge, 1995.
[152] S. Ptasinska, S. Denifl, S. Gohlke, P. Scheier, E. Illenberger, and T. D. Mark.
Ang. Chem. Int. Ed., 45:1893, 2006.
[153] S. Ptasinska, S. Denifl, P. Scheier, and T. D. Mark. J. Chem. Phys., 120:8505,
2004.
[154] C. Puzzarini, R. Tarroni, P. Palmieri, S. Carter, and L. Dore. Mol. Phys., 87:879,
1996.
[155] R. Renner. Z. Phys., 92:172, 1934.
[156] T. N. Rescigno, D. A. Byrum, W. A. Isaacs, and C. W. McCurdy. Phys. Rev. A,
60:2186, 1999.
[157] J. D. Robinson and R. J. Foltynowiczand E. R. Grant. J. Chem. Phys., 112:1701,
2000.

165
[158] B. R. Rowe, J. C. Gomet, A. Canosa, C. Rebrion, and J. B. A. Mitchell. J. Chem.
Phys., 96:1105, 1992.
[159] L. Sanche. Mass. Spectrom. Rev., 21:349, 2002.
[160] R. Santra and L. S. Cederbaum. Phys. Rep., 368:1, 2002.
[161] R. Santra and L. S. Cederbaum. Phys. Rep., 368:1, 2002.
[162] R. Santra, J. M. Shainline, and C. H. Greene. Phys. Rev. A, 71:032703, 2005.
[163] B. K. Sarpal and J. Tennyson. J. Phys. B, 25:L49, 1992.
[164] A.M. Scheer, K. Aflatooni, G. A. Gallup, and P. D. Burrow. Phys. Rev. Lett.,
92:068102, 2004.
[165] L. I. Schiff. Quantum mechanics. McGraw-Hill, 1968.
[166] B. I. Schneider and T. N. Rescigno. Phys. Rev. A, 37:3749, 1988.
[167] B. I. Schneider, T. N. Rescigno, B. H. Lengsfield, and C. W. McCurdy. Phys.
Rev. Lett., 66:2728, 1991.
[168] Scientific
ComputingGroup,
UniversityofBasel.
http://www.computational.unibas.ch/cs/scicomp/software/pardiso/.

Pardiso.

[169] M. J. Seaton. Rep. Prog. Phys., 46:167, 1983.
[170] M. J. Seaton. Comp. Phys. Comm., 146:225, 2002.
[171] T. C. Shen, C. Wang, G. C. Abeln, J. R. Tucker, J. W. Lyding, P. Avouris, andR. E. Walkup. Science, 268:1590, 1995.
[172] J. Shertzer and J. Botero. Phys. Rev. A, 49:3673, 1994.
[173] J. C. Slater. Quantum Theory of Molecules and Solids, volume IV. McGraw-Hill,
New York, 1974.
[174] F. T. Smith. Phys. Rev., 118:349, 1960.
[175] E. Snyder. unpublished, 2004.
[176] T. Sommerfeld. J. Phys. Chem. A, 108:9150, 2004.
[177] T. Sommerfeld and T. Posset. J. Chem. Phys,, 119:7714, 2003.
[178] O. Sueoka, C. Makochekanwa, H. Tanino, and M. Kimura.
72:042705, 2005.

Phys. Rev. A,[179] O. Sueoka and S. Mori. J. Phys. B, 19:4035, 1986.
[180] J. Q. Sun and C. D. Lin. Phys. Rev. A, 46:5489, 1992.
[181] A. Szabo and N. S. Ostlund. Modern Quantum Chemistry. Dover, Mineola, NY,
1996.

166
[182] C. Szmytkowski, A. Zecca, G. Karwasz, S. Oss, K. Maciag, B. Marinkovic, R. S.
Brusa, and R. Grisenti. J. Phys. B, 20:5817, 1987.
[183] H. Takagi. In S. L. Guberman, editor, Dissociative Recombination of MolecularIons With electrons, page 177. Kluwer Academic / Plenum Publishers, New York,2003.
[184] D. Talbi, F. Pauzat, and Y. Ellinger. Chem. Phys., 126:291, 1988.
[185] K. Tanaka and E. R. Davidson. J. Chem. Phys., 70:2904, 1978.
[186] M. Tashiro and S. Kato. J. Chem. Phys., 117:2053, 2002.
[187] J. Tennyson, C. J. Noble, and S. Salvini. J. Phys. B, 17:905, 1984.
[188] J. Tennyson and B. T. Sutcliffe. J. Chem. Phys., 77:4061, 1982.
[189] O. I. Tolstikhin, V. N. Ostrovsky, and H. Nakamura. Phys. Rev. Lett., 79:2026,
1997.
[190] O. I. Tolstikhin, V. N. Ostrovsky, and H. Nakamura. Phys. Rev. A, 58:2077, 1998.
[191] S. Tonzani. Comp. Phys. Comm., submitted, 2006.
[192] S. Tonzani and C. H. Greene. J. Chem. Phys., 122:014111, 2005.
[193] S. Tonzani and C. H. Greene. unpublished, 2005.
[194] S. Tonzani and C. H. Greene. J. Chem. Phys., 124:054312, 2006.
[195] S. Tonzani and C. H. Greene. J. Chem. Phys., submitted.
[196] C. S. Trevisan, A. E. Orel, and T. N. Rescigno. J. Phys. B, 39:L255, 2006.
[197] V. Vorsa, S. Nandi, P. J. Campagnola, M. Larsson, and W. C. Lineberger. J.
Chem. Phys,, 106:1402, 1996.
[198] A. Wasserman, N. T. Maitra, and K. Burke. J. Chem. Phys., 122:144103, 2005.
[199] S. Watanabe and C. H. Greene. Phys. Rev. A, 22:158, 1980.
[200] D. K. Watson, R. R. Lucchese, V. McKoy, and T. N. Rescigno. Phys. Rev. A,
21:738, 1980.
[201] C. A. Weatherford, M. Dong, and B. C. Saha. Int. J. Quant. Chem., 65:591, 1997.
[202] C. A. Weatherford and A. Temkin. Phys. Rev. A, 49:2580, 1994.
[203] H. J. Werner, P. J. Knowles, R. Lindh, M. Schutz, et al. MOLPRO, a package of
ab initio programs. Birmingham, UK, 2003. http://www.molpro.net.
[204] R. Wester, U. Hechtfischer, L. Knoll, M. Lange, J. Levin, M. Scheffel, D. Schwalm,
A. Wolf, A. Baer, Z. Vager, D. Zajfman, M. Mladenovic, and S. Schmatz. J. Chem.
Phys., 116:7000, 2002.

167[205] C. Winstead, P.G. Hipes, M. A. P. Lima, and V. McKoy. J. Chem. Phys., 94:5455,
1991.
[206] C. Winstead and V. McKoy. Phys. Rev. A, 73:012711, 2006.
[207] L. Yang, H. Agren, V. Carravetta, O. Vahtras, L. Karlsson, B. Wannberg, D. M. P.
Holland, and M. A. McDonald. J. Electr. Spectr., 94:163, 1998.
[208] D. R. Yarkony. Rev. Mod. Phys., 68:985, 1996.
[209] A. Zecca, C. Perazzolli, and M. J. Brunger. J. Phys. B, 38:2079, 2005.
[210] Y. Zheng, P. Cloutier, D. J. Hunting, J. R. Wagner, and L. Sanche. J. Chem.
Phys., 124:064710, 2006.
[211] E. J. Zuckerman, E. E. Mayer, R. J. Foltynowicz, J. D. Robinson, S. H. Jen,
M. C. Konopk, T. Sanford, H. G. Hedderich, I-C. Chen, and E. R. Grant. J.
Chem. Phys., 113:5372, 2000.

Appendix A

Finite element matrices

Starting from Eq. 2.8 we define the matrices Γ and Λ in our finite element
basis by transforming first to spherical coordinates (the box is spherical and the grid isalso defined in spherical coordinates), and then to rescaled coordinates, which are thevariables of the local polynomials. In the rescaled variables each sector is transformed
to a cube, in which the range of each variable is from 0 to 1. The basis set is composedof cubic Hermite polynomials in each dimension. In the rescaled coordinates systemthey are:ψ10 (ξ) = 1 − 3ξ 2 + 2ξ 3

(A.1)

ψ20 (ξ) = ξ(ξ − 1)2(A.2)

ψ11 (ξ) = ξ 2 (3 − 2ξ)(A.3)

ψ21 (ξ) = ξ 2 (ξ − 1).

(A.4)

The wavefunction inside each sector can therefore be expanded asu(ξ1 , ξ2 , ξ3 ) =

X

(lmn)

ψil (ξ1 )ψjm (ξ2 )ψkn (ξ3 )Cnode(A.5)

i,j,k,l,m,nwhere i, j, k can be 1 if the polynomial has nonzero value at some node or 2 if it hasnonzero derivative, whereas l, m, n can assume values of 0 if that node is the first for

169

Finite element basis
1

1
2
3
4

0.8
0.6
0.4
0.2
0
−0.2
0

0.2

0.4

0.6

0.8

1

ξ

6
5

ξ3

8
7

4

2

ξ1
1

ξ2

3

Figure A.1: Top: finite element basis functions from Eqs. A.1-A.4, the four polynomials
are plotted in the same order as in the text. Bottom: nodal structure for one cubic sector
in the rescaled variables.

170
the variable of the polynomial in the sector or 1 if it is the last; ξi are the local rescaled
variables . The nodal structure of each element is represented in Fig. A.1, together(lmn)with the basis functions. The coefficients Cnode are the values of the wavefunction andits derivatives at the node, and they are to be determined solving Eq. 2.8. If we define
ak,p = xk,p,i+1 − xk,p,i(A.6)

xk,p = ak,p ξk + xk,p,i(A.7)

where k indexes the spherical coordinates and p the sectors in which they are defined,
xk,p,i and xk,p,i+1 are the initial and final points for the variable xk in sector p, theexpressions for the matrices become:
Γij =

Z "X
3

#

F (xk ) ∂ui ∂uj+ 2ui (U − E)uj ar aθ aφ r2 sin2 θdξ1 dξ2 dξ3
a
a
∂ξ
∂ξ
k
k
k
k
k=1
Λmn =

Z

∗
Ylm
(θ, φ)Yl0 m0 (θ, φ) sin θ dθ dφ = δll0 δmm0

(A.8)

(A.9)

where F (xk ) is a spherical coordinates scale factor, and it is 1 if xk = r and 1/r 2 and1/(r 2 sin2 θ) for θ and φ respectively. Imposing function and derivative continuity foru(ξ1 , ξ2 , ξ3 ) amounts to require that the indices of the same node across neighboring
sectors be the same. This in turn leads to having to perform a sum of the integrals in
Eq. A.8 when evaluating the matrix element at a node, across all sectors that sharethat node. The continuity conditions we apply are very similar to the ones in Eqs. 1.16,
1.18 of Ref. [21], but adapted to a three dimensional case.

Appendix B

Related publications• “Electron-molecule scattering calculations in a 3D finite element R-matrix approach” S. Tonzani and C. H. Greene, J. Chem. Phys., 122, 014111 (2005)• “Dissociative recombination of HCO+ ” . Larsson, S. Tonzani, R. Santra and C.
H. Greene, J. Phys. Conference Series, 4, 148 (2005)• “Low energy electron scattering from DNA and RNA bases: shape resonancesand radiation damage” S. Tonzani and C. H. Greene, J. Chem. Phys., 124,
054312 (2006) and Virtual Journal of Biological Research 11 (4) (2006)• “Renner-Teller effects in HCO+ dissociative recombination” I. Mikhaylov, V.
Kokoouline, A. Larson, S. Tonzani and C. H. Greene, Phys. Rev. A, 74, 032707(2006).
• “FERM3D: A finite element R-matrix electron molecule scattering code” S.
Tonzani, Computer Physics Communications (in press).
• “Radiation damage to DNA: electron scattering from the backbone subunits”
S. Tonzani and C. H. Greene, J. Chem. Phys., 125, 094504 (2006)