Document [original]

Theoretical Modeling and Simulation of

Electron-Phonon Scattering Processes in

Molecular Electronic Devices

Alessio Gagliardi

Theoretical Modeling and Simulation of

Electron-Phonon Scattering Processes in

Molecular Electronic Devices

Dissertation zur Erlangung des akademischen Grades

Doktor der Naturwissenschaften (Dr. rer. nat.)

vorgelegt dem

Department Physik der Fakult¨

at f¨

ur Naturwissenschaften

Universit¨

at Paderborn

Dipl. Ing. Alessio Gagliardi

Paderborn, 2007

Dem Department Physik der Fakult¨

at f¨

ur Naturwissenschaften als Dissertation vorgelegt.

Eingereicht am: 4.5.2007

Tag der m¨

undlichen Pr¨

ufung: 29.6.2007

Promotionskommission

Vorsitzender Prof. Dr. Artur Zrenner

Erstgutachter Prof. Dr. Thomas Frauenheim

Zweitgutachter Prof. Dr. Wolf Gero Schmidt

Drittgutachter Prof. Dr. Aldo Di Carlo

Beisitzer Prof. Dr. Torsten Meier

Chapter 6reproduced with permission from A. Gagliardi, A. Pecchia, A. Di Carlo, S.

American Chemical Society.

Chapter 7reproduced in part with permission from, A. Gagliardi, G. C. Solomon, A.

Pecchia, J. R. Reimers, A. Di Carlo, Th. Frauenheim, N. S. Hush, Journal of Chemical

Chapter 8reproduced in part with permission from A. Gagliardi, G. C. Solomon, A.

Pecchia, J. R. Reimers, A. Di Carlo, Th. Frauenheim, N. S. Hush, Nano Letters, Vol. 6,

Chapter 9reproduced in part with permission from A. Gagliardi, G. C. Solomon, A.

Pecchia, J. R. Reimers, A. Di Carlo, Th. Frauenheim, N. S. Hush, Physical Review B,

Archiv

Elektronische Dissertationen und Habilitationen der Universit¨

at Paderborn

http://www.ub.upb.de/ebibliothek/hochschulschriften/ediss

Alessio Gagliardi, Theoretical Modeling and Simulation of Electron-Phonon Scattering Processes

in Molecular Electronic Devices. PhD Thesis (in English), Department of Physics, Faculty of

Science, University of Paderborn, Germany (2007).

Abstract

The entire field of molecular electronics relies on a very simple concept: a single organic molecule,

or a thin layer of organic molecules, are used as a bridge between two metallic contacts. The

bias applied induces a current from one contact to the other. The basic idea is that the current is

dominated by the chemical properties of the molecule between the contacts and that this electric

characteristic can be easily tuned by an external perturbation in order to get a switch. However, de-

spite the intense investigation many issues are still waiting to be fully addressed. One of them is the

dissipation induced in the molecule by the coupling between the vibrational modes and the charge

carriers. Due to the size of the device this interaction must be described quantum mechanically via

an electron-phonon scattering process. The quantitative evaluation of dissipation in molecular elec-

tronics is fundamental for the entire field, especially for future introduction of such kind of devices

in the electronic industry.

The author of this work and Dr. Alessandro Pecchia of the University of Rome “Tor Vergata”

have implemented a code based on non-equilibrium Green’s functions formalism to investigate

electron-phonon scattering in molecular devices. The first application presented is the evaluation

of the dissipation in octanedithiols sandwiched between two gold contacts showing which are the

most relevant modes in scattering electrons. The electron-phonon scattering is much more than the

study of the “quantum resistance”. The inelastically scattered electrons in fact can be directly used

as probes to investigate the geometry and the electronic structure of the device. This spectroscopy

technique is called inelastic electron tunneling spectroscopy (IETS). Despite the IETS is an old

spectroscopy technique, its application in the molecular electronic field is recent and many features

need to be fully understood.

The electron-phonon code has been used to investigate the properties of IETS, first simulating

directly IETS in octane-thiol systems, later on addressing the problem of selection rules in this

spectroscopy technique. It is well known in fact that the selectivity of IETS is rather different from

the one of other more conventional techniques, like infrared or Raman. However the lack of a

precise knowledge of selection rules in IETS renders a clear assignment of IET spectra difficult.

The definition of selection rules requires two steps. First, the definition of the point group of

the system. It is obvious that, also if the molecule in the gas phase shows high symmetry, this

will be reduced by the presence of the electrodes and the distortions of the structure. Despite that,

the molecule remains the dominant feature of IETS measurements. We can univocally identify the

point group of the device just taking the point group of the unperturbed molecule and reducing it to

a subgroup due to the effects of the electrodes. The second issue is to define a basis set to which the

point group is applied to. It is well known that canonical molecular orbitals are not the best choice

due to the strong interference effects between them during transport. However, the current can

be decomposed in a sum of independent “channels” obtained diagonalizing directly the operator

involved in the transmission and not the Hamiltonian. The choice of this set of channels is not

unique, but all of them permit to avoid interference.

These tools have been applied to investigate IETS in a benzenedithiols molecule between model

gold contacts. The interpretation of the current in terms of independent channels provides not only

information about selection rules, but also in which region electrons are scattered, leading to a deep

comprehension of where the dissipation occurs and which chemical substitutions in the device can

improve its performances.

Alessio Gagliardi, Theoretische Modellierung und Simulation von Elektron–Phonon Streuprozessen

in Molekularelektronischen Bauteilen. Dissertation (in englischer Sprache), Department Physik,

Fakult¨

at f¨

ur Naturwissenschaften, Universit¨

at Paderborn (2007).

Kurzfassung

Das Gebiet der Molekularelektronik beruht auf einem einfachen Konzept: Ein einzelnes oder eine d¨

unne

Schicht organischer Molek¨

ule verbindet zwei metallische Kontakte. Das Anlegen einer Spannung induziert

einen Strom von einem Kontakt zum anderen. Die Grundidee ist, dass der Strom wesentlich von den chemis-

chen Eigenschaften des Molek¨

uls zwischen den Kontakten abh¨

angt und dass diese elektronische Charakter-

istik leicht mit Hilfe einer ¨

außeren St¨

orung gestimmt werden kann, um einen Schalter zu realisieren. Trotz

intensiver Untersuchungen warten jedoch viele Aspekte noch immer auf ihre eingehende Behandlung. Einer

davon ist die Dissipation, welche im Molek¨

ul durch die Kopplung zwischen den vibronischen Moden und

den Ladungstr¨

agern verursacht wird. Durch die Gr¨

oße des Bauteils muss diese Wechselwirkung quanten-

mechanisch durch einen Elektron–Photon Streuprozess beschrieben werden. Die quantitative Auswertung

der Dissipation ist fundamental f¨

ur das gesamte Feld der Molekularelektronik, insbesondere in Hinblick auf

die industrielle Anwendungen.

Der Author dieser Arbeit und Dr. Alessandro Pecchia von der Universit¨

at von Rom “Tor Vergata” haben

einen Code implementiert, der auf dem Nichtgleichgewichts–Green–Funktionen–Formalismus basiert, um

Elektron–Phonon Streuung in molekularen Bauteilen zu untersuchen. Die erste Anwendung ist die Berech-

nung der Dissipation in Oktanthiolat eingebettet zwischen zwei Goldkontakten, die zeigt, welches die rele-

vantesten Moden bei der Elektronenstreuung sind, eine fundamentale Frage f¨

ur die Stabilit¨

at des Bauteils.

Die inelastische Streuung umfasst jedoch weitaus mehr als eine Untersuchung des “Quantenwiderstands” des

Bauteils. Die inelastisch gestreuten Elektronen k¨

onnen direkt als Sensor f¨

ur die Geometrie und Elektronen-

struktur des Bauteils verwendet werden. Diese spektroskopische Technik wird als Inelastische Elektronen–

Tunnel–Spektroskopie (IETS) bezeichnet. Obwohl bereits recht alt, ist die IETS auf dem Gebiet der Moleku-

larelektronik neu, und viele Eigenschaften m¨

ussen erst noch verstanden werden.

Der Elektron–Photon Code wurde verwendet, um die Eigenschaften der IETS zu untersuchen, zun¨

achst

durch direkte Simulation des IETS in Oktandithiol–Systemen, sp¨

ater, um das Problem der Auswahlregeln

anzugehen. Tats¨

achlich ist es bekannt, dass die Selektivit¨

at der IETS sehr verschieden von der anderer,

konventionellerer Techniken wie Infrarot oder Raman ist. Das Fehlen einer genauen Kenntnis ¨

uber die

Auswahlregeln in der IETS gestaltet eine klare Zuordnung der Spektren schwierig. Die Definition von

Auswahlregeln erfordert zwei Schritte. Zuerst muss die Punktgruppe des Systems festgelegt werden. Auch

wenn das Molek¨

ul in der Gasphase hochsymmetrisch ist, ist es offensichtlich, dass die Symmetrie durch

Anwesenheit der Elektroden und strukturelle Verzerrungen verringert wird. Dennoch bleibt das Molek¨

ul das

wesentliche Merkmal der IETS Messungen. Unter Ber¨

ucksichtigung des wesentlichen Effekts der Elektro-

den, k¨

onnen wir die Punktgruppe des Systems eindeutig definieren als diejenige des ungest¨

orten Molek¨

uls,

abz¨

uglich aller End–End–Symmetrieoperationen. Das zweite Problem besteht darin, einen geeigneten Basis-

satz zu finden. Es ist wohl bekannt, dass kanonische Molek¨

ulorbitale nicht die beste Wahl sind, wegen der

starken Interferenzeffekte zwischen ihnen w¨

ahrend des Transports. Der Strom kann jedoch in eine Summe

voneinander unabh¨

angiger “Kan¨

ale” zerlegt werden, die man durch Diagonalisierung der Matrix des Trans-

missionsoperators (nicht des Hamiltonians) erh¨

alt. Die Wahl dieser Menge von Kan¨

alen ist nicht eindeutig,

jede Wahl erm¨

oglicht es jedoch, Interferenz zu vermeiden. Die Interpretation des Stromes mit Hilfe un-

abh¨

angiger Kan¨

ale liefert auch Informationen dar¨

uber, wo die Elektronen gestreut werden, was zu einem

tiefer gehenden Verst¨

andnis f¨

uhrt, wo die Dissipation stattfindet und welche chemischen Substitutionen am

Bauteil die Leistungsf¨

ahigkeit des Bauteils erh¨

ohen k¨

onnen.

Contents

Introduction 1

1 Molecular Electronics: a Brief Overview 5

1.1 Molecular electronics experimental techniques . . . . . . . . . . . . . . . . 6

1.1.1 Break junction experiments . . . . . . . . . . . . . . . . . . . . . 6

1.1.2 Molecular monolayer devices . . . . . . . . . . . . . . . . . . . . 8

1.1.3 Nanopores .............................. 10

1.1.4 Electromigration experiments . . . . . . . . . . . . . . . . . . . . 11

1.1.5 Scanning tunneling microscope experiments . . . . . . . . . . . . . 12

1.2 Inelastic electron tunneling spectroscopy . . . . . . . . . . . . . . . . . . . 13

2 Modeling Molecular Electronic Devices 17

2.1 The general problem of molecular conduction . . . . . . . . . . . . . . . . 18

2.2 The Hamiltonian of the system . . . . . . . . . . . . . . . . . . . . . . . . 20

2.3 From the Hamiltonian to the current: Meir-Wingreen equation . . . . . . . 22

2.4 Lifetimesofinterest.............................. 26

3 Non-Equilibrium Green’s Functions 29

3.1 Threerepresentations............................. 30

3.1.1 Schr¨

odinger representation . . . . . . . . . . . . . . . . . . . . . . 30

3.1.2 Heisenberg representation . . . . . . . . . . . . . . . . . . . . . . 31

3.1.3 Interaction representation . . . . . . . . . . . . . . . . . . . . . . . 31

3.2 S-Matrix.................................... 32

3.3 Equilibrium Green’s functions . . . . . . . . . . . . . . . . . . . . . . . . 33

3.4 Wick’stheorem................................ 34

3.5 Feynman’sdiagrams ............................. 36

3.6 Dyson’sequation ............................... 39

3.7 Time-loop S-matrix:NEGF.......................... 41

3.7.1 Dyson’s equations for NEGF . . . . . . . . . . . . . . . . . . . . . 44

4 Density Functional Based Tight-Binding 47

4.1 DFT and Kohn-Sham formulation . . . . . . . . . . . . . . . . . . . . . . 48

CONTENTS

4.1.1 The Kohn-Sham equations . . . . . . . . . . . . . . . . . . . . . . 49

4.2 DFTB: method and approximations . . . . . . . . . . . . . . . . . . . . . 51

4.2.1 Pseudo-atomic starting density . . . . . . . . . . . . . . . . . . . . 52

4.2.2 Tight-binding integrals and the two-centre approximation . . . . . . 53

4.2.3 Repulsivepotential.......................... 54

4.2.4 Second-order correction . . . . . . . . . . . . . . . . . . . . . . . 55

4.2.5 DFTB secular equation . . . . . . . . . . . . . . . . . . . . . . . . 57

4.3 Disadvantages of DFT in transport simulations . . . . . . . . . . . . . . . 58

5 The Electron-Phonon Code 61

5.1 Approximations in the electron-phonon code . . . . . . . . . . . . . . . . . 62

5.2 The scheme of the device and the open boundary conditions . . . . . . . . 64

5.3 The electron-phonon self-energies . . . . . . . . . . . . . . . . . . . . . . 66

5.4 Computation of the electron-phonon couplings . . . . . . . . . . . . . . . 68

5.5 The flowchart of the code . . . . . . . . . . . . . . . . . . . . . . . . . . . 69

6 Power Dissipation in Molecular Electronic Devices 73

6.1 Dissipation in alkanethiols . . . . . . . . . . . . . . . . . . . . . . . . . . 74

6.1.1 Geometry and vibrational modes . . . . . . . . . . . . . . . . . . . 74

6.1.2 Power dissipation in the molecule . . . . . . . . . . . . . . . . . . 75

7 Simulation of Inelastic Electron Tunneling Spectroscopy in Alkanethiols 81

7.1 IETSapproximation.............................. 81

7.2 The choice of binding site . . . . . . . . . . . . . . . . . . . . . . . . . . . 82

7.2.1 Atomic partitioning of the system . . . . . . . . . . . . . . . . . . 83

7.2.2 Correction of the DFTB vibrational frequencies . . . . . . . . . . . 84

7.3 Discussion of the IETS simulations . . . . . . . . . . . . . . . . . . . . . . 84

7.3.1 Nature of the orbitals controlling conduction and IETS . . . . . . . 88

7.3.2 Nature of the vibrations that produce IETS . . . . . . . . . . . . . 89

7.4 Conclusions.................................. 91

8 The Role of Symmetry and Channels in Conduction 93

8.1 The definition of the symmetry of conduction . . . . . . . . . . . . . . . . 94

8.1.1 Application to molecular system . . . . . . . . . . . . . . . . . . . 95

8.2 The definition of channels in elastic transport . . . . . . . . . . . . . . . . 99

8.2.1 B¨

uttikerchannels...........................100

8.2.2 Γ-channels..............................103

8.3 Conclusions..................................106

9 Propensity Rules in IETS 109

9.1 Theoreticalformalism.............................109

9.2 From Γ-channels to A-channels . . . . . . . . . . . . . . . . . . . . . . . 111

9.3 Propensity rules in IETS . . . . . . . . . . . . . . . . . . . . . . . . . . . 118

9.4 Conclusions..................................121

CONTENTS

vii

Conclusion 123

A Atomic Units 127

B Surface Green’s Function: Decimation Technique 129

Bibliography 133

OwnPublications..................................145

Personal Contributions 147

Acknowledgments 148

Introduction

[...] I don’t know how to miniaturize a computer on a small scale in a practice

way, but I do now that computing machines are very large; they fill rooms.

Why can’t we make them very small, make them of little wires, little elements

and by little, I mean little. [...]

[...] But there is plenty of room to make them smaller. There is nothing that I

can see in the physical laws that says the computer elements cannot be made

enormously smaller than they are now [...]

[...] It is not an attempt to violate any laws; it is something, in principle, that

can be done; but, in practice, it has not been done because we are too big [...]

Richard P. Feynman There is a plenty of room at the bottom, 1959.

The idea of engineering at molecular scale has always been a dream of scientists. In his

speech quoted at the beginning of this introduction, Feynman, more than fifty years ago,

remarked that “there is a plenty of room at the bottom” to miniaturize our devices up to a

very microscopic scale.

In the last 30 years, the silicon metal oxide semiconductor (MOS) field-effect transistor

(FET) has been the foundation of integrated circuits [1]. In the 70s, integrated circuits con-

tained thousands of MOSFETs, each having dimensions of tens of micrometres. Today’s

chips contain almost one billion MOSFETs, each having physical dimensions of tens of

nanometres. The exponential increase of device density (Moore’s Law [2]) and scaling of

device dimension into the nanotechnology regime (see Fig. (1)) have resulted in the vast

improvements observed in numerous electronic devices, from PCs to control systems found

in airplanes.

Although substantial innovation was required to realize this rate of dimensional scaling,

until recently very little has changed in the materials and design of the basic MOSFET. In

the future, a variety of new materials and device structures will be required to continue

MOSFET scaling [3,4,5,6]. Furthermore, as silicon MOSFET technology approaches

its limits [7], entirely new device structures and computational paradigms will be required

to replace and improve traditional MOSFETs. These possible emerging technologies span

from transistors made of silicon nanowires to devices made of nanoscale molecules.

Introduction

Figure 1: The Moore’s Law in a diagram: the exponential growth of the number of transistors on

a single chip in the last thirty years.

Independently the future of electronics will be based on nanoscale planar silicon MOS-

FETs or an emerging technology which totally replaces the MOSFET, the electronic prop-

erties of all of these nanodevices are extremely susceptible to small perturbations in proper-

ties such as dimension, structure, roughness and defects, which means there is a significant

need for precise metrology. Furthermore, the insertion of a variety of new materials into

the conventional MOSFET and radically new materials and devices for potential emerg-

ing replacement technologies (for example, molecular and spin electronics) present further

challenges for metrology. Although the devices and materials being considered for future

MOS technologies and beyond are broad, the overarching challenges to metrology remain.

The ability to measure the physical, chemical and electronic properties that control final

device electrical characteristics will be crucial to the research and development of future

electronic devices.

In this rush to get smaller and faster devices, molecular electronic represents the ulti-

mate frontier. The potential for molecular electronics originates from its nanoscale dimen-

sion, the possibility of synthesizing very specific electronic properties, and the promise

of fabrication through self-assembly. Although numerous molecular devices have been

demonstrated, there is still significant controversy concerning the repeatability of measure-

ments and the understanding of charge transport [8,9,10]. Many problems are related

to the difficulties to investigate the geometry and the electronic structure of the molecule

connected to the metal electrodes. The final device structure dramatically changes the

intended electronic properties of the material of interest. The interaction between the elec-

trodes and molecules is critical to understand the transport properties of the final device.

Forming ohmic but non-interacting metallic contacts to molecules is a specific problem. A

new set of experimental techniques, such as inelastic tunneling spectroscopy, which probes

Introduction

vibrational modes of the molecules through measurement of the device current-voltage be-

havior, have been developed to determine which molecular species and bonds are present

in the final device structure [10]. New effects can arise creating a completely new kind

of electronics, such as Kondo effect [11], negative differential resistance [12] or Coulomb

blockade [11]. All these new problems, techniques and effects represent a continuous chal-

lenge for researchers.

The silicon electronics has been a technological and economic engine for over thirty

years. Metrology has been an important infrastructure enabling the exponential increase

of device density and exponential decrease of device dimension. New nanoscale materials

and devices will be required to sustain this technological revolution but the extreme sen-

sitivity of the electronic properties of these devices to their nanoscale physical properties

and the insertion of radically new materials present significant challenges for their devel-

opment. The ability to measure the physical, chemical and electronic properties of these

new materials and devices will become even more critical in continuing the advance of the

MOS and finding technologies to extend or replace it.

Concluding, “there is a plenty of space at the bottom” and we are just scratching the

surface.

Chapter 1

Molecular Electronics: a Brief Overview

In this Chapter is given a brief introduction to some of the most important experimental

techniques which are gaining popularity in the molecular electronic field. The techniques

presented are not the only in existence and also for the techniques enlisted here are many

variants not mentioned. However, a very detailed discussion of all methods to create molec-

ular devices is beyond the scope of this work. The techniques described have been selected

not only for their importance, but also for the relevance they have in the present work.

The setup of molecular electronic devices, despite the wide spectrum of variation that

have been developed in the recent years, can be summarized in few words: an organic or

inorganic molecule, or a thin layer of molecules, bridging two contacts. The geometry

can be very different from one device to the other as much as the molecules used and the

materials of the electrodes. Moreover a gate can or cannot be present. This has produced

a wide range of possibilities for this kind of devices, nevertheless the main characteristic

of the device remains the molecule (or molecules) used to create the bridge. In fact the

very basic idea of molecular electronics is that the electronic properties of the bridging

molecules can be used directly to create a switch, and in electronics a switch is the first

step to define a logic circuit.

The first experiments for molecular electronics concerned mainly measurements of

conductivity, because this is the first, the easiest and, in many respects, the most impor-

tant quantity we want to know about an hypothetic molecular device. Nevertheless, in the

following years, other aspects have become more and more important. The discovery of

very peculiar effects related to the strong confinement of the electrons in the molecular

device region, important spin-spin correlation processes and many other properties related

to the quantum nature of the physics at molecular scale have started to be investigated

experimentally and theoretically. Moreover another important issue has become very pop-

ular: the dissipation in these devices. It is obvious that if a current is flowing in this tiny

molecular bridge a certain degradation of the flux of current must be expected.

However, the magnitude and the quality of this degradation, that we can address as

a sort of “quantum resistance”, are still under investigation. Even more, the dissipation

is certainly a fundamental ingredient to address the stability of electronic devices, in fact

Chapter 1. Molecular Electronics: a Brief Overview

only a complete understanding of the exchange of energy between the environment and

the molecule can give an insight into the mechanisms which underlie the conduction in

the device. Finally, the electrons inelastically scattered can be directly used as probes to

investigate the device itself as in the Inelastic Electron Tunneling Spectroscopy (IETS). The

last part of this Chapter is devoted to describe the IETS due to the relevance this technique

has in the present work.

Of all the experimental setups described in this Chapter Nanopores and Molecular

Break-Junction are certainly the most important for the present work. Finally a remark:

every technique is offering particular advantages and drawbacks, but, so far, none of them

can be viewed as a “standard model” for experimental setup and all of them are far away

from being implemented on an industrial scale to yield commercial devices. Contacting

and manipulating few molecules between bulk contacts on the nanometer scale in a reli-

able and reproducible fashion still constitutes a scientific and engineering challenge.

1.1 Molecular electronics experimental techniques

1.1.1 Break junction experiments

The use of Mechanically Controlled Break-Junctions (MCBJ) [14,15] for contacting single

molecules is a further development of single-atom contact conductance experiments [16,

17]. The gap to introduce the molecules which will form the bridge is created by elongating

a metal wire by a mechanical stress. Shortly before breaking, the contact between the

two halves reduces to one single atom. Further elongation yields an atomic chain of gold

atoms until the wire finally breaks. The design principal of MCBJ device is simple: a thin

gold film is structured, using e-beam lithography, on an elastic substrate. The gold film is

designed to have a constriction. By underetching this constriction, a freestanding bridge is

created. The sample is placed in a bending mechanism as depicted in Fig. (1.1).

The very precise control of the bending is allowed by a three point geometry. To avoid

contaminations, the bending mechanism is placed in a vacuum chamber. The pushing rod

(a piezoelectric element) is moved upward until the bridge breaks. The upward movement

is converted, by the bending of the substrate, in a variation of the distance between the two

halves of the metallic sample. The control of this movement is extremely precise, up to 0.1

A. Now, the contacts are opened and closed repeatedly until quantized conductance values

are observed, corresponding to the formation of a sharp atomic point contacts.

Then the chamber is filled with an inert gas (like argon) and a solution of the molecule

(with terminal protective groups attached to the thiol groups) to be measured is applied

on the opened contacts. The terminal group prevents the molecules from forming clusters

and only attach in contacts with the gold forming a stable covalent bond between sulfur

and the gold. The thiol group is the most popular choice because the strong covalent bond

it forms with the gold atoms. After a while the solvent will evaporate and some of the

molecules will loose their protection terminal group and the sulfur will chemically bind to

one electrode.

1.1. Molecular electronics experimental techniques

Figure 1.1: 1)Schematic of the measurement process. (A) The gold wire of the break junction

before breaking and tip formation. (B) After addition of the molecule (in this case benzene-1,4-

dithiol in a tetrahydrofuran (THF) solution), layers form on the gold wire surfaces. (C) Mechanical

breakage of the wire in solution produces two opposing gold contacts that are covered by a layer

of molecule each. (D) After the solvent is evaporated, the gold contacts are slowly moved together

until the onset of conductance is achieved. 2)A schematic of the experimental setup with (a) the

bending beam, (b) the counter support, (c) the notched gold wire, (d) the glue contacts, (e) the

piezo element, (f) the glass tube containing the solution. 3)A schematic of a benzene-1,4-dithiolate

layer of molecule between proximal gold electrodes formed in a MBJE. The thiolate is normally H

terminated after deposition; end groups denoted as X can be either H or Au, with the Au potentially

arising from a previous contacts/retraction event. These molecules remain nearly perpendicular to

the Au surface, making other molecular orientations unlikely (from [13]).

To remove molecules that are not chemically bound to the gold, the junction is rinsed

with a solvent. After evacuation, the bridge is slowly closed with bias voltage applied.

Due to the small distance between the two electrodes, a high electric field is present, in

which the molecules orient themselves toward the opposite electrode. At one point, the

first molecules will reach the second electrode and chemically bind to it. This situation can

be recognized by a plateau in the current. Now, current-voltage curves can be recorded.

Break-junction offer the best control for single-molecules experiments [13,18,19,20,21],

nevertheless the main drawback of the technique is the impossibility to introduce a gate in

the device due to the large diameter of the breaking gold wire. Every gate electric field is

in fact completely screened by the metal.

Chapter 1. Molecular Electronics: a Brief Overview

Figure 1.2: A) Langmuir-Blodgett (LB) films: a layer of molecules with an hydrophilic group is

laying on the surface of water. The solid substrate is merged in the water (a). When the layer is

compressed the molecules start to climb to the solid substrate (b). The step can be repeated to

get multi-layers systems (c). B) A self-assembled monolayer structure (SAM) formed by molecules

with a thiol group in solution with a metallic substrate. The molecules, after a while, bind to the

substrate.

1.1.2 Molecular monolayer devices

The main difference between this approach and the former is that more than one molecule

is used to make the device. The principal advantage relates to the simplicity to produce well

structured layers. The substrate of the layer act as one of the electrode while the second is

then evaporated on top, yielding a film device, or it is a Scanning Tunneling Microscope

(STM) tip. The main difficulty in film devices is the deposition of the second electrode

that can easily damage the monolayer. Although working devices are achieved in most

cases, the actual microscopic structure of the measured device (structural integrity of the

monolayer, defects) is still not clear [22]. Several devices based on sandwiching molecular

films between two electrodes have been reported [23,24,25]. Different techniques can be

applied to produce a molecular monolayer: Langmuir-Blodgett films and Self-Assembled

Monolayers. Due to the simplicity of their formation these devices are very promising

1.1. Molecular electronics experimental techniques

possibilities for the future of molecular electronics, especially for hybrid devices which

combine with more standard electronics.

Langmuir-Blodgett films (LB): Langmuir-Blodgett films consist of mono-molecular

layers stacked sequentially onto a solid substrate. A solid substrate, covered by a hy-

drophilic material, is lowered into the water, breaking through the LB film which stays

on the surface of the liquid. The molecules of the LB must contain a hydrophilic group,

e.g. acid or alcohol group, and a hydrophobic group, usually an aliphatic chain, to allow

film formation on the surface of the water. In water, all the molecules will then align in

the same direction with the hydrophilic end at the water side. A continuous monolayer,

which can later be transferred on the substrate, is formed by compression of the film using

a moveable barrier. The substrate is slowly removed from the bath while at the same time

the barrier is shifted to keep the monolayer intact. This way, a well-ordered monolayer of

the organic molecule is formed on the surface of the substrate (see Fig. (1.2, a)). The big

disadvantage of this method is that the molecules are not chemically bound to the substrate,

as a hydrophilic group is required on the substrate side of the molecule.

Self-Assembled Monolayers (SAMs): Molecular self-assembly is a very well known

strategy for nano-fabrication that involves designing molecules and supramolecular enti-

ties so that shape-complementarity causes them to aggregate into desired structures. Self-

assembly has a number of advantages as a strategy: first, it carries out many of the most dif-

ficult steps in nano-fabrication, those involving atomic-level modification of structure, us-

ing the very highly developed techniques of synthetic chemistry. Second, it draws from the

enormous wealth of examples in biology for inspiration: self-assembly is one of the most

important strategies used in biology for the development of complex, functional structures.

Third, it can incorporate biological structures directly as components in the final systems.

Fourth, because it requires that the target structures be the thermodynamically most stable

ones open to the system, it tends to produce structures that are relatively defect-free and

self-healing.

Although there are countless examples of self-assembly all around us the basic rules

that govern these assemblies are not understood in useful detail, and self-assembling pro-

cesses cannot, in general, be designed and carried out “to order”. Many of the ideas that are

crucial to the development of this area: “molecular shape”, the interplay between enthalpy

and entropy, the nature of non-covalent forces that connect the particles in self-assembled

molecular aggregates, are simply not yet under the control of investigators. In the typi-

cal SAM setup for molecular electronics the monolayer forms from molecules in solution

on the substrate itself, due to the properties of the molecules [23]. These systems – self-

assembled monolayers (SAMs) – are reasonably well understood, and increasingly useful

technologically. The crucial dimension in SAMs is the thickness perpendicular to the plane

of the monolayer: this dimension, and the composition along this axis, can be controlled

very simply at the scale of 0.1 nm by controlling the structures of the molecules making

up the monolayer. After a layer has formed, the substrate is removed from the solution

(see Fig. (1.2, b)). The second electrode can then be evaporated on top of the molecular

layer. Many measurements using SAM devices have been reported [24,26]. For SAMs

Chapter 1. Molecular Electronics: a Brief Overview

the technique is simpler than LB and chemical bonds to the substrate can and have to be

achieved.

1.1.3 Nanopores

Figure 1.3: Fabrication of the heterostructure. (a)Cross section of a silicon wafer showing the

bowl-shaped pore etched in suspended SiN membrane with a diameter about 300 ˚

A. (b)Au-Ti

top electrode/self-assembled monolayer/Au bottom electrode sandwich structure in the nanopore

(see [27]).

The idea of nanopore links directly to the one of SAM. The idea is to have a small hole

in which deposit a smaller number of molecules than in SAMs and LBs techniques. The

self-assembling properties of the molecules and the confinement of the hole lead the system

to form a well structured layer. The first measurement of conduction with nanopores was

presented by Zhou and Reed [27] at the same time of their publication of break-junction

experiment. The hole is etched into a freestanding silicon nitride layer by reactive ion

etching. After the etching, gold is evaporated onto the hole from underneath. Then a SAM

is allowed to form on the evaporated gold. Finally, titanium and gold are evaporated from

above to yield the molecular device (Fig. (1.3)). The main advantage of this structure is its

stability and the reduction of size if compared to more traditional monolayers devices. The

big disadvantage is the deposition of the second electrode, the effects of that are still not

clear as for SAMs, moreover it is not clear which is the interaction between the titanium

and the molecules.

1.1. Molecular electronics experimental techniques

Figure 1.4: Kondo effect in a single molecule device. Experiment by J. Park et al. [11] using

electromigration contacts. Left: molecules investigated, consisting of a Co ion trapped in a cage of

benzene rings. Right: observed Kondo peak in the shorter molecule. Displayed is the differential

conductance over bias voltage in V. Different traces correspond to different temperatures. The inset

shows the decrease of the peak height with temperature.

1.1.4 Electromigration experiments

Devices fabricated by electromigration offer the possibility to study the influence of a gate

electrode. Usually, a 10 to 15 nm wide wire is fabricated on a SiO2insulating layer using

e-beam lithography [28]. A doped Si substrate underneath the SiO2can be used as a gate

electrode. After adding a solution of the molecules to be investigated onto the wire, the

single-molecule contact is then created by breaking the wire through electromigration: at

low temperature the applied bias voltage is increased to large values until the wire breaks.

This produces a gap of about 1 to 2 nanometers, across which in about 10% of the samples

a molecule is found, so that the current-voltage characteristic can be recorded. Using the

electromigration technique, J. Park et al. [11] managed to observe the Kondo effect (see

Fig. (1.4)) and the Coulomb blockade effect in single molecules consisting of a cobalt ion

trapped in a polypyridyl cage. For the longer, less conducting molecule, a conductance gap

for small bias voltage was observed. In dependence of the gate voltage different threshold

bias values were necessary until the current increased, corresponding to Coulomb blockade

behavior. 2Dplots of conductance in dependence of applied bias and gate voltage yielded

Coulomb diamonds. For the shorter molecule, where the cobalt-ion acted as an “impurity

spin” a Kondo resonance was reported: at zero bias, a peak in the difference conductance

has been reported, that has a logarithmic temperature dependence with temperature be-

tween 3 K and 20 K.

As a second test, a magnetic field was applied which resulted in a split of the peak,

as would be expected from Kondo physics. A second group, H. Park et al. [29], also

observed a Kondo resonance in single molecules using electromigration techniques to con-

tact the molecules. There are many other experiments that make use of the electromigration

technique to contact single molecules [29,30,31].

Chapter 1. Molecular Electronics: a Brief Overview

1.1.5 Scanning tunneling microscope experiments

Figure 1.5: (A)Conductance of a gold contact formed between a gold STM tip and a gold sub-

strate decreases in quantum steps near multiples of G0= (2e2/h)as the tip is pulled away from

the substrate. (B)A corresponding conductance histogram constructed from 1000 conductance

curves as shown in (A) shows well-defined peaks near 1 G0, 2 G0, and 3 G0due to conductance

quantization. (C)When the contact shown in (A) is completely broken, corresponding to the col-

lapse of the last quantum step, a new series of conductance steps appears if molecules such as 4,4’

bipyridine are present in the solution. These steps are due to the formation of the stable molecular

junction between the tip and the substrate electrodes. (D)A conductance histogram obtained from

100 measurements as shown in (C) shows peaks near 1 ×, 2 ×, and 3 ×0.01 G0that are ascribed

to one, two, and three molecules, respectively. (Eand F) In the absence of molecules, no such steps

or peaks are observed within the same conductance range, see [32].

Using monolayer film devices, it is also possible to conduct measurements on single

molecules. A Scanning Tunneling Microscope (STM) tip is used as second electrode.

Often films consisting of a matrix of spacer molecules plus the molecules to be measured

are used. By using these smaller, insulating molecules (usually based on alkyl chains) as

spacers, the molecules of interested are separated spatially from each other enhancing the

probability to contact a single molecule with the STM tip. As di-thiol SAMs (with sulfur

groups on both ends) are difficult to produce, the molecules can usually only be chemically

bound to the substrate, the second contact then is a tunneling contact. In consequence,

STM based single-molecule measurements always have strongly asymmetric coupling to

the electrodes, which makes interpretation of the experimental data more complicated.

The first demonstration of a single-molecule measurement based on STM techniques

was done by Bumm and Tour [33]. Here, longer molecules consisting of several benzene

rings in a SAM with shorter, insulating molecules as spacers, were contacted using the

STM tip. The molecules were covalently bound to the substrate via thiol groups, whereas

the contact to the STM tip was a tunneling contact. Since then, many other STM based

conductance measurements have been reported [23,32,34]. In 2003, Tao and coworkers

1.2. Inelastic electron tunneling spectroscopy

developed a method that allows the repeated formation of molecular junctions using a gold

STM tip in fast succession [32]. The STM tip is moved into and out of contact with a

gold substrate in a solution containing the sample molecules. When the tip is pulled out

of contact, a stepwise decrease in the conduction is observed, (Fig. (1.5)). These steps are

the usual conduction steps observed in metallic point contacts. When the tip is pulled out

further, additional plateaus appear, corresponding to a resistance two orders of magnitude

higher. By statistical analysis it was possible to attribute these steps to the measurement of

a distinct number of molecules between the contacts.

