## 2 Materials and Methods  
All DFT calculations, in this study, were carried out using the Cambridge Sequential Total Energy Package [36]. More specifically, an ultrasoft pseudopotential [37] plane wave approach was adopted in the calculations, based on the Perdew-Burke-Ernzerhof generalized gradient approximation [38] for estimating the exchange-correlation functional. The energy cut-off for the plane wave basis sets was 400 eV, and a Monkhorst-Pack scheme [39] (which was used for \(k\)-point sampling of the Brillouin zone) of \(2\times 2\times 1\)\(k\)-points was used for structural optimization processes and self-consistent  
Figure 1: Schematic illustration of a silicon-vacancy (SiV) center in bulk diamond, with Si positioned in a split-vacancy site, showing \(D_{3d}\) symmetry alongside the symmetry axis of the [111] direction. The original two C atoms are shown as white spheres and illustrated using the letter V (vacancy). energy calculations of diamond slabs. For geometry optimization, the convergence tolerances for energy, maximum force, and maximum displacement were \(1.0\times 10^{-5}\) eV/atom, 0.03 eV/A, and 0.001 A, respectively. The lattice constant of diamond was calculated as being 3.558 A (where a \(6\times 6\times 6\) \(k\)-points mesh was used for an eight-atom primitive diamond lattice), which was in excellent agreement with the experimental value [40]. All parameters were selected following careful test calculations, ensuring the convergence of subsequent results. The H-, O- and N-terminated diamond (001) surfaces were modeled as slabs with sizes of roughly 10.06 \(\times\) 10.06 \(\times\) 24.27 \(A^{3}\) and comprised 14 C layers with 224 C atoms in total when perfect, alongside a vacuum layer thickness of 11 A. Here, the thick vacuum layer of 11 A was large enough to ensure that the electric field decays were negligible within the vacuum, thus preventing artificial charge transfer between the two polar ends of the slab and its replica. Meanwhile, 14 layers of C, with the dangling bonds of the lowest layer of C atoms saturated with hydrogen atoms, were tested and found to be adequate for simulating bulk diamond. These H atoms were, together with the lowest layer of C atoms, fixed in space during the geometry optimization process and during subsequent calculations. The remaining atoms in the slabs were allowed to relax freely using the Broyden-Fletcher-Goldfarb-Shanno algorithm [41]. The formation energy \(E_{f}[X^{q}]\) of defect or impurity \(X\) in a charge state of \(q\) with respect to the perfect H-, O-, and N-terminated surface or bulk diamond, respectively, is defined as [42]  
\[E_{f}[X^{q}]\ =\ E_{tot}^{defect}[X,\ q]-\ E_{tot}^{perfect}\ -\ \sum n_{i}\mu_{i}\ -\ q[E_{F}\ +\ E_{V}+\Delta V], \tag{1}\]  
where \(E_{tot}^{defect}[X,\ q]\) is the total energy of the slab or bulk diamond containing defects in the supercell, \(E_{tot}^{perfect}\) is the total energy of the perfect H-, O-, or N-terminated diamond (001) surface or bulk diamond, \(n_{i}\) denotes the total number of atoms of element \(i\) (host atoms C or impurity silicon atoms (Si) removed from (\(n_{i}<0\)) or added to (\(n_{i}>0\)) in the diamond model), and \(\mu_{i}\) refers to the chemical potential of element \(i\). The value of \(\mu_{C}\) is the energy of an individual C atom, obtained from the total energy of a perfect 216-atom diamond supercell, whereas \(\mu_{Si}\) was deduced from cubic silicon carbide (SiC). \(E_{F}\) is the Fermi level with respect to the valence band maximum (VBM) (\(E_{V}\)), and \(\Delta V\) is the necessary correction for aligning the electrostatic potentials between the perfect surface and a surface with defects in a charged state. The final term is zero, because, in our work, we primarily investigated the \(E_{f}\) of a diamond surface and of bulk diamond with defects in a neutral charge state. A negative value of formation energy corresponded to an easier formation process with lower energetic barriers.