System-size convergence of point defect properties: The case of the silicon vacancy

We present a comprehensive study of the vacancy in bulk silicon in all its charge states from $2+$ to $2-$, using a supercell approach within plane-wave density-functional theory, and systematically quantify the various contributions to the well-known finite size errors associated with calculating formation energies and stable charge state transition levels of isolated defects with periodic boundary conditions. Furthermore, we find that transition levels converge faster with respect to supercell size when only the $\Gamma$-point is sampled in the Brillouin zone, as opposed to a dense k-point sampling. This arises from the fact that defect level at the $\Gamma$-point quickly converges to a fixed value which correctly describes the bonding at the defect center. Our calculated transition levels with 1000-atom supercells and $\Gamma$-point only sampling are in good agreement with available experimental results. We also demonstrate two simple and accurate approaches for calculating the valence band offsets that are required for computing formation energies of charged defects, one based on a potential averaging scheme and the other using maximally-localized Wannier functions (MLWFs). Finally, we show that MLWFs provide a clear description of the nature of the electronic bonding at the defect center that verifies the canonical Watkins model.

Figure 1

Potential alignment correction for ${\mathrm{v}}^{2-}$ using the MLWF and real-space Voronoi cell methods. The values are adjusted with respect to the bulk supercell. The dashed horizontal lines show the value for the furthest point from the defect center.

Figure 2

Dispersion of the unrelaxed defect level. The symbols show the position of the level at $\Gamma$, and the bars show the extent of the dispersion obtained from our dense k-point sampling. The positions are given with respect to the bulk VBM.

Figure 3

Relaxation effects for the neutral vacancy. Squares show the defect formation energy for the unrelaxed lattice and circles for the relaxed one. Labeled below is the point group symmetry of the defect center after relaxation. The last column shows the relaxed formation energy from previous studies: (a) Ref. 6 ($N=216$, $k_{\mathrm{MP}}=2$, LDA); (b) Ref. 3 ($N=216$, $k_{\mathrm{MP}}=2$, GGA); (c) Ref. 4 ($N=256$, $k_{\mathrm{MP}}=2$, GGA); (d) Ref. 7 ($N=1000$, $k_{\mathrm{MP}}=3$, LDA); (e) Ref. 7 ($N=1000$, $k_{\mathrm{MP}}=3$, GGA). Our calculations were performed with an ultrasoft (norm-conserving) pseudopotential and the LDA functional, with a converged plane-wave energy cutoff of 400 eV (800 eV) and a k-point grid of $k_{\mathrm{MP}}=3$ ($k_{\mathrm{MP}}=4$).

Figure 4

Formation energy of the different charge states of the vacancy as a function of the electron chemical potential (plotted relative to the VBM). The energy range shown covers the DFT band gap for silicon. The thermodynamically stable charge state at each point is highlighted in bold, and the circles indicate the level position for the stable transitions. The dashed vertical lines show the available experimental values for the transition levels. The results shown are for the 256-atom bcc supercell with a $\Gamma$-point sampling.

Figure 5

Watkins’ LCAO model for the silicon vacancy, deduced from electron paramagnetic resonance (EPR) studies of the defect center (Ref. 36), and later confirmed by electron-nuclear double resonance (ENDOR) measurements (Refs. 37 and 38). The four orbitals associated with the neighboring atoms of the defect site are the dangling bonds resulting from the removal of the central silicon atom. The figure shows the predicted ionic configuration and point group symmetry for different charge states of the vacancy due to Jahn-Teller distortion (Ref. 39). For the ${\mathrm{D}}_{2d}$ and ${\mathrm{C}}_{2v}$ configurations, the first charge state listed in the figure refers to the spin-polarized case with only one electron in the highest occupied orbital, and the second charge state refers to the non-spin-polarized case with two electrons.

Figure 6

Contour-surface plots of the MLWFs most strongly associated with the defect center for different charge states of the vacancy, calculated using the 256-atom bcc supercell and a norm-conserving pseudopotential; for computational efficiency, we only include the $\Gamma$-point wave function in the wannierization procedure, since we expect a negligible difference in the qualitative features of the resulting Wannier functions. Each separate closed surface is an individual Wannier function. The amplitude $A$ of the contour shown in each figure is given in the caption in terms of the primitive cell volume $v$. In general, the Wannier functions also have small negative components on the backbonds (not pictured). $f$ is the electronic occupancy of a single orbital. The wannierization procedure is performed using the wannier90 (Ref. 24) code (version 1.2).