Ab-initio calculations of shallow dopant qubits in silicon from pseudopotential and all-electron mixed approach

Obtaining an accurate first-principle description of the electronic properties of dopant qubits is critical for engineering and optimizing high-performance quantum computing. However, density functional theory (DFT) has had limited success in providing a full quantitative description of these dopants due to their large wavefunction extent. Here, we build on recent advances in DFT to evaluate phosphorus dopants in silicon on a lattice comprised of 4096 atoms with hybrid functionals on a pseudopotential and all-electron mixed approach. Remarkable agreement is achieved with experimental measurements including: the electron-nuclear hyperfine coupling (115.5 MHz) and its electric field response (−2.65 × 10−3 μm2/V2), the binding energy (46.07 meV), excited valley-orbital energies of 1sT2 (37.22 meV) and 1sE (35.87 meV) states, and super-hyperfine couplings of the proximal shells of the silicon lattice. This quantitative description of spin and orbital properties of phosphorus dopant simultaneously from a single theoretical framework will help as a predictive tool for the design of qubits. Modelling the quantum transport properties of qubit arrays and the electronic properties of dopant qubits is computationally challenging yet crucial for device optimization in quantum computing. Here, the authors compare different DFT-based methods to describe the properties of shallow donor-based qubits in silicon.

T he exponential progress of microelectronics over the last decade has been dramatically dependent on the excellent host material, silicon. The weak spin-orbit coupling and the existence of zero nuclear spin isotopes have made silicon a prime semiconductor platform for the emerging field of quantum information technology 1 . Shallow hydrogenic impurities (e.g., phosphorus, bismuth, and arsenic) 2 offer natural Coulomb confinement to single electrons and provide access to single electron and nuclear spins for encoding quantum information 3 . The relatively large electronic wavefunctions of the donors allow gate control of orbital and spin properties. This has led to the demonstration of high fidelity single 4 and two-qubit logic gates 5,6 with exceptional spin coherence times 7,8 and all-electrical spin readout and initialization 9 . Advanced lithographic techniques have been developed to deterministically place one to few phosphorous donors in a plane of silicon with a lattice constant precision 10,11 , leading to the fabrication of single-atom transistors 12 , few-atom thick nanowires 13 , and low-noise atomic spin qubits 4,5 .
Accurate modeling of the spin and orbital properties of shallow dopants is crucial for designing and optimizing high-performance spin qubits 14 . Two dominant approaches in the community are the multi-valley effective mass technique and the Slater-Koster tight-binding (TB) method, as these enable calculations at the length scale of 10-100 nm. While effective mass theory is a continuum minimal-basis approach lacking atomic resolution, Slater-Koster TB theory is an atomistic full Brillouin Zone method 15 . However, both these methods are semi-empirical and single-electron in nature, relying on different levels of simplifications and parameters from ab-initio methods. The ability to extend the simulations to a larger number of electrons and compute properties without fitting parameters can be achieved by ab-initio density functional theory (DFT). DFT has been extensively applied to elucidate the electronic properties of semiconductors and is currently the only approach to obtain absolute values of Fermi contact hyperfine (HF) parameters between electron and nuclear spins in multi-electron systems 16 . Although DFT has succeeded in providing accurate hyperfine parameters for some highly localized deep impurity systems 17 , it is extremely challenging to simulate relatively delocalized shallow dopant states due to the system size limitations and excessive computational burden 18 .
While DFT calculations on shallow dopants have improved over the years, past works were able to calculate either the hyperfine coupling or binding energy with accuracy, but not both. Swift et al. 19 employed pseudopotential and extrapolation approaches with HSE and PBE functionals working in tandem, achieving outstanding accuracy in binding energy calculations of bismuth and arsenic donors in Si. However, the predefined frozen core region of Kohn-Sham orbitals results in an incomplete description of the hyperfine coupling. Smith et al. 20 achieved success in the binding energy calculations directly from the wavefunctions extracted from the supercells containing 10648 atoms and the impurity potential of the phosphorus donor using an empirical model. And Gerstmann 21 accurately reproduced the isotropic hyperfine and super-hyperfine parameters of shallow donors in Si using Green's function approach. However, the empirical correction included in these methods lacks the capacity to provide an accurate evaluation of the hyperfine coupling and orbital splitting of the donor simultaneously from a single theoretical framework. In addition, electric field dependency of the hyperfine coupling, critical for quantum computing, has never been obtained from ab-initio calculations.
In the present work, these obstacles are overcome by using pseudopotential (PSP) and all-electron (AE) mixed Gaussian type localized basis sets combined with hybrid functional to conduct large-scale DFT calculations containing up to thousands of atoms (4096 atoms; 4 × 4 × 4 nm 3 ). A host of important spin and orbital properties of a single phosphorus donor in silicon (Si:P), such as hyperfine coupling and its electric field dependency, superhyperfine couplings at silicon lattice points within a few shells of the impurity, and excited orbital (valley-orbit) splitting are accurately evaluated simultaneously. To obtain consistently accurate values of all these quantities together, the wavefunction, and hence the electron density needs to be accurately described not only close to the donor nucleus but also over a large part of the silicon lattice extending over a Bohr radius of 1-2 nm. Figure 1 shows a schematic of a single P dopant as a spin qubit in a silicon host, with the corresponding charge density difference map obtained from DFT showing the electronic interactions close to the donor atom. We focus on P donors due to their technological significance in quantum computing 2 , but the presented methods also apply to other shallow Group V dopants in silicon. For comparison, three separate approaches are used, namely, (1) the Pseudopotential approach (PSP) in VASP 22 , (2) the All-Electron approach (AE) in CP2K 23 , and (3) the mixed Pseudopotential and All-Electron approach (PSP-AE) in CP2K. The three-pronged approaches enable us to contrast the shortcomings of each technique and help us to assess computational cost versus accuracy. We find that the use of proper basis functions, state-ofthe-art functionals, and extrapolation method can yield quantitatively accurate results without significantly increasing the computational burden.
Recently papers have emphazied the importance of using allelectron level basis sets in spin-Hamiltonian calculations for quantitative predictions of spin dynamics 24,25 . As an AE approach based on Gaussian type localized basis sets, CP2K solves the inherent problems induced by the PSP approximation such as the disregard of the core spin-polarization effect and the contribution of the exchange-correlation potential in the vicinity of the nucleus 26 . Hence, it provides high computational efficiency and superior accuracy for core-electron related properties such as hyperfine coupling, nuclear magnetic resonance, and g-factor [27][28][29][30] . Additionally, hybrid functional is employed, which can further improve the accuracy of the wavefunction amplitude and the electron localization at the donor nucleus 19,31 . With the combination of AE approach Fig. 1 Illustration of single P donor in Si host in a quantum information processor. A-gates above the donors control the hyperfine coupling and hence the resonance frequency of the nuclear spin qubits; J-gates between the donors control the electron-mediated coupling between adjacent nuclear spins. The black background Si region below the barrier is used to reflect a more bulk-like scenario for qubits embedded into Si at depths >10 nm. The embedded charge density difference map illustrates the electronic interaction around the donors (isosurface value of 1 × 10 −3 e/ Å 3 ), where cyan area represents electron depletion and yellow area represents electron accumulation. and hybrid functional, a value of 115.5 MHz HF constant is predicted in the present work which is in excellent agreement with the experimental value of 117 MHz 32 . The corresponding binding energy 1 s(A 1 ) of the ground state and energies of excited 1 s(T 2 ) and 1 s(E) states are found to be 46.07 meV, 37.22 meV, and 35.87 meV, respectively, in good agreement with the experimental values 33,34 . Furthermore, the Stark shift of the hyperfine coupling is of utmost importance for electrically tuning electron-nuclear resonance frequency of donor spin qubits and it is computed from DFT in quantitative agreement with measurements. Finally, we show that the computational burden of an AE method can be largely alleviated by using a PSP-AE mixed method which can achieve the same accuracy as the AE approach.
The method introduced in the present work, is currently, to the best of our knowledge, the only approach that can simultaneously reproduce all the important properties of shallow donors in predictive accuracy (e.g., anisotropic/isotropic hyperfine, superhyperfine interactions, and their electric field dependency, as well as binding energy calculations). Therefore, the present work represents significant progress in ab-initio calculations of shallow spin defects and paves the way for predictive exploration of similar quantum defects in a host of emerging materials.

