ab inito local vibrational modes of light impurities in silicon

We have developed a formulation of density functional perturbation theory for the calculation of vibrational frequencies in molecules and solids, which uses numerical atomic orbitals as a basis set for the electronic states. The (harmonic) dynamical matrix is extracted directly from the first order change in the density matrix with respect to infinitesimal atomic displacements from the equilibrium configuration. We have applied this method to study the vibrational properties of a number of hydrogen-related complexes and light impurities in silicon. The diagonalization of the dynamical matrix provides the vibrational modes and frequencies, including the local vibrational modes (LVMs) associated with the defects. In addition to tests on simple molecules, results for interstitial hydrogen, hydrogen dimers, vacancy-hydrogen and self-interstitial-hydrogen complexes, the boron-hydrogen pair, substitutional C, and several O-related defects in c-Si are presented. The average error relative to experiment for the aprox.60 predicted LVMs is about 2% with most highly harmonic modes being extremely close and the more anharmonic ones within 5-6% of the measured values.


I. INTRODUCTION
The knowledge of the structures of impurities and defects is an essential prerequisite for understanding the electrical and optical changes that these complexes induce in semiconductors such as crystalline silicon. 1,2 The presence of light impurities such as H, B, C, or O results in the appearance of infrared ͑IR͒ or Raman-active local vibrational modes ͑LVM's͒, usually well isolated from the frequency range of the phonons of the host material. The observation of LVM's, coupled with isotope substitutions and uniaxial stress measurements, provide precious information about the type and number of impurity atoms involved and the symmetry of the defect. However, these data are rarely sufficient to identify the defect unambiguously.
Since the early days of Stein, 3 a large number of vibrational modes have been identified through the interplay of experiment and theory. The calculation of LVM's at the ab initio level provides a critical link between theory and experiment. This is particularly true in the case of hydrogen, since it binds covalently in the immediate vicinity of many impurities and defects, thus giving rise to a number of LVM's in the range ϳ800 to ϳ2200 cm Ϫ1 . Other common impurities in Si which produce LVM's are B, C, and O, but any element lighter than Si can in principle be observed by LVM spectroscopy.
The computation of systematically accurate vibrational frequencies is a challenge for first-principles theory, given their sensitivity on the details of bonding geometry and electronic structure. Various approaches have been used to calculate LVM's, from semiempirical models, 4 to ab initio Hartree-Fock 5 and density-functional theory. [6][7][8] Typical accuracies in the calculated vibrational modes for light impurities in silicon in former works are within 3% and 10% of the experimental data, which means in some cases a deviation of over 100 cm Ϫ1 .
In most cases, frequencies are calculated in the spirit of the frozen phonon approximation. One computes the total energy of the system in the equilibrium configuration ͑that in which the forces acting on the atoms are zero͒, and then for small displacements of selected atoms ͑either individually or in the direction of a normal mode, if this is known͒. The actual value of the atomic displacements ͑typically a few hundredths of an Å͒ are parameters chosen by the user. One can either fit the energy vs displacement to a polynomial, and extract a specific vibrational mode, 6,9 or compute the dynamical matrix by finite differences. 7 When a few specific modes are all that is needed, only the movement of the atoms involved in those modes is considered. In these methods, it is not possible to isolate the harmonic contributions completely from the anharmonic ones, since finite displacements always involve some anharmonic effects. For this reason, the frequencies obtained in this approach are sometimes referred to as quasiharmonic. 10 One can also calculate vibrational properties from constant-temperature molecular-dynamics simulations, for instance by extracting selected frequencies from the velocityvelocity autocorrelation function 11 or by using more sophisticated spectral estimators, like the multiple signal classification ͑MUSIC͒ algorithm. 12,13 This is computationally exhausting, since long molecular dynamics runs are required, but potentially very accurate. 13 This also allows the calculation of frequencies as a function of temperature.
However, a calculation of vibrational frequencies does not necessarily require the actual displacement of the atoms, as in the methods described above. Linear-response theory ͑in particular through the application of perturbation theory in density-functional theory͒ was thoroughly used in the past 14 to compute the response of the system to infinitesimal atomic displacements, and from that, the vibrational frequencies in the harmonic approximation. This can be done with only a knowledge of the electronic solution in the equilibrium con-figuration. The advantage of this approach is that anharmonic effects are eliminated, and that no reference is needed to explicit finite atomic displacements. In addition, this approach allows one to compute phonons with an arbitrary q vector in crystalline systems, without having to consider a supercell commensurate with the periodicity of the phonon, as is required in the frozen-phonon and molecular-dynamics approaches.
We prosent a method, based on density-functional perturbation theory ͑DFPT͒, to compute vibrational frequencies in the harmonic approximation. We use a basis set of numerical atomic orbitals to expand the electronic wave functions, which makes the method computationally very efficient, and allows us to calculate systems with a large number of atoms. We then apply it to make a systematic study of a number of defect centers in silicon, involving light impurities and their complexes with intrinsic defects ͑vacancies and self-interstitials͒. In most cases, a comparison of the calculated and measured vibrational frequencies is very favorable, improving on the results obtained by other approaches.
The outline of this paper is as follows. We first discuss the theoretical method and the model used to describe the defects. Then we compare the vibrational properties obtained for a variety of complexes with experimental data as well as other first-principles calculations in the literature. Finally, we discuss the results.

