0lopscience
¡opscience.iop.org
Home Search Collections Journals About Contact us My IOPscience
The treatment of electronic excitations in atomistic models of radiation damage in metals
This article has been downloaded from IOPscience. Please scroll down to see the full text article.
2010 Rep. Prog. Phys. 73 116501
(http://iopscience.iop.org/0034-4885/73/11/116501)
View the table of contents for this issue, or go to the journal homepage for more
Download details: IP Address: 129.174.21.5
The article was downloaded on 24/01/2013 at 01:33
Please note that terms and conditions apply.
Rep. Prog. Phys. 73 (2010) 116501 (40pp) doi:10.1088/0034-4885/73/11/116501
The treatment of electronic excitations in atomistic models of radiation damage in metals
C P Race1, D R Mason1, M W Finnis12, W M C Foulkes1, A P Horsfield2 and A P Sutton1
1 Department of Physics, Imperial College, London SW7 2AZ, UK
2 Department of Materials, Imperial College, London SW7 2AZ, UK
Received 31 December 2009
Published 26 October 2010
Online at stacks.iop.org/RoPP/73/116501
Abstract
Atomistic simulations are a primary means of understanding the damage done to metallic materials by high energy particulate radiation. In many situations the electrons in a target material are known to exert a strong influence on the rate and type of damage. The dynamic exchange of energy between electrons and ions can act to damp the ionic motion, to inhibit the production of defects or to quench in damage, depending on the situation. Finding ways to incorporate these electronic effects into atomistic simulations of radiation damage is a topic of current major interest, driven by materials science challenges in diverse areas such as energy production and device manufacture.
In this review, we discuss the range of approaches that have been used to tackle these challenges. We compare augmented classical models of various kinds and consider recent work applying semi-classical techniques to allow the explicit incorporation of quantum mechanical electrons within atomistic simulations of radiation damage. We also outline the body of theoretical work on stopping power and electron-phonon coupling used to inform efforts to incorporate electronic effects in atomistic simulations and to evaluate their performance.
(Some figures in this article are in colour only in the electronic version)
This article was invited by M-Y Chou.
Contents
1. Introduction
1.1. The radiation damage cascade
1.2. Transfer of energy between ions and electrons in radiation damage
2. Models of electronic stopping power
2.1. Models of the stopping of fast, light particles
2.2. Models of the stopping of fast, heavy particles
2.3. Models of the stopping of slow, heavy particles
2.4. Empirical models of electronic stopping power
3. Models of electron-phonon coupling
3.1. The importance of electron-phonon coupling
3.2. Two-temperature models
3.3. The electron-phonon drag
3.4. Models of electron-phonon coupling
3.5. Estimates of electron-phonon coupling
4. Implicit incorporation of electronic effects in atomistic simulations 20
4.1. Experimental benchmarks 21
4.2. Damage functions 21
4.3. Binary collision models 22
4.4. Molecular dynamics with drag 23
4.5. Electrons as a heat bath 24
5. Explicit treatment of electrons 25
5.1. Electrons as an inhomogeneous heat bath 25
5.2. Adding quantum mechanical electrons 28
5.3. Time-dependent density functional theory 29
5.4. Time-dependent tight-binding 31
5.5. Correlated electron ion dynamics (CEID) 34
6. Concluding remarks 35 Acknowledgment 36 References 36
0034-4885/10/116501+40$90.00
© 2010 IOP Publishing Ltd Printed in the UK & the USA
1. Introduction
The study of the damage caused to materials by high energy particles is of huge practical importance. Beams of energetic particles are used as experimental probes, as therapeutic tools and as a means of materials modification in manufacturing processes. Unintended radiation damage to materials in many applications demands a detailed understanding of the physical processes at work.
Structural materials in fission and fusion reactor environments are subject to intense bombardment by high energy neutrons. Given that commercial viability demands long plant lifetimes, it is likely that future reactor materials will have to cope with lifetime dose rates of over 100 displacements per atom (dpa) [1]. Such high levels of atomic disturbance within an often carefully designed and fine-scaled microstructure can significantly influence the choice of material. The damage caused by irradiation can promote large long-term changes in bulk material properties, some of which, such as embrittlement, can threaten dangerous failure.
Characterizing this damage and understanding the mechanisms by which it is caused is central to efforts to improve the lifetime of materials in hostile environments. Experiment can only take us so far. Typical damage distributions have length scales that are accessible to the highest resolution experimental probes, but they are formed in times of no more than a few picoseconds; too fast to measure directly.
Filling this experimental gap, a century of theoretical work on the subject of radiation damage has built an impressive array of models, treating the interactions of one type of ion with another and treating in turn the interactions of those ions with the system of electrons in which they sit. In the last half century, these models have been incorporated into computer simulations of radiation damage events, which, as computing power has increased, have become more sophisticated, attacking ever higher energies and larger length and time scales. Today, computing resources allow simulations of radiation damage to advance to the next stage. They will move beyond the realm of classical molecular dynamics and towards a full quantum mechanical treatment of the metallic system.
A current challenge, acknowledged by recent efforts in the field, is to provide a better understanding of the effect of the interaction between ions and electrons on the outcome of radiation damage events [2] (see [3] for a general review of such effects in radiation damage). The electrons in a metal contribute to the potential in which the ions move and their state of excitation will affect that potential. Classical simulations assume the validity of the Born-Oppenheimer approximation: that the electrons will respond instantaneously to the motion of the ions. In fact their response time is finite and this will manifest itself increasingly strongly at higher ionic velocities. The electrons will also behave as a heat bath, in thermal contact with the system of ions. This heat bath will influence the production and healing of defects in the ionic system and the comparatively high thermal conductivity of the electrons will provide a means of enhanced energy transport away from areas of ionic disturbance.
In this paper we will consider the variety of ways in which electronic effects can be incorporated into atomistic models of radiation damage in metals, by which we mean simulation methods that follow the evolution of an explicitly represented set of ions. We will review the well established and much used analytical models and their application in dynamic simulations to yield an implicit treatment of energy loss to electrons. But we will also cover more recent simulation work, in which electrons are treated explicitly, and examine to what extent the results validate previous conclusions. And we will look to the future, identifying some techniques that improvements in computational power may soon render applicable to the simulation of radiation damage events.
We do not aim to present a comprehensive review of collision and stopping theory. Detailed discussions of these subjects are available elsewhere (see, for example, Bohr's review [4], the discussion of the work of Lindhard and co-workers by Ziegler et al [5] and the more recent book by Sigmund [6]). Instead we will consider only those aspects of stopping theory relevant to a discussion of atomistic simulation in metals, highlighting the essential aspects of the underlying physics that are included within or excluded from the various approximations made.
In the rest of this introduction we will discuss the evolution of a typical radiation damage event, considering where electronic effects might make themselves felt, and then introduce the basic theoretical framework on which much of the discussion in the literature of energy exchange between electrons and ions is based. In sections 2 and 3 we will describe the various analytical approaches that have been applied to energy transfer in different regimes of ionic motion. In sections 4 and 5 we will discuss the history, the state of the art and the future of atomistic simulations of radiation damage in metals.
1.1. The radiation damage cascade
Before delving into the nature of the interactions between the various particles involved in radiation damage, we will describe a typical event, introducing some of the terminology prevalent throughout the literature along the way. Much of our discussion will focus on the notion of a radiation damage collision cascade developing within a solid block of some material. The various ultimate causes of damage are brought together in a single phenomenological framework by considering the cascade to be initiated when an ion within the material, known as the primary knock-on atom (PKA), is set in motion, typically with a very high energy. This PKA is often an ion of the target material itself, accelerated by a collision with an intruding particle such as a neutron or a-particle, but it might also be the recoiling product of a radioactive decay process taking place, for example, within the storage medium for radioactive waste.
The distribution of kinetic energy of the PKA (its spectrum) and the subsequent pattern of development of the collision cascade depend on the initiating event [1]. Collisions with the 14 MeV neutrons emerging from a deuterium-tritium fusion reaction will produce recoil energies in iron of up to
1MeV with half of all recoils being above 10keV. The somewhat slower neutrons emerging from nuclear fission reactors produce recoil energies of up to several hundred kiloelectronvolts. The recoiling 234U nucleus from the decay of 238Pu will have an energy of around 100 keV.
A PKA with energy ~ 100 keV and above will have a very small cross-section for interaction with the nuclei of the solid and can travel large distances without undergoing a significant collision in a process known as channelling. Experiments in which 40 keV radioactive 125Xe ions are implanted into crystalline tungsten along channelling axes [7] show significant penetration at depths of up to 10-6 m. Comparison with results for amorphous tungsten, in which no penetration exceeds 0.1 x 10-6m, suggests that channelling plays a key role in determining the spatial distribution of damage. A channelling ion will progressively lose energy to the electrons of the solid and ultimately the channelling behaviour will end with a collision between a pair of ions. Depending on the energies and trajectories of the ions emerging from the collision one or both of them might be able to travel a significant distance before experiencing a further collision and a subcascade can form. However, an ion which has been significantly slowed, either by the collision or by energy loss to the electrons, and has a kinetic energy in the keV range, will have a large cross-section for interaction with other ions and will undergo collisions with every ion that it encounters. At this stage a displacement spike will form.
Before we discuss the evolution of the displacement spike in detail, we will consider the important role played by electrons in the processes discussed above. At the high ion velocities involved in channelling, kinetic energy is lost mainly to the electrons3, whilst the ions surrounding the channel are disturbed only slightly by their fleeting interaction with the passing projectile. However, experimental results show significant damage, referred to as an ion track, to the material surrounding the path of a channelling ion. The transfer of the energy required to displace these ions must be mediated by the electronic subsystem. The mechanism involved is the subject of some controversy. Two plausible models have been put forward.
The first is the so-called thermal spike model in which the electrons in the vicinity of the ion track undergo a large degree of excitation and subsequently lose energy to the ions in the same region, causing local melting. The plausibility of the thermal spike model depends on the competition between the rate at which the energy is conducted away from the track region within the electronic subsystem and the rate at which it is transferred to the ions. A model, due to Duffy et al, that is capable of investigating this balance [8,9] is discussed in section 4.5 on page 69.
The second, alternative view is the Coulomb explosion model in which the channelling ion is thought to ballistically eject electrons from the adjacent ions, causing a build-up of spatial charge in the track region [10]. It is the Coulombic repulsion between the ions that is responsible for the damage. The plausibility of a Coulomb explosion depends on the
3 The stopping power theories discussed in section 2 suggest that for a 500 keV Fe ion in iron, ~70% of the energy loss is to electrons.
ability of the target material to sustain a spatial charge distribution for time scales comparable to those required for ionic rearrangement to take place.
Ryazanov etal [11] have compared calculations of heating in copper for the two models of damage formation and found that only the Coulomb explosion can give sufficient heating to yield a molten region around the track. Fleischer et al [12] have shown that the extent of radiation damage in dielectrics correlates better with the rate of ionization by the PKA than with its rate of energy loss, again lending support to the Coulomb explosion model. Itoh and Stoneham [13] point to evidence from work on dichalcogenides that shows track formation does not occur for electrical conductivity exceeding 105 ^-1cm-1.
A third possibility, not widely considered in the literature, mirrors the thermal spike model, but without the requirement that the excitation energy of the electrons be passed to the ions in the form of thermal excitation. In fact, as discussed in section 5.4 on page 89 and in [14], a highly excited electronic subsystem implies significant weakening of the attractive bonding interaction between ions and the resulting outward pressure on the surrounding lattice may well be sufficient (and sufficiently long-lived) to provide a mechanism for damage formation. Distinguishing between these models requires a quantum mechanical model of electronic stopping and electron-phonon coupling.
Now, returning to our discussion of the evolution of a collision cascade, we will consider the displacement spike [15], first of all from the point of view purely of the interionic interactions. We can regard the displacement spike as being initiated by an energetic ion, moving sufficiently slowly that it interacts strongly with the surrounding ions, intruding into an undisturbed region of the target material, which for the sake of example we will take to be a crystalline lattice. A sequence of collisions takes place over a time scale of 1-10 ps in which the majority of ions over a region of 10-100 nm in size are displaced from their equilibrium lattice sites. This initial period of disruption is known as the displacement phase of the displacement spike and is followed by a relaxation phase during which the energy is rapidly repartitioned amongst the ions to yield a hot (and potentially molten) region, often referred to as a thermal spike. Finally, there follows a cooling phase in which the excited region grows and cools. By the end of the cooling phase, several hundred picoseconds after the initial PKA impact, the lattice will have healed itself to a large extent, many of the interstitial and vacancy defects formed in the displacement phase having recombined, and some final damage state will have formed. The nature of the damage will depend on the target material and the energy of the cascade and will be characterized by some particular population of defect types showing a particular tendency towards clustering. Broadly speaking, we tend to observe a core region with an elevated concentration of vacancies (the depleted zone) surrounded by a shell containing interstitial defects.
An important mechanism for damage production is the replacement collision sequence (RCS) in which an energetic ion collides with another ion at a low angle of incidence to a line of atoms in the crystal. A sequence of collisions ensues,
with each ion replacing the next along the line, carrying the resulting interstitial atom far enough from the corresponding vacancy to inhibit recombination. Successful formation of an RCS is highly dependent on the alignment of the ions and is thus less likely at high temperatures. It also requires a low ionic kinetic energy and even an RCS along a close-packed direction will require energies of less than 100 eV [16]. The lattice surrounding an RCS will tend to 'breathe' slightly as the disturbance passes, helping to steer the collisions and also gradually slowing their progress.
We expect that the electrons should play an important role in the evolution of the displacement spike. When the spike first forms, all of the excess energy is contained in the ionic subsystem and so although at low velocities (below ~500keV) the dominant mechanism for energy loss from an ion is to other ions, we expect a net transfer of energy into the electronic subsystem during the displacement and relaxation phases. During the cooling phase the electrons will function as a thermal bath in which the final defect distribution establishes itself. The dynamics of energy exchange with the ions, characterized by the electron-phonon coupling, as well as the high electronic thermal conductivity, will determine the effect of the electrons in quenching in, annealing out or inhibiting the production of defects.
1.2. Transfer of energy between ions and electrons in radiation damage
A key aim of the simulation of radiation damage collision cascades in metals is to establish analytical tools and methods of simulation able to predict the damage distribution at the end of the cooling phase. Experimental testing of the effects of 20 years of exposure in a high flux environment is difficult to achieve. The defect populations caused by irradiation constitute the initial conditions for the long-term microstructural evolution of the metal. The models reviewed here form the first level in a multiscale hierarchical description of radiation damage spanning time and length scales from the electronic to the geological.
Fundamentally important in the history of the development of the theory and simulation of radiation damage has been the evolution of the means by which energy exchange, between ions and between ions and electrons, is treated. The earliest theories undertook a wholly statistical treatment of cascade development. Collisions between particles were treated within various approximations and with increasing sophistication, ultimately yielding a body of theory able to make strikingly accurate predictions of particle implantation ranges and defect concentrations.
With the increasing availability of computer resources, work began on the direct simulation of cascade evolution. The earliest efforts married the established theory with Monte Carlo techniques to gain more accurate information about the results of radiation damage (see section 4.3). In such simulations, an explicit model of the ionic distribution makes its first appearance, although the ion-ion interactions are treated in a simplified way. More recent work has used molecular dynamics (MD) to give a completely explicit picture of ions
evolving under Newton's laws, under the influence of various force models. However, electrons within the metal have until recently been treated only implicitly, in the simplest case via their contribution to the interionic potential and in the most sophisticated case additionally as a viscous medium and energy sink.
Given the important role of electrons in providing a heat sink or reservoir and a means of energy redistribution, and given the likelihood that in many situations the electronic contribution to the interionic forces will violate the Born-Oppenheimer approximation implicitly assumed in the potentials of classical MD, we might reasonably conclude that the next stage in the evolving effort to model radiation damage should focus on improvements in the treatment of electrons. Some of the most recent work begins to address this challenge.
The full dynamics of a system of interacting nuclei and electrons are encapsulated within the time-dependent Schrodinger equation for the system. The solution of this equation being intractable for large systems, a given simulation method can be viewed as an approximation, to some lesser or greater extent, of the full many-body quantum mechanical dynamics. Broadly speaking, there are then two routes to incorporating electronic effects within such a simulation. One possibility is for the degree of approximation to be such as to retain the physics necessary to give rise to the electronic effects as a direct consequence of the dynamics of the model. That is to say that the model includes an explicit description of quantum mechanical electrons. Examples of such simulation methods are discussed in sections 5.3-5.5. Alternatively, the approximation to the system dynamics may be such as to exclude electronic effects, which are then added back in via some augmentation of the model intended to capture particular phenomena. Examples of this class of model are covered in sections 4 and 5.1.
The large body of theoretical and experimental work dedicated to understanding and measuring the effects of energy exchange between ions and electrons in radiation damage is thus highly relevant, either in evaluating models of the first class or in informing the phenomenological additions to models of the second class. In sections 2 and 3 we will undertake a brief review of the most significant material.
Traditionally, theorists have divided the process of energy transfer between ions and electrons into two regimes: the electronic stopping regime and the electron-phonon coupling regime. These regimes can be defined in terms of the predominant mode of ionic motion. In the electronic stopping regime, the ions are assumed to move ballistically, undergoing collisions with one another. In the electron-phonon coupling regime, the ions are assumed to oscillate around their equilibrium positions. It is important to note that the boundary between these two modes of behaviour is ill-defined and that the physics of energy exchange is the same in both cases [17]. In fact, the distinction is largely a practical one, with each regime having yielded a different theoretical treatment: for instance, extrapolation of electronic stopping theory down to ionic energies approaching those associated with phonon-like behaviour is unsuccessful.
The electronic stopping regime is characterized by high ionic energies, so that the dominant processes lead to a net
transfer of energy from ions to electrons. Ubiquitous within the literature concerned with such processes is the concept of electronic stopping power, defined as the rate of energy loss to the electrons along an ion's path and representable as a retarding force on the ion's motion. At high ionic velocities, these forces can significantly affect ion trajectories and in all events will act to dampen the evolution of an ionic system. Various models that aim to predict the rate of this energy loss are considered in section 2.
In the electron-phonon coupling regime, when ionic energies are much lower, energy transfer in both directions between electrons and ions will be significant and the electronic system will function predominantly as a heat reservoir. Theoretical models of electron-phonon coupling are considered in section 3.
2. Models of electronic stopping power
The problem of how a ballistic particle loses energy to the material through which it moves has attracted the attention of many researchers over the last century. The sheer volume of literature generated in a field whose importance spans from fundamental physics research to large-scale industrial application precludes a truly comprehensive treatment in any review paper4. Here we will give only a brief treatment, focussing on those aspects of the theory that have a bearing on atomistic simulations of radiation damage, either in helping us to understand the physics of the processes at work or in informing the construction of atomistic models.
We will begin by considering an energetic particle penetrating a stopping medium, assumed to consist of a collection of target particles. We will restrict our discussion mainly to the case of a solid target made up of ions with electrons in bound states and a gas of valence electrons. The penetrating particle, the projectile can lose energy via mechanisms falling into five categories:
(i) changes in the internal state of the target ions (electronic excitation and ionization) or excitations of the electron gas,
(ii) changes in the internal state of the projectile (electronic excitation, ionization and electron capture),
(iii) transfer of energy to the motion of the target ions (or to the generation of phonons),
(iv) emission of radiation (e.g. Bremsstrahlung and Cerenkov radiation), and
(v) chemical or nuclear reactions.
The general picture is clearly very complicated and efforts to understand the energy loss of a projectile are classified based on which mechanisms are significant. A standard classification scheme in the literature [20] divides projectiles by their atomic number Z1 into light (Z1 < 2), heavier or intermediate (3 < Z1 < 18) and heavy (Z1 > 19) ions. Projectiles are further classified by their kinetic energy per atomic mass unit E/W (effectively a measure of the velocity squared) into
4 Many excellent sources of further information exist, including the work of Bohr [4], Northcliffe [18], Seitz and Kohler [19], Ziegler et al [5] and Sigmund [6].
in" 10
Figure 1. The classification of electronic stopping behaviour into various regimes based on projectile atomic number and kinetic energy per atomic mass unit. A sample of applications is indicated on the chart according to their corresponding regimes.
fast (E/W > 10MeV), intermediate (100keV < E/W < 10MeV) and slow (1keV < E/W < 100keV). This classification scheme is illustrated in figure 1 along with an indication of to which region various applications of the theory correspond.
For atomistic simulations, the key concept in the literature on the energy loss of particles in matter is that of the stopping power, S, defined as the rate of loss of projectile kinetic energy E per unit length x along its path,
dE dx '
In fact, S has the dimensions of a force, but whilst the term 'stopping force' is now finding its way into the literature, 'stopping power' remains the dominant name. Because the energy transferred into the stopping medium can be attributed partly to the centre of mass motion of the target ions and partly to electronic excitations, it is common to split the stopping power into corresponding nuclear (n) and electronic (e) parts,
Each term will behave differently as a function of projectile velocity, with nuclear stopping dominating at low speeds and electronic stopping at higher speeds. This behaviour is most directly illuminated if we consider a target medium made up of free electrons or nuclei and an intruding point charge projectile of mass m1 and charge q1 moving with speed v. The energy transferred to a stationary target particle of mass m2 and charge q2 in a collision with impact parameter b is given by the Rutherford formula5,
2ffi2ff22
m2V2b2 \ 1 + (q1qi/^v2b)2
where the reduced mass ^ = m1m2/(m1 + m2) and where b is defined as the initial distance between the projectile and target
5 In all appearances of Coulomb's law, Gaussian units will be used.
keV/amu
S — Sn + Se
perpendicular to the projectile velocity. If we assume that the projectile path is unperturbed and the target particle does not move during the collision then (3) reduces to
2q2 gf v2b2
Assuming a target particle number density n and integrating over the range of valid impact parameters gives a stopping power,
S(v) = 2nn
4nqjq2 , bmax
-n ln ■
dbbT(b)
2.1. Models of the stopping of fast, light particles
The earliest theories of electronic stopping, from the 1910s, were due to Thomson [21] and Darwin [22]. They treated a point charge projectile interacting with free electrons of mass me in a target medium to arrive via the energy transfer formula (4) at a stopping power,
4nZ2 e4 S(ZU V) = .....2 Z2naLfree,
me v 1
Lfree = 2 ln
At large velocities, the behaviour of this stopping power will be dominated by the prefactor to the logarithm and so the fraction 1 /m2 will determine the relative contribution of nuclear and electronic energy losses. This result, that the electronic stopping power is the dominant process at high projectile velocity, remains generally true in more complex treatments.
Because atomistic models treat energy exchange between a projectile and the nuclei automatically, via some explicit model for the internuclear or interionic interaction, and because our key concern in this review is the treatment of electronic excitations, we will henceforth focus our discussion on theories of electronic stopping.
Over the decades, many different approaches to determining the electronic stopping power have been tried: classical and quantum mechanical, treating the target medium as a continuum and in a binary collision approximation, and varying between a first principles basis and completely empirical fitting. The aim is to produce models that can accurately predict electronic stopping powers for arbitrary combinations of projectile and target over given energy ranges. Which approach works best depends on the type and velocity of the projectile and the nature of the stopping medium.
Broadly speaking, the physics is at its simplest, and theory at its most successful (and abundant), in the case of fast, light projectiles. In this situation, the projectile can be regarded as a point charge, stripped of all electrons and losing energy only to excitation of target electrons and, at relativistic velocities, via radiative processes. For heavier particles it becomes necessary to consider bound electronic states and the additional energy loss mechanisms that appear as a result. And for slow particles, matters become even more complicated. Screening of the interaction of the projectile, now a complex compound object, with the target is significant and few simplifying assumptions can be made.
Atomistic simulations of radiation damage in metals cover all of the above regimes, dealing with channelling and sputtering of high energy incident particles and also with the evolution of displacement cascades down to thermal ionic velocities. We will therefore give a brief overview of examples from the whole range of stopping theory, beginning with the simplest case and discussing complexities as they arise.
where na is the number density of target atoms and Z2 is their atomic number. Quantities of the type Lfree, following the commonly occurring prefactor in (6), are known as stopping numbers. Application of the formula requires limits on the energy transfer, Tmax and Tmin. The maximum energy transfer is naturally set by considering a 'head-on' collision,
4«1 me 1
---m1v
(«1 + me)2 2
but in order to prevent divergence of S due to the long range nature of the Coulomb interaction, we also need to set a minimum energy transfer Tmin. Both Thomson and Darwin relied on an artificial limit corresponding to some maximum impact parameter (e.g. the atomic radius in Darwin's case).
This deficiency can be removed by taking into account the interaction of target electrons with target nuclei. Theories due to Bohr [23] andBethe [24] do just this and continue to be used (with various refinements) up to the present day.
2.1.1. Bohr formula. Bohr's theory of electronic stopping considers the energy transfer from a point charge projectile to classical electrons harmonically bound to the target atoms with angular frequencies wj. The interaction of the projectile with the electrons is treated under the assumption that the strength of the Coulomb interaction does not vary significantly over the range of motion of the electron. This is known as the dipole approximation because it amounts to neglecting all but the leading order term in electron position in the Fourier transform of the electric field at the electron—a long wavelength limit. Equivalently we can calculate the response of the electron to a force whose strength is calculated under the assumption that the projectile path is undisturbed and the electron remains at rest throughout the collision (the momentum approximation). Bohr's theory is thus a classical perturbation treatment and so is appropriate only for weak interactions. The final result for the stopping power is
4nZ2 e4 S(Z1, v) =-— naLBohr,
El Cmev3
where the constant C = 1.1229 and the values of fj give
the relative contributions of different frequencies &>j subject to j = 1.
Several important concepts arise in the derivation of Bohr's theory. The first is that when considering the binding of electrons, a natural upper limit to the impact parameter emerges, countering the long range nature of the Coulomb interaction. This corresponds to the situation in which the collision takes place so slowly that the harmonic electron moves appreciably during the interaction and no energy is transferred. This maximum impact parameter, bmax, then corresponds to an adiabatic radius,
where b/v is a measure of the time scale of the collision (the collision time). We can also recognize this limit as the distance at which the Coulomb interaction becomes significant compared with the binding forces on the electron.
Second, the energy transfer to a bound electron in the dipole approximation diverges for small impact parameters. Bohr removed this divergence by recognizing that for close collisions, such that the collision time is much smaller than the period of electronic motion,
b 1 - « —.
T(b) =
2Z2e4 mev2b2
1 + (Z1e2/mev2b)2
Z1e2 nev2 b*
— « —,
and so we have a validity condition,
Third, the treatment is classical and so relies on being able to describe the projectile as a well-confined wave packet throughout the collision. A well-collimated projectile beam with a spread in transverse momentum Spi will have a corresponding uncertainty in impact parameter Sb ~ h/2Spi. This will imply a spread in the transverse component of the momentum transferred in the collision of Sp2 ~ (2\Zie2\/b2v)Sb (from (4)). Minimizing ((Spi)2 + (Sp2)2)i/2 as a function of Sb and assuming that for the classical approximation to be valid this uncertainty in the transverse momentum must be much smaller than the total momentum transfer yields the condition
2|Z1 e2| h v
Thus, Bohr's classical treatment seems to become more valid with decreasing velocity below some high threshold. Unfortunately, the assumption is also made that the target electron is at rest during the collision, i.e.
v » vo,
where v0 = e2 /h = c/137 is the Bohr velocity. Combining these two criteria,
the effect of binding can be ignored. By treating collisions below some threshold impact parameter b* as being Coulomb collisions between free particles, Bohr avoided imposing an artificial lower limit on b.
However, this derivation implies a restriction on the applicability of the theory. To achieve a smooth join between the two limiting treatments of the collisions, interactions with bound electrons in the dipole approximation at large b and Coulomb interactions between free particles at small b, both must be valid at b*. The full Rutherford formula for the energy transfer in a collision between free particles is in this case
v0 « v « 2Z1 v0,
where the term Zie2/mev2b accounts for the effect of the deviation of the paths of the projectile and the position of the electron under the Coulomb force, i.e. the deviation from the assumptions used in the treatment of distant collisions, which must be small at b*. Hence,
For the close collisions to be unaffected by the binding forces, the collision time must be sufficiently short:
Z1e2Mj
we find that the Bohr theory should only be valid over a small range of high velocities above v0 for a heavy projectile ion.
The Bohr formula exhibits several features characteristic of experimental data for stopping powers (see figure 2). At high energies the stopping power drops away as i/v2. This is a consequence of the fundamental physics captured in the Rutherford scattering formula. The impulse on a target electron will be proportional to the collision time b/v and so the energy transfer will vary as the square of this. The Bohr formula also exhibits a peak in the stopping power at approximately the correct energy. This is known as the Bragg peak and its location is important experimentally because the depth resolution of experimental probes is maximized along with the stopping power at this energy.
At energies below the Bragg peak, the Bohr stopping power falls off, though much more rapidly than the experimental data. This fall-off is due to the effect of the binding of the electrons. At low velocities, the adiabatic radius v/wj becomes small and most collisions no longer transfer any energy.
2.1.2. Bethe formula. Bethe6 derived an alternative stopping power formula using quantum mechanical perturbation theory. In his derivation, the Coulomb interaction of an incoming particle, represented as a plane wave, with a target atom is treated in the first Born approximation [27] to give a stopping number,
'2me v2\
^Bethe = fj ln j
6 The original derivation is in German [24], but a thorough English language treatment can be found in Sigmund [6].
2mev2=h®o
Bohr Bethe Experimental
v=2ZiVo
10 10 10 10 10 10" Projectile Energy (keV/amu)
Figure 2. Sample electronic stopping power data for oxygen projectiles in a gold target. Experimental data (red crosses) from the database of Paul [25] are shown, along with stopping powers calculated via the Bohr (green long dashed line) and Bethe (blue short dashed line) theories. Various velocity thresholds discussed in the main text are indicated. The average excitation energy in the theoretical expressions, ln I = ^j fj ln(hwj), is calculated using a commonly used scaling relation I & Z2 x (10 eV) [26].
where wj is the angular frequency corresponding to the jth excitation of electrons in the target atom. fj is a generalized oscillator strength [27] in the limit of low momentum transfer and is proportional to the square of the matrix element of the projectile potential coupling the initial and final states in the j th excitation.
The derivation of the Bethe formula again relies on a split between close and distant collisions, in this case made at a particular momentum transfer hq*. If we consider the incoming particles as a plane wave with wave-vector k being scattered into a state with wave-vector k' whilst the target atom undergoes an excitation of energy hw0, where w0 is a typical excitation frequency, then for energy conservation we require
h «o +
■ (krl — k2) = 0.
h «0 +
h2q 2 2m 1
— hq ■ v = 0.
< h«0.
\h «0 2me
An example of the behaviour of the Bethe theory is shown in figure 2.
2.1.3. Corrections to the Bohr and Bethe Formulae. For most of the key simplifying assumptions used in the derivation of the Bohr (8) and Bethe (18) stopping numbers, there exist refinements and correcting terms. For a detailed discussion of these corrections see Ziegler's review [28].
Bloch derived a revised model for the electronic stopping number that tends to the Bohr model at lower velocity and the Bethe model at higher velocity. Bloch's model can be regarded as a correction to the Bethe stopping number to account for the failure of the Born approximation to correctly treat close collisions [29] (which are effectively Coulomb collisions between free particles, unaffected by electronic binding forces over the short time scale of the collision). These collisions are more important at lower projectile velocities when the effect of more distant collisions is removed by a smaller adiabatic radius. Conversely, the Bohr theory explicitly treats close collisions as free-Coulomb interactions, in which limit the classical and quantum mechanical treatments coincide, but fails because of its classical nature to give a good account of binding in distant collisions. The Bloch theory can thus also be seen as a correction to the Bohr theory in respect of more distant collisions and we can write the Bloch stopping number in two forms [6],
- 1.202
LBohr —
2 V Z1V0
The momentum transfer hq is given by q = k — k' and so, writing the projectile velocity v = hk/mi, we have
q * must clearly be larger than the minimum q at which (20) can be satisfied, so
q> W0. (21)
In order to apply the dipole approximation for distant collisions, the kinetic energy transferred to an electron must be small compared with a typical excitation frequency,
h2(q *)2
Together conditions (21) and (22) restrict the validity of the Bethe theory to the velocity range,
in the limit of high v and low Z1e.
The Bethe and Bohr theories both assume that the target electron is initially at rest. Corrections to account for electronic motion have been given for the Bohr model [30] and the Bethe model [31] and are known as shell corrections. A Thomas-Fermi model [32] for the target atom suggests that typical electron velocities will be Z^/3v0, giving a threshold for the importance of shell effects, v < Z^3vo.
A difference in the ranges of positive and negative pions observed by Smith et al [33] led to explorations of terms in the stopping cross-section of order Z3 and above, corresponding to higher order terms in the perturbation expansions in the Bethe and Bohr theories. Further observations by Barkas et al [34] and Andersen [35] led to such corrections being called the Barkas or Barkas-Andersen effect. The Barkas effect is also referred to as the polarization effect because the underlying cause of the difference in stopping between positively and negatively charged projectiles is the distortion of the electron density in the target ions. This alters the electron density experienced by the projectile, giving an enhanced stopping power for positively charged particles. Ashley et al [36] produced a theoretical treatment for the Bohr model and a simple treatment by Lindhard [37] gives an estimate of the size of the effect in an amended Bethe stopping number,
4^Z?e4 / 2Z1e2«o\ 2mev2
1 I 1 + 1 0 1e
in" 10
102 103 keV/amu
Figure 3. The projectile kinetic energy and atomic number regimes of electronic stopping theory in an iron target (Z2 = 26) showing order of magnitude thresholds for various effects and corrections. Bohr's classical threshold v < 2Zj v0 is shown, along with the thresholds for screening (v < Zi/iv0) and Barkas (polarization) (v < (Zj Z2)1/3v0) effects. The velocity at which shell effects become important (v < Z^/3 v0) and the Bohr velocity, v0, below which the projectile ion will have very low charge with many bound states are also indicated. (After Sigmund [39].)
where the second term in parentheses is known as the Barkas-Andersen parameter and takes the same form as the dimensionless parameter in the validity criterion (14) for the Bohr formula (8). I is a logarithmic average excitation energy
ln I = T.jfj ln(hj
Figure 3 shows the thresholds in projectile charge and kinetic energy per atomic mass unit at which various factors in the stopping power become important. Grüner et al [38] find that for a 1 MeV amu-1 nickel ion in a carbon target, the Bohr theory gives a stopping power 42% smaller than experimental results. A Barkas term over-corrects to give a number 28% too large and further addition of a shell correction gives a result 14% too large. In this case and in many others, the correction terms are significant and must be taken into account to obtain good predictions of stopping power.
2.2. Models of the stopping of fast, heavy particles
We have seen in section 2.1.3 that relaxation of the simplifying assumptions in the theories of Bohr and Bethe necessitates a series of corrections of considerable complexity. It is important to note that this complexity arises even in the simplest case where we treat projectile particles not only as point charges, but also as having very low charge (Zi < 2). The theories of Bohr, Bethe and Bloch are essentially perturbative, going to second order in the projectile charge Z1. Even with higher order corrections such as the Barkas term in Z3 and the so-called Bloch correction in Z4, these theories rapidly become inadequate for the treatment of heavy ions (Z1 > 2) even at very high velocities. The increasing strength of the inter-particle interactions means that the fundamental assumption in any perturbation theory, that the unperturbed evolution (be it the projectile path or the initial quantum state of the projectile-target system) is a good approximation to the perturbed
evolution, is invalidated. A number of models developed in the last decade and designed to treat the stopping of heavy ions non-perturbatively will be discussed below.
The physics of heavy ion stopping is fundamentally different from that of light particles. A highly charged nucleus will attract to itself a charge-compensating cloud of the electrons of the target medium. There will thus be a screening of the interaction of the projectile and target particles and a reduced stopping power. The intruding nucleus may also carry with it its own electrons in bound states so that the projectile becomes a compound object. Internal excitations of the projectile ion will become possible, enabling new mechanisms for energy loss7. The charge of the projectile will no longer simply be Zi and may change as the projectile traverses the stopping medium. Charge-changing processes of electron capture and ionization open further mechanisms of energy exchange.
2.2.1. The effective charge of the projectile. The possibility of bound electronic states of a projectile particle gives rise to the concept of an effective charge of the intruding ion, written Zj e. It is tempting to model the stopping of heavy ions by simply replacing Z1 with Zj in simpler perturbative stopping theories. This approach is not without some success (see section 2.4), but it is not immediately persuasive from a physical point of view.
There is no reason to assume that the charge state of a projectile ion will be constant during its flight and though it may acquire a well-defined mean value in a steady state, this will be the result of repeated charge-changing processes. Assuming a fluctuating number of bound electrons Nbound, so that at any time the effective charge can be written Zj e — (Zi — ^bound)e, we will denote the mean effective charge by (Zj)e. The earliest considerations of effective charge adopted an empirical definition. Writing a stopping power S(Z1, Z2, v) dependent on the atomic numbers of the projectile and target and on the projectile velocity, Northcliffe [18,40] defined an empirical effective charge (Zj)emp
= S(Z1,Z2,v)/S(1,Z2,v),
in which the effective charge is related to the proton stopping power in the same target. This definition is informed by the Z2 dependence in the stopping equation (5). Yarlagadda et al [41] compared experimental stopping data for protons and for carbon and iodine projectiles in a selection of targets up to Z2 — 79 (Au) to show that the quantity (Zj)emp/Z1 is independent of Z2 to within 10%. This suggests some validity for the effective charge concept and that, depending on the accuracy that we require, the treatment of projectile and target may be separable.
Various theoretical definitions of effective charge have also been considered in the literature. These often take the form of stripping models in which electrons are considered to be progressively stripped from the projectile ion as a function
7 An enhancement of the stopping force attributable to the electrons associated with the projectile is sometimes referred to as an anti-screening effect.
of its velocity based on some criterion related to the orbital velocities or potential energies of the electrons. A Thomas-Fermi model of the atom, in which electron orbital velocities are assumed to be ve -stripping criterion [39],
Z\/3v0, gives rise to a much used
1 ) strip
= (1 - exp(-v/ZLnV0))Z1,
in which electrons of the target medium, approaching with a relative velocity ~v, are assumed to knock slower electrons out of states bound to the projectile nucleus. Yarlagadda et al [41] used a similar model and obtained good agreement with empirical stopping data for ions heavier than chlorine. The less good agreement for lighter ions is improved by including corrections in the scaling relation analogous to the Bloch correction and the Z3 Barkas term.
Brandt and Kitagawa [42] note that when a projectile nucleus carries with it some distribution of bound electronic charge, an electron of the target medium will experience a projectile charge that depends on its impact parameter with the projectile nucleus. Thus, even for a given net charge (Z1 - Nbound)e, the stopping power will still have a Z1 dependence, increasing with increasing atomic number. To quantify this effect, they calculate a stopping power for a projectile of atomic number Z1 and fixed ionization Q = (Z1 — Nbound)/Z1 using a Lenz-Jensen model [5] of the projectile charge density and a Lindhard type model (see section 2.3.2) with an approximate form for the dielectric function of a free-electron gas target of number density ne. Writing the stopping power as S(Z1, Q, ne), the effect of the charge distribution is then revealed in the quantity
IS(Z1, Q,n) S(Z1, 1,ne) .
If the charge distribution had no effect then we would expect to see the stopping power vary as Q2 (see (38)). If the target electrons penetrate the electron distribution of the projectile they will experience a higher than expected charge and the quantity (29) will be greater than 1. Experimental data for ions of boron through to fluorine of fixed charge in (111) channels in gold reveal the expected effect, that, for example, S(Ns+) > S(C5+) > S(B5+). The qualitative variation of the stopping power with Z1 and Q is well modelled by the theory of Brandt and Kitagawa [42]. In addition, experimental data for the stopping of nitrogen in gold, carbon and aluminium targets over a velocity range 0.7v0 < v < 1.5v0 show the ability of the theory, when combined with a simple stripping model for Q as a function of v, to capture the behaviour of the stopping power of light ions that proved problematic in the theory of Yarlagadda et al [41]. However, the resolution of the data is not good enough to reveal if variations due to the target element (modelled via the free-electron gas density ne) are well captured.
The ultimate aim of effective charge theory is to provide estimates of heavy ion stopping powers. This is done either by multiplying some reference stopping power (normally taken to be that of a proton or alpha particle in the same target medium, for which comprehensive data are available) by some velocity
dependent effective charge or using the effective charge as an input to a theoretical model of stopping power. Attempts to formulate such methods will be hampered by the fact that the mean stopping power over the fluctuating charge state of a projectile ion traversing a stopping medium, (s(Q)), will not in general be equal to the stopping power at the mean charge, s((Q)). More importantly, as pointed out by Sigmund [43], success depends on the stopping of heavy ions and the stopping of the projectile in the reference data set being governed by the same physical processes. This is unlikely to be the case; we have no reason to expect that a singly ionized gold ion will behave in the same way as a proton. Experimental validation of effective charge theories [41,42] tends to focus on the high velocity regime, where projectiles will be highly ionized and poorly screened. Such almost bare particle stopping is most susceptible to an effective charge treatment, but at lower velocities the effects of screening and bound states will become significant. As Sigmund discusses [39], any attempt to incorporate these effects in an effective charge model is likely to result in a model so clumsy that a direct calculation of heavy ion stopping powers, ignoring any reference data set, will be at least as straightforward.
2.2.2. Non-perturbative models of heavy ion stopping. As discussed above, a key difference between light and heavy ion stopping is the strength of the interaction between the projectile and the electrons of the target medium. In the latter case, the use of perturbation theories such as Bohr's and Bethe's is questionable and new approaches are needed.
Over the last decade, many researchers have been active in developing stopping power theories applicable to fast, heavy projectiles. Various models have arisen, which show good agreement with experimental data, and there is a significant volume of relevant literature. Again, we will not attempt a comprehensive survey and the interested reader should consult the discussion and references in Sigmund's book [6]. Instead we will discuss the general character of a sample of theories of fast, heavy ion stopping, considering the inputs upon which they rely and the results they produce and briefly discussing the way that they work.
The broad aim of all the models that we will consider is to take a set of input data, derived either from experiment or from other theoretical calculations, and from this derive a prediction of electronic stopping power as a function of projectile and target type and projectile kinetic energy. They aim to achieve this without the use of any adjustable parameters. Beyond these common traits, the models that have been developed are highly diverse; some are quantum mechanical, others classical; some are based on straight calculations, others require a set of simulations to be performed.
The binary theory of Sigmund and Schinner [44, 45] adopts Bohr's classical approach to calculating the scattering of bound target electrons by a projectile. A non-perturbative treatment is made possible by exploiting an approximate equivalence between the interactions of a bare projectile with a bound electron and of a screened projectile with a free
electron8. Information concerning the binding frequency of the j th electron is then encoded in the screened potential,
V(r) = e-ra'/v,
where v/ojj can be identified as an adiabatic radius for the j th electron. Each electron is treated separately and the individual stopping cross-sections summed together. In common with other non-perturbative treatments, the binary theory treats the Barkas effect automatically. The screening effect of bound projectile electrons (as opposed to the screening in (30) used to represent the effect of the binding of electrons in target ions) with a radius ascr is incorporated into the potential as
V(r) = -
Abound e
-rœj/v
(Z1 - Abound)e
4 = ( ^ )2+4- •
a2 \ v / at.
have a contribution of -13% (an accelerating effect) due to polarization of the projectile electrons by the ionized target. Comparisons of simulations of a 1 MeV amu-1 Ni projectile in solid carbon give agreement with experiment for the steady state charge and steady state stopping power to within 2.8% and 3.7%, respectively.
Non-perturbative quantum mechanical schemes have also been proposed. Motivated by a desire for economical stopping calculations with the accuracy of full quantum mechanical results, Grande and Schiwietz [46] developed the perturbative convolution approximation (PCA). This calculates the energy transfer AE from a projectile to the electrons of a target atom as a function of impact parameter b via an integral over the electron density p of the target,
AE(b, v) =
d2r±K(b - r±,v)p(r±,z), (32)
Further refinements allow for quantum mechanical corrections at higher velocities (comparable to the Bloch correction to the Bohr formula (8)), for shell effects and for projectile excitation (treated approximately by repeating the calculation and exchanging the roles of projectile and target). Implementation of the basic binary theory requires data on the atomic binding frequencies of the target electrons and their relative strengths (occupations) fj, along with knowledge of the effective charge (Z1 - Nbound)e.
The full model compares well with experimental data (from the database of Paul [25], discussed in section 2.4) for a range of projectiles (3 < Z1 < 18) in N, Al, Ni and C targets over a range of energies from 1 keV amu-1 to 100 MeV amu-1.
A second classical scheme, due to Griiner et al [38], employs the classical trajectory Monte Carlo method to treat a small number of target nuclei and the projectile nucleus, along with their associated electrons. The evolution of the system of electrons and nuclei is calculated using classical equations of motion for a large number of statistically sampled starting conditions.
In contrast to the binary theory, the classical trajectory method does not require a charge state as input, the only input parameters required being the nuclear charges and the electron orbital binding energies and their occupations. The recovery of details of the projectile charge as an output from the simulation is a stated aim of the technique and recognizes the important influence of charge on the stopping power. This aspect of the model better reflects reality, in which projectiles tend to reach a steady state that is independent of the input charge state. In addition, because all the constituent particles are treated explicitly, information about the relative contribution of various processes to the slowing of the projectile nucleus can be recovered. In a simulation of a 1 MeV amu-1 Ni ion in a gaseous Ar target [38], 80% of the energy loss is attributable to target ionization, 12% to target excitation and 20% to electron capture by the projectile. Projectile excitation is found to
8 There is an intuitive sense to this equivalence, in that the primary effect of the electron binding is to reduce the energy transfer at larger impact parameters.
where the convolution integrals are carried out in cylindrical polar coordinates about a z-axis through the target nucleus and parallel to the projectile velocity. K(b,v) gives the energy transfer to an electron at impact parameter b from a projectile with velocity v. The form of K(b, v), given in [46], is such as to treat close collisions as free and distant collisions within the dipole approximation. Inputs required for the model are the projectile screening function, the electron density of the target and the oscillator frequencies and oscillator strengths of the target.
Schiwietz and Grande [47] have also developed an extension to the PCA called the unitary convolution approximation (UCA), which implements a Bloch-like correction and extends the applicability of the scheme to high Z\. Comparisons of AE as a function of b and of Z\ with intensive quantum mechanical calculations [48] show good and very good agreement, respectively. Data for the stopping of oxygen in Al and Si across an energy range from 0.1 to 100 MeV amu-1 show good agreement with experiment [49].
Key experimental studies by Blazevic et al [50,51] provide important data for testing of the non-perturbative stopping models described above. By separating out the initial charge states of Ne ions prior to passage through thin carbon films and measuring their final charge state and energy, it is possible to determine the cross-sections for charge-changing processes and the charge dependent stopping power of carbon for neon. Blazevic et al [50] compare their results for charge dependent stopping against calculations within the binary theory of Sigmund and Schinner and the UCA method of Schiwietz and Grande amongst others. The UCA agrees almost to within experimental error. The binary theory scales less strongly with projectile charge than the experimental results.
We have considered the above theories as examples of the sorts of models now being applied to calculate stopping powers of fast, heavy ions. None of these models is ab initio in character: they take information about the excitation spectrum of the target as input and focus on the problem of calculating the collision dynamics and how they will stimulate the given excitations. Most models (that of Grüner et al being the exception) also require further input in the form of an effective
charge of the projectile. The effective charge concept is ill-defined and, at least in principle, much physics could be secreted within an appropriate velocity and target dependent form for an effective charge function. Add in the inherent complexity of many of the modern stopping models and it becomes clear that they are theories in the sense that they provide testable predictions rather than models designed to yield physical insight. As such, their success should be evaluated in comparison with that of various empirical fitting and interpolation models to be discussed below (see section 2.4).
Data on the performance of the models against experimental results demonstrate that the models can successfully calculate electronic stopping powers across an energy range from tens of keV amu—1 up to tens of MeV amu—1 for projectile ions in the 'light' and 'heavier' ranges (as defined in figure 1) in a wide variety of targets. The predictions are good to within experimental error, which can be as good as 2% at 50 MeV amu—1, but as bad as 20% at 5 keV amu—1 [20].
Given the current limits to the precision of experimental data, and given the nature of the stopping models as predictive tools, the problem of predicting stopping powers of fast and intermediate velocity ions could be considered to be solved, at least for the time being. The practical justification for further refining stopping power predictions for fast ions is not immediately clear. In fact, the applications shown in figure 1 suggest a need for more work on the stopping of slower and heavier projectiles, where the physics is more complicated and the theoretical and experimental literature much more sparse.
2.3. Models of the stopping of slow, heavy particles
The non-perturbative stopping models aim to provide reliable predictions of stopping power in cases where the projectile charge is large and the effects of screening must be taken into account. A key issue is that of the degree of ionization of the projectile ion and all the models except that of Gruiner et al require this as an input parameter. Whilst this is less of a difficulty at high velocities, where stripping models of effective charge should work reasonably well, as projectile velocities fall below the Bohr velocity, v0 = e2/h, we see the appearance of so-called Z1-structure. These oscillations in the behaviour of the stopping power with Z1 are due to the atomic structure of the projectile and demand a different theoretical treatment.
In this section we will examine models designed to predict the stopping powers of slow, heavy ions, with kinetic energies significantly below 1 MeV amu—1. We will begin by outlining three older models, which are still very much in use today, before going on to consider some more recent treatments. The older models fall into two categories: those which consider the energy loss as arising from the inelastic interaction of the projectile and a target atom in binary collisions and those which consider losses into a continuous electronic system.
2.3.1. Binary models of slow particle stopping. Into the first category fall two much used models due to Firsov [52] and Lindhard and Scharff [53]. Firsov ascribed the energy loss
from a projectile during a binary collision with a target atom to the exchange of electrons between the atoms. During the collision, the atoms are considered to form a quasimolecule. At any given separation, the motion of electrons within the quasimolecule will cause some electrons previously associated with the target atom to become associated with the projectile. This change of identity of an electron will be accompanied by a change in its momentum, proportional to the relative velocity of the two atoms. Using a Thomas-Fermi model for the colliding atoms, Firsov [52] obtained the following expression for the energy loss in a collision between a stationary target and a projectile moving with velocity v,
AE = 0.35
(Z1 + Z2)5/3
a0 1+0.16(Z1 + Z2)1/3rmin/fl0'
where rmin is the distance of closest approach and a0 is the Bohr radius h2/mee2.
Lindhard and Scharff [53] suggested an alternative formula, again based on a Thomas-Fermi model, but did not publish a derivation. Tilinin [54] has derived a more general result by considering the scattering of the electrons of the target atom by the screened field of the projectile. Use of the Thomas-Fermi model to obtain quantitative results gives an electronic stopping power,
8nnee a00Z1Z2
(zf + z22/3)—3/2
x(E,Z1/Z2)-, v0
where v0 is the Bohr velocity (v0 = e2/h) and the function t is the result of an integral over the electronic densities experienced during the collision. In Lindhard and Scharff's original formula [53], t is replaced by a constant ~ 1 — 2, used to obtain an improved fit to empirical data. Indeed, provided v is not too small, then t ~ 1 and Tilinin's formula (34) approaches that of Lindhard and Scharff. The key feature of (34) is the proportionality of the stopping power to projectile velocity, a feature shared with the stopping power implied by Firsov's formula (33).
Because the models of Firsov and Lindhard and Scharff both rely on a simple Thomas-Fermi model for atomic structure, they predict that stopping power will have a simple monotonic dependence on the atomic numbers of projectile and target atoms. They thus neglect both Z1-structure and an additional fluctuation in stopping power with target atomic number known as Z2-structure.
2.3.2. Electron gas models of slow particle stopping. The second class of models of slow particle stopping treats the electrons of the target medium as an electron gas. Fermi and Teller [55] gave an early treatment of the stopping power of a free-electron gas, pointing out that, since the maximum energy transfer from the projectile to a target electron corresponds to the case of a head-on collision, only electrons with velocities within v of the Fermi velocity vF can take part in stopping. Their final result is
2me e4 v S = —-—^— ln 3nh3
again giving a stopping power proportional to the projectile velocity.
Lindhard [56] gave a more general treatment by considering the force acting on a projectile due to the change in the electronic density distribution caused by the electric field of the projectile. The response of the electronic system is assumed to be characterized by a frequency and wavevector dependent dielectric constant e(q, m). This corresponds to allowing the total potential t) at any point x and time t to depend on the potential due to the projectile ^Proj(x', t') at points x' and times t', via the integral equation
f' dt /
J— 'X> J
dx'e 1(x-x',t-t/)0Proj(x',t/), (36)
e(q, œ)
0Proj (q, œ).
2e2 fqv f œ dœ
v2 Jo Jo
q " [ e(q, œ)
™ 1.0
—;-0 -B..............B
.................1.................m
250 œ T3
where e 1 is a linear operator. In Fourier space, this relationship becomes
The history dependence introduced by (36) allows for a finite response time of the electron gas to the potential due to the projectile. The centre of the screening cloud around an intruding charge will tend to lag behind, giving rise to a retarding force on the charge9. Lindhard [56] derived an expression for this retarding force,
for a general stopping medium characterized by e(q, a>) and where indicates the taking of the imaginary part.
Various alternative derivations of the Lindhard result have appeared in the literature. In particular, Ritchie [58] treated the problem of finding the induced charge density directly within first-order perturbation theory, without the need to introduce a classical electric field. The results obtained are equivalent to Lindhard's results using the Lindhard expression for the dielectric function of a free-electron gas [56].
Lindhard [56] gives results for the limiting cases of low and high projectile velocities. These results are also derived by Lindhard and Winther [59] alongside further discussion of the nature of the excitations of the electron gas. At high but non-relativistic velocities, the stopping power reduces to the Bethe formula (18). At low velocities the Fermi-Teller formula (35) is recovered.
2.3.3. Non-linear calculations of electron gas stopping. At electron densities typical of metals, and at lower projectile velocities, the perturbative approach behind Lindhard's stopping power formula is no longer appropriate and a nonlinear theory must be used. Various approaches have been tried and the review by Echenique et al [57] contains detailed discussions of both linear and non-linear treatments of stopping by an electron gas.
9 Note that at higher velocities v > V0, this simple picture of a retarded screening response is greatly complicated by the appearance of strong oscillations in the induced charge density. A full discussion of these so-called wake effects can be found in the review by Echenique et al [57].
Figure 4. Calculated results for a proton in a free-electron gas of varying density (parametrized by the one-electron radius rs = (3/4nne)1/3). The screening length is indicated for linear response (orange vertical crosses) and non-linear DFT calculations (red diagonal crosses). Distances are given in atomic units, a0 = 0.529 A. On the right-hand axis are plotted the ratio of the charge density at the proton position to the background density for linear response (purple squares) and non-linear DFT calculations (blue circles). (Data from Almbladh et al [60].)
The deficiency of a linear treatment was highlighted by Almbladh et al [60] who calculated the screening of a stationary proton in a free-electron gas using Kohn-Sham density functional theory (DFT) [61]. They compared DFT calculations of the screening length and the relative charge density at the proton position with linear calculations using the Kohn-Sham form of the dielectric function. The results are presented in figure 4 from which it can be seen that the linear theory underestimates the extent to which charge piles up around the proton and overestimates the variation in the screening length (which varies little in the non-linear calculation). At the lowest density considered, the results are similar to those for a 1s atomic orbital about the proton, indicating the importance of bound states at lower electron densities. One of the key features that emerges in non-linear treatments is the oscillatory Z1 -structure in the stopping power at low velocities.
Echenique et al [62] calculated stopping powers for ions with v < vF, the Fermi velocity, in an electron gas within time-dependent DFT (see section 5.3 for a brief discussion of the relevant theory). Figure 5 shows the behaviour of the stopping power for hydrogen and helium nuclei calculated in the nonlinear DFT and the linear response theory. In both cases the non-linear DFT calculations show a more rapid decrease in stopping power with decreasing electron gas density, consistent with the formation of bound atomic-like states that screen out interactions with the gas. Conversely, high electron gas densities screen the nuclei so efficiently that bound states can no longer develop and the linear and non-linear results converge. A final feature of note is that the stopping power for He lies below that for H at low electron gas densities. This is because the higher nuclear charge of He is more effective at producing bound states and thereby screening the stopping interaction. Such a feature can never emerge in a linear theory.
Z --r 1 1 II-
; "^S^B
H-iv \ ^
z >>> 1 i v. \ ^ -
Figure 5. Stopping powers as a function of the one-electron radius rs = (3/4nne)1/3. Results for linear response theory are shown for hydrogen (curve A) and helium (curve B). Results for non-linear DFT are shown by curves D (hydrogen) and E (helium). (Reprinted with permission from Echenique et al [62]. Copyright 1986 American Physical Society.)
Echenique et al used their calculated stopping powers along with the theoretical equivalent of the empirical definition (27) to calculate the effective charge on a projectile as a function of atomic number Z1 at several electron gas densities. Clear oscillations are present (see figure 6), which at low density follow the pattern expected from consideration of atomic structure, with deep troughs occurring at Z1 = 2, 10, 18, corresponding to stable filled shells. At higher densities, more effective screening means that higher nuclear charges are necessary in order to form bound states and the pattern shifts upwards in atomic number. Calculations for Z1 dependence of the stopping power of ions in (110) and (111) channels in silicon up to Z1 = 20 show good agreement with experiment.
Arista [63] has presented results for velocity and Z1-dependent stopping based on calculating the transport cross-section from the phase shifts for free electrons scattering off a screened projectile ion. A variety of screening functions are used and are adjusted in order to satisfy the condition that the total charge of the screening cloud cancels the residual charge (Z1 - Nbound)e of the ion. Calculations of the stopping power of bare ions for Z1 = 1 and Z1 = 7 show that, whereas for the lower charge the results of Arista's model are very similar to the linear dielectric results, for the higher charge the non-linear theory gives a significantly (60%) reduced stopping power. Again, this is typical of the overestimation of stopping power that results from perturbative treatments of highly charged projectiles.
Figure 6. The effective charge ZJ of point charge projectiles as a function of atomic number Z1 in a free-electron gas from non-linear DFT. Results for a variety of electron densities labelled by the one-electron radius rs = (3/4nne)1/3 are shown. (Reprinted with permission from Echenique et al [62]. Copyright 1986 American Physical Society.)
■ àiP~~~~ ^ooooooç, v=2
r ffl v=1°oo0°^ 1
.............. \v=0.2yr ............................
Figure 7. The Z1 -structure in the stopping power for bare ions in a free-electron gas of density equivalent to that in carbon at several projectile velocities (given in units of v0). (Reprinted with permission from Arista [63]. Copyright 2002 Elsevier.)
Arista's calculations of the Z1 dependence of stopping powers at various velocities show strong Z1-structure for v < v0, which is all but washed out with rising velocity by v & 2v0 (see figure 7). The positions of the minima shift upwards with increasing velocity as it becomes harder to form bound states. Comparison of calculations of stopping power for a carbon target at a projectile velocity of v = 0.8v0 with experimental data shows good qualitative agreement with the positions of the peaks and troughs in the Z1-structure up to Z1 = 40 and good quantitative agreement of stopping powers up to Z1 ~ 20.
2.3.4. Atomic environment dependence of electronic stopping. Beyond considerations of projectile and target type and projectile velocity, we might also expect the stopping power felt by a projectile to depend on the positions of the surrounding target ions. Electron gas models like those discussed above clearly exclude all such effects. Models that treat stopping as a consequence of energy transfer in binary encounters between the projectile and target ions can potentially yield some form of environment dependence through the impact parameter dependence of the energy transfer: indeed, the provision of such information is a stated aim for the binary theory and the UCA discussed in section 2.2.2. But these treatments will miss environmental dependence of a many-atom nature, such as we would expect to become increasingly important at lower velocities.
Campillo et al [64] and Pitarke and Campillo [65] have derived an expression for a position dependent stopping power within the linear response approximation using a response function calculated with time-dependent DFT (see section 5.3). Comparison of stopping powers calculated for aluminium with equivalent calculations for a free-electron gas at the same density show an enhancement in the stopping power of 7% in the former case at low velocities (v < v0). This enhancement can be attributed to band-structure effects and, given that aluminium is free-electron-like, can be interpreted as a lower bound for the importance of such effects in metals generally.
Campillo et al [64] calculated the stopping power of protons at v = 0.2v0 along (100) and (1 1 1) channels in aluminium as a function of impact parameter with the channel walls. They found variations of up to 20% around the average stopping power.
Another model that yields an environment-dependent stopping power is that of Dorado and Flores [66], who used a linear combination of atomic orbitals (LCAO) [67] model to study helium projectiles in alkali metals. Calculations of the stopping power of He in a (10 0) direction in sodium as a function of the position of the trajectory in the plane perpendicular to the channel show variations of up to 100% around the average value. This variation is much stronger than that found by Campillo et al. We should note that no direct experimental validation is possible in either case.
2.4. Empirical models of electronic stopping power
If our primary aim in modelling electronic stopping is to predict stopping powers for arbitrary projectile energies and arbitrary combinations of projectile and target, then the most direct approach is to interpolate and extrapolate from the range of available experimental data. In following this route we lose much opportunity for physical insight, but as a practical tool for providing inputs to other calculations the approach has merit and is still much used (for example in the molecular dynamics simulations considered in sections 4 and 5).
The potential success of empirical fitting methods is indicated by useful scaling relations in the behaviour of experimental stopping power data. Figure 8 illustrates how data sets for different targets and projectiles can be superimposed via some simple transformations, leaving the
B in Au Be in Au Li in Au He in Au H in Au H in Fe H in Al
1—1-1—H-
100 1000 Energy (keV)
10000 100000
: + ; X i 1 i i i | i i i | i B in Au Be in Au i i I i < < i < ' i: (b) :
. □ Li in Au
o He in Au
: A H in Au
; v H in Fe
■ « H in Al djaJMfa,
« ÄgpP......
: oi^A ■ A A ........... V ..........
Energy/Z
100 4/3
1000 (keV/amu)
Figure 8. A demonstration of how simple scaling relationships (b) can capture much of the behaviour of the electronic stopping power (a) for a variety of projectile and target combinations. The scaling of stopping power by 1 /Z2 is informed by the prefactor in the fast particle stopping theories (6), (8) and (18) and the normalization of the particle velocity by 1 /Z^3 is suggested by the Thomas-Fermi scaling of electronic velocities. (Data are from the database of Paul [25].)
fitting process to capture the underlying shape of the stopping power curve as a function of projectile energy and any residual structure in the data.
Two popular fitting schemes are those of Ziegler et al [5], embodied in the SRIM code10 [68] and of Paul and Schinner [69,70], implemented in the MSTAR code [25]. The approach of Ziegler et al is based on scaling proton stopping powers, S(Z1 = 1,Z2, v), by an effective charge fraction y such that
S(Z1,Z2 ,v) = Z2y2 S(Z1 = 1,Z2 ,v).
The proton stopping powers are calculated using the local density approximation of Lindhard and Scharff [71], which writes the stopping power of a target medium as an integral over the electron density of the medium, ne (x),
S(Z1 = 1, Z2, v) = j dx n(x)S(Z1 = 1, ne(x)), (40)
where S(Z\ = 1,ne) is the stopping power for a proton of a free-electron gas of constant density ne calculated using
10 SRIM stands for 'the stopping and range of ions in matter'.
dielectric stopping theory. Hartree-Fock models of the target atoms are used to determine ne(x). Empirical fitting finds its way into the model via an adjustment to the calculated proton stopping powers in the form of a fitted empirical factor varying between 1.0 and 1.2. A further empirical fitting function is then used to determine the state of ionization of the projectile and this provides the basis for a calculation of the effective charge fraction y from theoretical considerations.
Paul and Schinner [69,70] take stopping data for helium as their experimental reference and fit the quantity,
Srel —
5(Z1,Z2 ,v)/Zl S(Z1 — 2, Z2,v)/(2)2 ,
as a function of Z1 and v. A three parameter fitting function of prescribed form is used to fit data at each value of Z1. The parameters themselves are then fitted as functions of Z1 to yield a universal fitting scheme.
In [72], Paul and Schinner compare experimental data for stopping of carbon projectiles in amorphous carbon targets with the prediction of various empirical fitting and theoretical models. The fitting models of Ziegler et al and Paul and Schinner give an understandably good match, whilst the binary theory of Sigmund and Schinner [44, 45] performs similarly well over the full projectile energy range from 1keVamu-1 to 100MeVamu-1. The UCA (unitary convolution approximation, see page 31) of Grande and Schiwietz [49] performs well at higher energies, but significantly underestimates the stopping power for projectile energies below 1 MeV amu-1.
In this section we have considered a variety of models for the energy transfer from ions to electrons when the ions are sufficiently energetic to be moving ballistically through a stopping medium. In this regime, the quantity of interest is the electronic stopping power and it is the aim of most models to provide reliable predictions thereof.
We have seen that the complexity of the physical problem to be solved varies considerably over the range of possible projectile energies, with many more factors to be taken into account at lower energy. The simpler problem of the stopping of fast particles (v > v0) of low charge (Z1 < 2) was tackled first and with most success by the research community and continues to attract the bulk of theoretical and experimental attention. In contrast, models of the stopping of slow ions are less well established. This is, at least in part, justified by the relative dominance of nuclear stopping at lower projectile velocities, but even in this regime we still expect a net energy transfer from ions to electrons in most situations of interest. We should note that the majority of practical applications shown in figure 1 lie squarely within the low velocity regime.
At very low velocities (v < 0.1 v0) the stopping power formalism may no longer provide an appropriate treatment in many cases. In radiation damage collision cascades in particular, a separation of electronic and nuclear losses into two force-like terms looks questionable and a direct simulation of the energy exchange processes, taking into account both ion and electron dynamics explicitly, may be called for. Some approaches of this type will be considered in section 5.
3. Models of electron-phonon coupling
The electronic stopping power theories, which model the interaction of ions and electrons via a force on the ions opposing their velocities, are derived assuming that each ion's motion can be treated independently. Either a single ion moving through some continuous stopping medium is considered or the ionic motion is regarded as a sequence of binary collision interactions, occurring independently of the surrounding ions. These approximations make sense when ions move with high velocities, interacting only fleetingly with surrounding ions except when a collision takes place and their interactions are then dominated by the strong binary interaction forces at close approach. Later on in a collision cascade, however, the initial PKA energy will be shared amongst many ions, all moving in complicated motion relative to one another. At some stage, the cascade may resemble a molten region of the target material. The forces on a given ion will not be dominated by those from any single neighbour and a many-atom treatment will be necessary. At still later times, the lattice of the target material will have largely healed and we must consider the collective motion of a thermally excited set of ions, representable as a superposition of phonons.
Radiation damage theory has traditionally treated the exchange of energy between ions and electrons in these later, lower energy stages of a cascade as a separate problem from that of electronic stopping. This electron-phonon coupling regime has given rise to its own body of theory, working with a different set of approximations, but it is important to realize that the physics of all the energy exchange processes is fully captured by the time-dependent Schrodinger equation (TDSE) for a set of interacting quantum mechanical nuclei and electrons. Indeed, as we will discuss in section 5, appropriate approximations to the TDSE can, in principle, give rise in a single simulation framework to all the phenomena treated by electronic stopping and electron-phonon coupling theories. However, the practical and theoretical challenges in implementing such a framework on the time and length scales of typical radiation damage phenomena are huge. Information and insight provided by theoretical models of the electron-phonon coupling regime thus remain important for informing and validating the various methodologies used to capture energy exchange processes in atomistic simulations at lower ion energies.
Our account of the theory of electron-phonon coupling will be briefer than that of the theory of electronic stopping given in section 2 for three reasons. First, whilst the study of electron-phonon coupling is an entire research field in itself, it is of much broader relevance than simply to radiation damage phenomena. Much of the literature is concerned with the concept of the electron-phonon coupling in general (see, for example, [73] or [74] for excellent accounts). Here we will focus on the work that directly addresses or informs atomistic simulations of radiation damage.
Second, as we shall see, experimental and theoretical estimates of numerical measures of electron-phonon coupling, generally calculated as a rate of energy transfer between the electronic and ionic subsystems per unit volume per
degree temperature difference, can vary over several orders of magnitude in a given material. This stands in stark contrast to the high precision with which electronic stopping powers can often be measured and reduces the direct usefulness of the data for parametrizing atomistic models. As discussed in sections 4.4 and 5.1, the large uncertainties often prompt researchers to explore a broad range of values for the electron-phonon coupling in their simulations, rather than relying on any one particular estimate from the literature.
Third, attempts to incorporate the effects of electron-phonon coupling into atomistic simulations are still at an early stage and even somewhat controversial. Thus far, only very simple, single-parameter models have been tried (sections 4 and 5.1) and it is far from clear how the true energy exchange in the events modelled in such simulations relates to the circumstances treated in experimental and theoretical studies of electron-phonon coupling. As a specific example, several simulation schemes incorporate electron-phonon coupling as a simple drag force (similar to a stopping power) assumed to act on lower energy ions. This force acts on ions individually, completely ignoring the collective nature of the motion that lies at the heart of theoretical treatments of electron-phonon coupling. Such drag forces are also assumed to apply up to cascade energies at which the target material may be molten and a description in terms of phonons may be inappropriate. When the term 'electron-phonon coupling' is applied in radiation damage theory it should be regarded as referring to the rate of exchange of energy between electrons and ions in general below some energy threshold, rather than specifically to the interaction of electrons and a well-defined phonon system.
For the reasons given above, we will focus our discussion on the few theoretical treatments of electron-phonon coupling that have been heavily cited in the radiation damage literature and have been important in the development of atomistic radiation damage models.
3.1. The importance of electron-phonon coupling in radiation damage
A short time after a thermal spike has formed, when almost all atoms over an extended region of size ~ 1000 A are significantly excited, but few atoms are moving ballistically through the target medium, we can regard our system as being composed of an excited ionic subsystem interacting with an electronic subsystem. The electrons will likely be initially much cooler than the ions and so will in the first instance act as a heat sink. Because of the relatively high electronic thermal conductivity, any energy transferred from the ions to the electrons can be rapidly transported away and the electronic subsystem provides an important mechanism for cooling the ions. Depending on the energies involved, and on the balance between the rate of energy exchange between ions and electrons (i.e. the strength of the electron-phonon coupling) and the electronic thermal conductivity, this cooling might be rapid enough to inhibit the production of defects early in the evolution of the displacement spike or it may act to quench in defects as the system returns to equilibrium.
This balance between the rates of energy exchange and energy transport was investigated by Flynn and Averback [75]. They considered a thermal spike formed by depositing an energy Q into a spherical region of the target of radius r. The energy per ion is thus Q(r0/r)3, where r0 is the Wigner-Seitz radius11. Writing the energy per atom in terms of an ionic temperature TI, i.e. as 3kB TI, we can consider the evolution of the thermal spike as it grows and cools to be described by
r(t) a
3kB Ti (t)
Energy exchange between ions and electrons will occur when electrons are scattered from state to state by imperfections in the crystal lattice, emitting or absorbing phonons. A measure of the rate of this scattering is the distance travelled by an electron between scattering events, the electron mean free path, Xmfp. In the harmonic approximation, when the ions are not much displaced from their equilibrium lattice sites, this can be written
rpTp Ti(ty
where A.^ <
where T0 is the ionic temperature at which Xmfp — r0. T0 thus encapsulates information about the strength of the electron-phonon coupling, with a low value corresponding to a strong coupling (a high rate of energy exchange). Values for T0 can be calculated from the electrical resistivity using formula [76] Amfp — (92 X 10-18 Qm2)(r0/pe(Tref)flo), where pe(Tref) is the electrical resistivity at a reference temperature Tref and a0 is the Bohr radius. Comparing the mean free path to the size of the thermal spike
Amfp _ / 3kBY/3 _T
Q ) T{/3(t)
Flynn and Averback point out that because Xmfp has the stronger dependence on ionic temperature, there will always be some point in time at which the thermal spike has cooled sufficiently that Xmfp > r .At this point, heating of the electrons by the lattice will be very inefficient. The critical factor then is T0. For a material with strong electron-phonon coupling and so a low value of To, heating of the electrons by the lattice will tend to remain effective for longer times and the electrons will tend to remain in equilibrium with that lattice throughout the development of the thermal spike.
To quantify the likely variation in the behaviour of real metals, Flynn and Averback present a picture of electrons diffusing out of the thermal spike, acquiring an energy kB©D with each scattering event, where ©D is the Debye temperature. If escape involves a random walk with on average (r/Xmfp)2 scattering events then the electrons will acquire a temperature Te — ©D(r/Xmfp)2. Equilibration of the electrons with the ions, Te — TI, will therefore be possible whilst the thermal spike remains above a critical temperature,
11 Defined as the radius of the spherical volume equivalent to the volume per atom in the solid, i.e. 3nr3 — 1/na for a number density of atoms na.
which has a very strong dependence on T0. Flynn and Averback [75] quote values for T0 of 4.5 x 104K and 1.5 x 104K for copper and nickel, respectively. These imply corresponding values for Tcrit of 2 x 105 and 300 K, suggesting very different behaviours for the two metals.
3.2. Two-temperature models
A particularly simple picture of the interacting electronic and ionic subsystems emerges if we assume that though they are out of equilibrium with one another, they are each internally in an equilibrium state, so that the state of the combined system can be parametrized by an electronic and an ionic temperature, Te and TI, respectively. Such two-temperature models have been widely explored and form part of some of the atomistic simulation schemes [8,9,77-80] to be discussed in section 5.
Assuming an initially hot ionic system, the two-temperature model will be valid provided the rate of thermalization of the electronic subsystem is significantly higher than the rate of energy transfer into the subsystem. Since electron-electron interaction time scales are typically of the order of a few hundred femtoseconds12, whereas electron-phonon interaction time scales vary between tens and hundreds of femtoseconds13 the approximation may be invalid in many cases. In fact, as we discuss in the context of some results from recent simulations [14] in section 5.4, the above condition can be considerably relaxed because the mode of excitation of the electronic subsystem is such that it remains close to equilibrium with a steadily rising temperature.
We can represent the evolution of a two-temperature model using coupled heat flow equations for the evolving spatial temperature distributions Te(x, t) and TI(x, t), d
Ce(Te) — Te(x,t) = Vxke (Te, TI)VxTe]- gp(Te, Ti)[Te -Ti], dt
ci(Ti) — Ti(x, t) = Vx[ki(Ti)VxTi] + gp(Te, Ti)[Te - Ti], dt
where ce and cI are the electronic and ionic heat capacities per unit volume and Ke and ki are the electronic and ionic thermal conductivities, respectively. gp(Te, TI) is the electron-phonon coupling (measured in Wm-3 K-1 or dimensionally equivalent units).
By extending the arguments of Flynn and Averback [75], Finnis et al [83] develop a scheme to estimate the value of gp in different metals. This scheme will be discussed in section 3.4. Their numerical solutions of (46) and (47) reveal the much more rapid cooling of the ionic subsystem in nickel when compared with copper (see figure 9) and suggest an important role for electron-phonon coupling in determining collision cascade dynamics.
12 For example, Del Fatti et al [81], in femtosecond laser experiments, find time scales of 350 fs in silver and 500 fs in gold.
13 The time scale of the electron-phonon interaction can be calculated as the ratio of the electronic heat capacity per unit volume and the electron-phonon coupling constant (defined below). Using experimental data from Qiu and Tien [82] gives values of 650 fs for copper (weakly coupled) and 64 fs for vanadium (strongly coupled).
1 1 .....1 ' 1 1 ' 1 1 ■ 1 - cu
- \ \ \
- \ V 1 1 . 1
100 200 lime (is!
400 1000 2000 4000 10.000
Figure 9. The evolution of the ionic temperature given by the coupled equations (46) and (47) for values of gp, the electron-phonon coupling constant, for Cu and Ni. (Reprinted with permission from Finnis et al [83]. Copyright 1991 American Physical Society.)
3.3. The electron-phonon drag
Finnis et al [83] consider how the two-temperature picture could be incorporated into an atomistic simulation. Noting that the ionic thermal conductivity will generally be small, they simplify (47) to
dTI gp
-7- = — (Te - TI).
This equation then gives the rate at which energy should be removed from (or injected into) the ions in an atomistic simulation due to their interaction with the electronic subsystem. Ideally, the required energy change of the ionic system would be effected by the excitation or de-excitation of the physically correct phonon modes, but such a process is computationally expensive. Instead, approximate methods such as uniformly scaling (up or down) the ionic velocities or applying a force parallel to the velocity are generally used. Energy transfer from electrons to ions is also often implemented using a stochastic force. Finnis et al adopt the use of a damping force, defined for the i th ion, with velocity vi as
Fi = -Pi vi. (49)
This force does work on the ion at a rate of -Pi v2, which can be equated with the rate of energy transfer due to the electron-phonon coupling 3kBdTi/dt = (gp/cI)(Te - TI), where we have introduced a temperature per ion, Ti. If we further equate the thermal energy associated with Ti with the ionic kinetic energy, 3kBTi = miv.2 for ions of mass mi, then we can write the drag coefficient Pi corresponding to the electron-phonon coupling as
gpmi Ti Te
Pl = Q (. (50)
In the case where the ions are hotter than the electrons Pi will be positive and will provide a drag force that removes energy
from the ions. In the opposite case of Te > Ti, a negative pi will provide an accelerating force. In a practical simulation scheme, Ti can be replaced by an average TI over some coarsegrained cell of the ionic system [8]. Finnis et al [83] point out that if one is calculating the 'ionic temperature' Ti of an ion it is necessary to deal with the singularity at Ti = mi v2/3kB = 0 for momentarily stationary ions. They do this by making the substitution T-1 ^ IT2 + (Te/20)2}-1/2, where the factor of 1/20 is chosen to be compatible with their simulation time-steps.
3.4. Models of electron-phonon coupling
Various analytical models for calculating the electron-phonon coupling gp have been proposed. One of the most physically transparent is due to Finnis et al [83] and arises from an extension of the analysis by Flynn and Averback [75]. We will consider their derivation and compare the result with other commonly cited analyses.
Finnis et al [83] once again consider an electron with mean free path Amfp = r0 T0 /TI acquiring energy kB©D in each of a series of collisions with lattice distortions. If the local electronic temperature is Te and the electronic density of states at the Fermi level is D(eF) then ~kBTeD(eF) electrons will be able to participate in the energy exchange, given the requirements of quantum mechanical exclusion. The scattering rate will be vF/Amfp, where vF is the Fermi velocity, and so the rate at which electrons acquire energy will be
dEe kB ©o D(sF)vF Ti Te
kB ©oD(sF )vFTe
corresponding to a damping coefficient,
3©pCeVFmi (Ti - Te n 2roToci
been used on a larger scale by other researchers [8,78,84-86] and we will discuss them in detail in sections 4 and 5.
Various more formal treatments of the electron-phonon coupling exist in the literature [17, 87, 88], but they all reduce to a form similar to (52), with different numerical prefactors depending on the extent to which they contain details of the true electronic and lattice structures. These treatments all begin by considering the general case of a quantum mechanical ionic subsystem characterized by phonons of energy h (q) and momentum hq populated according to a set of occupation numbers N(q, s), where s indexes the phonon branch. This system is coupled to a quantum mechanical electronic subsystem of electronic states of energy ev (k) and momentum hk with occupations f(ev(k)), where v is a bandindex. At equilibrium the f(e(k, v)) will have a Fermi-Dirac distribution and the N(q, s) a Bose-Einstein distribution. If we consider the case of an electron scattering from state (k, v) to state (k', v') with the emission or absorption of a phonon (q,s), then we must ensure energy conservation ev> (k') -ev(k) = hQ.s(q) and momentum conservation k' — k = q. A Fermi's golden rule (FGR) analysis [32] (first-order time-dependent perturbation theory) treating the lattice distortion due to the phonons as a perturbation gives an expression for the rate of the transfer of energy hm and momentum q from the ions to the electrons [17],
s k,v kk
xS(sv(k') - Sv(k) - hQs(q))\Vkvk>v>(q,s)\2
x f(ev(k)[1 - f(ev' (k' )] N (q, s)
Finnis et al then turn this into a net rate of energy transfer by letting TI ^ (TI - Te). Whilst this gives an energy transfer that vanishes as required when TI = Te, there is no physical argument to support the substitution. TI appears in (51) as the temperature determining the rate of scattering of electrons by ions and we would not expect this scattering rate to go to zero at equilibrium between electrons and ions. What we would expect to go to zero is the net average energy transfer associated with such collisions, and in (51) the temperature determining this is ©D. The above comments not withstanding, we will continue to follow the analysis given by Finnis et al and write for the electron-phonon coupling,
- f(sv,(k')[1 - f(ev(k)][N(q, s) + 1] ,
where Vkvk>v> (q, s) is the coupling matrix element between electronic states (k, v) and (k', v') due to the lattice distortion by phonon mode (q, s). The first term in braces corresponds to the stimulated absorption of a phonon and the second term to stimulated and spontaneous emission of a phonon. The rate of energy absorption by the electrons can then be written as
dEe dt
= y drnhrn).
The task now is to find suitable approximations to formulae (54) and (55). Kaganov et al [88] approximate the electronic system as a free-electron gas and find in the case of a hot electron-ion system not too far out of equilibrium, Ti > ©d, - Te\ ^ Ti, an electron-phonon coupling constant,
nme UeV,
6Te Te
where, again, ce and cI are the electronic and ionic heat capacities per unit volume and the density of states D(eF) has been subsumed into ce = (n2/3)kBD(eF)Te.
Finnis et al [83] undertake some preliminary simulations in copper and nickel using their damping scheme. Their results are inconclusive about the effect of electron-phonon coupling on the residual defect population, but similar schemes have
where Te — Amfp/vF is the electron-phonon scattering time and vs is the speed of sound in the lattice. Equation (56) has the same form as the result (52) due to Finnis et al [83] and the two formulae are made equivalent by the transformation,
2 X1/3 T2
Ti ©d'
Table 1. A sample of experimental and theoretical estimates of the electron-phonon coupling gp from the literature. All data given in units of 1016 Wm-3 K-1. See main text for a discussion of the trends and variation in the values.
Theoretical Experimental
Equation (52) Equation (56)
Finnis Gao et al Wang et al Qiu and Qiu and
et al [83] [84] [89] Tien [82] Tien [82]
V 4803a 183 648 523 ± 37
Cr 179 45 42 ± 5
Fe 1815.0 119
Ni 3164.1 1714.5 107
Cu 81.9 40.1 (36.4a) 12.7 14 4.8 ± 0.7
Ag 9.4a 3.34 3.1 2.8
Au 14.2a 2.3 2.6 2.8 ± 0.5
W 27.6 27 26 ± 3
a Calculated by the present authors using the approach of Gao et al
where z is the number of valence electrons per atom in the free-electron gas of Kaganov et al [88].
An analysis by Koponen [17] yields an expression for the energy transfer rate
^ = 4nD(£F) i doja2F(o)hoj)2{N(Tl) - N(Te)}, dt J
where N(T) is a Bose-Einstein distribution at temperature T for the phonon occupations. The derivation of (58) is too involved to admit a simple discussion, but the important point is that the so-called spectral function a2F(co) is an experimentally measurable property of a metal, encapsulating the information about the electronic and phonon densities of states and the coupling matrix explicitly present in (54).
3.5. Estimates of electron-phonon coupling
As a demonstration of the uncertainty involved in selecting a numerical value of the electron-phonon coupling for use in an atomistic simulation, we give a selection of literature estimates of gp in table 1. Included are values calculated by Finnis etal [83] using (52) and values calculated using the same formula by Gao et al [84]. These differ by approximately a factor of 2 due to the choice of values for T0.
Also included are theoretical values calculated by Wang et al [89] and Qiu and Tien [82] using formula (56) due to Kaganov et al [88]. The large differences in the values for chromium and vanadium are due to differing assumptions about the number of valence electrons contributing to the free-electron gas density: Wang et al assume ne = na; Qiu and Tien assume a variable ratio 0.5 < ne/na < 2.0.
The large difference between the values for nickel calculated using (52) and (56) is due to the presence of the electronic density of states at the Fermi level in the former. This is particularly high in nickel (compared, say, with copper) and thus implies a large electron-phonon coupling. This band-structure dependent effect is absent from the free-electron based formula (56).
Experimental values of gp were tabulated by Qiu and Tien [82], derived from short pulse laser heating experiments
(see [90,91] for examples, [92] for a review) in which excess energy is injected into the electronic system by a laser pulse of ~100 fs duration and the system monitored as it relaxes to equilibrium. These values show much better agreement with calculations using (56) than with those using (52), despite the presence of more band-structure dependent effects in the latter formula.
Overall, the values of electron-phonon coupling obtained by each method show a variation of several orders of magnitude between gold (weakly coupled) and vanadium (strongly coupled). Such variation implies the possibility of strong material dependent effects on the later stages of cascade development.
The relative variation in the experimental values between copper and vanadium is best reproduced by formula (52) as implemented by Finnis et al [83] and Gao et al [84]. Formula (56) used by Wang et al [89] is only able to capture the variation if the number of valence electrons is varied. Qiu and Tien [82] use such a variation, but the justification for, for example, setting ne = 2na for vanadium, but ne = 0.5na for its neighbour chromium is not clear.
Equally importantly, for any given material, the estimated values of gp vary by over an order of magnitude and it is not clear which value should be adopted as a parameter in any particular simulation scheme.
4. Implicit incorporation of electronic effects in atomistic simulations
In the following two sections we will examine the various methods which exist for incorporating the effect of the electronic system into simulations of radiation damage. Such methods fall into two broad classes depending on how the electrons are modelled and we will treat each in turn. To date, the bulk of effort has focussed on methods from the first of these classes—those which treat the electrons implicitly within a simulation of ion behaviour. A second class of methods which incorporate an explicit model of the electrons has recently begun to be explored and offers much opportunity for future work. These methods are considered in detail in section 5.
4.1. Experimental benchmarks
The main aim of most simulations of radiation damage is to establish as accurately as possible the primary state of damage at the end of the relaxation phase and the factors that influence it. In considering the validity of the various approaches we will need to bear in mind the extent to which experimental verification is possible. We will not consider the methods used to obtain such experimental data (see Averback and Diaz de la Rubia [15] for a review), rather only the types of information available. These are the following.
Range distributions. The perpendicular distance from the surface at which an intruding ion comes to rest is known as its range and the standard deviation of the range is referred to as the straggling. Large quantities of experimental data on both quantities are available, but as Ziegler et al [5] point out different experimental techniques involve different biasing factors and often use slightly different definitions of range and straggling. The data available for model fitting thus exhibit broad scatter in many cases.
Mixing. The ionic mixing parameter is defined as
<x2> = J2[Ri(t) - Ri(0)]2
where {Ri(t)} are the ionic positions at some time t during a cascade, {Ri(0)} are the corresponding initial positions, na is the ionic number density and Tc is known as the cascade energy and is defined as that portion of the initial PKA energy available to cause atomic displacements (i.e. it excludes energy lost to the electrons). fIM is normally given in units of A5 eV-1 and uncertainty in the experimental measurements is around 10-20% [93].
Displacement threshold energy. The displacement threshold energy (or just threshold energy), Ed, is defined as the energy required to make a stable Frenkel pair. In crystals this energy will vary with direction, typically being significantly lower along close-packed directions. For a few materials the direction-dependent threshold energy surface has been mapped experimentally (see King et al [94] for an example in copper). Average threshold energies, measured using polycrystalline samples, are available for a broad range of materials [15].
The threshold energies, and particularly the threshold energy surfaces when available, are used for calibrating the short range repulsive portion of interatomic potentials in MD. After such a fitting to the data for copper in King et al [94], Foreman et al [95] found a good match to the shape of the experimental threshold energy surface from their simulations. King and Benedek [96] obtain good qualitative agreement with an unfitted Born-Mayer potential.
Damage function. The damage function v(T) is the number of defects produced by a PKA of energy T. Damage functions are often expressed in terms of a damage efficiency, fFP, defined via
v(T) = fo (T)vkp (T), (60)
where vKP is the Kinchin-Pease damage function (see (61) and discussion in section 4.2). Wallner et al [97] give experimental data for which they claim a maximum error of 15%, though this excludes further possible errors involved in their calculations; they rely, for instance, on data for the incremental resistivity due to an additional Frenkel pair, which are known only to within ±25% and are not expected to be additive at higher damage concentrations.
Clustering and defect distributions. Various methods exist for mapping out the distribution of defects in experimentally irradiated samples. Though all have their limitations (they might be sensitive to the effects of the surface of a small sample for instance) they can be used to build up a picture of typical defect distributions in different materials. The distributions and data on clustering provide support for the ability of classical MD to reproduce the correct broad patterns of behaviour.
Lengths of replacement collision sequences. The length of a typical replacement collision sequence (RCS) in a cascade simulation will be sensitive to the details of the interatomic potential and so experimental measurements of RCS length can be useful for calibration. Such measurements are not possible directly, but RCS lengths can be inferred by various means. As an example, Wei and Seidman [98] give a distribution of RCS length resulting from 20 keV self-irradiation of tungsten.
4.2. Damage functions
Before the availability of general purpose computers flexible enough to perform simulations of radiation damage events, predictions of damage relied on analytical approaches. We can view many of these approaches as cascade simulations in thought-experiment form, allowing expressions for the damage function, and range and straggling distributions to be derived.
The most well known of these expressions is the so-called Kinchin-Pease damage function, vKP (T). Kinchin and Pease [99] considered the case of a PKA of initial kinetic energy T initiating a series of collisions within a solid. The cascade is imagined to develop via a sequence of rounds of collisions, with each round doubling the number of ions involved in the cascade and the initial energy being shared via hard sphere collisions. The probability distribution for the energy transfer in such collisions is uniform between zero and the maximum possible for the projectile to target mass ratio: a simple form permitting an analytical solution for the energy distribution of all the atoms in the cascade in a given round of collisions. In a given round, any ions with energy between one and two times the displacement threshold energy, Ed, will be displaced themselves, but will not go on to displace further atoms. Thus the total number of atoms displaced in the cascade can be calculated as a sum over the energy distribution in all collision rounds. The final expression is
vkp (t) =
Various extensions and refinements to this expression have been proposed, but the most widely adopted is the Norgett-Robinson-Torrens (NRT) expression, also sometimes referred to as the modified Kinchin-Pease damage function. An obvious deficiency of (61) is the lack of any account of electronic energy loss. In fact, in some applications vKP is allowed to saturate above some energy 7^rit to take account of the fact that at high energies the energy loss from an ion is predominantly to electrons, but this modification is crude. Instead, Norgett et al [100] propose a revised form of the damage function
VNRT (T) = —, (62)
where T is that portion of the PKA kinetic energy not lost to electronic excitations (calculated using the universal functions of the LSS theory, see Lindhard et al [101]). The factor k is the displacement efficiency, designed to take account of the possibility of recombination of some of the Frenkel pairs formed in the cascade. Norgett et al propose a universal value of k = 0.8 based on a series of simulations carried out by Robinson and Torrens [102].
4.3. Binary collision models
Expressions like the Kinchin-Pease (61) and NRT (62) damage functions result from considering radiation damage events as a series of isolated binary collisions. An alternative to addressing such a process analytically is to carry out large simulations within the same approximation and directly observe the damage produced. This approach is called the binary collision approximation (BCA).
Simulations in the BCA begin with a PKA with some initial kinetic energy T0 moving through a simulation cell containing other ions. This PKA moves in a straight line at constant speed until it comes within range of another ion, at which point it is considered to undergo a collision with that ion. The velocities of the projectile and the target ion after the collision are calculated using the theory of simple scattering under an assumed potential. If the kinetic energy transferred from projectile to target is Tt and the projectile is left with kinetic energy Tp then there are four possible outcomes from the collision:
(i) Tt > Ed, Tp > Ecut : the target atom joins the cascade and both atoms go on to undergo further collisions.
(ii) T > Ed, Tp < Ecut: the target atom joins the cascade and the projectile replaces the target at its lattice site.
(iii) T < Ed, Tp > Ecut: the target atom remains on its lattice site and the projectile proceeds on a modified trajectory to undergo further collisions
(iv) T < Ed, Tp < Ecut: the target atom remains on its lattice site and the projectile becomes an interstitial atom.
The cut-off energy Ecut is a simulation parameter chosen to improve the results and need not take the same value as the displacement threshold energy Ed.
Various different features have been incorporated into the basic BCA scheme since the earliest simulations by Yoshida [103]. A clear distinction exists between models of amorphous
targets (see the early work of Oen et al [104,105]) and more computationally demanding treatments of crystalline lattices. In the former case, target ions are generated for collisions with cascade ions based on a random spread around some mean free-flight path. In the latter case, a search for collision partners must be made within the lattice surrounding the cascade ion path. Early simulations using the BCA in ordered lattices [106-108] provided the first means of exploring the anomalies discovered in experimental range distributions in polycrystalline targets [109] and attributed to ion channelling.
An important feature of BCA simulations is their relatively low computational cost. Because they treat ionic collisions in a very simple way and because the description of the target material is effectively generated 'on the fly' in the vicinity of cascade ion paths, it is possible to simulate large numbers of cascade events up to very high PKA energies and to generate good statistics for the damage distribution. This feature of the BCA has ensured its continued use alongside more realistic, but much more computationally costly, MD methods, even up to the present day. Well developed codes have been widely used, including the SRIM code [68] for simulations of amorphous targets and Marlowe [102] for the treatment of crystalline targets. Results from the latter were used to determine the displacement efficiency, k in the NRT expression (62).
However, the BCA has some major flaws. Whilst simulations can provide defect distributions if the initial and final positions of the cascade ions are recorded [110,111], there is no inbuilt mechanism for recombination of interstitials and vacancies to take place. The simplest way to allow for such healing is to define a recombination radius within which Frenkel pairs will be eliminated. Beeler [112] has also explored the possibility of defect annealing by allowing for a phase of diffusive defect motion following the conclusion of the cascade. More generally, we cannot be confident that the cascade evolution, and therefore the damage distribution, is realistic, because of the approximation of treating only binary interactions. This approximation gets worse at the lower ion velocities at which the final damage state is formed. This problem is highlighted by differences in the BCA literature in the way that the cascade evolution is followed. Which collision should be considered next at any particular time? When should multiple collisions be deemed to occur simultaneously? And how should cascade evolution be allowed to interact with existing defects given that the 'true' chronology of the cascade is uncertain?
From the point of view of this review, a key issue is the incorporation of electronic effects within the BCA. Many of the earliest simulations made no attempt to take into account the loss of ionic kinetic energy to the electrons. Later schemes such as those in MARLOWE and SRIM incorporate an energy loss into the calculation of the projectile and target velocities following a collision. The collisions are still treated elastically within the centre of mass frame, but an energy loss is calculated using the models of Firsov [52] or Lindhard and Scharff [53] (see section 2.3.1), which then influences the final trajectories.
It is not immediately obvious that such a means of incorporating electronic effects will give the correct cascade behaviour. In particular, when the predominant mode
of evolution involves glancing collisions (in the case of channelling) or very efficient energy transfer from projectile to target (in the case of an RCS) we might expect a more sophisticated treatment to be necessary. Whilst the outcome of BCA simulations agrees well with experimental results for ion range and straggling distributions, such statistics afford only a crude view of cascade behaviour; it is not clear to what extent such agreement is the result of the fitting of the values of the many empirical parameters within the simulation schemes.
4.4. Molecular dynamics with drag
An obvious improvement over the BCA would be to include in the cascade simulation an explicit representation of the ions and the forces between them. Classical molecular dynamics (MD) implements just such an improvement and, contingent on a good choice for the interionic force model, allows for much more realistic simulations of radiation damage and for the study of more complex phenomena. Computational power has been sufficient for theMD simulation of collision cascades since the 1960s (see Gibson etal [113] for an early example and various reviews [114,115]). In the early simulations, electronic effects were ignored.
Since the 1990s various models have explored the use of a viscous drag force, opposed to the ion velocity and proportional to its magnitude, to represent the effects of energy transfer to the electrons in both the electronic stopping and the electron-phonon coupling regimes. Such a force is clearly consistent with the slow particle stopping models discussed in sections 2.3.1 and 2.3.2 for the electronic stopping power and also conforms to the model proposed by Finnis et al [83] for the electron-phonon coupling (see section 3.3). We are thus considering ions moving under an equation of motion,
MiR = Fi - PiRi, (63)
where Fi is the force on ion i of mass Mi due to the other ions under the chosen interatomic potential. Ri, R and R indicate the ion position and its time derivatives and pi is the drag coefficient. In the simplest models pi is chosen to be a constant for all ions.
Nordlund et al have made extensive studies [93,116118] of the effect of electronic stopping power on the primary damage state. Their MD simulations use embedded atom model (EAM) potentials [119], adjusted to give correct melting temperatures and joined at distances of close approach to a repulsive potential (chosen according to the specification of Ziegler et al [5]) to improve the handling of energetic collisions. The electronic stopping power is incorporated as a frictional force, proportional to ion speed with constants of proportionality drawn from the SRIM code [5] and assumed to act only on ions whose kinetic energy Ti exceeds 10 eV:
Pi = P Ti > 10 keV,
i i (64)
= 0 Ti < 10 keV.
No clear reason for the kinetic energy cut-off is offered, except that the Lindhard form for the stopping power is only considered valid down to this velocity range.
In addition to this, a kinetic energy cut-off in the damping force will ensure that the long-term steady state of the ionic system has a finite temperature rather than all the ions being damped to zero velocity. However, the approach to this steady state will be unphysical. Rapid repartitioning of the initial PKA energy TpKA amongst all the Na ions will probably take the average kinetic energy per ion, ~TpKA/2Na (assuming equipartition between kinetic and potential energy), below the threshold within the time scale of the simulation. Statistical fluctuations will then repeatedly take some ions above the threshold until the ultimate steady state, in which the average total kinetic energy is 5 eV is reached. The temperature of this steady state, TI = 2(5eV)/3kBNa, though it would take an extremely long time to reach in practice, is thus system size dependent.
A common theme in the publications of Nordlund et al is a conclusion that the electron-phonon coupling, not taken into account in their model, has only a minor effect on the primary damage state. In [116] and [93] a comparison of the mixing parameter from simulations of ion beam mixing with data from experiments in several metals and semiconductors is carried out. This shows that much of the difference in behaviour between pairs of metals expected to have very different electron-phonon coupling strengths (e.g. Cu and Ni, and Au and Pt) is accounted for by the MD model without the effect of electron-phonon coupling. Nordlund et al conclude that the interionic potential (and hence the melting point and elastic properties) has the dominant effect on the extent of mixing and its variation between materials, and therefore electron-phonon coupling strength must play only a minor role. A detailed MD simulation study of the effect on the primary damage state of material properties and choice of potential is to be found in [93]. Similar conclusions are drawn by Zhong et al [118], who compare defect yields in the self-bombardment of tungsten between simulation and experiment, claiming agreement within ~30% for a model omitting electron-phonon coupling effects.
However, even if an MD model without electron-phonon coupling were to replicate the experimental results perfectly, this would not necessarily indicate the unimportance of electron-phonon coupling. Such models tend to contain various parameters, whose calibration is ultimately based on experimental data. When comparison of simulations is made only with high-level characteristics of real cascades, it is conceivable that the fitted values for these parameters might compensate for the lack of electron-phonon coupling to yield broadly correct behaviour. The details of the cascade evolution might nevertheless be inaccurately captured.
Before we discuss the above results further, we will consider the complementary work of Bacon et al [84,85]. In [84] an MD model of a -iron using a Finnis-Sinclair potential [120] is augmented with a frictional force to represent the effect of electron-phonon coupling. The simulations focus on the evolution of the cascade once a molten region has formed and the effect of the coupling strength on the defect yield is analysed. A range of values of coupling strength are tested, varying between zero and a value ten times that found for Ni in [87], and including the values for Cu and Ni calculated by
Figure 10. The evolution of the radius of the molten zone in cascade simulations in a-Fe. Increasing the electron-phonon coupling strength dramatically increases the rate of cooling of the cascade. (Reprinted with permission from Gao et al [84]. Copyright 1998 IOP Publishing.)
when the coupling strength is varied from zero to the value associated with nickel.
As pointed out by Nordlund et al [93], there is still considerable uncertainty as to how electron-phonon coupling should be treated in MD simulations of radiation damage. If we add to this the scarcity of experimental data against which to validate coupling models and the uncertainties in those data, then there is clearly a need for more work in this area. More realistic treatments of the interaction between electrons and ions could help to resolve this issue (see section 5).
4.5. Electrons as a heat bath
In addition to providing a mechanism for ionic energy loss, we also expect that the electrons in a metal will function as a heat bath, exchanging energy with the ions and enhancing the rate of energy transport away from the cascade region. Caro and Victoria [121] identify two key problems in MD simulations of radiation damage, which remain no matter how accurate the interatomic potentials become. These are the treatment of inelastic scattering and the need for a description of electronic thermal conductivity (see also [2]). In order to better capture the effect of electron-phonon coupling, they propose a treatment of the electrons as a Langevin heatbath, such that the ions obey the modified equations of motion,
MiRi = Fi + m(t) - ßiRi,
where Pi is a drag coefficient and ni(t) is a stochastic force, distributed with probability P(n)
<V) = 0, <V(t) • n(t')) = 2ßikBTeS(t - t'),
P(n) = (2tf <n2>)-1/2exp(-n2/2<n2>).
Figure 11. The number of Frenkel pairs produced in cascades with various PKA energies as a function of electron-phonon coupling strength. (Reprinted with permission from Gao et al [84]. Copyright 1998 IOP Publishing.)
Finnis etal [83]. The strength of the electron-phonon coupling is observed to have a dramatic effect on the rate at which the molten zone shrinks (see figure 10) and a less strong, though still significant effect, on the number of Frenkel pairs produced (figure 11).
The apparently contradictory conclusions concerning the importance of the electron-phonon coupling strength can be reconciled if we note that the claim by Nordlund et al of agreement to within ~30% between experiment and MD excluding coupling [118] is overstated. In the case of a 20 keV impact the simulations produced 123 vacancies compared with 81 measured via field ion microscopy (FIM). At 30 keV the data were 63 vacancies in the simulations to 125 in the experiments. As a percentage of the experimental results, these figures suggest a discrepancy of ~50%, leaving ample room for a significant contribution from electron-phonon coupling. Indeed, this is consistent with the change in the Frenkel pair yield observed by Gao et al [84] (and illustrated in figure 11)
Caro and Victoria note that a wide range of theoretical treatments yield a stopping power proportional to velocity, so that equation (65) could, in theory, describe both the electronic stopping and the electron-phonon coupling regimes, provided some means were found to account for the difference of up to several orders of magnitude in the value of ß between the two (see section 3).
The method proposed by Caro and Victoria adopts the density dependent stopping power formalism of Ziegler et al [5] (see (40) on page 42), assuming that the higher average electronic density experienced by an ion moving ballistically in the stopping power regime compared with an ion oscillating about its equilibrium position can correctly account for the variation in damping. Such a model can be efficiently implemented within MD simulations that make use of EAM potentials in which a measure of local electron density is readily available. The form of ß is determined empirically to match the linear response theory of Kitagawa and Ohtsuki [122] at high density and results derived from density functional theory by Echenique et al [123] at low density. The best fit is found to be
ßi = A log10(ap//3+ b),
where pi is the electron density experienced by an ion at Ri and A, a and b are constants whose values are given in [121]. Caro
and Victoria emphasize that there is no physical justification for this density dependence of the damping; it is simply a best fit to other models across a range of p.
Pronnecke et al [124] used the Caro and Victoria model in a set of simulations of 55 296 atoms in copper. They simulated initial cascade energies of 2.5 and 5.0 keV with electron-phonon coupling included and excluded in each case, making four simulations in total. Even the weak coupling in copper is found to reduce the duration of the cascade and decrease the extent of mixing significantly.
An acknowledged deficiency of the original model in [121] is the handling of the dependence of the stochastic force n(t) on electronic temperature. Under Langevin dynamics the fluctuation dissipation theorem states that at equilibrium n a vT and so in general some model for the evolution of the electronic temperature distribution is required. To avoid this need, Caro and Victoria assume that the rate of heat transport by the electrons and the strength of the electron-phonon coupling are such that the electronic system functions as a perfect heat sink, remaining at the target ambient temperature throughout any simulation. This is the case with the simulations in [124]. Unfortunately this assumption precludes any electronic heating, excluding the possibility of the electrons acting to anneal out defects and making it unlikely that the model will correctly describe cascades in metals with strong electron-phonon coupling.
Rectifying these deficiencies requires a change in the way that simulations deal with the electronic system, moving beyond merely attempting to capture its effect on the ions and towards an explicit model of its evolution.
5. Explicit treatment of electrons
Our consideration of the work of Caro and Victoria [121] in section 4.5 serves to highlight the limitations of attempting to treat electrons only implicitly via their effect on ion behaviour. Only relatively simple phenomena can be captured and often only with the sacrifice of physical content. Just as the explicit treatment of ion trajectories in classical MD simulations revealed a richness of behaviour uncaptured within the binary collision approximation, so we might expect that the full role played by electrons in radiation damage will be revealed only by their explicit treatment.
In this section we will introduce a variety of models that explicitly model the behaviour of electrons within simulations of radiation damage in metals. The range of work that we consider spans the state of the art in handling electrons as classical degrees of freedom within large-scale cascade simulations and recent successes in exploring the effect of quantum mechanical electrons on ion dynamics. We will also examine some further techniques that may soon be rendered applicable to radiation damage problems by advances in theory and computational resources.
5.1. Electrons as an inhomogeneous heat bath
Duffy and Rutherford [8,9] have developed an extension of the model of Caro and Victoria, which better captures the effect
of electron-ion interactions on ion dynamics by including a representation of the electrons as an inhomogeneous heat bath. The ions obey a Langevin equation of motion as in (65),
Mi Ri = Fi + Vi(t ) - Pi Ri,
but in this case the mean squared magnitude of the stochastic force ni(t) is varied throughout the simulation cell with the local electronic temperature.
The form of the damping coefficient Pi is also different from that used in [121]. The effect of electron-phonon coupling is represented by a constant Pp applied to all ions and the electronic stopping power is modelled via a second, usually much larger, constant Ps applied to ions moving faster than some threshold velocity vt,
ßi = ßp + ßs = ßp
R i > v, Ri < vt.
Values for are taken from the SRIM code [5] (fis/Mi = 1 ps-1 for bcc iron) and a variety of values for is explored (0.05 ps-1 < Pp/Mi < 30 ps-1). vt is set to correspond to an ionic kinetic energy of twice the cohesive energy.
The electronic temperature distribution is coarse-grained into cells of around 340 ions and evolved according to a heat diffusion equation,
Ce dTe = V(KeVTe) - gp (Te - Ti) + gsT/, at
where ce and Ke are the electronic heat capacity and thermal conductivity. Ti is an ion temperature defined as an average over the Ncell atoms of the coarse-graining cell as
3 1 —kBTl = -
2 B 1 NceH
and T{ is the equivalent average over only those N'cell atoms with Ri > vt. The second term on the right-hand side in (69) is thus the usual source term corresponding to energy exchange between ions and electrons via the forces -PpR-i and ni. The third term on the right-hand side in (69) is an additional source term corresponding to the electronic stopping force -Ps R. The values of gp and gs are chosen to maintain energy conservation in the exchange between ions and electrons.
Duffy and Rutherford refer to this model as an inhomogeneous Langevin thermostat. They give a full description in [8] and present preliminary results for a small number of cascade simulations. In [9] they add a further refinement by enlarging the spatial extent of the electronic system to ten times that of the ionic MD simulation, thereby enhancing the ability of the electrons to transport energy out of the cascade region. They apply a 300 K thermostat to the boundaries of the enlarged system.
Rutherford and Duffy use this enhanced model to simulate 10 keV cascades in iron, modelled using the magnetic interionic potentials of Dudarev and Derlet [125]. They explore the effect of electron-phonon coupling on cascade development by varying Pp across the large range of values
♦ Inhomogeneous □ Homogeneous
X(ps1)
Figure 12. The number of stable defect pairs at the conclusion of the cascade simulations as a function of electron-phonon coupling. The plotted measure is related to the damping coefficient by X = ep/Mi. The inhomogeneous thermostat, which allows for feedback of energy from the heated electron to the ions, results in prolonged annealing and a reduced defect yield. Note that the electronic stopping power (fis/Mi = 1 ps-1) is only applied in the inhomogeneous case. (Reprinted with permission from Rutherford and Duffy [9]. Copyright 2007 IOP Publishing.)
10000 9000 8000 7000 6000 5000 4000 3000 2000 1000 0
X (Ps-1)
Figure 13. The peak number of (mostly unstable) defects measured during the cascade evolution as a function of electron-phonon coupling (x = ep/Mi). At high coupling, the use of an inhomogeneous evolving electronic temperature distribution is seen to enhance the degree of ionic mixing. The cascade is prolonged by the feedback of energy from electrons to ions. As in figure 12, no simple interpretation of the data for x ^ 1 ps-1 is possible. (Reprinted with permission from Rutherford and Duffy [9]. Copyright 2007 IOP Publishing.)
found in the literature (see section 3.5). The overall trend for the effect of electron-phonon coupling on the number of stable defect pairs formed is for an increase at moderate values of fip due to freezing in the distribution of point defects and for a decrease at high fip because the removal of energy from the cascade reduces the size of the cascade and hence the number of point defects created.
A comparison between the use of a spatially varying electronic temperature and simulations using a homogeneous Langevin thermostat at 300 K (the perfect heat sink of Caro and Victoria [121]) shows that the former tends to reduce the stable defect yield across the full range of ep (see figure 12). Comparison of the maximum number of defects (most of which will heal) formed during the simulations shows that at moderate and high ep, inhomogeneity in the thermostat tends to increase the size and duration of the thermal spike (see figure 13). These effects can be understood as the result of the elevated electronic temperatures developed in the cascade region; the feedback of energy to the ions prolongs the thermal spike, increasing the peak defect yield, but also allowing for prolonged annealing, reducing the final stable defect yield.
We note that care should be taken in interpreting figures 12 and 13 for lower values of the coupling (fip/Mi < 1ps-1). In the simulations using a homogeneous thermostat, the electronic stopping power is omitted (в = 0) meaning that the effect of allowing the electrons to heat up is properly disentangled from that of a higher average damping only when
в » A.
In [9] the evolution of the electronic and ionic temperature distributions during cascade simulations is plotted and the effect of the electron-ion interaction can be seen in the dramatically faster decrease in the ionic temperature for high Pp.
Duffy et al [78] find a further application for their model in simulating the formation of damage tracks around high energy channelling ions. They take an MD cell of around one million ions, coupled as in [9] to a much larger electronic system. The electronic temperature in a central cylindrical region 2.52 nm across is elevated to 7.5 x 104 K to represent the energy that would have been transferred from an ion traversing the region and experiencing a stopping power of 10 keV nm-1. These simulations thus examine the thermal spike model of channelling ion damage and the authors explore the sensitivity of the defect yield not only to the value of the electron-phonon coupling, but also to values of the electronic specific heat and thermal conductivity. The defect yield is highly sensitive to ep, with higher coupling resulting in a greater number of stable defects. Increasing the electronic thermal conductivity Ke from the experimental value for Fe rapidly reduces the stable defect yield by reducing the size of the molten region. This result is consistent with experimental data for the track formation behaviour of a variety of metals (see [78]). However, the very low sensitivity to electronic specific heat found in the simulations is in contradiction to the experimental evidence.
Care should be taken in the comparison of these simulations with experiments on different metals because the same Fe potential is used in each case. The simulations therefore fail to take account of the variation of other material properties that will affect the experimental results: melting point in particular would be expected to have a large effect on defect yield. Also, the heat capacity and stopping power are not independent material properties. The heat capacity will be proportional to the electronic density of states at the Fermi level, D(eF), and the stopping power, being dependent on the number of available electronic transitions from occupied to unoccupied states, will tend to be proportional to [D(eF)]2.
Ivanov and Zhigilei have developed a similar model for an MD simulation of ions coupled to an evolving electronic temperature distribution. They apply their model to the study of melting in nickel films induced by short laser pulses [126] and the velocity of the melt front in nickel and gold [127]. The main difference between the work of Duffy et al and that of Ivanov and Zhigilei is the manner by which energy is returned from hot electrons to cooler ions. In the latter case this is implemented by making the damping term proportional to the difference between the local electronic and ionic temperatures (averaged over a coarse-graining cell),
Pi a (Ti - Te). (71)
It is thus possible for the drag force to accelerate ions in the presence of hotter electrons.
Duvenbeck et al [77, 79, 80] have developed a model in which the work done by a drag force on the ions acts as a source of heat in a diffusional model of electronic temperature. The aim of the work is to explore the behaviour of the electronic temperature at the surface in simulations of sputtering events. However, since the authors do not include a mechanism for return of energy from electrons to ions, rejecting this as a second-order effect, their model represents only a partial coupling of electron and ion evolution.
The model presented by Duffy etal is a significant advance over earlier classical MD simulations with simple Langevin dynamics at a fixed thermostat temperature. Even though the techniques discussed in the rest of this section provide a more realistic model of the electronic system, available computational resources will restrict simulations of large collision cascades (of hundreds of thousands of atoms or more) to using some sort of augmented MD model for the next decade at least. The inhomogeneous Langevin thermostat, then, can be seen as an important first step towards a molecular dynamics scheme optimized to capture the effects of the exchange of energy between electrons and ions. Initially, the primary role of more sophisticated models will be to identify and calibrate the necessary refinements to such a classical scheme.
We therefore lay out the key remaining gaps in the model of Duffy et al as a useful guide for future work in the field:
Thermal properties of electrons. As Duffy et al acknowledge, their treatment of the electronic thermal conductivity and heat capacity is only at the simplest level, with a fixed room-temperature value for the former and a tanh(Te) form for the latter, saturating at 3kB per atom for high Te. Many processes of radiation damage are such that the electrons and ions are strongly out of equilibrium, both within and between the subsystems, for significant periods of time and so more sophisticated treatments of heat transport might be required.
Ivanov and Zhigilei [126,127] implement a temperature dependent thermal conductivity of a form that tends to Ke ~ Te/TI at low temperatures and to Ke ~ Te5/2, characteristic of a low density plasma, at high Te.
Energy transfer from ions to electrons. Various ad hoc forms for an electronic damping of ionic motion have been proposed, but all have in common the fact that the damping force is
proportional to and directly opposed to the ionic velocity. Such forms are justified by an appeal to a large body of supporting literature, both experimental and theoretical. However, the theoretical work is carried out only within approximations that would be expected to yield a linear damping result and experiments are severely restricted in what they can measure. Neither has anything to say about the magnitude or importance of deviations from perfect linearity and opposition to velocity of the force on the ions due to the electrons.
Furthermore, most schemes make a distinction between the electronic stopping power regime and the electron-phonon coupling regime, between which a large difference in damping is expected. Bacon et al [84,85] and Nordlund et al [93,116118] choose to ignore one or the other regime entirely, Caro and Victoria propose an ad hoc fitted form based on local electronic density to accommodate both regimes and Duffy et al incorporate each via a separate constant force. We reiterate that the same physics of electron-ion interactions underlies energy transfer in both regimes [17]. The distinction between electronic stopping and electron-phonon coupling is normally drawn in terms of the correlation between the motion of different ions (small in the former case, large in the latter). But a molecular dynamics simulation, with its explicit treatment of ionic motion, contains full information about such correlation. This means that it should be possible to find a scheme that accounts for variation in the damping force in a physical way. The investigation of radiation damage with more advanced models should help to inform the design of such a scheme, incorporating a local environment and velocity dependent, anisotropic electronic damping force.
Ensuring the correct dynamics. Several potential problems arise with the inclusion of a spatially varying electronic temperature. The first of these is the correct handling of the boundary conditions on the electronic subsystem. In the first incarnation of their model [8], Duffy et al chose the same size for the electron and ion systems. All the energy transferred into the electrons was therefore trapped in a small region, resulting in anomalously high electronic temperatures which were only slowly reduced by the 300 K thermostat on the boundary. In later work [9] the electronic system is greatly enlarged, allowing for diffusive transport of energy out of the ionic MD region, but the system is still finite in extent and bounded by a simple 300 K thermostat. In fact, Duffy et al identify three phases in the evolution of the maximum electronic temperature: a rapid rise, followed by a rapid falloff dominated by electronic heat diffusion, followed by a final slower decay. It is not clear, but it is possible that the change in behaviour between the fast and slower decay is simply due to the establishment of a fairly uniform electronic temperature, with further cooling dominated by the thermostat.
A second issue is the handling of the temperatures outside the region of the MD simulation. Duffy et al implicitly assume that the ion temperature remains equal to the electronic temperature in the embedding region. In reality the ions will tend to remain cooler than the electrons, but will gradually heat up. In effect they provide a spatially varying bath for the electronic temperature in the embedding region and another heat sink, in addition to the boundary thermostat.
5.2. Adding quantum mechanical electrons
In the rest of this review we will consider a variety of simulation methodologies able to better capture the complex physics of radiation damage in metals. These methodologies vary widely in their complexity, but none of them can be claimed to provide a means for directly simulating radiation damage cascades over realistic time and length scales. Advances in computational resources will gradually render models in each step up the ladder of complexity applicable to direct simulation, but in the meantime the models should primarily be viewed as exploratory tools. This is not to belittle their importance. Whilst the inhomogeneous Langevin thermostat of Duffy et al captures much of the important physics in a phenomenological way, there remain many unanswered questions about the interaction of ions and electrons. Only by answering these questions can we hope to build more realistic molecular dynamics models. And an excellent way of answering them is by undertaking simulations in a framework of more sophisticated physics.
A non-relativistic quantum mechanical system of Nn nuclei at positions R = {R1, R2,..., RNn} and Ne electrons at positions r = {r1, r2,...,rNe} is properly described at time t by a many-body wavefunction $>(R, r, t). This wavefunction will evolve under a Hamiltonian H (R, r), which we assume to have no explicit time dependence, according to the time-dependent Schrodinger equation,
H(R, r)$(R, r, t) = ih — &(R, r, t).
H(R, r) incorporates terms corresponding to the kinetic energy of the nuclei, Tn (R), the kinetic energy of the electrons, Te (r), and to the electrostatic interaction energy of the nuclei with each other, Vnn(R), the electrons with each other Vee(r) and the nuclei with the electrons Vne(R, r),
H (R, r) = Tn (R) + Te(r)+Vnn (R)+Vee(r)+Vne(R, r). (73)
Solution of this full problem is computationally impossible for any system of appreciable size and progress can be made only by making simplifying approximations.
Most practical simulation schemes begin by making the semi-classical approximation, choosing to treat the relatively massive nuclei as classical particles. This then reduces the problem to one of calculating the evolution of the quantum mechanical electronic system, described in general by a many-body wavefunction ^(r, t), under an electronic Hamiltonian,
Of course, solution of the full many-body electronic problem still remains a challenge formidable to the point of computational intractability, and so the rest of the task of implementing a semi-classical simulation scheme consists in finding a suitable approximation to ^(r, t) and its evolution,
ih — V(r,t) = he(r; R(t))ty(r,t). (76)
Before we consider some specific examples of such schemes, we will briefly examine the dynamics of a semi-classical system. To make our analysis transparent we will choose to work in the basis of many-electron instantaneous eigenstates of the electronic Hamiltonian hie. Working from now on in the Dirac notation we will denote the eigenstates by {| fi(t)}} and their eigenvalues by {£¿(0} so that we have
he(r, R(t))l fj(r; R(t))} = ej| fj(r; R(t))}. (77)
We now expand our many-electron wavefunction in terms of these eigenstates,
| V(r, t)} = £aj(t)e-(l/h)/£jdt| fj(r; R(t))}, (78)
where the phase-factor exp[-(i/h )/ej dt ] is inserted for algebraic convenience. The expansion coefficients are given
aj(t) = e(l/h)iejdt{fj(r\ R(t))\^(r, t)). (79)
Directed by the expression for the Hellmann-Feynman force (75), we consider the quantity
( fi \y Rh e \ fj ) = VRejSij + ei (Vr ft \ fj ) + ej {ft \ Vr^j ),
Where we have used the fact that he is Hermitian. By the orthogonality of the eigenstates, (ft \ fj) = Stj, we have Vr ( ( ft \ fj )) = 0 ^ ( Vr ft \ fj ) = -(ft \ VRfj ) and so we can write the Hellmann-Feynman force as
Fe = -J2\at\2V Ret
+ a*iaj(et - ej)e
(i/h)f(e; -ej)dt
(ft \ Vrf j). (81)
Inserting expansion (78) into the time-dependent Schrodinger equation (76) andpremultiplying by ( ft(r; R(t)) \ yields the time evolution of the expansion coefficients,
he(r; R(t)) = t(r) + Vee(r) + Ke(r; R(t)) + Vnn(R), (74) dtat(t) = ^ ■ (ft \ VRfj), (82)
where the coordinates of the classical nuclei now enter as parameters in the electronic dynamics and determine the time dependence of the Hamiltonian. Vnn (R) is the interaction energy of the classical ions and is included to cancel a singularity in the electronic energy. In turn, the nuclei will respond to forces due to the electrons derived from the Hellmann-Feynman theorem [128],
Fe(R,t) = - j dr [^(r,t)]+VRhe(r; R(t))^(r,t). (75)
where we have made use of the chain rule, dfj/dt = (VRfj) ■ R.
If we consider starting our evolution at time t0 with the electronic system in its ground state, say \ ^(r, t0) ) = \ f1(r; R(t0)) ), such that at(t0) = St1, then (82) gives a picture of the evolution in which the nuclear motion stimulates transitions to excited states at a rate determined by the non-adiabatic coupling vectors (ft \ VRfj). These coupling vectors also appear in (81 ), where we can see that the electronic
Figure 14. A schematic representation of the variation of the energies E1j2 and occupation probabilities |a1212 of a pair of energy eigenstates as a function of an ionic coordinate R. We are assuming that R is changing sufficiently quickly that non-adiabatic transitions within the system are stimulated. As the probabilities evolve, the potential energy surface Eeh traversed by the Ehrenfest system is a weighted average, representative of neither eigenstate. (After [129].)
forces on the nuclei can be decomposed into motion on a set of potential energy surfaces feW} corresponding to the instantaneous eigenstates {| fi}} occupied according to the evolving expansion coefficients {ai(t)} and a set of non-adiabatic forces given by the second term in (81). These non-adiabatic forces act in the direction of the non-adiabatic coupling vectors and the ionic motion will do work against them equal to the electronic excitation energy implied by the changing expansion coefficients.
5.2.1. Surface hopping versus Ehrenfest dynamics. Most semi-classical simulation schemes implement the electronic evolution (82) with some approximate treatment of the electronic wavefunction. The nuclei are evolved under forces which include those given by (81). Such a scheme is referred to as Ehrenfest dynamics and is a mean-field approach in the sense that the ions move on a linear combination of the potential energy surfaces corresponding to the electronic energy eigenstates. A possible flaw in Ehrenfest dynamics becomes apparent if we consider the case of an excited eigenstate that has a very different shape to the ground state. We can regard the lai(t)|2 as representing the probability that at time t the electronic system has been excited into the ith eigenstate. Ehrenfest dynamics will correspond to motion on some weighted average potential energy surface (illustrated schematically in figure 14), which is representative of neither the ground state nor the excited state.
We could instead consider a dynamical evolution in which the nuclei all remain on a potential energy surface corresponding to a single electronic eigenstate. The evolving values of {|ai (t) |2} are then used to determine the probability that the electronic system should undergo a discrete 'hop' to a different eigenstate, changing the potential energy surface on
which the nuclei move. This is the surface hopping method of Tully and Preston [130].
Chief amongst the issues that a practical surface hopping algorithm must address are, first, the mechanism for adjusting the nuclear energy in order to conserve total energy when a hop takes place and, second, how (and how frequently) to select which hops are made. A reasonably strong consensus has been reached on the second issue, with most surface hopping implementations adopting the so-called fewest-switches approach [131]. The first issue is more complex, but there are strong arguments for adjusting the kinetic energies of the nuclei by altering the components of the velocities parallel to the non-adiabatic coupling vectors, since this is the direction in which the non-adiabatic forces, corresponding to the electronic excitations, lie. Problems can arise when there is a finite probability of an electronic transition whose energy exceeds the kinetic energy of the ions parallel to the non-adiabatic coupling vectors. Often such transitions are simply disallowed by the algorithm.
Several features of the surface hopping approach make it unsuitable as a means of simulating radiation damage in metals. First, such simulations will tend to involve large numbers of atoms and so calculating the electronic forces on the nuclei (as gradients of the eigenstate potential energy surfaces) may become prohibitively expensive depending on the model of the electronic system used. This fact, along with the stochastic nature of the method, which means that it is necessary to follow a large number of trajectories in order to adequately sample the hopping probabilities, restricts the successful use of surface hopping to simulations of small systems (e.g. molecular reactions [132] and the interaction of atomic clusters [133]). Second, the large number of electrons involved will lead to a dense spectrum of electronic eigenstates meaning that a large number of hopping probabilities must be calculated and sampled, again increasing the computational overhead.
Since in a metallic system we would not expect to see large differences in the shape of the potential energy surfaces corresponding to electronic eigenstates close in energy, the Ehrenfest dynamics method, with its mean-field treatment, currently represents a better compromise between accurately representing the evolution of coupled quantum mechanical electrons and classical nuclei and limiting the computational overhead. We consider several examples of Ehrenfest dynamics simulations of radiation damage events below, each adopting a different model of the electronic system.
5.3. Time-dependent density functional theory
Time-independent density functional theory (DFT) is well established as a method for calculating the ground state electronic energies of a large variety of systems [61] to within a few tens of an electron volt per atom. The time-dependent version of the formalism (TD-DFT) is a more recent innovation [134] and allows the calculation of certain excited state properties, in contrast to standard DFT. Access to information about excitations of the electronic system makes possible the linear response calculations discussed in section 2.3.3 and also
opens up the possibility of determining the evolution of a system of electrons under a time-dependent Hamiltonian.
The formalism of TD-DFT is broadly analogous to that of standard DFT. It can be shown that, for a given initial many-body electron state = ^(r,t = t0), the electronic densities p(x, t) and p'(x, t) corresponding to evolution of under external potentials V(x,t) and V'(x,t), respectively, will be different provided the potentials differ by more than a function of time only (here x is a position vector in three-dimensional space and the electron density corresponding to ^(r, t) is given by£ ¡J dr S(x — ri)|^(r, t)|2). This is the Runge-Gross theorem [134], the time-dependent version of the Hohenberg-Kohn theorem [135]. Once a particular choice is made for then a given external potential determines a unique density p(x,t) and it is possible to write any observable of the electronic system as a functional of p(x, t).
It can then be established, that for a system of Ne interacting electrons evolving under a potential Vint(x,t), there exists a unique potential Veff (x, t) that, when acting on a fictitious system of Ne non-interacting electrons, yields the same time-dependent density. The evolution of this fictitious system is governed by the time-dependent Kohn-Sham equations,
d /V2 Ik 7ttfi(x, t) = I--— + Veff (x, t) ) fi(x, t),
~1-1-1-1-r
p —» Li6 F
for a set of Ne single-particle wavefunctions {fi(x, t)}. We thus have a way of calculating the evolution of a set of interacting electrons, provided we have a way of determining the correct Veff (x, t).
More information about the theory and application of TD-DFT can be found in the many reviews on the subject (see, for example, [136,137]).
The application of TD-DFT to direct simulation of radiation damage events is hampered by the computational complexity of the method. First, the calculation of the electronic forces on the nuclei requires that time-dependent one-electron orbitals be calculated to high precision, which in turn demands many basis states per atom. Second, calculation of the electronic Hamiltonian is inherently time consuming, requiring many three-dimensional spatial integrals over the basis states (six-dimensional in cases where non-local pseudopotentials are employed in a local basis).
We are aware of only one example in the literature of an atomistic radiation damage simulation using TD-DFT. Pruneda et al [138] examined the case of protons and anti-protons penetrating the insulator lithium fluoride. Experimental data for the electronic stopping power of LiF for various penetrating particles suggest the existence of a threshold effect [139141]. Below a certain velocity (0.1 v0 for protons in LiF) the electronic stopping power drops nearly to zero, violating the standard dE/dx a v behaviour. This threshold effect is attributed to the existence of the band gap, which implies a minimum energy for excitations in the target material, in contrast to the case of metallic targets.
Pruneda et al [138] undertake a series of simulations in a 4 x 4 x 4 unit cell lattice of LiF (128 atoms) in which all the atoms are frozen at their perfect lattice sites and a proton or anti-proton is constrained to move down the centre of a
0.1 0.2 0.3 0.4 0.5 0. v (a.u.)
Figure 15. Electronic stopping power dE/dx as a function of particle velocity from time-dependent DFT simulations of channelling in LiF. The results for protons are shown as filled circles, anti-protons as empty circles. The crosses show results for protons when extra basis states are added along the channelling particle's path. The inset and other data are discussed in [138]. (Reprinted with permission from Pruneda et al [138]. Copyright 2007 American Physical Society.)
[110] channel at a fixed velocity. The simulations use the SIESTA TD-DFT code [142] with the adiabatic local density approximation (ALDA) to the exchange-correlation energy. The results of the simulations for projectile velocities up to 0.6 v0 are shown in figure 15 and a threshold effect at around the velocity suggested by experiment is clearly evident. If we regard the channelling particle as a periodic perturbation to the extended electronic system, with a frequency determined by the passage of the particle from one cell to the equivalent point in the next [143], then we can understand the threshold velocity as corresponding, via the passing frequency, to the minimum possible excitation in the electronic system, in this case the band gap.
Above the threshold, the stopping power at a given velocity will be dependent on the number of transitions within the electronic system with energies corresponding to the frequencies characteristic of the ionic motion. Such frequencies will often be very low—an iron atom with a kinetic energy of 100 eV typical of a collision cascade will have a velocity v = 0.19Afs—1 and, assuming a characteristic atomic separation of d = 2 A, we obtain a characteristic excitation energy of 2nkv/d = 0.4 eV. The number of transitions available will be dependent on the size of the system simulated (for low frequency excitations it will be proportional to the square of the electronic density of states at the Fermi level) and the number of bands and £-points, and so we must be wary of finite system size effects when drawing quantitative conclusions from simulations in small systems. These system size effects will be a particular problem in any attempt to use TD-DFT to study radiation damage, because the computational complexity of the method will restrict simulations to very small
cells. The mere fact that TD-DFT is able to accurately predict excited state energies is not enough to ensure accuracy in the calculated values of quantities which, like stopping power, involve transitions between states in a possibly extremely sparse eigenspectrum. So, whilst the accuracy of the calculated threshold velocity might be very good, depending as it does on a specific feature of the band structure, the calculated stopping powers above the threshold might be unreliable.
Pruneda et al state that their 128 atom super-cell with calculations at a single k-point was chosen after 'convergence tests', but also make it clear that the values found for the stopping power were not converged. They attribute a residual factor of ~2 discrepancy between their calculated stopping power and experimental values (after making an adjustment for the channelling geometry) partly to a 'finite-basis saturation effect at high velocities'. They experiment in the case of proton stopping by augmenting the basis set of the target with extra hydrogenic basis states every 0.5 A along the projectile's path and find a ~75% enhancement in the stopping power. These results correspond better to experiment, but the justification for adding the extra basis states is unclear. Indeed, any augmentation of a finite-basis set, so long as the added basis states were coupled to the ionic motion, might be expected to increase the stopping power by increasing the number of excitations through which the electrons could absorb ionic energy. We discuss finite system size effects further in section 5.4.
Pruneda et al also point out that errors in the electronic structure will yield discrepancies in the calculated stopping powers. Whilst these discrepancies will be dwarfed by finite system size effects in this case, they will be an ongoing problem until improved models for the fully time-dependent exchange-correlation energy are available.
Despite the foregoing discussion, there is no fundamental reason why simulations of small collision cascades (up to a few hundred moving atoms) should not be undertaken with TD-DFT at the present time. Improvements in computational power will gradually increase the size of systems which can be modelled and such simulations could yield further valuable information about band-structure dependent effects in radiation damage.
5.4. Time-dependent tight-binding (TDTB)
Time-dependent density functional theory is highly restricted in the size of systems that it can simulate, because of the computational complexity of the underlying description of the electronic structure. By choosing a more approximate model we can trade off the accuracy of the electronic structure for an ability to simulate much larger systems for much longer times. Semi-empirical tight-binding (SETB) models offer just this compromise.
Semi-empirical tight-binding can be regarded as an approximation to the local orbital implementation of density functional theory14 in which the real space integrals needed to calculate the Hamiltonian matrix elements and the overlap of basis states are replaced with parametrized functions (for
14 This latter type of model is sometimes referred to as ab initio tight-binding.
details see [119,144,145]). The effects of the exchange-correlation term are implicitly included in these parameters. The parameters of SETB models are then derived either from experimental data or via recourse to more sophisticated models, such as DFT. A time-dependent form of SETB has been derived by Todorov [146] using a Lagrangian formalism.
Recent work by Mason et al [147] has explored the effect of electronic excitations in simulations of radiation damage in metals using a simplified tight-binding model. Whilst there exist tight-binding models with parametrizations that reproduce electronic structures with high accuracy [148], such models require large numbers of basis states per atom. This makes them too unwieldy for use in large simulations and also hampers attempts to interpret the dynamical evolution of the electronic system. At this stage in the effort to understand the quantum mechanical nature of non-adiabatic energy exchange between ions and electrons in radiation damage processes it will be helpful to focus one prong of attack on capturing the behaviour of simple systems on the largest possible length and time scales. Mason et al justify their choice of simplified TDTB dynamics as being the simplest possible way to introduce an explicit system of quantum mechanical electrons into an atomistic simulation. In any case, band-structure dependent effects will still appear in the system dynamics, and whilst the tight-binding band structure is not an accurate representation of any real metal, the dependence of aspects of the dynamical evolution on well-defined features of the band structure can still be quantified and understood.
Mason et al [147] adopt the single s-orbital tight-binding model of Sutton et al [149], selecting a parametrization that reproduces the mechanical and structural properties of copper along with a metallic electronic structure. The model uses a basis of orthogonal local orbitals {| a }} atoms at positions R = {R1, R2,..., RN consistent Hamiltonian for the electronic system is written as
centred on the Na The non-self-
h e(R(t)) = £| a ) Y (|Ra(t) - Rß(t)\) < ß |,
where he is now a one-electron Hamiltonian and where Y(\x\) are hopping integrals varying as an inverse power of the internuclear separation |x|. The electronic system is represented by a single-particle density matrix p(t), which is initialized at time t0 to be
P(to) = J2 f(£i(to); T)\MR(A>)) X fi(R(to)) \, (85)
where {| fi (R(t))}} are instantaneous eigenstates of h e(t) with eigenvalues {£¿(0} and f(e; i, T) = {exp[(e - i)/kBTe] + 1}-1 is a Fermi-Dirac distribution at temperature Te, chemical potential The evolution of this density matrix is given by the Liouville equation, equivalent to the time-dependent Schrodinger equation,
h — p (t) = [h e (R(t)),p(t)]. at
Electronic forces on the ions are given by the Hellmann-Feynman theorem (the independent electron equivalent of (75))
1.E+01 1.E+00 1.E-01 1.E-02 1.E-03 1.E-04
1.E-05
к ♦
0.1 1 10 driving frequency Q (rad/fs )
4— T=10 к ■ T=100 К
* T=1000 к л T=1e4K
—4— T=1e5 К ■ - T=1 еб К
natural frequency
1.6E-13
1.6E-14
1.6E-15
1.6E-16 j?
1.6E-17
1.6E-18
1.6E-19
Figure 16. The damping coefficient f corresponding to the average drag force on a single oscillating ion in a perfect lattice. f is shown as a function of driving frequency for different initial electronic temperatures. At low frequencies, finite-size effects are evident in the erratic behaviour of the plots. Above Q ~ 100 rad fs—1 the oscillator energy exceeds the finite band-width of the TB model, no electronic transitions can be stimulated and the damping force goes to zero. (Reprinted with permission from Mason et al [147]. Copyright 2007 IOP Publishing.)
to give an ionic equation of motion, MR (t) = -Уд УЮп (R(t)) - Tr (p (t) Vnh e (R(t))\, (87)
where Vion is a repulsive ion-ion interaction and it is assumed (for notational simplicity only) that all the ions have the same mass MI. Tr{-} indicates a trace over a complete basis of the electronic system, easily computable in the local orbital basis. Now (86) and (87), along with the definition of the Hamiltonian (84), form a closed set of equations for the dynamics of the semi-classical system of ions and electrons.
The time-evolved density matrix p(t) encapsulates all the information about the quantum mechanical electronic system and allows all quantities of interest to be calculated. In addition, by constructing an adiabatic density matrix p0(t), which represents the state the electronic system would be in if the ions had traversed their paths infinitely slowly, we can calculate the irreversible transfer of energy from ions to electrons,
AE = Tr (p(t) - po(t))he(R(t))
Mason et al [147] have explored the energy transfer response of a 1120 atom perfect lattice super-cell of their model to a single forced oscillating ion as a function of the frequency, position and direction of the oscillator and of the electronic temperature. The energy transfer in each case gives a measure of the coefficient f of the drag force on the oscillating ion (cf the discussion in sections 4.4 and 3.3, in particular (63)). The results of these oscillator simulations are reproduced in figures 16 and 17 and reveal that in general the damping coefficient is a tensor quantity dependent on the local atomic environment, the velocity of the moving ion and the electronic temperature.
Reference [147] also contains a detailed treatment of the effect of finite system size on the calculated energy transfer. A time-dependent perturbation theory analysis shows that, in any
given system, the irreversible energy transfer to the electrons is only well behaved for a finite time. The key issue is the degree of resolution of the energy transitions stimulated by the oscillating ion. The time-energy uncertainty relation suggests this should vary as ~ T/t and at some point the resolution will become so fine that the oscillator can no longer stimulate transitions within the discrete spectrum of the finite system. This problem is a general one, which needs to be carefully addressed in any semi-classical simulation in a finite system.
Despite the complexity of the most general form of the damping coefficient в demonstrated by Mason et al [147], we might still hope that a simple velocity independent damping force (analogous to a stopping power linearly proportional to velocity) would do a good job of representing the average effect of electrons on the ion dynamics in collision cascades. Le Page et al [150] examine this possibility by undertaking a series of 240 cascade simulations with the model of Mason et al [147]. Simulation cells of 2016 atoms with periodic boundary conditions are used to simulate cascades with a range of PKA energies from 100 eV to 1 keV in 24 different directions. By calculating the work done as an integral over the ionic trajectories for different models of the damping force and comparing it with the irreversible energy transfer (88), le Page et al test the validity of various MD-with-damping simulation schemes. Figure 18 shows the results and it is clear that they offer strong support for the use of a simple velocity-independent damping coefficient as a first approximation to the true non-adiabatic force. Importantly, the model of Caro and Victoria [121] (see section 4.5) gives an improvement over a simple damping in certain circumstances. The density dependence of their model goes some way to capturing the factor of 2 increase found in the damping of replacement collision sequences in the simulations of le Page et al [150]. The results offer no support for the use of a lower kinetic energy cut-off for the application of the damping force of the type implemented by Nordlund et al [93,117] (see section 4.4).
Figure 17. As figure 16 except that the oscillating ion is moved away from its lattice site to the positions indicated and oscillated at a fixed frequency, Q = 1 rad fs-1. The damping for an off-site oscillator is anisotropic and is plotted for a selection of directions of oscillation. A rich structure emerges, with a variation in f of an order of magnitude with direction and position. (Reprinted with permission from Mason et al [147]. Copyright 2007 IOP Publishing.)
1000eV
<1 "a œ
<1 "a œ
scaled AEehr(t)
scaled AEehr(t)
Figure 18. Plots of the irreversible energy transfer from ions to electrons in semi-classical (TDTB) simulations of collision cascades. Data from 24 simulations at PKA energies of 100 ev (left panel) and 1 keV (right panel) are shown. In each plot, for each simulation, the energy calculated according to one of three classical damping models is plotted against the energy transfer from the TDTB simulations. Top row: simple fixed damping model with no cut-off. Middle row: fixed damping applied above a 10 eV ionic kinetic energy cut-off (cf the discussion of the work of Nordlund et al [93,116-118] in section 4.4). Bottom row: the local electron density dependent model of Caro and Victoria [121] (see section 4.5). (Reprinted with permission from le Page et al [150]. Copyright 2009 IOP Publishing.)
An advantage of the approach of Mason et al is that the nature of the electronic excitations stimulated by ionic motion can be analysed in detail. Race et al [14] apply the method to study several hundred cascades in 2016 atom super-cells with PKA energies from 1 keV up to 50 keV. They find that the spectrum of electronic excitations is well described by a Fermi-Dirac distribution at an elevated temperature (see figure 19), even though the electron dynamics in the simulations do
not include the direct electron-electron interactions able to thermalize a non-equilibrium distribution15. This property of the excitation spectrum can be attributed to the characteristic
15 In reality the electron-ion interactions would also act to thermalize the electronic system. However Ehrenfest dynamics does not include the effect of spontaneous phonon emission (see section 5.5) and so in this case the electron-ion interaction cannot produce thermalization of the electrons at a fixed electronic energy.
Occupations Original distribution Best-fit distribution
JT 0.6 n a
-3 -2 Energy (ev)
Figure 19. The occupations of the instantaneous eigenstates around the Fermi level 225 fs into a sample TDTB simulation of a 2016 atom cascade. The excitations are seen to be well modelled by a best-fit thermal function at Te = 6055 K despite the lack of thermalizing electron-electron interactions. The initial temperature of Te = 300 K is shown for comparison. (Reprinted with permission from Race et al [14]. Copyright 2009 IOP Publishing.)
spectrum of the ionic motion. Even for a 50keV cascade, most of the transitions stimulated are small on the scale of the width of the Fermi-Dirac distribution and so the evolution of the excitation spectrum takes the form of a one-dimensional diffusion in energy space. This strong evidence for a well-defined electronic temperature justifies an assumption made in two temperature models (see sections 3.2 and 5.1) often used to simplify the development of excitation-dependent potentials for use in MD simulations of hot materials [151]. It will also greatly simplify efforts to devise and implement an augmented MD scheme that correctly treats energy exchange between ions and electrons.
Information about the effect of electronic excitations on the interionic forces is also available and in principle the direction and magnitude of the non-adiabatic force (given by the second term in (81)) can be analysed in detail. Studies by Race et al [14] of the effect of the electronic excitation on the average bond strength in evolving cascades show that significant changes may occur in situations where electronic temperatures above 10 000 K develop (see figure 20). Classical MD simulations of such processes might then need to make use of electronic temperature-dependent potentials.
Whilst the simulations so far undertaken by Mason et al [14,147,150] have been of low energy cascades, there is no reason that their model or one similar should not be applied to much higher energy phenomena. One possibility would be to carry out simulations like those of Pruneda et al [138] (discussed in section 5.3) of high energy channelling ions, but in the much bigger system sizes made possible by the simpler treatment of electronic structure. Such large-scale semi-classical simulations are currently the only way to search for excitation phenomena dependent both on the detailed atomistic evolution and on the quantum mechanical nature of electrons.
Extrapolated data —I-Perfect crystal — Simulation data
Temperature (K)
100000
Figure 20. The percentage reduction in the average magnitude of the attractive electronic force with rising temperature in a TDTB model. Data shown are the following: (red) squares, reduction in the force found in cascade simulations up to 50 keV PKA energy, plotted against the best-fit temperature to the excited electronic system. (Purple) circles, the calculated reduction in the force based on the temperature dependence of the bond orders in a perfect crystal. (Blue) crosses, the change in force found for an electronic temperature in a static system distorted to reflect the ionic positions during cascade simulations. (Reprinted with permission from Race et al [14]. Copyright 2009 IOP Publishing.)
5.5. Correlated electron ion dynamics (CEID)
The approximation of the ions as classical particles, inherent in Ehrenfest dynamics (ED), prevents the equilibration of the electronic and ionic subsystems. ED accurately reproduces the transfer of energy from excited ions into cooler electrons. However, it fails to produce thermal equilibrium between electrons and ions as the spontaneous emission of phonons is suppressed by the mean-field approximation [152-157]. This asymmetry exists because each ion is treated explicitly (so their fluctuations are visible to the electrons which can thus identify the ionic temperature) while the electrons are experienced by the ions as a structureless fluid whose temperature cannot be identified through the forces. This is a completely general property of the Ehrenfest approximation, independent of the level of description of the electrons.
The inability of ED to reproduce spontaneous phonon emission and so give the correct long-term behaviour of a system of ions and electrons is not always a problem in simulations of radiation damage. In the early stages of collision cascades and in phenomena such as channelling, the initial conditions of the combined system place a large excess of energy in the motion of the ions. The predominant mode of energy exchange is therefore from ions to electrons, exactly the process that ED can model well [156].
However, many key open questions in radiation damage centre on the electron-phonon coupling and the effects of energy exchange in both directions as the electrons and ions approach equilibrium. Modelling the effect of the electrons in quenching in or inhibiting the formation of defects in the displacement spike is critically dependent on a correct
Current computational capacity
Millions of atoms with statistics
Millions of atoms
Thousands of atoms
Hundreds of atoms
A few QM atoms
Binary collision / Continuum model
Non-perturbative representation
Direct simulation
Explicit ionic evolution
Implicit electronic effects
Explicit quantum electrons
Accurate electronic structure
Quantum mechanical Ions
Section 5.5
Figure 21. An overview of the theories and models covered in this review. The types of model are represented in terms of their incremental physical content along with an indication of the size of system that can be modelled with current computational resources. Augmented classical molecular dynamics models represent a current 'best compromise', able to handle simulations of millions of atoms, but having the potential to include much of the key physics of energy exchange between ions and electrons. Such models are informed by analytical stopping power and electron-phonon coupling theories and by the results of simulations with more complex models.
treatment of the full electron-phonon coupling, including spontaneous emission of phonons.
Correlated electron-ion dynamics (CEID) is a systematic method for extending ED so as to reintroduce the correct flow of energy from the excited electrons back to the ions. This is achieved by including small quantum fluctuations in the ionic trajectories through a low order moment expansion [153,154]. This allows the ions to probe the response of the electrons to small changes in their trajectory from which information about the internal state of the electrons can be determined. The moments correspond to powers of the instantaneous ionic positions and momenta relative to the mean values for the trajectories. The first moment gives us heating [153], while the second moment is needed to introduce the scattering of electrons by heated ions [154]. In order to capture strong non-adiabatic effects the formalism has had to be reworked to make the expansion stable. This has been achieved by the introduction of an efficient basis set expansion [158].
The ability of CEID to describe correctly the transfer of energy between electrons and ions over a wide range of electron-phonon coupling strengths has now been demonstrated: it has been applied to heating in current carrying nanowires [153,154,159,160] and simple two level systems [158]. Thus in principle it can be applied to problems in radiation damage, including the equilibration stage after the initial cascade
formation. Just as for ED, we anticipate that CEID will be used to inform much simpler and computationally efficient models for large-scale simulations. An intermediate step might be to use CEID to support simple corrections to ED simulations to enable them to reach thermal equilibrium. However, the complexity of CEID calculations and their memory requirements are much larger than for ED, which currently makes even reference simulations very expensive. Reducing the computational burden is an ongoing research topic.
6. Concluding remarks
Our aim in this review has been to outline and explain the various ways in which electronic effects can be incorporated within atomistic simulations of radiation damage in metals (see figure 21). We began by discussing the theoretical treatments of electronic stopping power and electron-phonon coupling. From our point of view these represent a means of understanding and quantifying the role of electrons in the various events that make up a radiation damage process. They tell us about how ionic collision dynamics are altered when excitation of the electrons is taken into account or how a ballistically moving ion will lose energy to the surrounding electron gas.
Moving on we considered various approaches to the simulation of radiation damage, in which a series of interactions between classical ions is allowed to unfold dynamically according to some approximate physical model. We discussed the various ways in which electronic effects have been incorporated into these simulations, always informed by stopping power and electron-phonon coupling theory. The electrons might make themselves felt as an inelastic loss in ionic collisions (as in the BCA) or they might manifest themselves as a viscous medium providing a drag force on the ions, or as a stochastic force, buffeting the ions and returning energy to them. We discussed the application of such models and what has been learnt from them—above all, that defect yields can be materially affected by the electrons.
Yet in discussing the full range of stopping power theory we saw how the effect of electrons becomes increasingly complex as we consider slower or heavier ions. With this increasing complexity the literature becomes increasingly sparse, in contrast to the comprehensive treatment of fast, light ions. The theories of slow ion stopping, all yielding a stopping force proportional to ion velocity, are not expected to remain valid down to the energies involved in many damage scenarios of technological importance. Somewhere in the murk of this low energy regime, only dimly lit by experimental data, analytical theory switches its attention to the electron-phonon coupling. The interactions of ions and electrons are treated within a different set of approximations to produce estimates of the rate of energy exchange that are, as we saw, inconclusive. The various different theoretical treatments conflict with one another and with experiment. This uncertainty, along with the expected low energy failure of stopping theory, demands new approaches to investigating electronic effects.
There are several reasons that we might expect simple theories to fail in the low energy regime that encompasses much of the evolution of a radiation damage cascade. When ions move more slowly it is harder to view the evolution of the cascade as a series of separable binary encounters. Equally we expect that nuclear and electronic energy losses will become correlated and that it will no longer be valid to treat the electrons as a homogeneous stopping medium. At lower velocities too, the excitations of the electronic system stimulated by the ionic motion will approach the energy scale of details in the electronic structure. These factors point the way towards the latest generation of atomistic models of radiation damage. We concluded our review by considering models that combine a dynamically evolving set of ions and an explicitly represented system of quantum mechanical electrons. Such models allow for the possibility of incorporating the effects of electronic structure and of the many-body nature of the cascade evolution. We presented the results of simulations undertaken to date within such a framework.
Finally, we gave some thought to the future of radiation damage simulations. As available computational resources improve, new techniques, incorporating some of the effects of the quantum mechanical nature of the ions, will become useful in investigating radiation damage phenomena. Such effects must be taken into account in a full description of the later stages of a collision cascade when hot electrons
may be in contact with cooler ions. At the same time the semi-classical techniques discussed above will become applicable to ever larger length and time scales, ultimately allowing for quantitative investigations via the simulation of full cascades within an accurate model of the electronic structure. In the short-term, however, the main role of more advanced techniques will be to inform the design of better classical models allowing the impact of electronic excitations on systems of millions of atoms to be more fully explored.
Acknowledgment
CPR and DRM were supported by EPSRC grant number EP/C524403.
References
[1] Greenwood L R 1994 Neutron interactions and atomic recoil
spectra J. Nucl. Mater. 216 29-44
[2] Stoneham A M 1997 Finding the gaps: problems in radiation
damage theory Radiat. Effects Defects Solids 142 191-203
[3] Stoneham A M 1990 Energy transfer between electrons and
ions in collision cascades in solids Nucl. Instrum. Methods Phys. Res. B 48 389-98
[4] Bohr N 1948 The penetration of atomic particles through
matter Mat. Fys. Med.—K. Dan. Vid. Sel. 18
[5] Ziegler J F, Biersack J P and Littmark U 1985 The Stopping
and Range of Ions in Solids (New York: Pergamon)
[6] Sigmund P 2006 Particle Penetration and Radiation
Effects—General Aspects and Stopping of Swift Point Charges (Berlin: Springer)
[7] Domeij B, Brown F, Davies J A, Piercy G R and
Kornelsen E V 1964 Anomalous penetration of heavy ions of kev energies in monocrystalline tungsten Phys. Rev. Lett. 12 363-6
[8] Duffy D M and Rutherford A M 2007 Including the effects of
electronic stopping and electron-ion interactions in radiation damage simulations J. Phys.: Condens. Matter 19 016207
[9] Rutherford A M and Duffy D M 2007 The effect of
electron-ion interactions on radiation damage simulations J. Phys.: Condens. Matter 19 496201
[10] Fleischer R L, Price P B and Walker R M 1965 Ion explosion
spike mechanism for formation of charged-particle tracks in solids J. Appl. Phys. 36 3645-52
[11] Ryazanov A I, Pavlov S A, Metelkin E V and Zhemerev A V
2005 Effect of coulomb explosion on track formation in metals irradiated by heavy ions J. Exp. Theor. Phys. 101 120-7
[12] Fleischer R L, Price P B, Walker R M and Hubbard E L 1967
Criterion for registration in dielectric track detectors Phys. Rev. 156 353-5
[13] Itoh N and Stoneham A M 2001 Materials Modification by
Electronic Excitation (Cambridge: Cambridge University Press)
[14] Race C P, Mason D R and Sutton A P 2009 Electronic
excitations and their effect on the interionic forces in simulations of radiation damage in metals J. Phys.: Condens. Matter 21 115702
[15] Averback R S and Diaz de la Rubia T 1997 Displacement
damage in irradiated metals and semiconductors Solid State Phys: Adv. Res. Appl. 51 281-402
[16] Thompson M W 1969 Defects and Radiation Damage in
Metals (Cambridge: Cambridge University Press)
[17] Koponen I 1993 Energy transfer between electrons and
ions in dense displacement cascades Phys. Rev. B 47 14011-9
[18] Northcliffe L C 1963 Passage of heavy ions through matter
Annu. Rev. Nucl. Sci. 13 67-102
[19] Seitz F and Koehler J S 1956 Displacement of atoms during
irradiation Solid State Phys.: Adv. Res. Appl. 2 305-448
[20] Baurichter A, Sigmund P and Sorensen A H 2002 Panel
discussion on stopping of heavy ions Nucl. Instrum. Methods Phys. Res. B 195 224-31
[21] Thomson J J 1912 Ionization by moving electrified particles
Phil. Mag. 23 449-57
[22] Darwin C G 1912 A theory of the absorption and scattering
of the a rays Phil. Mag. 13 901-20
[23] Bohr N 1912 On the theory of the decrease of velocity of
moving electrified particles on passing through matter Phil. Mag. 23 10-31
[24] Bethe H 1930 Zur theorie des durchgangs schneller
korpuskularstrahlen durch materie Ann. Phys., Lpz. 397 325-400 (in German)
[25] Paul H Stopping power for light ions
http://www.exphys.uni-linz.ac.at/stopping/
[26] Bloch F 1933 Bremsvermögen von atomen mit mehreren
elektronen Z. Phys. Hadrons Nuclei 81 363-76
[27] Schiff L I 1955 Quantum Mechanics (New York:
McGraw-Hill)
[28] Ziegler J F 1999 Stopping of energetic light ions in elemental
matter J. Appl. Phys. 85 1249-72
[29] Lindhard J and Sorensen A H 1996 Relativistic theory of
stopping for heavy ions Phys. Rev. A 53 2443-56
[30] Sigmund P 2000 Shell correction in bohr stopping theory Eur.
Phys. J. D—At. Mol. Opt. Plasma Phys. 12 111-6
[31] Fano U 1963 Penetration of protons, alpha particles, and
mesons Annu. Rev. Nucl. Sci. 13
[32] Messiah A 1999 Quantum Mechanics (New York: Dover)
[33] Smith F M, Birnbaum W and Barkas W H 1953
Measurements of meson masses and related quantities Phys. Rev. 91 765-6
[34] Barkas W H, Dyer J N and Heckman H H 1963 Resolution of
the a --mass anomaly Phys. Rev. Lett. 11 26-8
[35] Andersen H H, Simonsen H and Srensen H 1969 An
experimental investigation of charge-dependent deviations from the bethe stopping power formula Nucl. Phys. A 125 171-5
[36] Ashley J C, Ritchie R H and Brandt W 1972 z? effect in the
stopping power of matter for charged particles Phys. Rev. B 5 2393-7
[37] Lindhard J 1976 The barkas effect - or z3, z4 -corrections to
stopping of swift charged particles Nucl. Instrum. Methods 132 1-5
[38] Grüner F, Bell F, Assmann W and Schubert M 2004
Integrated approach to the electronic interaction of swift heavy ions with solids and gases Phys. Rev. Lett. 93 213201
[39] Sigmund P 2004 Stopping of Heavy Ions—A Theoretical
Approach (Berlin: Springer)
[40] Northcliffe L C 1960 Energy loss and effective charge of
heavy ions in aluminum Phys. Rev. 120 1744-57
[41] Yarlagadda B S, Robinson J E and Brandt W 1978
Effective-charge theory and the electronic stopping power of solids Phys. Rev. B 17 3473-83
[42] Brandt W and Kitagawa M 1982 Effective stopping-power
charges of swift ions in condensed matter Phys. Rev. B 25 5631-7
[43] Sigmund P 1997 Charge-dependent electronic stopping
of swift nonrelativistic heavy ions Phys. Rev. A 56 3781-93
[44] Sigmund P and Schinner A 2000 Binary stopping theory for
swift heavy ions Eur. Phys. J. D 12 425-34
[45] Sigmund P and Schinner A 2002 Binary theory of electronic
stopping Nucl. Instrum. Methods Phys. Res. B 195 64-90
[46] Grande P L and Schiwietz G 1998 Impact-parameter
dependence of the electronic energy loss of fast ions Phys. Rev. A 58 3796-801
[47] Schiwietz G and Grande P L 1999 A unitary convolution
approximation for the impact-parameter dependent electronic energy loss Nucl. Instrum. Methods Phys. Res. B 153 1-9
[48] Schiwietz G 1990 Coupled-channel calculation of stopping
powers for intermediate-energy light ions penetrating atomic h and he targets Phys. Rev. A 42 296-306
[49] Grande P L and Schiwietz G 2002 The unitary convolution
approximation for heavy ions Nucl. Instrum. Methods Phys. Res. B 195 55-63
[50] Blazevic A, Bohlen H G and von Oertzen W 2000
Charge-state changing processes for Ne ions passing through thin carbon foils Phys. Rev. A 61 32901
[51] Blazevic A, Bohlen H G and von Oertzen W 2002 Stopping
power of swift neon ions in dependence on the charge state in the non-equilibrium regime Nucl. Instrum. Methods Phys. Res. B 190 64-8
[52] Firsov O B 1959 A quantitative interpretation of the mean
electronic excitation energy in atomic collisions Sov. Phys.—JETP 36 1076
[53] Lindhard J and Scharff M 1961 Energy dissipation by ions in
the keV region Phys. Rev. 124 128-30
[54] Tilinin I S 1995 Quasiclassical expression for inelastic energy
losses in atomic particle collisions below the Bohr velocity Phys. Rev. A 51 3058-65
[55] Fermi E and Teller E 1947 The capture of negative mesotrons
in matter Phys. Rev. 72 399-408
[56] Lindhard J 1954 On the properties of a gas of charged
particles Mat. Fys. Medd.—K. Dan. Vid. Sel. 28
[57] Echenique P M, Flores F and Ritchie R H 1990 Dynamic
screening of ions in condensed matter Solid State Phys.: Adv. Res. Appl. 43 229-308
[58] Ritchie R H 1959 Interaction of charged particles
with a degenerate fermi-dirac electron gas Phys. Rev. 114 644-54
[59] Lindhard J and Winther A 1964 Stopping power of electron
gas and equipartition rule Mat. Fys. Medd.—K. Dan. Vid. Sel. 34 1-21
[60] Almbladh C O, von Barth U, Popovic Z D and Stott M J 1976
Screening of a proton in an electron gas Phys. Rev. B 14 2250-4
[61] Kohanoff J 2006 Electronic Structure Calculations for Solids
and Molecules: Theory and Computational Methods (Cambridge: Cambridge University Press)
[62] Echenique P M, Nieminen R M, Ashley J C and Ritchie R H
1986 Nonlinear stopping power of an electron gas for slow ions Phys. Rev. A 33 897-904
[63] Arista N R 2002 Energy loss of ions in solids: non-linear
calculations for slow and swift ions Nucl. Instrum. Methods Phys. Res. B 195 91-105
[64] Campillo I, Pitarke J M and Eguiluz A G 1998 Electronic
stopping power of aluminum crystal Phys. Rev. B 58 10307-14
[65] Pitarke J M and Campillo I 2000 Band structure effects on
the interaction of charged particles with solids Nucl. Instrum. Methods Phys. Res. B 164-165 147-60
[66] Dorado J J and Flores F 1993 Molecular-orbital theory for
the stopping power of atoms in the low-velocity regime: the case of helium in alkali metals Phys. Rev. A 47 3062-72
[67] Ashcroft N W and Mermin N D 1976 Solid State Phys.
(London: Thomson Learning)
[68] Ziegler J F SRIM—the stopping and range of ions in matter
http://www.srim.org
[69] Paul H and Schinner A 2001 An empirical approach to the
stopping power of solids and gases for ions from 3Li to 18Ar Nucl. Instrum. Methods Phys. Res. B 179 299-315
[70] Paul H and Schinner A 2002 An empirical approach to the
stopping power of solids and gases for ions from 3Li to 18Ar: II Nucl. Instrum. Methods Phys. Res. B 195166-74
[71] Lindhard J and Scharff M 1953 Energy loss in matter by fast
particles of low charge Mat. Fys. Med.—K. Dan. Vid. Sel. 27 1-31
[72] Paul H and Schinner A 2003 Judging the reliability of
stopping power tables and programs for heavy ions Nucl. Instrum. Methods Phys. Res. B 209 252-8
[73] Grimvall G 1981 The Electron-Phonon Interaction in Metals
(Oxford: North-Holland)
[74] Ziman J M 2001 Electrons and Phonons (Oxford: Oxford
University Press)
[75] Flynn C P and Averback R S 1988 Electron-phonon
interactions in energetic displacement cascades Phys. Rev.B 38 7118-20
[76] Ashcroft N W and Mermin N D 1976 Solid State Physics
(London: Thomson Learning) p 52
[77] Duvenbeck A and Wucher A 2005 Low-energy electronic
excitation in atomic collision cascades: a nonlinear transport model Phys. Rev. B 72 165408
[78] Duffy D M, Itoh N, Rutherford A M and Stoneham A M
2008 Making tracks in metals J. Phys.: Condens. Matter 20 082201
[79] Duvenbeck A, Sroubek Z and Wucher A 2005 Electronic
excitation in atomic collision cascades Nucl. Instrum. Methods Phys. Res. B 228 325-9
[80] Duvenbeck A, Sroubek F, Sroubek Z and Wucher A 2004
Computer simulation of low-energy electronic excitations in atomic collision cascades Nucl. Instrum. Methods Phys. Res. B 225 464-77
[81] Del Fatti N, Voisin C, Achermann M, Tzortzakis S,
Christofilos D and Vallee F 2000 Nonequilibrium electron dynamics in noble metals Phys. Rev. B 61 16956-66
[82] Qiu T Q and Tien C L 1992 Short-pulse laser heating on
metals Int. J. Heat Mass Transfer 35 719-26
[83] Finnis M W, Agnew P and Foreman A J E 1991 Thermal
excitation of electrons in energetic displacement cascades Phys. Rev. B 44 567-74
[84] Gao F, Bacon D J, Flewitt P E J and Lewis T A 1998 The
effects of electron-phonon coupling on defect production by displacement cascades in a-iron Modelling Simul. Mater. Sci. Eng. 6 543-56
[85] Kapinos V G and Bacon D J 1994 Influence of ion-electron
interaction on the formation mechanism of depleted zones in displacement cascades in metals Phys. Rev. B 50 13194-203
[86] Kapinos V G and Bacon D J 1996 Effect of melting and
electron-phonon coupling on the collapse of depleted zones in copper, nickel, and a-iron Phys. Rev. B 53 8287-95
[87] Allen P B 1987 Theory of thermal relaxation of electrons in
metals Phys. Rev. Lett. 59 1460-3
[88] Kaganov M I, Lifshitz I M and Tanatorov L V 1956
Relaxation between electrons and the crystalline lattice J. Exp. Theor. Phys. 4 173-8
[89] Wang Z G, Dufour C, Paumier E and Toulemonde M 1994
The se sensitivity of metals under swift-heavy-ion irradiation: a transient thermal process J. Phys.: Condens. Matter 6 6733-50
[90] Brorson S D, Kazeroonian A, Moodera J S, Face D W,
Cheng T K, Ippen E P, Dresselhaus M S and Dresselhaus G 1990 Femtosecond room-temperature measurement of the electron-phonon coupling constant
k in metallic superconductors Phys. Rev. Lett. 64 2172-5
[91] Elsayed-Ali H E, Norris T B, Pessot M A and Mourou G A
1987 Time-resolved observation of electron-phonon relaxation in copper Phys. Rev. Lett. 58 1212-5
[92] Bennemann K H 2004 Ultrafast dynamics in solids J. Phys.:
Condens. Matter 16 R995-R1056
[93] Nordlund K, Ghaly M, Averback R S, Caturla M, Diaz de la
Rubia T and Tarus J 1998 Defect production in collision cascades in elemental semiconductors and fcc metals Phys. Rev. B 57 7556-70
[94] Wayne King E, Merkle K L and Meshii M 1983 Threshold
energy surface and frenkel pair resistivity for Cu J. Nucl. Mater. 117 12-25
[95] Foreman A J E, English C A and Phythian W J 1992
Molecular dynamics calculations of displacement threshold energies and replacement collision sequences in copper using a many-body potential Phil. Mag. A 66 655-69
[96] Wayne King E and Benedek R 1981 Computer simulation
study of the displacement threshold-energy surface in Cu Phys. Rev. B 23 6335-9
[97] Wallner G, Anand M S, Greenwood L R, Kirk M A,
Mansell W and Waschkowski W 1988 Defect production rates in metals by reactor neutron irradiation at 4.6 K J. Nucl. Mater. 152 146-53
[98] Wei C-Y and Seidman D N 1981 The spatial distribution of
self-interstitial atoms around depleted zones in tungsten ion-irradiated at 10 K Phil. Mag. A 43 1419-39
[99] Kinchin G H and Pease R S 1955 The displacement of atoms
in solids by radiation Rep. Prog. Phys. 18 1-51
[100] Norgett M J, Robinson M T and Torrens I M 1975 A
proposed method of calculating displacement dose rates Nucl. Eng. Des. 33 50-4
[101] Lindhard J, Scharff M and Schott H E 1963 Range concepts
and heavy ion ranges (notes on atomic collisions, ii). Mat. Fys. Med.—K. Dan. Vid. Sel. 33 1-42
[102] Robinson M T and Torrens I M 1974 Computer simulation
of atomic-displacement cascades in solids in the binary-collision approximation Phys. Rev. B 9 5008-24
[103] Masayuki Yoshida 1961 Distribution of interstitials and
vacancies produced by an incident and fast neutron J. Phys. Soc. Japan 16 44-50
[104] Oen O S, Holmes D K and Robinson M T 1963 Ranges of
energetic atoms in solids J. Appl. Phys. 34 302-12
[105] Oen O S and Robinson M T 1964 Monte Carlo range
calculations for a Thomas-Fermi potential J. Appl. Phys. 35 2515-21
[106] Beeler J R and Besco D G 1963 Range and damage effects of
tunnel trajectories in a wurtzite structure J. Appl. Phys. 34 2873-8
[107] Robinson M T and Oen O S 1963 The channeling of
energetic atoms in crystal lattices Appl. Phys. Lett. 2 30-2
[108] Robinson M T and Ordean S O 1963 Computer studies of the
slowing down of energetic atoms in crystals Phys. Rev. 132 2385-98
[109] Piercy G R, Brown F, Davies J A and McCargo M 1963
Experimental evidence for the increase of heavy ion ranges by channeling in crystalline structure Phys. Rev. Lett. 10 399-400
[110] Beeler J R 1966 Vacancy and interstitial cluster production in
neutron-irradiated a-iron J. Appl. Phys. 37 3000-9
[111] Beeler J R 1964 Primary damage state in neutron-irradiated
iron J. Appl. Phys. 35 2226-36
[112] Beeler J R 1966 Displacement spikes in cubic metals: I.
a-iron, copper, and tungsten Phys. Rev. 150 470-87
[113] Gibson J B, Goland A N, Milgram M and Vineyard G H 1960
Dynamics of radiation damage Phys. Rev. 120 1229-53
120 121
Malerba L 2006 Molecular dynamics simulation of [137]
displacement cascades in Fe: a critical review J. Nucl. Mater. 351 28-38 [138]
Diaz de la Rubia T 1996 Irradiation-induced defect production in elemental metals and semiconductors: a review of recent molecular dynamics studies Annu. Rev. [139] Mater. Sci. 26 613-49 Nordlund K, Ghaly M and Averback R S 1998 Mechanisms of ion beam mixing in metals and semiconductors J. Appl. Phys. 83 1238-46 [140]
Nordlund K, Wei L, Zhong Y and Averback R S 1998 Role of electron-phonon coupling on collision cascade development in Ni, Pd and Pt Phys. Rev. B 57 R13965-8 Zhong Y, Nordlund K, Ghaly M and Averback R S 1998 [141]
Defect production in tungsten: a comparison between field-ion microscopy and molecular-dynamics simulations Phys. Rev. B 58 2361-4 Finnis M W 2003 Interatomic Forces in Condensed Matter [142]
(Oxford: Oxford University Press) Finnis M W and Sinclair J E 1984 A simple empirical n-body
potential for transition metals Phil. Mag. A 50 45-55 Caro A and Victoria M 1989 Ion-electron interaction [143]
in molecular-dynamics cascades Phys. Rev. A 40 2287-91 [144]
Kitagawa M and Ohtsuki Y H 1974 Inelastic scattering of
slow ions in channeling Phys. Rev. B 9 4719-23 Echenique P M, Nieminen R M and Ritchie R H 1981 [145]
Density functional calculation of stopping power of an electron gas for slow ions Solid State Commun. 37 779-81 [146]
Pronnecke S, Caro A, Victoria M, Diaz de la Rubia T and
Guinan M W 1991 The effect of electronic energy loss on [147] the dynamics of thermal spikes in Cu J. Mater. Res. 6 483-91
Dudarev S L and Derlet P M 2005 A 'magnetic' interatomic
potential for molecular dynamics simulations J. Phys.: [148]
Condens. Matter 17 7097-118 Ivanov D S and Zhigilei L V 2003 Combined
atomistic-continuum modeling of short-pulse laser
melting and disintegration of metal films Phys. Rev. B [149]
68 64114
Ivanov D S and Zhigilei L V 2007 Kinetic limit of heterogeneous melting in metals Phys. Rev. Lett. 98 195701 [150]
Feynman R P 1939 Forces in molecules Phys. Rev. 56 340-3
Doltsinis N L and Marx D 2002 First principles molecular
dynamics involving excited states and nonadiabatic [151]
transitions J. Theor. Comput. Chem. 1 319-49 Tully J C and Preston R K 1971 Trajectory surface hopping approach to nonadiabatic molecular collisions: the reaction of H+ with D2 J. Chem. Phys. 55 562-72 Tully J C 1990 Molecular dynamics with electronic [152]
transitions J. Chem. Phys. 93 1061-71 Tully J C 1974 Collisions of F(2P1/2) with H2 J. Chem. Phys. 60 3042-50
Faeder J, Delaney N, Maslen P E and Parson R 1997 Charge
flow and solvent dynamics in the photodissociation of [153]
cluster ions: a nonadiabatic molecular dynamics study of I- • Arn Chem. Phys. Lett. 270 196-205 Erich Runge and Gross E K U 1984 Density-functional
theory for time-dependent systems Phys. Rev. Lett. [154]
52 997
Hohenberg P and Kohn W 1964 Inhomogeneous electron gas
Phys. Rev. 136 B864-71 Gross E K U, Dobson J F and Petersilka M 1996 Density [155]
Functional Theory of Time-Dependent Phenomena (Topics in Current Chemistry vol 181) (Berlin: Springer)
Joubert D (ed) 1997 Density Functionals: Theory and
Applications (Berlin: Springer) Pruneda J M, Sanchez-Portal D, Arnau A, Juaristi J I and Artacho E 2007 Electronic stopping power in LiF from first principles Phys. Rev. Lett. 99 235501 Auth C, Mertens A, Winter H and Borisov A 1998 Threshold in the stopping of slow protons scattered from the surface of a wide-band-gap insulator Phys. Rev. Lett. 81 4831-4
Draxler M, Chenakin S P, Markin S N and Bauer P 2005 Apparent velocity threshold in the electronic stopping of slow hydrogen ions in LiF Phys. Rev. Lett. 95 113201
Markin S N, Primetzhofer D and Bauer P 2009 Vanishing electronic energy loss of very slow light ions in insulators with large band gaps Phys. Rev. Lett. 103 113201
Tsolakidis A, Sanchez-Portal D and Martin R M 2002 Calculation of the optical response of atomic clusters using time-dependent density functional theory and local orbitals Phys. Rev. B 66 235416 Artacho E 2007 Electronic stopping in insulators: a simple
model J. Phys.: Condens. Matter 19 275211 Matthew W, Foulkes C and Haydock R 1989 Tight-binding models and density-functional theory Phys. Rev. B 39 12520-36
Sutton A P, Finnis M W, Pettifor D G and Ohta Y 1988 The tight-binding bond model J. Phys. C: Solid State Phys. 21 35-66
Todorov T N 2001 Time-dependent tight binding J. Phys.:
Conden. Matter 13 10125-48 Mason D R, le Page J, Race C P, Foulkes W M C,
Finnis M W and Sutton A P 2007 Electronic damping of atomic dynamics in irradiation damage of metals J. Phys.: Condens. Matter 19 436209 Papaconstantopoulos D A and Mehl M J 2003 The
Slater-Koster tight-binding method: a computationally efficient and accurate approach J. Phys.: Condens. Matter 15R413-40
Sutton A P M, Todorov T N, Cawkwell M J and Hoekstra J 2001 A simple model of atomic interactions in noble metals based explicitly on electronic structure Phil. Mag. A 81 1833-48 le Page J, Mason D R, Race C P and Foulkes W M C 2009 How good is damped molecular dynamics as a method to simulate radiation damage in metals? New J. Phys. 11 013004
Khakshouri S, Alfe D and Duffy D M 2008 Development of an electron-temperature-dependent interatomic potential for molecular dynamics simulation of tungsten under electronic excitation Phys. Rev. B 78 224304
Horsfield A P, Bowler D R, Fisher A J, Todorov T N and Montgomery M J 2004 Power dissipation in nanoscale conductors: classical, semi-classical and quantum dynamics J. Phys.: Condens. Matter 16 3609-22
Horsfield A P, Bowler D R, Fisher A J, Todorov T N and Sanchez C 2004 Beyond ehrenfest: correlated non-adiabatic molecular dynamics J. Phys.: Condens. Matter 16 8251-66 Horsfield A P, Bowler D R, Fisher A J, Todorov T N and Sanchez C G 2005 Correlated electron-ion dynamics: the excitation of atomic motion by energetic electrons J. Phys.: Condens. Matter 17 4793-812 Horsfield A P, Bowler D R, Ness H, Sanchez C G,
Todorov T N and Fisher A J 2006 The transfer of energy between electrons and ions in solids Rep. Prog. Phys. 69 1195-234
[156] le Page J, Mason D R and Foulkes W M C 2008 The
ehrenfest approximation for electrons coupled to a phonon system J. Phys.: Condens. Matter 20 1 25212
[157] le Page J 2008 The transfer of energy between electrons and
ions in solids PhD Thesis University of London
[158] Stella L, Meister M, Fisher A J and Horsfield A P 2007
Robust non-adiabatic molecular dynamics for metals and insulator J. Phys. Chem. 127 214104
[159] McEniry E J, Bowler D R, Dundas D, Horsfield A P,
Sanchez C G and Todorov T N 2007 Dynamical simulation of inelastic quantum transport J. Phys.: Condens. Matter 19 196201
[160] McEniry E J, Frederiksen T, Todorov T N, Dundas D and
Horsfield A P 2008 Inelastic quantum transport in nanostructures: the self-consistent born approximation and correlated electron-ion dynamics Phys. Rev. B 78 35446