1.2 Inelastic electron tunneling spectroscopy

IETS is an all electronic spectroscopy that has been extensively reviewed [35,36]. By mea-

suring currents and voltages across a metal-adsorbate-metal (M-A-M’) device in a special

way, one is able to extract vibrational and electronic spectroscopic information about the

metals (magnons and phonons), the insulator, and the adsorbate. It has been successfully

applied to problems in surface chemistry and catalysis, adhesion and corrosion, molecular

vibrational and electronic spectroscopy. The motivation for tunneling studies is based on

several unusual properties of the electron-molecule scattering process. Despite the com-

plexity of the experimental setup, the molecule must be connected to two electrodes, there

are some remarkable advantages (see [35]):

•Optically forbidden (IR and/or Raman) fundamental transitions may be observed. In

addition, spin and dipole forbidden electronic transitions are also allowed .

•Overtone and combination bands are very weak or absent. Thus, virtually every band

is due to a fundamental transition.

•Very small amounts of compound can provide a complete spectrum.

•IETS offers the opportunity of studying molecules in buried interfaces thereby pro-

viding unique insights into the effects on molecular states produced by intimate con-

tact with metallic and insulating layers.

•IETS provides a high sensitivity method for studying chemisorption and catalysis.

•It is an all solid-state spectroscopy. In principle the sample and spectrometer could

be built on a single chip. When this feature is combined with the further development

of porous electrodes, IETS becomes uniquely situated for sensor applications.

The IETS is based on electron tunneling. By electron tunneling I mean the motion

of electrons from one classically allowed region to another through a region where the

electron is classically forbidden to exist (see Fig. (1.6)). The barrier region is one in which

the potential energy, U, is greater than the total classical energy, E. If the particle moves

from one allowed region to another, it must tunnel through the potential barrier. Of course,

Chapter 1. Molecular Electronics: a Brief Overview

there is no such prohibition in quantum mechanics. An electron of energy Eimpinging on

the potential barrier from the left in Fig. (1.6) has an exponentially damped probability of

penetrating the barrier a distance z. If Uand Eare both large (several electron volts) and

dis very small (of the order of 1 nm), there is a finite probability of the incident electron

emerging into the classically allowed region III.

Figure 1.6: Schematic representation of the elastic tunneling process.

The tunneling can be divided as the sum of two fluxes: the first related to elastic tun-

neling when the electron flows from one contact to the other without any interaction with

the barrier, the second is the inelastic tunneling when the electron undergoes an interaction

with the degrees of freedom inside the barrier (like the modes of vibrations of the molecular

region).

A schematic view of what happen for IETS is shown in Fig. (1.7, a). In the crude model,

the free electron theory is used to represent the metals and the insulator plus adsorbate is

treated as a vacuum space. The electron states are filled up to the Fermi energy, EF, and the

work function of the metal is given by

. It is assumed that we are at very low temperature

so that the occupation in the contacts follow a real step-like function. The hatched regions

in Fig. (1.7, a) represents completely filled metal orbitals. The HOMO and LUMO of the

molecular region are also shown in Fig. (1.7, a), as are a few schematic ground vibrational

levels. The interaction between the top metal and the molecule can produce a shift in

molecular energy levels and its level spacing as well as the vibrational spectrum. As the

two metals and insulator are brought together, electrons flow via tunneling until the Fermi

surfaces are matched. If the left-hand metal electrode is taken as the reference, the right-

hand electrode develops a net charge, of sign and magnitude determined by the difference

in work functions. Once the Fermi energies are matched, no net current flows and the

approximate barrier height is

∼(

M+

M0)/2. In the case of metals with widely different

1.2. Inelastic electron tunneling spectroscopy

work functions there is an internal potential generated that can amount to millions of volts

per centimeter.

Figure 1.7: Energy diagram for an Al-Al2O3-anthracene-Pb tunnel diode showing elastic and

inelastic tunneling processes. (a) The hatched regions represent the filled states of the top and

bottom metal electrodes. The area in the center represents the anthracene coated alumina. The

HOMO (

) and LUMO (

∗) orbital energies and a few vibrational levels are indicated. The case

shown is where the bias energy (eV) is just sufficient to allow inelastic tunneling with excitation of

the first vibrational level, eV = h

. Energy loss (equilibration) for the tunneling electron occurs

through a cascade process in the M’ electrode. (b) The I-V curve, conductance-V curve, and the

IETS spectrum that would result from both elastic processes and the first inelastic channel as shown

in the barrier diagram.

When a small bias is applied a tunneling current starts to flow. At low temperature

the absorption of phonons is much smaller than the emission from the electrons which is

the dominating process. Until the bias is smaller than the smallest vibrational frequency

no scattering is allowed because after the emission of a phonon the electron will end up

below the Fermi energy of the second electrode. This means that only the elastic current

can pass. This elastic tunneling current is the ever-present source of the background signal

with which all tunneling spectroscopists must deal. Since IETS depends on derivatives of

Iwith respect to V(see Fig. (1.7, b)) it would first appear that the elastic current would

have little effect on the IETS. Unfortunately, the number of electrons tunneling elastically

is usually orders of magnitude larger than the number utilizing inelastic channels. Thus,

relatively small deviations of the elastic current produce large background changes. This

background is depicted as linear in Fig. (1.7, b), but can become nonlinear and even non-

monotonic.

Chapter 1. Molecular Electronics: a Brief Overview

When the bias is large enough, in addition to elastic tunneling, there are other tunneling

mechanisms which may contribute to the current: the inelastic current as shown in Fig. (1.7,

b). The moving electronic charge interacts with the time-varying molecular vibrations to

induce excitation of the molecule in the barrier with concomitant loss of energy by the

electron. From an energetically point of view the applied voltage is less than h

/e, the

inelastic channel is closed because the final states are already filled (Pauli principle). At

V=h

/ethe inelastic channel opens. Further increases in Vresult in additional possible

final states with an associated increase in current due to this channel. As is depicted in

Fig. (1.7, b), there is a change in the slope in the I-V curve at V=h

/e. If one measures

the conductance, dI/dV, the opening of the inelastic channel is signaled by a step. Plotting

d2I/dV2versus Vproduces a peak at V=h

/e. Both vibrational and electronic transitions

may be observed as peaks in the d2I/dV2versus Vplots.

The width of the peaks in IETS depends upon the sharpness of the onset of the inelastic

process, which in turn depends upon the thermal distribution of electron energies about EF.

Thus, the IETS line width depends on temperature and is about 3.5 cm−1K−1. Because of

this, vibrational IETS is most often performed below 5 K. Electronic transitions are usually

much broader than vibrational transitions. In its simplest form, an IET spectrum is a plot

of d2I/dV2versus V. It turns out that using (d2I/dV2)/(dI/dV)as the y-axis provides

spectra having flatter baselines and is most appropriate for high bias work [37].

There are very few hard selection rules in IETS. Instead, there are some selection pref-

erences (referred in literature as propensity rules [38]). In IETS we see some bands that

appear in IR or Raman spectra and others that are optically forbidden. Thus, IETS is an ex-

cellent complementary vibrational spectroscopy for high symmetry molecules where many

of the modes are optically forbidden.

In principle, the spectrum in IETS is usually assigned by comparison with IR, Raman,

and HREELS results for monolayers of the molecule in question. However, the absence of

simple selection rules means that IETS may exhibit markedly different spectra with charac-

teristics that are difficult to predict. There have been several theoretical approaches devel-

oped to simulate IETS [39,40,41,42,43,44,45,46], approaches with varying strengths

and weaknesses. The theoretically most advanced ones have been applied to model sys-

tems [43] while many of the earlier ones made significant approximations in order to treat

real molecules with ease.

Chapter 2

Modeling Molecular Electronic Devices

The main aim of every theoretical transport model is of course the description of the flow

of current in the device region, in our case focusing on the dissipation induced by electron-

phonon interaction (vibronic coupling). In Chapter 1some of the most popular experi-

mental setups for molecular electronic devices were presented. The reason for that brief

introduction was to show the kind of systems that must be simulated. Following those

setups, in this second Chapter, a theoretical approach to model the device is discussed.

In the first part of this Chapter is shown how the device is partitioned in different re-

gions. This partitioning is related to the way the bias drops through the entire device. First

we have the proper molecular bridge connected to two leads, the leads are the interface

between the molecular bridge and the bulk contacts. The contacts are bulk regions which

act as two semi-infinite reservoirs of charges. The current flowing in the molecular bridge

is a flux of charges going from one contact to the other. This flux is influenced by many as-

pects of the system under investigation. First of all, the electronic structure of the molecule

and of the two leads, which is related to the geometry of the system and the material used

to make the device. The interface between molecule and contacts should also define the

charge transfer that occurs between the two also at zero bias and, eventually, the relax-

ation of the molecule when sandwiched between the two electrodes. Two other important

aspects influence the total current: the bias applied, and the electron-phonon interaction

which changes the distribution of electrons in the molecule and determine the power dissi-

pated in the bridge.

In the second part of this Chapter is presented the Meir-Wingreen equation [47], for

the calculation of the steady-state current including the effect of inelastic scattering. This

general equation can be simplified to a complete elastic transport equation, the Landauer-

B¨

uttiker equation [48,49], if the inelastic scattering interaction is neglected.

During this Chapter extensive use of the formalism called second quantization is made.

It is a particular convenient formalism for transport problems. For the non expert reader I

suggest, from the literature, a brief list of excellent sources about the topic [50,51,52].

Chapter 2. Modeling Molecular Electronic Devices

2.1 The general problem of molecular conduction

All the experimental setups show very similar features: the presence of two contacts and

a molecule, or a layer of molecules, which bridges them, sometimes also a gate is present.

In Fig. (2.1) a schematic view of a typical molecular electronic device is shown.

Figure 2.1: A schematic view of a typical molecular electronic device. In the picture different

regions can be recognized: the molecular (M) and the surfaces or leads regions (L) which together

form the device region (D). The two semi-infinite contacts (CLand CR)act as charge reservoirs.

To better simulate such devices, the entire system is usually split in different parts. For

the simplest geometry, we need at least two electrodes, (CLand CR, but the formalism and

the experiments can be generalized to a larger number of electrodes) and a device region

(D). The contacts (CLand CR) represent two semi-infinite reservoirs of charge carriers.

Electrons within the contacts are assumed to be in thermodynamic equilibrium and their

electron distributions can be correctly described by means of two Fermi functions with a

definite chemical potential (

Land

R), see Fig. (2.2). If no bias is applied the two chemi-

cal potentials are the same. The contacts are assumed to be perfect, defect-free crystalline

structures, in the sense that they can always be constructed starting from a defined unit

cell. Apart from this assumption, no constraint is imposed on the dimensionality, trans-

verse size and transverse shape of them. They can be one-dimensional atomic chains,

two-dimensional planar atomic strips or three-dimensional wires of arbitrary cross-section.

The device region (D) is obviously the most important part of the system. It is the

region in which the main processes happen. It can be formed by every collection of atoms

between the contacts. When a bias is applied the reservoirs are assumed to remain locally

in thermodynamic equilibrium. This means that the only effect of the bias is to change the

two chemical potentials that now have different values. All the voltage drop is instead in

the device region, which is in non-equilibrium condition (see Fig. (2.2)). The main con-

sequence of this is that both contacts try to shift the chemical potential of the device to its

own level. The contact with the higher chemical potential tries to charge the device while

the other contact tries to discharge it. The final effect is a net rate of electrons (or holes)

flowing through the device from one contact to the other, that means a current. Due to the

fact that the device region must include all the part of the system where there is a significant

2.1. The general problem of molecular conduction

Figure 2.2: The voltage drop in the different regions when an external bias is applied.

Land

are the two chemical potentials of the left and the right contact respectively. The drop of the bias

occurs practically entirely in the molecular region and to a small extent in the lead regions also.

voltage drop it cannot reduce, usually, to only the molecular part. Some layers of the con-

tacts must be included in the device, they are called leads. The system represented by the

molecule plus the leads is often called extended molecule in the literature. The size of the

leads depends mostly on the screening of the bias induced by the material of the contacts,

so that different materials can result in different leads. The partitioning between device and

contacts is an important issue, however, for many qualitative simulations we can effectively

restrict the device region to the molecular bridge only, without large consequences.

The switching of the bias induces a current. After a first transient regime, where the

current shows a clear time dependence, it reaches rapidly a steady-state. Despite the impor-

tance of the transient regime for the stability and the switching properties of an electronic

device, the steady state regime is the only one investigated in the present work, also because

is the one which is usually measured in experiments.

Theoretically the first challenge is to map the open boundary conditions represented

by the contacts into the finite size device region. The problem is solved defining two po-

tentials, called self-energies, which describe the incoming and outgoing of charges to and

from the device. This simplification allows to reduce the dimensionality of our system.

Instead of including the entire semi-infinite contacts in the simulated model I need just the

device region and some chunks of the two contacts. From these chunks the self-energies

can be computed (see Chapter 5and Appendix Bfor a complete description of the evalua-

tion of the self-energies). They are non local, non hermitian, energy dependent potentials.

The non hermitianity is the fundamental property of these potentials. Including them in the

Schr¨

odinger equation we define for the device region a new set of eigenvalues and eigen-

states. The eigenvalues are complex numbers, where the real part is the energy of the level

while the imaginary part is the inverse of a lifetime. The lifetime takes into account the

time spent by the charge in the device, before it is absorbed by the contacts.

The eigenstates represent a set of transmission states, connected through the self-energies

to the external environment; they are called scattering states. Each of these scattering states

is carrying a finite amount of current due to its occupancy. In case of molecular devices

Chapter 2. Modeling Molecular Electronic Devices

these scattering states give a simple explanation of the strong non-linearity which is nor-

mally observed in current-voltage characteristics. When in fact one of these states enters in

the energy window between the two chemical potentials of the contacts another channel for

the current is switched on and more current can flow. The knowledge of scattering states

completely covers the problem of describing quantum conductance in the case of elastic

transport, that means when inelastic scattering processes are not relevant.

In case scattering processes (e.g. electron-phonon) are present the problem is compli-

cated by the fact that emission or absorption of phonons can change the channel in which

one charge propagates and scatters it into another channel. The calculation of scattering

states can be achieved by several techniques, such as transfer matrix (TM) and Green’s

function (GF) approaches [53,54]. Although the computational cost of the TM method is

generally smaller with respect to the GF technique, the latter approach is more powerful. In

fact it can be extended and generalized to a statistical Non-Equilibrium Green’s Function

(NEGF) theory [50,55] accounting for the correct population of the states with the inclu-

sion of scattering mechanisms (electron-electron, electron-phonon, etc). A more detailed

treatment of the NEGF method is presented in Chapter 3.

2.2 The Hamiltonian of the system

The general Hamiltonian for the entire system must consider the presence of the device and

the contacts regions and, also, of the phonons in the molecule which are connected to an

external bath in thermodynamic equilibrium. The external bath of phonons represents the

thermal vibrations of the contacts. With this in mind it is possible to write the Hamiltonian

of the entire system, using second quantization formalism:

H=ˆ

HD+ˆ

VDC +ˆ

HC(2.1)

where:

HD=

∑

i=1

iˆ

d†

iˆ

di+∑

ωα

ˆa†

ˆa

+∑

i,j;

i j ˆ

d†

iˆ

dj(2.2)

HC=∑

k∈L,R

kˆc†

kˆck+∑

=L,R

ωβ

b†

(2.3)

VDC =∑

i;k∈L,R³Vki ˆc†

kˆ

di+V∗

ki ˆ

d†

iˆck´+∑

;

=L,R

αβ

(2.4)

Eqn. (2.1) is the sum of three terms: the contributions to the Hamiltonian of the device,

HD, of the contacts, ˆ

HC, and the coupling between the two, ˆ

VDC. In Eqn. (2.2), Eqn. (2.3)

and Eqn. (2.4), ˆa( ˆa†) and ˆ

b(ˆ

b†) are annihilation (creation) operators for the phonons in the

device and in the external bath respectively. The ˆ

d(ˆ

d†) and ˆc(ˆc†) are similar operators for

electrons in the device and in the contacts. The creation and annihilation operators satisfy

2.2. The Hamiltonian of the system

the usual commutation rules:

[ˆai,ˆa†

j] =

i j (2.5)

[ˆci,ˆc†

j]+=

i j (2.6)

where the first relation is for boson operators and the second, the anticommutation, is for

fermion operators. Finally, Qa

and Qb

are vibration coordinate operators:

=ˆa

+ˆa†

(2.7)

=ˆ

+ˆ

b†

(2.8)

In the equation

iis the energy state of the device central region,

ωα

(

ωβ

) are the energies

of phonons in (and out) of the device, M

i j is the electron-phonon coupling for the

mode of vibration,

kis the energy of one electron in the contacts, Vki is the contact-device

coupling and finally U

αβ

is the phonon-phonon coupling between phonon in the device and

in the external bath. Eqn. (2.2) is the sum of three terms related to the electronic structure,

the phonon population and the interaction between phonon and electrons in the device.

Eqn. (2.3) is the same for the two contacts, but because the contacts are in equilibrium, there

is no electron-phonon term here. Finally, the last equation, (2.4), describes the coupling

between contacts and device (first term) and the coupling between phonons in the bath and

in the device (second term).

The model is characterized by other parameters which are fundamental for the evalua-

tion of the current. These quantities are:

ΓL,R

i=2

∑

i j;k∈L,R|Vik|2

(

i−

k)(2.9)

it represents the molecule-contact coupling by its effect on the lifetime broadening on a

molecular level iand,

L,R

=∑

;

∈L,R|U

αβ

(

ωα

−

ωβ

)(2.10)

is similarly the lifetime broadening induced by the relaxation of phonons of the device in

the external bath. The central quantity in the electron-phonon interaction is obviously the

vibronic coupling M

i j .

In the formulation implemented in my code the last term of Eqn. (2.4) is neglected.

The phonon population is assumed to be in perfect equilibrium with the environment. This

means that all the U

αβ

coefficients are zero and that there is no relaxation of the phonons

from the device region into the contacts. This approximation reduces the applicability of

the code to low temperature simulations only, because relaxation processes are particularly

important at higher temperature. Nevertheless many experiments in molecular electronics

Chapter 2. Modeling Molecular Electronic Devices

are performed at liquid helium temperature (4.2 K); however, an extension of the formalism

including electron-phonon relaxation has been implemented recently, see [56].

2.3 From the Hamiltonian to the current: Meir-Wingreen

equation

Following the approach of Caroli et al. [57] it is possible to relate the Hamiltonian of

Eqn. (2.1) to the current. From the continuity equation it is possible to define the current

flowing through the left electrode:

IL=−eh˙

NLi=−ie

hh[ˆ

H,ˆ

NL]i(2.11)

where ˆ

His the Hamiltonian of the system and ˆ

NLis the number operator for the device

which couples with the left electrode only. It is defined as,

NL=∑

d†

kˆ

dk(2.12)

The only term in the Hamiltonian which does not commute with the number operator is

the interaction term between the device and the left contact (first term in Eqn. (2.4)). After

solving Eqn. (2.11) we get

IL=2ie

h∑

i;k∈L³Vkihˆc†

kˆ

dii−V∗

kihˆ

d†

iˆcki´(2.13)

where the factor 2 comes from the spin degeneracy. If we analyze in more details Eqn. (2.13)

we find that the total current through the left interface is the difference between two fluxes:

the outgoing, hˆc†

kˆ

dii, and the incoming, hˆ

d†

iˆcki, electrons. The direction of the flux of

charges can be easily checked considering the kind of operators involved, i.e., for the in-

coming flux hˆ

d†

iˆckiwe have that one electron is annihilated in the left contact, ˆck, and one

is created in the device, ˆ

d†

i. We can reformulate Eqn. 2.13 and terms like hˆ

d†

iˆckiby intro-

ducing the so called lesser Green’s function, which is a time dependent quantity, G<

ik(t,t0)

defined as

ik(t,t0) = ihˆc†

kˆ

diie−iE

h(t−t0)(2.14)

The G<

ik depends only on the difference of tand t0because we are computing the steady-

state current. We can define t−t0=

and Fourier transform the function G<

ik:

2.3. From the Hamiltonian to the current: Meir-Wingreen equation

ik(

) = Z+∞

−∞

ik(E)eiE

ik(E) = Z+∞

−∞

ik(

)e−iE

h(2.15)

(2.16)

The G<

ik, apart from a phase factor, represents the correlation, in energy or in time, between

one electron in the state iof the device and another one in the state kof the left contact.

The two fluxes of the current are related to these correlation functions as follows:

hˆc†

kˆ

dii=−iG<

ik(0) = −iZ+∞

−∞

ik(E)(2.17)

hˆ

d†

iˆcki=−iG<

ki(0) = −iZ+∞

−∞

ki(E).(2.18)

Substituting the latter relations in Eqn. (2.13) for the current we get,

IL=2e

h∑

i;k∈LZ∞

−∞

£VkiG<

ik(E)−V∗

kiG<

ki(E)¤(2.19)

The latter formula permits to introduce the Green’s function formalism in the transport

calculations. The main drawback of Eqn. (2.19) is that the quantities G<

ik and G<

ki depend

on operators which apply both to the device region, ˆ

d†

iand ˆ

di, and the left contact, ˆckand

ˆc†

k. What would be more desirable for the calculation is the possibility to decompose the

G<functions in combination of terms which depend on operators in only one part of the

system. This can be done using the so called Dyson’s equation. A rigorous derivation

of the Dyson’s equation is not given here, but in Chapter 3. In order to apply the Dyson’s

equation we must first define two new propagators, the retarded and advanced propagators,

Grand Ga:

im(

) = −i

(

)h[ˆ

di,ˆ

d†

m]+ie−iE

(2.20)

im(

) = i

(−

)h[ˆ

di,ˆ

d†

m]+ieiE

(2.21)

The retarded function represents a dynamic propagators of one particle from state ito state

m, while the advanced propagator describes the propagation of a charge vacancy, or a hole,

from state ito state m. The retarded, advanced and lesser Green’s functions are related by

Chapter 2. Modeling Molecular Electronic Devices

the following equations, see [57,47]:

ik(E) = ∑

V∗

km[Gr

im(E)g<

kk(E)+G<

im(E)ga

kk(E)] (2.22)

ki(E) = ∑

Vkm[gr

kk(E)G<

mi(E)+g<

kk(E)Ga

mi(E)].(2.23)

where capital letters were used for propagators which operate on the device region only

and small letters for propagators which operate on the left contact. In these last relations

we can see that, apart for the Vkm coefficient which couples both device and contact, all the

propagators depend on only one part of the system. Inserting Eqs. (2.22) and (2.23) in the

current expression and interchanging the indices mand iin the last term we obtain

IL=2e

h∑

i,m;k∈LZ∞

−∞

VkiV∗

km[g<

kk{Gr

im −Ga

im}+G<

im{gr

kk −ga

kk}].(2.24)

Because the contacts are free of any interaction their Green’s functions assume a very

simple form:

kk =2

inL

(E−

k)(2.25)

where nLis the Fermi distribution function of the left contact, and

kk −gr

kk =2

(E−

k).(2.26)

It is possible to use a generalization of ΓL

i(Eqn. (2.9)) to put all the reference about the

contact-device interface in one function,

ΓL

mi =2

∑

k∈L

VkiV∗

(E−

k)(2.27)

Substituting Eqn. (2.25), Eqn. (2.26) and Eqn. (2.27) in the current formula we get,

IL=2ie

hZ+∞

−∞

dE ∑

i,m

ΓL

mi[nL(Gr

im −Ga

im)+ G<

im](2.28)

which in matrix notation becomes

IL=2ie

hZ+∞

−∞

Tr[ΓLnL(Gr−Ga)+ΓLG<]dE.(2.29)

The total current is the difference between the charges entering through the left contact and

leaving through the right one

I=1

2(IL−IR)(2.30)

2.3. From the Hamiltonian to the current: Meir-Wingreen equation

so we get finally for the current between the left and right junctions

I=ie

hZ+∞

−∞

Tr[(nLΓL−nRΓR)(Gr−Ga)+(ΓL−ΓR)G<]dE.(2.31)

In a steady state the flux between the left junction must be equal in magnitude and just

opposite in sign: IL=−IR. This means that the total current is just equal to the left current.

Moreover, we can define two new quantities which represents the injection (extraction) of

charges through the left contact

Σ<

L=inLΓL(2.32)

Σ>

L=−i(1−nL)ΓL(2.33)

Substituting the former in the current equation and using the identity (see Chapter 3)G>−

G<=Gr−Ga, we get, after some manipulations, the Meir-Wingreen equation

I=2e

hZ+∞

−∞

Tr[Σ<

LG>−Σ>

LG<]dE.(2.34)

Further simplifications can be done to the Meir-Wingreen equation if there is no scattering

in the device region, in that case the formula of the current reduces to the famous Landauer-

B¨

uttiker equation. First of all, without scattering, we have the following, see [53],

G<=inLGrΓLGa+inRGrΓRGa(2.35)

Gr−Ga=−iGr(ΓL+ΓR)Ga(2.36)

and substituting them in the Meir-Wingreen equation,

I=2e

hZ+∞

−∞

Tr[ΓLGrΓRGa](nL−nR)dE.(2.37)

In the Landauer-B¨

uttiker equation the current is split in two parts: a transmission factor

which depends on the coupling with the contacts and the propagators and an occupation

factor which depends on the distribution of electrons in the contacts only. In other words,

in elastic transport the occupation is dominated by the contacts and the dynamics mostly

by the device. In the inelastic case, described by Meir-Wingreen, the current is just the

difference between two fluxes (in and out of the device region), where the occupation

information and the dynamics information are strongly correlated due to the scattering.

Chapter 2. Modeling Molecular Electronic Devices

2.4 Lifetimes of interest

The equation displayed in (2.34) is the basis of the non-equilibrium evaluation of the cur-

rent. However, despite the Meir-Wingreen evaluates a steady state current, time has still a

central role.

The Hamiltonian and the Meir-Wingreen formulas are in fact related to many pa-

rameters which are connected to the lifetime of all the processes involved in transport

through scattering and some energy parameters. These lifetimes are usually very different

in timescale and so the effect of one process on the full transport can change dramatically

depending on the system under investigation. The final result is that the effect of electron-

phonon interaction in molecular conduction can produce a large set of different behaviors

observed experimentally.

The first important energy parameter is ∆E, the energy difference between the Fermi

energy of the electrode and the closest molecular orbital (the highest occupied molecular

orbital (HOMO) or the lowest unoccupied molecular orbital (LUMO)). The second pa-

rameter is Γ, the broadening of molecular orbital levels induced by the coupling with the

contacts. If the device behaves like a molecular chain, there is another important value, the

bandwidth VBof the “conduction band” of the chain. Finally, the energy value kBT(where

kBis the Boltzman constant and Tis the absolute temperature) is quite important giving

the thermal effects to the electronic and phonon populations.

All these energy parameters are related to some time constants which are fundamental

for the transport. First of all there is the transmission time (¯

h/Γk) which is the time an

electron in the state kspends in the device. Second, we have the relaxation time for phonons

(¯

γα

) which is the time needed by a phonon in the device to relax in the thermal bath into

the contacts. This contribution is neglected due to the fact that the system is constrained

to be always in thermal equilibrium with the bath. Finally ¯

h/VBis the lifetime spent by an

electron in one molecular site if the device behaves like a wire.

Another important time scale, and that is less obvious, is the so called dephasing time

(

) which describes the time needed by the modes of vibration to destroy completely the

phase coherence of the electron. When

¿¯

h/Γkthen the transport reduces to a collection

of hopping from one site of the device to the other. On the other hand if the dephasing

does not occur, the transport can be treated using elastic transport and adding the electron-

phonon contribution as a perturbation.

The relation between all these energy parameters and time scales defines which model

is more appropriate for the system under investigation. The comparison between dephasing

time and transmission time decides if the transport is fundamentally an elastic coherent

tunneling or it is more a collection of jumps. When either ∆Eor Γkare much larger

compared to the electron-phonon interaction the transport becomes fundamentally elastic

(especially if ∆E≈0, that means at resonance). In an intermediate case we can use a

perturbation correction to introduce the inelastic part of the current induced by the device

vibrations. If instead Γkis small (weak coupling regime) the device behaves fundamentally

as a dot and also if the electron-phonon coupling is small there is time for the formation

of polaron-like states inside the device. The transmission is basically a hopping between

the device and the two contacts assisted by vibrations. These different cases can be easily

2.4. Lifetimes of interest

tested by a comparison between the different quantities,

¯¯¯¯¯

p(∆E)2+(Γk/2)2¯¯¯¯¯

=∆.(2.38)

The two limits (∆¿1 and ∆À1) are called strong and weak electron-molecule coupling

regime respectively. In the present work we investigate the first case only.

The discussion of course is not confined to only these parameters, that anyway remain

the most important. Other aspects can play a fundamental role. One for all is the possibility

that strong many body effects can change the picture. When in fact a charge enters the

device there is a distortion. If the transmission time is long enough and the distortion is

localized, this creates a quasi-particle state, a polaron, which propagates in the system. The

different nature of the polaron, basically an electron dressed by a collection of phonons,

can change dramatically the conduction. Moreover, the change induced in the device can

change the time and energy scale previously described changing eventually the type of

conduction.

Chapter 3

Non-Equilibrium Green’s Functions

The main result of Chapter 2is the Meir-Wingreen equation (Eqn. (2.34)), which gen-

eralizes the Landauer-B¨

uttiker equation for elastic transport. The Landauer formula is

in fact contained in the Meir-Wingreen as a special case, when scattering processes are

neglected. We can expect that Landauer is a good approximation of transport for short

devices, strongly bound to the contacts, at low temperature. In these conditions, in fact, the

probability that scattering processes occur is reduced. However, also in that case, inelastic

scattering can play an important role. In Chapter 2, for the evaluation of the current, a

collection of propagators were introduced, related by a set of so called Dyson’s equations.

These propagators are related to the statistical population of electrons in the system like the

lesser (G<(E)) Green’s functions or to their dynamics like retarded (Gr(E)) and advanced

(Ga(E)) propagators.

Their evaluation requires the Non Equilibrium Green’s Function (NEGF) formalism.

NEGF methods represent the state of the art for molecular electronics transport calcula-

tions. This is based on two facts: first, the open boundary conditions of the problem (the

connection between the device and the semi-infinite contacts) can be easily handled in the

calculations. Second, scattering processes can be included in the formalism following a

general and systematic way based on a perturbation expansion of the Green’s functions.

The NEGF formalism relies on the works of Kadanoff, Baym and Keldysh [55,58] mainly.

Many body calculations with NEGF are often done for simulating the behavior of sys-

tems at zero temperature. Of course, real experiments are never performed at zero tem-

perature, although, in molecular electronics, many of them are often at low temperature.

However, many quantities are not very sensitive to temperature and the effect of the latter

can be included as a perturbation to a full zero temperature approach. Furthermore, the

zero temperature property of a system is an important conceptual quantity, the ground state

of an interacting system. A system is often described as its ground state plus its excitations,

and the ground state may be deduced from a zero temperature calculation.

The starting point of every Green’s function calculation is the Hamiltonian of the prob-

lem ˆ

H. It is presumed that one is trying to solve a Hamiltonian which cannot be solved

Chapter 3. Non-Equilibrium Green’s Functions

exactly. The usual approach is to set

H=ˆ

H0+ˆ

V(3.1)

where ˆ

H0is a Hamiltonian which may be solved exactly. In our approach the ˆ

H0can be

split in the sum of three contributions:

H0=ˆ

0+ˆ

ΣL+ˆ

ΣR(3.2)

Where the first term is a single particle Hamiltonian for the device region represented by

a Density Functional Theory (DFT) Hamiltonian, the description of the particular DFT

method applied in our case is the topic of Chapter 4. The other two terms (ˆ

ΣL,R) are the

contributions from the open boundary conditions (the way of computing them is described

in Chapter 5and Appendix B). The last term ˆ

Vrepresents all the remaining parts of ˆ

H, in

our case the effects of the electron-phonon interaction.

3.1 Three representations

In this brief paragraph the three most common representations for operators and wave

functions are presented. They are the Schr¨

odinger, the Heisenberg and the interaction

representation. The paragraph is just a short refresh, for more details the reader is referred

to other books which fully treat the topic [52,51]. All the operators are marked by a hat

symbol (Schr¨

odinger and Heisenberg representations), except interaction representation

for which a tilde is used. From now and through all this work, apart from a few sections,

atomic units will be used for which ¯

h=me=e= 1 (see Appendix A).

3.1.1 Schr¨

odinger representation

Elementary quantum mechanics is usually presented in Schr¨

odinger representation, which

is based on the formula (in atomic units):

∂ψ

(t)

∂

t=ˆ

(t)(3.3)

which has the formal operator solution:

(t) = e−iˆ

(0)(3.4)

The use of this formula requires two assumptions:

1. The wave functions

(t)are time dependent, even if this time dependence is a simple

factor of exp(−iEt).

2. Operators, such as the Hamiltonian ˆ

H, are independent of time.

3.1. Three representations

3.1.2 Heisenberg representation

It is possible to solve quantum mechanical problems in another way which gives exactly

the same answers. The Heisenberg representation has the following properties:

1. The wave functions are independent of time.

2. The operators are time dependent, and this dependence is given by,

O(t) = eiˆ

Ht ˆ

O(0)e−iˆ

Ht (3.5)

or, equivalently, one is trying to solve the equation of motion which is derived from

a time derivative:

∂

O(t)

∂

t= [ ˆ

O(t),ˆ

H].(3.6)

3.1.3 Interaction representation

The interaction representation is the third way of recasting all quantum mechanics. Here

both the wave functions and the operators are time dependent. The Hamiltonian is sepa-

rated in two parts,

H=ˆ

H0+ˆ

V(3.7)

where ˆ

H0is the unperturbed (easy) part, while the ˆ

Vcontains the interactions. In interaction

representation,

1. Operators have a time dependence

O(t) = eiˆ

H0tˆ

Oe−iˆ

H0t.(3.8)

2. Wave functions, ˜

have a time dependence

(t) = eiˆ

H0te−iˆ

(0) = e−iˆ

(0) = ˆ

U(t)

(0).(3.9)

It is assumed that [ˆ

H0,ˆ

V]6=0. The main result of this division is that the time dependence

of the wave functions is related to the ˜

Vonly, the complicate part, while all the operators

depend on the unperturbed ˆ

H0. It is important to note that ˆ

H0is the same in Schr¨

odinger

and interaction representations.

The operator ˆ

U(t), in Eqn. (3.9), reduces to unity for t=0. It can be represented as

a series expanding the exponential operators. We define first an operator ˆ

Tcalled time-

ordering operator particular useful in this case:

T[˜

A(t1)˜

B(t2)] =

(t1−t2)˜

A(t1)˜

B(t2)±

(t2−t1)˜

B(t2)˜

A(t1)(3.10)

Chapter 3. Non-Equilibrium Green’s Functions

where

is the usual step function and the upper (lower) sign is for Boson (Fermion) oper-

ators. The function of ˆ

Tis to time order the operators following their statistic (symmetric

or antisymmetric). The operator ˆ

Ucan be expanded as

U(t) = 1+

∞

∑

n=1

(−i)n

n!Zt

dt1Zt

dt2...Zt

dtnˆ

T[˜

V(t1)˜

V(t2)... ˜

V(tn)] (3.11)

U(t) = ˆ

Texp·−iZt

dt1˜

V(t1)¸(3.12)

where the last form is a short way to represent the infinite series.

3.2 S-Matrix

We define a new operator called S-Matrix, ˆ

S(t,t0), which propagates the wave function

(t0)in time:

(t) = ˆ

S(t,t0)˜

(t0) = ˆ

S(t,t0)ˆ

U(t0)˜

(0).(3.13)

This leads to the conclusion that the operator ˆ

S(t,t0)has the following formula:

S(t,t0) = ˆ

U(t)ˆ

U†(t0)(3.14)

From the previous section and Eqn. (3.14) we can define a list of properties for the ˆ

operator:

1. ˆ

S(t,t) = 1

2. ˆ

S†(t,t0) = ˆ

S(t0,t)

3. ˆ

S(t,t0)ˆ

S(t0,t00) = ˆ

S(t,t00)

Furthermore, from the properties of the ˆ

U(t)operator it is possible to define a simple

exponential form for the ˆ

Soperator, similar to the expansion for ˆ

U(t):

S(t,t0) = ˆ

Texp·−iZt

t0dt1˜

V(t1)¸(3.15)

The wave function ˜

(0)≡