A. Ground-state description
In this work, we use the fully self-consistent ab initio code SIESTA. 15,16 The electronic energy is obtained from density-functional theory ͑DFT͒ 17,18 within the local-density approximation. The exchange-correlation potential is that of Ceperley and Alder 19 as parametrized by Perdew and Zunger. 20 Norm-conserving pseudopotentials 21 in the Kleinman-Bylander form 22 are used to remove the core electrons from the calculations.
The valence-electron wave functions are described with numerical linear combinations of atomic orbitals of the Sankey type, 23 but generalized to be arbitrarily complete with the inclusion of multiple-zeta orbitals and polarization states. 24 These orbitals are numerical solutions of a free atom with the appropriate pseudopotential, and are strictly zero beyond some cutoff radius. This makes the calculation of the ground-state Kohn-Sham Hamiltonian scale linearly with the number of atoms, 15,16 allowing calculations in very large systems with a modest computational cost.
In the present work, the basis sets include single-zeta ͑SZ͒, double-zeta ͑DZ͒, double-zeta plus polarization ͑DZP͒, and triple-zeta plus polarization ͑TZP͒. A DZP basis includes two sets of s and p's plus one set of d's on Si, O, C or B, and two s's and one set of p's on H. The radial cutoff of the atomic orbitals was determined as described in Ref. 24, with an energy shift of 0.5 eV, and a split norm of 0.15 for all the species except H, for which the split norm was 0.5. The charge density is projected on a real space grid with an equivalent cutoff of 90-150 Ry to calculate the exchangecorrelation and Hartree potentials. The host crystal is repre-sented by a periodic supercell of 64 host atoms, and the k-point sampling is reduced to the ⌫ point. This restriction appears to be quite sufficient for a calculation of vibrational spectra. Tests have been performed for selected defects in a 128-host-atom cell, and the results are within a few wave numbers of those obtained in the 64 host atoms cell.
In order to determine the equilibrium structure of the defects studied, we have relaxed all the atomic coordinates with a conjugate gradient algorithm, reaching a tolerance in the forces of F max Ͻ0.01 eV/Å . The dynamical matrix for the whole cell is computed ͑see below͒ from this ground state, and its eigenfrequencies and eigenmodes are obtained.

