- Email: [email protected]
Ab initio molecular dynamics simulations on structural change of supercooled liquid Si at different temperatures from 1700 to 1100 K
Physica B 406 (2011) 3991–3996
Contents lists available at ScienceDirect
Physica B journal homepage: www.elsevier.com/locate/physb
Review
Ab initio molecular dynamics simulations on structural change of supercooled liquid Si at different temperatures from 1700 to 1100 K Yubing Wang a,b,n, Gang Zhao c, Changsong Liu b, Zhengang Zhu b a
Institute of Intelligent Machines, Chinese Academy of Sciences, Hefei 230031, PR China Key Laboratory of Materials Physics, Institute of Solid State Physics, Chinese Academy of Sciences, Post Office 1129, Hefei 230031, PR China c Department of Physics and Electronic Engineering, Ludong University, Hongqi Road, No. 186, Yantai 264025, PR China b
a r t i c l e i n f o
abstract
Article history: Received 2 February 2010 Received in revised form 30 July 2011 Accepted 31 July 2011 Available online 9 August 2011
Using ab initio molecular dynamics simulations, the local atomic structure and electronic properties of supercooled liquid Si (l-Si) at different temperatures from 1700 to 1100 K were studied. Our calculated coordination numbers present no obvious change in the temperature range investigated. Our results indicate that the structure of supercooled l-Si may be well described as a combined local atomic configuration of white-tin and diamond type structures. Upon cooling from 1700 to 1100 K, the tetrahedral white-tin type ordering collapses gradually toward the tetrahedral diamond-type structure. No drastic change behavior is observed in our work. & 2011 Elsevier B.V. All rights reserved.
Keywords: Temperature-induced structural change White-tin type structure Diamond type structure
Contents 1. 2. 3. 4.
Introduction . . . . . . . . . . Computational methods Results and discussion . . Conclusions . . . . . . . . . . Acknowledgments . . . . . References . . . . . . . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
1. Introduction Silicon is one of the most intensively investigated materials for many years due to its great technological importance. Silicon single crystals widely used for ultralarge-scale integrated circuit devices mainly are grown from melts by the Czochralski method, and rapid development of technology requires upgrading of the quality and scaling up of the size of Si devices. Therefore the fundamental physical properties, including the atomic-scale structure of liquid Silicon (l-Si), are very crucial for understanding and controlling the crystal-growth processes [1]. l-Si has several intriguing properties but its structure is still ambiguous particularly in the supercooled region, which is
n Corresponding author at: Key Laboratory of Materials Physics, Institute of Solid State Physics, Chinese Academy of Sciences, Post Office 1129, Hefei 230031, PR China. Tel.: þ 86 551 5591429; fax: þ86 661 5591434. E-mail address: [email protected] (Y. Wang).
0921-4526/$ - see front matter & 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.physb.2011.07.062
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. 3991 . 3992 . 3992 . 3995 . 3995 . 3996
important for the growth of a single crystal. Upon melting the crystal structure with the coordination number of 4 develops into a more compact liquid structure characterized by a coordination number exceeding 6, which suggests that the covalent-bonding character still remains in the liquid state [2–4]. Moreover, l-Si is metallic while the covalent-bonding crystal and amorphous phase are semiconducting [2–4]. In the supercooled liquid state, the controversy mainly focuses on the existence of a first-order liquid–liquid (L–L) phase transition. Since experimental studies on l-Si have been hampered in deeply supercooled regions, several computer simulations have been performed simultaneously. For example x-ray measurements and ab initio molecular dynamics simulations by N. Jakse et al. have been performed to study the atomic structure of l-Si for a deep supercooling of 230 K below the melting point (1685 K) and show conclusively that the coordination number decreases as the liquid is supercooled, indicative of an L–L phase transition [5,6]. Recently they reported new results of first principles molecular dynamics simulations,
3992
Y. Wang et al. / Physica B 406 (2011) 3991–3996
which show the increase of tetrahedral local structures upon cooling and confirm their early speculations of a liquid–liquid phase transition [7]. By using orbital-free ab initio molecular dynamics simulations of Si at different temperatures from 1550 to 1100 K, M. Colakogullari et al. have advanced similar viewpoints [8]. A molecular-dynamics study by Sastry and Angell also demonstrated that the Stillinger–Weber (SW) potential yields a discontinuous L–L transition at around 1060 K [9]. Kimura et al. have employed electromagnetic levitation to obtain a supercooling of 290 K below the melting point; however, the results show no discontinuous behavior in this temperature region [10]. Meanwhile, using the technique of electrostatic levitation together with high-energy x-ray diffraction, Kim et al. obtained high quality structural data more deeply into the supercooled regime 316 K below the melting point [11]. In their results, no change of coordination number is observed and the existence of the L–L phase transition is in doubt [11]. The discrepancies between these results described above show that none of the groups found convincing evidence of a drastic change, which left the question of an L–L transition on supercooled l-Si still unsolved [12,13]. To obtain microscopic atomic structure and electronic properties of supercooled l-Si and investigate possible structural changes, we carried out ab initio molecular dynamics simulations on deeply supercooled region from 1700 to 1100 K. The paper is organized as follows: in Section 2, we describe the method of our simulations; the calculated results and the related discussions are presented in Section 3; finally a short summary is given in Section 4.
2. Computational methods Our simulations were performed with the framework of the density function theory (DFT) [14]. We used the Vienna ab initio simulation package (VASP) [15], and employed ultrasoft pseudopotential (USPP) of the Vanderbilt type with the generalized gradient approximation (GGA) formulated by Perdew and Wang (PW91) to the exchange-correlation energy [16–18]. The 80-atom system was put in a simple cubic box with periodic boundary conditions. Only the G point was used to sample the Birllouin zone of the supercell. The electronic wave functions were expanded in the plane-wave basis set, with an energy cutoff of 250 eV. Our canonical ensemble simulations have been performed at 1100, 1300, 1500 and 1700 K with a Nose´ thermostat to control temperature [19]. The experimental densities determined by Watanabe et al. were used at each temperature, as listed in Table 1 [20]. The Verlet algorithm was used to integrate Newton’s equations of motion and the time step of ion motion was 3 fs [21]. The Kohn–Sham energy function was minimized by the preconditioned conjugate-gradient method [22]. The initial atomic configuration adopted was a random distribution of 80 atoms on the grid, which was constructed by dividing the supercell into 5 5 5 square segments. Then the system was heated up to 2000 K by rescaling the ionic velocities. After equilibration for 6 ps at this temperature, we gradually reduced the temperature to 1700 K. For other temperatures, we repeated this procedure and only changed the final temperature in 1500, 1300 and 1100 K. At each temperature, the physical Table 1 Experimental densities used in our simulations [18] and the calculated coordination numbers at four temperatures. T (K)
r0 (at/A˚ 3)
1100 0.05665
1300 0.05645
1500 0.05616
1700 0.05569
N
5.93 (4.84)
6.08 (4.64)
6.01 (4.6)
5.96 (4.54)
properties of interest were obtained by averaging over 6 ps after an equilibration taking 3 ps.
3. Results and discussion Firstly, the mean square displacement (MSD) was calculated to investigate the information of microscopic atomic motion and analyze the character of four supercooled liquid states. The timedependent mean-square displacement (MSD) of supercooled l-Si is defined in the following way: * + 2 N 1 X ! ! 2 /DrðtÞ S ¼ ð1Þ r i ðt þt0 Þ r i ðt0 Þ , N i¼1 ! where the summation goes over all N atoms of the supercell, r i gives the coordinates of atom i, t0 is an arbitrary origin of time, and the angular brackets denote an average over all possible time origins. In the liquid state the behavior of the MSD in the limit of long time is linear with time, which is correlated to the diffusion coefficient D by the so-called Einstein relationship ! /D r ðtÞ2 S-6Dt þ C,
ð2Þ
where D is the diffusion coefficient and C is a constant. The MSD for Si atoms at four different temperatures from 1100 to 1700 K are shown in Fig. 1(a). MSD always increases with time and this is consistent with the generalized Einstein model in which the atoms in a cage, formed by its surrounding atoms, oscillate about a center, which itself is undergoing Brownian motion. This confirms well that our sample is in liquid state. Using Eq. (2), we have calculated the diffusion coefficient D from the slopes of the fitting lines in Fig. 1(a). In Fig. 1(b) we can find that the natural logarithm of D approximately has a linear relationship with the reciprocal of temperature (1/T). It means the temperature-dependent diffusion coefficient of Si atom follows the Arrhenius relationship D ¼ D0 expðED =kB TÞ,
ð3Þ
where ED is called the activation energy, D0 is the preexponential factor, and kB is the Boltzmann constant. By a linear fit, we can calculate that ED and D0 are 0.25 eV and 2.83 10 4 cm2 s 1, respectively. We can notice that the calculated values of D are smaller than the corresponding values in Ref. [23], indicating that Si atoms in the present work are relatively less diffusive in liquid state and this also confirms that our samples are supercooled liquids. The pair correlation function, g(r), plays a central role in the physics of liquids, because it is directly measurable and various properties can be estimated from it when coupled with an appropriate theory. This quantity is used to describe threedimensional spatial distribution of particles from one-dimensional point of view, and is proportional to the probability of finding a particle at a distance r from a reference particle. We can calculate g(r) from the atomic coordinates as * + XX ! ! 1 gðrÞ ¼ dð r r ij Þ , ð4Þ r0 N i j a i where N is the total number of atoms, r0 is the average number ! ! ! ! ! density of atoms, and r ij ¼ r j r i where r i and r j are the position vectors of atoms i and j, respectively. Using the atomic coordinates in the molecular dynamics simulations, the pair correlation functions calculated at four temperatures are shown in Fig. 2. At 1100 K the pair correlation function is characterized by a low peak maximum, which follows the usual first peak corresponding to the nearest neighbors, located at around 0.245 and 0.39 nm. With increasing temperature, the first peak goes lower and broader while the second low peak decreases obviously
Y. Wang et al. / Physica B 406 (2011) 3991–3996
3993
Fig. 1. (a) Variation of the mean square displacement (MSD) with time t for a Si atom in supercooled liquid Si and (b) the diffusion coefficient DSi (cm2/s) of a Si atom in supercooled liquid Si as a function of temperature T.
where r0 is the average number density of atoms and ! ! ! ! ! Q ¼ 9 Q 9 ¼ 9 k k0 9 with k and k0 the wave vector of incident and scattered waves, respectively. Since the calculation of S(Q) by Eq. (5) on a small simulation box may lead to significant errors, a direct calculation from the atomic coordinates in the simulations was also performed [2,25,26]: * + ! ! 1 XX expði Q U r ij Þ N dQ ,0 , ð6Þ SðQ Þ ¼ N i jai
Fig. 2. Pair correlation functions of supercooled liquid Si at four temperatures.
and almost disappears finally. As the positions of the first and second peaks of the pair-distribution function in the diamond type structure are 0.235 and 0.384 nm, respectively, the shortrange structure at 1100 K is supposed to be characterized by this structure as in the solid state [24]. The white-tin type structure at high temperatures is more densely packed, in which the first and second nearest neighbors reside on 0.253 and 0.266 nm far away, respectively [24]. The positions of the first two nearest neighbors are very close so that these two peaks may merge into one broad peak, so we may infer that the white-tin type structure exists at high temperatures. Then we investigate the temperature dependence of the atomic structure in the reciprocal space. In one-component liquids, the structure factor S(Q) can be obtained by the Fourier transformation of pair-correlation function g(r): SðQ Þ ¼ 1 þ 4pr0
Z
½ gðrÞ1
sin Qr 2 r dr, Qr
ð5Þ
! where N is the total number of atoms and r ij is the interatomic distance between i and j. The calculated results at four temperatures and the experimental data at 1557 K (dash line) are shown in Fig. 3. The calculated results by the two methods are in good agreement, indicating that errors by the Fourier transformation in our simulation box are negligible. The comparison between the simulations and experiments proves to be satisfactory and suggests that our results are acceptable on the whole. In our results, a large temperature dependence of the structure factors can be observed in this range. At 1700 K, one can easily find two peaks: the first one (P1) is quite broad ranging from 26 to 34 nm 1, and the second peak (P2) is at around 56 nm 1. Upon cooling, a shoulder on the high-Q side of the first peak grows in intensity while its position shifts towards the right slightly. Finally, it tends to be more like two decomposed small peaks (P1a and P1b) at 1100 K, at around 26.1 and 34.2 nm 1, respectively. The second peak (P2) seems insensitive to the variation of temperature. With reference to Ref. [26] using the reverse Monte Carlo (RMC) method Watanabe et al. have assumed the modified white-tin structure to calculate S(Q), which correctly reproduced the shoulder on the right side of the first peak. The white-tin type atomic structure, in which the tetrahedra are compressed in one direction and elongated along the other two axes so that each atom has four nearest neighbors and the other two slightly far away, has been proposed by some groups to model the short-range of l-Si [27,28], and can also be evidenced in the present work at high temperatures. Given the pair correlation functions we have calculated the coordination number N, i.e., the average number of near-neighbor atoms around a Si atom, according to Z Rcutoff N¼ 4pr 2 r0 gðrÞdr, ð7Þ 0
3994
Y. Wang et al. / Physica B 406 (2011) 3991–3996
Fig. 4. Bond angle distribution function g3(y) at four different temperatures. The cutoff distance is 0.300 nm. The two dotted lines denote the location of peaks.
Fig. 3. Structure factors of supercooled liquid Si calculated by the Fourier transformation (solid line) and direct calculation (gray line and circle) at four different temperatures, together with the experimental data at 1557 K (dash line).
where r0 is the number density of the atom and Rcutoff is the position of the first minimum of the pair correlation functions, ˚ The calculated results are listed in here chosen to be 3.00 A. Table 1, from which one can easily find that our results show constant coordination numbers of approximately 6 over the entire temperature region. We notice that the first minimum at 1700 K is poorly defined, so we also calculate the coordination number by the symmetric option according to Z Rmax N¼2 4pr 2 r0 gðrÞ dr, ð8Þ 0
where Rmax is the position of the first maximum of the pair ˚ The calculated correlation functions, here chosen to be 2.50 A. results by symmetric option are also given in Table 1, as listed in parentheses, which are also insensitive to the variation of temperature but are smaller than the results calculated from the first minimum method. The temperature independence in our calculated coordination number is in conflict with the results of several experimental and other simulation studies based on the SW potential [5,6,11], suggesting that the existence of a liquid– liquid phase transitions is not substantiated in our work. Additional information on the local short-range order in supercooled l-Si can be obtained from the triplet correlation. This is usually measured by the bond-angle distribution function g3(y), in which y is the angle formed by a pair of vectors drawn from a central atom to any other two atoms within a sphere of cutoff radius Rcutoff. Our calculated g3(y) with the cutoff distance of ˚ shown in Fig. 4, is rather broad with two peaks at around 3.00 A, 561 and 961. Upon cooling from 1700 to 1300 K, these two peaks are almost independent of the change of temperature; however, as temperature is lowered to 1100 K, the two peaks both display obvious changes: the 561 peak at 1100 K goes lower; the 961 one
goes higher and the trough between these two peaks goes deeper. We have also calculated bond-angle distribution functions for coordination number of the central atom of 4, 5, 6, 7, and 8, in order to further investigate the origin of the peaks, which are shown in Fig. 5. For the fourfold coordination the bond angles are intensively distributed around the tetrahedral angle, indicating that the atomic configurations constructed by sp3 bonding is still retained in the supercooled liquid state. For the fivefold coordination the main peak moves towards the small angle, located around 981, and a new peak at around 591 appears. For the sixfold and sevenfold coordinations the 591 peak grows and the 981 peak goes weaker, in which the latter peak moves to the small angle. The 591 peak is considerably strong and the 981 peak evolves to have a broad shoulder in the eightfold coordination. It should be noted that peaks in the bond angle distribution with 4, 5, 6, and 7 nearest neighbors at 1100 K present different properties compared with those at other temperatures: the first peak is weaker while the second peak is stronger and narrower. The bond angles in white-tin type structure associated with four nearest and two next-nearest neighbors are 74.61, 94.01, 105.41, 149.31, and 1801 [27,28]. Our calculated bond angle distribution functions with six coordination numbers present small peaks at around the second broad peak (941) at higher temperatures, which may also serve as a position proof of the white-tin type structure. The microscopic atomic structure is correlated with the electronic structure. Here we have studied the electronic density of state (DOS), i.e., the DOS is composed of angular-momentumresolved contributions. By projecting all the wavefunctions in a sphere of radius R around atom i onto the spherical harmonic (l, m), we obtained the (l, m) angular momentum component of the atom i. The calculated DOS of supercooled l-Si is shown in Fig. 6, from which we can observe an obvious dip at the Fermi level EF and upon cooling the dip becomes deeper. Based on the above-described structural parameters in the present work, we assume that these structural peculiarities may originate from a portion of crystalline Si existing in the supercooled liquid state. We could find geometrical characteristic evidence of both normal form and high-pressure modification solid Si, which have diamond short-range structure and white-tin type structure, respectively. At higher temperatures, the structure of supercooled l-Si may be well described as a combined local atomic configuration with the characteristic features of white-tin ordering predominating over diamond-type ordering. Upon cooling, the fundamental white-tin structure is continuously
Y. Wang et al. / Physica B 406 (2011) 3991–3996
3995
Fig. 5. Bond angle distribution functions g3(y) for atoms with different coordination numbers at four different temperatures. The cutoff distance is 0.300 nm. The dotted lines denote the location of peaks.
previous claims of a liquid–liquid phase transition is found in our simulations.
4. Conclusions
Fig. 6. Total densities of states for supercooled liquid Si at four different temperatures.
deformed to the diamond structure; i.e., the diamond local atomic configuration is enhanced. Our calculated coordination numbers present no obvious change in the temperature region investigated, which suggests that no structural evidence of
In summary, ab initio molecular dynamics simulations were performed on supercooled liquid Si by using the cell size of 80 atoms with decreasing temperature ranging from 1700 to 1100 K. The dynamical, structural, and electronic properties and their temperature dependence have been studied. The mean square displacement, structural factors, pair correlation functions, coordination number, bond-angle distribution functions, and electronic density of state at each temperature were calculated. By analysis, we conclude that the structure of supercooled l-Si may be well described as a combined local atomic configuration of white-tin and diamond type structures. Upon cooling from 1700 to 1100 K, the tetrahedral white-tin type ordering collapses gradually toward the tetrahedral diamond-type structure. Our calculated coordination numbers show no obvious change in the whole temperature region, which suggests that no discontinuous behavior is observed in our simulations.
Acknowledgments This work is supported by the National Science Foundation of China (Grant nos. 10674135, 10874182, and 50803066), and by
3996
Y. Wang et al. / Physica B 406 (2011) 3991–3996
the Center for Computational Science, Hefei Institutes of Physical Sciences. References [1] S. Kimura, J. Crystallogr. Soc. Jpn. 17 (1990) 264. [2] Y. Waseda, The Structure of Non-Crystalline Materials: Liquid and Amorphous Solids, McGraw-Hill, New York, 1980. [3] Y. Waseda, K. Suziki, Z. Phys. B 20 (1975) 339. [4] J.P. Gabathuler, S. Steeb, Z. Naturforsch. A34 (1979) 1314. [5] S. Ansell, S. Krishnan, J.J. Felten, D.L. Price, J. Phys.: Condens. Matter 10 (1998) L73. [6] N. Jakse, L. Hennet, D.L. Price, S. Krishnan, T. Key, E. Artacho, B. Glorieux, A. Pasturel, M.L. Saboungi, Appl. Phys. Lett. 83 (2003) 4734. [7] N. Jakse, A. Pasturel, Phys. Rev. Lett. 99 (2007) 205702. [8] M. C - olakogullari, S. Dalgic-, L.E. Gonza´lez, D.J. Gona´lez, Eur. Phys. J. Spec. Top. 196 (2011) 45. [9] S. Sastry, C.A. Angell, Nat. Mater. 2 (2003) 739. [10] H. Kimura, M. Watanabe, K. Izumi, T. Hibiya, D.H. Moritz, T. Schenk, K.R. Bauchspieb, S. Schneider, K. Funakoshi, M. Hanfland, Appl. Phys. Lett. 78 (2001) 604.
[11] T.H. Kim, G.W. Lee, B. Sieve, A.K. Gangopadhyay, R.W Hyers, T.J. Rathz, J.R. Rogers, D.S. Robinson, K.F. Kelton, A.I. Goldman, Phys. Rev. Lett. 95 (2005) 085501. [12] N. Jakse, A. Pasturel, J. Chem. Phys. 129 (2008) 104503. [13] S.N. Yannopoulos, J. Chem. Phys. 130 (2009) 247102. [14] W. Kohn, L.J. Sham, Phys. Rev. 140 (1965) A1133. ¨ [15] G. Kresse, J. Furthmuller, Phys. Rev. B 54 (1996) 11169. [16] K. Laasonen, A. Pasquarello, R. Car, C. Lee, D. Vanderbilt, Phys. Rev. B 47 (1993) 10142. [17] D. Vanderbilt, Phys. Rev. B 41 (1990) 7892. [18] J.P. Perdew, Y. Wang, Phys. Rev. B 44 (1991) 13298. [19] S. Nose´, J. Chem. Phys. 81 (1984) 511. [20] M. Watanabe, M. Adachi, T. Morishita, K. Higuchi, H. Kobakazu, H. Fukuyama, Faraday Discuss. Chem. Soc. 136 (2007) 279. [21] L. Verlet, Phys. Rev. 159 (1967) 98. [22] G. Kresse, J. Hafner, Phys. Rev. B 49 (1994) 14251. [23] G.X. Li, C.S. Liu, Z.G. Zhu, Phys. Rev. B 71 (2005) 094209. [24] Y. Waseda, K. Suzuki, Z. Phys. B 20 (1975) 339. [25] G. Zhao, C.S. Liu, Y.N. Wu, E.G. Jia, Z.G. Zhu, Phys. Rev. B 74 (2006) 184202. [26] M. Watanabe, K. Higuchi, A. Mizuno, K. Nagashio, K. Kuribayashi, Y. Katayama, J. Cryst. Growth 294 (2006) 16. [27] H. Ogawa, Y. Waseda, Z. Naturforsch. 49a (1994) 987. [28] V. Petkov, J. Phys.: Condens. Matter 7 (1995) 5745.
Contents lists available at ScienceDirect
Physica B journal homepage: www.elsevier.com/locate/physb
Review
Ab initio molecular dynamics simulations on structural change of supercooled liquid Si at different temperatures from 1700 to 1100 K Yubing Wang a,b,n, Gang Zhao c, Changsong Liu b, Zhengang Zhu b a
Institute of Intelligent Machines, Chinese Academy of Sciences, Hefei 230031, PR China Key Laboratory of Materials Physics, Institute of Solid State Physics, Chinese Academy of Sciences, Post Office 1129, Hefei 230031, PR China c Department of Physics and Electronic Engineering, Ludong University, Hongqi Road, No. 186, Yantai 264025, PR China b
a r t i c l e i n f o
abstract
Article history: Received 2 February 2010 Received in revised form 30 July 2011 Accepted 31 July 2011 Available online 9 August 2011
Using ab initio molecular dynamics simulations, the local atomic structure and electronic properties of supercooled liquid Si (l-Si) at different temperatures from 1700 to 1100 K were studied. Our calculated coordination numbers present no obvious change in the temperature range investigated. Our results indicate that the structure of supercooled l-Si may be well described as a combined local atomic configuration of white-tin and diamond type structures. Upon cooling from 1700 to 1100 K, the tetrahedral white-tin type ordering collapses gradually toward the tetrahedral diamond-type structure. No drastic change behavior is observed in our work. & 2011 Elsevier B.V. All rights reserved.
Keywords: Temperature-induced structural change White-tin type structure Diamond type structure
Contents 1. 2. 3. 4.
Introduction . . . . . . . . . . Computational methods Results and discussion . . Conclusions . . . . . . . . . . Acknowledgments . . . . . References . . . . . . . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
1. Introduction Silicon is one of the most intensively investigated materials for many years due to its great technological importance. Silicon single crystals widely used for ultralarge-scale integrated circuit devices mainly are grown from melts by the Czochralski method, and rapid development of technology requires upgrading of the quality and scaling up of the size of Si devices. Therefore the fundamental physical properties, including the atomic-scale structure of liquid Silicon (l-Si), are very crucial for understanding and controlling the crystal-growth processes [1]. l-Si has several intriguing properties but its structure is still ambiguous particularly in the supercooled region, which is
n Corresponding author at: Key Laboratory of Materials Physics, Institute of Solid State Physics, Chinese Academy of Sciences, Post Office 1129, Hefei 230031, PR China. Tel.: þ 86 551 5591429; fax: þ86 661 5591434. E-mail address: [email protected] (Y. Wang).
0921-4526/$ - see front matter & 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.physb.2011.07.062
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. 3991 . 3992 . 3992 . 3995 . 3995 . 3996
important for the growth of a single crystal. Upon melting the crystal structure with the coordination number of 4 develops into a more compact liquid structure characterized by a coordination number exceeding 6, which suggests that the covalent-bonding character still remains in the liquid state [2–4]. Moreover, l-Si is metallic while the covalent-bonding crystal and amorphous phase are semiconducting [2–4]. In the supercooled liquid state, the controversy mainly focuses on the existence of a first-order liquid–liquid (L–L) phase transition. Since experimental studies on l-Si have been hampered in deeply supercooled regions, several computer simulations have been performed simultaneously. For example x-ray measurements and ab initio molecular dynamics simulations by N. Jakse et al. have been performed to study the atomic structure of l-Si for a deep supercooling of 230 K below the melting point (1685 K) and show conclusively that the coordination number decreases as the liquid is supercooled, indicative of an L–L phase transition [5,6]. Recently they reported new results of first principles molecular dynamics simulations,
3992
Y. Wang et al. / Physica B 406 (2011) 3991–3996
which show the increase of tetrahedral local structures upon cooling and confirm their early speculations of a liquid–liquid phase transition [7]. By using orbital-free ab initio molecular dynamics simulations of Si at different temperatures from 1550 to 1100 K, M. Colakogullari et al. have advanced similar viewpoints [8]. A molecular-dynamics study by Sastry and Angell also demonstrated that the Stillinger–Weber (SW) potential yields a discontinuous L–L transition at around 1060 K [9]. Kimura et al. have employed electromagnetic levitation to obtain a supercooling of 290 K below the melting point; however, the results show no discontinuous behavior in this temperature region [10]. Meanwhile, using the technique of electrostatic levitation together with high-energy x-ray diffraction, Kim et al. obtained high quality structural data more deeply into the supercooled regime 316 K below the melting point [11]. In their results, no change of coordination number is observed and the existence of the L–L phase transition is in doubt [11]. The discrepancies between these results described above show that none of the groups found convincing evidence of a drastic change, which left the question of an L–L transition on supercooled l-Si still unsolved [12,13]. To obtain microscopic atomic structure and electronic properties of supercooled l-Si and investigate possible structural changes, we carried out ab initio molecular dynamics simulations on deeply supercooled region from 1700 to 1100 K. The paper is organized as follows: in Section 2, we describe the method of our simulations; the calculated results and the related discussions are presented in Section 3; finally a short summary is given in Section 4.
2. Computational methods Our simulations were performed with the framework of the density function theory (DFT) [14]. We used the Vienna ab initio simulation package (VASP) [15], and employed ultrasoft pseudopotential (USPP) of the Vanderbilt type with the generalized gradient approximation (GGA) formulated by Perdew and Wang (PW91) to the exchange-correlation energy [16–18]. The 80-atom system was put in a simple cubic box with periodic boundary conditions. Only the G point was used to sample the Birllouin zone of the supercell. The electronic wave functions were expanded in the plane-wave basis set, with an energy cutoff of 250 eV. Our canonical ensemble simulations have been performed at 1100, 1300, 1500 and 1700 K with a Nose´ thermostat to control temperature [19]. The experimental densities determined by Watanabe et al. were used at each temperature, as listed in Table 1 [20]. The Verlet algorithm was used to integrate Newton’s equations of motion and the time step of ion motion was 3 fs [21]. The Kohn–Sham energy function was minimized by the preconditioned conjugate-gradient method [22]. The initial atomic configuration adopted was a random distribution of 80 atoms on the grid, which was constructed by dividing the supercell into 5 5 5 square segments. Then the system was heated up to 2000 K by rescaling the ionic velocities. After equilibration for 6 ps at this temperature, we gradually reduced the temperature to 1700 K. For other temperatures, we repeated this procedure and only changed the final temperature in 1500, 1300 and 1100 K. At each temperature, the physical Table 1 Experimental densities used in our simulations [18] and the calculated coordination numbers at four temperatures. T (K)
r0 (at/A˚ 3)
1100 0.05665
1300 0.05645
1500 0.05616
1700 0.05569
N
5.93 (4.84)
6.08 (4.64)
6.01 (4.6)
5.96 (4.54)
properties of interest were obtained by averaging over 6 ps after an equilibration taking 3 ps.
3. Results and discussion Firstly, the mean square displacement (MSD) was calculated to investigate the information of microscopic atomic motion and analyze the character of four supercooled liquid states. The timedependent mean-square displacement (MSD) of supercooled l-Si is defined in the following way: * + 2 N 1 X ! ! 2 /DrðtÞ S ¼ ð1Þ r i ðt þt0 Þ r i ðt0 Þ , N i¼1 ! where the summation goes over all N atoms of the supercell, r i gives the coordinates of atom i, t0 is an arbitrary origin of time, and the angular brackets denote an average over all possible time origins. In the liquid state the behavior of the MSD in the limit of long time is linear with time, which is correlated to the diffusion coefficient D by the so-called Einstein relationship ! /D r ðtÞ2 S-6Dt þ C,
ð2Þ
where D is the diffusion coefficient and C is a constant. The MSD for Si atoms at four different temperatures from 1100 to 1700 K are shown in Fig. 1(a). MSD always increases with time and this is consistent with the generalized Einstein model in which the atoms in a cage, formed by its surrounding atoms, oscillate about a center, which itself is undergoing Brownian motion. This confirms well that our sample is in liquid state. Using Eq. (2), we have calculated the diffusion coefficient D from the slopes of the fitting lines in Fig. 1(a). In Fig. 1(b) we can find that the natural logarithm of D approximately has a linear relationship with the reciprocal of temperature (1/T). It means the temperature-dependent diffusion coefficient of Si atom follows the Arrhenius relationship D ¼ D0 expðED =kB TÞ,
ð3Þ
where ED is called the activation energy, D0 is the preexponential factor, and kB is the Boltzmann constant. By a linear fit, we can calculate that ED and D0 are 0.25 eV and 2.83 10 4 cm2 s 1, respectively. We can notice that the calculated values of D are smaller than the corresponding values in Ref. [23], indicating that Si atoms in the present work are relatively less diffusive in liquid state and this also confirms that our samples are supercooled liquids. The pair correlation function, g(r), plays a central role in the physics of liquids, because it is directly measurable and various properties can be estimated from it when coupled with an appropriate theory. This quantity is used to describe threedimensional spatial distribution of particles from one-dimensional point of view, and is proportional to the probability of finding a particle at a distance r from a reference particle. We can calculate g(r) from the atomic coordinates as * + XX ! ! 1 gðrÞ ¼ dð r r ij Þ , ð4Þ r0 N i j a i where N is the total number of atoms, r0 is the average number ! ! ! ! ! density of atoms, and r ij ¼ r j r i where r i and r j are the position vectors of atoms i and j, respectively. Using the atomic coordinates in the molecular dynamics simulations, the pair correlation functions calculated at four temperatures are shown in Fig. 2. At 1100 K the pair correlation function is characterized by a low peak maximum, which follows the usual first peak corresponding to the nearest neighbors, located at around 0.245 and 0.39 nm. With increasing temperature, the first peak goes lower and broader while the second low peak decreases obviously
Y. Wang et al. / Physica B 406 (2011) 3991–3996
3993
Fig. 1. (a) Variation of the mean square displacement (MSD) with time t for a Si atom in supercooled liquid Si and (b) the diffusion coefficient DSi (cm2/s) of a Si atom in supercooled liquid Si as a function of temperature T.
where r0 is the average number density of atoms and ! ! ! ! ! Q ¼ 9 Q 9 ¼ 9 k k0 9 with k and k0 the wave vector of incident and scattered waves, respectively. Since the calculation of S(Q) by Eq. (5) on a small simulation box may lead to significant errors, a direct calculation from the atomic coordinates in the simulations was also performed [2,25,26]: * + ! ! 1 XX expði Q U r ij Þ N dQ ,0 , ð6Þ SðQ Þ ¼ N i jai
Fig. 2. Pair correlation functions of supercooled liquid Si at four temperatures.
and almost disappears finally. As the positions of the first and second peaks of the pair-distribution function in the diamond type structure are 0.235 and 0.384 nm, respectively, the shortrange structure at 1100 K is supposed to be characterized by this structure as in the solid state [24]. The white-tin type structure at high temperatures is more densely packed, in which the first and second nearest neighbors reside on 0.253 and 0.266 nm far away, respectively [24]. The positions of the first two nearest neighbors are very close so that these two peaks may merge into one broad peak, so we may infer that the white-tin type structure exists at high temperatures. Then we investigate the temperature dependence of the atomic structure in the reciprocal space. In one-component liquids, the structure factor S(Q) can be obtained by the Fourier transformation of pair-correlation function g(r): SðQ Þ ¼ 1 þ 4pr0
Z
½ gðrÞ1
sin Qr 2 r dr, Qr
ð5Þ
! where N is the total number of atoms and r ij is the interatomic distance between i and j. The calculated results at four temperatures and the experimental data at 1557 K (dash line) are shown in Fig. 3. The calculated results by the two methods are in good agreement, indicating that errors by the Fourier transformation in our simulation box are negligible. The comparison between the simulations and experiments proves to be satisfactory and suggests that our results are acceptable on the whole. In our results, a large temperature dependence of the structure factors can be observed in this range. At 1700 K, one can easily find two peaks: the first one (P1) is quite broad ranging from 26 to 34 nm 1, and the second peak (P2) is at around 56 nm 1. Upon cooling, a shoulder on the high-Q side of the first peak grows in intensity while its position shifts towards the right slightly. Finally, it tends to be more like two decomposed small peaks (P1a and P1b) at 1100 K, at around 26.1 and 34.2 nm 1, respectively. The second peak (P2) seems insensitive to the variation of temperature. With reference to Ref. [26] using the reverse Monte Carlo (RMC) method Watanabe et al. have assumed the modified white-tin structure to calculate S(Q), which correctly reproduced the shoulder on the right side of the first peak. The white-tin type atomic structure, in which the tetrahedra are compressed in one direction and elongated along the other two axes so that each atom has four nearest neighbors and the other two slightly far away, has been proposed by some groups to model the short-range of l-Si [27,28], and can also be evidenced in the present work at high temperatures. Given the pair correlation functions we have calculated the coordination number N, i.e., the average number of near-neighbor atoms around a Si atom, according to Z Rcutoff N¼ 4pr 2 r0 gðrÞdr, ð7Þ 0
3994
Y. Wang et al. / Physica B 406 (2011) 3991–3996
Fig. 4. Bond angle distribution function g3(y) at four different temperatures. The cutoff distance is 0.300 nm. The two dotted lines denote the location of peaks.
Fig. 3. Structure factors of supercooled liquid Si calculated by the Fourier transformation (solid line) and direct calculation (gray line and circle) at four different temperatures, together with the experimental data at 1557 K (dash line).
where r0 is the number density of the atom and Rcutoff is the position of the first minimum of the pair correlation functions, ˚ The calculated results are listed in here chosen to be 3.00 A. Table 1, from which one can easily find that our results show constant coordination numbers of approximately 6 over the entire temperature region. We notice that the first minimum at 1700 K is poorly defined, so we also calculate the coordination number by the symmetric option according to Z Rmax N¼2 4pr 2 r0 gðrÞ dr, ð8Þ 0
where Rmax is the position of the first maximum of the pair ˚ The calculated correlation functions, here chosen to be 2.50 A. results by symmetric option are also given in Table 1, as listed in parentheses, which are also insensitive to the variation of temperature but are smaller than the results calculated from the first minimum method. The temperature independence in our calculated coordination number is in conflict with the results of several experimental and other simulation studies based on the SW potential [5,6,11], suggesting that the existence of a liquid– liquid phase transitions is not substantiated in our work. Additional information on the local short-range order in supercooled l-Si can be obtained from the triplet correlation. This is usually measured by the bond-angle distribution function g3(y), in which y is the angle formed by a pair of vectors drawn from a central atom to any other two atoms within a sphere of cutoff radius Rcutoff. Our calculated g3(y) with the cutoff distance of ˚ shown in Fig. 4, is rather broad with two peaks at around 3.00 A, 561 and 961. Upon cooling from 1700 to 1300 K, these two peaks are almost independent of the change of temperature; however, as temperature is lowered to 1100 K, the two peaks both display obvious changes: the 561 peak at 1100 K goes lower; the 961 one
goes higher and the trough between these two peaks goes deeper. We have also calculated bond-angle distribution functions for coordination number of the central atom of 4, 5, 6, 7, and 8, in order to further investigate the origin of the peaks, which are shown in Fig. 5. For the fourfold coordination the bond angles are intensively distributed around the tetrahedral angle, indicating that the atomic configurations constructed by sp3 bonding is still retained in the supercooled liquid state. For the fivefold coordination the main peak moves towards the small angle, located around 981, and a new peak at around 591 appears. For the sixfold and sevenfold coordinations the 591 peak grows and the 981 peak goes weaker, in which the latter peak moves to the small angle. The 591 peak is considerably strong and the 981 peak evolves to have a broad shoulder in the eightfold coordination. It should be noted that peaks in the bond angle distribution with 4, 5, 6, and 7 nearest neighbors at 1100 K present different properties compared with those at other temperatures: the first peak is weaker while the second peak is stronger and narrower. The bond angles in white-tin type structure associated with four nearest and two next-nearest neighbors are 74.61, 94.01, 105.41, 149.31, and 1801 [27,28]. Our calculated bond angle distribution functions with six coordination numbers present small peaks at around the second broad peak (941) at higher temperatures, which may also serve as a position proof of the white-tin type structure. The microscopic atomic structure is correlated with the electronic structure. Here we have studied the electronic density of state (DOS), i.e., the DOS is composed of angular-momentumresolved contributions. By projecting all the wavefunctions in a sphere of radius R around atom i onto the spherical harmonic (l, m), we obtained the (l, m) angular momentum component of the atom i. The calculated DOS of supercooled l-Si is shown in Fig. 6, from which we can observe an obvious dip at the Fermi level EF and upon cooling the dip becomes deeper. Based on the above-described structural parameters in the present work, we assume that these structural peculiarities may originate from a portion of crystalline Si existing in the supercooled liquid state. We could find geometrical characteristic evidence of both normal form and high-pressure modification solid Si, which have diamond short-range structure and white-tin type structure, respectively. At higher temperatures, the structure of supercooled l-Si may be well described as a combined local atomic configuration with the characteristic features of white-tin ordering predominating over diamond-type ordering. Upon cooling, the fundamental white-tin structure is continuously
Y. Wang et al. / Physica B 406 (2011) 3991–3996
3995
Fig. 5. Bond angle distribution functions g3(y) for atoms with different coordination numbers at four different temperatures. The cutoff distance is 0.300 nm. The dotted lines denote the location of peaks.
previous claims of a liquid–liquid phase transition is found in our simulations.
4. Conclusions
Fig. 6. Total densities of states for supercooled liquid Si at four different temperatures.
deformed to the diamond structure; i.e., the diamond local atomic configuration is enhanced. Our calculated coordination numbers present no obvious change in the temperature region investigated, which suggests that no structural evidence of
In summary, ab initio molecular dynamics simulations were performed on supercooled liquid Si by using the cell size of 80 atoms with decreasing temperature ranging from 1700 to 1100 K. The dynamical, structural, and electronic properties and their temperature dependence have been studied. The mean square displacement, structural factors, pair correlation functions, coordination number, bond-angle distribution functions, and electronic density of state at each temperature were calculated. By analysis, we conclude that the structure of supercooled l-Si may be well described as a combined local atomic configuration of white-tin and diamond type structures. Upon cooling from 1700 to 1100 K, the tetrahedral white-tin type ordering collapses gradually toward the tetrahedral diamond-type structure. Our calculated coordination numbers show no obvious change in the whole temperature region, which suggests that no discontinuous behavior is observed in our simulations.
Acknowledgments This work is supported by the National Science Foundation of China (Grant nos. 10674135, 10874182, and 50803066), and by
3996
Y. Wang et al. / Physica B 406 (2011) 3991–3996
the Center for Computational Science, Hefei Institutes of Physical Sciences. References [1] S. Kimura, J. Crystallogr. Soc. Jpn. 17 (1990) 264. [2] Y. Waseda, The Structure of Non-Crystalline Materials: Liquid and Amorphous Solids, McGraw-Hill, New York, 1980. [3] Y. Waseda, K. Suziki, Z. Phys. B 20 (1975) 339. [4] J.P. Gabathuler, S. Steeb, Z. Naturforsch. A34 (1979) 1314. [5] S. Ansell, S. Krishnan, J.J. Felten, D.L. Price, J. Phys.: Condens. Matter 10 (1998) L73. [6] N. Jakse, L. Hennet, D.L. Price, S. Krishnan, T. Key, E. Artacho, B. Glorieux, A. Pasturel, M.L. Saboungi, Appl. Phys. Lett. 83 (2003) 4734. [7] N. Jakse, A. Pasturel, Phys. Rev. Lett. 99 (2007) 205702. [8] M. C - olakogullari, S. Dalgic-, L.E. Gonza´lez, D.J. Gona´lez, Eur. Phys. J. Spec. Top. 196 (2011) 45. [9] S. Sastry, C.A. Angell, Nat. Mater. 2 (2003) 739. [10] H. Kimura, M. Watanabe, K. Izumi, T. Hibiya, D.H. Moritz, T. Schenk, K.R. Bauchspieb, S. Schneider, K. Funakoshi, M. Hanfland, Appl. Phys. Lett. 78 (2001) 604.
[11] T.H. Kim, G.W. Lee, B. Sieve, A.K. Gangopadhyay, R.W Hyers, T.J. Rathz, J.R. Rogers, D.S. Robinson, K.F. Kelton, A.I. Goldman, Phys. Rev. Lett. 95 (2005) 085501. [12] N. Jakse, A. Pasturel, J. Chem. Phys. 129 (2008) 104503. [13] S.N. Yannopoulos, J. Chem. Phys. 130 (2009) 247102. [14] W. Kohn, L.J. Sham, Phys. Rev. 140 (1965) A1133. ¨ [15] G. Kresse, J. Furthmuller, Phys. Rev. B 54 (1996) 11169. [16] K. Laasonen, A. Pasquarello, R. Car, C. Lee, D. Vanderbilt, Phys. Rev. B 47 (1993) 10142. [17] D. Vanderbilt, Phys. Rev. B 41 (1990) 7892. [18] J.P. Perdew, Y. Wang, Phys. Rev. B 44 (1991) 13298. [19] S. Nose´, J. Chem. Phys. 81 (1984) 511. [20] M. Watanabe, M. Adachi, T. Morishita, K. Higuchi, H. Kobakazu, H. Fukuyama, Faraday Discuss. Chem. Soc. 136 (2007) 279. [21] L. Verlet, Phys. Rev. 159 (1967) 98. [22] G. Kresse, J. Hafner, Phys. Rev. B 49 (1994) 14251. [23] G.X. Li, C.S. Liu, Z.G. Zhu, Phys. Rev. B 71 (2005) 094209. [24] Y. Waseda, K. Suzuki, Z. Phys. B 20 (1975) 339. [25] G. Zhao, C.S. Liu, Y.N. Wu, E.G. Jia, Z.G. Zhu, Phys. Rev. B 74 (2006) 184202. [26] M. Watanabe, K. Higuchi, A. Mizuno, K. Nagashio, K. Kuribayashi, Y. Katayama, J. Cryst. Growth 294 (2006) 16. [27] H. Ogawa, Y. Waseda, Z. Naturforsch. 49a (1994) 987. [28] V. Petkov, J. Phys.: Condens. Matter 7 (1995) 5745.