Molecular-dynamics simulation of the structure and diffusion ...

PHYSICAL REVIEW B

VOLUME 61, NUMBER 14

1 APRIL 2000-II

Molecular-dynamics simulation of the structure and diffusion properties of liquid silicon Z. G. Zhu and C. S. Liu Laboratory of Internal Friction and Defects in Solids, Institute of Solid State Physics, Chinese Academy of Sciences, P.O. Box 1129, Hefei 230031, People’s Republic of China 共Received 26 October 1999兲 Using molecular-dynamics methods and the activation-relaxation technique, we have investigated the inherent structure and diffusion properties of liquid silicon. With increasing density, the 52° and 60° peaks 共attributed to long bonds兲 in the bond-angle distribution functions decrease in height, while the spread main peak 共mainly related to the bonds containing covalent character兲 increases and moves towards the tetrahedral angle. The change in density does not give rise to a clear change in the diffusion constants. With the change of temperature, the diffusion coefficients obtained from the average mean-square displacement can be fitted by Arrhenius equation. The fit yields an activation energy of 0.92 eV and a pre-exponential factor of 30.8 ⫻10⫺3 cm2 s⫺1. However, the activation energy, which is determined from the activation-relaxation technique using a Metropolis accept-reject criterion with a fictitious temperature of 0.5 eV, is in the range of 0.22 to 1.0 eV and shows a steep increase at low temperature. The very large pre-exponential factor suggests that the interatomic forces obtained from the Tersoff potential are very strong. The information obtained in this paper is consistent to some extent with the recent experimental results of some physical properties of liquid silicon.

I. INTRODUCTION

Silicon is one of the most complex elemental liquids. Crystalline and amorphous silicon form covalent bond while liquid silicon has metallic character. Upon melting, the density of silicon increases by ⬃10%,1 and the average coordination of liquid silicon increases from 4 in the crystal to over 6.2,3 The experimentally determined atomic structure factor of liquid silicon exhibits peculiar features, and strongly differs from that characteristic of simple liquid metals. Although the conductivity of liquid silicon exceeds that of the solid by a factor of 20,4 the low coordination of liquid silicon indicates that silicon appears not to be a ‘‘free-electron’’ metal, but contains ‘‘covalent character.’’ Recently, Kimura and co-workers have observed the anomalous temperature dependence of physical properties of liquid silicon in the temperature range just above the melting point, such as the density,5 the viscosity,6 the surface tension,7 and the electrical resistivity.8 For instance, when temperatures are lower than about 1430 °C, the density and the viscosity show a steep increase. The high-temperature x-ray-diffraction study of melt silicon by Waseda et al.9 indicate that, the structure factors have a characteristic small hump on the higher wave side of the first peak and such specific feature becomes slightly obscure as the temperature increases. A thorough understanding of these behaviors is technologically important because the anomalous temperature range is involved in the crystal growth from molten silicon. Molecular-dynamics 共MD兲 simulations provide a useful technique to analyze atomic structures of the liquid phase. Although ab initio simulations for atomic-scale studies have made great progress, they are still difficult to apply for largescale systems. Therefore many workers have developed several empirical interatomic potentials for silicon.10 Recently, the potential proposed by Tersoff11 has been applied to study the validity of this potential for liquid silicon:12,13 The pair0163-1829/2000/61共14兲/9322共5兲/$15.00

PRB 61