B. Linear-response theory
An implementation of DFPT was developed to compute the electronic response to infinitesimal atomic displacements. As is well known from the ''2nϩ1'' theorem in quantum mechanics, 25 the first-order change of the electronic wave function in a perturbative expansion allows the computation of the second-order change in the energies. This implies that only the properties of the unperturbed ground state are needed to obtain the linear response of the system. The extension of this theorem to DFT is that a knowledge of the first-order change in the electronic density variationally determines the second-order change in the energy. In this way, we analytically obtain the dynamical matrix from the gradient of the density relative to atomic displacement and, from this, the vibrational properties, without physically displacing any atom.
Here we briefly describe the key points of our formulation. In the Appendix, additional technicalities are presented. A complete report will be published elsewhere. 26 The change in the electronic wave function is obtained by solving the first-order perturbation expansion of the Schrödinger equation ͑the Sternheimer equation 27 where 0,i are the ground-state electronic wave functions and ␦ i the first order perturbation of i when an atom is displaced ͑if we consider atom ␣, this would be ‫ץ‬ ␣ i ). As we expand these wave functions in terms of atomic orbitals centered at R , the derivatives can be directly written in terms of the derivatives of the atomic orbitals:

͑3͒
Only the orbitals centered on atom ␣ appear in the last term which is equal to Ϫc i •ٌ (rϪR ). The change in the coefficients, ‫ץ‬ ␣ c i , is obtained from Eq. ͑1͒. ‫ץ‬ i is then used to compute the perturbation in the electronic density This allows a computation of the dynamical matrix, by an explicit derivation of the forces on all the atoms ␤ in the system ͑the expressions of which can be found in the Appendix, for our approach͒ with respect to the infinitesimal displacement of one of them (␣): We remark that only terms up to first order in the electronic wave functions appear in the resulting formulas. Note that only linear effects in the wave function are obtained in this method, which is consistent with the harmonic approximation implicitly assumed in the diagonalization of the dynamical matrix. Thus one expects to obtain high-quality frequencies for the vibrational modes that are harmonic, but the frequencies of modes involving large anharmonic contributions will be less accurate. Although phonons with an arbitrary q vector can be obtained in the present approach, here we only calculate vibrations which are periodic with the simulation supercell ͑i.e., qϭ⌫͒, since we are interested in LVM's.

