Defects in silicon: the role of vibrational entropy

Solid State Communications 128 (2003) 181–185 www.elsevier.com/locate/ssc

Defects in silicon: the role of vibrational entropy M. Sanati, S.K. Estreicher* Physics Department, Texas Tech University, Lubbock, TX 79409-1051, USA Received 7 August 2003; accepted 8 August 2003 by M. Cardona

Abstract Vibrational free energies are calculated from first-principles in the same Si periodic supercells routinely used to perform defect calculations. The specific heat, vibrational entropy, and zero-point energy obtained in defect-free cells are very close to the measured values. The importance of the vibrational part of the free energy is studied in the case of two defect problems: the relative energies of the H2 and Hp2 dimers and the binding energy of a copper pair. In both cases, the vibrational entropy term causes total energy differences to change by about 0.2 eV between 0 and 800 K. We also comment on the rotational entropy in the case of H2 and the configurational entropy in the case of the Cu pair. These examples illustrate the importance of extending first-principles calculations of defects in semiconductors to include free energy contributions. q 2003 Elsevier Ltd. All rights reserved. PACS: 61.72.Bb; 63.20.Mt; 71.55.Cn Keywords: A. Semiconductors; C. Impurities in semiconductors; C. Point defects; D. Thermodynamic properties

First-principles calculations of properties of defects in semiconductors have evolved tremendously in the past decades [1]. Today’s most common approach involves periodic supercells to represent the host crystal (64 host atoms is a typical size) and first-principles densityfunctional theory for the electronic structure. The total electronic energy Eelec: includes the nuclear repulsion energies (see e.g. [2,3]) and is calculated accurately. The results include reliable defect geometries, increasingly accurate (within 2 – 5% of experiment) local vibrational modes (LVMs), and the formation, binding, and various activation energies are typically within a few tenths of an eV of the measured values. The uncertainties in energy differences are associated with the size of the cell, basis set, k-point sampling, choice of exchange-correlation potential, total zero-point energy (ZPE) contributions, and other factors which depend on the defect under study and/or the user. Molecular-dynamic simulations allow real-time simulations of fast processes at finite temperatures. * Corresponding author. Tel.: þ1-806-7423723; fax: þ 1-8067421181. E-mail address: [email protected] (S.K. Estreicher). 0038-1098/$ - see front matter q 2003 Elsevier Ltd. All rights reserved. doi:10.1016/j.ssc.2003.08.005

However, this temperature is defined from the kinetic energy of the nuclei and the free energy is not included. However, even at T ¼ 0 K; the total energy is Eelec: plus the free energy. The Gibbs free energy G (constant pressure) is most relevant to experimentalists but we focus here on the Helmholtz free energy F (constant volume). Working at constant volume rather than constant pressure is appropriate up to several hundred degrees Celsius in Si because the thermal expansion coefficient is very small, 4.68 £ 1026 K21 at room temperature (RT), and the phonon frequencies shift very slowly with T: At RT in Si, the difference between the constant-pressure and constantvolume specific heats [4] CP 2 CV is 0.0165 J/mol K, a correction of only 0.08% to CP ¼ 20 J=mol K: Thus, in the case of semiconductors such as Si, calculating the vibrational spectra at T ¼ 0 K and ignoring the temperature-dependence of the lattice constant is justified and greatly simplifies the calculations. The Helmholtz free energy F ¼ U 2 TS contains the internal energy U and the entropy S; which may include electronic, configurational, vibrational, and rotational contributions. In insulators and moderately-doped semiconductors, the electronic entropy can be ignored up to high

182

M. Sanati, S.K. Estreicher / Solid State Communications 128 (2003) 181–185

