Density functional theory calculations have been carried out using both plane-wave and LCAO basis sets. For the PW basis set, the Vienna ab initio
simulation package (vasp
] software was used with projector augmented wave
] pseudo-potentials for Si and P. Due to the nature of the PW basis set, there exists a simple relationship between the cut-off energy and basis set completeness. For the structures considered in this work, the calculations were found to be converged for PW cut-offs of 450 eV.
Localised basis set calculations were performed using the Spanish Initiative for Electronic Simulations with Thousands of Atoms (siesta
] software. In this case, the P and Si ionic cores were represented by norm-conserving Troullier-Martins pseudo-potentials
]. The Kohn-Sham orbitals were expanded in the default single-ζ
polarized (SZP) or double-ζ
polarized (DZP) basis sets, which consist of 9 and 13 basis functions per atom, respectively. Both the SZP and DZP sets contain s
-, and d
-type functions. These calculations were found to be converged for a mesh grid energy cut-off of 300 Ry. In all cases, the generalized gradient approximation PBE
] exchange-correlation functional was used.
The lattice parameter for bulk Si was calculated using an eight-atom cell and found to be converged for all methods with a 12 × 12 × 12 Monkhorst-Pack (MP) k
]. The resulting values are presented in Table
and were used in all subsequent calculations.
Eight-atom cubic unit cell equilibrium lattice parameters for different methods used in this work
In modelling δ
-doped Si:P, as used in another work
], we adopted a tetragonal supercell description of the system, akin to those of other works
]. In accordance with the experiment, we inserted the P layer in a monatomic (001) plane as one atom in four to achieve 25% doping. This will henceforth be referred to as 1/4 monolayer (ML) doping. In this case, the smallest repeating in-plane unit had 4 atoms/ML (to achieve one in four dopings) and was a square with the sides parallel to the  and
10] directions. The square had a side length
), where a
is the simple cubic lattice constant of bulk silicon. The phosphorus layers had to be separated by a considerable amount of silicon due to the large Bohr radius of the hydrogen-like orbital introduced by P in Si (approximately 2.5 nm). Carter et al.
] showed that this far exceeded the sub-nanometre cell side length. If desired, cells with a lower in-plane density of dopants may be constructed by lengthening the cell in the x
directions, such that more Si atoms occupy the doped monolayer in the cell - though this would significantly increase the computational cost of such a calculation.
(001) Planar slice of the c(2×2) structure at the 1/4 ML doped monolayer. One of the Si sites has been replaced by a P atom (shown in dark gray). The periodic boundaries are shown in black.
A collection of tetragonal cells comprising 4, 8, 16, 32, 40, 60, 80, 120, 160 and 200 monolayers was constructed, having four atomic sites per monolayer and oriented with faces in the , [
10], and  directions (see Figure
). Cells used in PW calculations began at 4 layers and ran to 80 layers; larger cells were not computationally tractable with this method. SZP and DZP models began at 40 layers to overlap with PW for the converging region and were then extended to their tractable limit (200 and 160 layers, respectively) to study convergence past the capability of PW.
Figure 2 Ball and stick model of a δ-doped Si:P layer viewed along the  direction. Thirty-two layers in the  direction are shown. Si atoms (small gray spheres), P atoms (large dark gray spheres), covalent bonds (gray sticks), repeating cell boundary (more ...)
For tetragonal cells, the k
-point sampling was set as a 9 × 9 × N
Γ-centred MP mesh as we have found that failing to include Γ in the mesh can lead to the anomalous placement of the Fermi level on band structure diagrams. N
varied from 12 to 1 as the cells became more elongated (see Appendix 1). We note that, as mentioned in the work of Carter et al.
], the large supercells involved required very gradual (<0.1%) mixing of the new density matrix with the prior step, leading to many hundreds of self-consistent cycles before convergence was achieved.
Although it has been previously found that relaxing the positions of the nuclei gave negligible differences (<0.005 Å) to the geometry
], this was for a 12-layer cell and may not have included enough space between periodic repetitions of the doping plane for the full effect to be seen. Whilst a 40-layer model was optimised in the work of Carter et al.
], this made use of a mixed atom pseudo-potential and is not explicitly comparable to the models presented here. We have performed a test relaxation on a 40-layer cell using the PW basis (vasp
). The maximum subsequent ionic displacement was 0.05 Å, with most being an order of magnitude smaller. The energy gained in relaxing the cell was less than 37 meV (or 230 μ
eV/atom). We therefore regarded any changes to the structure as negligibly small, confirming the results of Carter et al.
], and proceeded without ionic relaxation.
Single-point energy calculations were carried out with both software programs; for vasp
, the electronic energy convergence criterion was set to 10−6
eV, and the tetrahedron method with Blöchl correction
] was used. For siesta
, a two-stage process was carried out: Fermi-Dirac electronic smearing of 300 K was applied in order to converge the density matrix within a tolerance of one part in 10−4
; the calculation was then restarted with the smearing of 0 K, and a new electronic energy tolerance criterion of 10−6
eV was applied (except for the 120- and 160-layer DZP models for which this was intractable; a tolerance of 10−4
eV was used in these cases). This two-stage process aided convergence as well as ensuring that the energy levels obtained were comparably accurate across methods. In addition, for each doped cell thus developed and studied, an undoped bulk Si cell of the same dimensions was constructed to aid in isolating those features primarily due to the doping.