III. RESULTS
Tests of the method for free SiH 4 , CO, CO 2 , and H 2 lead to very good agreement with experiment ͑see Table I͒. Using an appropriate basis set is essential to reproduce accurate frequencies. In most of the cases a DZ set gives good values, but in some configurations a more complete basis is required. This is particularly true for bending modes. Oxygen ''likes'' to have polarization orbitals, and thus the frequencies are better when these are included. Note that the vibrational frequencies are more sensitive to the basis set size than the structural properties such as bond lengths. In general, the largest improvements in the frequencies correspond to going from SZ to DZ, then DZ to DZP, but the TZ basis produce only marginal improvements.
A number of defects containing light impurities in c-Si are now considered. These are H at a bond-center ͑BC͒ site, H dimers (H 2 and H 2 *), the hydrogenated vacancy (VH n , nϭ1,2,3,4) and the saturated divacancy V 2 H 6 , the selfinterstitial-hydrogen IH 2 complex, the ͕B,H͖ pair, substitutional C, interstitial O, and two charge states of the A-center ͑oxygen-vacancy complex͒. These embrace a range of Si-X bonding configurations, with XϭH, C, O and B. In most cases, we have tried different basis sets in order to check or improve the accuracy of our calculations when the chemical properties are particularly complex. We show frequencies for the DZ and DZP basis sets. Typically, the calculations of the whole dynamical matrix for each of the systems considered here take around one week of CPU in a 733-MHz Pentium-III processor. Although the diagonalization of the dynamical matrix gives all the (⌫ point͒ modes in the cell, we present here only the stretching and some bending LVM's. We also obtain the eigenvectors which we use to determine the symmetry of the corresponding vibrational modes. Our results are compared with available experimental data, and with theoretical frequencies obtained by other authors using firstprinciples DFT.

A. H BC , H 2 , and H 2 *
At the BC site, 29 hydrogen exists in the ϩ1 and 0 charge states. In the latter case, the odd electron does not participate in the bonding but resides in a nonbonding orbital primarily localized on the two Si atoms adjacent to the proton. Chemically, this is a three-center, two-electron bond, very much like the type of H bonding occurring in boron hydrides. The bond is somewhat compressed because optimal relaxation ͑no Si second-nearest neighbors͒ would likely result in a longer Si-H-Si bond and a frequency lower than observed. This suggests that as H moves along the trigonal axis, it tends to form something like Si . . . H-Si, then Si-H-Si, then Si-H . . . Si, a process which is highly anharmonic. The calculated frequency is given in Table II.
Raman 40 and IR 33,41 measurements of H 2 in silicon reveal a considerable softening of the stretching mode with respect to the frequency of H 2 in the gas phase. A number of calculations ͑for a review, see Ref. 39͒ found the molecule to be stable at the tetrahedral interstitial ͑T͒ site. The electron affinity of the Si atoms surrounding H 2 is at least partly responsible for a small charge transfer from H 2 to its Si neighbors, which results in a weakening of the H-H bond. Even though the H-H stretch mode is not expected to be fully harmonic, our calculated frequency is close to the experimental one ͑Table II͒. Note that the errors relative to experiment in the D 2 and HD frequencies are very different than the error in the H 2 frequency. This is also a clear feature of these frequencies when they are calculated dynamically from the v-v autocorrelation function. 42 The trigonal H 2 * defect 8 consists of one hydrogen atom close to the antibonding ͑AB͒ site, and the other near the BC site ͑see Fig. 1͒. The two H atoms are inequivalent. The Si-H AB bond length is slightly longer than the Si-H BC one ͑we obtain 1.580 and 1.510 Å, respectively͒ which gives rise to different stretch frequencies for the two atoms: 2135 cm Ϫ1 for H BC and 1750 cm Ϫ1 for H AB . We also obtain two degenerate wagging modes, associated with the H BC atom, at 560 cm Ϫ1 , and two wagging modes with 839 and 843 cm Ϫ1 , related to H AB . The latter should be degenerate, but small inaccuracies in the atomic relaxations render them at slightly different frequencies.

B. Hydrogen and native defects
A considerable number of IR and Raman lines are related to H-intrinsic defect complexes. It has been noted 43,44 that vibrational modes above 2000 cm Ϫ1 are mainly related to H in vacancies, while those lines below 2000 cm Ϫ1 are predominant for H-self-interstitial systems or H at AB sites. A large number of geometrical configurations may lead to very similar vibrational lines, making it difficult to identify these defects unambiguously. As noted by other groups, 2,45-47 the stretching frequencies in VH n (nϭ1,2,3,4) systems increases with n due to the repulsive H-H interaction. Thus, the highest IR line is that of VH 4 . Si-H bonds point toward the center of the vacancy along the trigonal axes.
In our calculations ͑Table III͒, VH has monoclinic symmetry, and the H oscillates parallel to the ͗111͘ direction. In orthorhombic VH 2 , two equivalent H's have stretching modes along the ͗100͘ and ͗001͘ directions. The frequencies for these modes are 2121 and 2144 cm Ϫ1 respectively. VH 3 has C 3v symmetry. The A singlet involves the movement of three H atoms toward the vacancy, while in the twofolddegenerate E mode one of the atoms moves in opposition. VH 4 has T d symmetry. In addition to the threefolddegenerate T 2 mode at 2205, we obtain an IR-inactive singlet A 1 mode at 2265 cm Ϫ1 .
The vibrational modes of V 2 H 6 are almost identical to those of VH 3 : The fully saturated divacancy behaves very much like two weakly coupled VH 3 complexes. The A 2 singlet at 2176 cm Ϫ1 induces a dipole along the ͗111͘ direction. In addition to this mode and the IR-active E doublet, we obtain two IR-inactive modes at 2186 and 2134 cm Ϫ1 .
The IH 2 complex 50 has two equivalent and weakly coupled hydrogen atoms, which yields two very similar stretching frequencies. Uniaxial stress measurements show that the two hydrogen atoms are equivalent. Our relaxed structure has an almost C 2v symmetry, with the A mode higher than the B mode, confirming early results. 50 The deviation from perfect symmetry is due to the finite tolerance  in the geometry optimization. This deviation is seen when comparing the IHD and IDH complexes: they should be identical, but we find their frequencies to be off by 2 cm Ϫ1 . Note that we reproduce the correct ordering for the bending modes of IH 2 at 737 and 732 cm Ϫ1 ͑Table IV͒.

C. Oxygen, carbon, and boron in silicon
Oxygen is a well-known impurity which is especially important in Czochralski-grown Si, and a considerable amount of effort was focused in understanding its properties. 51 We have computed the LVM frequencies for interstitial oxygen (O i ) and two charge states of the vacancy-oxygen complex (A center͒. The results are in Table V. Frequencies for O i were computed with the oxygen placed at the BC site, where the probability of finding this delocalized atom is at a maximum, and the classical harmonic potential can better describe the local modes. The IR-active A 2u mode corresponds to the asymmetric-stretching mode, while the A 1g is the symmetric one. The E u mode involves the movement of nearest silicon atoms with no participation of the oxygen. 52 Finally, Table VI shows the triply degenerate mode of substitutional carbon as well as the LVM's associated with a ͕B,H͖ complex. 2

IV. CONCLUSIONS
We have presented a development of DFPT, using localized atomic wave functions as a basis set, and applied it to the study of LVM's for light impurities in silicon. In contrast to other methods, the dynamical matrix is computed analytically without actually displacing any atom from its equilibrium position. The calculations are based on the ground-state density matrix, as computed with the SIESTA package.
Tests of the method for free molecules (SiH 4 , H 2 , CO, and CO 2 ) show that this approach is highly accurate in situations where the anharmonic contributions are small. Note that the frequencies are obtained at Tϭ0 K, while experimental data are obtained at low, but nonzero, temperatures.
We have used a variety of basis set sizes to describe the The average error relative to experiment of the ϳ60 calculated modes is about 2%. In situations where large anharmonic contributions are present, the accuracy of the method decreases somewhat ͑5-6 %͒. This occurs, for example, when H is close to a BC position. However, in most cases the calculated frequencies are in remarkable agreement ͑0-2 %͒ with experimental data, implying that this perturbative approach is totally justified and that the ground-state density matrix calculated with SIESTA is very reliable.
where n i is the occupation of state i, and , loop over the orbitals in the basis set. When we consider a perturbed system, we can express the change in the electronic density in a similar way, defining a density matrix ‫ץ‬ ␣ :

͑A2͒
The matrix ‫ץ‬ ␣ involves derivatives of the coefficients c i , which are obtained from the solution of the Eq. ͑1͒. Its final expression 26 is

͑A3͒
where ‫ץ‬ ␣ H and ‫ץ‬ ␣ S denote the change in the Hamiltonian and overlap matrices H ϭ͗ ͉Ĥ ͉ ͘ and S ϭ͗ ͉ ͘. Once the change in the density matrix ‫ץ‬ ␣ is known, we derive the expression of the forces in SIESTA with respect to atomic displacement ͑Ref. 16͒: Here H 0 stands for the non-self-consistent contribution of the Hamiltonian ͑kinetic plus nonlocal pseudopotential terms͒; the matrix E ϭ ͚ i ⑀ i c i * c i is defined in a similar way to , V NA is the screened neutral-atom potential; ␦(r) is the difference (r)Ϫ atom (r), where atom is a sum of localized atomic densities; V H ␦ is the Hartree potential of ␦; and V xc is the exchange-correlation potential.
Note that, in ‫ץ‬ ␣ F ␤ ϭϪ‫ץ‬ 2 E/‫ץ‬ ␣ ‫ץ‬ ␤ , just first-order derivatives of will appear, and that the knowledge of ‫ץ‬ ␣ would enable the computation of a whole row of the dynamical matrix. More details of the computation of these matrix elements will be published elsewhere. 26 *Corresponding author: Email address: pruneda@icmab.es