Results
Hyperfine coupling. The Hamiltonian describing the hyperfine interaction between an electronic spin S and a nuclear spin I is known as: where A I iso is the isotropic or the Fermi contact hyperfine coupling at the site of the nuclear spin and is calculated by 35 : where μ 0 is the magnetic susceptibility of free space. γ e and γ I are the electron gyromagnetic ratio and the nuclear gyromagnetic ratio of the nucleus at R I , respectively. δ T ðrÞ is a smeared out δ function to take scalar-relativistic effects into account 35 and ρ s is the spin density. For all-electron basis sets approach, the relativistic effects onto the A I iso are ignored and the last term becomes ρ s R I À Á . If the nuclear spin site is the impurity site itself, we obtain the hyperfine coupling of the donor; while if the site is a silicon site in the vicinity of the donor, occupied by a 29 Si isotope of 28 Si, we obtain the super-hyperfine couplings.
Pseudopotential approximation approach (VASP-PSP). In order to predict the HF constants from DFT results, a very recently proposed extrapolating method is employed in the present work with minor modification 19 . Briefly, since the error induced by the exponential tail of wavefunctions inversely scales with the supercell size, the calculated HF constants (A I iso ) of a variety of supercell sizes obtained from the PBE functional are plotted as a function of 1/N ( Fig. 2a and Supplementary Table S1), where N represents the number of atoms in the supercell. A fitting function is then obtained based on the polynomial least-square method starting from the supercell size of n = 4 (N = 512) to n = 6 (N = 1728), which gives an expression of A I iso = 7.47x 2 + 72.8x + 43.6. Subsequently, an intercept value of 43.6 MHz is achieved by the extrapolation of N⟶∞ (i.e., x → 0), corresponding to the predicted HF constant.
In contrast to the experimental HF constant of 117 MHz 32 , the predicted value is significantly underestimated owing to the excess delocalization of wavefunction inherent in PBE. The use of hybrid functional can improve the description of localization. As reported, the fitting functions obtained from PBE and HSE functionals provide exactly the same slope value 19 . Hence the fitting function A I iso from PBE is subsequently employed to fit HF results obtained from HSE at the supercell size of n = 5 (i.e., A I iso = 136.9 MHz and x = 1) (Fig. 2a). This gives an intercept value of 56.7 MHz (Supplementary Table S2 a slight improvement over PBE. However, the HF coupling is still considerably underestimated due to the inherent problems of the PSP approximation method discussed in the introduction section. It is worth noting that due to the excessive computational burden and lower computational efficiency of large-scale systems, a limited supercell size in the VASP-PSP and the previously reported linear-extrapolation 19 approaches results in a weak convergence of the fitting function (i.e., the linear coefficient is dominant in A I iso ), leading to a significant underestimation of the HF constants. This poor convergence issue and underestimation of the HF constants can be considerably eliminated by the employment of all-electron basis sets and larger supercell sizes as described in the following CP2K-AE and CP2K-PSP-AE approaches, providing a significant improvement in evaluating the core-electron dominant properties of shallow donors.
All-electron approach (CP2K-AE). To eliminate the inherent problems induced by the PSP approach, GAPW method implemented in CP2K using an AE basis set is employed to calculate the Fermi contact HF parameter. Due to the higher computational efficiency, a larger supercell size is evaluated for the PBE approach, containing 4096 (n = 8) atoms. The calculated HF parameters for the donor are presented in Supplementary Table S1. In contrast to VASP-PSP approach, the use of AE scheme provides higher HF constants for different supercells due to the more localized nature of Gaussian type basis sets and improved description of the core spinpolarization effect. The same extrapolating method as the VASP-PSP approach is used to predict the HF constants. With the increase in supercell size, the quadratic coefficient of the polynomial fitting function becomes more dominant. The fitting function starting from the supercell size of n = 5 (N = 1000) to n = 8 (N = 4096) for PBE results gives an expression of A I iso = 20.22x 2 + 15.3x + 110.0 (Fig. 2b). This is subsequently employed to extrapolate HF results obtained from HSE at n = 5 (i.e., A I iso = 151.0 MHz and x = 1), giving an intercept value of 115.5 MHz (Supplementary Table S2), corresponding to the final prediction of HF constant, in excellent agreement with the experimental value of 117 MHz 32 .
Pseudopotential approximation and all-electron mixed approach (CP2K-PSP-AE). Although the AE approach significantly improves the inherent delocalization problem induced by the PSP method, the consideration of all electrons for each atom considerably increases the computational cost, making large-scale system calculations expensive. Therefore, a PSP and AE mixed approach is employed in the present work to efficiently and accurately evaluate the hyperfine coupling interaction for largescale systems. To achieve consistent results in contrast to pob-DZVP AE approach, pseudopotentials from the same family (DZVP-GTH-PBE) were used. The lattice parameter, the bond length of Si-Si, and the angle of Si-Si-Si of the geometry optimized from DZVP-GTH-PBE pseudopotential using PBE are 5.479 Å, 2.373 Å, and 109.5°, consistent with the results from pob-DZVP AE approach using PBE (5.486 Å, 2.375 Å, and 109.5°, respectively).
The CP2K-PSP-AE scheme provides similar HF constants compared with CP2K-AE approach (Supplementary Table S1). The same extrapolating method starting from the supercell size of n = 5 (N = 1000) to n = 8 (N = 4096) is used, providing a fitting function of A I iso = 19.49x 2 + 18.9x + 102.3 (Fig. 2c). A subsequent extrapolation on HSE results n = 5 gives an intercept value of 112.6 MHz (Supplementary Table S2), in good agreement with the experimental value and the HF constant predicted from CP2K-AE approach. Furthermore, a comparison of the computational cost of CP2K-PSP-AE mixed approach with CP2K-AE approach summarized in Supplementary Table S3 indicates a significant improvement in the computational speed (4.5-5 times faster), permitting accurate spin-density description with higher computational efficiency.
Super-hyperfine coupling. In addition to the calculation of HF parameters, super-hyperfine (SHF) interactions are evaluated. In contrast to the HF parameters which describe the wavefunction distribution on the central donor nucleus, SHF interactions depend on the wavefunction distribution in the silicon lattice sites surrounding the donor which can be occupied by the spin ½ isotope 29 Si. Five shells of silicon atoms closest to the donor center are investigated, namely (1,1,1), (2,2,0), (1,1, 3), (0,0,4), and (3,3,1) in units of a 0 /4 with a 0 being the lattice constant.
The SHF parameters are calculated in a similar way to the HF parameter through the extrapolation approach. The fitting functions are obtained starting from the supercell size of n = 4 (N = 512) to n = 6 (N = 1728) for VASP approach (Supplementary Figure S1) or n = 8 (N = 4096) for CP2K approach (Supplementary Figs. S2 and S3). These fitting functions are subsequently used to fit the SHF results obtained from HSE at the supercell size of n = 5, giving the final prediction of SHF constants. As shown in Fig. 3, the 4 th -nearest-neighbor shell (0,0,4) gives the highest SHF constant owing to the Kohn-Luttinger oscillations 36 , consistent with the experimental results 32 and data from Green's function calculation (LMTO-GF) 21 . A better agreement with Green's functional method is obtained from AE approach.
It's worth noting although the differences between VASP and CP2K results might be partially induced by the different optimized geometries, the consideration of the core electrons (i.e., PSP vs. AE basis sets) is the dominant factor that influences the HF and SHF calculations. These results prove that a combination of hybrid functional and AE approach can significantly improve the delocalization problem induced by PSP and PBE functional and accurately describe the spin-density of the donor state. The anisotropic SHF constants of the five shells of silicon atoms closest to the donor center are subsequently analyzed. Since most of the anisotropic parts (kHz) are at least 3-4 orders of magnitude smaller compared with the isotropic (MHz) SHF constants, the extrapolation method is not applicable in this case due to the sensitivity and significant fluctuation of the anisotropic parts as a function of supercell size. Hence, the anisotropic constants calculated from the largest supercells (n = 6 for VASP and n = 8 for CP2K) were used as the final prediction. As shown in Supplementary Table S4, both VASP and CP2K results exhibit good agreement with the experimental value 37 . The only exception is the (0,0,4) results obtained from VASP which provide different values of |B zz |(23 kHz) and |B xy |(49 kHz) components, contradictory to the experimental measurements (41 kHz for both |B zz | and |B xy |). In this sense, CP2K results show better consistency with the experimental values, where an anisotropy of these values could not be resolved.
Electric-field dependent hyperfine coupling. Different from previous studies that investigate the influence of strain-induced internal field on the hyperfine coupling 19,36,38 , an external static electric field is employed in the present work. The presence of an external electric field (ε) can distort the shape of the donor wavefunction and pull the wavefunction away from the donor site. This will result in a decrease in the electric field dependent hyperfine parameter A(ε) which is proportional to jΨðε;R I Þj 2 , whereR I represents the donor site. Experimental measurements 39 have shown a quadratic dependence of the hyperfine coupling with electric field, which is parameterized as 40 : where 4Aε ð Þ represents the change of hyperfine parameter in the presence of electric fieldε (i.e., Aε ð Þ À A 0 ). η 2 and η 1 are the coefficients of the quadratic and linear terms respectively. For bulk donors, η 2 was measured to be −3.7 × 10 −3 μm 2 /V 2 for Si:Sb and η 1 was found to be negligibly small. Since 4Aε ð Þ is a relative value, there is no need to employ all-electron basis sets for all atoms; hence only CP2K-PSP-AE method is used. Detailed calculation parameters are described in the Methods section.  39 . Although the calculated quadratic coefficient is larger than the experimental measurement owing to the limitation of supercell size, the trend shows that it is approaching the experimental value with the increase of supercell size. In contrast, the linear coefficient η 1 is relatively stable and close to zero for all supercells.
The effect of electric field on the isotropic and anisotropic parts of HF and SHF interactions is further analyzed. In this part, AE basis sets were used for P donor and Si atoms located in (1,1,1) and (0,0,4) shells, while the rest Si atoms were treated with PSP in order to reduce the computational burden. (1,1,1) and (0,0,4) shells were chosen to reflect the effect of electric field on the inner and outer 29 Si atoms. Compared with the isotropic parts, the value of the anisotropic parts are at least two orders of magnitudes smaller except the (1,1,1) shell (Supplementary  Table S6), consistent with the experimental results 37 . And the change of the anisotropic parts under the effect of electric field is at least one order of magnitude smaller compared with the isotropic parts. The anisotropic parts of the HF constant of P donor are increased under electric field, due to the symmetrybreaking effect, except the |B xy | component. In terms of the anisotropic parts of the four 29 Si atoms located in inner (1,1,1) shell, only the diagonal components (i.e., |B xx |, |B yy |, and |B zz |) are increased under electric field, while the non-diagonal components show an opposite trend. For the anisotropic parts of the six 29 Si atoms located in outer (0,0,4) shell, the effect of electric field can be summarized as two cases: (1) For Si1 and Si2 atoms which are along the electric field z-direction and located on the opposite sides of the P donor (i.e., (0,0,4) and (0,0,−4)), the diagonal components (i.e., |B xx |, |B yy |, and |B zz |) and one non-diagonal component (|B xy |) are decreased under electric field, while the rest two non-diagonal components (i.e., |B xz |, |B yz |) are increased; (2) In contrast, the anisotropic parts of the rest four Si atoms showing a different trend, with one diagonal component decreased while all the rest diagonal and non-diagonal components increased. It's worth noting that the |B zz | (19 kHz) and |B xy | (27 kHz) components of (0,0,4) shell in this case are different, contradictory to the experimental measurements, showing the same issue as previously discussed in the anisotropic VASP SHF results (Supplementary Table S4). This further proves the importance of using all-electron basis sets for all atoms to correctly describe the core-electron contributions to the spinpolarization effect. Binding energy calculation. Binding energies (E) are calculated from the Kohn-Sham eigenvalues as described in a previous work 19 : where ϵ cb represents the eigenvalue of the conduction band minimum (CBM) in bulk Si and ϵ donor represents the resultant eigenvalue in P doped Si. ∇V, which is calculated by the Freysoldt method 41,42 , is the correction term and is used to align the Kohn-Sham levels between the bulk Si and Si:P system (Supplementary Table S7).
Pseudopotential approximation approach (VASP-PSP). The final binding energy is predicted in a similar way to the calculation of HF constant by fitting the binding energy results from PBE as a function of the supercell size starting from n = 4 (N = 512) to n = 6 (N = 1728) using the linear least-square method. This provides a fitting function of E b = 0.0312x + 0.0325 (Fig. 5a). A subsequent extrapolation to N⟶∞ (i.e., x⟶0) gives an intercept value of 32.5 meV, representing the prediction of binding energy from the PBE functional. However, the delocalization of the wavefunction and the underestimation of the exchange-splitting effect (δ ex ) from PBE functional result in an underestimation of the binding energy compared with the experimental value of 45.59 meV 33 . Therefore, HSE functional is employed and the underestimation of δ ex in PBE is considered (Supplementary Table S7). δ ex represents the difference between the spin-up and spin-down eigenvalues of the donor state 19 and it also shows a clear trend inversely scaling with the supercell size (N) (Fig. 5b), giving slope values of 0.0086 (s PBE;δ ex ) and 0.0304 (s HSE;δ ex ) for PBE and HSE functionals, respectively. The final HSE binding energy slope (s HSE ) is calculated by adding half the HSE exchange-splitting slope ( 1 2 s HSE;δ ex ) on the spin-averaged PBE slope (s PBE À 1 2 s PBE;δ ex ) 19 : This gives a slope value of 0.0421 (s HSE ), which is subsequently used to extrapolate the binding energy calculated from HSE functional at the supercell size of n = 5 to the N⟶∞ limit (i.e., E b = 0.088 and x = 1) ( Table 1). The resultant intercept value corresponds to the final prediction of the binding energy, with the value of 46.07 meV, in excellent agreement with the experimental value (45.59 meV).
The same method was employed to calculate the higher 1 s(T 2 ) and 1 s(E) orbital states, giving the final prediction of 37.22 meV and 35.87 meV, respectively (Table 1 and Supplementary Figure S4), showing good agreement with the experimental value of 33.88 meV and 32.54 meV 34 . The splitting of these states is due to the so-called valley-orbit splitting originating from strong donor core potential. These low-manifold states in the donor spectra are responsible for several properties of donor spin qubits including the spin-lattice relaxation times 43 . These results indicate that HSE has the capacity of reducing the error induced by the PBE underestimation of binding energy, while the size limitation of HSE functional could be corrected by the PBE scaling, achieving significant agreement with experiments.
All-electron approach (CP2K-AE). Unfortunately, binding energy cannot be calculated currently in this method because the output of electrostatic potential which is required for calculating the correction term (e∇V) is incompatible with GAPW method in CP2K.
Pseudopotential approximation approach (CP2K-PSP). Unlike GAPW approach, electrostatic potential can be obtained with GPW method in CP2K. Therefore, the binding energy is calculated in the same method described in the VASP-PSP approach. While performing the hybrid functional calculation using the pseudopotential approximation approach in CP2K (CP2K-PSP), the auxiliary density matrix method (ADMM) was employed to reduce the computational burden and accelerate the calculation speed 44 . The calculated results are summarized in Table 1. The same extrapolating method is used by linearly fitting the PBE results as a function of the supercell size starting from n = 4 (N = 512) to n = 8 (N = 4096), This provides a fitting function of E b = 0.0345x + 0.0319 (Fig. 5c) in the prediction of Fermi contact HF, its electric field dependency, and SHF has been achieved with all-electron level analysis and hybrid functional working in tandem. This highlights the importance of using all-electron level analysis to provide satisfactory elucidation of shallow donor systems. Results also reveal that a combination of pseudopotential approximation and allelectron level analysis can permit an accurate description of spindensity at the donor nucleus and generate all calculated parameters in a considerably more efficient way. Additionally, calculations of energies of the ground and excited states have exhibited outstanding agreement with the experimental value, demonstrating the predictive nature of the extrapolating approach for shallow donor systems. In the qubit architectures, the dopants are typically embedded into silicon at depths >10 nm (Fig. 1), which is larger than several Bohr radii of the dopants. Hence, the bulk-like scenario is more appropriate for these qubits. While the method introduced in the present work is focused on evaluating the properties of bulk systems, further consideration of surface effects for dopants buried a few nm below the barrier may be also achievable, owing to the superior computational efficiency of the PSP-AE mixed approach and its capacity of investigating large-scale systems containing thousands of atoms.

Methods
Pseudopotential approximation approach (VASP-PSP). The first set of calculations were conducted with the projector augmented-wave (PAW) method implemented in the Vienna Ab-initio Simulation Package (VASP) 22  most accurately reproduce the band gap (1.14 eV) and lattice constant (5.473 Å). A matrix diagonalization method was applied for wavefunction optimization. The convergence threshold on energy for self-consistent field was set to 10 −5 Hartree (Ha). Double-ζ valence with polarization quality (pob-DZVP) for Si 52 and 6-311G** Gaussian type basis set for P 53 , respectively, were used, together with a 450 Ry cutoff energy for the auxiliary plane-waves grid. The largest supercell corresponded to n = 8 containing 4096 atoms. A single k point was employed for all supercells. The geometry was optimized using Broyden-Fletcher-Goldfarb-Shanno (BFGS) [54][55][56][57] and Limited-memory BFGS 58 optimizers, with the convergence criteria of 3 × 10 −4 Ha/Bohr and 1.5 × 10 −3 Bohr for the root mean square (RMS) force and RMS atomic displacement, respectively.
Pseudopotential approximation and all-electron mixed approach (CP2K-PSP-AE). The third set of calculations included two steps for PBE functional calculations: (1) The geometry optimization was performed with gaussian and planewaves (GPW) scheme in CP2K/Quickstep 59 in combination with Goedecker −Teter−Hutter pseudopotentials 60,61 and (2) Double-ζ valence with polarization quality (pob-DZVP) for Si 52 and 6-311G** Gaussian type basis set for P 53 , respectively, were used to perform a single-point calculation based on the optimized geometry. For hybrid functional calculations, the procedure was the same as that used in the second set of calculations (i.e., using AE scheme for both geometry optimization and single-point calculation), owing to the vastly lower computational cost of HSE in AE scheme in contrast to that in PSP approach. All the rest of the parameters and convergence criteria were the same as that used in the second set of calculations.
Electric field-dependent hyperfine parameter calculation. A static electric field was applied to evaluate its influence on hyperfine coupling. A variety of electric field values was investigated, including: 0, 2E-07, 4E-07, 6E-07, 8E-07, 10E-07, and 12E-07 (units are in a.u.). In order to consider the symmetry-breaking effect induced by the electric field, geometry relaxation of ligand atoms was performed after applying an electric field with the GAPW scheme 51 . Goedecker−Teter −Hutter pseudopotential for Si 60,61 and modified pcH-2 Gaussian type basis set for P were used 62 , respectively. The rest of the parameters and convergence criteria are the same as those used in the second set of calculations.

Data availability
The data that support the findings of this study are available in the NOMAD repository at https://doi.org/10.17172/NOMAD/2022.06.25-1. All other data are available from the corresponding authors on reasonable request.