correlation function produced by the Tersoff potential is consistent with that obtained by the Stillinger-Weber potential and ab initio calculations except for a slight quantitative difference; the bond-angle distribution function of the Tersoff liquid is entirely different from that of the Stillinger-Weber liquid, but it is much similar to that of ab initio simulations; The Tersoff potential yields a similar tendency to velocity autocorrelation function obtained by ab initio method, whereas the oscillation period is shorter than that of ab initio calculations. Using this potential, we have recently carried out MD simulations of the local inherent structure of liquid silicon at different temperatures.14 Our result shows that the 50°–60° peak in bond angle distributions decomposes into two peaks located at 52° and 60°, respectively, and a new peak at 75° appears. Similar to some physical properties of liquid silicon, it has been observed that the anomalous temperature dependence of the height of 52° peak and the probability of the covalent bonds whose bond angle is greater than 67°. This anomalous behavior may play an important role on the anomalous feature of some physical properties in liquid silicon such as electrical resistivity. However, how does the anomalous change in density with temperature give rise to the change in the structure? The anomalous increase in viscosity with decreasing temperature implies that the property of diffusion will be unusual. The activationrelaxation technique 共ART兲 has been successfully applied to amorphous semiconductors and metallic glasses.15,16 Can we study the diffusion property of liquid silicon by using this technique or can we find some information to help us understand the diffusion properties of liquid silicon? In the present work, we performed MD simulations using Tersoff potential and applied the ART in order to investigate the inherent structures with varied density and diffusion properties of liquid silicon. This paper is organized as follows. The computational methods are presented in Sec. II. In Sec. III we provide the calculated results and discussions. Conclusive remarks are summarized in Sec. IV. 9322

©2000 The American Physical Society

PRB 61

MOLECULAR-DYNAMICS SIMULATION OF THE . . .

9323

II. COMPUTATIONAL METHODS

The same computational procedure as described in our previous work14 was employed. The Tersoff potential11 was employed to model the silicon atomic interactions. MD simulations under constant volume and constant temperature conditions were carried out for N⫽512 atoms in a cubic box with periodic boundary conditions. The Newtonian equations of motion were integrated using a velocity-Verlet algorithm with a time step ⌬t⫽2⫻10⫺3 ps. After the equilibrium liquid state was obtained, another run of 12 000 time steps at the given temperature or density was performed to obtain 60 configurations and to determine the diffusion coefficient from the atomic average mean-square displacement. These configurations were relaxed to their closest local minimum by a conjugate gradient energy minimization and then we employed the obtained configurations to analyze the inherent structure following Stillinger idea.17 In calculation, we fix the density at 2.57 g/cm3, and change the temperatures to 2800, 3000, 3200, 3400, 3600, and 3800 K, or fix the temperature at 3000 K, and adjust densities to 2.61, 2.59, 2.57, 2.55, 2.53 g/cm3. The activation-relaxation technique 共ART兲 presented by Barkema and Mousseau allows the system to evolve following well-defined paths between local energy minima 共for details of the ART see Ref. 18兲. This technique proceeds in two steps, activation and relaxation. That is, the local energy minimum is first moved to a nearby saddle point, then the system is pushed over the saddle point and relaxed to a new minimum. We start from the obtained 60 inherent structural configurations. And then we apply ART iteratively. ART moves a local minimum to a saddle point using a Metropolis accept-reject criterion with a fictitious temperature of 0.5 eV. An event in the ART is determined as a move from a local energy minimum to another nearby minimum. If an ART move is accepted, it is a successful event. At a certain temperature, 1200 successful events are performed. So the average local saddle energy and the average local minimum energy can be determined from these successful ART events, and then we can determine the average energy barrier or the activation energy. III. RESULTS AND DISCUSSIONS

The triplet correlation, which is particularly important because liquid silicon retains some covalent bonding effect and has directional forces, is conveniently measured by the bondangle distribution function g 3 ( ␪ ). Here ␪ is the angle between two vectors that join a central atom with two neighbors. The previous molecular-dynamics simulations indicate that14 the height of the 52° peak exhibits an anomalous temperature dependence, at first increases and then decreases with increasing temperature, and the probability of those covalent bonds with bond angle greater than 67° also shows an anomalous decrease at a certain temperature. They qualitatively give a better explanation of why the electrical resistivity of liquid silicon shows a local minimum in the temperature range 1450–1500 °C. However, we cannot well interpret the anomalous behavior of other physical properties such as the viscosity. Here we should point out that the anomalous temperature range of the density, the viscosity, and the surface tension is almost identical, i.e., about 1430 °C, while the

FIG. 1. Bond-angle distribution function g 3 ( ␪ ) at five different densities. The cutoff distance is 3.1 Å. Inset: the relative probability P(52°) of bond angle 52° 共⫾3°兲 found in the system at five different densities.