(0)is equivalent in both interaction and Heisenberg form and

it represents the ground state at time zero. The main problem of this formulation is that we

do not know this function because we cannot solve the exact Hamiltonian ˆ

H, even for its

lowest eigenvalue and eigenvector. In order to overcome this problem we can try to connect

the exact ground state of ˆ

Hwith the ground state of ˆ

H0,

0. The relationship between the

3.3. Equilibrium Green’s functions

two ground states ˜

(0)and

0at zero temperature was established by Gell-Mann and Low

[59]:

(0) = ˆ

S(0,−∞)

0(3.16)

Thus we are really asserting that ˜

(t)is equal to

0for (t→−∞). The traditional argument

is that one starts in the dim past with a wave function

0which does not contain the effects

of the interaction ˜

V. The operator ˆ

Sbrings this wave function adiabatically up to the present

t=0. The wave function ˜

(t)for (t→+∞) is related to the non interacting ground state

by just an arbitrary phase factor exp(iL):

(∞) = ˆ

S(∞,−∞)

0=eiL

0.(3.17)

3.3 Equilibrium Green’s functions

In this work we are interested in the propagation of two particles only: electrons and

phonons. So we will focus our attention to the Green’s functions for them. They rep-

resent the two typical cases (Fermionic and Bosonic particles) of propagators because the

basic method used for computing them can be applied, more or less, for all other particles

or excitations. At zero temperature and for steady state (important for the time dependence)

the electron Green’s function is defined as follows:

,t−t0) = −ih|ˆ

T[ˆc

(t)ˆc†

(t0)]|i (3.18)

where

is the set of quantum numbers for the problem of interest. The ˆc(ˆc†) are the usual

annihilation (creation) operators.

At zero temperature the state |imust be the ground state. If we have chosen the Hamil-

tonian of the problem to be ˆ

H, then |iis the ground state of ˆ

H. The problem arises because

the exact ground state of ˆ

His not known. The ground state |i is independent of time and

the creation and annihilation operators are in Heisenberg representation.

For t>t0we have, using the ˆ

Toperator:

,t>t0) = −ih|ˆc

(t)ˆc†

(t0)|i (3.19)

Here one takes the real ground state, and at a time t0one creates an excitation

. At a later

time tone destroys the same excitation. If

is an eigenstate of ˆ

Hthen this state would

propagate with a simple exponential time dependence:

,t>t0) = −ie−i(t−t0)(

ελ

−

0)(3.20)

For the other time arrangement, t<t0, we have that, on the contrary, a particle is removed

from the system. So the Green’s function describes the propagation of a hole:

,t<t0) = ih|ˆc†

(t)ˆc

(t0)|i (3.21)

Chapter 3. Non-Equilibrium Green’s Functions

So for t>t0the Green’s function is related to the retarded Green’s function and for t<t0

to the advanced propagator.

The problem in evaluating the Green’s function arises because both the ground state

and the operators depend on the full Hamiltonian. The idea is to isolate the effect of the ˆ

operator, so to say of the interaction part. This can be done switching into the interaction

representation. Substituting Eqn. (3.16) in the Green’s function equation and after some

manipulation, we get the following:

,t−t0) = −ih

0|ˆ

T[˜c

(t)˜c†

(t0)ˆ

S(∞,−∞)]|

0|ˆ

S(∞,−∞)|

0i(3.22)

The advantage of this representation is that now all the information about the interaction

is in the ˆ

Soperator. On the other hand the creation and annihilation operators and the

ground state are the not-interacting ones and do not contain any information about ˜

V. The

Toperator automatically handle the right sorting of all the operators between the square

brackets splitting the S-matrix in the proper way, following the properties of the ˆ

Soperator

enlisted in the previous section.

A Green’s function can also be defined for the special case where the interactions ˜

V=0

and hence the ˆ

S-matrix is unity. This Green’s function plays a special role in the formalism,

and we designate it by G0:

G0(

,t−t0) = −ih

0|ˆ

T[˜c

(t)˜c†

(t0)]|

0i(3.23)

G0is often called the unperturbed Green’s function.

3.4 Wick’s theorem

The Green’s function is evaluated by expanding the S-matrix in a series (Eqn. 3.15) such

as:

,t−t0) =

∞

∑

n=0

(−i)n+1

n!Z+∞

−∞

dt1...Z+∞

−∞

dtnh

0|ˆ

T[˜c

(t)˜

V(t1)˜

V(t2)... ˜

V(tn)˜c†

(t0)]|

0|ˆ

S(∞,−∞)|

0i(3.24)

Let us, for the moment, ignore the denominator. We are more interested how evaluate

time-ordered brackets like:

0|ˆ

T[˜c

(t)˜

V(t1)˜

V(t2)... ˜

V(tn)˜c†

(t0)]|

0i(3.25)

Our interaction is represented by the electron-phonon operator like the last term in Eqn. (2.2)

in Chapter 2. This means that our brackets contain simultaneously operators for phonons

(˜

) and for electrons (˜c

, ˜c†

). This is not a problem, it is possible to separate them

3.4. Wick’s theorem

because they commute. Taking, for example, a term like the following,

0|ˆ

T[˜c

(t)˜c†

1(t1)˜

(t1)˜c

2(t2)˜c†

3(t3)˜

(t2)]|

0i(3.26)

It can be separated as a product of terms which contain only operators of one type:

0|ˆ

T[˜c

(t)˜c†

1(t1)˜c

2(t2)˜c†

3(t3)]|

0ih

0|ˆ

T[˜

(t1)˜

(t2)]|

0i(3.27)

In the rest of the paragraph the terms for electron operators are the only ones considered

due to the fact that the brackets with phonon operators behave exactly in the same manner

except for the effect of the ˆ

Toperator (see Eqn. (3.10)).

In the electron factor of the time ordered bracket there are two creation and two an-

nihilation operators. It is rather complicate to evaluate this bracket. In order to solve the

problem we can notice at first that all brackets like Eqn. (3.27) always contain the same

number of creation and annihilation operators. Thus one is always trying to evaluate the

product of ncreation operators and nannihilation operators averaged by the non interact-

ing ground state |

0i. The effect of a creation operator c†

p(t)is to put an electron into the

state

p. In order to bring back the system into the non interacting ground state one of the

annihilation operators c

p(t)must destroy the state

pbefore the bra h

0|. For example,

the four operator term discussed before equals zero unless

1and

3or unless

3and

2. Another important help comes from Wick’s theorem. It states that a

time-ordered product of operators can be split in a sum of products of brackets which con-

tain only a pair of operators (pairing). Of all the possible arrangements in pairings only few

are not zero. Applying the theorem to the bracket in Eqn. (3.27) we obtain the following,

0|ˆ

T[˜c

(t)˜c†

1(t1)˜

(t1)˜c

2(t2)˜c†

3(t3)˜

(t2)]|

0i=

δλλ

δλ

δαβ

0|ˆ

T[˜c

(t)˜c†

1(t1)]|

0ih

0|ˆ

T[˜c

2(t2)˜c†

3(t3)]|

0i×

× h

0|ˆ

T[˜

(t1)˜

(t2)]|

0i−

δλλ

δλ

δαβ

0|ˆ

T[˜c

(t)˜c†

3(t3)]|

0i×

× h

0|ˆ

T[˜c

2(t2)˜c†

1(t1)]|

0ih

0|ˆ

T[˜

(t1)˜

(t2)]|

0i(3.28)

Note that there is a time-ordering operator ˆ

Tin each pairing bracket. For 3 creation and

annihilation operators there are six possible pairings; for noperators there are n! possible

pairings.

A few simple rules should be kept in mind when making these pairings. The first

is that a sign change occurs (only for Fermion operators) each time the position of two

neighboring operators are interchanged.

The second rule is a method of treating the “time ordering” of two operators which

occur at the same time, such as

0|ˆ

T[˜c†

(t)˜c

1(t)]|

0i(3.29)

Chapter 3. Non-Equilibrium Green’s Functions

In these cases the annihilation operator always goes to the right,

δλλ

0|˜c†

(t)˜c

1(t)|

0i=

δλλ

1nF(3.30)

where nFis just the average occupation probability for electrons the

state and it is inde-

pendent of time. When two electron operators have different time arguments in a pairing,

we put the creation operator to the right:

0|ˆ

T[˜c†

(t)˜c

1(t1)]|

0i=

δλλ

0|ˆ

T[˜c

(t)˜c†

1(t1)]|

0i(3.31)

This term can be immediately identified as the unperturbed Green’s function iG0(

,t−t1).

Our previous examples (Eqn. (3.27)) can also be written in terms of Green’s functions:

0|ˆ

T[˜c

(t)˜c†

1(t1)˜

(t1)˜c

2(t2)˜c†

3(t3)˜

(t2)]|

δλλ

δλ

δαβ

iG0(

,t−t1)iG0(

2,t2−t3)D

0(t1−t2)+

−

δλλ

δλ

δαβ

iG0(

,t−t3)iG0(

2,t2−t1)D

0(t1−t2)(3.32)

where D0is the analogous of the unperturbed propagator for phonons. In summary, Wick’s

theorem tells us that a time-ordering bracket may be evaluated by expanding it into all

possible pairings and that each of these pairings will be an unperturbed propagators or a

number operator nF. In this way we can rewrite all the terms in the perturbation expan-

sion of Gas terms which depend only on the unperturbed propagators and equilibrium

occupation probabilities.

3.5 Feynman’s diagrams

Feynman introduced the idea of representing the kind of terms in Eqn. (3.32) by diagrams.

These diagrams, are extremely useful for providing an insight into the physical process

which these terms represent. They can be drawn both for the Green’s function depending

on time G(

,t)as well as for functions which are Fourier-transformed and depend on

energy G(

,E).

The diagrams in the time space are drawn by representing the electron Green’s function

G0by a solid line which goes from t0to t, as shown in Fig. (3.1). An arrow is often in-

cluded to represent the time direction. The phonon Green’s function D

0is represented by a

wiggly line. Usually double lines are used to represent the full propagators and single lines

for unperturbed propagators. The phonon propagator have not a directional arrow because

the ˜

operator for phonons is a linear combination of both annihilation and creation oper-

ators. This mean that a phonon propagator represents simultaneously an emission (forward

arrow) and an absorption (backward arrow). The factor nFis represented as a closed loop

in time.

By using these pictorial elements, we can construct the diagrams for all the terms like

3.5. Feynman’s diagrams

Figure 3.1: Diagrams for the three basic quantities: the unperturbed propagator for electrons

(solid line) and for phonons (wiggly line) and the closed loop for the electron distribution nF(Fermi

function).

Eqn. (3.32) in our expansion of the full propagators (Eqn. (3.24)). For example, for the

second order term, expanded with Wick’s theorem, we get six contributions. These six

terms are shown in Fig. (3.2) in their diagram representations.

In Fig. (3.2) diagrams like (c),(d), and (e)are usually zero. In fact in solid state

of physics, for periodic systems, such kind of contributions represent a translation of the

crystal or a permanent strain, and neither of these is meant to be in Hamiltonian. Despite

that for molecular electronics, for which the region with electron-phonon interaction is

only the device region and the system is not invariant after a translation, terms like those

give a non vanishing contribution (see [60]). However, the two terms of primary interest

remain (a)and (b). They look alike except for the labeling of the integration variables so

they give the same contribution.

Finally the terms (d)and (f)show an interesting characteristic: the diagrams are

formed by two distinct parts, they are disconnected diagrams. There is a relevant theo-

rem, called Linked Cluster Theorem, which states that the denominator of the formula of

the Green’s function, in Eqn. (3.22), h

0|S(t,t0)|

0i, cancels with all the disconnected con-

tribution in the numerator. This means that of all the diagrams we must take care of only

those terms which produce connected diagrams. This is an important simplification of the

Green’s function expansion: not only it reduces the number of terms we must evaluate for

every order, but also delete some of the most problematic terms from the series, in fact the

disconnected diagrams are often affected by upsetting divergencies.

The evaluation of diagrams is really useful and very simple, much simpler than taking

care of the analytic counterpart. Moreover there is a dual relationship between all the non

zero pairings of terms in the series of the full propagator and the diagrams set. In principle

Chapter 3. Non-Equilibrium Green’s Functions

Figure 3.2: The six non null contributions for the second order expansion of the full propagator, in

diagram form.

to evaluate the Green’s function one should expand the S-matrix, decide which terms are

connected and which are zero, and finally obtain the steady state Green’s function by a

Fourier transform of the sum of the relevant diagrams. But this laborious procedure can be

avoided because the diagrams can be written down directly by following few simple rules.

These rules are enlisted here for a norder diagram:

1. Draw nwiggly parallel lines and two external dots;

2. Join all wiggly lines and external dots with solid lines and define on these solid lines

a direction with an arrow. Join them in all linked distinct topological way;

3. For every wiggly line: from each vertex point one solid line must enter and one must

leave;

4. The external dots are connected to only one solid line; for one of the dot it must be

an entering line and for the other a leaving line.

Once we have all the distinct linked diagrams they can be translated in analytic form with

another simple set of rules:

1. If there are many modes of vibration the diagram must be preceded by a sum over

all the modes for every wiggly line;

2. For each solid line, introduce an electron free propagator, G0;

3. For each wiggly line insert a free phonon propagator, D

0and multiply for the proper

electron-phonon coupling M

λλ

1for every vertex of the wiggly line;

3.6. Dyson’s equation

4. Conserve energy and quantum numbers at each vertex. Thus both electron lines and

phonon lines have their variables labeled to conform with this rule.

5. Sum (or integrate) over internal degrees of freedom: quantum number, energy, and

so on. All labels are internal except the two connected to the fixed external dots;

6. Finally, multiply the result by the factor:

P=(i)m

)4m(−1)F(2S+1)F(3.33)

where Fis the number of closed fermion loops and Sis the spin index (usually for

electrons (2S+1) = 2). The index mis chosen as follows:

•For electron diagrams, mis the number of internal phonon wiggly lines.

•For phonon diagrams, mis one-half the number of vertices.

3.6 Dyson’s equation

Up to now the formalism shows how to describe the full many-body Green’s function in

a series of terms (via the Wick’s theorem) which depend on unperturbed propagators only.

The zero order Green’s functions are supposed to be known exactly because they depend

on the ˆ

H0part of the Hamiltonian.

Moreover the technique of diagrams, introduced in the previous section, permits to

simplify dramatically the evaluation of the series. The diagrams in fact can be drawn

following very simple rules and represent a straightforward picture to understand the kind

of contribution that every term gives to the final sum in a physical sense. On top of that

from the connectivity of the diagrams we can see immediately which diagrams give a final

contribution to the sum and which instead cancel with the denominator (Linked Cluster

Theorem). It is possible now to reformulate the entire many body Green’s function in a

very nice formula called Dyson’s equation. The presentation of the Dyson’s equation and

of the so called self-energy in the present work is not rigorous and I refer the interested

reader to other books [50,51,52,53].

First of all, from the drawing rules for diagrams, we know that every diagram has a

pair of external propagators which connect the wiggly lines with the external dots. Let’s

consider two external electron propagators the variables of which are not summed and

integrated. These two propagators are linked by a complicate series of other electron and

phonon propagators. This means that we can rewrite the series for the full propagator

(Eqn. (3.24)) as shown in Fig. (3.3, a). The labels are in energy after a Fourier transform.

The new formulation shows that the perturbation theory for the Green’s function can

be written as:

,E) = G0(

,E)+G0(

,E)˜

Σph(E)G0(

,E)(3.34)

Chapter 3. Non-Equilibrium Green’s Functions

Figure 3.3: The diagrammatic representations of Eqn. (3.34) and Eqn. (3.35).

where ˜

Σph represents all the contributions between square brackets in Fig. (3.3, a). ˜

Σph

can be rearranged in a nicer way. First of all we can define two classes of diagrams:

the first class is represented by the so called reducible diagrams. A reducible diagram

is a diagram that can be split in other two more elementary diagrams just removing one

electron propagator (a solid line). An example of a reducible contribution is represented by

the third diagram inside the bracket in Fig. (3.3, a). All diagrams which are not reducible

are called irreducible (see for example the first two diagrams in the same figure). The sum

over square brackets for ˜

Σph can be rearranged tacking only the irreducible diagrams into

account. If we do so the Fourier transform for the full Green’s function becomes:

,E) = G0(

,E)+G0(

,E)Σph(E)G(

,E)(3.35)

The difference between Eqn. (3.34) and Eqn. (3.35) is that now Σph contains only the

irreducible diagrams. This equation is called Dyson’s equation. In diagram form it looks

like in Fig. (3.3, b). Eqn. (3.35) can be rearranged as follows

,E) = (G−1

,E)−Σph(E))−1.(3.36)

The Dyson’s equation is the fundamental equation for the solution of the propagator as

much as the Meir-Wingreen equation is for the current. It states that the exact Green’s

function is obtained from Eqn. (3.35) by just calculating the self-energy Σph. The self-

energy is a summation of an infinite number of distinct diagrams. This method is only

useful if we can approximate Σph by the lowest few irreducible terms. Σph(E)is a sort of

potential that includes all the effects of the phonons over the propagation of the electrons.

Except in a few rare cases, it is impossible to get Σph exactly, and one must be content with

an approximate result.

If the approximate result is not a very good approximation, one should not try to solve

the problem in this fashion, in fact the evaluation cost for the irreducible diagrams in-

creases very quickly. Basically one should realize that Dyson’s equation is useful only in

weak coupling theory, where the perturbation is sufficiently weak that an adequate approx-

imation is obtained with a few terms in Σph(E). An important issue is the convergence of

Eqn. (3.36) which depends fundamentally on the shape of the self-energy. The presence of

3.7. Time-loop

-matrix: NEGF

poles in the Green’s function tell us that this argument is a very delicate one and that the

behavior of the Green’s function obtained via perturbation expansion is not easy. A way to

avoid bad divergences relies on the so called renormalization technique which is beyond

the scope of this work (see [51]).

One aspect of this result deserves a special mention. The self-energy has real and

imaginary parts which switches sign depending on whether E>

or E<

where

is the

chemical potential of the system. The same happens if we compute a self-energy for the

phonon propagators:

Im(Σph(E)) <0E>

Im(Σph(E)) >0E<

(3.37)

For the phonon propagator we can find another Dyson’s equation with a different self-

energy, Π(E), which is sometimes called polarization operator. This name is quite de-

scriptive, since the self-energy effects arise from the phonons causing polarization in the

medium. In our methodology the self-energy of the phonon is neglected.

The real and imaginary parts of the self-energies Σph(E)and Π(E)each have inter-

pretations. The imaginary part Im(Σph(E)) and Im(Π(E)) are interpreted as causing the

damping of the particle motion. They are related to the finite mean free path of the excita-

tion or to the energy uncertainty. The real parts are actual energy shifts of the excitation,

which may also change its dynamical motion. The excitation may alter its effective mass

or group velocity because of the self-energy contributions.

3.7 Time-loop S-matrix: NEGF

The formulation of the previous section has been done for the Equilibrium Green’s function

in which the system is supposed to remain in the equilibrium ground state at the beginning

and at the end of the process (that means, the system starts in the ground state for t=−∞

and returns in the same state for t= +∞).

Under this assumption, all the interaction processes included in the self-energy cannot

represent real transition of the system because it is forced to go back to the ground state,

they are usually called virtual processes. This formalism is all we need to compute the

Landauer-B¨

uttiker equation for transport. In the other cases we must move into the Non-

equilibrium formalism. Fortunately for us the machinery developed for the equilibrium

case remains practically intact. This means that the Green’s function machinery can be ap-

plied, after some generalizations, also for transport problems in which scattering processes

are present.

The S-matrix we have defined in Eqn. (3.15) depends on real time taken over the in-

terval (−∞,+∞). Now, assuming that also for a real interacting system we can assume the

adiabatic switching on of the interaction part of the Hamiltonian ˆ

V, the system is in its

ground state for t→−∞. The problem is that this is not true for t→+∞and the state for

Chapter 3. Non-Equilibrium Green’s Functions

that must be defined carefully.

Schwinger suggested another method of handling the asymptotic limit t→+∞. He

proposed that the time integral in the S-matrix should have two pieces: one going from

(−∞,

) while the second going from (

,−∞). Finally, we take the limit

→+∞. The

integration is a time loop which starts and ends at t=−∞. The new Green’s function

for non equilibrium system is defined over this new loop time variable and retains all the

properties and satisfies all the theorems of the equilibrium propagators discussed before

(see Fig. (3.4)).

Now we can define four different Green’s functions over the time loop axis, depending

on which branch the annihilation and creation operators are defined. If we manage also

the phonon operators we have then to handle about eight different Green’s functions. For-

tunately, there are some relations between these different propagators so that we need to

compute only two of them for the electrons and two for the phonons. In the case of equi-

librium Green’s function a further relation reduces the number of independent propagators

to one for the electrons and one for the phonons (see [61]).

Limiting to the case of electrons, the reason why there are four different propagators is

related to the S-matrix in the contour time:

S(−∞,−∞) = ˆ

TCtexp·−iZCt

dt1˜

V(t1)¸(3.38)

The integration path is the time-loop shown in Fig. (3.4). The variable t1goes (−∞,

)

and then (

,−∞). The operator ˆ

TCtorders along the entire loop, with earliest values t1in

C+

toccurring first. In expanding the S-matrix, we will encounter Green’s functions of the

form:

,t1−t2) = −ih|ˆ

TCt[˜c

(t1)˜c†

(t2)]|i (3.39)

Figure 3.4: The complex time contour for the NEGF formalism. The upper (time ordered) and

lower (anti-time ordered) branches are denoted with C+

tand C−

trespectively.

Depending where the two operators are positioned over the contour we can define the

four possibilities

1. Both on the upper branch C+

t: time ordered Green’s function Gt.

2. Both on the lower branch C−

t: anti-time ordered Green’s function G¯

3. The annihilation operator on the upper branch and the the creation operator on the

lower branch: greater Green’s function G>.

3.7. Time-loop

-matrix: NEGF

4. The annihilation operator on the lower branch and the the creation operator on the

upper branch: lesser Green’s function G<.

The analytic formulation of these four propagators is the following:

G>(t1,t2) = −ih|˜c

(t1)˜c†

(t2)|i (3.40)

G<(t1,t2) = ih|˜c

(t2)˜c†

(t1)|i (3.41)

Gt(t1,t2) =

(t1−t2)G>+

(t2−t1)G<(3.42)

G¯

t(t1,t2) =

(t2−t1)G>+

(t1−t2)G<(3.43)

Gr(t1,t2) = Gt−G<=G>−G¯

t(3.44)

Ga(t1,t2) = Gt−G>=G<−G¯

t(3.45)

where we have added two more propagators, the retarded (Gr) and advanced (Ga), obtained

as a linear combination of the previous four. They are particular important for their ana-

lytical properties, the retarded propagator has poles only in the lower half of the complex

plane, whilst the advanced in the upper part. Moreover, it can be easily demonstrate two

fundamental identities between the lesser, greater, retarded and advanced propagators:

G>(E)−G<(E) = Gr(E)−Ga(E)(3.46)

Ga(E)=(Gr(E))†.(3.47)

Just to give an example we write the unperturbed propagators for a free electron in time

domain and, after a Fourier transform, in energy domain:

;t,0) = −i(1−nF)e−i

ελ

t→−2

i(1−nF)

(E−

ελ

)(3.48)

;t,0) = inFe−i

ελ

t→2

inF

(E−

ελ

)(3.49)

;t,0) = −i[

(t)−nF]e−i

ελ

t→1

E−

ελ

δλ

(3.50)

G¯

tf(

;t,0) = −i[

(−t)−nF]e−i

ελ

t→−1

E−

ελ

−i

δλ

(3.51)

;t,0) = −i

(t)e−i

ελ

t→1

E−

ελ

(3.52)

;t,0) = i

(−t)e−i

ελ

t→1

E−

ελ

−i

(3.53)

where

and

δλ

are infinitesimal quantity,

is always positive while

δλ

is negative for

occupied states and positive for unoccupied states.

The same rules apply for phonon propagators:

Chapter 3. Non-Equilibrium Green’s Functions

(t1,t2) = −ih| ˜

(t1)˜

(t2)|i (3.54)

(t1,t2) = −ih| ˜

(t2)˜

(t1)|i (3.55)

(t1,t2) =

(t1−t2)D>

(t2−t1)D<

(3.56)

D¯

(t1,t2) =

(t2−t1)D>

(t1−t2)D<

(3.57)

(t1,t2) = Dt

−D<

(t1−t2)[D>

−D<

](3.58)

(t1,t2) = Dt

−D>

=−

(t2−t1)[D>

−D<

](3.59)

For the free phonon propagators we get,

;t,0) = −i[(N

+1)e−i

ωα

t+N

ωα

→ −2

i[(N

+1)

(E−

ωα

)+N

αδ

(E+

ωα

)] (3.60)

;t,0) = −i[(N

+1)ei

ωα

t+N

e−i

ωα

→ −2

i[(N

+1)

(E+

ωα

)+N

αδ

(E−

ωα

)] (3.61)

;t,0) = −2

(t)sin(

ωα

t)→2

ωα

E2−

(3.62)

;t,0) = −2

(−t)sin(

ωα

t)→2

ωα

E2−

−i

(3.63)

The phonon Green’s functions in equilibrium are expressed in terms of the phonon occu-

pation number N

, which equals to a Bose-Einstein distribution in thermal equilibrium at

finite temperature.

3.7.1 Dyson’s equations for NEGF

Keldysh was the first to reformulate the NEGF formalism in a matrix notation which pre-

serves the intuitive picture of the loop time contour, but, at the some time, results to be a

better tool to handle the many propagators. The advantage of the matrix notation is that we

can reformulate all the propagators in one matrix 2×2 and we can extract all the properties

of the different propagators just by the usual matrix algebra:

G=µGt−G<

G>−G¯

t¶(3.64)

where Gsatisfies the Dyson’s equation

G=G0+G0ΣphG.(3.65)

3.7. Time-loop

-matrix: NEGF

This matrix notation, as already suggested by Eqn. (3.65), can be extended to all the other

quantities, like self-energies and so on. Applying a rotation in the Keldysh space, it is

possible to get another particularly useful matrix representation in which, instead of the

time ordered and anti-time ordered propagators, the retarded and advanced propagators

appear:

G0=µGrG>+G<

0Ga¶.(3.66)

Using the Dyson’s equation in Eqn. (3.65) we get the Dyson’s equations for the retarded,

advanced and lesser and greater quantities:

Gr,a=Gr,a

0+Gr,a

0Σr,aGr,a(3.67)

G<,> = [1+GrΣr]G<,>

0[1+GaΣa]+Gr

0Σ<,>

ph Ga(3.68)

Iterating the matrix representation we can obtain the formula for operators Xwhich de-

pends on two or more operators (X=AB or X=ABC). A summary of the results, using

the so called Langreth rules [61], is presented in Table (3.1).

Contour time Real time

X=RCtAB X<=Rt[ArB<+A<Ba]

Xr=RtArBr

Y=RCtABC Y<=Rt[ArBrC<+ArB<Ca+A<BaCa]

Yr=Rt[ArBrCr]

X(t,t0) = A(t,t0)B(t,t0)X<(t,t0) = A<(t,t0)B<(t,t0)

Xr(t,t0) = A<(t,t0)Br(t,t0) + Ar(t,t0)B<(t,t0)+Ar(t,t0)Br(t,t0)

Y(t,t0) = A(t,t0)B(t0,t)Y<(t,t0) = A<(t,t0)B>(t0,t)

Yr(t,t0) = A<(t,t0)Ba(t,t0) + Ar(t,t0)B<(t0,t)

Table 3.1: Langreth rules for analytical continuation of propagators and self-energies from the

contour time back to real time. The rules for advanced and greater quantities can be obtained

starting from the one for retarded and lesser, just exchanging all Frwith Faand all F<with F>

and viceversa.

Chapter 4

Density Functional Based Tight-Binding

In Chapter 3was shown that a full Non-Equilibrium Green’s Function (NEGF) formalism

can be developed via a perturbation approach starting from zero order building blocks. The

zero order Green’s functions, both for phonons and electrons, are supposed to be exactly

solvable. This assumption comes from the fact that the Hamiltonian ˆ

Hcan be partitioned

into two parts, ˆ

H0and ˆ

V.ˆ

H0is a reasonable approximation of the full Hamiltonian that

can be solved exactly and that contain both a single particle Hamiltonian ˆ

0and two other

terms for the open boundary conditions.

In many NEGF calculations for molecular transport simulations the starting point is

a Density Functional Theory (DFT) calculation. In fact, despite some important draw-

backs of the method, DFT is still the best compromise between complexity and accuracy

for molecular electronic simulations. The main advantage of DFT compared to other ab-

initio methods, is the possibility to include a large number of atoms from the contacts for

the extended molecule regions and to handle the open boundary conditions. Moreover,

the electron-electron correlation can be included immediately via an exchange-correlation

functional without any other effort at the Green’s function level. However, it is well known

that the correlation included in DFT, especially at the lowest level of approximation, the so

called Local Density Approximation (LDA), gives a poor description of the electronic mu-

tual interaction. This leads, in many cases of interest, to a weak evaluation of the electronic

spectrum of the molecule and to quantitative and sometimes qualitative errors in the current

calculated, compared to experiments. One dramatic example is represented by the under-

estimation of the energy gap between the Highest Occupied Molecular Orbital (HOMO)

and the Lowest Unoccupied Molecular Orbital (LUMO) which produces an overestimation

of the magnitude of the current (sometimes of a couple of orders of magnitude). At the end

of the present Chapter these problems will be addressed in a more specific way.

The present work is based on a code the development of which started in the early

1990, a DFT based Tight-Binding called DFTB. Recently, a new version of the code used

