Thin Solid Films 332 (1998) 109±112
Thermal diffusion in molecular dynamics simulations of thin ®lm diamond deposition B.A. Pailthorpe a, b, c,*, D. Mitchell b, N.S. Bordes a, b a
Department of Applied Physics, University of Sydney, Sydney, NSW 2006, Australia b Sydney VisLab, University of Sydney, Sydney, NSW 2006, Australia c San Diego Supercomputer Center, P.O. Box 85608, San Diego, CA 92138-5608, USA
Abstract Molecular dynamics simulations of carbon atom depositions are used to investigate energy diffusion from the impact zone. A modi®ed Stillinger±Weber potential models the carbon interactions for both sp2 and sp3 bonding. Simulations were performed on 50 eV carbon atom depositions onto the (111) surface of a 3:8 £ 3:4 £ 1:0 nm diamond slab containing 2816 atoms in 11 layers of 256 atoms each. The bottom layer was thermostated to 300 K. At every 100th simulation time step (27 fs), the average local kinetic energy, and hence local temperature, is calculated. To do this the substrate is divided into a set of 15 concentric hemispherical zones, each of thickness one atomic diameter (0.14 nm) and centered on the impact point. A 50-eV incident atom heats the local impact zone above 10 000 K. After the initial large transient (200 fs) the impact zone has cooled below 3000 K, then near 1000 K by 1 ps. Thereafter the temperature pro®le decays approximately as described by diffusion theory, perturbed by atomic scale ¯uctuations. A continuum model of classical energy transfer is provided by the traditional thermal diffusion equation. The results show that continuum diffusion theory describes well energy diffusion in low energy atomic deposition processes, at distance and time scales larger than 1.5 nm and 1-2 ps, beyond which the energy decays essentially exponentially. q 1998 Published by Elsevier Science S.A. All rights reserved. Keywords: Simulations; Diamond; Stillinger±Weber; Thermal diffusion
1. Introduction There has been signi®cant interest [1±3] in the growth of synthetic diamond ®lms by low energy carbon plasmas and ion beams directed onto low temperature substrates. A wide variety of experimental techniques [4] are providing valuable data on the bonding states, structure and properties of diamond ®lms deposited by many techniques. One route to theoretical understanding is provided by molecular dynamics (MD) computer simulations which provide insight into atomic scale processes at growing diamond surfaces and detailed structural descriptions of the material. The work described here compares computational simulation studies to traditional continuum diffusion theories applied to energy dissipation in growing diamond thin ®lms, and so builds a bridge to a wider class of theoretical methods that could be applied to furthering our mechanistic and structural understanding of diamond thin ®lms. In order to progress MD studies a description of interatomic interactions is required. We have developed a Stillin* Corresponding author Tel.:161 2 93513863; fax:161 2 9351 7726; email:
[email protected].
ger±Weber (SW) class of interatomic potential [5] suitable for modeling both sp 3 and sp 2 interactions in diamond-like carbon [6]. Our potential function was derived originally from Hartree±Fock calculations for small carbon clusters in tetrahedral [6,7] and trigonal con®gurations [8,9]. It has been applied to modeling bulk diamond, carbon clusters and thin ®lm forms of high density amorphous carbon [7,9±12]. The SW potential is an empirical inter-atomic potential which was developed originally for silicon [5] and is a sum over physically intuitive two-body bond stretching and three-body bond bending interactions. This potential has been applied successfully to the study of bulk [13±15] and surface [16±18] properties of silicon, to sputtering of silicon surfaces [19], to small silicon clusters [20,21] and to bulk germanium [22]. Another [23,24] interatomic potential, developed by Tersoff, which includes the effects of coordination number has improved the description of silicon clusters and surfaces, has been used with some success to simulate growth of amorphous carbon ®lms [25], and has been extended by Brenner [26] to study carbon ®lms grown from hydrocarbon species. A tight binding MD scheme, using an approximate localized wave function, has been applied to amorphous carbons [27,28], but fails to predict
0040-6090/98/$ - see front matter q 1998 Published by Elsevier Science S.A. All rights reserved. PII S0 040-6090(98)010 13-X
110
B.A. Pailthorpe et al. / Thin Solid Films 332 (1998) 109±112
Fig. 1. Schematic illustration of the simulated system, comprising a diamond (111) substrate at 300 K, an incident 50 eV neutral carbon atom and hemispherical zones for temperature calculations.
the observed sp 3 fractions and structures [29,30]. More fundamentally, the Car±Parrinello ab initio molecular dynamics (CPMD) technique has been used to investigate bulk molten and quenched carbon [31±33] and high density tetrahedrally coordinated amorphous carbon (taC) [29,30], and has provided unique and valuable insights into the bonding states and structures present. However, with presently available computer capacity it does not appear feasible yet to apply the Car±Parrinello technique to surface and thin ®lm studies because of the reduced symmetry. The present study continues with the semi-classical approach and is investigating carbon ®lms using a combined sp 3-sp 2 SW potential. We note parenthetically that other [34], essentially thermodynamic, theories explore the dependence of sp 3 fraction on deposition rates. 2. Theory and procedure The functional form of the SW potential has been reported widely [5±9,11±22] so it will not be repeated here. Our previously presented results show that, for the diamond con®guration, the parameters of the SW potential most appropriate to carbon are: A 5:298336, B 0:503230, a 1:9, l 1:10, gd 24:0, sd 0:1378 nm and ed 3:214 eV. Here sd corresponds to the atomic diameter and e d to the pair sp 3 bonding energy; these values
of s d and e d form the basis for reduced units herein. For sp 2 bonding in graphitic structures our group has calculated and tested SW potential parameters, ®nding A, B and a values as above, and gg 15:0, s g 0.1237 nm and e g 4.299 eV [8]. Being ¯exible between sp 3 and sp 2 bonding, the mixed potential can describe better the under-coordinated states observed in small carbon clusters, molten carbon and surface deposition processes [12]. Ultimately full CPMD simulations [29,30] are required but are computationally prohibitive for the moment. The molecular dynamics technique is widely known [35] and has been previously described for the present application [7,9,11]. Minimum-image periodic boundary conditions, neighbor listing and the Verlet algorithm are used with integration time steps of 0.27 fs. After the deposition of each atom the system dynamics is followed for times of 1±5 ps. The starting con®guration of atoms derives from a diamond lattice at 0 K, with a density appropriate to that of natural diamond, 3.52 g/cm 3. The temperature of the system was increased up to 300 K by randomly assigning and gradually rescaling all atomic velocities. The diamond (111) surface was modeled by 11 layers of 256 atoms each, periodically repeated in the x±y plane. The fundamental cell of 2816 atoms was contained in a box of dimensions 3:8 £ 3:4 £ 1:0 nm. The bottom layer of atoms was rigidly ®xed (to an in®nitely massive substrate) to preserve the integrity of the model ®lm. The next layer was `thermostated' to the substrate temperature ^25% by periodic rescaling of its atomic velocities which were modulated additionally by a random phase to eliminate spurious phonon re¯ections. The procedure corresponds to rates of heat removal to the substrate and to deposition rates which are orders of magnitude greater than is experimentally realized. For instance a typical cooling rate of 10 15 K/s was imposed on the model thin ®lm, compared to laser glazing rates [36] of 10 12 K/s. That said, comparisons to CPMD simulations of bulk carbon quenching processes [37] indicate that the taC metastable con®gurations are actually frozen in by such rapid quenching, and are annealed out by longer times at elevated temperatures. The outer layer of atoms presented one dangling bond per atom to the vacuum. No hydrogen is involved in the experimental fabrication of taC ®lms via vacuum arc techniques [38±40]. Surface reconstruction was not expected nor observed. Energy transfer during diamond thin ®lm growth was investigated by introducing a single carbon atom at a height of two atomic diameters above the surface, traveling at normal incidence towards the free surface with the prescribed incident kinetic energy of 50 eV. The dynamics of the entire system was then monitored. The simulation was run until all the energy imparted to the ®lm by the incoming atom was dissipated by the thermostating process described above. The incident atom was aimed at the centre of the cell, approximately over a sub-surface atom which presented no dangling bonds
B.A. Pailthorpe et al. / Thin Solid Films 332 (1998) 109±112
Fig. 2. The pro®le of local temperatures in a diamond substrate maintained at 300 K, as a function of time and radius from the impact zone of a 50 eV neutral carbon atom. Predictions of continuum diffusion theory are shown on the left and results of MD simulations on the right.
3. Results The average atomic kinetic energy, at every 100th MD time step (i.e. 27 fs), is calculated in 15 hemispherical
111
shells, akin to onion skin layers, of approximately one atomic diameter, s d spacing. The geometry of the system is illustrated in Fig. 1. Each shell typically contains at least 5±10 atoms, aside from the central impact zone, with shells at larger radii containing many more atoms. The average kinetic energy in each shell is then used to calculate P a local temperature via the kinetic theory relationship: mv2 3=2 NkT. The continuum theory predictions are calculated using the standard 3D diffusion equation [41] solved in a semi-in®nite half space representing the substrate. Using an average of known values of thermal conductivity for diamond-like carbons [42,43], here taken to be 1000 (w/ mK), theoretical temperature pro®les, T(r,t) are calculated and plotted for comparison with the MD simulations. The energy delivered by the incoming atom to the substrate propagates, and ultimately diffuses, away from the initial impact site in two distinctive regimes: (1) an initial thermal spike of 40 000 K is `instantly' (within a few tens of atomic collisions) transformed into a molten zone of some ®ve atoms at approximately 10 000 K, which rapidly cools to 3000 K within about 0.2 ps; and (2) a longer time decay, over some 2 ps, from 3000 K to the ambient temperature of 300 K. These results are presented as 3D surface plots of the temperature pro®le, T(r,t). These color plots are generated using the scienti®c visualization software package, AVS [44] and are presented here in grayscale. The original 3D color plots are available at our web site (www.vislab.usyd.edu.au/research/). Note that the binned (histogram) MD data is interpolated and smoothed by the software (AVS `®eld to mesh' module). Fig. 2 shows the full decay pro®le, including the transient phase from a time of 0.08 ps after the atomic impact, but excludes the initial impact event which is dominated by collisional processes. Within 5±10 atomic diameters of the impact zone large scale ¯uctuations dominate T(r,t), even for long times, approaching 3 ps. Fig. 3 shows the long time decay from 1.13 ps onwards. Again there are large scale temperature ¯uctuations within the ®rst 5±10 atomic diameters. Beyond that radius (r . 1:5 nm) such ¯uctuations die off rapidly. Then T(r,t) follows the trends described by the continuum diffusion theory perturbed by decaying ¯uctuations. 4. Discussion and conclusions
Fig. 3. The pro®le of local temperatures, calculated as in Fig. 1, plotted only for longer times (@t . 1.13 ps) to show decaying ¯uctuations.
The present calculations provide a bridge between numerical simulations of atomic processes in thin ®lm carbon structures, and traditional continuum theories. The approach is semi-classical, rigorous and fully three-dimensional and can treat interfaces and free surfaces approximately. However, it is limited by the classical assumptions. Full quantum (CPMD) simulations are not yet practical for these systems. The MD simulations reported herein describe energy diffusion after single neutral carbon atom impacts on a model diamond (111) surface at
112
B.A. Pailthorpe et al. / Thin Solid Films 332 (1998) 109±112
300 K using a modi®ed Stillinger±Weber potential, parameterized for both sp 3 and sp 2 bonding in carbon. The results show that continuum diffusion theory describes well energy diffusion in low energy atomic deposition processes, at distance and time scales larger than 1.5 nm and 1±2 ps, beyond which the energy decays essentially exponentially Color versions of Figs. 1±3, along with QuickTime animations (or a videotape) of our previous simulations of thin ®lm diamond growth are available also at www.vislab.usyd.edu.au/research/. Acknowledgements The valuable discussions and collaborations of Dr. D.R. McKenzie and P. Guan are gratefully acknowledged. Financial support was provided by the Australian Research Council and by His Royal Highness Prince Nawaf Bin Abdul Aziz of the Kingdom of Saudi Arabia through the Science Foundation for Physics. The San Diego Supercomputer Center hosts a Senior Fellowship for BAP and provides computational resources which are gratefully acknowledged. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10]
J.C. Angus, C.C. Hayman, Science, 241 (1988) 913. W. A Yarbrough, R. Messier, Science 247 (1990) 688. D.E. Koshl Jr., Science 250 (1990) 1637. J. Robertson, Adv. Phys. 35 (1986) 318. F.H. Stillinger, T.A. Weber, Phys. Rev. B 31 (1985) 5262. P. Mahon, B.A. Pailthorpe, G. Bacskay, Phil. Mag. B 63 (1991) 1419. B.A. Pailthorpe, P. Mahon, Thin Solid Films 193/194 (1990) 34. A Bensan, Thesis, School of Physics, University of Sydney, 1990. E.G. Gerstner, B.A. Pailthorpe, J. Non-Cryst. Solids 189 (1995) 258. D.R. Mckenzie, D. Muller, B.A. Pailthorpe, Phys. Rev. Lett. 67 (1991) 773. [11] B.A. Pailthorpe, J. Appl. Phys. 70 (1991) 543. [12] P. Guan, Ph.D. Thesis, University of Sydney, 1997.
[13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43] [44]
M.D.Kluge, J.R Ray and A. Rahman, J. Chem. Phys., 85 34 (1986). M.D. Kluge, J.R. Ray, A. Rahman, Phys. Rev. B 36 (1987) 4232. W.D. Luedtke, U. Landman, Phys. Rev. B 40 (1989) 1164. F.F. Abraham, J.Q. Broughton, Phys. Rev. Lett. 56 (1986) 734. M. Schneider, I.K. Schuller, A. Rahman, Phys. Rev. B 36 (1987) 1340. E.T. Gawlinski, J.D. Gunton, Phys. Rev. B 36 (1987) 4774. R. Smith, D.E. Harrison Jr., B. Garrison, Phys. Rev. B 40 (1989) 93. E. Blaisten-Borojas, D. Levesque, Phys. Rev. B 34 (1986) 3910. B.P. Feuston, R.K. Kalia, P Vashista, Phys. Rev. B 36 (1987) 6222. K. Ding, H.C. Anderson, Phys. Rev. B 34 (1986) 6987. J. Tersoff, Phys. Rev. Lett. 56 (1986) 632. J. Tersoff, Phys. Rev. Lett. 61 (1988) 2879. H.-P. Kaukonen, R.M. Nieminen, Phys. Rev. Lett. 68 (1992) 620. D.W. Brenner, Phys. Rev. B 42 (1990) 9458. C.Z. Wang, K.M. Ho, Phys. Rev.Lett. 71 (1993) 1184. C.Z. Wang, K.M. Ho, Phys. Rev.Lett. 72 (1994) 2667. N.A. Marks, D.R. McKenzie, B.A. Pailthorpe, M. Bernasconi, M. Parrinello, Phys. Rev. Lett. 76 (1996) 768. N.A. Marks, D.R. McKenzie, B.A. Pailthorpe, M. Bernasconi, M. Parrinello, Phys. Rev. B 54 (1996) 9703. J. Galli, R.M. Martin, R. Car, M. Parrinello, Phys. Rev. Lett. 60 (1988) 204. J. Galli, R.M. Martin, R. Car, M. Parrinello, Phys. Rev. Lett. 63 (1989) 988. J. Galli, R.M. Martin, R. Car, M. Parrinello, Science 250 (1990) 1547. Y. Yin, D.R. McKenzie, Thin Solid Films 280 (1996) 95. B.J. Berne, Statistical Mechanics, Plenum, New York, 1977. R.Zallen, The physics of Amorphous Solids, Wiley, New York, 1983. N.A. Marks, Phys. Rev. B 65 (1997) 2441. S.D. Berger, D.R. McKenzie, P.J. Martin, Phil. Mag. Lett. 57 (1988) 285. D.R. Mckenzie, D. Muller, B.A. Pailthorpe, et al., Diamond Relat. Mater. 1 (1991) 51. D.R. McKenzie, Y. Yin, N.A. Marks, C.A. Davis, B.A. Pailthorpe, G.A.J. Amaratunga, V.S. Veerasamy, Diamond Relat. Mater. 3 (1994) 353. J. Crank, The Mathematics of Diffusion, Oxford University Press, New York, 1975. R.C. Weast (ed.), Handbook of Chemistry and Physics, CRC Press, Bocs Raton, FL, 1975. C.J. Morath, et. al, J. Appl. Phys. 76 (1994) 2636. Advanced Visual Systems, Inc. Waltham, MA.