temperatures since few electrons populate the excited states. The configurational entropy is proportional to the logarithm of the number of possible configurations and must be calculated for each defect separately. We show here that accurate vibrational free energies can be calculated from first-principles in the same periodic supercells commonly used for electronic structure calculations at T ¼ 0 K: The contribution of the rotational and configurational entropies are discussed. Since Einstein’s pioneering calculation of the specific heat of elemental solids [5], many authors have calculated phonon dispersion curves, phonon free energies, and thermodynamic properties of defect-free crystals. This includes first-principles studies similar to the ones we use here [6 –9]. As concerning defects, empirical calculations of free energies have been performed [10], and the formation free energy of the Si vacancy has been obtained from firstprinciples by thermodynamic integration [11]. The approach proposed here involves first-principles calculations of the entire dynamical matrix of a periodic supercell, with or without a defect. The phonon density of states (DoS) (including any LVMs) is then used to calculated the vibrational part of the free energy. Our calculations are based on the SIESTA [2,3] method (Spanish Initiative for Electronic Simulations with Thousands of Atoms), which expands on Sankey’s ‘ab initio tight-binding’ ideas [12,13]. The basis sets for the valence states are linear combinations of (numerical) atomic orbitals. Double-zeta basis sets (with d orbitals for Cu) are used here. The (local) exchange-correlation potential is that of Ceperley – Alder [14] as parameterized by Perdew – Zunger [15]. Norm-conserving pseudopotentials in the Kleiman – Bylander form [16] remove the core regions from the calculations. The charge density is projected on a real-space grid with an equivalent cutoff of 150Ry to calculate the exchange-correlation and Hartree potentials. The k-point sampling is limited to k ¼ 0: The (harmonic) dynamical matrices (DMs) are calculated from derivatives of the density matrix relative to nuclear coordinates using linear response theory [17,18]. This approach has been used to predict [18] accurate LVMs for light impurities and identify ‘pseudo’-LVMs [19] in Si. In the Si64 and Si128 cells, the G phonon is calculated at 540 and 535 cm21, respectively. This should be compared to the low-temperature harmonic [20] Si frequency, 531 cm21. The DM gives all the normal-mode frequencies of the cell which in turn leads to the phonon DoS gðvÞ: In perfectcrystal studies, the DoS is normally calculated in a primitive cell with many k-points to sample the Brillouin zone. However, defect calculations require the use of the largest possible cell and therefore very few k-points: In this paper, we present evidence that accurate perfect-crystal ZPE, specific heat, and vibrational entropy can be obtained from the DM of the Si64 cell with k ¼ 0: The method is expanded to calculate the temperature dependence of the total energy difference between the interstitial H2 molecule [21] and the

Hp2 complex [22] and the binding energy of the substitutional – interstitial {Cus,Cui} pair [19] in Si. The vibrational part of the Helmholtz free energy (of the cell) is given by Refs. [17,23] ð1 Fvib: ¼ kB T ln{sinhðx=2Þ}gðvÞdv ð1Þ 0

where x ¼ "v=kB T kB Ðis the Boltzmann constant, and if N is the number of atoms, 1 0 gðvÞ; dv ¼ 3N is the total number of modes. In the perfect cell, the integral extends up to the G phonon and gðv . vG Þ ¼ 0: In a cell containing a defect, gðvÞ has continuous and discrete parts g0 ðvÞ þ SLVM dðv 2 vLVM Þ where the sum includes all the modes above the G 0 phonon. The perturbed quantities are primed. Once Fvib: is calculated, the vibrational entropy Svib: and the vibrational contribution to the specific heat at constant volume CV are obtained from standard thermodynamics relations. Note that Fvib: reduces to the total ZPE when T ¼ 0 K: In the harmonic approximation, !   ›Fvib: ›2 Fvib: Svib: ¼ 2 ; CV ¼ 2T : ð2Þ ›T V ›T 2 V If one starts with an accurate DM, the set of 192 (or 384) normal-mode frequencies calculated in the Si64 (or Si128) cell suffices to produce ZPE, Svib: and CV in good agreement with experiment. Even if one takes gðvÞ to be a simple sum of d-functions over all the discrete frequencies, one obtains values within 5% of the measured ones. For example, the ZPE, Svib: at 300 K and CV at 800 K are 6.3 kJ/mol, 17.9 J/mol K, and 24.0 J/mol K, respectively, as compared to the measured 6.0, 18.8, 25.6, respectively. However, better results can be obtained by averaging the set of discrete frequencies to generate a continuous gðvÞ: The averaging process must stop when nearing the calculated G phonon, which must remain at its initial (pre-averaging) calculated value. One can average the frequency spectrum using Lorentzians or Gaussians and vary their width to produce thermodynamic quantities more accurate than those obtained with d-functions; or sharper phonon DoS plots. We obtain satisfactory results with a two-step process, with typical results within 2% of the experimental data. First, we build a histogram by averaging, at each wavenumber, the number of modes within a window Dv: The window width is of the order of a typical error bar on the calculated modes, 30 cm21. The histogram is then smoothed using Gaussians with 7 cm21 width. The results depend very little on the details of the averaging, which is highly desirable since we repeat the same process when dealing with defects. For example, if we vary Dv by a factor of 2 (from 30 to 60 cm21), the calculated Debye temperature, entropy (at 300 K), and specific heat (at 300 K) vary from 668.6 to 667.3 K, 18.20 to 18.40 J/mol K, and 19.47 to 19.52 J/mol K, respectively. The results are similarly independent of the Gaussian width. The key is to start with a set of accurate normal-mode frequencies, especially at the high-frequency end.