in this work, not specifically for transport, has been released under the name DFTB+ (see at

the website http://www.dftb-plus.info). However, several other implementations of DFTB

exist. The description of DFTB is briefly treated in this Chapter, nevertheless, we refer the

Chapter 4. Density Functional Based Tight-Binding

interested reader to the wide literature about the method of DFTB [62,63,64] and more in

general DFT [65,66] for a more comprehensive documentation.

4.1 DFT and Kohn-Sham formulation

Compared to other methods, DFT recasts the problem of calculating the electronic structure

of a system in terms of the electronic density. The spin-free density in terms of the wave

function can be formally defined as follows:

n(r) = NZΨ∗(r,r2,...,rN)Ψ(r,r2,...,rN)dr2...drN(4.1)

where Ψ({ri})is the multi-particle wave function. The concept of employing the electron

density in an energy functional dates back to Thomas-Fermi theory [67,68]. The very

foundation of DFT are two rigorous statements. The first is a lemma linking the density

and the potential of a Hamiltonian. The second is a theorem about the existence of a uni-

versal functional of the electron density, universal in the sense that does not depend on

the external potential. This functional can be minimized, with a contribution which comes

from the external potential, in a variational sense. The minimum is achieved for the ground

states density and gives the ground state energy.

Lemma (Basic Lemma of Hohenberg-Kohn) The ground state density n(r)of a bound

system of a fixed number of interacting electrons in an external potential Vext(r)determines

this potential uniquely (to within an additive constant).

The very first step for the study of molecules and solids is the calculation of the total

energy of a given arrangement of atoms; eventually, we want to find, within suitable bound-

ary conditions, the atomic arrangement that is lowest in total energy and hence the most

stable. Besides contributions from the internuclear interactions, this task requires knowl-

edge of the energy of the electronic system. The energy of the system is related to the wave

function and so to the Schr¨

odinger equation. The lowest eigenvalue represents the ground

state of the system. For a physical system the problem can be rearranged starting from the

density of electrons. From the Lemma we know that the density fixes the external potential

and, therefore, the Hamiltonian. On the other hand the Hamiltonian defines completely

the spectrum of wave functions. This means that there is a fundamental connection be-

tween the electron density and the ground state of the system, Ψ0({ri}). Hence, the latter

is strictly a functional of n(r), albeit not a trivial one. This means that the energy of the

system can be expressed as a pure functional of the density:

E[n] = ZVext(r)n(r)dr+F[n](4.2)

The functional is split in two parts: the first term depends only on the external potential and

the density, the second contains all the information about the kinetic and Coulomb inter-

action between the electrons. The functional F[n]is universal in the sense that it does not

4.1. DFT and Kohn-Sham formulation

refer to an external potential. A minimization step, in the space of all densities for which

the functional F[n]is defined, is required to find the density nGS of the system and so the

ground state energy. This variational principle is supported by the following theorem.

Theorem (Hohenberg-Kohn variational principle) For all densities the functional E[n]

has a minimum, which it assumes at the ground state density nGS.

min(E[n]) = E[nGS](4.3)

The electronic ground state energy of a system of interacting electrons is thus expressed

in terms of their spatial charge density instead of the many-particle wave function. The

main problem of DFT is to find the shape of the universal functional F[n]. The statements

on which DFT is based demonstrate only the existence of this functional, but they do not

suggest how we can build it.

4.1.1 The Kohn-Sham equations

In 1965 Kohn and Sham [69] proposed a new reformulation of the energy functional,

Eqn. (4.2), in order to find an analytical solution for it. The functional F[n]is split into

parts:

F[n] = T0[n]+EH[n]+Exc[n](4.4)

There are three components. The first term, T0[n], is the kinetic energy functional for

a fictitious system of non-interacting electrons producing the same density as n(r). The

second term, EH[n], is the so called Hartree energy, arising classically from the mutual

Coulomb repulsion of all electrons:

EH[n] = 1

2Z Z n(r)n(r0)

|r−r0|drdr0(4.5)

The last term, Exc[n], called the exchange-correlation functional, is a correction term, which

accounts for all many-body effects in F[n]. The treatment of this term decides upon the vi-

ability of any DFT implementation. All terms but the last of Eqn. (4.4) are known analyt-

ically. The functional form of the correction term Exc[n]is unknown and must be approx-

imated. Practical applications of the DFT are classified according to the approximations

taken for the exchange-correlation functional Exc[n].

The first attempt to give a shape to the exchange-correlation functional has been the

LDA,

Exc[n] = Z

xc[n(r)]n(r)dr(4.6)

where

xc[n]is the exchange-correlation energy per electron in a uniform electron gas of

density n(r) = const. The energy

xc[n]is a function of only the local density value, and

no longer a functional of the global density distribution and its form is well known.

Chapter 4. Density Functional Based Tight-Binding

The T0[n]term can be described defining a set of single particles orbitals which are

orthogonal and normalized and that are related to the real density of the system as follows:

n(r) =

∞

∑

i=1

ni|Φi(r)|2(4.7)

hΦi|Φii=Z|Φi(r)|2dr=1 (4.8)

where niare the occupation numbers of these orbitals. The occupation numbers are con-

strained to be a certain set of values: integer numbers (0 and 1 for spin-unrestricted solution

and 0,1 and 2 for spin-restricted solutions) or real numbers following a Fermi-Dirac distri-

bution in case of small thermal excitation. The kinetic term becomes:

T0[{Φi(r)}] =

∞

∑

i=1

niZΦ∗

i(r)µ−4

2¶Φi(r)dr(4.9)

The entire functional, E[n], can be expressed as a functional of the new orbitals. The

minimum of the energy functional with respect to the density becomes a variation respect

the orbitals, with the constraint that the orbitals must remain orthogonal and normalized:

Φ∗

i(r)(E[n(r)]+

∞

∑

i=1

i·1−Z|Φi(r)|2dr¸)=0 (4.10)

where the

iare the Lagrange parameters. Solving Eqn. 4.10 we get at the end a Schr¨

odinger-

like equation,

·−4

2+Vext(r)+VH([n(r)],r)+Vxc([n(r)],r)¸Φi(r) =

iΦi(r),(4.11)

this equation plus Eqn. (4.7) are called the Kohn-Sham equations. The three terms Vext,VH

and Vxc are the variational counterpart of the energy functional. Their sum is a potentials

which describe the external field, the classical electron-electron repulsion and the many-

body effects (exchange due to Pauli exclusion principle and correlation) respectively. The

picture of Kohn-Sham equations is a set of Nnot interacting particles merged in a effec-

tive external potential which is the sum of these three terms. We can be tempted to give a

physical interpretation to the eigenvalues (Lagrange parameters) of Eqn. (4.11)

i. Janak

has demonstrated [70] that, in fact, the eigenvalue of the highest occupied orbital is fun-

damentally the first ionization energy changed in sign. The interpretation of the others

eigenvalues is more controversial.

An analytic solution of the Kohn-Sham equations is rather difficult. Instead, the un-

known functions Φi(r)are sought numerically employing tools from mathematical physics.

Two classes of approaches will yield a solution: real-space grid methods and basis func-

tion expansions. Either approach will transform the differential equation for the unknown

functions into a set of algebraic ones for unknown expansion coefficients of a reference

4.2. DFTB: method and approximations

basis of functions. The general expansion of the one-electron wave functions reads:

Φi(r) =

∑

νφν

(r)(4.12)

Once the basis is chosen, depending to the kind of problem, the coefficients Ci

are obtained

solving a general secular equation:

∑

¡h

µν

−

µν

¢=0 (4.13)

where

µν

φµ

|ˆ

φν

i=Z

∗

·−4

2+VH+Vxc([n(r)],r)¸

φν

dr(4.14)

µν

φµ

φν

i=Z

∗

µφν

dr(4.15)

Since the orbital solution depends on the density, and in turn, by Eqn. (4.12), on the orbitals

themselves, Eqn. (4.14) needs to be solved self-consistently starting from a guess density,

n0(r).

4.2 DFTB: method and approximations

The Kohn-Sham equation are implemented in the DFTB method with some approxima-

tions. These mostly take the form of separations of global quantities like potentials and

wave functions into atomic contributions. This separation principle allows the effective

treatment of compound atomic systems by building as much as possible upon preparatory

work, which is performed beforehand on isolated subsystems of atoms and atom pairs.

This previous work is used to evaluate the computational expensive reference Hamiltonian

and other contributions, like the nuclei repulsion, into simple tables.

The total-energy expression of the DFTB approach follows the one that can be derived

from full Kohn-Sham equations,

Etot[n0+∆n] =

occ

∑

i=1

nihΦi|ˆ

h0|Φii+Erep[n0]+E2[n0,∆n]

=EBS +Erep +E2(4.16)

where EBS is the trace of the eigenstates of the system over a Hamiltonian ˆ

h0which depends

on the guess density n0only:

h0=−4

2+Ve f f [n0,r].(4.17)

Chapter 4. Density Functional Based Tight-Binding

The second term, Erep constitutes a repulsion energy similar to standard tight-binding the-

ory, which subsumes the double counting terms of the reference Hamiltonian, as well as

the nuclear repulsion. Finally, the last term, E2, describes atomic charge fluctuations and

is subject to a self-consistent treatment.

The main approximations to obtain the DFTB method from the stationary principle of

DFT are in turn:

•superposition of pseudo-atomic densities as starting density;

•usually a minimal-basis, valence-only Linear Combination of Molecular Orbital

(LCAO) wave functions, is used;

•two-centre Hamiltonian (neglect of crystal-field and three-centre terms);

•repulsive pair potential for the double-counting and inter-nuclear energies;

•for the second-order correction:

1. monopole approximation and extrapolation of

Vxc[n]/

2. Mulliken charges;

3. Integral approximations for the γmatrix;

4.2.1 Pseudo-atomic starting density

In DFTB, the starting density is chosen as superposition of slightly compressed densities

of neutral atoms,

n0(r) = ∑

0(ra);ra=r−Ra(4.18)

The densities of free atoms are too diffuse to be a good initial guess in compound systems.

The densities are the result of a self-consistent LDA or General Gradient Approximation

(GGA) calculation of pseudo-atoms, i.e., atoms placed within a weak parabolic constriction

potential as expressed in a modified Kohn-Sham equation:

·−4

2+Vat

e f f [na

0]+µr

r0¶m¸

psat

(r) =

psat

νψ

psat

(r)(4.19)

The constriction potential is characterized by its exponent mand range r0. The exponent

was shown to have rather small influence on the final results, so that m=2 is usually

used. For the range parameter, a number of calculations has lead to results with optimal

transferability in covalent systems (barring 3dtransition metals) using r0≈1.85rcov, where

rcov is the covalent radius of the given element. The pseudo atomic wavefunctions are used

for the LCAO basis. These wavefunctions are represented by Slater-type orbitals (STO)

characterized by coefficients ai j and exponents

ψν

(r)≡

lmn(r) = ∑

i=1∑

j=0

ai jrl+je−

irYlm ³r

r´.(4.20)

4.2. DFTB: method and approximations

4.2.2 Tight-binding integrals and the two-centre approximation

The pseudo-atomic are used as basis set for the expansion of the eigenfunctions

Φi(r) = ∑

a∑

[a]

νψν

(ra).(4.21)

The reference Hamiltonian can be represented on this basis, giving matrix elements de-

noted h0

µν

and non-orthogonal overlap elements S

µν

ψµ

|ˆ

h0|

ψν

i(4.22)

µν

ψµ

ψν

i=Z

∗

µψν

dr.(4.23)

These integrals are calculated immediately and are tabulated as function of distance be-

tween the two centres. The integrals at general difference vectors as needed in the system,

∆Rab =Ra−Rb, are transformed using projection relations [71]. Due to symmetry, only

few integrals between basis functions remain nonzero for angular momenta up to l=2.

In order to achieve the two-centre representation for the Hamiltonian matrix elements,

the effective Kohn-Sham potential is formally decomposed into atomic-like contributions.

Due to the non-linearity of the exchange-correlation potential, there are two ways to do this

in practice, either to sum atomic potentials or atomic densities:

Ve f f [n0,r]≈½∑cV0

c[n0

c,rc]potential superposition

Ve f f [∑cn0

c]density superposition (4.24)

Depending on the centres involved for the basis functions and the potential, the Hamilto-

nian matrix elements fall into a number of categories, which are summarized in Table (9.1).

Type Classification Centres Status

(A) onsite-terms a=b=cretained

(B) crystal-field terms a=b6=cneglected

(D) three-centre terms a6=b6=cand b6=cneglected

Table 4.1: Integral types in the DFTB Hamiltonian h0

µν

. Centres a and b denote orbital centres,

with basis functions

∈a,

∈b and c is the potential centre.

The neglect of three-centre terms provides the largest practical simplification because

on the one hand the handling for three centres is more involved than for two centres and

on the other hand there are many more combinations. On the other hand, the crystal-

field integrals are relatively simple. However, the integral neglections work only in concert

because there is error-cancelation of a considerable degree. The justification of this process

is complex and was discussed originally by Seifert et al. [72] and later reviewed in ref. [63].

Chapter 4. Density Functional Based Tight-Binding

Moreover, using a valence-only basis assures the core-valence orthogonalisation between

different centres.

Among the retained integrals, type (A) represents on-site energies

εν

of single atoms,

νν

εν

(4.25)

which are obtained in the first step of the atom calculation. In the other remaining integral

type (C), the potentials and densities of two distinct atoms are to be combined:

µν

ψµ

|−4

2+½V0

a[ra]+V0

b[rb]

Ve f f [n0

a+n0

b]¾|

ψν

i;a[

],b[

].(4.26)

4.2.3 Repulsive potential

The next term that must be evaluated is the so called repulsive term which includes the

double-counting terms between EBS, evaluated at the input density, and E2, the self-con-

sistent correction, in the energy functional. Moreover, it includes also the inter-nuclear

repulsion

Enuc =1

Nat

∑

a,b6=a

ZaZb

|Ra−Rb|.(4.27)

The DFTB Method sum of short-ranged repulsive pair potentials:

Erep[n0,{Ra}]≈1

2∑

a∑

a6=b

Vab

rep(|Ra−Rb|).(4.28)

This approximation is justified by two observations:

•With n0represented by neutral atomic fragments, there are no long-range Coulomb

interactions in the combined electrostatic Hartree and nuclear contributions to the

double counting energy, EH[∑ana

0] + Enuc(Ra), due to mutual screening. Further-

more, since the atomic starting densities are spherically symmetric, the Hartree inte-

grals for atom pairs (4.5) depend on internuclear distance only. The same is trivially

the case for the nuclear repulsion. Therefore, these contributions together are repre-

sentable by short-range pair potentials without loss of accuracy.

•Contributions due to exchange and correlation are not separable into pair potential

form per se because of the non-linearity of the xc-functional. However, a cluster ex-

pansion [63] allows to extract two-body components. Its higher order terms involve

the overlap of the densities of three centres, which is negligible for the compressed

starting densities used here. The remaining two-body terms may again be represented

by pair-potentials.

All monomer contributions are contained within the atomic orbital energies

εν

to ensure

that the repulsive potential actually goes to zero in the dissociation limit. The repulsive

4.2. DFTB: method and approximations

potential is obtained from self-contained ab initio calculations of the energy of a set of

reference molecules for a range of a typical bond length. This is the most delicate term in

the DFTB method, because still there is not a systematic way to get a set of molecules and

compounds that can assure the transferability of the values obtained. Most conveniently,

the reference molecules are dimers, but also methane-like structures. For solid state calcu-

lations, as in gold crystals which are widely used in this work, crystal reference structures

may be used. For each reference structure j, the energy difference,

rep(r) = Ej

sc f (r)−Ej

BS(r)(4.29)

is calculated as function of distance r. As a consequence of the concepts and approxima-

tions taken thus far, the set of repulsive energies has the following important properties:

•Usually ˜

Erep >0. For small distances, i.e., close atoms, Erep has a steep repulsive

slope indicative of strong Pauli repulsion of the electron shells and the ion-ion repul-

sion.

•It decays rapidly to zero between typical first- and second-neighbor distances. This is

another indication that the pair-potential representation of the double-counting terms

embodied within Erep, is valid.

•Most importantly, ˜

rep for different molecules are close to each other. This property

is indicative of the degree of transferability of the method, i.e., its applicability to a

wide range of atomic structures.

It was found that ˜

Erep(r)for high-coordinated crystal phases such as fcc and bcc often does

not coincide with the curve of other structures. Such phases cannot be treated with the

present method. All the crystalline gold electrodes used in the simulations in the last four

Chapters (6-9) were in fact not allowed to relax using DFTB. For practical calculations,

the set of repulsive energies are represented numerically by splines for each atom type

combination.

4.2.4 Second-order correction

The last term in the DFTB total-energy expression to discuss is the second-order correc-

tion E2[n0,∆n]. This term becomes important in the simulation of heteroatomic molecules

and polar semiconductors where chemical bonding is influenced considerably by charge

transfer effects and long-range Coulomb interactions.

In line with previous procedures, the charge fluctuations ∆nare decomposed into atomic

contributions which are expected to decay rapidly with increasing distance from their cen-

tre. The second-order term then reads,

E2[n0,∆n] = 1

2∑

a,bZ Z ·1

|r−r0|+

Vxc([n(r),r])

n(r0)¸n0

∆na∆nbdrdr0(4.30)

Chapter 4. Density Functional Based Tight-Binding

The term is expected to give a small contribution of the total energy, so that integrations

can be approximated. First, one expresses ∆naas a mere monopole contribution:

∆na(r) = ∑

lmFa

lm(ra)Ylm µra

ra¶≈∆qaFa

00(ra)Y00;ra=r−Ra.(4.31)

Thus, Eqn. (4.30) takes a simple matrix form:

E2[n0,∆n]≈1

∑

a,b

ab∆qa∆qb(4.32)

where γmatrix is:

ab =1

Z Z ·1

|r−r0|+

Vxc[n]

n(r0)¸n0

00(ra)Fb

00(r0

b)drdr0.(4.33)

The γmatrix contains the dependency of E2on ionic positions. This matrix, except for

the Vxc part, is well known from the CNDO-formalism by Pople, Santry and Segal [73].

The Fi

00 term is the normalized radial dependence of the density fluctuation on atom i. If

the latter are assumed to have a fixed form, spherically symmetric, to be weighted by ∆qa,

the only geometry parameter in

ab will be the interatomic distance R=|Ra−Rb|. The

limit R→0 means coinciding atomic centres aand b. In this case,

aa equals a Hubbard-

like parameter Uafor the atom. The Hubbard parameter is related to the chemical hard-

ness

≈Ua/2, which is a measure of the ionization potential and electron affinity of the

atom. Neglecting the influence of the environment and employing Janak’s theorem [70],

Uais obtained non-empirically at the DFT level during the pseudo-atom calculation as the

derivative of the HOMO of the free atom with respect to its occupation number:

aa =Ua=

∂

2Eat

∂

at ¯¯¯¯q=q0

∂ε

HOMO

∂

nHOMO

(4.34)

In the limit of large interatomic distances,

ab reduces to a 1/Rdependency, since Vxc

interactions vanish in this case. To obtain a continuous transition between the limits of

small and large interatomic distances an interpolation formula is used [74]. Assuming a

normalized Slater-type function with a range parameter

00(r) =

e−

a|r−Ra|(4.35)

the Coulomb integral of two such charge densities reads, after some manipulations (for

integrals of this type, see Pople and Beveridge [75]):

ab(

b,R) = 1

R−S(

b,R)(4.36)

where

4.2. DFTB: method and approximations

b,R) = e−

aRÃ

b−

a)2−

b−3

b−

a)3!+

+e−

bRµ

a−

b)2−

a−3

a−

b)3¶;

a6=

b(4.37)

,R) = e−

Rµ48+33

R+9

R2+

48R¶;

(4.38)

In the limit of short distances

ab should equal the Hubbard parameter for an isolated atom.

This being a known quantity, one obtains a one-to-one relation between the hitherto un-

specified charge fluctuation range

aof an atom and its Hubbard parameter:

aa =Ua=5

a(4.39)

Since Eqn. (4.36) was derived from pure Coulomb interactions, yet the Hubbard param-

eter incorporates xc contributions, this step may seem inconsistent. To reconcile the ap-

proaches, the following should be considered: appreciable deviations from the point charge

behavior of 1/Roccur only close to atoms, namely, within radii of typical bond lengths.

Adjusting the inner limit to include xc interactions by way of equating it to the ab-initio

DFT value, one corrects mainly the on-site terms of

ab.

To finally evaluate the second-order energy contributions, one needs the atomic charge

deviations ∆qa. These are obtained from Mulliken charges. Given real-valued eigenvectors

in an atomic basis with overlap matrix S, the Mulliken charges qaand charge deviations

∆qaon atoms are:

qa=

occ

∑

i=1

niqi

occ

∑

i=1

ni∑

[a]

∑

µν

=∑

[a]

∑

µν

∆qa=qa−q0

a(4.40)

where q0

aare the charges of the respective neutral atoms.

4.2.5 DFTB secular equation

With the approximations discussed above the DFTB total-energy reads:

Etot =

occ

∑

i=1

ni∑

∑

µν

Nat

∑

ab∆qa∆qb+Erep[Ra](4.41)

Chapter 4. Density Functional Based Tight-Binding

Given a set of atomic coordinates {Ra}and the resulting matrices h0,Sand γ, the LCAO

coefficients which minimize the DFTB total energy (4.41) are found by the variation prin-

ciple subject to orbital normalization:

∂

"Etot +

occ

∑

i=1

ni˜

iÃ1−∑

∑

µν

!#=0 (4.42)

The procedure is quite similar to the standard tight binding case. The resulting secular

equation reads:

∑

µν

−˜

µν

) = 0 (4.43)

where, as a consequence of the minimization of E2alongside EBS the original Hamiltonian

matrix elements are augmented as follows:

µν

=h0

µν

∑

(

ac +

bc)(qc−q0

∈a;

∈b; (4.44)

As for the generic Kohn-Sham system, Eqn. (4.43) is a generalized eigenvalue problem

which is solved using standard libraries. Since the charges in the augmented Hamiltonian

depend upon the coefficients of the solution, the equation has to be solved self-consistently.

The iteration is driven by some mixing algorithm for the charges {qa}.

Forces on atoms can be obtained analytically without resorting to repeated energy cal-

culations for finite displacements. The essential result of the DFTB calculation is available

in form of the total energy and the electronic eigenvectors. Other physical observables

like charge distributions or vibrational properties may be extracted, the latter not without

considerable numerical effort.

4.3 Disadvantages of DFT in transport simulations

There is a number of significant limitations known to be associated with modern density

functionals. These include [76]: poor treatment of dispersive forces, poor treatment of

covalent bond breakage, poor treatment of systems involving partial electron removal, and

poor treatment of conjugated

systems.

Firstly, dispersive forces are very important for the determination of the geometrical

structure of interfaces, being the primary contributor to physisorptive processes including

the strong

-stacking interactions that can occur between gold and aromatic molecules, the

interactions of nitrogen bases with gold, and contributions along with covalent bonding to

thiol-gold interactions. While these forces are treated sporadically by DFT, they are also

excluded from the basic DFTB hamiltonian. However, DFTB+ now includes provision for

the inclusion of these forces.

Secondly, in general, covalent bond breakage is treated poorly by modern DFT func-

4.3. Disadvantages of DFT in transport simulations

tionals, with their poor description of long-range electron correlation leading to contami-

nation of supposedly radical-like reaction products with ionic-like structures [77,78]. This

failure has immediate consequences for the determination of realistic strained electrode-

molecule junction structures as dynamic bond breaking and forming processes are likely

to be important. However, it has more profound effects owing to the way in which Green’s

function-based methods perceive the current carrying process [79]. Effectively, they per-

ceive conduction as occurring between the two macroscopic contacts through tunneling.

All that is important is the tunneling probability, effectively obtained from the calculated

energy gap between the electrode-localized tunneling orbitals. The process by which this

energy gap arises is irrelevant: it could be via direct through-space interactions of the elec-

trodes or via superexchange [80] through molecule-assisted pathways. Viewed in this way,

through-molecule conduction is analogous to through-space conduction between two elec-

trodes with broken inter-electrode covalent bonds, allowing errors associated with covalent

bond breakage to enter calculations in a profound way. Covalent bond breakage is an open-

shell problem for which valid solutions may be obtained using spin-unrestricted means if

just one bond is involved [81], a situation not typical in molecular electronics applications.

Inappropriate application of the Kohn-Sham theorem leads to error up to 20% in the

underestimation of the magnitude of excited state energies, energies interpreted by Green’s

function codes as depicting the strength of tunneling, leading to large overestimation of

the conductivity. Also, as the LDA functional in DFT provides realistic treatments of the

metal electrodes but gives a poor description of the molecular component. Only hybrid

functionals, however, include the long-range exchange interactions that are often critical

in electron transfer problems. Hence the development of improved density functionals is a

priority.

The third limitation of DFT is related to the so called self-interaction energy error.

This error arises from the approximated form of the exchange potential in DFT. In exact

exchange potential a charge in a certain state cannot interact with itself because the term

of self-interaction in the Hartree energy cancels exactly with the same term in exchange

potential. In the case of DFT this cancelation is not exact and this means that a certain self-

interaction survives. This is particularly problematic for localized states [66]. Practically,

this results in low energies for charge transfer states, a type of state that can contribute to

through-molecule conductivity, and an incorrectly positioned molecular highest-occupied

molecular energy (HOMO) level that is typically too high in energy by 3-4 eV. In molec-

ular electronics applications, this leads to errors in the line-up of the orbital bands of the

electrodes and the molecule [82].

Lastly, we consider limitations of standard DFT in treating conjugated

systems.

Many molecules used in molecular electronics applications, especially the highly conduc-

tive ones, such as 1,4-benzenedithiol, are of this type. Molecular conduction involves

partial oxidation or reduction processes of the molecule, processes whose energies are

identified in single-particle implementations of Green’s function kinetics theories with

molecular-orbital energy differences. While such an identification follows naturally from

Hartree-Fock based approaches through application of Koopmans’theorem, in DFT the

calculated orbital-energy differences actually reflect the much smaller energies associated

with optical transitions [82]. Also DFTB is subject to errors of this type [83]. For standard

Chapter 4. Density Functional Based Tight-Binding

DFT, the interplay between the self-interaction correction (HOMO error) and the unde-

sired physical interpretation of the band gap is complex. Fortunately, the two errors cancel

to some extent for the band line-up of the molecular lowest-energy unoccupied molecular

orbital (LUMO).

However, as the self-interaction error also results in significant underestimation of the

actual excited-state optical transition energies, both effects reinforce each other in making

molecules appear metal-like. Polyacetylene, for example, is predicted by modern density

functionals to have a ground state of at least triplet spin multiplicity rather than being a

wide band-gap semiconductor [84], with associated errors in calculated polarizabilities for

the actual closed-shell ground state [85]. Porphyrin and chlorophyll molecules, ubiquitous

for their roles in naturally occurring systems involving through-molecule conduction, are

also poorly described [84].

Chapter 5

The Electron-Phonon Code

The present Chapter is focused on the implementation of the Green’s function code for the

electron-phonon interaction based on DFTB. The code was implemented by the author of

the present work and Dr. Alessandro Pecchia, of the Univesrity of Rome “Tor Vergata”.

It belongs to the so called gDFTB family, see [86], a set of transport codes applying equi-

librium and non equilibrium Green’s function formalism and based on DFTB. The codes

already developed, beside the electron-phonon, are a gDFTB version for Landauer-B¨

uttiker

evaluation of the elastic current and another for a full treatment of the bias in the device

with the solution of the Poisson equation in NEGF. The first code is already contained in

both the electron-phonon and bias codes as a special case when the interaction, in one case,

or the Poisson solving, in the other case, are neglected. A further step will be the merging

of the bias and electron-phonon codes in a final program. The aim of this final gDFTB

code is to be merged in a full multiscale approach to transport called Tibercad. This final

package will contain modules to simulate transport in electronic devices at different scales,

from micrometers down to nanoscale, and interlace them. The gDFTB is the atomistic part.

The implementation of an electron-phonon calculation requires a certain number of

steps. First of all the definition of the system simulated, that is the definition of the par-

titioning in device and contact regions. The second step is the relaxation of the geometry

between the two electrodes and the calculation of the zero order DFTB Hamiltonian ˆ

The third step is the calculation of the normal mode of vibrations and the evaluation of

the electron-phonon coupling matrices elements M

µν

, where the index

runs over all the

normal modes and the indices

and

over the local basis set. This last step is very deli-

cate because the electron-phonon interaction depends strongly on the couplings. To obtain

a good set of couplings we need a fully relaxed geometry, the relaxation of the structure

must in fact reach a global minimum in order to really catch a good description of the

vibrational spectrum.

Finally the calculation of the current can be performed. The zero order Hamiltonian

is used to compute the retarded self-energies needed by the code, for the open boundary

conditions (Σr

Land Σr

R) and for the electron-phonon interaction (Σr

ph). Indeed, in the NEGF,

we need to compute also the equivalent quantities for the lesser and greater quantities. All

Chapter 5. The Electron-Phonon Code

the matrices can be later on used to compute the current, via the Meir-Wingreen equation,

and the power dissipated in the device.

The Chapter is organized as follows: first, a brief description of the main approxi-

mations used in the code is given. Then, the treatment of the open boundary conditions

follows. The next sections are devoted to describe in more details the calculation of the

electron-phonon self-energies and the vibronic couplings. Finally, the last part of the Chap-

ter describes in details the flowchart of the code.

5.1 Approximations in the electron-phonon code

This short paragraph is focused on the approximations implemented in the electron-phonon

code. A full description of the code and its potentiality cannot dispense to address also

the main simplifications introduced in the calculations of the quantities needed for the

evaluation of the current.

The set of approximations related to the DFTB Hamiltonian have been already de-

scribed in Chapter 4and they will be not repeated here. Further approximations appear at

the NEGF level. First of all the calculation of the electron-phonon self-energy is performed

at the first level of perturbation, called Born approximation (BA). The diagrammatic rep-

resentation of that term is shown in Fig. (5.2, a).

A better description of interaction can be performed allowing a self-consistent Born ap-

proximation (SCBA) in which the propagators of the self-energy is renormalized Fig. (5.2,

b). In this way we are including other diagrams in the self-energy. However, the difference

in the results between SCBA and BA is usually pretty small, within few percentage. The

real advantage of SCBA is that it assures the conservation of the charges in the current.

This conservation is in fact slightly broken in the mere BA and this can become a problem

in case of high current amount, when for example one of the states of the device is in reso-

nance. However, also an SCBA is not able to describe scattering processes if we are in the

strong electron-phonon couple regime. The entire perturbation approach to the self-energy,

in fact, fails badly in the latter case.

Problems can arise also in the weak coupling regime if the device region is weakly

bounded to the contacts (weak contact couplings). In that case it is possible that the electron

lifetime in the device is long enough for the formation of complex deformations of the

device itself, like polaron states. Also in that case the perturbation approach is not able

to recover the correct quantities. Fortunately, in many molecular electronics experimental

setups, like in break-junction experiments, the molecule is strongly bounded to the contacts

via thiolate groups and the lifetime of the electrons is small enough and we can assume that

we are in the weak electron-phonon coupling regime.

The phonon Green’s functions are the free phonon propagators enlisted in Eqs. (3.59).

The latter approximation means that the phonon population is regarded as a simple set of

Einstein oscillators with infinite lifetime. Moreover, we consider them in thermodynamic

equilibrium. The two main consequences are that the phonon population for every mode,

, is described by a simple Bose-Einstein distribution and, second, that the relaxation of

the phonons in the contacts is neglected.

5.1. Approximations in the electron-phonon code

Another approximation is related to the real part of the electron-phonon self-energy. As

explained the self-energy contains a real and an imaginary part where the former represents

a shift of the energy levels and the latter a lifetime of the particle induced by the electron-

phonon interaction. The real and imaginary parts for the retarded Green’s function are

correlated by the Kramers-Kronig relations which states the following:

Re(Σr

ph) = 1

PZ+∞

−∞

dE0Im(Σr

ph)

E0−E.(5.1)

In the calculations, this real part is neglected because the effects of the shift in energy

induced by the contacts is much larger than the one induced by the electron-phonon inter-

action.

The evaluation of the normal modes, vibrational vectors and frequencies, is performed

getting the relaxed geometry after a relaxation of the system. Later on the equilibrium po-

tential for the distribution of atoms is perturbed up to the second order defining a parabolic

approximation for the potential in the limit of small vibrations, the so called harmonic ap-

proximation. In the evaluation of the spectrum is fundamental to define which atoms are

allowed to vibrate. It is obvious that, until we do not allow the relaxation of the phonons

from the device to the contacts, the knowledge of the phonon dispersion function for the

bulks of the contacts is not required. What is more critic is the decision of which part of

the device can vibrate, that means if also some of the leads atoms, the one closest to the

interface between the molecule and the leads, should vibrate. In fact the vibrations of the

surface and the molecule together can have important effects on the total conduction, like

conduction modulation, where the charging and discharging of the molecule induced by

vibrations can change dramatically the conduction [87]. In all the calculations presented

in Chapters 6-9only the molecules were allowed to vibrate. This approximation can be

particularly questionable when the atoms in the leads have more or less the same mass of

the atoms in the molecular region (sulphur, carbon mainly). In that case it is obvious that

a vibration of the molecule can easily stimulate a vibration of the contacts. However, in

many experimental setups of interest, the contacts are made of gold that means that the

ratio between the masses of gold atoms and the sulphur junction atoms of the thiol group is

very large. It is possible to assume that the atoms in the gold are frozen in their equilibrium

positions.

Finally some words must be spent about the modeling of the external bias. The effect

of the voltage is equivalent to the one of an external electric field. The bias drops linearly in

the device regions, the effect of polarization induced by charge arrangements is neglected.

This effect, that in return changes the bias drop profile, can be included only by a full

treatment of the Poisson equation.

Chapter 5. The Electron-Phonon Code

5.2 The scheme of the device and the open boundary con-

ditions

The real system which is simulated in the code looks like the one depicted in Fig. (5.1).

The device region is formed by the molecule (M) plus the two leads (L). For metal con-

tacts, due to the good charge screening, the leads reduce to a couple of layers of atoms.

Periodic boundary conditions must be used to ensure that the principal layers (PL) of the

electrodes behave like in a bulk. The periodic boundary conditions in the direction of the

current produce problems when a real bias is applied between the two contacts. So, the pe-

riodicity in that direction must be used only for the evaluation of the contacts self-energy

and discarded after for the calculation with the bias applied.

Figure 5.1: A view of a typical molecular electronic device. In the picture different regions can be

recognized, the molecular (M) and the surfaces or leads regions (L) which together form the device

region (D). The two semi-infinite contacts (CLand CR), the charge reservoirs, are represented by

two chunks of the real electrodes. This chunks are formed by two principal layers (PL) each. The

PLs are the basis to evaluate the self-energies for the open boundary conditions.

The partitioning of the system as shown in the figure is reflected in the total DFTB

Hamiltonian matrix as:

0=



HDJDL JDR

J†

DL HL0

J†

DR 0HR

(5.2)

where the three diagonal blocks are the Hamiltonian parts for the device, left contact and

right contact respectively and the off-diagonal terms are the coupling between the three.

5.2. The scheme of the device and the open boundary conditions

Clearly, there is no coupling between the two contacts directly, this ensures that the self-

energies for the contacts are two independent contributions.

The task to evaluate the self-energies is fundamental in transport calculations. The two

self-energies represent a sort of equivalent potential, principally located in the junctions

between the device and the contacts, which record all the perturbation effects induced

by the semi-infinite contacts. To evaluate the self-energies we can expand the Dyson’s

equation for the zeroth order retarded propagator in matrix form

(

S−H)Gr

0=1(5.3)

where

= (E+i

)and

is an infinitesimal positive quantity. We can expand the latter

equation





SD−HD

SDL −JDL

SDR −JDR

S†

DL −J†

SL−HL0

S†

DR −J†

DR 0

SR−HR

˙





D,0Gr

DL,0Gr

DR,0

LD,0Gr

L,0Gr

LR,0

RD,0Gr

RL,0Gr

R,0



=1(5.4)

where we have emphasized the block structure of the total Green’s function matrix for the

entire system. However, we are interested in the device block Gr

D,0only. To get that part

we can start from the propagators of the isolated contacts and device, assuming that we can

ideally split the system in three non interacting parts:

L= (

SL−HL)−1(5.5)

R= (

SR−HR)−1(5.6)

D= (

SD−HD)−1(5.7)

The computation of the two free propagators for the semi-infinite electrodes is explained in

more details in Appendix B. Expanding Eqn. (5.4) for the device block and manipulating

them we get the three equations:

LD,0=−gr

S†

DL −J†

DL)Gr

D,0(5.8)

RD,0=−gr

S†

DR −J†

DR)Gr

D,0(5.9)

(

SD−HD)Gr

D,0+(

SDL −JDL)Gr

LD,0+(

SDR −JDR)Gr

RD,0=1(5.10)

Chapter 5. The Electron-Phonon Code

Substituting Eqn. (5.8) and Eqn. (5.9) in Eqn. (5.10) we get at the end the following:

D,0= [ESD−HD−(

SDL −JDL)gr

S†

DL −J†

DL)−(

SDR −JDR)gr

S†

DR −J†

DR)]−1

= [ESD−HD−Σr

L(E)−Σr

R(E)]−1.(5.11)

The latter is the zero order Green’s function with the self-energy for the open boundary

conditions: Σr

L(R). The self-energies depend on the coupling between the device region

and the contacts and on the free propagators of the contacts. The advantage of using a

linear combination of atomic orbitals basis set is particular evident here where we need

that the coupling blocks do not extend too much, but are quite localized in the interface

between contacts and device.

The two PLs in Fig. (5.1) are used to evaluate the free propagators of the contacts

for the evaluation of the self-energies and after are not take into account anymore in the

proceeding of the calculations. The self-energies in fact map the semi-infinite contacts on

the device region. So fort all the matrices implemented for the calculation of the current

have the dimension of the device region only.

5.3 The electron-phonon self-energies

For the electron-phonon interaction we can use the perturbation expansion starting from

the zeroth order propagators including only the open boundary conditions. The approxi-

mation implemented here is the BA, in its consistent or not consistent form. The diagram

representations are shown in Fig. (5.2, a), BA, and Fig. (5.2, b), SCBA.

Figure 5.2: The diagrammatic representations of the self-energy for electron-phonon in Born ap-

proximation (a) or in self-consistent Born approximation (b).

In analytic form the lesser and greater self-energies are the following

Σ<,>

ph (E) = ∑

Z∞

−∞

G<,>

0(E−E0)M

D<,>

,0(E0)dE0(5.12)

5.3. The electron-phonon self-energies

for BA and

Σ<,>

ph (E) = ∑

Z∞

−∞

G<,>(E−E0)M

D<,>

,0(E0)dE0(5.13)

for SCBA. The only difference in the two self-energies is that in one case we use the zeroth

order propagators, and in the second case the full renormalized G<(>)defined as

G<,> =Gr[Σ<,>

L+Σ<,>

R+Σ<,>

ph ]Ga.(5.14)

The evaluation of the integration in Eqn. (5.12) and Eqn. (5.13) is quite expensive compu-

tationally. However, following the approximation that our phonon population is described

by a collection of free oscillators, we obtain a simple formula for the phonon lesser and

greater propagators,

,0=−2

i[(N

+1)

(E−

ωα

)+N

αδ

(E+

ωα

)] (5.15)

,0= [(N

+1)

(E+

ωα

)+N

αδ

(E−

ωα

)].(5.16)

Under those approximations Eqn. (5.12) reduces to,

Σ>

ph =∑

+1)M