anomalous temperature range of the electrical resistivity is close to 1500 共see Fig. 3 in Ref. 8兲. According to Kimura and co-workers’ experiments, this difference in the anomalous temperature range does not originate from the experimental uncertainty. We have calculated the bond-angle distribution functions from the inherent structures at five different densities, as shown in Figs. 1 and 2 with the cutoff distance of 3.1 and 2.9 Å, respectively. As can be seen in Figs. 1 and 2 and two insets, the characteristic features are the same as previously reported. That is, the 50°–60° peak decomposes into two peaks that are located at 52° and 58° and a small new peak appears to be located at 75°. With the choice of smaller cutoff distance, the 52° peak disappears in the bond-angle distribution functions and the 58° peak is now located at 60°. It can also be seen that, with the increase in density, the height of both 52° and 60° peaks decreases and so does the height of 75° peak; while for the spread main peak, its height is increased and its location moves towards the tetrahedral angle. A bond-angle distribution function, calculated with choice of a smaller cutoff distance close to the covalent bond length, nearly shows a single peak close to the tetrahedral angle 共this is in agreement with that reported in Ref. 19兲, and the peak height increases with increasing density. These results indicate that the bonds with relatively long bond length decrease and the bonds containing covalent character increase with increasing density. The diffusion coefficient has been calculated from the average mean-square displacement. Table I provides the diffusion coefficients at different temperatures or densities. The diffusion coefficient obtained in this work is much larger than that (1.0– 2.1⫻10⫺4 cm2 s⫺1兲 calculated from ab initio

Z. G. ZHU AND C. S. LIU

9324

FIG. 2. Bond-angle distribution function g 3 ( ␪ ) at five different densities. The cutoff distance is 2.9 Å. Inset: the relative probability P(60°) of bond angle 60° 共⫾7°兲 found in the system at five different densities.

simulations19,20 and other semiempirical or empirical simulations21–24 共see Table I in Ref. 20兲. As can be seen from Table I, the change in density does not lead to a clear change in the diffusion constant. In general, however, the increase in the density gives rise to decrease in the diffusion coefficient. We think it is attributed to the small change in the density and the large pre-exponential of the diffusion coefficient 共see below兲. From Table I it can be also seen that the diffusion coefficient increases greatly with increasing temperature. Figure 3 shows the natural logarithm of diffusion coefficients versus the reciprocal of temperatures, along with a curve fit obtained assuming Arrhenius behavior for the liquid. The fit is quite satisfactory and yields a value of activation energy of 0.92 eV and a value of the preexponential of 30.8⫻10⫺3 cm2 s⫺1. The present activation energy may be compared with the Stillinger-Weber values of 0.56 eV, however, the pre-exponential is much larger than the Stillinger-Weber values of 3.5⫻10⫺3 cm2 s⫺1. 25,26 The very large pre-exponential and the large activation energy suggest that the interatomic forces obtained from Tersoff potential are stronger than those of the actual system, as the well-known high melting point of Tersoff potential does, and the same suggestion has been proposed in some papers.13,27

PRB 61

FIG. 3. Temperature dependence of the diffusion coefficient, the dotted line is a curve fit obtained assuming Arrhenius behavior.

The fairly large pre-exponential may contribute to why the diffusion coefficient obtained in this work is much larger than those calculated from ab initio simulations and other empirical simulations, and also may be responsible for why the increase in the density does not result in a clear decrease in the diffusion constant. We have calculated the activation energies or the average energy barriers at six different temperatures by determining the difference between the average local saddle energy and the local minimum energy from 1200 successful ART events at a certain temperature, which are presented in Fig. 4. As can be seen from this figure, the average energy barrier 共1.03 eV兲 at 2800 K is much larger than those 共0.2–0.3 eV兲 at other five relative high temperatures. In order to corroborate this behavior of the average energy barrier, we have performed another 1200 successful ART events at 2800, 2850, 2900, and 2950 K, respectively. The obtained average energy barriers are also presented in Fig. 4. From this figure, we can clear see that the average energy barrier shows a much steep increase at low temperatures with decreasing temperature. There exists the difference between the average energy barriers at some temperatures 共obtained from the ART兲 and the activation energy 共fitted from the diffusion coefficients兲. Note that 共i兲 the magnitude of the average energy barrier is

TABLE I. Diffusion coefficients for liquid silicon at different temperatures and different densities. Temperature 共K兲

3000

Density 共g/cm3兲 2.53 2.55 2.57 2.59 2.61 Diffusion coefficient (10⫺4 cm2 s⫺1) 8.8 8.3 9.0 9.5 8.7

2800 3000 3200 3400 3600 3800 6.5

9.0

2.57 10.5 13.1 15.5 18.3

PRB 61

MOLECULAR-DYNAMICS SIMULATION OF THE . . .

9325

crease with decreasing temperature. So the viscosity and the surface tension increase on the entire trend with the decrease in temperature; in the meantime, the increases of density enhance these changes in the viscosity and the surface tension. The steep increase of the viscosity at low temperature may be attributed to the abnormal increase of the average energy barrier and the size of the ART events. It should be stressed that the activation energy and the size of the ART events exhibit a much steep increase with the temperature, whereas the similar characteristic of the local structures is not observed. More extensive study is essential. Since any global arrangement of the particles depends on the distribution of the clusters, the global arrangement of the atoms in liquid silicon must depend on the temperature. It should be important to find how the global arrangement of the atoms in liquid silicon change with the temperature in future studies. IV. CONCLUSIONS

FIG. 4. Activation energy obtained by the activation-relaxation technique versus the temperature.

related to the choice of a Metropolis accept-reject criterion; 共ii兲 the activation energy is fitted assuming Arrhenius behavior for the liquid, and it is an averaged result for six temperatures. We also examine the size of an event, which is defined by the number of atoms that move more than a certain threshold displace. The number of atoms displaced from a minimum to a saddle point or from the saddle point to the new minimum shows the similar temperature dependence as the activation energy. Based on the present obtained and previous results, we try to make the following analysis of some physical properties. The anomalous increase or decrease in the bonds contributed to the 52° peak and the covalent bonds with the temperature gives rise to the anomalous change in the conduction electrons, so the anomalous behavior of electrical resistivity takes place. However, it does not lead to the anomalous behavior of viscosity in the corresponding temperature range because the average energy barrier or activation energy almost does not change. The dynamics in the materials, to a first approximation, is determined by the activation energy, i.e., the energy needed to bring a configuration from a local minimum to a nearby saddle point. Except for the anomalous behavior of those above bonds, on the whole trend of change, those bonds that correspond to long bonds and contribute to the 52° and 60° peaks decrease with decreasing temperature, while those bonds with covalent characters in-

1

Y. Waseda, The Structure of Non-Crystalline Materials: Liquid and Amorphous Solids 共McGraw-Hill, New York, 1980兲. 2 Y. Waseda and K. Suziki, Z. Phys. B 20, 339 共1975兲. 3 J. P. Gabathuler and S. Steeb, Z. Naturforsch. A 34, 1314 共1979兲. 4 V. M. Glazov, S. N. Chizhevskaya, and N. N. Glagoleva, Liquid

We have applied molecular-dynamics methods and the activation-relaxation technique to investigate the inherent structure and diffusion properties of liquid silicon. The bonds contributed to the 52° and 60° peaks in the bond-angle distribution functions decrease with increasing density; while the height of the spread main peak that is mainly attributed to the covalent bonds increases and its location moves towards the tetrahedral angle. With the change of temperature, the diffusion coefficients obtained from the average mean-square displacement can be fitted with Arrhenius equation, which yields a value of activation energy of 0.92 eV and a value of a value of the pre-exponential of 30.8⫻10⫺3 cm2 s⫺1. In contrast with this, the activation energy determined from the activation-relaxation technique is in the range of 0.22–1.0 eV and shows a steep increase at the low temperature, and so does the size of the successful activation-relaxation events. The very large pre-exponential and the large activation energy suggest that the interatomic forces obtained from Tersoff potential are stronger than those of the actual system. The steep increase of the activation energy and the size of the activation-relaxation events at low temperature may contribute to the anomalous behavior of the viscosity. ACKNOWLEDGMENTS

We are greatly indebted to Normand Mousseau for sending the code of the activation-relaxation technique and his help with this technique. C. S. Liu thanks D. Y. Sun and Junchao Xia for their useful discussions. This work was supported by the National Natural Science of China 共Grant No. 19874067兲 and the Foundation of Chinese Academy of Science 共Grant No. KJ952-J1-412兲.

Semiconductors 共Plenum, New York, 1980兲. H. Sasaki, E. Tokizaki, K. Terashima, and S. Kimura, Jpn. J. Appl. Phys., Part 1 33, 3803 共1994兲; H. Sasaki, E. Tokizaki, S. Takeda, K. Terashima, and S. Kimura, ibid. 34, 482 共1995兲. 6 H. Sasaki, E. Tokizaki, X. Huang, K. Terashima, and S. Kimura, 5

9326

Z. G. ZHU AND C. S. LIU

Jpn. J. Appl. Phys., Part 1 34, 3432 共1995兲. H. Sasaki, Y. Anzai, X. Huang, K. Terashima, and S. Kimura, Jpn. J. Appl. Phys., Part 1 34, 414 共1995兲. 8 H. Sasaki, A. Ikari, K. Terashima, and S. Kimura, Jpn. J. Appl. Phys., Part 1 34, 3426 共1995兲. 9 Y. Waseda, K. Shinoda, K. Sugiyama, S. Takeda, K. Terashima, and J. M. Toguri, Jpn. J. Appl. Phys., Part 1 34, 4124 共1995兲. 10 H. Balamane, T. Halicioglu, and W. A. Tiller, Phys. Rev. B 46, 2250 共1992兲, and references therein. 11 J. Tersoff, Phys. Rev. B 38, 9902 共1988兲; 39, 5566 共1989兲. 12 M. Ishimaru, K. Yoshida, and T. Motooka, Phys. Rev. B 53, 7176 共1996兲. 13 M. Ishimaru, K. Yoshida, T. Kumamoto, and T. Motooka, Phys. Rev. B 54, 4638 共1996兲. 14 C. S. Liu, Z. G. Zhu, J. Xia, and D. Y. Sun, Phys. Rev. B 60, 3194 共1999兲. 15 G. T. Barkema and N. Mousseau, Phys. Rev. Lett. 77, 4358 共1996兲; 81, 1865 共1998兲. 16 N. Mousseau and L. J. Lewis, Phys. Rev. Lett. 78, 1484 共1997兲; Phys. Rev. B 56, 9461 共1997兲. 7

PRB 61

F. H. Stillinger and T. A. Weber, Phys. Rev. A 25, 978 共1982兲; 28, 2408 共1983兲; F. H. Stillinger, Science 267, 1935 共1995兲. 18 N. Mousseau and G. T. Barkema, Phys. Rev. E 57, 2419 共1998兲. 19 I. Stich, R. Car, and M. Parrinelo, Phys. Rev. Lett. 63, 2240 共1989兲; Phys. Rev. B 44, 4262 共1991兲. 20 J. R. Chelikowsky, N. Troullier, and N. Binggeli, Phys. Rev. B 49, 114 共1994兲. 21 G. Servalli and L. Colombo, Europhys. Lett. 22, 107 共1993兲. 22 C. Z. Wang, C. T. Chan, and K. M. Ho, Phys. Rev. Lett. 66, 189 共1991兲. 23 R. Virkkumen, K. Laasonen, and R. M. Nieminen, J. Phys.: Condens. Matter 3, 7455 共1991兲. 24 P. Allen and J. Broughton, J. Phys. Chem. 91, 4964 共1987兲. 25 M. H. Grabow, G. H. Gilmer, and A. F. Bakker, in Atomic Scale Calculations in Materials Science, edited by J. Tersoff, D. Vanderbilt, and V. Vitek, MRS Symposia Proceedings No. 141 共Materials Research Society, Pittsburgh, 1989兲, p. 349. 26 S. J. Cook and P. Clancy, Phys. Rev. B 47, 7686 共1993兲. 27 K. Moriguchi and A. Shintani, Jpn. J. Appl. Phys., Part 1 37, 414 共1998兲. 17