9. SCF and impurity state analysis of P-in-Si donor system

9.1. Introduction

In this tutorial, we demonstrate how to perform a self-consistent field (SCF) calculation and impurity state analysis for a phosphorus donor in bulk silicon using density functional theory (DFT), following the methodology proposed in Songqi Jia et al., Appl. Phys. Lett. 125, 184001 (2024). The goal is to accurately extract the impurity states introduced by the dopant atom, which requires particular attention to the size and resolution of the simulation domain.

Atomistic simulations of this kind require extremely large supercells. In our case, a system of over 10,000 atoms is used to faithfully capture the spatial extent and symmetry of the impurity state’s wavefunction. This is necessary because silicon’s conduction electrons have a large effective Bohr radius—about 1.8 nanometers (J. S. Smith et al., Sci. Rep. 7, 6010 (2017)). If smaller supercells are used, the impurity wavefunctions become artificially confined by periodic boundaries, leading to unphysical localization and inaccurate energy levels. By performing SCF calculations on sufficiently large systems, we ensure that the impurity states are converged with respect to supercell size and can be used as physically meaningful quantities.

These converged SCF impurity states serve as a well-defined single-particle basis for many-body approaches, such as exact diagonalization, to compute critical quantum transport properties like the charging (addition) energy of a donor quantum dot. While this tutorial focuses exclusively on the SCF portion, the results presented here provide the essential foundation for more advanced quantum modeling and simulation.

We will:

1. Perform an SCF calculation for a large silicon supercell with a single phosphorus dopant.

2. Identify the impurity states from the Kohn-Sham spectrum.

3. Verify that the impurity states have their largest LCAO coefficient on the dopant atom.

4. Extract the six near-degenerate impurity states (A₁, T₂, E) and analyze their energies.

5. Visualize the wavefunction of the A₁ impurity state.

6. Compute the charge density of the two-electron state at the SCF level.

Note: This is a non-quantum demo, intended to showcase SCF and wavefunction-level information only.

9.2. System description

The system under investigation is a cubic silicon supercell consisting of 10,648 atoms, in which a single phosphorus (P) atom replaces a silicon (Si) atom at the center. This creates a model of a donor quantum dot embedded in a bulk Si environment. The supercell has a side length of approximately 59.7 Å, providing over 3 nm of buffer between the impurity and the periodic boundaries— a necessary condition for minimizing finite-size effects on the localized impurity states.

Fig. 9.2.2 Atomic structure of the P-in-Si donor system. Left: Full 10,648-atom supercell with a cubic lattice of side length 59.7 Å. Right: Zoomed-in view showing the phosphorus impurity atom (yellow) at the center, surrounded by its nearest-neighbor silicon atoms.

The impurity state introduced by the P atom is expected to be hydrogenic in character, with the lowest-energy 1s-like state split into an \(A_1\) singlet, \(T_2\) triplet, and \(E\) doublet due to the valley-orbit interaction in silicon. These states are localized around the impurity and extend over a few nanometers, requiring a supercell much larger than this spatial extent to avoid unphysical wavefunction overlap with its periodic images.

To accurately model this system, we use the following computational setup:

• Exchange-correlation functional: LDA (local density approximation)

• Pseudopotentials: Troullier-Martins norm-conserving

• Basis set: Linear Combination of Atomic Orbitals (LCAO), double-zeta polarized (DZP)

• K-point sampling: \(\Gamma\)-point only

• Real-space mesh: \(206^3\) points (for solving the Hartree potential)

This setup enables us to achieve convergence of the impurity states and faithfully capture their spatial and energetic features.

9.3. SCF calculation

First, we calculate the ground state energy of the phosphorus-doped silicon system using a self-consistent field (SCF) method within DFT. The LCAO method is enabled to make the large system size tractable while still providing a high-quality description of the electronic structure.

The input file for the RESCU SCF calculation is structured as follows:

% Enable LCAO method
LCAO.status = true;

% Define calculation type and output settings
info.calculationType = 'self-consistent';
info.savepath = 'p_in_si';

% Define atomic structure
atom.xyz             = 'atoms_10000.xyz';
element(1).species   = 'Si';
element(1).path      = 'Si_TM_LDA.mat';
element(2).species   = 'P';
element(2).path      = 'P_TM_LDA.mat';

% Define simulation domain and resolution
domain.latvec        = [[59.7377 0 0]; [0 59.7377 0]; [0 0 59.7377]];
domain.lowres        = 0.55;

% Define exchange-correlation functionals
functional.list      = {'XC_LDA_X', 'XC_LDA_C_PW'};
functional.libxc     = true;

% Define k-point sampling and mixing settings
kpoint.gridn         = [1, 1, 1];
kpoint.sigma         = 0.01;
mixing.type          = 'density';
mixing.method        = 'pulay';
mixing.tol           = [1e-4, 1e-4];
mixing.maxhistory    = 6;

% Set simulation memory and eigensolver settings
eigensolver.maxit    = 15;

% Define units
units.length         = 'Angstrom';
units.energy         = 'eV';

% Use spin-degenerate calculation
spin.type            = 'degenerate';

% Define differential operator
diffop.method        = 'fft';

% Set SCF iteration cap
option.maxSCFiteration = 200;

% Disable GPU usage
gpu.status           = 0;

This configuration enables SCF convergence in large supercells (over 10,000 atoms) with a balance of numerical stability and speed.

After running the SCF, we obtain converged Kohn-Sham (KS) energies and orbitals, which are then used to identify and analyze the impurity states.

9.4. Verification of impurity states

After the SCF calculation is completed, we analyze the KS eigenstates to identify the localized impurity states introduced by the phosphorus atom. In silicon, a hydrogenic donor state is expected to form near the conduction band minimum (CBM), exhibiting a characteristic splitting pattern due to valley-orbit interaction.

To extract these states, we run a band structure calculation that uses the converged charge density from the SCF step and outputs the eigenstates near the CBM, as well as their corresponding wavefunctions and LCAO expansion coefficients.

The RESCU input file is as follows:

% Define calculation type and charge density input
info.calculationType = 'band-structure';
info.savepath        = 'bs_p_in_si';
rho.in               = 'p_in_si';

% Enable LCAO method
LCAO.status          = true;

% Define atomic species
element(1).species   = 'Si';
element(1).path      = 'Si_TM_LDA.mat';
element(2).species   = 'P';
element(2).path      = 'P_TM_LDA.mat';

% Define k-point sampling (only Gamma point)
kpoint.kdirect       = [0, 0, 0];

% Set range of bands to output
eigensolver.bandi    = [21296, 21315];

% Save wavefunctions and LCAO coefficients
option.saveWavefunction = true;
option.saveLCAO         = true;

This configuration extracts the KS eigenstates near the conduction band minimum and saves their wavefunctions and LCAO expansion data for post-processing.

From the KS energy spectrum, we locate six eigenstates just above the CBM:

• One singlet state labeled \(A_1\)

• Three degenerate triplet states labeled \(T_2\)

• Two degenerate doublet states labeled \(E\)

These six states are tabulated below, with their energies referenced to the Fermi level:

Table 9.4.1 Identified impurity state energies

Impurity state	Energy (eV)
\(A_1\)	0.015
\(T_2\) (1)	0.025
\(T_2\) (2)	0.025
\(T_2\) (3)	0.025
\(E\) (1)	0.026
\(E\) (2)	0.026

This energy structure is consistent with the hydrogenic model for a donor in silicon. To confirm these are truly impurity states, we analyze the spatial distribution and LCAO composition of the wavefunctions.

Each KS state \(\Psi_n(r)\) is expressed as a linear combination of atomic orbitals:

\[\Psi_n(r) = \sum_{i} C_{ni} \phi_i(r),\]

where \(\phi_i(r)\) is the atomic orbital on site \(i\), and \(C_{ni}\) is the LCAO expansion coefficient. For a true impurity state, the largest value of \(|C_{ni}|^2\) should be located on the phosphorus atom.

By inspecting the coefficients for the six impurity-related KS states, we observe that they all have their maximum orbital weight on the P atom, confirming their localized character. These states are also spatially confined within approximately 2-3 nm around the dopant site, as expected based on the effective Bohr radius in silicon.

This verification ensures that the impurity states are properly isolated from the bulk states and are suitable for further physical interpretation or use in many-body basis construction.

9.5. Wavefunction visualization

To further validate the impurity character of the donor states, we visualize the spatial distribution of the lowest-energy KS wavefunction, corresponding to the \(A_1\) state. This state is expected to be spherically symmetric and strongly localized around the phosphorus atom, with an extent of approximately 2-3 nm.

The wavefunction \(\Psi_{A_1}(r)\) is obtained from the RESCU band structure calculation, where the following options were enabled:

option.saveWavefunction = true;
option.saveLCAO         = true;

This ensures that the real-space representation of the wavefunction is saved and can be post-processed using visualization tools.

We plot the isosurface of \(|\Psi_{A_1}(r)|^2\) in two orthogonal views. The left panel shows the projection along the (001) direction, and the right panel shows the (011) view, obtained by rotating the supercell 45° around the x-axis. An arrow indicates the rotation between the two perspectives.

Fig. 9.5.1 Isosurface plot of the \(A_1\) state wavefunction probability density \(|\Psi_{A_1}(r)|^2\). The isosurface threshold is set to \(7.3 \times 10^{-13}\) atomic units. Left: (001) view. Right: (011) view after a 45° rotation about the x-axis. The wavefunction is centered at the phosphorus atom and exhibits exponential decay in all directions.

This visualization confirms that the \(A_1\) impurity state is spatially confined and localized at the donor site, with a shape and decay profile consistent with a shallow hydrogenic state in silicon. The results validate that the KS states are converged and physically meaningful for further analysis.

9.6. Charge density calculation

In addition to the single-particle wavefunction, we analyze the real-space charge density associated with the impurity states. While the full many-body charge density includes electron-electron interactions, here we focus on both the non-interacting KS picture and the interacting two-electron ground state.

The total non-interacting charge density \(\rho(r)\) is computed directly from the occupied KS orbitals:

\[\rho(r) = \sum_{i}^{\text{occ}} |\Psi_i(r)|^2,\]

where \(\Psi_i(r)\) denotes the KS eigenstate of orbital \(i\). In our case, we focus on the lowest-energy impurity state (\(A_1\)) filled with two opposite-spin electrons.

To account for Coulomb interaction between the two electrons, we construct a many-body Hamiltonian using the KS single-particle states as a basis and solve it by exact diagonalization. The many-body Hamiltonian is:

(9.6.1)\[\hat{H} = \sum_{i r} \epsilon_i \hat{c}^{\dagger}_{i r} \hat{c}_{i r} + \frac{1}{2} \sum_{i j k l, r r'} V_{i j k l} \hat{c}^{\dagger}_{i r} \hat{c}^{\dagger}_{j r'} \hat{c}_{k r'} \hat{c}_{l r}\]

Here, \(\epsilon_i\) are the KS eigenvalues, \(\hat{c}_{i r}\) and \(\hat{c}^\dagger_{i r}\) are annihilation and creation operators for a single-particle state \(i\) with spin \(r\), and \(V_{ijkl}\) are the two-body Coulomb integrals:

(9.6.2)\[V_{i j k l} = \int d^3 r_1 \, d^3 r_2 \, \psi_i^*(r_1) \psi_j^*(r_2) \frac{e^2}{4 \pi \varepsilon |r_1 - r_2|} \psi_k(r_2) \psi_l(r_1)\]

The dielectric constant \(\varepsilon\) is taken as the bulk silicon value. This approach captures the electrostatic repulsion between electrons confined to the localized impurity states.

The exact diagonalization of the Hamiltonian in (9.6.1) is carried out by the QTCAD spin-qubit simulation package [QTCAD], which uses the finite element method (FEM) for spatial discretization. KS wavefunctions are interpolated from the original RESCU real-space grid onto a FEM mesh before computing \(V_{ijkl}\).

The many-body charge density of the two-electron ground state, referred to in the original paper as the two-electron charge density, is then evaluated as:

(9.6.3)\[\langle \hat{\rho}(r) \rangle_{\Psi} = \sum_{j j'} \psi_j^*(r) \psi_{j'}(r) \langle \Psi | \hat{c}^\dagger_j \hat{c}_{j'} | \Psi \rangle\]

Here, \(\psi_j(r)\) are the KS impurity basis states, and \(|\Psi\rangle\) is the ground state of the interacting two-electron system.

The three views below visualize this many-body charge density using a combination of cross-sectional slice, isosurface, and combined renderings. These images correspond to what is referred to as the two-electron charge density in the original publication.

Fig. 9.6.1 Cross-sectional x–y slice of the two-electron charge density \(\langle \hat{\rho}(r) \rangle_{\Psi}\) through the plane containing the phosphorus impurity atom. The charge is sharply localized and decays radially with smooth symmetry.

Fig. 9.6.2 Isosurface of the two-electron charge density at a fixed threshold (\(1.365 \times 10^{27}\) m⁻³), showing tetrahedral extension toward the four nearest-neighbor Si atoms.

Fig. 9.6.3 Combined slice and isosurface rendering of the two-electron charge density. This composite view shows both the internal distribution and outer decay profile of the donor-bound electrons.

These results confirm that the donor-bound electron pair remains tightly localized around the impurity atom and follows the expected hydrogenic profile. This charge density is critical for understanding the electrostatics of donor quantum dots in quantum information applications.

9.7. Conclusion

In this tutorial, we demonstrated how to perform a self-consistent field (SCF) calculation for a phosphorus donor in bulk silicon using a large supercell containing more than 10,000 atoms. Using the RESCU DFT solver with the linear combination of atomic orbitals (LCAO) method, we successfully converged the electronic structure of the donor system at the \(\Gamma\)-point.

We then extracted and verified the six donor impurity states near the conduction band minimum, consisting of a singlet \(A_1\), triplet \(T_2\), and doublet \(E\) levels. These states were identified through energy, spatial localization, and LCAO coefficient analysis — with the largest wavefunction weights centered on the phosphorus atom.

The wavefunction of the \(A_1\) state was visualized in real space and shown to be highly localized with hydrogenic symmetry, consistent with expectations for a shallow donor in silicon. The SCF-derived charge density of the doubly occupied impurity state further confirmed this localization, showing a peak near the dopant and tetrahedral decay toward the nearest neighbors.

While this tutorial does not include many-body interactions or quantum diagonalization, the results presented here form the essential single-particle foundation for such calculations. The converged KS states and real-space orbitals can be used as an input basis for more advanced simulations, such as computing charging energies and spin-dependent transport properties of donor quantum dots.

This workflow highlights the feasibility and importance of large-scale DFT calculations for accurately capturing impurity states in realistic atomistic models, and provides a stepping stone for future quantum device modeling.

References

• 19. Jia et al., Atomistic first-principles modeling of single donor spin-qubit, Appl. Phys. Lett. 125, 184001 (2024).

• 10. 19. Smith et al., Ab initio calculation of energy levels for phosphorus donors in silicon, Sci. Rep. 7, 6010 (2017).

• QTCAD: Spin-qubit simulation package, https://docs.nanoacademic.com/qtcad/