and other phenomena
and other phenomena
Proefschrift
ter verkrijging van de graad van doctor
aan de Technische Universiteit Delft,
op gezag van de Rector Magnificus prof. dr. ir. J.T. Fokkema,
voorzitter van het College voor Promoties,
in het openbaar te verdedigen op 5 oktober 2004 om 10.30 uur,
door
Theodorus Petrus Cornelis KLAVER
materiaalkundig ingenieur,
geboren te Den Haag
Samenstelling promotiecommissie:
Rector Magnificus voorzitter
Prof. dr. B. J. Thijsse Technische Universiteit Delft, promotor Prof. dr. S. W. de Leeuw Technische Universiteit Delft
Prof. dr. G. J. Kramer Technische Universiteit Eindhoven Prof. E. van der Giessen Rijksuniversiteit Groningen
Prof. dr. K. Albe Technische Universität Darmstadt
Prof. dr. J. Fransaer Katholieke Universiteit Leuven
Dr. ir. J. Sietsma Technische Universiteit Delft
Published and distributed by: DUP Science DUP Science is an imprint of
Delft University Press P. O. Box 98 2600MG Delft The Netherlands Telephone: +31 15 2785678 Telefax: +31 15 2785706 E-mail: info@library.tudelft.nl ISBN 90-407-2526-8
Keywords: simulation, film, deposition Copyright © 2004 by Peter Klaver
All rights reserved. No part of the material protected by this copyright notice may be reproduced or utilized in any form or by any means, electronic or mechanical, including photocopying , recording or by any information storage and retrieval system, without written permission from the publisher: Delft University Press
Printed in The Netherlands
This work is part of the research program of the Stichting voor Fundamenteel Onderzoek der Materie (Foundation for Fundamental Research of Matter), and was made possible by financial support from the Nederlandse Organisatie voor Wetenschappelijk Onderzoek NWO (Netherlands Organization for Scientific Research).
1. Introduction 1
1.1 Scope of this thesis 2
1.2 Computational details 4
2. Thin Ta films: growth, stability, and diffusion studied by Molecular 11
Dynamics simulations
2.1 Introduction 12
2.2 Computational details 12
2.3 Results 13
2.3.1 The β-Ta structure 13
2.3.2 Thermal PVD film growth 15
2.3.3 Ion bombardment 16
2.3.4 Annealing 17
2.3.5 Vacancy and adatom diffusion 19
2.3.6 DFT verification of EAM results 22
2.4 Conclusions 22
References 23
3. Molecular Dynamics study of Cu thin film deposition on β-Ta 25
3.1 Introduction 26
3.2 Computational details 26
3.3 Results 27
3.3.1 Cu deposition on a β-Ta (001) hexagonal surface 27
3.3.2 Cu deposition on a β-Ta (001) intermediate surface 28
3.3.3 Cu deposition on a β-Ta (001) grown surface 29
3.4 Conclusions 30
References 30
4. Molecular Dynamics simulations of Cu/Ta and Ta/Cu thin film growth 31
4.1 Introduction 32
4.2 Simulation details 32
4.2.1 Computational details 32
4.2.2 Simulation temperature 33
4.2.3 DFT verification of the empirical Ta-Cu potential 34
4.3 Results 35
4.3.1 Cu deposition on bcc Ta substrates 35
4.3.1.1 Cu deposition an a flat bcc (100) substrate 35
4.3.1.2 Cu deposition an a flat bcc (110) substrate 36
4.3.1.3 Cu deposition on a rough bcc (110) substrate 37
4.3.2 Comparison between Cu deposition on bcc Ta and β-Ta substrates 38
4.3.3 Ta deposition on Cu (111) substrates 39
4.4 Conclusions 41
Appendix 42
References 42
5. The interaction of N with atomically dispersed Ti, V, Cr, Mo, and Ni in 45
ferritic steel
5.1 Introduction 46
5.2 Calculations and results 47
5.3 Discussion 50
5.4 Conclusions 51
6.3 Results 56
6.4 Discussion 58
6.5 Conclusions 59
References 59
7. Constructing a Density Functional Theory-based empirical Si potential 61
7.1 Introduction 62
7.2 Comparing empirical potentials to Density Functional Theory data 62
7.3 Towards a new MEAM silicon potential 63
References 66
Appendix A. Cu and Ta properties 68
Summary 69
Samenvatting 72
Dankwoord 77
Curriculum Vitae 79
1.1 Scope of this thesis
The rate of progress in the computer industry shows no signs of slowing down. Almost three decades after it was first expressed in 1975, Moore’s law is still valid. Keeping up with Moore’s law requires constant improvements and innovations in chip design and manufacturing. One way to improve chip performance is to reduce the feature sizes. This means that the current carrying interconnects in chips will get thinner and will be required to handle higher and higher current densities. The typical minimum feature size in current state-of-the-art chips is 90 nm and the corresponding maximum current density is already more than 106 A/cm2. The feature size is
expected to decrease to 45 nm by 2010 and the current density is expected to approximately double accordingly [1, 2, 3].
When current densities exceed a certain value the currents cause electromigration, i. e. the momentum transfer from the ‘electronic wind’ causes the atoms in the interconnect to be displaced. This unwanted material transport can lead to device failure through various
mechanisms, e. g. void formation or stress buildup [4]. Aluminium (Al) has been the standard interconnect metal for decades, up to the end of the 1990s. It has good processing properties and a reasonably low electrical resistance (2.7 µΩ-cm). However, copper (Cu) has still better electrical properties1
. It has a lower electrical resistivity (1.7 µΩ-cm) and a resistivity to electromigration which is more than an order of magnitude higher than for Al [5, 6, 7, 8, 9, 10, 11]. The reason why Cu has not always been used in the first place is the high rate at which it diffuses into silicon (Si) and silicon oxide (SiO2) [12, 13, 14, 15, 16]. If Cu diffuses into Si it causes deep level traps which
are detrimental to the Si semiconducting properties [17, 18, 19]. The relentless pace of progress in the computer industry requires the transition from Al to Cu interconnects, but this transition can only be made if the Cu can be kept from diffusing into Si. For this reason diffusion barriers are required which must at the very least possess two somewhat contradictory properties: they must be reasonably chemically inert to Cu, Si and SiO2, yet they must show good adhesion to these
materials. Tantalum (Ta) is one of a number of materials which meet these requirements [19, 20, 21, 22] and indeed it has already been applied as a diffusion barrier material in chips with Cu interconnects. However, the level of understanding and control over the use of Cu interconnects is not yet at the level of Al interconnects. More research is needed in order to push Cu interconnect technology to the same optimal level as has been done for Al.
Even if Cu interconnect technology is pushed to its extreme, there will be limits to what can be gained by doing so. As interconnects get thinner their electrical resistance increases due to surface scattering. So even though the resistance to electromigration of Cu is more than an order of magnitude higher than for Al, this does not mean that interconnect thicknesses can be reduced proportionally. Even if surface scattering would not a problem there would be other limitations. For instance, if current interconnects would be reduced in thickness by a few dozens times, the resulting interconnects would be no more than ~1 nm thick. This means that quantum confinement effects start to play a role. Therefore it is expected that by 2010 a more drastic solution will have to be found, like replacing metal interconnects by nanotubes or dispense with electrical signals altogether and switch to optical propagation or microwave interconnects [1, 2, 3]. Still, Cu interconnects should help to keep up the pace of progress in the computer industry to the end of this decade.
The extent in which electromigration manifests itself in Cu depends on the Cu microstructure. This thesis covers the microstructural evolution of evaporated thin Cu/Ta and Ta/Cu films. More specifically, phenomena like growth mode, surface roughness evolution, diffusion, nucleation, grain boundaries, grain growth, creation of stacking faults and dislocations, interface mixing, epitaxy, phase stability, and texture development are studied. The paragraphs above would seem to indicate a very practical, even urgent relevance of this research. They are indeed effective
propaganda to get research proposals approved in which ‘applicability’ or ‘social relevance’ is one of the criteria for approval. But to quote Richard Feynman: ’Science is like sex: sometimes
1
something useful comes out, but that is not the reason we are doing it’. The real reason for doing this research is that the Ta/Cu system offers various angles to investigate interesting phenomena, like the metastable β-Ta crystal structure with its inequivalent lattice positions, hetero-epitaxy between two elements with very different atomic sizes, crystal structures and melting points, and non-equilibrium intermixing in a system of thermodynamically immiscible elements.
Atomistic simulations were chosen as the method of investigation. Atomistic simulations offer the possibility to perform a detailed atomic-scale analysis of the results. They also offer the possibility to control and manipulate atoms, e. g. an atom in an energy minimum can be given increasing amounts of kinetic energy to see how much is required to overcome an energy barrier and move to another minimum. Through this control a wide range of experiments can be simulated quite easily. There are many different kinds of atomistic simulations, each with its strengths and weaknesses. The results in this thesis are obtained from an empirical force field method (the Embedded Atom Method, EAM) and from Density Functional Theory (DFT).
The EAM is a fast N-body calculation method which is well suited for large metallic systems. The computational power of a handful of todays modern desktop computers is sufficient to calculate the time evolution of ~105 atoms for ~10-8 s. The great disadvantage of EAM potentials
(and other empirical potentials) is that they require parameter fitting to relevant data for each possible pair of interacting elements before they can be used. This means that any simulation results that involve situations that are far removed from the data in the database should be considered interpolations or extrapolations. This of course includes any newly discovered phenomena. Also, the fitting of the empirical parameters can be a very time-consuming task. Moreover, if no relevant experimental data is available, it would first require the calculation of an ab-initio data set. Ab-initio calculations like DFT are much more widely applicable (without sacrificing accuracy) than empirical calculations. However, DFT calculations are orders of magnitude more demanding in terms of computer resources. Computers that can handle EAM calculations of ~105 atoms for ~10-8 s can only handle DFT calculations of ~100 atoms for ~10-11 s
in the same time. In order to combine the strong points of both methods, large-scale calculations in this thesis (film deposition, ion bombardment, annealing deposited films) are carried out with EAM potentials, and small Cu/Ta test systems are calculated both with EAM and DFT potentials. The calculation of small systems with both methods allows one to determine how much EAM results deviate from the more accurate values and to estimate if any of the large-scale EAM results are unreliable. In chapter 1.2 the computational methods are presented.
In chapter 2 the results of Ta deposition (deposition energy < 1 eV/atom) on bcc Ta and β-Ta substrates are reported. The β-Ta films are annealed and bombarded with 250 eV Ar ions to investigate the stability of the β-Ta crystal structure. Both crystal structures are simulated with an EAM potential that has been fitted to bcc crystal properties. DFT calculations show that several properties of the β-Ta crystal structure are accurately reproduced by this bcc-fitted EAM potential. Apart from being useed in large scale calculations, the EAM potentials are also used to study diffusion of single adatoms on atomically flat surfaces and vacancies, including the somewhat complex vacancy diffusion in β-Ta.
The deposited β-Ta films can be used for subsequent Cu deposition. In chapter 3 the results of Cu deposition on these films as well as on perfectly flat β-Ta substrates are reported. Cu deposition on bcc Ta substrates has also been simulated and the results of these simulations are reported in chapter 4. Part of the Cu/Ta work consists of DFT calculations to verify the accuracy of the Cu-Ta potential. Implications of the outcome of the DFT calculations for the EAM results of chapters 3 and 4 are discussed. Also, a comparison between the resulting microstructures of Cu films on bcc Ta and β-Ta substrates is made. Finally, chapter 4 briefly discusses Ta deposition on Cu
substrates.
DFT calculations are a widely applicable technique. Three DFT investigations unrelated to Cu/Ta are also reported in this thesis.
The first DFT study deals with the interaction of nitrogen (N) interstitials with atomically dispersed alloy atoms in steels.
Hardness experiments can give an excellent indication of the presence of N in steels. N concentration profiles can be also calculated. There is some disagreement between an often-used model by Sun and Bell and experimental hardness data. The model assumes that N is either completely immobilised in nitrides or can diffuse freely through the lattice. Trapping of N by single alloy atoms is not accounted for. If such trapping takes place it would lower the N diffusion coefficient and this would lead to better agreement between the calculated concentration profiles and experimental data. Therefore the question is: is there any trapping energy which ties N interstitials to single alloy atoms and if so, how large is this energy? The DFT calculations that help to answer this question for Ti, V, Cr, Ni, and Mo are reported in chapter 5.
Chapter 6 reports the results of DFT calculations of aluminium (Al) alloys with interstitials. The question underlying this investigation is somewhat similar to that concerning N trapping in steels. In this case the energy by which alloy atoms are tied into interstitials is investigated. It is
reasonable to assume that it is energetically favourable to have small atoms in interstitial positions. Having a higher than average concentration of alloy atoms in interstitial positions can lead to concentration fluctuations and even segregation. This can be understood as follows: interstitials in metals are generally more mobile than vacancies. If alloy atoms are strongly tied into interstitials they can travel with the interstitials over long distances. Eventually the interstitials will end up in interstitial sink sites or link up to form platelets in the bulk. At these sites an alloy enrichment will appear, or segregation may even occur. Other parts of the alloy will be depleted of alloy element atoms. For interstitial-mediated segregation to occur it is essential that the alloy atoms are tied into the interstitial with sufficient energy. The energies by which alloy atoms are tied into interstitials in Al have been calculated for Mg, Si, Cr, Mn, Fe, Co, Ni, and Cu.
DFT calculations can be used to verify empirical potentials. Even better than verifying afterwards is to create potentials based on relevant DFT data before running the calculations with empirical potentials. Chapter 7 describes the first results of an effort to do this for Si. The reason for doing this is that four existing empirical Si potentials (Modified Embedded Atom Method, Stillinger-Weber, Tersoff-III, and a potential by Lenosky et al) all have serious shortcomings when applied to situations that deviate strongly from the equilibrium lattice.
DFT calculations were carried out to determine the cohesion energies of Si in various crystal structures, the unrelaxed vacancy formation energies for these crystal structures, and the energies of several small clusters of atoms. These energies are used as fitting data in an effort to construct a Modified Embedded Atom Method Si potential.
1.2 Computational details 1.2.1 Embedded Atom Method
For large-scale calculations the EAM is used. In the EAM the total energy U of a system of N atoms is € U = Fi i=1 N
∑
(ρi(ri)) +1 2 j=1, j≠iφij(rij) N∑
i=1 N∑
(1)in which Fi is the embedding function, i.e. the function that describes how the electron density
part of the potential energy of atom i at position ri depends on the electron density ρi at that
position, and φij is the pair potential contribution of two atoms i and j separated by a distance rij.
The electron density experienced by an atom only depends on the distances to other atoms. The value at a certain position is a lineair superposition of the densities contributed by different atoms,
€
ρi = fj(rij) j=1: j≠i
N
∑
, (2)where fj is the electron density contributed by atom j at the position of atom i.
1.2.2 Ta-Ta interaction
Johnson and Oh have introduced a functional form for EAM potentials for bcc metals [23]. The embedding function, electron density around each atom, and the pair potential are given by
€ F(ρ) = − EC − E1V UF
(
)
1 − ln ρ ρe n ρ ρe n , (3) € f (r) = fe r1e r β , (4) € φ(r) = K3 r r1e −1 3 + K2 r r1e −1 2 + K1 r r1e −1 + K0, (5) in which € n = 9ΩB −15ΩG β2 EC − E1V UF(
)
, (6) € K3 = − 15ΩG 3A + 2 3 4A − 1 2S , (7) € K2 = −15ΩG 3A + 2 15 16A − 3 4S − 9 8AS + 7 8 , (8) € K1 = −15ΩG 7 8− 3 4S , (9) € K0 = −15ΩG 3A + 2 201 112A − 27 28S − 87 56AS + 187 168 −1 7E1V UF, (10) in which € S =r2e r1e, (11) € A = 2C44 C11− C12 , (12) € G =3C44+ C11− C12 5 , (13) € B =C11+ 2C12 3 . (14)In these equations EC is the equilibrium cohesive energy per atom,
€
E1VUF is the unrelaxed monovacancy formation energy, ρe is the (dimensionless) equilibrium electron density, fe is a
dimensionless element-specific weight factor that cancels out for single-element systems, r1e and
r2e are the first and second equilibrium neighbour distances, β is a constant (of little influence, a
value of 6 was chosen by Johnson and Oh), Ω is the atomic volume, B is the bulk modulus, G is the Voigt average shear modulus, A is the anisotropy factor, and C11, C12, and C44 are the cubic
elastic constants.
To implement eqns. (1)-(14), a few small alterations have been made. Firstly, the radial interaction range around each atom has been limited to the first and second neighbour atoms in the equilibrium Ta bcc crystal only. This was achieved by replacing the electron density function in the interval between
rs = r2e + 0.1(r2e - r3e) and rc = r2e + 0.5(r2e - r3e) by a smooth cut-off function,
€ fc(r) =−2 f (rs) + (rs− rc) f '(rs) (rs− rc) 3 (r − rc) 3 +3 f (rs) − (rs − rc) f '(rs) (rs− rc) 2 (r − rc) 2. (15)
This function has the properties fc(rs) = f(rs), fc’(rs) = f’(rs), fc(rc) = fc’(rc) = 0 . Therefore the
corrected function smoothly connects to the original function at rs and then gradually goes to zero
at rc (3.98 Å). The pair potential has been adapted in a similar manner.
Secondly, the embedding function of eqn. (3) is very steep near ρ = 0, the derivative tends to -∞. As a result of discretization errors this singularity can produce unreliable results. To alleviate this problem, the embedding function has been multiplied by the function
€ h(ρ) = 1 − e −1 2 ρ ρs 2 , (16)
in which ρs was given the value 0.5. This fixes the problem with the singularity. There is no
theoretical basis for this, but only the low-density part of the embedding function is significantly affected. The embedding energy for low densities is small, consequently modifying it by a small fraction has no significant impact.
Table I shows the input values that were used in the creation of the Ta-Ta EAM potential. Since the Johnson-Oh EAM was derived for low energy situations, different potentials are used for short distances. For distances smaller than r1e (2.86 Å) the Johnson-Oh pair
potential is stiffened. Instead of φ(r) Johnson and Oh have proposed to use φS(r), given by
€ φS
( )
r =φ( )
r + 4.5 1 + 4 A − 0.1 (
φ( )
r −φ( )
r1e)
r r1e −1 2 (17)For even smaller distances (< 1.8 Å) the Screened Coulomb pair potential with Molière weight factors and Firsov
screening length [24] is used. This purely repulsive potential has the form
€ φ(r) = Z1Z2e 2 4πε0r γ r a , (18)
Table I. Values used for constructing the EAM potential for bcc Ta. The values are the lattice parameter a, atomic volume Ω, the bulk modulus B, the Voigt average shearmodulus G, the anisotropy ratio A, the cohesive energy EC, and the unrelaxed
monovacancy formation energy
€ E1VUF. quantity a (Å) 3.3026 ΩB (eV) 21.66 ΩG (eV) 7.91 A 1.57 EC (eV) 8.089 E1VUF (eV) 2.95
in which € γ r a = cie−di r a 1 3
∑
, ci =γ(0) = 1 1 3∑
, (19) and € a = 9π 2 128 1 3 aB Z1 1 2 + Z 2 1 2 −2 3 , (20)in which Z1 and Z2 are the atom numbers of the interacting atoms at distance r, e is the elementary
charge, ε0 is the vacuum permittivity, γ is the screening function, a is the screening length, ci and di
are constants, and aB is the Bohr radius. The Molière values for ci and di are given in table II.
Since the EAM potential is based on a combination of a pair potential and an electron density and the Screened Coulomb potential is a pure pair potential, a method has to be devised to smoothly go from one potential to the other. One way to do this is to let the electron density function fall off to zero when the distance from the nucleus decreases towards the range where the Screened Coulomb pair potential takes effect. This is done by multiplying f(r) with a Fermi-Dirac-like function,
€ g(r) = 1 e rz−r Δrz + 1 , (21)
in which rz is a constant that determines where the function is equal to 1/2 and Δrz is a constant
that determines how steeply the function rises from 0 to 1. For rz and Δrz the values of 1.8 Å and
0.025 Å were chosen. In this way the electron density and hence the embedding energy are practically zero near the atom core. The stiffened Johnson-Oh and Screened Coulomb pair potentials have been smoothly spline-connected in a small r-range around 1.8 Å, where they intersect and their slopes are almost the same.
1.2.3 Ta-Ar interaction
For the Ta-Ar interaction the Firsov-Molière Screened Coulomb pair potential has been used with the values in table II. No attractive potential was added to this purely repulsive potential. This is justified because Ar atoms are practically always in high-energy situations where neglecting a few meV of binding energy has no influence.
1.2.4 Cu-Cu interaction
Johnson and Oh have introduced the following functional form for EAM potentials for fcc metals [25,26]: € F(ρ) = a ρ ρe n + b ρ ρe (22) € f (r) = fold
( )
r − fc( )
r (23) € φ(r) =φold( )
r −φc( )
r (24)Table II: Molière constants in the screening function.
i ci di
1 0.35 0.3 2 0.55 1.2 3 0.10 6.0
where € fold(r) = fee −β r r 1e −1 (25) € φold(r) = φee −γ r r1e−1 (26) € fc
( )
r = fold( )
rc + g r( )
f 'old( )
rc / g' r( )
c (27) € φc( )
r =φold( )
rc + g r( )
φ'old( )
rc / g' r( )
c (28) € g r( )
= 1 − e δ r r1e− rc r1e (29)Eqns. (23)-(29) already include a smooth function cut-off, no adaptation as in the case of the Ta-Ta interaction is necessary. The embedding function in eqn. (22) has no singularity at ρ = 0, so multiplying it with eqn. (16) is also not necessary. Finally, Cu films were never bombarded with Ar ions, so it was not necessary to stiffen the pair potential or replace the interaction in the core region by a Screened Coulomb pair potential.
Table III lists the parameters used in the construction of the Cu-Cu potential. These values were either generally sensible choices made by Johnson et al (for those parameters which do not have too great an influence), or were fitted to cohesion energy, vacancy formation energy, atomic volume, and Voigt average bulk and shear moduli data. The cut-off radius includes the first, second, and third-nearest
neighbours in the equilibrium Cu fcc lattice.
The atom size appears in the potential not only through the atomic volume to which a number of parameters are fitted, it is also present through r1e (2.55 Å) which in turn influences the equilibrium electron
density ρe. Cu has a rather large thermal expansion coefficient at room
temperature (16.5*10-6/K) which would give rise to thermal stresses if Cu is simulated at high
temperatures but at room temperature lattice spacing (which is often the case in this thesis). Therefore different Cu potentials with different values for r1e were constructed for high
temperatures. The 0 K value for r1e was chosen to be slightly smaller, such that at the higher
simulation temperatures the Cu lattice has the experimental room temperature lattice constant. It was verified that after making small changes to r1e the potential still accurately reproduces the
other properties to which the potential was fitted.
The value for fe was determined to reproduce the heat of formation of a 50-50% ordered CuTa
alloy, see next section.
1.2.5 Cu-Ta interaction
Johnson has presented a model, based on the heat of solution, for a pair potential to describe the interaction between unlike elements in binary alloys [27]. In its original form, the pair potential between atoms of elements a and b is given by
€ φab(r) = 1 2 fb
( )
r fa( )
r φaa( )
r + fa( )
r fb( )
r φbb( )
r . (30)Table III. Parameter values used in the construction of the Cu-Cu potential.
parameter β 5.0 γ 8.5 δ 20.0 rc (Å) 4.85 φe (eV) 0.36952 a -4.0956 b -1.6979 n 0.44217
To avoid the singularities that arise when the different elements have different radial cut-off radii (4.85 Å for Cu, 3.99 Å for Ta), we slightly adapted Johnsons functional form for Cu-Ta to
€ φCuTa
( )
r =1 2 fTa2( )
r φCu( )
r + fCu 2 r( )
φTa( )
r 1 2(
fCu( )
r + fTa( )
r)
2 . (31)The functional form of eq. (31) has only one free parameter, i.e. the ratio of the weight factors fe,Cu
and fe,Ta. This parameter was fitted to the heat of formation of a fully relaxed, CsCl-type, ordered
50-50 atom% Ta-Cu alloy. The (positive) heat of formation of such an alloy is 0.03 eV/atom according to Miedema [28]. A value of 0.8924 for fe,Cu/fe,Ta (fe,Cu = 0.8924, fe,Ta =1) results in this
heat of formation.
1.2.6 CAMELION
The potentials described above are used with the CAMELION code, developed at Delft
University. CAMELION is a fast, efficient Fortran Molecular Dynamics code that runs on both single-cpu and parallel computer systems. CAMELION calculates energies and forces for a system of atoms described by EAM potentials. Energies and forces due to the pair potentials and electron densities are linearly interpolated from values tabulated at equal intervals Δr2 = 0.036 Å2.
Embedding energies and forces are linearly interpolated from values tabulated at Δρ = 0.01. The Velocity-Verlet algorithm [29] is used to advance the atoms with time. The time-step Δt is variable and is chosen such that in each timestep no atom will move more than a specified distance Δrmax:
€ Δrmax
(
)
2 < v2( )
Δt 2 + (a • v) Δt( )
3+ 1 4a 2 Δt( )
4, (32)where v is the velocity and a is acceleration. The distance estimate includes the acceleration because in certain situations this part far outweighs the velocity part. In the present simulations, Δrmax = 0.02 Å was chosen. Periodic boundary conditions are used to simulate ‘large’ systems.
Temperature is controlled through a Berendsen thermostat [30], typically with a 0.018 ps time constant. Atoms can be introduced into the system under any angle with a selectable kinetic energy and are allowed to leave the system (sputtering). Any number of atoms can be fixed to their equilibrium positions by a harmonic force. This prevents systems from drifting away by the momentum of incoming atoms, yet it interferes less with the natural evolution of the system than fixed atoms. In most of the simulations in this thesis the atoms in the bottom monolayer are held fixed with a force constant of 19.4 N/m. Configurations are stored in text files. Although not used for work in this thesis, CAMELION can also use Modified Embedded Atom Method [31]
potentials, do Dynamic Monte Carlo calculations [32], and use a Berendsen barostat for pressure control [30].
1.2.7 VASP
For small systems DFT calculations were done with the Vienna Ab-initio Simulation Package (VASP) [33, 34]. VASP is a plane-wave code that uses either Vanderbilt [35, 36] type ultra-soft pseudopotentials (USPP) or the Projector Augmented Wave (PAW) method [37] in the Local Density Approximation (LDA) or Generalised Gradient Approximation (GGA). VASP is a fast DFT code which has many built-in options for automated structural relaxation and ab-initio Molecular Dynamics. The Fortran90 source code is available to users and can be compiled on many platforms, both for single-cpu and parallel execution.
References
[1] International Technology Roadmap for Semiconductors, 2001 edition, Interconnect, http://public.itrs.net/
[2] International Technology Roadmap for Semiconductors, 2002 update, Interconnect, http://public.itrs.net/
[3] International Technology Roadmap for Semiconductors, 2003 Edition, Executive Summary, http://public.itrs.net/
[4] J. A. Nucci, A. Straub, E. Bischoff, E. Arzt, and C. A. Volkert, J. Mater. Res. 17 (2002) 2727
[5] R. Rosenberg, D. C. Edelstein, C.-K. Hu, and K. P. Rodbell, Annu. Rev. Mater. Sci. 30 (2000) 229
[6] D. Wiess, O. Kraft, and E. Arzt, J. Mater. Res. 17 (2002) 1363
[7] C.-L. Lin, P.-S. Chen, and M.-C. Chen, J. Vac. Sci. Technol. B 20 (2002) 1111 [8] C.H. Lin, J.P. Chu, T. Mahalingam, and T.N. Lin, J. Mater. Res. 18 (2003) 1429 [9] L. Chen, N. Magtoto, B. Ekstrom, and J. Kelber, Thin Solid Films 376 (2000) 115 [10] H. Ono, T. Nakano, and T. Ohta, Appl. Phys. Lett. 64 (1994) 1511
[11] K.-M. Yin et al, Thin Solid Films 388 (2001) 15
[12] R. P. Vinci, S. A. Forrest, and J. C. Bravman, J. Mater. Res. 17 (2002) 1863 [13] K.-H. Min, K.-C. Chun, and K.-B. Kim, J. Vac. Sci. Technol. B 14 (1996) 3263
[14] K. Radhakrishnan, Ng Geok Ing, and R. Gopalakrishnan, Mat. Sci. Eng. B 57 (1999) 224 [15] E. Kolawa, J. S. Chen, J. S. Reid, P. J. Pokela, and M.-A. Nicolet, J. Appl. Phys. 70 (1991)
1369
[16] L. C. Lane, T. C. Nason, G.-R. Yang, T.-M. Lu, and H. Bakhru, J. App. Phys 69 (1991) 6719
[17] A. E. Kaloyeros and E. Eisenbraun, Annu. Rev. Mater. Sci. 30 (2000) 363
[18] A. E. Kaloyeros, X. Chen, S. Lane, and H. L. Frisch, J. Mater. Res. 15 (2000) 2800 [19] K. Holloway, P. M. Fryer, C. Cabral, J. M. E. Harper, P. J. Bailey, and K. H. Kelleher,
J. Appl. Phys. 71 (1992) 5433
[20] D. Fischer, O. Meissner, B. Bendjus, J. Schreiber, M. Stavrev, and C. Wenzel, Surf. Int. An. 25 (1997) 522
[21] M. Stavrev, D. Fischer, C. Wenzel, K. Drescher, and N. Mattern, Thin Solid Films, 307 (1997) 79
[22] L. A. Clevenger, A. Mutscheller, J. M. E. Harper, C. Cabral, and K. Barmak, J. App. Phys 72 (1992) 4918
[23] R.A. Johnson, D.J. Oh, J. Mater. Res. 4 (1989) 1195
[24] W. Eckstein, Computer Simulation of Ion-Solid Interactions, Springer, Heidelberg, 1991, p.40.
[25] R.A. Johnson, D.J. Oh, J. Mater. Res. 3 no. 3 (1988) 471.
[26] R.A. Johnson, D.J. Oh, Embedded atom method for close-packed metals, in ‘Atomistic Simulations of Materials: Beyond Pair Potentials’ eds. Vaclav Vitec, David J. Srolovitz, Jan 1989, p233-238.
[27] R. A. Johnson, Phys. Rev. B 41 (1990) 9717
[28] F. R. de Boer, R. Boom, W. C. M. Mattens, A. R. Miedema, and A. K. Niessen, Cohesion in metals, ed. F. R. de Boer and D. Pettifor (North Holland Physics Publishing 1989) p. 544 [29] W. C. Swope, H. C. Andersen, P. H. Berens, and K. R. Wilson, J. Chem. Phys. 81 (1984)
3684
[30] H. J. C. Berendsen et al, J. Chem. Phys. 81 (1984) 3684 [31] M. I. Baskes, Phys. Rev. B 46 (1992) 2727
[32] G. Dereli, Molec. Simul. 8 (1992) 351
[33] G. Kresse, J. Hafner, Phys. Rev. B 47 (1993) RC558 [34] G. Kresse, J. Furthmuller, Phys.Rev . B 54 (1996) 11169 [35] D. Vanderbilt, Phys. Rev. B 41 (1990) 7892
[36] G. Kresse, J. Hafner, J. Phys. Condens. Mat. 6 (1994) 8245 [37] G. Kresse en D. Joubert, Phys. Rev. B 59 (1999) 1758.
Thin Ta films: growth, stability, and diffusion,
studied by Molecular Dynamics simulations
Abstract
We report results of classical Molecular Dynamics simulations of bcc and β-Ta thin films. Thermal PVD film growth, surface roughness, argon ion bombardment, phase stability and transformation, vacancy and adatom diffusion, and thermal relaxation kinetics are discussed. Distinct differences between the two structures are observed, including a complex vacancy diffusion mechanism in β-Ta. Embedded Atom Method potentials, which were fitted to bcc properties, have been used to model the Ta-Ta interactions. In order to verify the application of these potentials to the more complex β-Ta structure, we have also performed Density Functional Theory calculations. Results and implications of these calculations are discussed.
2.1 Introduction
The unabating need to produce faster computers is forcing chipmakers to replace the traditional aluminium interconnects with copper ones. The electrical properties of copper are superior to those of aluminium (higher electrical conductivity, much higher resistance to electromigration). However, copper diffuses into silicon and silicon oxide quite rapidly, degrading the device properties. Therefore diffusion barriers must be used to keep copper from interacting with silicon and silicon oxide. Tantalum is one of a number of metals and compounds under consideration for their good barrier properties (a small diffusion coefficient for copper and the somewhat
contradictory combination of good adhesion and low chemical affinity). Understanding the growth and the stability of a Ta thin film is therefore directly related to its technological importance. In addition, Ta is interesting from a more fundamental point of view because of the metastable β-phase in which it sometimes appears [1] in thin films. This is not only an academic issue. The difference between films grown in the β-phase and those in the bcc equilibrium phase has an influence on the Cu seed layer that is usually deposited on top of the Ta layer in the subsequent IC-processing step. Microstructure and surface roughness of the Ta layer will affect adhesion, epitaxy, diffusion, and interface conductivity.
The purpose of this work is to obtain atomic-level information on the evolution of Ta films. More specifically, we study the growth of β and bcc Ta films, and we investigate diffusion and stability under Ar bombardment and thermal annealing, using classical molecular-dynamics simulations. The conditions under which Ta films are deposited in IC manufacturing are not exactly copied in this work. Rather, we have chosen to investigate the fundamental atomic
processes by clearly distinguishing between deposition, ion bombardment, and thermal annealing. As far as we know, β-Ta has not been studied by molecular-dynamics methods before.
Part of this paper covers Density Functional Theory results which we have calculated and used, among other things, as a means to do a computational verification of the β-Ta structure and the classical molecular dynamics results.
2.2 Computational details
The classical simulations were performed using the CAMELION Molecular Dynamics code developed in Delft [2], which employs Embedded Atom Method (EAM) potentials in the Johnson-Oh form [3] to model the Ta-Ta interactions. The Ta-Ar and Ar-Ar interactions were modelled by the Firsov-Molière Screened Coulomb potential [4]. Details are described in [2].
The Johnson-Oh EAM potential for Ta has been fitted to properties of the perfect bcc lattice (cohesion energy, lattice constant) and to properties that represent small distortions of the lattice (vacancy formation energy, elastic constants). Applying this potential to situations which deviate strongly from the perfect bcc lattice, like the metastable β-Ta phase, should therefore be
considered a large extrapolation. It is not obvious that these situations are well described by the same EAM parameter set. In order to verify the validity of the EAM results, we have performed ab-initio calculations using the program VASP (Vienna Ab-initio Simulation Package), developed at the Institut für Materialphysik of the Universität Wien [5,6]. VASP is a plane-wave Density Functional Theory program that uses ultra-soft Vanderbilt type pseudopotentials [7,8]. More details are given in sec. 3.6.
Several conditions in the classical simulations have no other particular significance than that they are the same as those used in thermal helium desorption experiments that are carried out in our group to study thermal and ion-beam assisted Ta film deposition: the slightly off-normal incidence of the Ta vapour flux (15°) and the normal incidence of the Ar ions. The direction of the incoming Ta atoms has an in-plane component making a 17° angle with the <100>-axis, which is intended as an “arbitrary” angle without a specific relation to the crystal structure. The Ta atoms arrive at the surface with a thermal energy corresponding to 2000 K, the Ar ion energy is 250 eV.
The deposition rate in the simulations (0.294 m/s) is ten orders of magnitude higher than in the experiments, in order to finish simulations in a feasible amount of time. In order to permit the
system nevertheless a realistic time period for thermal relaxation in the simulations, we have adopted the simple method of giving the atoms in the film and substrate a higher mean kinetic energy than in the experiments. Running the simulation at a higher temperature thus compensates to first approximation for the high deposition rate. This approach is motivated by the Arrhenius model
N = ν0te
−Q / kT, (1)
where N denotes the expected number of activated atomic processes of activation energy Q that occur in time t at temperature T , and ν0 is the effective attempt frequency, which we will show
later to be approximately 5 × 1013 s−1. Eq. (1) shows that in a typical room temperature experiment
taking 10 s, processes up to 0.84 eV can be considered as being thermally activated, because these will occur 3 or more times (the number 3 is somewhat arbitrarily chosen). Holding the film at 1000 K during a 1 ns simulation has the same effect, which is why 1000 K was chosen as the standard simulation temperature. Note that this reasoning does not give the full energetic picture: as each atom condensates on the film, an additional 8 eV condensation energy becomes locally available.
In order to study phase transformation in detail, each atom is given a crystal symmetry index according to the shape of its coordination polyhedron. For example, an atom is marked as bcc-like, when its coordination polyhedron has a closer resemblance to that of an atom in a perfect bcc crystal than to an atom in any of the alternative structures considered (fcc, hcp, icosahedral, diamond). Additionally, a certain minimum resemblance criterion must be met, or else the atom’s local symmetry is marked “undefined”. The method of comparison is based on the spherical harmonic rotational invariants of ref. [9]. Configurations are quickly cooled to 0 K before applying the rotational invariants algorithm, in order to exclude thermal vibrations.
Where applicable, the surface roughness of a film is defined as the standard deviation of the height-coordinates of the surface atoms. Surface atoms are those atoms which lack five or more nearest or next nearest neighbours.
2.3 Results
2.3.1 The β-Ta structure
In the literature there is some uncertainty about the metastable tetragonal β-Ta structure, which is reported to be pseudo-hexagonal, like β-uranium, in a somewhat older paper by Mosely and Seabrook [10] (see fig. 1). Some impurity stabilized Ta structures have been reported [11,12], which are also referred to as β-Ta by some, but these will not be considered in this paper. We have not found any other proposed structure for β-Ta than that by Mosely and Seabrook, nor a (more recent) confirmation. It has been entered into the JCPDS XRD-data index and is therefore assumed to be correct by most people. However, Donohue and Einspahr [13] have cast some doubt on the exact nature of the β-U structure [10], although the structures that they have compared differ only slightly. In view of this uncertainty an EAM structure determination was undertaken.
After relaxing the atomic positions, the cell volume, and the c/a ratio of the Moseley and Seabrook β-Ta cell, the EAM potential was found to stabilise a β-phase with a 1.018 nm × 1.018 nm × 0.533 nm unit cell, which is within 0.4 % of the experimental results [10], which include stresses. All atomic equilibrium positions ended up within 0.003 nm of their starting values. The average atomic cohesion energy of the β-phase increased by only 2 meV to 8.01 eV, 0.08 eV less than the value for the bcc phase (throughout this paper, the cohesion energy of an atom is defined as the negative of its potential energy after rapid cooling to 0 K). A difference of 0.08 eV may be somewhat large for a structure which is still stable in experiments at up to 1000 K. On the other hand, a phase transformation from β to bcc Ta was never observed in any simulation, as will be
discussed later. So even if the energy difference is somewhat large, it certainly does not result in any unduly strong tendency to form bcc Ta. This shows that the Johnson-Oh EAM potential can be used to describe the β-phase for our purposes. This is consistent with recent work of Mishin et al. on copper, who conclude that EAM potentials can be quite useful for non-equilibrium structures [14]. A relaxation of one of the alternative β-Ta structures proposed by Donohue and Einspahr showed that the atoms relaxed to the same positions as when starting from the Mosely and Seabrook structure. These results, although encouraging, should be treated with some caution. Creating a β-Mo structure also
resulted in a stable structure (average atomic cohesion energy 0.05 eV less than for bcc Mo), while to our knowledge β-Mo has never been observed experimentally. Therefore a small difference in cohesion energy predicted by EAM potentials is no guarantee of an actual metastable phase.
The various inequivalent atomic positions in β-Ta have substantially different individual cohesion energies. These are shown in fig. 2 along with the different relaxed vacancy formation energies at these positions. For comparison, in bcc Ta the cohesion energy is 8.09 eV and the relaxed vacancy formation energy is 2.86 eV. The different vacancy formation energies in β-Ta have important consequences for vacancy diffusion, as will be discussed in section 3.5.
Surface reconstruction was not observed on the β-Ta (001) surface, neither when the outer plane was a hexagonal plane, nor when it was an 'intermediate' plane. At 1000 K an intermediate
Fig. 1. Tetragonal β-Ta unit cell according to Moseley and Seabrook [10], seen from the <001> (left) and <010> directions (right). Cell dimensions are 1.019 nm × 1.019 nm × 0.531 nm. The atoms can be thought of as forming two (pseudo-) hexagonal planes with different orientations, while eight atoms lie halfway between the hexagonal planes. The grey atoms form the hexagonal plane at z = 0, the light atoms the hexagonal plane at z = c/2, and the black atoms lie in the 'intermediate' planes at z = c/4 and 3c/4 (the c-axis is taken in the z-direction).
type 1: Ec = 7.74 eV, Ev = 2.27 eV type 2: Ec = 7.79 eV, Ev = 2.47 eV type 4: Ec = 8.31 eV, Ev = 3.47 eV type 3: Ec = 8.09 eV, Ev = 2.99 eV
All atoms in intermediate planes have Ec = 8.06 eV, Ev = 2.81 eV 1 2 2 2 2 4 4 3 3 3 3
Fig. 2. View of a β-Ta crystal along the <001> direction with atoms shaded according to their cohesion energy Ec (negative of the potential energy). The top plane is a hexagonal plane. Also
indicated are the relaxed vacancy formation energies Ev at the different atomic sites. The atom
marked '1' is the atom in the centre of the unit cell in Fig. 1. Type-3 atoms have been marked with small white circles. A few hexagons have been drawn to indicate the orientation of the hexagonal planes.
outer plane does become disordered, because the low-coordinated atoms of the intermediate plane (with fewer than 8 neighbours) can easily move to other local potential energy minima in the same plane. Disordering is an appropriate term here, because the positions that these atoms assume are sometimes not an epitaxial extension of the substrate. This may have consequences for the adhesion of a copper overlayer.
The (001) surface energy is 1.92 J/m2 with a hexagonal plane as outer plane and 2.00 J/m2 with
an 'intermediate' plane as outer plane.
2.3.2 Thermal PVD film growth
PVD films nominally 2.5 nm thick have been grown on single-crystal bcc Ta (110) and β-Ta (001) substrates of 18 nm × 18 nm, held at 300 and 1000 K. As discussed above, the 1000 K simulations are intended to model a room-temperature experiment, while the 300 K simulations allow us to study deposition under conditions where the atomic jumps of activation energies higher than approximately 0.25 eV are frozen in (Eq. (1)). Fig. 3 shows the surface roughness during deposition. It is seen that increasing the deposition temperature from 300 K to 1000 K, i.e. activating the processes between 0.25 eV and 0.83 eV, only causes a modest roughness reduction. This is not surprising, since on the surface of a grown Ta film the majority of possible atomic jumps and rearrangements have activation energies above approximately 0.9 eV. This we will show later. The occasional low activation energy, like the 0.44 eV activation energy that will turn out to govern adatom diffusion on bcc (110) surfaces, has little influence on a grown surface, which consists of dozens of small surfaces with different orientations. The lower roughness at higher simulation temperature is caused mainly by a smaller number of low-coordinated surface atoms; no differences in the larger scale morphology of the films are observed. Fig. 3 shows that the surface roughness should be above ≈ 0.25 nm before the small temperature effect on the roughness becomes noticeable. When the roughness increases, the number of available low-energy rearrangement modes is seen to increase.
For comparison the roughness development during Mo growth on a bcc (110) Mo substrate is also shown in Fig. 3. The data were taken from earlier simulations, performed under the same
conditions [15]. The Mo film was found to develop a structure with much more open space than Ta (the density being 15-20 % below the bulk value), resulting in a pronounced columnar structure. It is not immediately obvious why a growing Ta film remains flatter than a Mo film. At room temperature, adatom mobilities on growing surfaces of both materials are negligible, and the cohesion energies are not much different. The main difference between Ta and Mo lies in the average shear modulus, which for Ta is only 54% of the value for Mo. Apparently the lower shear modulus allows easier atomic rearrangements during the condensation of new atoms.
The rising roughness curves in Fig. 3 indicate the onset of surface patterning. Fig. 4 shows the top and side views of a β-Ta and a bcc β-Ta film after 2.5 nm growth. It can be seen that there are clear
differences between the two films. The β-Ta film grows with an average lateral
Fig. 3. Surface roughness as a function of deposition thickness for aβ-Ta (001) film grown at 1000 K (bottom curve), aβ-Ta (001) film grown at 300 K (lower solid curve), a bcc Ta (110) film grown at 1000 K (dashed curve), a bcc Ta (110) film grown at 300 K (dotted curve), and a (110) Mo film grown at 350 K (upper solid curve).
feature size of approximately 3 nm, without a clear orientational preference. The bcc film is noticeably patterned, showing elongated features having a typical size of 4 nm in the <100> direction and 2 nm in the <110> direction. Apart from this, in the <110> direction the elevations are separated by much deeper depressions than in the <100> direction. This directional patterning was also found for Mo and is contributed to a preference of the material to form (100) surface nanofacets during growth. It should be further noted that the substrates used in this work are sufficiently large, since fig. 4 shows no signs of periodicities enforced by the boundary conditions.
The Ta films grew perfectly epitaxially. Only the β-Ta films grown at 300 K did contain some non-epitaxial clusters of atoms. However, given the perfect epitaxy found for the 1000 K
simulations (the disordering of an outer intermediate plane is not a significant effect in a growing film), the non-perfect epitaxy at 300 K is certainly just the result of the limited mobility of low-coordinated adatoms.
2.3.3 Ion bombardment
After deposition the β-Ta film grown at 1000 K was bombarded with a dose of up to 3 × 1018/m2
of 250 eV Ar ions, to investigate the stability of the β-Ta structure against ion bombardment. This had a noticeable influence on the film microstructure. Fig. 5 (a) shows the atoms marked bcc in the film during several stages of the bombardment. A number of surface atoms are spuriously
<100> <010> <010> <001> <100> <001> <110> <100> <100> <001> <001> <110>
Fig. 4. Views of 2.5 nm thick films grown at 1000 K on (left) β-Ta (001) and (right) bcc Ta (001) substrates. In the plan view, along the c-axis, atoms are shaded according to their z-coordinate. In the side views, atoms are shaded according to their cohesion energy. Note the orientational preference of the surface pattern of the bcc film, which is absent for the β-Ta film.
determined as bcc coordinated because they lack some neighbours. These atoms should be ignored; they are only shown to indicate the position of the surface. Down to approximately 4 nm below the surface, the Ar ions create 22-atom clusters of which 16 atoms are marked as bcc-like. Such a cluster of atoms is shown enlarged in fig. 5 (b). It can be seen that four heavily distorted bcc cells have formed. Such clusters of bcc-like cells always have one of two very particular orientations to the surrounding β-Ta lattice. A few smaller pockets were also marked as bcc, but
closer inspection showed that there were no visually discernable bcc clusters present at those sites (which did however, clearly deviate from β-Ta). The bcc marking in those cases is due to certain non-ideally tuned parameters within the rotational invariants algorithm. After the
bombardment, 2% of the atoms in the top 4 nm of the film became bcc-coordinated. The total energy input of the bombardment divided by the number of bcc-coordinated atoms it creates, is 170 eV per atom. It is tempting to conclude that on average 170 eV Ar impact energy is needed for the β→bcc transformation of a single Ta atom. However, although numerically correct, this does not describe the formation mechanism. What we observed is that an incoming argon ion either creates a typical 22-atom bcc cluster or does nothing at all in terms of phase transformation. With few exceptions, smaller or bigger bcc pockets do not appear. In combination with the fact that the 22-atom clusters have only a few well-defined orientations with respect to the β-Ta lattice, the conclusion seems justified that the 22-atom bcc cluster is a particularly stable cluster which, once nucleated, fits well within the β-structure. When it originates, the increased lattice stresses counteract the nucleation of new bcc clusters in its proximity. This is supported by the observations that neighbouring ion impacts never produced a larger cluster and that the cluster production rate decreased during the
bombardment.
Ar ions are incorporated in the film at substitutional sites, and after the 3 × 1018/m2 dose the
surface concentration had reached a value of 4 × 1017/m2, corresponding to an average impurity
concentration of 1.8 × 10−3 in the top 4 nm of the film. At the end of the bombardment the Ar
concentration had not yet reached its equilibrium steady state value. The concentration of empty vacancies at the end of the bombardment was roughly 10-4.
2.3.4 Annealing
To investigate the thermal stability of β-Ta films, the β-Ta film grown at 1000 K has been annealed for 23 ns at 800, 1000, and 1200 K. Fig. 6 (a) shows the average potential energy of the film atoms versus time; to convert to zero-temperature values, 3kT/2 was subtracted from each
(b)
Fig. 5. Local atomic crystal symmetry in aβ-Ta film during 250 eV Ar bombardment. In (a) the atoms shown are those which are determined as bcc-like coordinated atoms. Left: film before the bombardment. A number of surface atoms are spuriously determined as bcc coordinated because they lack some neighbours. These atoms should be ignored; they are only shown to indicate the position of the surface. Middle: film after a 6 × 1017
/m2
Ar dose. Right: film after a 3 × 1018/m2 Ar dose.
In (b) an Ar-created 22-atom cluster is shown, forming four heavily distorted bcc cells. Note the different distances between pairs of atoms. The orientation with respect to the β-Ta lattice can be best imagined by noting that the front and back planes of the cluster lie parallel to β-(001) planes. The two large atoms shown in black are atoms in the intermediate planes adjacent to the bcc cluster. By imagining such atoms at all four corners of the diamond formed by the four unit cells, two diamonds of atoms appear in the intermediate planes, which are the diamonds of intermediate atoms shown on the left in fig. 1.
4.0 nm (a)
curve prior to plotting. The results clearly indicate that the film thermally relaxes. According to the simple model of Eq. (1) atomic processes with activation energies up to 0.89, 1.11, and 1.33 eV are expected to take place three times at 800, 1000, and 1200 K, respectively. However, it is questionable whether these numbers are the best starting point for an analysis of long-term behaviour on a surface of complex geometry. What one observes in the simulations is a multitude of atomic processes that slowly die out over time; examples are diffusing adatoms being captured by more stable clusters, and clusters that rearrange to more stable configurations. A better
approach is to assume that the observable in question (the potential energy) relaxes via a large number of relaxation modes, each characterized by its own relaxation time τ. Because the thermal activation of atomic vibrations is the driving force of the relaxation, the temperature dependence of each τ is given by an expression similar to Eq. (1),
τ = 1 ν0
eQ / kT
. (2)
When the collection of available relaxation modes is described as an activation energy spectrum
P(Q, t), and the nonequilibrium part of the observable, EP(t)−EP(∞), is taken to be the integral over
this spectrum, the observable will relax according to
EP(t) − EP(∞) = P(Q, 0)e−ν0te− Q/ kTdQ
∫
, (3)provided that each relaxation mode decays exponentially, P(Q, t) = P(Q, 0)exp(−t/τ). The second factor in the integrand of Eq. (3) is essentially a step function of Q (in fact it is slightly
broadened), where the step occurs at the value Qs = kT ln ν0t. In other words, as time proceeds the
initial spectrum P(Q, 0) becomes increasingly exhausted, because Qs increases. As a result, the
potential energy steadily decreases,
EP(t) − EP(∞) = P(Q, 0)
kT lnν0t
∞
∫
dQ. (4)To verify this approach, we must be able to find a value of ν0 that causes all isotherms to fall on a
single curve when they are plotted as a function of kT ln ν0t. Apart from the slight broadening of
the step function, this curve would then be equal to the negative of the initial spectrum P(Q, 0).
Fig. 6. (a) Average atomic potential energy as a function of time for three β-Ta films, annealed at 800, 1000 and 1200 K. To convert to zero-temperature values, 3kT/2 was subtracted from each isotherm prior to plotting. (b) Same isotherms plotted as a function of Qs = kT ln ν0t, with ν0 = 1 × 10
Fig. 6 (b) shows that, using the value ν0 = 1 × 10
14 s−1, this analysis leads to a quite consistent
result, since a single curve is clearly obtained. The value of ν0 lies well within the range of other
frequency values for this system, as we will see. The shape of the energy spectrum shows that, as was mentioned earlier, activation energies below ≈ 0.9 eV are almost absent on an as-grown Ta surface, and that higher activation energies become increasingly important for thermal stability issues. There are no signs that the spectrum falls off to zero within the energy range probed. However, from the viewpoint of room temperature applications the curve in Fig. 6 (b) extends far enough. The analysis predicts that at 300 K the spectrum will be activated up to Qs = 0.89, 1.19,
and 1.34 eV, respectively, for film lifetimes of 10 s, 10 d, and 10 y.
A β → bcc transformation was never observed during thermal annealing. This is not surprising considering the impact energies needed in the ion bombardment to induce this transformation. Overall surface roughness was reduced by only a few hundredth of nm, showing that one needs to overcome barriers beyond 1.5 eV to decrease roughness significantly.
The ion-bombarded β-Ta film has also been annealed, at 800 K. Annealing reduced the number of atoms marked as bcc-coordinated, see fig. 7. The initial strong decrease is due to the disappearance of the small clusters mentioned earlier. Because of the short time in which these small clusters disappear, their accumulation during the bombardment was probably not realistic. With a realistic ion impact rate, there would have been time for the small clusters to disappear before the next impact. The slow decrease at later annealing time is predominantly due to surface atoms marked as bcc which move to “non-bcc” positions. Importantly, none of the typical 22-atom bcc clusters disappeared or formed aggregates. This again underlines their special embedding in the β-Ta lattice. We conclude that their presence is a realistic effect of the ion bombardment.
2.3.5 Vacancy and adatom diffusion
Apart from the 'macroscopic' relaxation behaviour addressed in the previous section, the diffusion behaviour in a number of well-defined systems containing a vacancy or an adatom has also been studied. It should be noted that for the β-Ta systems, simulation temperatures were often several hundred degrees higher than the reported β → bcc transformation temperature (800-1000 K, depending on deposition conditions [16, 17]). Such high temperatures were needed in order to observe an adequate number of atomic jumps during a feasible simulation time. However, the phase transformation was not observed in the simulations, although the simulation times were long enough. Very probably these test systems were too small to undergo the phase transformation in a simulation box of fixed size (and in some cases without a free surface).
Only a few systems obeyed simple, single-activation-energy behaviour. For these systems migration activation energies and prefactors are listed in Table I, along with the scarce
experimental data. The table shows that the frequency value 5 × 1013 s−1 used in Eq. (1) and the
value 1 × 1014 s−1 resulting from Eq. (4) lie well within the range of values obtained for these
single processes. The values for vacancy diffusion in bcc Ta agree well with the experimental data. As expected, the activation energy is somewhat higher than the value 1.2 eV found from simulations on Mo [15]; the ratio of the two values is almost the same as the ratio of the cohesion energies of Ta and Mo (1.19). However, contrary to Mo adatom diffusion on Mo(100), which was
Fig. 7. Number of bcc atoms in the ion-bombarded film during annealing at 800 K.
found to take place via a simple hopping mechanism with an activation energy of more than 1.4 eV, Ta adatom diffusion on bcc Ta(100) is observed to occur predominantly through an exchange mechanism with an activation energy of 1.4 eV. At this moment, we do not fully understand why this difference occurs. Adatoms on bcc Ta (110) are very mobile (Q = 0.44 eV) and migrate through hopping. For comparison, from the Mo simulations a value of 0.55 eV was found for Mo adatoms on Mo (110). Somewhat surprisingly this value is a little greater than that for Ta.
Diffusion on β-Ta (001) surfaces is not governed by a single adatom migration energy. As the temperature increases above 1400 K, surface atoms are observed to break free from their positions in a hexagonal outer plane, often in groups, within 71 ns simulation time. The atoms then move to several different local potential energy minima, which are not always epitaxial extensions to the substrate. Adatoms appear to have hardly any more mobility on such surfaces than any other atom in the top plane. The same is true for adatoms on intermediate outer planes.
The most interesting diffusion phenomena were observed during monovacancy diffusion in bulk β-Ta. Because of the different vacancy formation energies, vacancies prefer to remain at those sites which have the lowest vacancy formation energy. Fig. 8 again shows the β-Ta structure coloured according to the atomic cohesion energies, but this time (compared to Fig. 2) a number of atoms has been removed. It can be seen that a vacancy in a position having a low vacancy formation energy (type-1 or type-2 in Fig. 2) is surrounded by atoms with high cohesion energies (the diamond-shaped rings of type-3 and type-4 atoms) on all sides in the (001) plane. The positions in the adjacent hexagonal planes opposite the vacancy also have high formation energies, as fig. 8 shows. Therefore, vacancies in type-1 and type-2 positions are essentially trapped by a surrounding potential energy 'cage'. A vacancy can move relatively easy among type-1 and type-2 positions inside the cage, which have either exactly or almost exactly the same vacancy formation energies. However, a vacancy movement over a longer distance always implies the motion of atoms from high to low cohesion energies. This creates an extra energy barrier for long distance vacancy migration, see Fig. 9. From this one could expect that vacancy diffusion would be almost completely frozen out
Table I. Prefactors ν0 and migration activation energies Q from simulations
and experiments. ν0 (Hz) simulations Q (eV) simulations ν0 (Hz) experimental Q (eV) experimen tal bcc Ta vacancy diffusion 2 × 1014 1.4 1.5 × 1015 (a) 1.45 (b)
bcc Ta (100) adatom diffusion 6 × 10 13 1.4 - -bcc Ta (110) adatom diffusion 5 × 10 12 0.44 -
-β-Ta vacancy jump between
intermediate planes 2 × 10
13 1.0 -
-(a) Ref [18].
(b) The number is obtained by subtracting the literature average of the values for the relaxed vacancy formation energy (2.95 eV) [19,20] from the value for the self-diffusion energy (4.4 eV) [18]. For comparison: the current EAM value for the relaxed vacancy formation energy is 2.86 eV.
Fig. 8. Two β-Ta (001) hexagonal planes with atoms shaded according to their cohesion energies. The intermediate plane between the two hexagonal planes is not shown. Dark atoms have high cohesion energies. Compared to Fig. 2, a cluster of atoms with low cohesion energies on the left top side has been removed in order to show the atoms with high cohesion energies in the next lower plane.
in β-Ta, limiting diffusion to those rare cases when a vacancy jumps from a position with a low formation energy into a position with a high formation energy, after which the vacancy would either quickly jump back – thereby annihilating the diffusion attempt – , or move another step in the same direction, to reclaim the invested energy difference. These rare cases have indeed been
seen. However, an interesting alternative mechanism was observed, by which the system is able to circumvent the extra energy barrier. The principle is straightforward. When two neighbouring atoms jump simultaneously, the first of the pair moving into an existing vacancy and the second one following, then the vacancy moves two places. If the vacancy thereby moves from a position with low vacancy formation energy to another position of low formation energy, the energy increase of the first atom is concurrently compensated by the movement of the second atom. This
principle is illustrated in Fig.10 (a). Such a mechanism was observed to take place in β-Ta in the <001> direction. Neighbouring atoms at type-2 and type-3 sites in adjacent hexagonal planes, the type-3 atom situated opposite a vacancy at a type-2 site in the next hexagonal plane, sometimes moved simultaneously as a dimer, allowing the vacancy to jump into a new cage two planes further. A snapshot from a dimer jump in progress is displayed in Fig. 10 (b). For reasons that are as yet undiscovered, type-2 and type-3 atoms in the same hexagonal planes never moved in pairs. Because it proved very difficult to disentangle the different processes from the data (dimer jumps, various in-cage jumps and the sporadic single vacancy jumps to adjacent hexagonal planes), the migration energies and frequency factors for diffusion involving the hexagonal planes could not be determined. Vacancy jumps from one intermediate plane to the next intermediate plane were also found. Since in this process only one type of jump is involved, it was possible to determine the activation energy and prefactor, see table I. Vacancy jumps from intermediate to hexagonal planes or vice versa were never observed.
The overall result is highly anisotropic diffusion behaviour. The dimer mechanism only allows one-dimensional movement, into the <001> direction, and the same is true for vacancy diffusion via intermediate planes. Vacancies can jump to other positions in the same cage, or in the same cage two planes further, but they cannot jump into a different cage, either in the same (001) plane
type 1 type 2 type 3 type 4 inter- medi-ate 1.4 eV 0.2 eV 0.3 eV 0.7 eV 1.2 eV 1.0 eV
Fig. 9. Potential energy barriers for monovacancy diffusion in bcc (left) and β-Ta (right, schematically). In bcc Ta there is just the barrier due to the periodic nature of the lattice. In β-Ta there are additional contributions due to the energy difference between inequivalent vacancy sites.
Fig. 10. Vacancy diffusion through a dimer mechanism. The <001> axis is displayed horizontally. Left: schematic representation of dimer diffusion, through which the high energy barriers shown in fig.9 are lowered (black = type-2 atom, white = type-3 atom). Right: snapshot from a simulation during a dimer jump. The two dimer atoms are the large atoms indicated by arrows. The dimer atoms move from one hexagonal plane to the adjacent hexagonal plane. As can be seen, the dimer is somewhat stretched during the jump.