M. Sanati, S.K. Estreicher / Solid State Communications 128 (2003) 181–185

Finally, the results obtained in Si128 are within 1% of the Si64 ones. The computational effort required to calculate DMs increases substantially with cell size. Therefore, the numbers cited in the tables and figures below are all from the Si64 cell. Table 1 compares the measured and calculated ZPE, that is Fvib: ðT ¼ 0Þ; as well as vibrational entropy and specific heat at 300 and 800 K. The calculated specific heat and vibrational entropy for the perfect cell are compared to experimental data in Fig. 1. Within linear response theory, the normal modes are purely harmonic and CV ðT ! 1Þ ¼ 3kB : In Fig. 1, the experimental data begin to deviate from theory around 800 or 900 K because of the growing contributions of anharmonicity (and electronic entropy). In the case of Si, all the calculated quantities are in very good agreement with the measured ones at least up to 800 K. Thus, accurate free energies, specific heats, and vibrational entropies can indeed be calculated from first principles, in real space, using the DM matrix of the Si64 cell. These calculations are easily extended to defect problems. First, we focus on the relative energy of the H2 and Hp2 dimers. H2 is the interstitial hydrogen molecule, which is not bound to the crystal. In Hp2, one Si – Si bond is replaced by two Si – H bonds along the same trigonal axis. Hp2 perturbs the phonon spectrum of Si much more than H2, and both defects have high-energy LVMs. Second, we calculate the binding energy of the substitutional/interstitial copper pair {Cus,Cui}. Even though Cu is heavier than Si and no LVM exists, each of the Cui, Cus, and {Cus,Cui} defects perturbs the phonon spectrum of the host crystal in different ways, leading to different vibrational energies as T increases. In all cases, the total ZPE contributions are fully included. Table 2 lists the contributions to Etot ¼ Eelec: þ Fvib: and DEtot at T ¼ 0; 300, and 800 K. The difference in total energy between the two dimers is plotted as a function of T in Fig. 2. The energy difference in favor of H2 increases with T: As shown in Table 2, most of this change is caused by the vibrational entropy term. The interstitial H2 molecule has been seen by Raman [25] as well as infra-red (IR) absorption spectroscopy [26, 27]. The Raman studies involved plasma exposure at about 200 8C and high concentrations of H in a thin layer. This Table 1 Calculated (in Si64) and measured (Expt.) ZPE (kJ/mol), vibrational entropy Svib: (J/mol K), and specific heat C (measured: CP ; calculated: CV ; in J/mol K), per atom

T (K) Calc. Expt. a b

Fig. 1. Calculated (black line) and measured (þs: Ref. [24], £ s: Ref. [4]) specific heat (top) and vibrational entropy (bottom) per atom vs. temperature. The calculated specific heat is CV and the measured one is CP :