G>(E−

ωα

G>(E+

ωα

(5.17)

Σ<

ph =∑

+1)M

G<(E+

ωα

G<(E−

ωα

.(5.18)

The retarded and advanced propagators should be also renormalized taking into account

that, until the matrices are symmetric, it is possible to use the identity

i(Σ>

ph −Σ<

ph) = i(Σr

ph −Σa

ph) = −2Im(ΣR

ph)(5.19)

the symmetry of the matrices in fact ensures that (Σr

ph)†= (Σr

ph)∗=Σa

ph. The iden-

tity (5.19) can be used to compute the imaginary part of the retarded self-energy for the

electron-phonon interaction. The real part is neglected. Neglecting the real part we are

neglecting all the polaronic effects induced by the interaction between the charge and its

environment. The shift in the device spectrum is in fact zero, and the only effect of the

vibration is a correction to the electron lifetime.

The lesser and greater self-energies for electron-phonon interaction represent, as for the

equivalent quantities of the contacts, a flux of injecting and outjecting of charges, respec-

tively. In fact the scattering processes can be seen as a third virtual contact, see Fig. (5.3),

connected to the external phonon bath, with which it exchanges electrons. For charge

conservation the total current through this third contact must be zero for every mode,

=2e

hZ+∞

−∞

dETr[Σ<

G>−Σ>

G<] = Z+∞

−∞

(E)dE =0.(5.20)

Chapter 5. The Electron-Phonon Code

where J

(E)is a density of current for energy unit.

Figure 5.3: The schematic representation of the device connected to the two real electrodes and

the third virtual contact which connects the device to the external bath.

In few words, this means that for every charge that enter in the virtual contact, another

one must come out. However, the energy of the two charges is not the same, but differs

of the value ±

ωα

depending if it is describing an absorption or an emission process. This

means that also if the current I

is null, the power dissipated is not:

eZ+∞

−∞

E J

(E)dE 6=0.(5.21)

To evaluate the power dissipated by every mode of vibration the Meir-Wingreen equation

is invoked,

hZ+∞

−∞

dE E ¡Tr[Σ<

G>−Σ>

G<]¢(5.22)

where Σ<,>

ph =∑

Σ<,>

5.4 Computation of the electron-phonon couplings

As in most studies of vibronic coupling effects in molecular electronic, the small amplitude

of vibrational motion is usually invoked to approximate the electron-phonon couplings, see

[88]. In this case we can reduce the vibronic coupling to the expansion up to the first order

of the Hamiltonian of the device respect the displacements of the ions along every normal

mode of vibration

i j =∑

N∑

s¯

2mN

ωα

hi|

∂

,N|jie

,N(5.23)

5.5. The flowchart of the code

where Nruns over all the nuclei,

over the three cartesian coordinates, mNis the mass

of the Nth nuclei,

ωα

the frequency and e

,Nthe proper mode eigenvector element. The

expansion of the Hamiltonian can be further developed:

hi|

∂

,N|ji=

∂

hi|ˆ

HD|ji

∂

,N−hi0|ˆ

HD|ji−hi|ˆ

HD|j0i(5.24)

where hi0| ≡

∂

hi|/

∂

,Nrepresents the change in basis orbitals with displacements, and

using the identity

1=∑

i j |ii(S−1)i jhj|(5.25)

where Sis the overlap matrix, we get

hi|

∂

,N|ji=

∂

hi|ˆ

HD|ji

∂

,N−∑

nl hi0|ni(S−1)nlhl|ˆ

HD|ji−∑

nl hi|ˆ

HD|ni(S−1)nlhl|j0i(5.26)

The terms hi0|nican be evaluated numerically computing the first derivative of the overlap

matrix elements

∂

Sin/

∂

,Nusing the unperturbed and deformed geometries. Substitut-

ing Eqn. 5.26 in Eqn. 5.23 we get the final form of the electron-phonon coupling matrix

elements:

i j =∑

N∑

s¯

2mN

ωα

∂

i j

∂

,N−∑

∂

Sin

∂

(S−1)nlHD

l j −∑

in(S−1)nl

∂

Sl j

∂

,N!e

(5.27)

The importance to retain all the matrix elements is related to the fact that the inelastic

current can be strongly affected by motion that changes the overlap between different sites.

For example when two neighbor atoms in a molecular structure are oriented so that they

lie along the tunneling direction, a modulation of their distance can dramatically affect the

inelastic current.

5.5 The flowchart of the code

The final implementation of the code is shown in the flowchart in Fig. (5.4). After the

relaxation of the geometry, the DFTB code is used to compute the Hamiltonian H0

0for

bare electrons in the whole system without phonons, the modes of vibration of the device

(frequencies and eigenvectors) and the overlap matrix S. Before starting the evaluation

of the current a calculation at zero bias is done to evaluate the electron-phonon coupling

matrices M

. The matrices are stored to be used later.

Chapter 5. The Electron-Phonon Code

Figure 5.4: Flowchart of the code.

Now, it is possible to start the real calculation for the current. The very first step is

the definition of a set of parameters, that means the chemical potential of the two contacts

Land

R, related to the Fermi energy of the contacts at equilibrium EF. Other important

parameters are the applied bias Vand the temperature of the charge carriers Teand the

temperature of the device Td. The former affects the electron distribution in the contacts,

which in return changes the injection of charges in the device. The latter fixes the phonon

population.

5.5. The flowchart of the code

Concerning the numerical calculation there are two important parameters to set. One

is the energy grid of linearly spaced points ∆Eused for the numerical sampling. It should

span a sufficiently large energy range while at the same time resolve the variations of all

Green’s functions, self-energies, etc. The second is a parameter ∆which sets a convergence

criterium for the iteration procedure toward self-consistency. Finally, only in case the full

I-V characteristic is computed, there is the voltage step ∆V. Every calculation computes a

single current-voltage point and the power dissipation at that voltage. Using a script it is

possible to repeat the calculation for many points and get the full I-V curve.

The self-energies for the contacts and the bias are further used to compute the zero

order retarded Green’s function,

[ES−H0

0−Σr

L(E)−Σr

R(E)]Gr

0=1.(5.28)

Using the identity Ga= (Gr)†, the zeroth order of the lesser (greater) Green’s function can

be computed:

G<,>

0=Gr

0[Σ<,>

L+Σ<,>

R]Ga

0(5.29)

where

Σ<

L(R)=i(nL(R))ΓL(R)=i(nL(R))(−2)Im(Σr

L(R))(5.30)

Σ>

L(R)=−i(1−nL(R))ΓL(R)(5.31)

and nL(R)is the Fermi function for the left (right) contact.

The zeroth order Green’s functions are used to compute the self-energies for the electron-

phonon interaction. At this stage lies the difference between BA and SCBA. If we are in

BA the electron-phonon self-energies are used to renormalize the propagators, Gr

0,G<

0and

0and the Meir-Wingreen equation is evaluated for that voltage point. In case we are in

the SCBA, the new propagators are used to recompute the electron-phonon self-energies

and the scheme is iterated until convergence is reached. Finally, the current, the power

dissipation and the density of states are evaluated,

I(V) = 2e

∑

i=1

Tr[Σ<

L(Ei)G>(Ei)−Σ>

L(Ei)G<(Ei)]∆E(5.32)

(V) = 2

∑

i=1

Ei³Tr[Σ<

ph,

(Ei)G>(Ei)−Σ>

ph,

(Ei)G<(Ei)]´∆E(5.33)

(E,V) = Tr[i(Gr−Ga)] = Tr[A(E,V)].(5.34)

In case the number of modes included in the calculation is set to zero, the code reduces to

evaluate the elastic current via a Landauer-B¨

uttiker equation.

Chapter 6

Power Dissipation at Low Temperature

in Molecular Electronic Devices

The first application of the electron-phonon code is the evaluation of the power dissipated

in molecular devices. Due to the approximations enlisted previously the simulations can

evaluate dissipation at low temperature only.

The dissipation in molecular electronic devices is certainly a very important issue for

the future of the field. In order to generalize the method it is fundamental to extend the

technique of the gDFTB code also for higher temperatures. It is in fact clear that these

kind of simulations will become fundamental in the future to test the stability of molecular

junctions at work temperatures inside real electronic circuits. Recently an extension of the

method, to take care of real non equilibrium population of phonons in the molecular region

due to the coupling with the contacts, has been developed by Pecchia et al. [56].

Intuitively we can think that electron-phonon interaction is one of the main inelastic

mechanism responsible for thermal relaxations of electrons in molecular junctions, its un-

derstanding is important for controlling not only the stability of the device, but also for

other fundamental aspects. In fact, knowing were the dissipation is the highest in the

molecule can suggest how to improve, with molecular substituents, the performance of the

device.

In this Chapter is investigated the dissipation in octanedithiol molecules sandwiched

between two gold electrodes.

Chapter 6. Power Dissipation in Molecular Electronic Devices

6.1 Dissipation in alkanethiols

6.1.1 Geometry and vibrational modes

Although a simple and clear understanding of transport mechanisms has not been reached

yet for most of the molecular compounds, transport in alkanethiols is essentially at a ma-

ture stage of experimental development. Such type of molecules align on Au surfaces via

covalent S-Au bonds to form regular and stable self-assembled monolayers (SAMs). They

are characterized by a large optical band gap ( >5 eV), making them very stable to pho-

todegradation. The same large gap is responsible for a very low conduction via tunneling

in Au-thiol-Au structures, giving good electrical stability. For the reasons above, it is not

surprising that experiments of conduction through these compounds have been reproduced

by many research groups [89,90]. This paragraph is focused on octanedithiol between two

(111) gold surfaces.

The gold atoms, composing the two contacts, are kept at fixed positions, corresponding

to an ideal fcc crystal with Au-Au separation of 2.884 ˚

A. First the octanedithiol saturated

with a hydrogen termination was relaxed on top of one gold surface comprising six atomic

layers. After this first step, the hydrogen was removed and the second Au surface was

put, taking care that the Au-S distance obtained at the first interface was reproduced at the

second interface. Then this system was relaxed again. Periodic boundary conditions were

used in all these calculations. The final relaxed structure is shown in Fig. (6.1). Although

the global minimum of the isolated octanedithiol is the one with the highest symmetry

configuration, this symmetry is broken in the relaxed structure of the molecule attached to

Au electrodes.

Figure 6.1: Diagram representing the relaxed atomic coordinates of the octanedithiol between Au

contacts.

The vibrational frequencies are slightly affected by this, since many degeneracies are

lifted. However, the differences in the whole spectrum of vibrations between the molecule

in gas phase and relaxed on gold are within 20%, which is inconsequential for the present

analysis. The S atom is found to form a bond with an energy minimum at the hollow

site of the Au(111) crystal. Actually, the exact minimal configuration is not completely

6.1. Dissipation in alkanethiols

understood, it is found with the sulfur atom slightly shifted from the hollow position, as

reported in other ab-initio DFT calculations [91,92]. The sulfur atom is found at 2.76 ˚

from the Au plane. After the relaxation of the molecular coordinates, obtained by imposing

the atomic forces to be smaller than 10−4a.u., the computation of the vibrational modes

and frequencies was done.

6.1.2 Power dissipation in the molecule

While crossing the system, the electrons interact with the molecular ionic vibrations from

which they can be inelastically scattered. In the present context the assumption was made

that the metal ions do not move; henceforth, the electron-phonon scattering within the leads

is neglected. This is a reasonable approximation given the very large difference in atomic

masses between the Au atoms and the organic elements. It is clear that this assumption

must be relaxed if we handle different kind of contacts, like silicon contacts.

Figure 6.2: (left) Tunneling probability as a function of injection energy across a molecule of

octanethiol. The Fermi energy corresponds to the Au contacts. (right) The computed I-V character-

istics of the system (solid line) compared to experimental results (see [90]).

The separation between the molecular and the contact regions is obtained by two leads

composed by three layers of gold atoms representing the metal surfaces that have been

treated as to be part of the molecular region. Because of the short range of the atomic

orbital interactions and since the Au atoms do not move, such three layers ensure that the

contact and electron-phonon self-energies act on orthogonal subspaces of the Hamiltonian

Chapter 6. Power Dissipation in Molecular Electronic Devices

matrix. In this calculation the Born approximation (BA) was applied for the calculation

of the phonon self-energies. After that a self-consistent solution (SCBA) was done for

comparison.

Figure 6.3: Calculated (a) Scattering rate and (b) power dissipated for each vibrational mode.

Since I am mainly interested on inelastic phonon emission, the approximation of the

phonon self-energy to be purely imaginary, neglecting the shift of the van-Hoove singulari-

ties (i.e., essentially neglecting energy shifts due to polaron-like states), is reasonable. This

assumption is valid since the molecular levels are far from the relevant range of injection

energies and therefore the small shift will not modify appreciably the tunneling current.

In any case, the BA is valid in the small electron-phonon interaction limit and far from

resonant situations where polaron effects should be taken into account. Using the mathe-

matical machinery of the NEGF, the tunneling current through an alkylthiolate molecule,

Au-S(CH2)8S-Au, has been computed. Applying the Meir-Wingreen equation to get the

current we obtain the curve shown in Fig. (6.2). The current computed is compared to

experimental data obtained by measuring the current through a SAM assembled within a

nanopore [90].

The simulation has been performed on a single molecule connecting the two con-

tacts, assuming that the SAM assembled in a lattice (√3×√3)R30◦, from the measured

6.1. Dissipation in alkanethiols

nanopore diameter of 45 nm, it is possible to estimate that approximately 10,000 molecules

are sampled in parallel. Accordingly, the experimental measurements are scaled by a factor

10−4in order to compare with the calculations. The order of magnitude of the tunneling

current is predicted very well. This is important, particularly in relation to the absolute

magnitude of the power dissipated in such system, as discussed later.

For each of the 78 modes, the contribution to the inelastic current at T=0 K and for

an applied bias of 2 Volts has been computed . At T=0, the inelastic part of the tunneling

probability is related to the net phonon emission rate by

−1

inel/e, where the

index

stands for the vibrational mode. This quantity is used as a measure of the electron-phonon

strength and is used to compare the different modes. A summary of such computations is

reported in Table (6.1) and in Fig. (6.3, a). The lowest vibrational modes correspond to

oscillations of the carbon atoms in the backbone plane, resembling the first harmonics of a

string, and rigid twist around the C-C bonds of large subunits involving two or more CH2

groups.

ω τ

−1

description

exp.

34.8 0.608 (CH2)2twist

106.9 0.127 (CH2)2twist

109.3 0.189 (CH2)2twist

178.8 0.157 (CH2)2twist

217.0 0.314 Au-S stretch 2251-

320.6 0.470 Au-S stretch -2551

687.5 0.383 S-C stretch 650-7062

693.6 0.462 S-C stretch

848.1 0.444 CH2rock 715-9251,3

1076.8 0.848 C-C stretch 10501

1136.2 0.871 C-C stretch 11202

1222.8 0.097 CH2wag 13303

1259.9 0.626 CH2twist 12651

1479.2 0.147 CH2twist

1518.2 0.082 CH2scissor 14551

1881.0 0.243 C-C + H swing

2976.5 0.002 C-H stretch 28601

3421.4 0.200 H-Au stretch

Table 6.1: Summary of the most important calculated frequencies (cm−1) of the octanedithiol,

phonon emission rate (ps−1), followed by a brief description of the mode and the typical experi-

mental value, obtained by (1) HREELS [93], (2) Raman [94]or (3) IR [95].

The two sulfur atoms remain practically fixed or slide slightly over the gold surface.

The frequencies of such modes are affected by large relative errors because they are sensi-

tive to small differences in the atomic pair potentials and to the relaxed geometry. Further-

more, these modes are not easily accessible experimentally. We find that modes associated

to internal twists of the C-C backbone involving the motion of (CH2)2subunits are more

Chapter 6. Power Dissipation in Molecular Electronic Devices

effective in scattering electrons. As expected, the Au-S stretching modes, which are found

in the range between 250 cm−1to 320 cm−1, in agreement with experiments [90], give a

large contribution to the electron-phonon scattering. Similarly, a large contribution is given

by the C-S stretch modes, found around 690 cm−1. The band of modes found between 1000

and 1100 cm−1corresponds to C-C stretch modes. These give the largest contribution to

the inelastic current. A large contribution comes also from the modes related to motions of

the CH2modes, particularly rocking and twisting. Wagging and scissoring modes, instead,

give a smaller contribution.

In the simulation there is also a non negligible signal from some vibrational modes

that have not been discussed in experimental papers, since they are special to the Au-

octanedithiol-Au system and are not seen in isolated molecules. In particular, two modes

around 1900 cm−1(0.23 eV) are associated with a slight rotation of the C-S bond, a stretch

of the C-C bonds involving the motion of the C at position close to sulfur, and a pronounced

swing of the H atom closest to the gold surface. Another mode affecting the inelastic

current is related to the oscillation of the hydrogens closest to the Au surfaces. Such modes

could be seen as a Au-H stretch and are found at a frequency of 3400 cm−1(0.42 eV).

Figure 6.4: (left) Tunneling probability as a function of injection energy across a molecule of

octanedithiol. The contribution of the inelastic (dashed) and elastic (solid) tunneling current are

shown. The black curves are the first-order BA, the red ones are the SCBA. (right) The SCBA I-V

characteristics of the system. The continuous line is the total current, the broken line is the elastic

component and the dotted line the inelastic component.

The 34 modes giving an inelastic rate greater than 108 ps−1have been included for the

calculation of the inelastic current. Figure (6.4) shows the elastic and inelastic contributions

6.1. Dissipation in alkanethiols

to the total tunneling current resolved in energy for the bias of 2.0 Vand the calculated

I-V characteristics of the system. The peaks in the transmission function correspond to

features of the surface density of states of gold. In Fig. (6.4) the first-order BA has been

compared with the SCBA. It is important to remark that for this system the first-order

Born approximation gives already an acceptable result and the self-consistent loop does not

introduce substantial changes in the transmission probability in the relevant energy range.

Differences are appreciably deep in the energy gap, where the transmission probability is

already small and the total contribution to the current is negligible.

At this point it is possible to compute the amount of power dissipated in the molecule

due to inelastic phonon emission. This calculation can be obtained by considering the

virtual contact current. The power dissipated can be easily calculated by Eqn. (5.21). This

quantity is the virtual contact current, simply representing the total current scattered from

its original energy to a new one after phonon emission (or absorption). The calculation

shows that the power dissipated in the octanedithiol isW=0.16 nWatt (1 eV/ns) at 2 Volt of

applied bias. It is interesting, however, to analyze the power dissipated in each vibrational

mode, as shown in Fig. (6.3, b). This is not directly proportional to the emission rate, since

the power dissipated depends on the phonon energy as well. As expected, the low frequency

modes contribute little to the dissipation, despite the relatively large scattering rates. The

modes giving the largest contribution are found in the band of C-C stretch modes, the C-S

stretch, and the CH2rocking modes. Considerable power is also absorbed by the modes at

1900 and 3400 cm−1, involving essentially movements of the hydrogens close to the Au

surfaces.

Chapter 7

Simulation of IETS in Alkanethiols

In this Chapter an application of the electron-phonon code to simulate Inelastic Electron

Tunneling Spectroscopy (IETS) measurements is presented. The systems under investi-

gation are self-assembled monolayers of alkanethiols between two gold electrodes [90].

At least three papers have investigated the IETS of alkanethiols on gold. First, a sophisti-

cated theoretical method [96] was applied to varying length alkanethiols [42] to give results

which reproduced experimentally observed odd-even effects with varying chain lengths.

Second, by making use of a model based on perturbation theory [45] it was shown that it is

possible to get informative predictions of IETS showing good agreement with experiment

[97]. Third, a number of different geometries for alkanethiol binding were considered with

an examination of the impact of this binding on the IETS [46]. Here, I continue from this

work and look further at the dependence on molecular structure and how this relates to the

understanding of IETS.

The Chapter is organized as follows: after an introduction about the calculation of the

geometry and of the vibrational frequencies, the IETS spectrum is discussed. The analysis

section is divided in four parts. First, it is shown the general features of the spectrum.

Second, it is looked at the unambiguous assignment of this spectrum that is achievable

through simulations. Third, it is examined the nature of the low-voltage conduction channel

and its control of the IETS signal. Finally it is shown the nature of the vibrational modes

that produce the IETS.

7.1 IETS approximation

IETS is defined as d2I/dV2. To obtain d2I/dV 2we should take, from the Meir-Wingreen,

I(V) = Z+∞

−∞

J(E,V)dE (7.1)

Chapter 7. Simulation of Inelastic Electron Tunneling Spectroscopy in Alkanethiols

and derive two times respect the bias. However, in IETS the bias applied is very small, so

that we can consider that the density of current is independent from the bias: J(E,V)≈

J(E). Moreover, for small electron-phonon coupling the integration limits reduce to the

chemical potentials of the two contacts. Considering one of the contact as the reference

one, we get:

I(V)≈ZEF+eV

J(E)dE =ZEF/e+V

EF/e

eJ(e

(7.2)

where it has been made a simple variable substitution in the last equation:e

=E, with e

the electron charge. The first derivative of the current respect to the bias becomes:

dI(V)

dV ≈eJ(eV)(7.3)

Iterating the latter for the second derivative and with a further variable substitution E=eV,

it is obtained

dI2(V)

dV2≈e2dJ

dE .(7.4)

This final relation has been used to generate the IETS. This approximation neglects the

voltage dependence of the electronic structure, a property considered to be minimal at the

low voltages considered.

7.2 The choice of binding site

The binding geometry of alkylthiol molecules chemisorbed (with terminal proton loss from

sulfur during the adsorption process) on gold surfaces has received considerable attention

in the literature [98,99,100,101] without any clear conclusion. Calculations performed

on flat Au(111) surfaces [102,103] suggest different geometries from those resulting from

thiol-induced surface reconstructions [82]. In this section I present results for a represen-

tative geometry which is consistent with recent experimental results [104].

In Chapter 6was already discussed a geometry for octanethiol between two gold con-

tacts. However, for the IETS simulation a new geometry has been computed. A different

relaxation scheme was adopted to preserve the symmetry in the junction between molecule

and gold at the two surfaces.

The geometry of the molecule in the junction is determined in two steps. First, an

optimized geometry is obtained for chemisorbed octanethiol or in the case of the other

chain lengths the monothiol of that alkane chemisorbed through the terminal sulfur to a top

site on a single Au(111) surface; while experimentally it is likely that sulfur does indeed

sit vertically above a gold atom, it is probable that it is some reconstructed surface instead

of a pure (111) surface. Periodic boundary conditions are used, with the chemisorbed

molecules sufficiently far apart to be considered isolated [a p(5×5)unit cell].

The geometry for the full electrode-molecule-electrode system was then generated by

7.2. The choice of binding site

symmetrizing about a point of inversion between the central bond of the alkane chain to

give a chemisorbed alkanedithiol bound to two cofacial Au(111) surfaces. Optimized ge-

ometries for octanedithiol in two slightly different local minima are shown in Fig. (7.1).

The lowest-energy structure has the octane backbone at a tilt angle of 28◦(from the Au

surface normal), reducing to 18◦for the other. These values differ significantly from what

has been reported experimentally [105,106,107] for condensed monolyers, supporting the

proposal that this property is controlled by interchain and solvent interactions. Detailed

structures for molecules sandwiched between two electrodes remain unknown.

Figure 7.1: Two optimized geometries of octanedithiol between the two electrodes. The top geom-

etry is the lower-energy minima of the two.

7.2.1 Atomic partitioning of the system

In the case of highly conducting conjugated molecules it is necessary to include a number

of gold atoms in the extended molecule in order to accurately calculate the conductance

of these systems. If no gold atoms are included in the extended molecule the self-energies

do not shift and broaden the molecular states sufficiently to accurately model the coherent

conduction channels that link the interface regions of the real semi-infinite systems. While

the inclusion of a large number of gold atoms in the junction region significantly enhances

the calculation of coherent transport, it does so by introducing a discrete approximation to

the notionally continuous distribution of the low-voltage interface to interface conduction

channels.

Chapter 7. Simulation of Inelastic Electron Tunneling Spectroscopy in Alkanethiols

While the effects of excluding the interface region may be more important for con-

duction through conjugated molecules, much smaller effects are expected for alkanethiols.

The smaller band gap HOMO-LUMO in conjugated systems means that the conductance

varies significantly depending upon how the system is partitioned as even a slight shift in

the molecular energy levels is important for the conductance at low voltage.

For alkanethiols, however, the significantly larger band gap means that the exact char-

acter of the molecular energy levels at low voltage is not as important. Further, the IETS

is, as would be expected, very sensitive to changes in the number of conducting channels

available. In the low-voltage regime considered here ( <0.4 V) the density of coherent

conduction channels should be voltage independent, as it only depends upon the gold den-

sity of states; it is the change in the number of accessible incoherent channels as the voltage

increases that gives the IETS signal. This will produce a large and broad feature underlying

the fine structure of the IETS signal. For this reason we do not include any gold atom in

the extended molecule.

7.2.2 Correction of the DFTB vibrational frequencies

The DFTB method is known to overestimate harmonic frequencies [62] as do other theo-

retical methods. In order to have sufficient accuracy in the calculated spectra to confirm

experimental assignments for observed peaks, it is necessary to correct the DFTB frequen-

cies with a scale factor

. This factor is generated for a gas-phase octanedithiol molecule

by the procedure described by Scott and Radom [108] as

=∑all

theor

exp

∑all

theor

i)2(7.5)

where

theor

iand

exp

iare theoretical harmonic and experimental frequencies, respectively.

In this case, however, due to the lack of experimentally assigned gas-phase frequencies

for the octanedithiol, the experimental frequencies were substituted with B3LYP6-31g(d)

frequencies calculated using GAUSSIAN 03 which had first been corrected with the rele-

vant scale factor. This gave a scale factor of 0.913 for the DFTB frequencies. Implicit in

this procedure is the assumption that the majority of the error in the calculated frequencies

arises due to intrinsic limitations in the electronic-structure methodology rather than with

the specific treatment of the interaction of the molecule with the gold surfaces. All fre-

quencies presented in this IETS spectrum, as well as all figures, have been corrected with

this scale factor.

7.3 Discussion of the IETS simulations

The observed and calculated IETS of chemisorbed octanedithiol are shown in Fig. (7.2)

as a function of the applied voltage, 0.1 Vcorresponds to 806.5 cm−1vibrational wave

number. The calculated spectrum is obtained from the low-energy structure shown in

7.3. Discussion of the IETS simulations

Fig. (7.1). Both spectra depict a variety of well-resolved peaks covering the vibrational

energy range, but they differ in terms of the spectral base line and the band shapes.

Figure 7.2: The calculated IETS for octaneditiol between two gold electrodes with the experimental

results from ref. [90] shown in the inset. Numbered peaks are assigned in Table (7.1).

A close inspection of the figure shows that the observed and calculated peaks occur

at similar voltages, and, while significant differences between the observed and calculated

intensities are apparent, the number of resolved bands is approximately the same. This

facilitates the assignment of the observed spectra, as discussed in detail later. The accuracy

of the calculated band positions and intensities depends most significantly on the accuracy

of the DFTB method used to calculate the electronic and nuclear structure. While in gen-

eral density-functional tight-binding calculations such as these perform very well for the

calculation of geometries and vibrational frequencies, latter improved by the use of the

scale factor, errors in calculated energy levels do arise. As the intensities of calculated

IETS bands depend significantly on the values for the energy levels, lower accuracy is, in

general, expected for this property.

Previously [90], the IET spectrum of octanedithiol has been assigned by comparison

of the observed bands with vibration frequencies obtained by IR spectroscopy, Raman

spectroscopy, and HREELS. This assignment is not a straightforward process as the other

techniques have specific selection rules, unlike IETS. It has been observed previously that

in IET spectra both IR-active modes and Raman-active modes can be seen as well as addi-

tional modes, although not all IR-active or Raman-active modes may be seen. The experi-

mental assignment is indicated in the inset of Fig. (7.2) and comprises eight firmly assigned

molecular vibrations as well as seven tentative assignments, marked “ ∗”, to vibrations in

the surrounding Si3N4matrix. For the calculated spectrum shown in Fig. (7.2), simple

assignments arise for all bands except for that labeled 1, and these are given in Table (7.1)

where they are compared to the experimental assignments.

Chapter 7. Simulation of Inelastic Electron Tunneling Spectroscopy in Alkanethiols

Calculated spectrum Experimental spectrum

Peak Volt. (V)

(cm−1)Assignment†Volt. (V)

(cm−1)Assignment

1 0.007 61 Many modes

2 0.018 151 C-C-C scissor 0.015 121 *

3 0.025 204 Au-S stretch 0.033 266 Au-S stretch

4 0.036 302 S-C-C scissor

5 0.051 424 C-C-C scissor 0.058 468 *

6 0.080 648 C-S stretch 0.080 645 C-S stretch

7 0.091 737 CH22 in-plane rock

8 0.102 825 CH2in-plane rock (all) 0.107 863 in-plane rock

9 0.120 972 CH2in-plane rock (all)

10 0.142 1149 C-C stretch (cent.) 0.133 1073 C-C stretch

11 0.150 1208 C-C stretch (ext.)

12 0.161 1296 CH2twist (ext.) 0.158 1274 CH2wag

13 0.179 1444 CH2scissor (ext.)

14 0.189 1525 CH2scissor (cent.)

15 0.321 2593 C-H stretch sym. (ext.)

16 0.333 2688 C-H stretch asym. (cent.)

17 0.343 2769 C-H stretch sym. (ext.)

18 0.348 2806 C-H stretch asym. (cent.) 0.357 2879 C-H str. asym.

Table 7.1: The assignment of the peaks observed in the calculated IET spectra shown in Fig. (7.2).

†: the labels (ext.), (cent.) and (all) stand for “extremities”, “central” and “all the molecule”

respectively and are related to the localization of the mode of vibration.

Of the seven bands originally tentatively assigned to the Si3N4matrix, five occur above

0.25 eV (2000 cm−1) in a region in which the molecule does not have vibrational fun-

damentals. Any molecular band in this region would thus have to be a combination or

overtone band. As such bands have never been observed to produce a strong signal in IET

spectra [109], they are explicitly excluded from the electron-phonon coupling and hence

the calculations provide no new information concerning the nature of these five observed

modes. However, a variety of molecular peaks are predicted in regions (0.01 and 0.07 V)

of the other two peaks that were not originally authoritatively assigned. The calculations

thus suggest new assignments of these peaks.

Of the eight bands positively identified experimentally, seven correspond to calculated

bands of the same nature, as indicated in Table (7.1); however, the calculations enhance

these assignments by providing more specific information. For example, the most intense

band in the observed and calculated spectra was originally assigned as a C-C stretch but the

calculations indicate that the various C-C stretches within octanedithiol should have differ-

ent frequencies and very different intensities, allowing the observed band to be assigned to

C-C stretch modes at the extremities of the alkane chain. Similarly, the calculations reveal

that the observed CH2stretch and rock modes are localized on the extremities. These con-

clusions are obtained following the examination of the IETS spectrum shown in Fig. (7.3)

calculated at much higher resolution than the experimental spectrum.

Due to the effect of modulation broadening, high-resolution experimental spectra with

7.3. Discussion of the IETS simulations

36 clearly identifiable modes will be difficult to obtain and hence calculations offer a long-

term solution to the assignment problem of IETS. While the calculations indicate that the

vibrational modes at the center and extremities of the molecule are distinguishable, they

also indicate that most of the major peaks shown in low resolution in Fig. (7.2) actually

arise from multiple modes as indicated in Fig. (7.3).

Figure 7.3: The calculated IET spectrum at high resolution showing the contribution of individual

modes to the broad peaks observed experimentally.

For all bands other than peak 1, all contributing modes are modes of the same type so

it is possible to assign the peaks. For peak 1, the contributions are from different types of

modes associated with translation of the carbon chain in the two directions parallel to the

gold surface, the rotation of the molecule about the chain direction, and some asymmetric

combinations of these modes about the center of the molecule.

The only region of the spectrum where the calculated assignment is not fully consistent

with the experimental assignment is part of the fingerprint region (1200-1500 cm−1) in

which a variety of different types of modes occur such as CH2wags and CH2twists that

often mix strongly with each other. The most prominent disparity is the assignment of

peak 12 as a CH2twist in the calculated spectrum compared with the CH2wag assigned at

1274 cm−1in the experimental spectrum. Throughout this region DFTB predicts a strong

Duschinsky rotation between the modes, and in addition DFTB and B3LYP calculations

differ significantly in their perception of the mixing for gas-phase octanedithiol.

While all calculations predict the IETS signal for the twist to be stronger than that for

the wag, I find that the calculated twist frequency is very sensitive to interface geometry,

leading to the possibility that structural irregularities wash out the signal from the twist.

Computationally, this makes a definitive assignment impossible without detailed knowl-

edge of the experimental structure, independent of the reliability of the method. Also,

Chapter 7. Simulation of Inelastic Electron Tunneling Spectroscopy in Alkanethiols

the arguments used to make the original experimental assignments of the IR and Raman

spectra may need to be reexamined.

7.3.1 Nature of the orbitals controlling conduction and IETS

IETS is controlled by the electron-phonon couplings M

µν

between atomic orbitals

and

in mode

as well as the Green’s function matrices lesser and greater that also act to

determine the coherent electron transport. Shown in Fig. (7.4, a) is the coherent transmis-

sion T(E)calculated using the Landauer formalism as in Eqn. (2.37) where the contact

self-energies are related to the ΓL(R)matrices.

Figure 7.4: The transmission (a) and IETS (b) for the normal system and the system where the

sulfur lone-pair and d orbitals have been shifted to very high energy, showing that these are the

dominant conduction channels at low voltage.

The peaks in this function correspond to conduction channels embodied within G<(>)

the general nature of such channels is quite complex and their effects on IETS will be de-

scribed in detail in Chapter 8. Here, I show that broad features can be identified with con-

duction involving the sulfur lone-pair and valence dorbitals by taking the DFTB-calculated

Hamiltonian matrix in the L¨

owdin orthogonalized atomic-orbital basis and artificially shift-

7.3. Discussion of the IETS simulations

ing these orbitals to remote energies. The modified coherent transmission is also shown in

Fig. (7.4, a) and it is evident that there is a dramatic lowering of the transmission near

-4 eV, attributable to the lone-pair orbitals, and 0 eV, attributable to the dorbitals. It is

hence clear that at the Fermi energy of the electrodes these orbitals constitute the primary

conduction channel(s).

As alkanethiols have very large band gaps (HOMO-LUMO gap) these primary conduc-

tion channels at low voltage are not resonant channels with a large amount of molecular

character, as they are for, say, 1,4-benzenedithiol, but rather are states that are primarily

localized on the electrodes with a small amount of molecular character. This is effectively

a conduction through the wings of the probability distribution functions representing the

molecular orbitals; these functions show significant broadening owing to their coupling

with the electrodes. Figure (7.4, b) shows the calculated IETS for both the original and

modified Hamiltonians, and it is evident that the majority of the IETS signal disappears

when the sulfur lone pair and dorbitals are shifted to very high energy.

Interestingly, the peak in the IETS resulting from the C-S stretch mode does not de-

crease proportionately with the rest of the spectrum when the orbitals are removed. The

reason for this is that this mode couples strongly with the high-energy transmission channel

at ca. -8 eV associated with the C-S

bond, the channel that produces the weak remnant

transmission at the Fermi energy depicted in Fig. (7.4, a) after the removal of the other

sulfur orbitals. Hence it is clear that the conduction channels that control the low voltage

coherent transport also control IETS. The signal obtained in IETS thus arises as a result

of energy loss from electrons in the low-voltage conduction channels as they excite the

vibrational modes of the molecule.

7.3.2 Nature of the vibrations that produce IETS

As the low-voltage conduction channel dominates IETS, the peaks of highest intensity will

be those associated with the vibrational modes that provide the greatest perturbation in

the region where this channel has electron density. As the main conduction channels for

octanedithiol have been determined to involve sulfur lone-pair and dorbitals as well as

some contribution for the C-S sigma system, vibrations localized on the S and terminal C

atoms are expected to be the most significant. It is this physical principle that underpins

the previously mentioned assignment of the dominant IETS modes as modes involving the

extremities of the bridging molecule.

This physical picture leads to the prediction that the IET spectra of alkanethiols will be

very sensitive to changes in the binding geometry and less so to changes in the backbone

of the chain. To test this hypothesis numerically, I calculate the IETS for alkanethiols of

modified chain length and different binding site. Figure (7.5) shows the calculated current-

voltage characteristics (left) and the corresponding IET spectra (right) for C4, C8, C12, and

C16 alkanethiols.

The geometry optimizations for these molecules were performed so as to achieve sim-

ilar binding configurations, and this is reflected in the similarity of the IETS features. The

current decreases approximately exponentially with increasing length as is expected and

as a result so does the magnitude of the IETS base line. Beyond the change in base line,

Chapter 7. Simulation of Inelastic Electron Tunneling Spectroscopy in Alkanethiols

Figure 7.5: The simulated current vs voltage and IETS for C4, C8, C12, and C16 length chains.

however, the differences are slight and simply reflect the subtle variation in the calculated

molecular binding geometry (for example, the slight rotation of the terminal CH2group).

This result shows that the vibrational modes in the central region of the molecule do not

contribute significantly to the IETS signal, as hypothesized.

This result is consistent with that found previously by Chen et al. [42] for the effect

of chain length on IETS characteristics. In that work the low-energy modes (below 50

meV) were predicted to be very much smaller in magnitude than the other features, in

contrast with the present results, however, the predicted length dependencies were similar.

As the treatment of the interface region in this previous work was quite different in that the

electrodes were treated with a jellium model, it can be concluded that the major conclusion

concerning the significance of the terminal vibrations is found to be robustly obtained.

As there are many minima on the potential-energy surface for binding octanedithiol

on gold it is possible to investigate two similar binding sites that feature the sulfur atoms

bound directly above an Au site. The first geometry is the low-energy structure that has

been used up to this point and the second is a higher-energy optimized geometry shown in

Fig. (7.1). In this second structure the tilt angle to the Au surface normal is now 18◦; and

there is a greater twist in the alkane chain; both geometries are shown together in Fig. (7.1)

for comparison.

Figure (7.6) shows the calculated IETS spectra for these two geometries, and these dis-

play qualitative differences that reflect the changes in the binding site in addition to the

changes found in the magnitude of the observed current (evidenced in the IETS by the

changes in the base line). The most significant change observed is the modulation of the

intensities of the C-S and C-H stretches and CH2in-plane rock modes between the two ge-

7.4. Conclusions

Figure 7.6: The calculated IETS for the two different binding sites corresponding to the two differ-

ent geometries shown in Fig. (7.1) The modulation of the C-S stretch and C-H stretch modes reflects

the changes in the geometry of the terminal CH2group with respect to the gold surface.

ometries. In the low-energy geometry we have the large C-S stretch peak corresponding to

sensitivity to this vibrational mode as the molecular component of the conduction channel

is predominantly sulfur density. Conversely, in the high-energy geometry with the greater

twist angle in the alkane chain a hydrogen atom of the first CH2group is now in close prox-

imity to the surface. This results in the molecular component of the conduction channel

having significant density on this hydrogen atom in addition to the sulfur, and in turn the

IETS signal shows sensitivity to the CH2in-plane rock and C-H stretch of the extremities.

7.4 Conclusions

The IETS of chemisorbed octanedithiol between two gold contacts shows to reproduce key

features of the experimental spectrum for this molecule, including the verification of the

major molecular vibrational modes involved. However, the calculations are shown to lead

to increased specificity in the assignments and to new assignments. The previous experi-

mental assignments were obtained by the comparison of the IETS results with those from

other spectroscopies, but the different selection rules and chemical environments make this

a difficult process. The calculations suggest that a number of the modes that were attributed

to possibly coming from the Si3N4matrix [90] may in fact have molecular character.

Further, the evidence of a mixing of modes with different symmetry in the fingerprint

region suggests that the detailed assignment of this region may prove illusory without de-

tailed knowledge of experimental structure. They show that for chemisorbed alkanethiols,

Chapter 7. Simulation of Inelastic Electron Tunneling Spectroscopy in Alkanethiols

IETS is dominated by vibrational modes localized at the ends of the molecule, providing

additional assignment information that is critical to the understanding of the nature of the

process. It may be possible to verify these features of the assignment if higher-resolution

spectra can be obtained, beyond the modulation-broadening limit, as such spectra are pre-

dicted to be able to discriminate in both frequency and intensity between central and ter-

minal vibrational modes.

The calculations indicate that there is a strong correlation between the conduction chan-

nels that dominate the low voltage coherent conductance of alkanethiols and those that

dominate IETS. These channels are shown to be dominated by the sulfur lone-pair and

other orbitals, leading to the result that IETS is most sensitive to vibrations localized at the

ends of the molecule. For conduction through conjugated molecules, the lower molecular

band gap results in greater participation of the center region of the molecule in the conduc-

tion channel, and so the dominance of terminal modes is not expected to be as pronounced.

For alkanethiols, however, the deduced sensitivity to the metal-molecule interface region

indicates that IETS spectra should show considerable variation with subtle changes in bind-

ing site and small variations with increasing chain length, and these predictions are verified

by the simulations.

Other calculations [46] of the IETS of octanedithiol have predicted significant differ-

ences in the relative intensities of the C-C and C-S modes to those shown in Fig. (7.2), most

likely as the result of the different interface geometries and treatments used. While opti-

mized structures were used in the calculations, many binding sites are possible on smooth

gold surfaces ([82]), and many irregular binding sites are likely to occur in real metal-

molecule junctions. Experimental spectra obtained in, say, nanopores involving thousands

of molecules are likely to depict an ensemble of possibly very different junctions, making

detailed interpretation difficult. The calculations are also approximate in that they do not

allow the gold surface atoms to relax, nor do they include gold surface vibrational modes.

Consequently, the predicted properties for terminal vibrational modes are not expected to

be as accurate as those for the central modes. Hence, while IETS are predicted to show

details revealing critical but largely unknown aspects of the alkanethiol-metal interface,

the design of computations that can authoritatively unravel these features poses a major

challenge.

The most striking feature of IETS results remains; this is the first technique to provide

conclusive evidence that current flows through a molecule, as opposed to a metal filament.

Chapter 8

The Role of Symmetry and Channels in

Conduction

The analysis of IETS in Chapter 7arises many questions. First, it was shown the diffi-

culty of a clear assignment of the peaks in the IET spectrum. This is related to two facts:

the uncertainty in the definition of the geometry of the molecule and the metallic surfaces

from experiments data and the different selection rules shown by IETS. Second, from the

analysis of the conduction in the octanedithiol molecule was clearly shown that only few

“channels” conduce all the current which is inelastically scattered by the most relevant

vibrational modes, and those channel can be partially related to some canonical molecu-

lar orbitals of the molecule. However, the canonical molecular orbitals are usually not a

suitable set of channels for transport problems, due to the presence of the open boundary

conditions, that means the coupling of the molecule with the contacts.

The present Chapter is focused on answering these two points. The Chapter is artic-

ulated as follows: in the first part it is shown how define a point group for the molecule

between the electrodes. It is clear that starting from the molecule in gas phase, the effects

of the contacts is to reduce the symmetry of the system in a subgroup of the original sym-

metry point group. The second part of the Chapter is devoted in the definition of channels

for transport. These channels decompose the total current in a sum of non interfering con-

tributions. Every channel can be labeled according to its irreducible representation from

the reduced point group.

During this Chapter it is made use of the Landauer-B¨

uttiker formalism, that means I

am considering elastic current only. In Chapter 9the formalism will be generalized to

include also scattering processes in the full Meir-Wingreen equation. This will lead to un

investigation of selection rules in IETS.

Chapter 8. The Role of Symmetry and Channels in Conduction

8.1 The definition of the symmetry of conduction

It is well known that one of the most powerful concepts to aid in understanding molecular

properties is symmetry. As researchers today expend more effort to utilize molecular prop-

erties in single molecule electronic devices, it is important to understand how symmetry

operates within these regimes.

Understanding the symmetry of a molecule carrying current is an important step in

molecular electronics for resolution of issues such as whether there can be symmetry se-

lection rules in IETS, or the means by which molecular orbitals influence conduction char-

acteristics. These issues in turn feed back to experimental design, providing guidance for

device characterization and synthetic targets. This section will show how symmetry oper-

ates within the standard conduction formalism of molecular electronics and from this illus-

trate the additional information that symmetry assignment provides concerning molecular

conduction.

In order to describe the symmetry of molecular conductance, it is necessary to under-

stand the relationship between the symmetry of each of the quantities within the Meir-

Wingreen equation in terms of the geometric point group of the entire system. Assuming

for simplicity the absence of an applied bias (V= 0), the first important consideration is that

the self-energies Σr

L(R)arising from the electrode-molecule interactions do not commute

with all of the symmetry operators that define the point group of the bridging molecule.

Moreover, the application of an electric field only increases the symmetry breaking.

First, a single electrode will couple strongly to one end of the molecule and weakly to

the other. This results in the individual Σr

Land Σr

Rmatrices not having an end-to-end sym-

metry even if such symmetry is present in the geometric point group of the full electrode-

molecule-electrode system. With end-to-end symmetry I mean that the symmetry of the

self-energies cannot contain all the symmetry operations of the device point group which

relate one end of the device to the other. Second, as Σr

L(R)embody the interactions between

the molecule and the left (right) electrode, the symmetry of these matrices embodies only

the symmetry elements that are common to both the molecule and the electrode junction.

These two levels of symmetry breaking provide a natural framework to define the sym-

metry of conductance. The symmetry of the transmission through the full system can be

defined exactly as the conductance point group; the symmetry operators that define this

point group will commute with all quantities appearing in the NEGF formalism. Hence

the conductance point group contains the common elements of the symmetries of Σr

Land

Σr

R, called the junction-conductance point groups, as well as those of the geometric point

group of the isolated molecule. Put another way, the conductance point group contains all

the symmetry of the two electrodes and the molecule and also all the symmetry breaking

interactions between the three.

While the conductance point group and the individual junction-conductance point groups

are important as they describe molecular conductance within the NEGF formalism exactly,

in practical applications the required precise junction geometries and electrode symmetries

are likely to be unknown and/or uncontrollable [110]. However, the dominant symmetry

breaking contribution in molecular conductance comes not from the details of the elec-

trode geometry but from the loss of end-to-end symmetry which is a universal feature of

8.1. The definition of the symmetry of conduction

all Σr

L(R)matrices. Consequently, I introduce a molecular-conductance point group com-

prising only the symmetry operators of the geometric point group of the molecule that do

not interchange its two ends.

This is of practical utility as it allows the dominant symmetry contributions to the trans-

mission through a molecule to be determined from knowledge of just the molecular geom-

etry alone without need for knowledge of the details of the electrode binding. The disad-

vantage is that the molecular-conductance point group may not be an exact descriptor of

the symmetry of transmission as it neglects the possible symmetry breaking contributions

arising from the internal electrode geometries. However, I will show later that in the case

of sulfur-gold binding this symmetry reduction either does not occur or is inconsequential

from a practical standpoint.

8.1.1 Application to molecular system

The approach makes it possible to separate the different symmetry contributions to the

conductance for a molecular system. To demonstrate this, calculations are performed for

conduction through a 1,4-benzenedithiol molecule chemisorbed on two gold electrodes.

Two binding configurations are considered: the bridge binding site where the sulfur atom

sits above two gold atoms and the hollow binding site where the sulfur sits above three

gold atoms. The molecule is defined as being the central C6H4S2moiety whose geometric

point group is D2h. The important difference between the two binding geometries is that

the differing symmetry in the binding gold atoms makes the conductance point group Csin

the case of the hollow binding site, whereas the bridge binding site leads to a conductance

point group of C2v. Different interface geometries thus result in different conductance point

groups. However, when the molecule is considered alone, the molecular-conductance point

group is insensitive to the interface structure: it is C2vindependent of binding site and hence

unambiguously determinable.

Figure (8.1) shows how the molecular-conductance point group is determined through

examination of the symmetry operators of the D2hpoint group of the C6H4S2fragment.

Those operators shown in red-the center of inversion, the

xy plane, and the C2xand C2y

axes are not part of the molecular-conductance point group of C2vas only the

xz and

planes and the C2zaxis may possibly commute with Σr

L(R).

It has been discussed previously that the impact of the electrodes on the molecular

electronic structure is treated more precisely, with consequent improvements in numerical

accuracy, when a number of electrode atoms are included with the molecule to form an

“extended molecule”. This inclusion has ramifications on the molecular-conductance point

group of the system, as now electrode atoms are included in the geometric point group

under consideration. If the shape of the unit cell is chosen with care, a bridge binding site

geometry will still have a molecular-conductance point group of C2v; however, a hollow

binding site geometry can only have Cssymmetry.

For the bridge binding site, the zero-voltage transmission can thus be separated into the

four different symmetry components, A1, A2, B1, and B2consisting of two

and two

representations, as shown in Fig. (8.2) for the molecule alone (black curve) and an extended

Chapter 8. The Role of Symmetry and Channels in Conduction

Figure 8.1: The symmetry operators that define the geometric point group of D2h, those shown

in green, are the operators present in the molecular-conductance point group of C2v, while those

shown in red reflect end-to-end symmetry and do not commute with Σr

L(R).

molecule with 25 gold atoms included on each side (red curve). This immediately yields

an intuitive picture of conduction in this system that is not otherwise demonstrable. At

very low energies (below -9 eV), conduction is dominated by transmission through the

system with A1symmetry; however, in the range from -9 to 4 eV that includes all accessible

energies, the dominant transmission is through the

system with B1symmetry. There is

nothing surprising in this result; it just provides a formal basis for the general assumptions

about the nature of the conductivity in conjugated systems of this kind. An interesting

feature to note is that while the region of high transmission through the

system does not

lie in an energy range that is accessible using gold electrodes whereas that for the

system

does, the transmission coefficients become large in both cases. This result accords with

what is known from intramolecular electron transfer experiments [111] where it was seen

that when systems are designed to make the

system energetically accessible then high

rates of electron transfer can be achieved over large distances.

Additionally, it is evident that the inclusion of electrode atoms in an extended molecule

simply acts as a perturbation on the symmetry separated transmissions without changing

their character. Hence including only the isolated molecule in the molecular component is

a reasonable approximation for revealing the character of symmetry separated components

of the transmission. This approximation facilitates the definition of a single conductance

symmetry for a particular molecule, irrespective of the binding site.

The bridge binding site considered above is a special case in which inclusion of the gold

atoms at the binding site does not lower the symmetry of the system. This results in the

molecular-conductance point group being the same as the conductance point group. Gener-

8.1. The definition of the symmetry of conduction

Figure 8.2: Symmetry separated components of the transmission through chemisorbed 1,4-

benzenedithiol bound to the bridge site with two gold contacts. In each case the black curve is

the transmission calculated without any gold atoms included in the extended molecule and the red

curve results from 25 gold atoms being included on each side.

ally, however, this is not the case. The simplest illustration of this is again the chemisorbed

1,4-benzenedithiol bound this time to the hollow site. In this case the symmetry of the full

system is now C2v, and although the geometric point group of the central C6H4S2moiety

is still D2h, the conductance point group (defined by the common elements of the junction-

conductance point group and the geometric point group) is no longer C2vbut Csas the C2

axis of the full system is not the electrode-molecule-electrode axis.

Figure (8.3) provides a comparison of the two ways of treating the symmetry of the

system where no gold atoms are included in the extended molecule. First are the two

symmetry components of the transmission through the molecule alone that arise, without

approximation, as the system can be described using its Csconductance point group (black

curve). These results are compared with those obtained approximately by neglecting the

terms in Σr

L(R)that break C2vsymmetry (green curve). While this approximation allows

the conductance to be classified in terms of the A1, A2, B1, and B2representations of the

C2vpoint group, the A1and B1as well as the A2and B2components may be summed and

compared to the true A0and A00 components, respectively.

Figure (8.3) shows that the results obtained using the molecular-conductance point

group are indistinguishable from the exact transmission at most energies. The important

practical result is that the symmetry reduction in systems bound through a terminal sul-

fur can be considered, to a very good approximation, to only be the loss of end-to-end

symmetry irrespective of the symmetry of the electrode-molecule junction.

Also shown in Fig. (8.3) are the results obtained using the Cspoint group with the

inclusion of 25 gold atoms per contact in the extended molecule. It is again evident that

the inclusion of gold atoms in the extended molecule perturbs the symmetry separated

Chapter 8. The Role of Symmetry and Channels in Conduction

Figure 8.3: Symmetry separated components of the transmission through chemisorbed 1,4-

benzenedithiol bound to the hollow site with two gold contacts. In each case the green curve is

the sum of the symmetry components when the conductance point group is approximated to be C2v,

while the black curve is the transmission when the system is treated without approximation as hav-

ing a conductance point group of Cs. Also shown in red is the transmission when 25 gold atoms are

included on each side to form an extended molecule.

components of the transmission without changing their character. The difference in this

instance is that with an extended molecule it is not possible to examine the hollow binding

site transmission in C2vsymmetry. Consequently, the level of detail that can be obtained in

symmetry analysis is reduced. The effect of any modification of the molecule which does

not affect all symmetry components equally, for example, chemical substitution, will be

most clearly seen in a higher symmetry representation.

The separation of transmission into symmetry components as shown in Figs. (8.2)

and (8.3) indicates that at no energy is the transmission purely through channels of just one

symmetry. There are energies at which the transmission is dominated by a single symmetry

component but other symmetry components are always present, albeit that their contribu-

tions are orders of magnitude less. This fact may be vital to interpretation of spectroscopic

properties associated with electron transmission, because channels that make small con-

tributions to the total current can, in principle, dominate effects such as IETS. Further,

a relationship between the peaks in each symmetry component of the transmission and

molecular orbitals of known symmetry can be identified leading to a useful chemical pic-

ture of molecular conductance that can aid synthetic design.

An alternative perspective on the significance of the conductance point group can be

gained by considering the nature of the wave vectors transmitted through the entire molecu-

lar component, whether this constitutes the real molecule or some extended one. In the bulk

electrode on either side, the wave vectors of the transmitted and reflected electrons can be

defined. From this point of view the molecule can be considered as a momentum filter only

8.2. The definition of channels in elastic transport

allowing certain wave vectors to be transmitted [112]. It is the conductance point group that

determines those linear combinations of wave vectors that will be allowed to pass through

the molecule from the set of all wave vectors incident on the electrode-molecule interface.

While all of the discussion so far has been concerned with systems where the scattering

processes are neglected, Σ<,>

ph =0, or when Σ<,>

ph commutes with all symmetry operators

that define the geometric point group of the molecular system, the applicability of symme-

try characterization is not limited to such situations. These conditions on Σ<,>

ph result in

it not modifying the conductance point group of the system from what would be expected

through analysis of Σr

L(R)alone. In the case where Σ<,>

ph is not equal to zero and does

not have the symmetry of the geometric point group of the system, it may act to further

reduce the conductance point group. If Σ<,>

ph has the same symmetry as Σr

L(R), this will

not be the case as no further symmetry breaking will be introduced. However, if Σ<,>

ph has

lower symmetry than Σr

L(R), it will become the dominant component in the symmetry re-

duction of the conduction process and the conductance point group will be determined also

by the symmetry of Σ<,>

ph . A symmetry breaking component arising from Σ<,>

ph may be

large, depending upon the effect being included, in which case analysis using molecular-

conductance point group may be a poor approximation to the transmission obtained using

the conductance point group.

8.2 The definition of channels in elastic transport

The idea of conduction channels is connected to the real promise of molecular electronics

which is not just in the small size of molecular components but also in tunable electrical

characteristics. If molecular properties can control conduction characteristics, then these

could be as varied as the immense range of chemical substitution allows. In the devices

developed today it is not always clear that it is the molecular characteristics, as opposed

to properties of the junction, that dominate the transmission [110,22]. However, increas-

ingly sophisticated experimental techniques are allowing more information to be obtained

about the conducting entity than just simple current-voltage characteristics. Recently, shot-

noise measurements on a single molecule were reported [113], which made it possible to

distinguish between proposed binding geometries because of the differing numbers of con-

duction channels that were active in each case. The conduction channels observed in shot-

noise measurements must have a molecular origin; however, it has not been clear how to

describe them in terms of the properties of the molecule.

Qualitatively, it is clear that resonant transport occurs at energies close to that of the

molecular orbitals of the molecular component of the system [114]. However, molecular

conduction has, so far, eluded quantitative description in chemical terms. It was shown

by B¨

uttiker that it is possible to separate transmission into components from individual

channels [115] whose transmission coefficient lies between zero and one. This has been

applied to molecular systems [116] and necked metallic wires [112], but it is not clear

how they may be quantitatively described in terms of molecular properties, for example

the molecular orbitals, of the system. In the case of single-atom metallic wires, it has been

100

Chapter 8. The Role of Symmetry and Channels in Conduction

shown both experimentally [17] and theoretically [117] that the number of channels relates

to the valence of the atom; however, this needs to be extended to describe how conduction

channels operate in molecular systems.

Here it will be shown how conduction channels can be linked to molecular properties.

First, it is shown that only a small number of conduction channels are expected because of

the nature of chemical binding of molecules to electrodes. Second, it is shown alternatively

how conduction channels can be understood in terms of the connections that they provide

between incoming and outgoing junction states.

8.2.1 B¨

uttiker channels

Using the ubiquitous Landauer-B¨

uttiker equation, the current (I) through a molecule can

be expressed as a function of bias (V) as

I(V) = 2e

hZ+∞

−∞

g(E,V)(nL−nR)dE (8.1)

The quantum conductance, g(E,V)(the sum of the transmission probabilities through all

available channels) is defined as usual:

g(E,V) = Tr[ΓL(E,V)Gr(E,V)ΓR(E,V)Ga(E,V)].(8.2)

The one-electron Hamiltonians, HD,HL, and HR, as well as the electrode-molecule

coupling matrices, JDL and JDR, are partitioned from the Hamiltonian of the full electrode-

molecule-electrode system H0

0, Eqn. (5.2). As shown by B¨

uttiker, the contributions from in-

dividual eigenchannels can be determined by diagonalizing T=ΓLGrΓRGafrom Eqn. (8.2)

to yield [115]:

g(E,V) = ∑

Ti(E,V)(8.3)

where Tiare the resulting eigenvalues and 0 ≤Ti≤1. In order to label the eigenvalues Ti

with an irreducible representation of the molecular-conductance point group of the system,

before diagonalizing T(E,V)we perform a transformation of the matrix in the symmetry

adapted basis. The matrix results block-diagonalized and every block refers to a single

irreducible representation of the point group:

T=





IR10 0

0IR20

0 0 ...



(8.4)

where IRistands for the ith irreducible representation of the symmetry group.

Performing this analysis for a 1,4-benzenedithiol illuminates a number of interesting

features. The molecule is chemisorbed at the bridge site of two gold electrodes as shown

in Fig. (8.1), with partitioning of H0

0into its HD,HL, and HRcomponents occurring across

8.2. The definition of channels in elastic transport

101

the gold-sulfur bonds. The first feature of this system is that there are very few eigen-

channels with relevant (above 10−14) transmission. Figure (8.4) shows the eigenchannel

transmissions for this system separated into the symmetry components of the molecular-

conductance point group.

Figure 8.4: Contributions to the transmission from individual conduction channels separated by

symmetry group for clarity. In each symmetry group, the primary conduction channel (dominant

transmission at the Fermi energy) is shown in red, the secondary channel is shown in blue, and only

in the A1case is a third channel visible in green.

The number of eigenchannels is formally related to the dimension of the molecule (for

the minimal basis set used in the system here we have 46 channels); however, the vast

majority of these channels do not contribute to the total transmission. The small number

of nonzero eigenchannel transmissions can be understood by considering the nature of the

chemical binding to the electrodes and the impact this has on the matrices that generate

the transmission. The main feature of the ΓL(R)(E,V)matrices is that they contain the

information about the bottleneck of the junction between the sulfur atoms and the contacts.

So, if the electrode only couples with a small number of molecular states then there will

be many near zero elements in JDL and JDR. As a result, rows and columns of near-

zero elements appear in ΓL(R)(E,V), which in turn leads to eigenchannels with negligible

transmission. This is an intuitive picture of transmission; the number of open eigenchannels

will not be greater than those supported by the smallest bottleneck in the system. It has long

been recognized to apply to necked junctions in atomic wires [117] but here it is shown to

be generally relevant to bridged-conductor conductivity.

Mathematically, this process is readily identified when the bottleneck occurs at the

sites at which HD,HL, and HRare partitioned, as described above. If the situation is

examined numerically, including 25 gold atoms on each side in the extended molecule.

This results in a modified transmission curve as the influence of the electrodes is treated

more precisely through the inclusion of explicit gold atoms in the molecular part compared

102

Chapter 8. The Role of Symmetry and Channels in Conduction

to when the effect of the electrodes is limited to the self energies, Σr

L(E,V)and Σr

R(E,V).

However, the number of eigenchannels found with transmission above 10−10 is invariant

to this enhancement. The internal coupling elements between the gold and sulfur atoms

within HDare the dominant factor controlling the number of open channels by effectively

reducing the rank of the product of g(E,V). This is an important result: it is the physical

and chemical nature of the system, not the choice of partitioning in the calculation, that

determines the number of conduction channels with significant transmission. The number

of channels with significant transmission is determined by the size of the bottleneck that

the connection through the sulfur allows and is not determined by the size of the matrix

HD, which determines the theoretical limit for the number of perceived eigenchannels.

However, through judicious choice of the partitioning, the nature of the conductivity may

be described using simple models [118,119] or molecular symmetry.

The number of eigenchannels with nonzero transmission, as well as the total transmis-

sion, is dependent on the size of the basis set used for the atoms at the bottleneck. In the

above calculations, the basis set used for the sulfur atom has s,p, and dfunctions. If the

dfunctions are removed, then the A2channel shown in Fig. (8.4) disappears completely

in addition to the changes in shape that come from the total transmission being modified.

A result similar to this was obtained in the work of Cuevas et al. [117] for tunneling

through atomic wires. Significantly, in that application the electronic coupling elements

were assumed to connect only nearest-neighbor atoms, resulting in the size of the bottle-

neck being exactly the dimension of the number of basis functions on a single metal atom.

For molecules in comprehensive treatments, it is not always reasonable to assume that only

nearest-neighbor coupling elements are nonzero. At the molecule-electrode junction, for

example, the accuracy of such an assumption will be dependent on the nature of the bind-

ing group. In the case of chemisorbed 1,4-benzenedithiol, direct Au-S and S-C couplings

are the dominant terms but an extensive calculation will also give small but nonzero Au-H

and Au-C couplings. This means that the number of open conduction channels will not be

precisely determined by the Au-S interaction.

Beyond the technicalities of the calculations, the significance of this result is that the

number of open conduction channels arises from a chemical property of the system. The

smallest bottleneck in the system can be controlled by chemical substitution, be it by multi-

ple binding groups between the molecule and the electrode or by different backbone struc-

tures for molecular wires. This understanding allows for appropriate synthetic targets to be

envisaged without any unforseen constraints in conductivity dominating the conductance

characteristics. Additionally, when the bottleneck occurs at the binding site, as was the case

for the example system above, the size of this bottleneck, that is, the number of observable

conduction channels, will change with changes in the binding geometry. For example, if

the molecule was tilted on the surface direct Au-H and Au-C coupling would be more sig-

nificant and may increase the number of conduction channels observed. Conversely, if the

terminal sulfur atom sat directly above a super-surface gold atom or at the apex of a single

atom on a gold tip, then this effect may reduce the coupling with the sulfur orbitals of

symmetry and thereby reduce the number of channels observed. This means that molecules

such as chemisorbed 1,4-benzenedithiol can be used as probes for the nature of thiol bind-

ing on gold surfaces as the number of channels observed in shot-noise measurements will

8.2. The definition of channels in elastic transport

103

be sensitive to the geometry of the junction region.

Understanding the cause of the small number of open eigenchannels does not, how-

ever, provide any further illumination of the relationship between conduction channels and

molecular orbitals. Indeed, a striking feature of Fig. (8.4) is that each conduction channel

spans a range that encompasses many molecular orbitals further complicating any quanti-

tative description of conduction through a single orbital. As has been reported elsewhere

[118], when the eigenvectors are computed in the atomic orbital basis they are difficult

to interpret. In fact, the details of the unusual character of the eigenvectors of g(E,V) are

important. The eigenvectors at all transmission energies Eare dominated by coefficients

on the sulfur atom at the charge-injecting electrode. The charge injecting electrode (L in

the notation used herein) is distinguished by the asymmetry in g(E,V), which allows the

total transmission to be calculated through the knowledge of the transmission across the

interface between the charge injecting electrode and the molecule. Where two channels

cross (akin to avoided crossings between molecular potential energy surfaces) the identity

of each channel can, in principle, be determined from the eigenvectors. However, there

are energies at which two channels meet but do not cross, a situation that cannot occur

for Hermitian operators, and the eigenvectors are not continuous over this region with the

meeting point coinciding with the discontinuity. Also, the eigenvectors are not guaranteed

to be orthogonal (zero overlap) and this point actually coincides with the maximum in the

overlap between the two eigenvectors. An example of where this unusual behavior occurs

is the meeting of two channels around -13 eV in the A1symmetry transmission. This re-

gion is shown in high resolution in Fig. (8.5) with the coefficients on the sulfur atom at the

charge-injecting electrode as a function of energy shown below for the primary and sec-

ondary conduction channels shown above, the discontinuity being readily apparent around

-13.15 eV, whereas the overlap peaks at 0.8, indicating nearly parallel eigenvectors. In

fact, the determinant of the eigenvector matrix, which is 1 for an orthogonal transforma-

tion, takes a maximum value as a function of energy of just 10−32, indicating a high degree

of eigenvector linear dependency.

8.2.2 Γ-channels

As mentioned previously, the eigenvectors of the B¨

uttiker channels reflect only the site of

charge injection and do not capture the full symmetry of through-molecule conductivity.

A means by which this symmetry can be perceived is through pre-diagonalization of the

two sets of molecule-electrode interactions ΓL(R)(E,V)as suggested by Troisi and Ratner

[38]. I ignore the orbital overlap during this process, introducing the (real diagonal) eigen-

value and (real orthogonal) eigenvector matrices ΓL(R)(E,V)and DL(R)(E,V), respectively,

through

ΓL(R)(E,V) = D†

L(R)(E,V)ΓL(R)(E,V)DL(R)(E,V)(8.5)

In order to assign symmetry to the eigenvalues a transformation of the two matrices in

the symmetry adapted basis must be performed before the diagonalization. If the two

junctions are symmetrically related, then ˆ

ΓL(E,V) = ˆ

ΓR(E,V) = ˆ

Γ(E,V)so the current

104

Chapter 8. The Role of Symmetry and Channels in Conduction

Figure 8.5: Unusual properties of the eigenvectors associated with the primary and secondary

channel of A1symmetry at the point where they meet but do not cross. The maximum in the over-

lap and the discontinuities in the coefficients on the sulfur at the charge injecting electrode atom

illustrate the unusual characteristics of the eigenvectors associated with conduction channels.

from Eqn. (8.2) can be written without approximation as

g(E,V) = Tr[ˆ

Γ(E,V)ˆ

Gr(E,V)ˆ

Γ(E,V)ˆ

Ga(E,V)]

=∑

i j

ii (E,V)ˆ

j j(E,V)|ˆ

i j(E,V)|2(8.6)

where

Gr(E,V) = D−1

L(E,V)Gr(E,V)DR(E,V)(8.7)

and

is the label for the irreducible representation to which belongs a certain eigenvalue.

This expression for the transmission suggests that the tunneling current can be viewed as

entering one electrode through the input junction channel i, being transferred through the

molecule with probability |ˆ

i j(E,V)|2and then leaving through the exit junction j. In this

8.2. The definition of channels in elastic transport

105

way, the total current can be considered to be decomposed into the sum of n2independent

channels, where nis the number of atomic orbitals of the molecule. This interpretation is

formally applicable whenever the eigenvalues ˆ

ii (E,V)are guaranteed to be all positive, as

is the case when the basis functions do not overlap. In general, the arguments used earlier

to interpret the B¨

uttiker conduction channels may again be applied, however, concluding

that most of the eigenvalues should be zero with just a few large and positive contribu-

tions facilitating current flow. Hence this interpretation, though only formally exact for

a restricted problem, is expected to be rather useful in characterizing through-molecule

conductivity.

For the case of 1,4-benzenedithiol symmetrically chemisorbed between two gold con-

tacts as studied earlier using full treatment of orbital overlap, numerical calculations reveal

only 15 junction eigenvalues ˆ

ii that are nonzero (i.e., ˆ

ii >10−16 within numerical pre-

cision); these are listed in Table (8.1). Of them, only one is negative corresponding to a

junction path of A1symmetry with ˆ

ii =−1.1×10−12. Hence, the interference contribu-

tions to Eqn. (8.6) are at least 12 orders of magnitude weaker than the direct contributions

and so the current may indeed be perceived in terms of simple contributions from each

junction linked through the molecule.

Figure 8.6: First three conduction channels obtained through the Troisi-Ratner method for each

system. In each case the primary channel is shown in green, the secondary channel is shown in

blue, and the tertiary channel is shown in yellow. Where visible, the red curve represents the total

transmission for each symmetry.

From Table (8.1), it can be seen that a significant number of the junction eigenval-

ues fall within 3 orders of magnitude of the most prolific junction channel. On the basis

of the very small number of B¨

uttiker channels depicted in Fig. (8.4) smaller number of

significant junction channels could be expected. The change is due to the neglect of the

orbital overlap in Eqn. (8.7), a feature required in order to obtain the desired |ˆ

Gr(E,V)|2

106

Chapter 8. The Role of Symmetry and Channels in Conduction

label A1A2B1B2

1 1.9462 ×10−11.2691 ×10−33.4940 ×10−28.1665 ×10−3

2 1.0196 ×10−24.4071 ×10−33.5221 ×10−4

3 1.5149 ×10−38.7315 ×10−84.7558 ×10−10

4 7.6416 ×10−59.4243 ×10−14

5 4.9025 ×10−8

6 3.6017 ×10−9

7 -1.1494 ×10−12

Table 8.1: Nonzero Eigenvalues of ˆ

ii at -5.0 eV Separated into Symmetry-Block contributions.

dependence of the molecular contribution in Eqn. (8.6). The contributions of the dominant

terms from the individual ˆ

ii (E,V)ˆ

j j(E,V)|ˆ

i j(E,V)|2terms in Eqn. (8.6) are shown in

Fig. (8.6). For each symmetry the contribution from the ˆ

11(E,V)ˆ

11(E,V)|ˆ

11(E,V)|2

term is in general the most important one. The energy dependence of these probabilities re-

flects largely the nature of the conducting molecule and is independent of the path that the

tunneling electrons takes through the entry and exit channels. The B¨

uttiker channels focus

on the conductivity bottlenecks and associated shot noise, whereas the Troisi-Ratner chan-

nels focus on the nature of both junctions and the net way in which the molecule passes

the current, but neither of these illuminates the means by which the molecule facilitates

conductivity.

8.3 Conclusions

The symmetry of molecular conductance is not simply that of the geometric point group of

the conducting molecule but instead can be described in terms of a molecular-conductance

point group and/or a full conductance point group. While the molecular-conductance point

group can always readily be determined from the geometric point group as is illustrated for

some molecules commonly used for molecular electronics applications in Fig. (8.7), the

conductance point group requires detailed knowledge of the possible symmetry breaking

contributions from the electrodes as well as any other effect included through additional

self-energies.

In the case of molecules bound to gold through a terminal sulfur, in the absence of

symmetry breaking terms in Σ<,>

ph , if the molecular-conductance point group is different

from the full conductance point group then the simpler description is shown to still provide

a good approximation to the transmission. It is thus possible to determine the symmetry

of conduction channels (and thereby understand the contribution that molecular symmetry

makes to conductance) and also understand how the selection rules of spectroscopies ap-

plied to molecules in this conducting state operate. Of course it is not always clear that

molecular properties will dominate single-molecule conduction; indeed, there is consider-

able evidence that the contacts may in some cases dominate the conduction characteristics.

8.3. Conclusions

107

Figure 8.7: The symmetry reduction between the geometric point group and the molecular-

conductance point group for some common molecules used in molecular electronics research. (*)

Precisely, the geometric point group is either C2vor C2hdepending upon the PR3rotation: The

C2hcase results in a molecular-conductance point group of Cs, whereas the C2vcase retains this

symmetry for the molecular-conductance point group.

Moreover, by separating molecular transmission in a variety of ways into non-interacting

components it is possible to isolate the roles of the junction region and the bridging molecule

in determining conductivity properties. First, the small number of rigorously defined eigen-

channels with nonzero transmission can be understood to arise from the size of the smallest

bottleneck in the system. This, in turn arises, from the chemical structure of the molecule or

the interaction of the binding group with the electrodes and consequently can be controlled

by synthetic modifications.

Further, it means that, with careful selection of target molecules, shot-noise measure-

ments can be used to probe the details of molecular binding geometry. Second, by pre-

diagonalizing Γ, conduction can be described in terms of molecular conductivity chan-

nels linking well defined junction states. Although this process is only approximate, the

interference terms found manifest in our model calculations of conductivity through 1,4-

benzenedithiol were 12 orders of magnitude smaller than the primary contributions, making

this an excellent practical method for understanding the role of the junction in the process.

108

Chapter 8. The Role of Symmetry and Channels in Conduction

Together these features make it possible to quantitatively describe molecular conduction in

terms of atomistic descriptors of the conducting system.

Chapter 9

Propensity Rules in Inelastic Electron

Tunneling Spectroscopy

This Chapter is devoted to the application of the theoretical machinery developed in Chap-

ter 8, molecular-conductance point group and channels, to investigate IETS. Already two

groups Troisi et al. [38,120] and Lorente et al. [39,40], have developed some intuitive

propensity rules that capture many of the features that lead to vibrational-mode selectivity

in IETS. The main aim of the present Chapter is to present a general approach to the inter-

pretation of IETS measurements, starting with fundamental principles implemented using

a generally applicable a priori computational scheme. Further it will be derived and justi-

fied a simplified form of the inelastic current. Finally, it will be shown how not only the

propensity rules Troisi and Ratner may be rigorously derived but also how these rules can

be extended. This can provide a thorough, yet easily, understood description of the general

phenomenon of IETS that may be readily applied in diverse practical applications.

9.1 Theoretical formalism

Starting from the ubiquitous Meir-Wingreen equation for the current:

I=2e

hZ+∞

−∞

Tr[Σ<

L(E)G>(E)−Σ>

L(E)G<(E)]dE (9.1)

and remembering that at Born Aprroximation (BA) level the electron-phonon self-energies

are

Σ<,>

ph (E) = ∑

Z∞

−∞

G<,>

0(E−E0)M

D<,>

,0(E0)dE0(9.2)

109

110

Chapter 9. Propensity Rules in IETS

with M

∂

−

∂

S−1HD−HDS−1

∂

(9.3)

for vibrational mode

with normal mode Q

. As IETS measurements are performed at

very low temperature (usually 4 K), we can simplify Eqn. (9.2) assuming temperature T

= 0. We also approximate the phonon population as a collection of Einstein oscillators in

thermodynamic equilibrium with the environment. Under this approximation, at T=0, the

population of phonon N

is set to zero and Eqn. (9.2) becomes

Σ<,>

ph (E) = M

G<,>

0(E±

ωα

(9.4)

where the upper (lower) sign is for the lesser (greater) self-energy and G<,>

0are the ze-

roth order lesser or greater Green’s functions obtained in the absence of electron-phonon

coupling. The current from Eqn. (9.1) thus becomes, expanding up to the second order in

I=2e

RµTr[ΓLGr

0ΓRGa

0]+∑

Tr[ΓLGr

0Σr

phGr

0ΓRGa

0]+Tr[ΓLGr

0ΓRGa

0Σa

phGa

0]dE¶

+∑

Θ(

L−

R−

ωα

)2e

R+

ωα

Tr[ΓLGr

0˜

ΓR˜

0]dE

=Iel +∑

Θ(

L−

R−

ωα

inel (9.5)

where the three terms in parenthesis are part of the coherent current including part of the

virtual electron-phonon scattering terms, whilst the latter term is the real inelastic current.

The inelastic current shows in fact a change in the limits of integrations due to the real

emission of a phonon.

inel is defined as

inel =2e

R+

ωα

dE (9.6)

with

=Tr[ΓLGr

0˜

ΓR˜

0].(9.7)

The use of the zeroth order Green’s functions in Eqs. (9.5) and (9.6) is a valid approx-

imation, because the electron-phonon interaction is small in the system considered here.

Numerical calculations, including the renormalization of the propagators lead to essen-

tially the same incoherent current, to within few percent of errors.

The first three terms in Eqn. (9.5) give a coherent contribution, describing the elastic

part of the current, while the last term is the inelastic component expressed as the sum of

9.2. From

-channels to A-channels

111

independent contributions, I

inel, from all of the vibrational modes of the molecule. In order

for phonon emission to take place, the applied bias must obey the condition

L−

ωα

which is emphasized by the Heaviside function in front of the inelastic component.

In the inelastic term the matrices with a tilde are computed at an energy shifted by the

energy of a phonon: E0=E−

ωα

. The shift comes from the conservation of energy and in

fact describes the lowering of the energy of the electron after the excitation of a phonon. We

hence see that the inelastic current depends on five basic quantities: the propagators Gr,a

the vibronic coupling matrix M

, the vibrational frequency

ωα

, the couplings between

the contacts and the molecule ΓL(R), and the Fermi energy EF. The main question is, of

course: what are the relationships between these quantities and how we can simplify the

picture in order to have a better insight into the physics involved ?

The application of the formalism developed in Chapter 8must certainly be applied. The

definition of the molecular conductance point group permits to assign univocally the sym-

metry to the molecule connected to the contacts and so the entire device part. This define to

which symmetry group belongs the device under investigation. After the assignment of the

molecular conductance point group is always possible to rotate all the matrices involved in

the inelastic transmission, which is the one relevant for the IETS, in the symmetry adapted

basis for the symmetry point group. This means that all the matrices are block diagonal-

ized and every block represents one irreducible representation of the molecular conduc-

tance point group. Diagonalizing the matrix we get that all the eigenvalues spectrum can

finally be assigned to a certain irreducible representation as the corresponding eigenvector.

Particular important is the transformation into the symmetry adapted basis of the electron-

phonon coupling matrices M

that, depending on the symmetry of the vibrational mode

they describe, can result in a block diagonal matrix in case of a totally symmetric mode or

in a matrix with off diagonal blocks which describe the mixing between channels which

belong to different irreducible representations induced by the vibrational mode. Since now

it is assumed that all the matrices are already rotated in the symmetry adapted basis corre-

sponding to the molecular conductance point group so that it is always possible to assign an

irreducible representation label to all the modes and the channels involved in the analysis

of the inelastic current at the base of the IETS signal.

9.2 From Γ-channels to A-channels

The specific system investigated in this study is the 1,4-benzenedithiol molecule chemi-

sorbed between two gold contacts already used in the previous Chapter, though the ap-

plications of the method presented are quite general. As already demonstrated, the sym-

metry properties of coherent transport through gold-thiol junctions are dominated by the

molecular symmetry and that low-dimensional models of the electrode capture most of the

essential features of observed IETS. So only two atoms are used to represent each lead, as

shown in Fig. (9.1).

The low-temperature IETS spectrum calculated using only the Born approximation

is shown and assigned in Fig. (9.2). We seek a basic understanding of how this spec-

trum arises and the relative propensities calculated for the modes. The chemisorbed 1,4-

112

Chapter 9. Propensity Rules in IETS

Figure 9.1: A chemisorbed 1,4-benzenedithiol molecule connected to two gold atoms on each side.

benzenedithiol molecule has d2hsymmetry, and the IETS-active modes are categorized in

Fig. (9.2) accordingly: the most prominent modes are in-plane totally symmetric agmodes

while the out-of-plane antisymmetric b2gand b3umodes are also significant. However,

the conduction process depicted by Eqn. (9.6) does not display d2hsymmetry as the ΓL(R)

matrices have only the symmetry of the Left and Right junctions, not the full molecular

or device symmetry, see Chapter 8. Analysis of Eqn. (9.6) thus must commence with this

reduced symmetry, the molecular conductance point group symmetry which, in this case,

is C2v; note that lower-case symbols are used throughout for the description of molecular

symmetry while upper-case symbols are used for the description of molecular-conductance

symmetry. We must demonstrate how IETS appears to take on the molecular symmetry

properties in spite of this limitation.

0 300 600 900 1200 1500 1800 2100 2400 2700 3000

Frequency (cm-1)

-1×10-3

1×10-3

2×10-3

3×10-3

4×10-3

d2 I / dV2 (A/V2)

0 0.1 0.2 0.3

Energy (eV)

ν30

b3u

ν29

b2g

ν26

ν24

b3u

ν20

agν19

b1g

ν15

ν9

ν8

Figure 9.2: The IET spectrum of our system. The broadening of the peaks has been introduced

empirically considering a phonon lifetime of 6.6 ×10−13 s, corresponding to a broadening of 2 ×

10−3eV.

The approach is based on the idea that both the elastic and inelastic current can be

expressed as the sum of a small number of essentially non interacting paths or channels

through the device. Using the same machinery used for elastic transport through the intro-

9.2. From

-channels to A-channels

113

duction of the transformation

ΓL(R)=D†

L(R)ΓL(R)DL(R),(9.8)

which reduces the electrode-molecule coupling matrices to diagonal form, a transformation

that captures the essence of the physical insight used previously by Troisi et al. [38,120]

in their proposed IETS propensity rules. After this transformation is applied we get the

following formula for the inelastic current of every mode:

inel =2e

R+

ωα

Ã∑

i j

ΓL,

ii ˆ

ΓR,

j j (E−

ωα

)|ˆ

Λi j|2!dE,(9.9)

where ˆ

Λ=ˆ

0ˆ

0(E−

ωα

)and

and

are the irreducible representation labels of the

left and right channels respectively.

The above junction channels will provide a good simple description of the (elastic or in-

elastic) conduction process whenever the bottleneck between the molecule and its contacts

is the most important physical element. In this depiction, |ˆ

Λi j|2determines the probability

that an electron or hole that enters the molecule from Left lead through channel i, with

symmetry

, is scattered inelastically out the Right lead through channel j, with symmetry

. In Fig. (9.3), upper section, it is plotted the amplitude ˆ

ΓL,

ii for an electron entering the

molecule in a window of energy of 5 eV encompassing the Fermi energy. The eigenvalues

show a very low dependence on energy as the s-band density of states of the gold con-

tacts is nearly energy independent. The graph shows also that effectively only one junction

channel, of symmetry A1in the molecular-conductance point group, dominates the process,

with the next most significant channel being of B1symmetry but an order of magnitude less

transmissive. This result suggests that the most intense IETS process is likely to involve

electrons or holes entering and exiting through the dominant A1Left and Right junction

channels, respectively; a process that is only possible when vibrations of A1symmetry

are involved. While this argument correctly predicts the propensity for agmodes apparent

in Fig. (9.2), the dominant process is actually found to involve the less transmissive B1

channels of the Left and Right junctions. This shows that the ˆ

Λmatrix, which contains

information about the chemical properties of the molecule, also plays an important role

in selecting which channel is the most important for the current, a role that arises as the

conduction in the molecule near the Fermi energy of chemisorbed 1,4-benzenedithiol is

dominated by the electronic

-system that embodies B1but not A1symmetry.

The promise of this approach points to the possibility that a new set of channels might

preserve all the nice characteristics of the previous ones, but allow for better insight into

the role of the molecular properties. My proposal lies in a second transformation based on

diagonalization of the matrices AR=Gr

0ΓRGa

0and (AL)∗=Ga

0ΓLGr

AR=C†

RARCR,

AL=C†

L(AL)∗CL.(9.10)

114

Chapter 9. Propensity Rules in IETS

-7 -6 -5 -4 -3 -2

10-4

10-3

10-2

10-1

100

Γ eigenvalues (a.u.)

-7 -6 -5 -4 -3 -2

Energy (eV)

10-1

100

101

102

103

104

A eigenvalues (a.u.)

Figure 9.3: The energy dependence around the Fermi energy EFfor the eigenvalues ˆ

ΓL,

ii of ΓL

(top) and AL,

ii of AL(bottom) partitioned into symmetry components.

The conjugation of ALarises because the Hamiltonian is a real matrix and therefore (Gr

0)†=

(Gr

0)∗=Ga

0. Moreover, due to the fact that ALis a Hermitian positive definite matrix, the

conjugation does not change its positive and real eigenvalues.

Although in this derivation I have used the unrenormalized propagators, the approach

is quite general since Eqn. (9.6) is valid also in the so called self-consistent Born approx-

imation (SCBA), provided all Gr,a

0are substituted by renormalized propagators Gr,a. The

electron-phonon self-energy, Σr,a

ph, do not reduce further the symmetry of (AL)∗and AR.

These matrices depict the coupling-weighted molecular density of states [53] and con-

tain information about not only the junctions but also the chemical properties of the molecule

as well. In Fig. (9.3), bottom section, there are plotted the eigenvalues AL,

ii to be compared

to those previously discussed for ˆ

ΓL,

ii . The first thing we can observe is that a strong en-

ergy dependence appears with the AL,

ii eigenvalues showing peaks close to the energies

of the molecular orbitals. The dominant eigenvalues at the Fermi energy are listed in Ta-

ble (9.1) and, as expected, by far the largest eigenvalue is found to be of B1symmetry,

named 1B1. Hence this new transformation provides improved insight into the physical

problem of IETS scattering. The inelastic conductances are then given by

=∑

i j

AL,

ii AR,

j j (E−

ωα

i jM

ji,(9.11)

where

=C†

CRand M

=C†

CL.(9.12)

9.2. From

-channels to A-channels

115

Differently from the Γ-transformation the Amatrices are positive definite and so their

rotation should preserve this property. A subtlety, however, is that due to the differing

energy dependencies of ΓLand ΓRapparent in Eqn. (9.5), the eigenvectors CLand CRare

not simply related to each other so that M

i j and M

ji become unrelated complex quantities

and hence the i j contributions in Eqn. (9.11) cannot strictly be interpreted as independent

channels.

Nevertheless, in Table (9.2) the major contributions to g

evaluated at the Fermi energy

from the double sum of Eqn. 9.11 are listed and at most two contributions account for at

least 95% of the IETS signal. Single significant contributions are found only for the case

of totally symmetric vibrations (agsymmetry in the full molecular point group), these

involving an electron entering through the 1B1channel of the Left lead, scattering off an ag

vibration and exiting through the 1B1channel of the Right lead. Otherwise, closely related

pairs of contributions are involved, with for example the intense b2gmode

29 activated by

a 50% contribution from an electron that enters the 1A1channel and exits the 1B1channel

combined with a 45% contribution from an electron that enters the 1B1channel and exits

the 1A1channel.

label AL,

ii (EF)

1A13.44

2A10.144

3A10.0128

1A20.0107

1B127.4

2B10.0705

1B22.26

Table 9.1: All significantly large eigenvalues in atomic units of the electrode-coupling weighted

density of states, AL,

ii at the Fermi energy EFfrom Fig. (9.3).

The analysis presented in Table (9.2) presumes that the total IETS signal I

inel from

Eqn. (9.5) can be approximated from the inelastic scattering conductance g

evaluated

only at the Fermi energy EF. To test this hypothesis, the relative values of the total inelas-

tic currents calculated using Eqn. (9.5) and its approximation are given in the table where

they are seen to be in very good agreement. The lack of symmetry in the leading terms

shown in the table arise due to the shift

ωα

that must be applied to the outgoing chan-

nel energies. Owing to the large values of the active frequencies compared to the energy

dependencies depicted in Fig. (9.4), such differences appear profound and indeed result

in factors of two differences between the displayed eigenvalues AL(EF)and the analogous

values of AR(EF−

ωα

). However, Troisi [38,120] has argued that such changes should not

qualitatively affect IETS propensity rules and hence in Table (9.2) an approximate analysis

of g

based on this further approximation is presented. Indeed, the relative propensities

evaluated using this crude approximation differ by at most 20% from the exact intensities

inel, indicating its usefulness. Further, this approximation leads to the re-expression of

116

Chapter 9. Propensity Rules in IETS

Eqn. (9.11) in terms of true channels as

=∑

i j

AL,

ii AR,

j j |M

i j|2(9.13)

and the dominant terms in these sums, now symmetrically disposed towards both leads, are

also shown in Table (9.2). Further, this approximation facilitates the rewriting of Eqn. (9.6)

in the simple form

Tr[ΓLGr

0ΓRGa

0] = Tr[ΓL

∂

QΓR

∂

Q](9.14)

using Gr,a

Gr,a

∂

Gr,a

∂

Q. This result is precisely that obtained by Troisi [38,120] using

a perturbation expansion of the elastic component of the current. Equation (9.14) puts this

very useful expression in the context of a general theory for the non-equilibrium transport

process and a hierarchy of numerical methods for obtaining the exact and approximate

IETS intensities.

9.2. From

-channels to A-channels

117

Vibration

inel g

(EF)exact g

(EF)approximate

ν ω

(eV [cm−1]) d2h(C2v)Total Total channel % channel % Total channel % channel %

30 0.014 (112) b3u(B1) 0.269 0.209 1A1→1B148 1B1→1A146 0.288 1A1→1B148 1B1→1A148

29 0.028 (228) b2g(B1) 0.189 0.173 1A1→1B150 1B1→1A145 0.229 1A1→1B148 1B1→1A148

26 0.043 (346) ag(A1) 0.199 0.175 1B1→1B197.9 0.225 1B1→1B198

25 0.045 (365) au(A2) 0.032 0.027 1B1→1B249 1B2→1B151 0.035 1B1→1B250 1B2→1B150

24 0.057 (462) b3u(B1) 0.023 0.024 1A1→1B157 1B1→1A142 0.029 1A1→1B150 1B1→1A150

23 0.072 (583) b1u(A1) 0.002 0.002 1B1→2B143 2B1→1B157 0.002 1B1→2B150 2B1→1B150

21 0.088 (713) b2g(B1) 0.014 0.008 1A1→1B154 1B1→1A146 0.01 1A1→1B150 1B1→1A150

20 0.098 (790) ag(A1) 0.127 0.103 1B1→1B197 0.121 1B1→1B196

19 0.102 (826) b1g(A2) 0.01 0.008 1B1→1B248 1B2→1B152 0.009 1B1→1B250 1B2→1B150

18 0.104 (836) b3u(B1) 0.008 0.004 1A1→1B155 1B1→1A145 0.005 1A1→1B150 1B1→1A150

17 0.125 (1009) b2g(B1) 0.004 0.004 1A1→1B162 1B1→1A138 0.004 1A1→1B150 1B1→1A150

15 0.138 (1114) ag(A1) 1.000 1.000 1B1→1B1100 1.000 1B1→1B198.7

13 0.161 (1303) b1u(A1) 0.004 0.004 1B1→2B149 2B1→1B151 0.003 1B1→2B150 2B1→1B150

9 0.205 (1656) ag(A1) 0.084 0.084 1B1→1B1100 0.073 1B1→1B199.6

8 0.212 (1712) ag(A1) 0.531 0.509 1B1→1B1100 0.452 1B1→1B199

Table 9.2: Properties of some significant IETS active vibrational modes

of frequency

and symmetry as specified in the molecular point group

(conductance point group), including the relative total intensity I

inel from Eqn. (9.5), the relative zero-voltage conductance g

from Eqn. (9.6)

evaluated at the Fermi energy EF, and the identity and contribution of its most significant channel contributor(s), and that as obtained approximately

through the replacement of AR(E−

ωα

)with AR(E).

118

Chapter 9. Propensity Rules in IETS

9.3 Propensity rules in IETS

Based on this simple formalism, the factors that lead to the propensity of different modes in

IETS are depicted graphically in Figure (9.4). The top part is a schematic plot of six charac-

teristic normal modes and their frequencies. The second and third rows show the modulus

of the eigenvectors of (AL)∗and ARassociated with the dominant channels through the

coupling-weighted molecular density that indicate how transported charges enter and leave

the molecular region. These Left and Right pairs of channels are coupled by a molecular

vibronic coupling term |M

i j|2whose origin in terms of interfering atomic contributions and

bond contributions is shown in the lower part of the figure.

Figure 9.4: Diagrammatic description of the origin of the IETS intensity for six characteristic vi-

brational modes. First row: depictions of the normal modes, including either arrows for in-plane

modes or open and closed circles for out-of-plane modes, along with the mode frequency in cm−1

and molecular symmetry. Second and third rows: circles indicating the absolute values of the

atomic contributions to the dominant eigenvectors CL

iof ALand CR

jof ARalong with the associated

eigenvalues AL,

ii and AR,

j j , in a.u., from Table (9.1). Fourth row: origin of the vibronic term M

i j

that couples these channels expressed in terms of atomic contributions (circles) and bond contribu-

tions (lines) colored blue and red to indicate constructively and destructively interfering processes,

respectively; also indicated is the total inelastic conductance at EFof the mode and the percentage

contribution arising from the indicated coupled channels (doubled for non-totally-symmetric modes

to account also for its symmetry-related counterpart), and the destructive interference indicator

To do this we can expand

i j|=Re(∑

E∑

i j

EF e−i

),(9.15)

9.3. Propensity rules in IETS

119

where Eand Frefer to different atoms,

is the phase of the complex number M

i j, and the

atomic (diagonal) and bond (off-diagonal) contributions

are defined using:

i j

EF =∑

m∈E∑

n∈F

(CL

im)∗M

mnCR

n j,(9.16)

where CL(R)

i j are the elements of the eigenvectors and M

mn are the vibronic coupling matrix

elements in the atomic-orbital basis. The off-diagonal elements of this coupling matrix

represent scattering from the bonds of the molecule while the diagonal elements represent

scattering off the atoms. Through chemical or other external modifications to the molecule,

Fig. (9.4) indicates how the IETS propensities of different modes can be manipulated.

To show how the lower frames of Fig. (9.4) depict the molecular properties, the dom-

inant contributions

i j

EF evaluated for the most intense IETS mode,

15, are given in Ta-

ble (9.3). The largest contribution (40 % of the total) arises from each of the two C-S

bonds. The upper frames in Fig. (9.4) indicate why a large contribution arises from the C-S

bonds: both the L and R channel eigenvectors contain components on these atoms, and the

vibrational mode acts to change the C-S bond length. In the lower frame of Fig. (9.4), the

thick and dark blue lines over the C-S bonds indicate this contribution graphically. From

Table (9.3), the next most important contribution is 19 % and arises from each of the C-Co

bonds, where Cois the carbon ortho to the linked carbon. In Fig. (9.4), the total intensity

of blue color for the C-S and C-Cobonds indicates their relative importance. Blue coloring

is used to indicate that the contributions from the C-S and C-Cobonds add constructively

to enhance the inelastic current, their net effect being to generate 2×40 % + 4×19 % =

158 % of the current. This exceeds 100 % as the remaining minor contributions to the cur-

rent, dominated as indicated in Table (9.3) by the atomic S (-8 % each) and C (-4 % each)

contributions, act destructively to reduce the current. Contributions that act destructively

are indicated in red in Fig. (9.4). The total color density shown in circles for the atomic

contributions is scaled to that of the bond vectors.

S C CoCo0

S -0.08

C 0.40 -0.03

Co-0.01 0.19 -004

Co0-0.01 0.02 -0.02 -0.04

Table 9.3: Normalized atomic (

i j

EE /ΣEF |

i j

EF |) and bond ( (

i j

EF +

i j

FE )/ΣEF |

i j

EF |) contributions

to the molecular vibrational couplings between input and output electron scattering channels, |M

i j|

from Eqn. (9.15), as visualized in Fig. (9.4); contributions are shown for the sulfur (S), connected

carbon (C), its ortho carbon (Co), and its symmetry-related neighbor (Co0).

For each mode depicted in Fig. (9.4), the impact of destructive interference between

electron-scattering pathways is quantified using the indicator

=1−|∑EF

i j

EF |

∑EF |

i j

EF |,(9.17)

120

Chapter 9. Propensity Rules in IETS

which takes on the value of

=0 to indicate only constructive interference to

≈1, in-

dicating strong destructive interference. Modes dominated by destructive interference are

likely to be highly sensitive to external modifications through modulation of the interfer-

ence. For the example of

15 considered previously, bond scattering is opposed by atomic

scattering so that

is quite large at 0.56. Of all of the vibrational modes, this mode, the

most intense mode, has one of the lowest values of

, however, indicating that interference

effects will in general limit the development of simple chemical models for IETS intensity.

From the analysis of the channels propensity rules can also be defined. The most active

modes are the agmodes. Based on the notion that the most intense modes will be those that

access the dominant 1B1paths through each of the junction-weighted densities of states,

the most active modes are expected to be of A1symmetry in the C2vconductance point

group. Such modes will have either agor b1usymmetry in the full molecular point group,

d2h. Indeed, agmodes are found to be the most active ones, but b1umodes are found to be

very weak. To understand this differentiation, we note that the transmission eigenvectors

CL,R

1B1associated with the Left and Right 1B1channels can be represented as a sum of terms

each with either agor b1usymmetry:

CL,R

1B1=1

√2(ψag±ψb1u)(9.18)

where the upper (lower) sign is for the left (right) eigenvector. Substituting Eqn. (9.18) into

Eqn. (9.12) splits the sum into four terms:

i j =1

2(ψag)†M

ψag−1

2(ψb1u)†M

ψb1u+(9.19)

2(ψag)†M

ψb1u−1

2(ψb1u)†M

ψag.

The product of the three elements in Eqn. (9.19) must be totally symmetric so only the

first two terms may be non zero for agmodes while only the last two terms are permissible

for b1umodes. For agmodes, the two allowed terms differ fundamentally in nature from

each other facilitating an allowed net contribution. However, for the b1umodes, the two

non-zero contributions exactly cancel each other, preventing inelastic scattering involving

the same input and output channels. Hence the inelastic scattering for the b1umode

shown in Fig. (9.4) involves two different B1channels, 1B1and 2B1, so that the b1umodes

thus behave in the same fashion as do all other non-totally-symmetric vibrations. Concep-

tual approaches that exploit the sparceness properties of ΓLand ΓRalone do lead to the

primary propensity rule favoring totally symmetric IETS excitations.

The propensity rules derived by transforming the Γs parallel these propensity rules

and can be thought of as arising through similar arguments and serve to identify the most

active modes in IETS. However, through the complete description of the junction-weighted

densities of states, it is possible to derive also propensity rules for the next most active

vibrations: these will be the ones associated with both the largest and the second largest

eigenvalues AL,

ii , in this case the 1B1and 1A1channels, respectively, which are coupled by

9.4. Conclusions

121

modes of b3uand b2gsymmetry. Indeed, modes of both of these symmetries are identified

as being quite prolific in Fig. (9.2) and in Table (9.2). Of note is the fact that these modes

are out-of-plane modes that do not have a components in the direction of the charge flow;

rather, they scatter electrons and holes between

and

channels.

Within each particular symmetry class, the most active modes are seen from Fig. (9.4)

to be those that involve atomic motion on atoms that have large amplitudes in the coupling-

weighted Left and Right channels. As a result, a wide variety of scattering paths through the

molecule can be invoked. As described earlier, the most intense agmode

15 from Fig. (9.4)

is an in-plane ring deformation mode that embodies some C-S stretching character, and all

parts of the molecule contribute to the scattering. Also shown in Fig. (9.4) is the nature and

scattering origin for the next most prolific mode, the agC-H bending mode

8. This mode

is coupled to some Co-Co0stretching character, and it is this that facilitates the inelastic

scattering from the input 1B1channel to the output 1B1channel.

For the intense non-totally symmetric modes that scatter charges between

and

channels, the most active modes are found to be those with significant C and S involvement.

The scattering is not generated symmetrically from each end of the molecule, as scattering

between input and output channels of the same type is constrained to be, but instead is

typically dominated by a single C-S bond. Typifying this behavior is

30, the b3umode

that couples the two most conductive junction channels 1A1and 1B1(see Table (9.1)).

Figure (9.4) shows that the input and output channels both have amplitude on the one C-S

bond while the out-of-plane vibration also has C-S amplitude, allowing the

mixing to

occur. This mode is indeed the third most active mode calculated for IETS (see Table (9.2)),

but, as the bottom frame of Fig. (9.4) indicates, the dominant scattering from the C-S bond

is strongly opposed by scattering from the individual C and S atoms. Were it not for this

interference, this vibrations would be even more prolific in the calculated IETS. The largest

destructive interference found was for the related b3umode

24 for which

= 0.91.

Finally, I consider the C-H out-of-plane modes

19 (b1g) and

25 (au) that have the

same form except for opposite end-to-end symmetry. These modes coupe the most prolific

and third-most prolific molecule-weighted junction channels 1B1and 1B2, but as can be

seen from the forms of the channels shown in Fig. (9.4), very few atoms are active in

both channels. Hence the scattering is intrinsically weak and contains contributions from

non intuitive non bonded 1,3- and 1,4- intermolecular interactions. While these modes

are predominantly C-H in character, they do provide scattering through their weak C-C

contamination.

9.4 Conclusions

In conclusion, two approaches has been investigated to analyze the IETS signal in a 1,4-

benzenedithiol molecule chemisorbed between two gold leads in the context of a full gen-

eral formalism for non-equilibrium elastic and inelastic conduction processes. In both ap-

proaches, the inelastic current is split into a small number of non interfering contributions

or channels. The first approach, based on the insight provided by Troisi et al. [38,120]

that only a few paths through the junctions are accessible, involves transformation of the

122

Chapter 9. Propensity Rules in IETS

molecule-electrode couplings ΓL(R)and reduces dramatically the complexity of the prob-

lem. However, detailed insight into the influence of the molecule on these paths is required

before the method can be put to practical use in the determination of propensity rules.

The second approach provides this insight automatically through the transformation of the

(AL)∗and ARmatrices that depicts the density of states of the molecule coupled to the

contacts.

The A-channel transformation allows the dominant channels for electron or hole con-

duction through the junctions, channels of different symmetries, to be identified, leading to

propensity rules based on the affect of the normal modes of vibration in scattering charges

between these channels. As for molecules such as 1,4-benzenedithiol with dominant

conduction character, the molecule-weighted junction channel 1B1is very much more

prolific than in any other channel. This leads to the first propensity rule that the totally

symmetric modes (ag) dominate IETS as only these can couple 1B1from the Left lead to

1B1in the Right lead. A similar propensity rule is expected for all molecules chemisorbed

to gold through sulfur links, independent of the actual identity of the dominant channel.

Weaker contributions to IETS are then identified from the nature of the next-most signif-

icant molecule-weighted junction channel eigenvalue(s), which for 1,4-benzenedithiol are

found to be 1A1and 1B2, leading to a propensity for b2gand b3umodes and b1gand au

modes, respectively, in its IETS.

Once these dominant channels are identified for modes of a particular symmetry type,

the most active modes can be determined by examination of how the normal mode affects

the atoms accessed by the appropriate molecule-weighted junction channels and the na-

ture of the junction channels: the same atoms must be involved in each channel, and the

vibration must perturb these atoms. For coupling between the 1B1and 1B2channels, few

atoms are involved in both channels and so the IETS is weak, but the overlap between the

1B1and 1A1channels is large and hence the IETS is strong. These intense IETS modes

involve the mixing of molecular

and

character through out-of-plane C-S distortions.

However, a more subtle feature acts to determine the final IETS intensities: the scattering

amplitude can be decomposed in terms of interfering contributions associated with scat-

tering from each atomic centre and from each bond in the molecule, and this interference

is in general large. A practical consequence of this is that chemical and other variations

are likely to modulate the IETS intensity associated with particular modes of vibration.

The analysis presented thus provides simple and effective a priori means by which a very

complex process involving no formal selection rules can be controlled and manipulated to

achieve desired outcomes.

Conclusion

My work was focused in the investigation of some of the effects induced by the vibrations

of the molecule to the conductivity in molecular electronic devices. As already mentioned,

this work is particular important in order to investigate the stability of molecular devices,

evaluating the dissipated power, and to simulate new techniques to probe the molecule, like

IETS.

I applied the NEGF machinery up to the first perturbation term to get a consistent

set of equations and evaluate the current including the vibrational effects. The code was

successful in handling systems at low temperature strongly bounded to the contacts. This is

the typical setup of many experiments using molecule on gold surfaces covalently bonded

via thiol groups, as in break-junction or self-assembled monolayer experiments.

The application of the code to alkanethiols, Chapter 6, has shown how the vibrational

degrees of freedom can play an important role in the final dissipation of the molecule. In

particular it has shown the importance in electron scattering of vibrational modes which

are localized in the junction region between the molecule and the contacts. This means

that to catch the real dissipation effects a good characterization of the surface cannot be

neglected. This demonstrate that if the molecule is characterizing the shape of the current-

voltage characteristic, the surface is fundamental to get the right magnitude of the current

and of the power dissipated, see also Di Ventra et. al. [121].

The same results were obtained in the IETS in alkanethiols, were the effect of sulphur-

gold interaction was the origin of one of the most intense peak in the spectrum. The IETS

has provided to be one of the most interesting topic for the application of the electron-

phonon code. In fact, it is a very flexible technique to investigate many aspects of molec-

ular devices which normally are beyond the possibility of more conventional spectroscopy

techniques in molecular electronics. In particular I used the simulation of IETS to in-

vestigate how the geometry configuration of the molecule between the metallic electrodes

and how the junction affects the spectrum. Particular important is this sensitivity to the

metal-molecule interface if molecular electronics shall move in the direction of high repro-

duceable devices for industry.

Moreover, from my simulation it is evident that the scattering of electrons in the molecule

is a complex interplay between the density of electrons, the external bias applied and the

localization of the most relevant vibrational modes. This suggests the possibility to decom-

123

124

Conclusion

pose the current into independent channels. Those channels can be labeled following the

reduced symmetry point group of the molecule between the electrodes. In case the junc-

tion reduces to a simple thiol group which bonds the molecule to the metal, the number of

dominant channels reduces drastically. The irreducible representations of those channels,

coupled with the symmetry of the modes, has permitted to cast a light to the complicate

propensity rules of IETS. Rules which depend on both the molecule and the interfaces with

the contacts.

The computational efficiency of DFTB and the flexibility of the code have permitted

to implement new features very easily during the proceeding of the PhD. In particular the

introduction of symmetry was fundamental to the investigation of selection rules in IETS.

This represents probably the largest strength of the method. However, other issues need to

be investigated:

•The strength of DFTB relies on the possibility to use a systematic method of parametriza-

tion for the zero order Hamiltonian. However apart for organic compounds still a re-

liable parametrization for organic elements is missing with metal, in particular with

gold.

•The inclusion of the charge polarization induced by the bias is fundamental for higher

voltages. The merging of the electron-phonon code with the Poisson solving gDFTB

is thus a mandatory future step.

•The correlation between electrons must be improved beyond the mean field approx-

imation implemented in DFTB. Many studies [122] show that the charge correlation

plays an important role in the alignment of the HOMO-LUMO gap and the Fermi

energy of the contacts, fundamental aspect for every quantitative evaluation of the

current. At the moment a scheme to include a quasiparticle correction, at the so

called GW level, is under development [123,124].

Concluding, this work was motivated by the large impact that molecular electronics has

mostly in the scientific community. Despite many years of research, the field has still

not reached a degree of maturity that permits the transition from the laboratories to the

industry. However, the idea of using a single molecule to characterize the entire device

is revolutionary for the impact it can have to the concept of electronics and in the entire

engineering field. It is certainly a step to get closer to the workings of nature, considering

that all the lifeforms on this planet are nothing else that extremely complex self-assembled

systems made of molecules.

The possibilities are very large, but the theoretical, as the experimental, work is ex-

tremely challenging. Already the definition of the system is rather complicate considering

that in the same system we have to face and describe at the same level of accuracy a bulk

region (the contacts), a surface (the leads) and a molecular region. This puts immediately

some constrains to the kind of theoretical methods that can be used, considering that ab-

initio quantum chemistry methods are at the moment unable to handle such large systems.

Since the first seminal paper about molecular electronics [125], many issues have been

addressed from the theoretical point of view, but many others need to be fully investigated.

Conclusion

125

Main fundamental points are still under debate: is a molecule mostly a wire or more a quan-

tum dot? Which is the effect of the contacts and especially the surface on the conductivity?

Which role do static and dynamic correlation play in conduction when considering also

effects such as Coulomb blockade? How can the spin-spin interaction change the entire

picture?

It is clear that only a full description of all these issues together can really address and

solve the problem of theoretical description of molecular electronic conductivity.

At the moment, the best compromise between all these issues remains the Green’s

function approach using DFT as starting point. However, it is now clear that the method

is unable to completely handle many fundamental characteristics of the molecular conduc-

tivity. There are basically two ways to go beyond a DFT description: changing directly the

DFT method as starting point for the Hamiltonian or refining the calculation with correc-

tions at the Green’s function level. This second approach has been followed in our group

in the gDFTB codes as, i. e., by the electron-electron correlation correction with GW.

However, the author suggests that in the future a more radical change in the approach

will be needed to really overcome the problems that cannot be fully described at the

DFT level of theory, many candidates for hybrid methodologies already exist, like time-

dependent DFT (TDDFT).

Appendix A

Atomic Units

The atomic units (a.u.) are a set of units specifically designed to fulfil two specific purposes:

first, present the results from atomistic calculations in digestible numbers and, second,

decouple the numerical calculations from specific values of fundamental constants. To see

how these units arise naturally let us consider the Schr¨

odinger equation for the hydrogen

atom. In SI (Systeme International d’Unites) units, we have

·−¯

2me

∇−e2

πε

0r¸

εφ

(A.1)

where ¯

his the Plank’s constant divided by 2

,meis the mass of the electron. To cast this

equation into dimensionless form we let x,yand z→

x0,

y0and

z0and obtain

·−¯

2me

2∇0−e2

πε

r0¸

εφ

0(A.2)

The constant in front of the kinetic and potential energy operators can then be factored,

provided we choose

such that

2=e2

πε

a(A.3)

where

ais the atomic unit energy called Hartree. Solving Eqn. (A.2) for

we find

πε

0¯

mee2=a0(A.4)

Thus

is just the Bohr radius a0which is the atomic unit of length called a Bohr. Finally,

since

a·−1

2∇0−1

r0¸

εφ

0(A.5)

127

128

Appendix A. Atomic Units

if we let

a, we obtain the dimensionless equation

·−1

2∇0−1

r0¸

0(A.6)

which is the Schr¨

odinger equation in atomic units.

Physical quantity Conversion factor X Value of X (SI)

Length a05.2918 ×10−11 m−1

Mass me9.1095 ×10−31 Kg

Charge e 1.6022 ×10−19 C

Energy

a4.3598 ×10−18 J

Action ¯

h1.0546 ×10−34 Js

Table A.1: Conversion of atomic units to SI units.

The solution of this equation for the ground state of the hydrogen atom yields an energy

0equal to -0.5 atomic units. Table (A.1) gives the conversion factor Xbetween atomic

units and SI units, such that the SI value of any quantity Qis related to its value in atomic

units Q0by

Q=XQ0.(A.7)

Conversion factors for a few other units, which are not related to SI but which are necessary

to read the existing literature, are so follows. One atomic unit of length equals 0.52918

Angstr¨

oms ( ˚

A) and one atomic unit of energy equals 27.211 electron volts (eV) or 627.51

kcal/mole.

Appendix B

Surface Green’s Function: Decimation

Technique

In this appendix the algorithm to evaluate the contacts self-energies is discussed. The

evaluation of the self-energies, from Eqn. (5.11) Chapter 5, is related to the computation of

the so called surface Green’s functions gr

Land gr

Rof the contacts. These Green’s functions

are the propagators of the two contacts assuming that they are disconnected from the bridge.

To evaluate them it is generally used a decimation algorithm, following the Guinea et al.

[126] scheme.

Figure B.1: The contact can be seen as a semi-infinite repetition of a fundamental supercell called

principal layer (PL). The size of the principal layer is chosen large enough in the way that only

adjacent PLs interact via the

blocks.

129

130

Appendix B. Surface Green’s Function: Decimation Technique

In Fig. (B.1) it is shown schematically the shape of one contact. We can define them

as the repetition of an infinite number of connected subsystem called principal layers (PL).

The PLs are defined in such a way that only nearest neighbor interact. The Dyson’s equa-

tion for the infinite system can be written in terms of block matrices as







S00 −H00

S01 −τ01 0 0 ···

S10 −τ10

S11 −H11

S12 −τ12 0

S21 −τ21

S22 −H22

S23 −τ23

0 0

S32 −τ32

S33 −H33

....







·





00 gr

01 gr

02 gr

02 ···

10 gr

11 gr

12 gr

20 gr

21 gr

22 gr

30 gr

31 gr

32 gr

....







=1(B.1)

where

= (E+i

)with

a positive infinitesimal quantity (in the calculations was used

10−6). The subindexes label the different PLs and the block gmn is the propagator connect-

ing the mth and nth layer. For an ideal system we can assume that the blocks of the PLs are

all equal (H00 =···=Hmm =HPL) as the coupling terms (τ10 =···=τnm =τ). Moreover

the coupling terms are related by the following relation: τmn =τ†

nm. The same relations

are for the overlap matrix S. We can separate the Dyson’s equations connecting the surface

layer, or 0th layer, with the rest, obtaining a set of linear equations

(

SPL −HPL)gr

00 −(

−τ)gr

10 =1

(

SPL −HPL)gr

10 −(

−τ)gr

20 −(

S†

−τ†)gr

00 =0

(

SPL −HPL)gr

20 −(

−τ)gr

30 −(

S†

−τ†)gr

10 =0

···

(

SPL −HPL)gr

m0−(

−τ)gr

m+1,0−(

S†

−τ†)gr

m−1,0=0

··· (B.2)

In the latter set of equations we can express all the blocks gr

m0, with modd, as the following:

m0= (

SPL −HPL)−1[(

−τ)gr

m+1,0+(

S†

−τ†)gr

m−1,0](B.3)

The new expression for every odd blocks of the Green’s function can be substituted in the

neighbor equations, getting the following

131

sgr

00 +τ1

1gr

20 =1

bgr

20 +τ1

1gr

40 +τ1

2gr

00 =0

bgr

60 +τ1

1gr

20 +τ1

2gr

20 =0

bgr

80 +τ1

1gr

20 +τ1

2gr

40 =0

···

sgr

m0+τ1

1gr

m+2,0+τ1

2gr

m−2,0=0

··· (B.4)

where

s= (

SPL −HPL)−(

−τ)(

SPL −HPL)−1(

S†

PL −τ†)

b= (

SPL −HPL)−(

SPL −τ)(

SPL −HPL)−1(

S†

−τ†)+

−(

S†

PL −τ†)(

SPL −HPL)−1(

−τ)

τ1

1=−(

−τ)(

SPL −HPL)−1(

−τ)

τ1

2=−(

S†

−τ†)(

SPL −HPL)−1(

S†

−τ†)(B.5)

We have actually recovered Eqn. (B.2) with new renormalized parameters. The first equa-

tion in Eqn. (B.4) connects now the renormalized 0th and 2nd layers and τ1and τ2are giving

the effective interaction between them, while Wsis the effective matrix of (

SPL −He f f )

for layer “0”, that means the layer in contact with the device.

We proceed in successive steps by eliminating the Green’s functions of the next alter-

nating layers, which will now be gr

2n,0for odd n. This procedure can be iterated, obtaining

at each step the corresponding renormalized Wb,Ws,τ1and τ2matrices relating to those

of the previous step in the following way:

s=Wi−1

s−τi−1(Wi−1

b)−1τi−1

s=Wi−1

b−τi−1(Wi−1

b)−1τi−1

2−τi−1

2(Wi−1

b)−1τi−1

τi

1=−τi−1

1(Wi−1

b)−1τi−1

τi

1=−τi−1

2(Wi−1

b)−1τi−1

2.(B.6)

After the ith step τi

1and τi

2give the effective interaction between layers “0” and 2ihaving

renormalized out those in between. We can thus expect both τ1and τ2to decrease for

any energy value with increasing number of steps, becoming quickly negligible as the

layers they effectively connect are just too far apart. The surface and bulk layers are then

practically decoupled, and we are left with the equations:

132

Appendix B. Surface Green’s Function: Decimation Technique

sgr

00 =1

bgr

2Nn,0=0(B.7)

where N is the number of steps necessary to achieve convergence in the iteration process.

Convergence is achieved when τN

1and τN

where

is a defined small threshold. The

Green’s function for the surface is therefore well described by

00 =¡WN

s¢−1(B.8)

In this example we have assumed that the chain of PLs grows in the right direction, which

is the case of the right contact. The same algorithm can be applied for the left contact,

where the PLs grows in the left direction, just exchanging τ1with τ2in the equations. In

order to run this algorithm we need only the Hamiltonian of two PLs and their coupling for

both the contacts.

Bibliography

[1] International Technology Roadmap for Semiconductors (Semiconductor Industry Associa-

tion), (2005).

[2] G. Moore. Cramming more components onto integrated circuits. Electronics, 38:114–117,

(1965).

[3] R. Chau, S. Datta, and A. Majumdar. Opportunities and challenges of iii-v nanoelectron-

ics for future high-speed, low-power logic applications. IEEE CSIC Symposium Technical

Digest, pages 17–20, (2005).

[4] G. D. Wilk, R. M. Wallace, and J. M. Anthony. High-kappa gate dielectrics: Current status

and materials properties considerations. J. Appl. Phys., 89:5243–5275, (2001).

[5] H. S. P. Wong, D. J. Frank, P. M. Solomon, C. H. J. Wann, and J. J. Welser. Nanoscale

CMOS. Proc. IEEE, 87:537–570, (1999).

[6] L. L. Chang, M. Ieong, and M. Yang. CMOS circuit performance enhancement by surface

orientation optimization. IEEE Trans. Electron. Dev., 51:1621–1627, (2004).

[7] D. A. Muller. A sound barrier for silicon? Nature Mater., 4:645–647, (2005).

[8] M. A. Reed. Molecular electronics: Back under control. Nature Mater., 3:286–287, (2004).

[9] C. A. Richter, D. R. Stewart, D. A. A. Ohlberg, and R. Stanley Williams. Electrical charac-

terization of Al/AlOx/molecule/Ti/Al devices. Appl. Phys. A, 80:1355–1362, (2005).

133

134

BIBLIOGRAPHY

[10] W. Y. Wang, T. H. Lee, and M. A. Reed. Electronic transport in molecular self-assembled

monolayer devices. Proc. IEEE, 93:1815–1824, (2005).

[11] J. Park, A. N. Pasupathy, J. I. Goldsmith, C. Chang, Y. Yaish, J. R. Petta, M. Rinkoski, J. P.

Sethne, H. D. Abruna, P. L. McEuen, and D. C. Ralph. Coulomb blockade and the Kondo

effect in single-atom transistors. Nature, 417:722, (2002).

[12] J. Chen and M. A. Reed. Electronic transport of molecular systems. Chem. Phys., 281:127,

(2002).

[13] M. A. Reed, C. Zhou, C. J. Muller, T. P. Burgin, and J. M. Tour. Conductance of a Molecular

Junction. Science, 278 (5336):252, (1997).

[14] C. Zhou, C. Muller, M. Deshpande, J. Sleight, and M. A. Reed. Microfabrication of a me-

chanically controllable break junction in silicon. App. Phys. Lett., 76:1160, (1995).

[15] J. van Ruitenbeek, A. Alvarez, I. Pineyro, C. Grahmann, P. Joyez, M. Devoret, D. Esteve,

and C. Urbina. Adjustable nanofabricated atomic size contacts. Rev. of Scien. Inst., 67:108,

(1996).

[16] C. Muller, J. van Ruitenbeek, and L. de Jongh. Experimental observation of the transition

from weak link to tunnel junction. Physica C, 191:485, (1992).

[17] E. Sheer, P. Joyez, D. Esteve, C. Urbina, and M. Devoret. Conduction channels transmissions

of atomic size aluminium contacts. Phys. Rev. Lett., 78:3535, (1997).

[18] J. Reichert, R. Ochs, D. Beckmann, H. Weber, M. Mayor, and H. von L ¨ohneysen. Driving

current through single organic molecules. Phys. Rev. Lett., 88:176804, (2002).

[19] H. Weber, J. Reichert, F. Weigend, R. Ochs, D. Beckmann, M. Mayor, R. Ahlrichs, and

H. von L ¨ohneysen. Electronic transport through single conjugated molecules. Chem. Phys.,

281:113–125, (2002).

[20] M. Mayor, C. von H ¨anisch, H. Weber, J. Reichert, and D. Beckmann. A trans-Platinum (II)

complex as a single-molecule insulator. Angw. Chem. Int. Ed., 41 (7):1183, (2002).

BIBLIOGRAPHY

135

[21] C. Kergueris, J. P. Bourgoin, S. Palacin, D. Esteve, C. Urbina, M. Magoga, and C. Joachim.

Electron transport through a metal-molecule-metal junction. Phys. Rev. B, 59:12505, (1999).

[22] C. Lau, D. Stewart, R. Williams, and M. Bockrath. Direct observation of nanoscale switching

centers in metal/molecule/metal structures. Nano Lett., 4:569, (2004).

[23] X. Cui, A. Primak, X. Zarate, J. Tomfohr, O. Sankey, A. Moore, T. Moore, D. Gust, G. Haris,

and S. Lindsay. Reproducible measurements of single-molecule conductivity. Science,

294:571, (2001).

[24] L. Patrone, S. Palacin, J. Bougoin, J. Lagoute, T. Zambelli, and S. Gauthier. Direct com-

parison of the electronic coupling efficiency of sulfur and selenium anchoring groups for

molecules adsorbed onto gold electrodes. Chem. Phys., 281:325–332, (2002).

[25] C. Collier, G. Mattersteig, E. Wong, Y. Luo, K. Beverly, J. Samapaio, F. Raymo, J. Stoddart,

and J. Health. A [2]Catenane-based solid state electronically reconfigurable switch. Science,

289:1172, (2000).

[26] Z. Donhauser, B. Mantooth, K. Kelly, L. Bumm, J. Monnell, J. Stapleton, D. Price Jr.,

A. Rawlett, D. Allara, J. Tour, and P. Weiss. Conductance switching in single molecules

through conformational changes. Science, 292:2303, (2001).

[27] C. Zhou, M. R. Deshpande, M. A. Reed, L. Jones II, and J. M. Tour. Nanoscale metal/self-

assembled monolayer/metal heterostructures. Appl. Phys. Lett., 71:611, (1997).

[28] Y. Kervennic, D. Vanmaekelbergh, L. Kouwenhoven, and H. Van der Zant. Planar nanocon-

tacts with atomically controlled separation. App. Phys. Lett., 83:3782, (2003).

[29] H. Park, A. Lim, A. Alivisatos, J. Park, and P. McEuen. Fabrication of metallic electrodes

with nanometer separation by electromigration. App. Phys. Lett., 75:301, (1999).

[30] M. Lambert, M. Goffman, J. Bourgoin, and P. Hesto. Fabrication and characterisation of

sub-3 nm gaps for single-cluster and single-molecule experiments. Nanotechnology, 14:772,

(2003).

136

BIBLIOGRAPHY

[31] S. Khondaker and Z. Yao. Fabrication of nanometer-spaced electrodes using gold nanoparti-

cles. App. Phys. Lett., 81:4613, (2002).

[32] B. Xu and N. J. Tao. Measurement of single-molecule resistance by repeated formation of

molecular junctions. Science, 301:1221, (2003).

[33] L. Bumm, J. Arnold, M. Cygan, T. Dunbar, T. Burgin, L. Jones II, D. Allara, J. Tour, and

P. Weiss. Are single molecular wires conducting? Science, 271:1705, (1996).

[34] X. Xiao, B. Xu, and N. Tao. Measurement of single molecule conductance: benzendithiol

and benzendimethanethiol. Nano Lett., 4:267, (2004).

[35] P. K. Hansma. Tunneling Spectroscopy. Plenum Press, New York, 1982.

[36] J. T. Yates and T. E. Madey. Vibrational Spectroscopy of Molecules on Surfaces. Plenum

Press, New York, 1987.

[37] K.W. Hipps and U. Mazur. Vibrational and Low-Lying Electronic Transitions in Tetraalky-

lammonium Salts of CoBr2−

4, CoCl2−

4, Co(CNS)2−

4As Observed by Raman, Infrared, and

Tunneling Spectroscopies. J. Phys. Chem., 91:5218, (1987).

[38] A. Troisi and M. A. Ratner. Molecular transport junctions: Propensity rules for inelastic

electron tunneling spectra. Nano Lett., 6:1784, (2006).

[39] N. Lorente, M. Persson, L. J. Lauhon, and W. Ho. Symmetry Selection Rules for Vibra-

tionally Inelastic Tunneling. Phys. Rev. Lett., 86:2593, (2001).

[40] N. Lorente and M. Persson. Theory of Single Molecule Vibrational Spectroscopy and Mi-

croscopy. Phys. Rev. Lett., 85:2997, (2000).

[41] H. Kato and H. Shinohara. Spatially dependent inelastic tunneling in a single metallo-

fullerene. Phys. Rev. Lett., 94:136802, (2005).

[42] Y.-C. Chen, M. Zwolak, and M. Di Ventra. Inelastic effects on the transport properties of

alkanethiols. Nano Lett., 5:621, (2005).

BIBLIOGRAPHY

137

[43] M. Galperin, M. A. Ratner, and A. Nitzan. On the line widths of vibrational features in

inelastic electron tunneling spectroscopy. Nano Lett., 4:1605, (2004).

[44] Tikhodeev and H. Ueba. Relation between inelastic electron tunneling and vibrational exci-

tation of single adsorbates on metal surfaces. Phys. Rev. B, 70:125414, (2004).

[45] A. Troisi, M. A. Ratner, and A. Nitzan. Vibronic effects in off-resonant molecular wire

conduction. J. Chem. Phys., 118:6072, (2003).

[46] J. Jiang, M. Kula, W. Lu, and Y. Luo. First-principles simulations of inelastic electron tun-

neling spectroscopy of molecular electronic devices. Nano Lett., 5:1551, (2005).

[47] Y. Meir and N. Wingreen. Landauer formula for the current through an interacting electron

region. Phys. Rev. Lett., 68:2512, (1992).

[48] M. B ¨uttiker. Four terminal phase-coherent conductance. Phys. Rev. Lett., 57:1761, (1986).

[49] M. B ¨uttiker, Y. Imry, R. Landauer, and S. Pinhas. Generalized many-channel conductance

formula with application to small rings. Phys. Rev. B, 31:6207, (1985).

[50] G. D. Mahan. Many-Particle Physics. Plenum, second edition, 1993.

[51] R. D. Mattuck. A Guide to Feynman Diagrams in the Many-Body Problem. Dover, second

edition, 1992.

[52] A. L. Fetter and J. D. Waleka. Quantum Theory of Many-Particle System. Dover, first edition,

2003.

[53] S. Datta. Electronic Transports in Mesoscopic Systems. Cambridge, 1995.

[54] T. N. Todorov. Tight-binding simulation of current-carrying nanostructures. J. Phys. Cond.

Matt., 14:3049, (2002).

[55] L. V. Keldysh. Diagram technique for non-equilibrium processes. Soviet Physics JETP,

20:1018, (1964).

138

BIBLIOGRAPHY

[56] A. Pecchia, G. Romano, and A. Di Carlo. Theory of heat dissipation in molecular electronics.

Phys. Rev. B, 75:035401, (2007).

[57] C. Caroli, R. Combescot, P. Nozieres, and D. Saint-James. Direct calculation of the tunneling

current. J. Phys. C: Solid St. Phys., 4:916, (1971).

[58] Leo P. Kadanoff and Gordon Baym. Quantum Statistical Mechanics. W. A. Benjamin, 1962.

[59] M. Gell-Mann and F. E. Low. Quantum Electrodynamics at small distances. Phys. Rev,

95:1300, (1954).

[60] T. Frederiksen, M. Paulsson, M. Brandbyge, and J. Antti-Pekka. Inelastic transport the-

ory from first-principles: methodology and applications for nanoscale devices. cond-matt,

(2006).

[61] Hartmut Haug and Anti-Pekka Jahuo. Quantum Kinetics in Transport and Optics of Semi-

conductors. Springer, 1996.

[62] T. Frauenheim, G. Seifert, M. Elstner, T. Niehaus, C. Koeler, M. Amkreutz, M. Sternberg,

Z. Hajnal, A. Di Carlo, and S. Suhai. Atomistic simulations of complex materials: ground-

state and excited-state properties. J. Phys.: Condensed Matter, 14:3015, (2002).

[63] T. Frauenheim, G. Seifert, M. Elstner, Z. Hajnal, G. Jungnickel, D. Porezag, S. Suhai, and

R. Scholz. A Self-Consistent Charge Density-Functional Based Tight-Binding Method for

Predictive Materials Simulations in Physics, Chemistry and Biology. Phys. Stat. Sol. (b), 217

(1):41, (2000).

[64] M. Elstner, D. Porezag, G. Jugnickel, J. Elsner, M. Haugk, T. Frauenheim, S. Suhai, and

G. Seifert. Self-consistent-charge density-functional tight-binding method for simulations of

complex materials properties. Phys. Rev. B, 58:7260, (1998).

[65] R. M. Dreizler and E. K. U. Gross. Density Functional Theory. Springer, Berlin and Heidel-

berg, 1990.

[66] R. O. Jones and O. Gunnarsson. The density functional formalism, its application and

prospects. Rev. Mod. Phys., 61:689, (1989).

BIBLIOGRAPHY

139

[67] L. H. Thomas. The calculation of atomic fields. Proc. Camb. Phil. Soc., 23:542, (1926).

[68] E. Fermi. Eine statistische Methode zur Bestimmung einiger Eigenschaften des Atoms und

ihre Anwendung auf die Theorie des periodischen Systems der Elemente. Z. Phys., 48:73,

(1928).

[69] W. Kohn. Nobel lecture: Electronic structure of matter-wave functions and density function-

als. Rev. Mod. Phys., 71:1253, (1999).

[70] J. F. Janak. Proof that

∂

iin density-functional theory. Phys. Rev. B, 18:7165,

(1978).

[71] J. C. Slater and G. F. Koster. Simplified LCAO method for the periodic potential problem.

Phys. Rev., 94:1498, (1954).

[72] G. Seifert, H. Eschrig, and W. Bieger. Eine approximative Variante des LCAO-X

Ver-

fahrens. Z. Phys. Chem., 267:529, (1986).

[73] J. A. Pople, D. P. Santry, and G. A. Segal. Approximate self-consistent molecular orbital

theory. I. Invariant procedures. J. Chem. Phys., 43:S129, (1965).

[74] M. Scholz and H.-J. Koehler. Quantenchemische Naeherungsverfahren und ihre Anwendung

in der organischen Chemie, vol. 3 of Quantenchemie-Ein Lehrgang. Deutscher Verlag der

Wissenschaften, 1981.

[75] J. A. Pople and D. L. Beveridge. Approximate Molecular Orbital Theory. McGraw Hill,

New York, 1970.

[76] J. R. Reimers, Z. L. Cai, A. Bili´c, N. S. Hush, and N. Y. Ann. Acad. Sci., 235:1006, (2003).

[77] R. Neumann, R. H. Nobes, and N. C. Handy. Molec. Phys., 87:1, (1996).

[78] Z.-L. Cai and J. R. Reimers. Application of time-dependent density-functional theory to the

3Σ−

ufirst excited state of H2.J. Chem. Phys., 112:527, (2000).

140

BIBLIOGRAPHY

[79] G. C. Solomon, J. R. Reimers, and N. S. Hush. Single molecule conductivity: The role of

junction-orbital degeneracy in the artificially high currents predicted by ab initio approaches.

J. Chem. Phys., 121:6615, (2004).

[80] J. Jortner, A. Nitzan, and M. A. Ratner. Lecture Notes in Physics: Introducing Molecular

Electronics.

[81] O. Gunnarsson and B. I. Lundqvist. Exchange and correlation in atoms, molecules, and

solids by the spin-density-functional formalism. Phys. Rev. B, 13:4274, (1976).

[82] A. Bili´c, J. R. Reimers, and N. S. Hush. The structure, energetics, and nature of the chemical

bonding of phenylthiol adsorbed on the Au(111) surface: Implications for density-functional

calculations of molecular-electronic conduction. J. Chem. Phys., 122:094708, (2005).

[83] T. A. Niehaus, S. Suhai, F. Della Sala, P. Lugli, M. Elstner, G. Seifert, and T. Frauenheim.

Tight-binding approach to time-dependent density-functional response theory. Phys. Rev. B,

63:085108/1, (2001).

[84] Z.-L. Cai, K. Sendt, and J. R. Reimers. Failure of density-functional theory and time-

dependent density-functional theory for large extended

systems. J. Chem. Phys., 117:5543,

(2002).

[85] S. J. A. van Gisbergen, P. R. T. Schipper, O. V. Gritensko, E. J. Baerends, J. D. Snijders,

B. Champagne, and B. Kirtman. Electric field dependence of the exchange-correlation po-

tential in molecular chains. Phys. Rev. Lett., 83:694, (1999).

[86] A. Pecchia and A. Di Carlo. Atomistic theory of transport in organic and inorganic nanos-

tructures. Rep. Prog. Phys., 67:1497, (2004).

[87] B. A. Mantooth, Z. J. Donauser, K. F. Kelly, and P. S. Weiss. Cross-correlation image tracking

for drift correction and adsorbate analysis. Rev. Scien. Instr., 73:313, (2002).

[88] M. Galperin, M. A. Ratner, and A. Nitzan. Molecular transport junctions: Vibrational effects.

J. Phys. Cond. Matt., 19:103201, (2007).

BIBLIOGRAPHY

141

[89] J. G. Kushmerick, J. Lazorcik, C. H. Patterson, R. Shashidhar, D. S. Seferos, and G. C.

Bazan. Vibronic contributions to charge transport across molecular junctions. Nano Lett.,

4:639, (2004).

[90] W. L. Wang, I. Kretzschmar, and M. Reed. Inelastic electron tunneling spectroscopy of an

alkanedithiol self-assembled monolayer. Nano Lett., 4:643, (2004).

[91] Y. Yourdshahyan, H. K. Zhang, and A. M. Rappe. n-alkyl thiol head-group interactions with

the Au(111) surface. Phys. Rev. B, 63:081405, (2001).

[92] J. Gottschalck and B. J. Hammer. A density functional theory study of the adsorption of

sulfur, mercapto, and methylthiolate on Au(111). J. Chem. Phys., 116:784, (2002).

[93] H. S. Kato, J. Noh, M. Hara, and M. Kawai. An HREELS Study of Alkanethiol Self-

Assembled Monolayers on Au(111). J. Phys. Chem. B, 106:9655, (2002).

[94] M. A. Bryant and J. E. Pemberton. Surface Raman Scattering of Self-Assembled Monolayers

Formed from 1-Alkanethiols: Behavior of Films at Au and Comparison to Films at Ag. J.

Am. Chem. Soc., 113:8284, (1991).

[95] C. Castiglioni, M. Gussoni, and G. Zerbi. Charge mobility in

-bonded molecules: The

infrared spectrum of polymethylene chains in the solid and liquid phases. J. Chem. Phys.,

95:7144, (1991).

[96] Y.-C. Chen, M. Zwolak, and M. Di Ventra. Inelastic current-voltage characteristics of atomic

and molecular junctions. Nano Lett., 4:1709, (2004).

[97] A. Troisi and M. A. Ratner. Modeling the inelastic electron tunneling spectra of molecular

wire junctions. Phys. Rev. B, 72:033408, (2005).

[98] W. Andreoni, A. Curioni, and H. Gronbeck. Int. J. Quantum Chem., 80:598, (2000).

[99] P. Fenter, F. Schreiber, L. Berman, G. Scoles, P. Eisenberger, and M. J. Bedzyk. On the

structure and evolution of the buried S/Au interface in self-assembled monolayers: X-ray

standing wave results. Surf. Sci., 412/413:213, (1998).

142

BIBLIOGRAPHY

[100] G. J. Kluth, C. Carraro, and R. Marboudian. Direct observation of sulfur dimers in alka-

nethiol self-assembled monolayers on Au(111). Phys. Rev. B, 59:449, (1999).

[101] M. S. Yeganeh, S. M. Dougal, R. S. Polizzotti, and P. Rabinowitz. Interfacial atomic structure

of a self-assembled alkyl thiol monolayer/Au(111): A sum frequency generation-study. Phys.

Rev. Lett., 74:1811, (1995).

[102] Y. Akinaga, T. Nakajima, and K. Hirao. A density functional study on the adsorption of

methanethiolate on the (111) surfaces of noble metals. J. Chem. Phys., 114:8555, (2001).

[103] Y. Yourdshahyan and A. M. Rappe. Structure and energetics of alkanethiol adsorption on the

Au(111) surface. J. Chem. Phys., 117:825, (2002).

[104] L. M. Molina and B. Hammer. Theoretical study of thiol-induced reconstructions on the

Au(111) surface. Chem. Phys. Lett., 360:264, (2002).

[105] C. D. Bain, E. B. Troughton, Y.-T. Tao, J. Evall, M. Whitesides, George, and R. G. Nuzzo.

Formation of monolayer films by the spontaneous assembly of organic thiols from solution

onto gold. J. Am. Chem. Soc., 111 (1):321, (1989).

[106] L. H. Dubois, B. R. Zegarski, and R. G. Nuzzo. Molecular ordering of organosulfur com-

pounds on Au(l11) and Au(l00): Adsorption from solution and in ultrahigh vacuum. J. Chem.

Phys., 98 (1):678, (1993).

[107] M. D. Porter, T. B. Bright, D. L. Allara, and C. E. D. Chidsey. Spontaneously organized

molecular assemblies. 4. Structural characterization of n-alkyl thiol monolayers on gold by

optical ellipsometry, infrared spectroscopy, and electrochemistry. J. Am. Chem. Soc., 109

(12):3559, (1987).

[108] A. P. Scott and L. Radom. Harmonic vibrational frequencies: An evaluation of Hartree-

Fock, Møller-Plesset, Quadratic Configuration Interaction, Density Functional Theory, and

Semiempirical Scale Factors. J. Phys. Chem., 100:16502, (1996).

[109] K. W. Hipps and U. Mazur. Handbook of Vibrational Spectrum. Wiley, New York, 2002.

BIBLIOGRAPHY

143

[110] H. Basch, R. Cohen, and M. A. Ratner. Interface geometry and molecular junction conduc-

tance: Geometric fluctuation and stochastic switching. Nano Lett., 5:1688, (2005).

[111] H. Oevering, J. W. Verhoeven, M. N. Paddon-Row, E. Cotsaris, and N. S. Hush. Long-range

exchange contribution to singlet-singlet energy transfer in a series of rigid bichromophoric

molecules. Chem. Phys. Lett., 8:117, (1985).

[112] M. Brandbyge, M. R. Sorensen, and K. W. Jacobsen. Conductance eigenchannels in

nanocontacts. Phys. Rev. B, 56:14956, (1997).

[113] D. Djukic and J. M. van Ruitenbeek. Shot noise measurements on a single molecule. Nano

Lett., 6:789, (2006).

[114] M. Kemp, A. Roitberg, V. Mujica, T. Wanta, and M. A. Ratner. Molecular wires: Extended

coupling and disorder effects. J. Phys. Chem., 100:8349, (1996).

[115] M. B ¨uttiker. Spatial variation of currents and fields due to localized scatterers in metallic

conduction. IBM J. Res. Dev., 32:63, (1988).

[116] J. Heurich, J. C. Cuevas, W. Wenzel, and G. Schon. Electrical transport through single-

molecule junctions: From molecular orbitals to conduction channels. Phys. Rev. Lett.,

88:256803, (2002).

[117] J. C. Cuevas, A. L. Yeyati, and A. Martin-Rodero. Microscopic origin of conducting channels

in metallic atomic-size contacts. Phys. Rev. Lett., 80:1066, (1998).

[118] D. Jacob and J. J. Palacios. Orbital eigenchannel analysis for ab initio quantum transport

calculations. Phys. Rev. B, 73:075429, (2006).

[119] A. Bagrets, N. Papanikolaou, and I. Mertig. Ab initio approach to the ballistic transport

through single atoms. Phys. Rev. B, 73:045428, (2006).

[120] A. Troisi and M. A. Ratner. Propensity rules for inelastic electron tunneling spectroscopy of

single-molecule transport junctions. J. Chem. Phys., 125, (2006).

144

BIBLIOGRAPHY

[121] M. Di Ventra, S. T. Pantelides, and N. D. Lang. First-principles calculation of transport

properties of a molecular device. Phys. Rev. Lett., 84:979, (2000).

[122] J. B. Neaton, M. S. Hybertsen, and S. G. Louie. Renormalization of molecular electronic

levels at metal-molecule interfaces. Phys. Rev. Lett., 97:216405, (2006).

[123] T. Niehaus, M. Rohlfing, F. Della Sala, A. Di Carlo, and T. Frauenheim. Quasiparticle

energies for large molecules: A tight-binding-based Green0s-function approach. Phys. Rev.

A, 71:022508, (2005).

[124] A. Gagliardi, T. A. Niehaus, T. Frauenheim, A. Pecchia, and A. Di Carlo. Quasiparticle

correction for electronic transport in molecular wires. J. Comp. Elect., 6:345, (2007).

[125] A. Aviram and M. A. Ratner. Molecular rectifier. Chem. Phys. Lett., 29:277, (1974).

[126] F. Guinea, C. Tejedor, F. Flores, and E. Louis. Effective two-dimensional Hamiltonian at

surfaces. Phys. Rev. B, 28:4397, (1983).

BIBLIOGRAPHY

145

Own Publications

1. The gDFTB Tool for Molecular Electronics,

A. Pecchia, L. Latessa, A. Gagliardi, Th. Frauenheim, A. Di Carlo,

Molecular and Nano Electronics: Analysis, Design and Simulation, Vol. 17, Editor J. Semi-

nario, Publisher Elsevier, 2006

2. A Priori Method for Propensity Rules for Inelastic Electron Tunneling Spectroscopy of Single-

Molecule Conduction,

A. Gagliardi, G. C. Solomon, A. Pecchia, Th. Frauenheim, A. Di Carlo, N. S. Hush and J. R.

Reimers,

Physical Review B, 2007, 75, 174306.

3. Understanding the Inelastic Electron-Tunneling Spectra of Alkanedithiols on Gold,

G. C. Solomon∗, A. Gagliardi∗, A. Pecchia, Th. Frauenheim, A. Di Carlo, J. R. Reimers and

N. S. Hush,

Journal of Chemical Physics, 2006, 124, 094704.

4. Molecular Origins of Conduction Channels Observed in Shot-Noise Measurements,

G. C. Solomon, A. Gagliardi, A. Pecchia, Th. Frauenheim, A. Di Carlo, J. R. Reimers and

N. S. Hush,

Nano Letters, 2006, 6, 2431.

5. The Symmetry of Single-Molecule Conduction,

G. C. Solomon, A. Gagliardi, A. Pecchia, Th. Frauenheim, A. Di Carlo, J. R. Reimers and

N. S. Hush,

Journal of Chemical Physics, 2006, 125, 184702.

6. Incoherent Electron-Phonon Scattering in Octanethiols,

A. Pecchia, A. Di Carlo, A. Gagliardi, S. Sanna, Th. Frauenheim and R. Gutierrez,

Nano Letters, 2004, 4, 2109.

146

BIBLIOGRAPHY

7. The Green’s Function Density-Functional Tight-Binding (gDFTB) Method for Molecular-

Electronic Conduction,

J. R. Reimers, G. C. Solomon, A. Gagliardi, A. Bili´c, N. S. Hush, Th. Frauenheim, A. Di

Carlo and A. Pecchia,

Journal of Physical Chemistry A, 2007, 111, 5692.

8. Simulations of Inelastic Tunnelling in Molecular Bridges,

A. Gagliardi, G. C. Solomon, A. Pecchia, A. Di Carlo, T. Frauenheim, J. R. Reimers, N. S.

Hush,

Non-Equilibrium Carrirer Dynamics in Semiconductor, Vol. 110, Editor M. Saraniti and U.

Ravaioli, Publisher Springer, 2006.

9. Atomistic Simulation of the Electronic Transport in Organic Nanostructures: Electron-Phonon

and Electron-Electron Interactions,

A. Pecchia, A. Di Carlo, A. Gagliardi, Th. Niehaus and Th. Frauenheim,

Journal of Computational Electronics, 4, 79, 2005.

10. Quasiparticle Correction for Electronic Transport in Molecular Wires,

A. Gagliardi, Th. A. Niehaus, Th. Frauenheim, A. Pecchia and A. Di Carlo,

Journal of Computational Electronics, 2007, 6, 345.

∗These authors contributed in a similar matter to publication.

Pieces of this thesis were published in (2), (3), (4), (5), (6) and (7).

Personal contributions

In the present work many topics were discussed. They were the results of collaborations with many

people and groups. In this last paragraph the main personal contributions of the author are enlisted.

The parts of the thesis which include new results are from Chapters 5to Chapter 9. The principal

contributions of the author are:

•Chapter 5: implementation of the Meir-Wingreen equation and the electron-phonon self-

energies in the Landauer-B¨

uttiker gDFTB code, Born approximation and part of the self

consistent Born approximation. Implementation of the calculation for the electron-phonon

couplings and the calculation of the power dissipated by virtual contacts.

•Chapter 6: Supervision for the code, analysis of the results with Dr. Alessandro Pecchia.

•Chapter 7: derivation of the simplified equation for the IETS. The contribution for all the

calculations, analysis and conclusions are equally shared with Dr. Gemma C. Solomon.

•Chapter 8: Supervision for the code, analysis and discussion of the results with Dr. Gemma

C. Solomon.

•Chapter 9: All the theoretical derivation, calculations, analysis and discussions performed by

the author.

148

Acknowledgements

I am most grateful to Prof. Dr. Frauenheim for his friendship, help and for all the freedom he gave

me in all these years. It is only thank to him if I was able to do so many collaboration and learning

so much.

I want to acknowledge also Prof. Jeffrey R. Reimers, of the University of Sydney, for his

invaluable help. I thank him also for his willingness and kindness. I am also indebted deeply to

Prof. Aldo Di Carlo of the University “Tor Vergata” in Rome. It is only thank to him if I had the

opportunity to start this PhD.

I want to express also my gratitude to some other people, first Dr. Alessandro Pecchia, for all

the interesting discussions about science and life in general. I am most grateful also to Dr. Gemma

C. Solomon for all the work done together and all the stimulating discussions we have in the last

three years. A nice friend in the right side of the world. I deep in debt with Marius Wanko and

Dr. Balint Aradi, especially for their revision of the manuscript of the thesis. I want to thank also

Simone Lange and Simone Sanna for their friendship and assistance.

I am also very grateful to Prof. Noel S. Hush of the University of Sydney for the discussion

concerning selection rules in inelastic electron tunneling spectroscopy. I would also express my

gratitude to the entire group of Bremen Material Computer Science and the component of the previ-

ous group in Paderborn, in particular Pia T¨

oelle, Dr. Thomas Niehaus and David Heringer. Thanks

to all of them for all the help provided in these three years. Finally I want to acknowledge the group

of Paderborn Institute for Scientific Computation (PaSco) for the support, the assistance and the

funding provided during my PhD. My thanks to all those whom, by mistake, I forgot to mention.

The last part of the acknowledgment is in Italian, for some very special persons. Innanzitutto

voglio ringraziare la mia famiglia, Enrica, mia madre, mio fratello, Silvia e Luciano. Persone

senza le quali questo lavoro non sarebbe sicuramente mai stato concluso. A queste persone devo

aggiungere un cane (Tango) e un gatto, anzi una Gatta (lei sa). Voglio inoltre ringraziare i miei

amici pi´

u stretti, Alessio, Andrea, Michele, Gabriele, Paolo, Mauro, Alessandro (Dottor) e Alberto.

Molta strada si ´

e fatta insieme e molta se ne far´

a ancora.

Per ultimi, ma non ultimi, due ringraziamneti speciali per due persone: una, purtroppo, partita e

149

l’altra appena arrivata. Mio padre e la piccola Giorgia. ´

E stato infatti grazie a mio padre se ho com-

inciato ad apprezzare le scienze e quindi una buona parte del merito (o della colpa) se io sono qui

e sua. Ricordo quando mi narr´

o il primo racconto di fantascienza, “Incontro con Rama” di Arthur

Clarke. Il suo riassunto pu´

o essere integralmente riportato su queste pagine: “Una astronave aliena

arriva nel nostro sistema solare e se ne va”. Mio padre era cos´

ı, ingegnere fino al midollo. Quando

chiesi maggiori dettagli mi regal´

o la sua edizione del romanzo. La seconda persona ´

e appena nata.

Questa Tesi ´

e anche dedicata a lei e spero che le porti fortuna.