situation is not discussed here. The samples in the IR work were hydrogenated by exposure to a gas near the melting point of Si then quenched. They contain H concentrations of about 1016 cm23. It has been shown [28] that hydrogen penetrates the samples in atomic form and that pairing occurs during the quench. Interstitial H2 molecules and {B,H} pairs are observed, but Hp2 is never seen [29]. Since its IR intensity is much larger than that of H2, even trace amounts of Hp2 should be visible in the IR spectra. The increase in the energy difference between H2 and Hp2 as a function of T provides a straightforward explanation. Indeed, with DE ¼ 0:285 eV at room T; the probability (Boltzmann factor) of Hp2 formation is less than 0.002%. When comparing the H2 and Hp2 pairs, one must also consider the contribution of the rotational entropy for H2. Molecular-dynamics simulations predicted that interstitial H2 [30] is a nearly free rotator in Si, and this was confirmed experimentally by the explicit observation of the ortho/para splitting by Raman spectroscopy [31]. Therefore, approximating H2 in Si by a free rotator is reasonable when Table 2 Contributions to Etot ¼ Eelec: þ Fvib: and total energy differences (in eV) at various temperatures for H2 and Hp2 in Si T (K) 0

ZPE

Svib:

Svib:

C

C

300

0 6.201 6.008a

300 18.395 18.820a

800 40.308 41.568b

300 19.522 19.992b

800 24.027 25.357b

800

Ref. [4]. Ref. [24].

183

Eelec: Fvib: ¼ ZPE TSvib: Fvib: Etot TSvib: Fvib: Etot

H2

Hp2

DEtot

26880.443 4.517 3.720 2.929 26877.514 21.804 27.463 26887.906

26880.199 4.545 3.685 2.970 26877.229 21.650 27.338 26887.537

20.244 20.272

DEtot , 0 means that H2 is more stable than Hp2.

20.285

20.369

184

M. Sanati, S.K. Estreicher / Solid State Communications 128 (2003) 181–185 Table 3 Contributions to Etot ¼ Eelec: þ Fvib: at 0, 300, and 800 K for the reaction Si63CusCui þ Si64 ! Si63Cus þ Si64Cui T (K) 0 300 800

Fig. 2. Total energy difference between the interstitial H2 molecule and the Hp2 pair in Si as a function of T:

calculating the rotational energy, which is done in a number of statistical physics texts [32]. At 500 K, TSrot: for H2 is about 0.09 eV. This is a small correction to the 9.88 eV contributed by TSvib: at that temperature (see Table 2). However, when considering the total energy difference between H2 and Hp2, the vibrational entropy contributions vastly compensate each other while TSrot: exists only for H2. Thus, in this case, the rotational entropy increases the slope in Fig. 2 by 0.09 eV at 500 K and reduces even further the probability of Hp2 formation in samples hydrogenated at high temperatures and quenched. More details about rotational entropy contributions will be discussed elsewhere [33]. Note that since the number of sites for H2 is half that for Hp2, the difference in configurational entropy between the two pairs is negligible. Our second example deals with the binding energy of the substitutional/interstitial copper pair in Si, {Cus,Cui}. The binding energy [19] Eb has been calculated to be 1.16 eV. This is the difference in electronic energies of the pair vs. isolated Cus and Cui species. This number decreases to 1.15 eV once the difference in total ZPE is included. However, the measured thermal dissociation energy [34] Ed is 1.02 ^ 0.07 eV. Since Ed ¼ Eb þ Em ; where Em is the migration barrier (0.18 ^ 0.02 eV [35]), the calculated value at T ¼ 0 K is not very close to the measured one. However, in the experiments, the dissociation does not occur at T ¼ 0 K but while the samples are annealed at various temperatures. The contributions to the total energy are given in Table 3 at 0, 300, and 800 K. The energy difference DE ¼ Etot {Si63 Cus Cui } þ Etot {Si64 } 2 Etot {Si63 Cus } 2 Etot {Si64 Cui } (where Etot ¼ Eelec: þ Fvib: ) is plotted as a function of T in Fig. 3. The binding energy is Eb ¼ 2DE: The calculations show that at a few hundred degree Celsius, the calculated Eb lies within the experimental range, even when the migration energy Em ¼ 0:18 eV is added to it. In this case, one must also consider the difference in configurational entropy between {Cus,Cui} and its bypro-

Fvib: ¼ ZPE TSvib: Fvib: TSvib: Fvib:

Si63CusCui

Si64

Si63Cus

Si64Cui

3.936 4.039 2.145 22.708 28.816

4.114 3.660 2.549 21.390 27.657

4.035 3.768 2.394 21.724 28.011

4.000 4.008 2.218 22.583 28.668

ducts Cus þ Cui. This is a straightforward calculation in the Si64 cell and gives Sconf: ¼ 22:4 £ 1024 eV=K: This adds 0.1 eV at 500 K and the slope of the curve in Fig. 3 is steeper as a result. More details about configurational entropy contributions will be discussed elsewhere [33]. Our key points can be summarized as follows. (1)

(2)

(3)

Accurate vibrational free energies (and therefore ZPEs, vibrational entropies, and specific heats) can be calculated from first-principles in real space in Si64 periodic supercells with k ¼ 0; if the complete DM of the cell is computed. Above 800 or 900 K, anharmonic and electronic entropy contributions cause the measured specific heat to begin to deviate from the one calculated in the defect-free cell. Our phonon DoS gðvÞ is obtained from the eigenvalues of the DM of the cell following a two-step averaging process. This improves the accuracy relative to measured quantities from about 5% (no averaging at all) to about 2%. Ongoing calculations in diamond and (wurtzite) GaN cells indicate that the same procedure can be applied hosts other than Si. When a defect is in the cell, the phonon DoS is perturbed. Its continuous part (up to the perturbed G phonon) is constructed in the same way and all the modes above the G phonon are added in a sum of d-functions.

Fig. 3. Binding energy of the {Cus,Cui} pair in Si as a function of T: The measured dissociation energy [34], 1.02 ^ 0.07 eV, is between the two horizontal dashed lines.

M. Sanati, S.K. Estreicher / Solid State Communications 128 (2003) 181–185

(4)

(5) (6) (7)

Applications to the relative energy of H2 and Hp2 and the binding energy of copper pairs show that the vibrational entropy term leads to contributions in total energy differences of the order of 0.20 eV at 800 K. The total energy difference curves start rather flat for the first 50–100 K, then slope away quasi-linearly with T: In some cases, rotational and/or configurational entropy contributions must be included. The total energy variations as a function of temperature depend on the defect in a non-trivial way. They can be larger for defects that have no high-frequency LVMs than for defects that do, as illustrated by the hydrogen pair and copper pair examples. This means that shortcuts (such as focusing only on the contributions of selected LVMs) is likely to produce incorrect results. The entire DM of the cell is needed in order to calculate the free energy of defects.

Acknowledgements This work is supported by the R.A. Welch Foundation, the National Renewable Energy Laboratory, and the Alexander von Humboldt Foundation. Many thanks to Texas Tech’s High Performance Computer Center for generous amounts of CPU time. Thanks are also due to M. Cardona for fruitful discussions and A.M. Stoneham for suggesting the problem.

References [1] S.K. Estreicher, Mater. Today 6 (6) (2003) 26. [2] D. Sa´nchez-Portal, P. Ordejo´n, E. Artacho, J.M. Soler, Int. J. Quantum Chem. 65 (1997) 453. [3] E. Artacho, D. Sa´nchez-Portal, P. Ordejo´n, A. Garcı´a, J.M. Soler, Phys. Status Solidi (b) 215 (1999) 809. [4] P. Flubacher, A.J. Leadbetter, J.A. Morrison, Philos. Mag. 4 (1959) 273 In this paper, ‘cal/g atom’ should read ‘cal/mol’ or be divided by atomic mass of Si. [5] A. Einstein, Ann. Phys. (Leipzig) 19 (1906) 289. [6] P. Giannozzi, S. deGironcoli, P. Pavone, S. Baroni, Phys. Rev. B 43 (1991) 7231. [7] O. Sugino, R. Car, Phys. Rev. Lett. 74 (1995) 1823.

185

[8] C. Lee, X. Gonze, Phys. Rev. B 56 (1997) 7321. [9] G. Kern, G. Kresse, J. Hafne, Phys. Rev. B 59 (1999) 8551. [10] R. LeSar, R. Najafabadi, D.J. Srolovitz, Phys. Rev. Lett. 63 (1989) 624. [11] E. Smargiassi, R. Car, Phys. Rev. B 53 (1996) 9760. [12] O.F. Sankey, D.J. Niklevski, Phys. Rev. B 40 (1989) 3979. O.F. Sankey, D.J. Niklevski, D.A. Drabold, J.D. Dow, Phys. Rev. B 41 (1990) 12750. [13] A.A. Demkov, J. Ortega, O.F. Sankey, M.P. Grumbach, Phys. Rev. B 52 (1995) 1618. [14] D.M. Ceperley, B.J. Alder, Phys. Rev. Lett. 45 (1980) 566. [15] S. Perdew, A. Zunger, Phys. Rev. B 32 (1981) 5048. [16] L. Kleiman, D.M. Bylander, Phys. Rev. Lett. 48 (1982) 1425. [17] X. Gonze, C. Lee, Phys. Rev. B 55 (1997) 10355. [18] J.M. Pruneda, S.K. Estreicher, J. Junquera, J. Ferrer, P. Ordejo´n, Phys. Rev. B 65 (2002) 075210. [19] S.K. Estreicher, D. West, J. Goss, S. Knack, J. Weber, Phys. Rev. Lett. 90 (2003) 035504. [20] F. Widulle, T. Ruf, M. Konuma, I. Silier, M. Cardona, W. Kriegseis, V.I. Ozhogin, Solid State Commun. 118 (2001) 1. [21] S.K. Estreicher, Acta Physica Polonica A 102 (2002) 403. ¨ berg, [22] J.D. Holbech, B. Bech Nielsen, R. Jones, P. Sitch, S. O Phys. Rev. Lett. 71 (1993) 875. [23] A.A. Maradudin, E.W. Montroll, G.H. Weiss, I.P. Ipatova, Theory of Lattice Dynamics in the Harmonic Approximation, Academic Press, New York, 1971. [24] I. Barin, Thermochemical data of pure substances, VCH, Weinheim, 1995. [25] A.W.R. Leitch, V. Alex, J. Weber, Phys. Rev. Lett. 81 (1998) 421. [26] R.E. Pritchard, M.J. Ashwin, J.H. Tucker, R.C. Newman, Phys. Rev. B 57 (1998) R15048. [27] J.A. Zhou, M. Stavola, Phys. Rev. Lett. 83 (1999) 1351. [28] R.C. Newman, R.E. Pritchard, J.H. Tucker, E.C. Lightowlers, Physica B 273–274 (1999) 164. [29] R.C. Newman, private communication. [30] S.K. Estreicher, K. Wells, P.A. Fedders, P. Ordejo´n, J. Phys.: Condens. Matter 13 (2001) 6271. [31] E.V. Lavrov, J. Weber, Phys. Rev. Lett. 89 (2002) 215501. [32] L.D. Landau, E.M. Lifshitz, Statistical Physics, 3rd ed., Pergamon Press, New York, 1980, Part 1, p. 137 ff. [33] M. Sanati, S.K. Estreicher, unpublished. [34] A.A. Istratov, H. Hieslmair, T. Heiser, C. Flink, E.R. Weber, Appl. Phys. Lett. 72 (1998) 474. [35] A.A. Istratov, C. Flink, H. Hieslmair, E.R. Weber, T. Heiser, Phys. Rev. Lett. 81 (1998) 1243.