Abstract
The differences in energy between electronic bands due to valley splitting are of paramount importance in interpreting transport spectroscopy experiments on stateoftheart quantum devices defined by scanning tunnelling microscope lithography. Using VASP, we develop a planewave density functional theory description of systems which is size limited due to computational tractability. Nonetheless, we provide valuable data for the benchmarking of empirical modelling techniques more capable of extending this discussion to confined disordered systems or actual devices. We then develop a less resourceintensive alternative via localised basis functions in SIESTA, retaining the physics of the planewave description, and extend this model beyond the capability of planewave methods to determine the ab initio valley splitting of wellisolated δlayers. In obtaining an agreement between planewave and localised methods, we show that valley splitting has been overestimated in previous ab initio calculations by more than 50%.
Keywords:
Density functional theory; Valley splitting; δDoped layers; Phosphorus in silicon; Basis sets; 73.22.f; 31.15.ae; 71.15.MbBackground
The study of the quantum properties of lowdimensional and doped structures is central to many nanotechnology applications [115]. Quantum devices in silicon have been the subject of concentrated recent interest, both experimental and theoretical, including the recent discussion of Ohm’s law at the nanoscale [16]. Efforts to make such devices have led to atomically precise fabrication methods which incorporate phosphorus atoms in a single monolayer of a silicon crystal [1720]. These dopant atoms can be arranged into arrays [21] or geometric patterns for wires [16,22] and associated tunnel junctions [23], gates, and quantum dots [24,25]  all of which are necessary components of a functioning device [26]. The patterns themselves define atomically abrupt regions of doped and undoped silicon. While silicon, bulkdoped silicon, and the physics of the phosphorus incorporation [27] are well understood, models of this quasitwodimensional phosphorus sheet are still in their initial stages. In particular, it is critical in many applications to understand the effect of this confinement on the conduction band valley degeneracy, inherent in the band structure of silicon. For example, the degeneracy of the valleys has the potential to cause decoherence in a spinbased quantum computer [28,29], and the degree of valley degeneracy lifting (valley splitting) defines the conduction properties of highly confined planar quantum dots [26].
The importance of understanding valley splitting in monolayer δdoped Si:P structures has led to a number of theoretical works in recent years, spanning several techniques, from pseudopotential theories via planar Wannier orbital bases [30], density functional theory (DFT) via linear combination of atomic orbital (LCAO) bases [31,32], to tightbinding models [3337] and effective mass theories (EMT) [3840]. We note that several of these papers are based upon the assumption that the effective masses of δdoped P in Si remain unchanged from bulkdoped values [38,39], an assumption which has been challenged [30,33]. Others assume doping over a multiatomic plane band [33,38] which no longer represents the state of the art in fabrication. There is currently little agreement between the valley splitting values obtained using these methods, with predictions ranging between 5 to 270 meV, depending on the calculational approach and the arrangement of dopant atoms within the δlayer. Density functional theory has been shown to be a useful tool in predicting how quantum confinement or doping perturbs the bulk electronic structure in silicon and diamondlike structures [4145]. The work of Carter et al. [31] represents the first attempt using DFT to model these devices by considering explicitly doped δlayers, using a localised basis set and the assumption that a basis set sufficient to describe bulk silicon will also adequately describe Pdoped Si. It might be expected, therefore, that the removal of the basis set assumption will lead to the best ab initio estimate of the valley splitting available, for a given arrangement of atoms. In the context of describing experimental devices, it is important to separate the effects of methodological choices, such as this, from more complicated effects due to physical realities, including disorder.
In this paper, we determine a consistent value of the valley splitting in explicitly δdoped structures by obtaining convergence between distinct DFT approaches in terms of basis set and system sizes. We perform a comparison of DFT techniques, involving localised numerical atomic orbitals and delocalised planewave (PW) basis sets. Convergence of results with regard to the amount of Si ‘cladding’ about the δdoped plane is studied. This corresponds to the normal criterion of supercell size, where periodic boundary conditions may introduce artificial interactions between replicated dopants in neighbouring cells. A benchmark is set via the delocalised basis for DFT models of δdoped Si:P against which the localised basis techniques are assessed. Implications for the type of modelling being undertaken are discussed, and the models extended beyond those tractable with planewave techniques. Using these calculations, we obtain converged values for properties such as band structures, energy levels, valley splitting, electronic densities of state and charge densities near the δdoped layer.
The paper is organised as follows: the ‘Methods’ section outlines the parameters used in our particular calculations; we present the results of our calculations in the ‘Results and discussion’ section and draw conclusions in the ‘Conclusions’ section. An elucidation of effects modifying the bulk band structure follows in Appendices 1 and 2 to provide a clear contrast to the properties deriving from the δdoping of the silicon discussed in the paper. The origin of valley splitting is discussed in Appendix 3.
Methods
Density functional theory calculations have been carried out using both planewave and LCAO basis sets. For the PW basis set, the Vienna ab initio simulation package (VASP) [46] software was used with projector augmented wave [46,47] pseudopotentials for Si and P. Due to the nature of the PW basis set, there exists a simple relationship between the cutoff energy and basis set completeness. For the structures considered in this work, the calculations were found to be converged for PW cutoffs of 450 eV.
Localised basis set calculations were performed using the Spanish Initiative for Electronic Simulations with Thousands of Atoms (SIESTA) [48] software. In this case, the P and Si ionic cores were represented by normconserving TroullierMartins pseudopotentials [49]. The KohnSham 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, p, and dtype functions. These calculations were found to be converged for a mesh grid energy cutoff of 300 Ry. In all cases, the generalized gradient approximation PBE [50] exchangecorrelation functional was used.
The lattice parameter for bulk Si was calculated using an eightatom cell and found to be converged for all methods with a 12 × 12 × 12 MonkhorstPack (MP) kpoint mesh [51]. The resulting values are presented in Table 1 and were used in all subsequent calculations.
Table 1. Eightatom cubic unit cell equilibrium lattice parameters for different methods used in this work
In modelling δdoped Si:P, as used in another work [26], we adopted a tetragonal supercell description of the system, akin to those of other works [30,31]. 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 inplane unit had 4 atoms/ML (to achieve one in four dopings) and was a square with the sides parallel to the [110] and 10] directions. The square had a side length (see Figure 1), 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 hydrogenlike orbital introduced by P in Si (approximately 2.5 nm). Carter et al. [31] showed that this far exceeded the subnanometre cell side length. If desired, cells with a lower inplane density of dopants may be constructed by lengthening the cell in the x and y 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.
Figure 1. (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 [110], [ 10], and [001] directions (see Figure 2). 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 [110] direction. Thirtytwo layers in the [001] direction are shown. Si atoms (small gray spheres), P atoms (large dark gray spheres), covalent bonds (gray sticks), repeating cell boundary (solid line).
For tetragonal cells, the kpoint 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. [32], the large supercells involved required very gradual (<0.1%) mixing of the new density matrix with the prior step, leading to many hundreds of selfconsistent 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 [31], this was for a 12layer cell and may not have included enough space between periodic repetitions of the doping plane for the full effect to be seen. Whilst a 40layer model was optimised in the work of Carter et al. [32], this made use of a mixed atom pseudopotential and is not explicitly comparable to the models presented here. We have performed a test relaxation on a 40layer 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. [31,32], and proceeded without ionic relaxation.
Singlepoint 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 [52] was used. For SIESTA, a twostage process was carried out: FermiDirac 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 160layer DZP models for which this was intractable; a tolerance of 10^{−4} eV was used in these cases). This twostage 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.
Results and discussion
Analysis of band structure
Once converged charge densities were obtained, band structures were calculated along the M–Γ–X highsymmetry pathway (as shown in Appendix 1), using at least 20 kpoints between highsymmetry points. For comparative purposes, the band structures have all been aligned at the valence band maximum (VBM).
Figure 3 contrasts the bulk and doped band structures for the 40layer PW calculation. DZP and SZP results are qualitatively similar on this scale, albeit with different band energies in the SZP model, and are omitted in the interest of clarity in the diagram. As discussed in Appendix 2, it is evident from the bulk values that the elongated cells have led to the folding of two conduction band minimum valleys towards the Γ point. Also visible is the difference that the doping potential makes to the system; what was the lowest unoccupied orbital (Γ_{1} band) in the bulk is now dragged down in energy by the extra ionic potential. It is of note that the region near Γ, corresponding to the k_{z }valleys which can be modelled as having different effective masses to the k_{x,y }valleys, [30] is brought lower than the region corresponding to the k_{x,y} valleys and is nondegenerate. The second (Γ_{2}) band behaves in a similar fashion. The third (δ) band appears to maintain a minimum away from the Γ point in the Σ_{TET} direction (which is equivalent to the Δ_{FCC }direction; see Appendix 1) but in a less parabolic fashion than the lower two; its minimum is similar to the value at Γ. This band is nondegenerate along this particular direction in kspace, but due to the supercell symmetry, it is actually fourfold degenerate, in contrast to the other bands.
Figure 3. Full band structure (colour online) of the 40layer tetragonal system calculated using PW (VASP). Bulk band structure (shaded gray background), doped band structure (solid black) and Fermi level (labelled solid red).
The Fermi level for the doped system is also shown, clearly being crossed by all three of these bands which are therefore able to act as open channels for conduction.
As mentioned above, the band structures are similar across all methods, but upon detailed inspection, important differences come to light. A closer look at the δ band shows a qualitative difference between the predictions using SZP (Figure 4c) and the PW and DZP results (Figure 4a,b): the models with a more complete basis predict the band minimum to occur in the Σ_{TET}(Δ_{FCC}) direction, below the value at Γ, while the SZP band structure shows the reverse  the minimum at Γ, a similar amount below a secondary minimum in the Σ_{TET }direction, a qualitative difference.
Figure 4. Band structure (colour online) of the 40layer tetragonal system zoomed in on the δ band. (a) PW (VASP), (b) DZP (SIESTA) and (c) SZP basis sets were used. Fermi level is shown by a solid horizontal red line.
The difference between the energies of the first two band minima (Γ_{1}− Γ_{2}, illustrated in Figure 5), or the valley splitting, from the PW and DZP calculations, agrees with each other to within ∼6 meV. Significantly, the value obtained using our SZP basis set differs by 52 meV, some 55% larger than the value obtained using the PW basis set. The importance of this discrepancy cannot be overstated; valley splitting is directly relatable to experimentally observable resonances in transport spectroscopy of devices made with this δdoping technology (see [26]).
Figure 5. Minimum band energies for tetragonal systems with 1/4 ML doping. (a) PW (VASP), (b) DZP (SIESTA) and (c) SZP (SIESTA) basis sets were used. Fermi level also shown where appropriate. Bold numbers indicate energy differences between band minima.
In the smallest cells (<16 layers), less than three bands are observed. This is likely due to the lack of cladding in the z direction, leading to a significant interaction between the dopant layers, raising the energy of each band. Whilst the absolute energy of each level still varies somewhat, even with over 100 layers incorporated, we find that the Γ_{1}–Γ_{2 }values are well converged with 80 layers of cladding for all methods (see Figure 5). Indeed, they may be considered reasonably converged even at the 40layer level (0.5 meV or less difference to the largest models considered). The differences between the energies of the second and third band minima (Γ_{2}–δ splittings) are also shown in Figure 5 and show good convergence (within 1 meV) for cells of 80 layers or larger.
The Fermi level follows a similar pattern to the Γ and δlevels. In particular, the gap between the Fermi level and Γ_{1 }level does not change by more than 1 meV from 60 to 160 layers.
Given that the properties of interest are the differences between the energy levels, rather than their absolute values (or position relative to the valence band), in the interest of computational efficiency, we observe that using the DZP basis with 80 layers of cladding is sufficient to achieve consistent, converged results.
Valley splitting
Table 2 summarises the valley splitting values of 1/4 ML Pdoped silicon obtained using different techniques, showing a large variation in the actual values. In order to make sense of these results, it is important to note two major factors that affect valley splitting: the doping method and the arrangement of phosphorus atoms in the δlayer. As the results from the work of Carter et al. [32] show, the use of implicit doping causes the valley splitting value to be much smaller than in an explicit case (∼7 meV vs. 120 meV). It is also shown that the use of random P coverage on the δlayer reduces the valley splitting value by only 40 to 50 meV compared with the fully ordered placement, leaving a large discrepancy between the valley splitting results from implicit and explicit doping. This large decrease in valley splitting due to implicit doping can be explained by the smearing of the doping layer in the direction normal to the δlayer, thereby decreasing the quantum confinement effect responsible for breaking the degeneracy in the system. Carter et al. [32] also shows that the arrangement of the phosphorus atoms in the δlayer strongly influences the valley splitting value. In particular, they showed that there is a difference of up to 220 meV between P doping along the [110] direction and along the [100] direction. It should be noted, however, that deterministic nearestneighbour donor placements are not yet physically realisable due to the P incorporation mechanism currently employed [27,53]. Similarly, the perfectly ordered arrangement discussed here is highly improbable, given the experimental limitations, but represents the ideal case from which effects such as disorder can be studied.
Table 2. Valley splitting values of 1/4 ML Pdoped silicon obtained using different techniques
Our results show that valley splitting is highly sensitive to the choice of basis set. Due to the nature of PW basis set, it is straightforward to improve its completeness by increasing the planewave cutoff energy. In this way, we establish the most accurate valley splitting value within the context of density functional theory. Using this benchmark value, we can then establish the validity and accuracy of other basis sets, which can be used to extend the system sizes to that beyond what is practical using a PW basis set. As seen in Table 2, the valley splitting value converges to 93 meV using 80layer cladding. The DZP localised basis set gives an excellent agreement at 99.5 meV using 80layer cladding (representing a 7% difference). On the other hand, our SZP localised basis set gave a value of 145 meV using the same amount of cladding. This represents a significant difference of 55% over the value obtained using PW basis set and demonstrates that SZP basis sets are unsuitable for accurate determination of valley splitting in these systems.
Density of states
The electronic density of states (eDOS) was calculated for each cell. Figure 6 compares the unscaled eDOS for bulk 80layer cells to that of doped cells varying from 40 to 80 layers. The bulk bandgap is visible, with the conduction band rising sharply to the right of the figure. The doped eDOS exhibits density in the bulk bandgap, although the features of the spectra differ slightly according to the basis set used.
Figure 6. Electronic densities of states for tetragonal systems with 0 and 1/4 ML doping. The DZP (SIESTA) basis set was used. The Fermi level is indicated by a solid vertical line with label, and 50meV smearing was applied for visualization purposes.
The Fermi energy exhibits convergence with respect to the amount of cladding, as reported above. It is also notable that the eDOS within the bandgap are nearly identical regardless of the cell length (in z). This indicates that layerlayer interactions are negligibly affecting the occupied states and, therefore, that the applied ‘cladding’ is sufficient to insulate against these effects.
Electronic width of the plane
In order to quantify the extent of the donorelectron distribution, we have integrated the local density of states between the VBM and Fermi level and have taken the planar average with respect to the zposition. Figure 7 shows the planar average of the donor electrons (a sum of both spinup and spindown channels) for the 80layer cell calculated using the DZP basis set. After removing the small oscillations related to the crystal lattice to focus on the physics of the δlayer, by Fourier transforming, a Lorentzian function was fitted to the distribution profile. (Initially, a threeparameter Gaussian fit similar to that used in [40] was tested, but the Lorentzian gave a better fit to the curve.)
Figure 7. Planar average of donorelectron density as a function of zposition for 1/4 MLdoped 80layer cell. The DZP basis set was used. The fitted Lorentzian function is also shown.
Table 3 summarises the maximum donorelectron density and the full width at half maximum (FWHM) for the 1/4 MLdoped cells, each calculated from the Lorentzian fit. Both of these properties are remarkably consistent with respect to the number of layers, indicating that they have converged sufficiently even at 40 layers.
Table 3. Calculated maximum donorelectron density,ρ_{max}, and FWHM
Our results differ from a previous DFT calculation [32] which cited an FWHM of 5.62 Å for a 1/4 MLdoped, 80layer cell calculated using the SZP basis set (and 10 × 10 × 1 kpoints). We note that those values were taken from the unfitted, untransformed donorelectron distribution and represent an approximately 15% underestimation in comparison with the DZP result. The peak height is not shown in the work of Carter et al. [32], but the value from another work [31] (1.7 × 10^{21} e/cm^{3}) is a factor of 0.44 smaller than the peak we observe here. This may be due, to some extent, to the larger width of the SZP model leading to an associated lowering of the peak density.
Conclusions
In this article, we have studied the valley splitting of the monolayer δdoped Si:P, using a density functional theory model with a planewave basis to establish firm grounds for comparison with less computationally intensive localisedbasis ab initio methods. We found that the description of these systems (by density functional theory, using SZP basis functions) overestimates the valley splitting by over 50%. We show that DZP basis sets are complete enough to deliver values within 10% of the planewave values and, due to their localised nature, are capable of calculating the properties of models twice as large as is tractable with planewave methods. These DZP models are converged with respect to size well before their tractable limit, which approaches that of SZP models.
Valley splittings are important in interpreting transport spectroscopy experiment data, where they relate to families of resonances, and in benchmarking other theoretical techniques more capable of actual device modelling. It is therefore pleasing to have an ab initio description of this effect which is fully converged with respect to basis completeness as well as the usual size effects and kpoint mesh density.
We have also studied the band structures with all three methods, finding that the DZP correctly determines the δband minima away from the Γ point, where the SZP method does not. We show that these minima occur in the Σ direction for the type of cell considered, not the δ direction as has been previously reported. Having established the DZP methodology as sufficient to describe the physics of these systems, we then calculated the electronic density of states and the electronic width of the δlayer. We found that previous SZP descriptions of these layers underestimate the width of the layers by almost 15%.
We have shown that the properties of interest of δdoped Si:P are well converged for 40layer supercells using a DZP description of the electronic density. We recommend the use of this amount of surrounding silicon, and technique, in any future DFT studies of these and similar systems  especially if interlayer interactions are to be minimised.
Appendix 1
Subtleties of bandstructure
Regardless of the type of calculation being undertaken, a band structure diagram is inherently linked to the type (shape and size) of cell being used to represent the system under consideration. For each of the 14 Bravais lattices available for threedimensional supercells, a particular Brillouin zone (BZ) with its own set of highsymmetry points exists in reciprocal space [54]. Similarly, each BZ has its own set of highsymmetry directions. Some of these BZs share a few highsymmetry point labels (or directions), such as X or L (δ or Σ), and they all contain Γ, but these points are not always located in the same place in reciprocal space.
A simple effect of this can be seen by increasing the size of a supercell. This has the result of shrinking the BZ and the coordinates of highsymmetry points on its boundary by a corresponding factor. Consider the conduction band minimum (CBM) found at the δ valley in the Si conduction band. This is commonly located at in the δ direction towards X (also Y, Z and their opposite directions). Should we increase the cell by a factor of 2, the BZ will shrink (BZ→BZ’), placing the valley outside the new BZ boundary (past X’); however, a valid solution in any BZ must be a solution in all BZs. This results in the phenomenon of band folding, whereby a band continuing past a BZ boundary reenters the BZ on the opposite side. Since the X direction in a facecentred cubic (FCC) BZ is sixfold symmetric, a solution near the opposite BZ boundary is also a solution near the one we are focussing on. This results in the appearance that the band continuing past the BZ boundary is ‘reflected’, or folded, back on itself into the first BZ. Since the new BZ boundary in this direction is now at , the location of the valley will be at , as mentioned in the work of Carter et al. [31]. Each further increase in the size of the supercell will result in more folding (and a denser band structure). Care is therefore required to distinguish between a new band and one which has been folded due to this effect when interpreting band structure.
Continuing with our example of silicon, whilst the classic band structure [55] is derived from the bulk Si primitive FCC cell (containing two atoms), it is often more convenient to use a simple cubic (SC) supercell (eight atoms) aligned with the 〈100〉 crystallographic directions. In this case, we experience some of the common labelling; the δ direction is defined in the same manner for both BZs, although we see band folding (in a similar manner to that discussed previously) due to the size difference of the reciprocal cells (see Figure 8). We also see a difference in that, although the Σ direction is consistent, the points at the BZ boundaries have different symmetries and, therefore, label (K_{FCC}, M_{SC}). (The L_{FCC} point and ⋀ _{FCC }direction have no equivalent for tetragonal cells, and hence, we do not consider band structure in that direction here).
Figure 8. Band structure and physical structure of FCC and SC cells. (a) Typical band structure of bulk Si for twoatom FCC (solid lines) and eightatom SC cells (dotted lines with squares), calculated using the VASP planewave method (see ‘Methods’ section). (b) Twoatom FCC cell. (c) Eightatom SC cell.
Consider now the δdoping case discussed in the ‘Methods’ section, where we wish to align our cell with the [110] and [ 10] directions (by rotating the cell 45° anticlockwise about z; this will also require a resizing of the cell in the plane to maintain periodicity  see Figure 9), to allow us to include precisely four atoms per monolayer (as required for the minimal representation of 1/4 ML doping). We now have a situation where the X_{TET }point in the new tetragonal BZ (see Figure 10) is no longer in the direction of the X_{SC} point in the simple cubic BZ, despite both X points being in the centre of a face of their BZ. Due to the rotation, what used to be the Δ_{SC} direction in the simple cubic BZ is now the Σ_{TET} direction (pointing towards M at the corner of the BZ in the k_{z }= 0 plane) in the tetragonal BZ. The tetragonal CBM, while physically still the same as the CBM in the FCC or simple cubic BZ, is not represented in the same fashion (see Figure 11).
Figure 9. Geometrical difference between the simple cubic and tetragonal cells. A (001) planar cut through an atomic monolayer is shown.
Figure 10. The Brillouin zone for a tetragonal cell. The M–Γ–X path used in this work is shown.
Figure 11. Band structure (colour online) diagram for tetragonal bulk Si structures with increasing number of layers. The VASP plane wave method was used (see ‘Methods’ section).
Appendix 2
Band folding in the z direction
Increasing the z dimension of the cell leads to successive folding points being introduced as the BZ shrinks along k_{z }(see Appendix 1). This has the effect of shifting the conduction band minima in the ± k_{z} directions closer and closer to the Γ point (see Figure 8a) and making the band structure extremely dense when plotting along k_{z}. This results in the value of the lowest unoccupied eigenstate at Γ being lowered as what were originally other sections of the band are successively mapped onto Γ, and after a sufficient number of folds, the value at Γ is indistinct from the original CBM value. The effects of this can be seen in Table 4, which describes increasingly elongated tetragonal cells of bulk Si. When we then plot the band structure in a different direction, e.g. along k_{x}, the translation of the minima from ± k_{z }onto the Γ point appear as a new band with twofold degeneracy. The degeneracy of the original band seems to drop from six to fourfold, in line with the reduced symmetry (we only explicitly calculate one, and the other three occur due to symmetry considerations). This is half of the origin of the ‘Γbands’ (more details are presented in Appendix 3). Once the k_{z} valleys are sited at Γ, parabolic dispersion corresponding to the transverse kinetic energy terms is observed along k_{x} and k_{y}, at least close to the band minimum (see Figure 11)  in contrast to the four ‘δbands’ whose dispersion (again parabolic) is governed by the longitudinal kinetic energy terms. The different curvatures are related to the different effective masses (transverse, longitudinal) of the silicon CBM. It should be noted that the bands are still degenerate in energy at this stage  their minima (and range) occur at (over) the same energy (energies) even though their projections onto the k_{x} axis are different.
Table 4. Energy levels of tetragonal bulk Si structures
All methods considered in Table 4 show the LUMO at Γ (folded in along ± k_{z}) approaching the CBM value as the amount of cladding increases; at 80 layers, the LUMO at Γ is within 1 meV of the CBM value. It is also of note that the PW indirect bandgap agrees well with the DZP value and less so with the SZP model. This is an indication that, although the behaviour of the LUMO with respect to the cell shape is well replicated, the SZP basis set is demonstrably incomplete. Conversely, pairwise comparisons between the PW and DZP results show agreement to within 5 meV.
It is important to distinguish effects indicating convergence with respect to cladding for doped cells (i.e. elimination of layerlayer interactions) from those mentioned previously derived from the shape and size of the supercell. Strictly, the convergence (with respect to the amount of encapsulating Si) of those results we wish to study in detail, such as the differences in energy between occupied levels in what was the bulk bandgap, provides the most appropriate measure of whether sufficient cladding has been applied.
Appendix 3
Valley splitting
Here, we discuss the origins of valley splitting, in the context of phosphorus donors in silicon. Following on from the discussion of Si band minima in Appendices 1 and 2, we have, via elongation of the supercell and consequent band folding, a situation where, instead of the sixfold degeneracy (due to the underlying symmetries of the Si crystal lattice), we see an apparent splitting of these states into two groups (6 → 2 + 4, or 2 Γ + 4 δ minima).
We now consider what happens in perfectly ordered δdoped monolayers, as per the main text. Here, we break the underlying Si crystal lattice symmetries by including foreign elements in the lattice. By placing the donors regularly (according to the original Si lattice pattern) in one [001] monolayer, we reduce the symmetry of the system to tetragonal, with the odd dimension being transverse to the plane of donors. This dimension can be periodic (as in the supercells described earlier), infinite (as in the EMT model of Drumm et al. [40]) or extremely long on the atomic scale (as the experiments are).
Immediately, therefore, we expect the same apparent 2 + 4 breaking of the original sixfold degenerate conduction band minima. Of course, as we have introduced phosphorus (which has one more electron and one more proton than silicon), this next band (still actually sixfold degenerate in bulk silicon) will be occupied and will now be influenced by the new potential. The subbands interact differently with the potential, thanks to the different curvatures in their dispersion relations and drop by different amounts into the bandgap. As discussed in detail in Drumm et al. [40], the filling of these subbands is partial rather than complete (or absent) and is governed by both the energy of their minima and their respective effective masses. We now have an actual breaking of the sixfold degeneracy into a true 2 + 4 system.
If we still look closer, we might expect these lower degeneracies to spontaneously break  nature, after all, is said to abhor degeneracy. Indeed, this does occur, but for this special case of δdoped Si:P, the effect is enhanced by the strong Vshaped potential about the monolayer due to the extra charge in the donor nuclei [40]. Consideration of odd and even solutions to the effective mass Schrödinger equation for this subband leads to their superposition(s) and subsequent energy difference. This is enhanced further in the KohnSham formalism, as evidenced in previous sections. (The four δ minima also split but on a farreduced scale not visible using current DFT techniques.) We thus expect, in the DFT picture, to see 6 →2 + 4→1 + 1 + 4 subband structure, namely the Γ_{1}, Γ_{2} and δ bands. The valley splitting which is the main focus of this paper is the energy difference between the Γ_{1} and Γ_{2} band minima due to the superposition of solutions.
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
DWD, SPR, and LCLH conceived the study. Density functional theory calculations were carried out by DWD, AB, and MCP. All authors contributed to the discussion of results and drafting of the the final manuscript. All authors read and approved the final manuscript.
Acknowledgements
The authors acknowledge funding by the ARC Discovery grant DP0986635. This research was undertaken on the NCI National Facility in Canberra, Australia, which is supported by the Australian Commonwealth Government. We thank Oliver Warschkow, Damien Carter and Nigel Marks for their feedback on our manuscript.
References

Shen G, Chen D: Onedimensional nanostructures and devices of IIV group semiconductors.
Nanoscale Res Lett 2009, 4(8):779788. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Dresselhaus MS, Chen G, Tang MY, Yang R, Lee H, Wang D, Ren Z, Fleurial JP, Gogna P: New directions for Lowdimensional thermoelectric materials.
Adv Mater 2007, 19:10431053. Publisher Full Text

Lu YH, Hong ZX, Feng YP, Russo SP: Roles of carbon in light emission of ZnO.
Appl Phys Lett 2010, 96(9):091914. Publisher Full Text

Zhao YS, Fu H, Peng A, Ma Y, Xiao D, Yao J: Lowdimensional nanomaterials based on small organic molecules: preparation and optoelectronic properties.
Adv Mater 2008, 20:28592876. Publisher Full Text

Drumm DW, Per MC, Russo SP, Hollenberg LCL: Thermodynamic stability of neutral Xe defects in diamond.

Tsu R: Superlattices: problems and new opportunities, nanosolids.
Nanoscale Res Lett 2011, 6:127. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Luo DS, Lin LH, Su YC, Wang YT, Peng ZF, Lo ST, Chen KY, Chang YH, Wu JY, Lin Y, Lin SD, Chen JC, Huang CF, Liang CT: A deltadoped quantum well system with additional modulation doping.
Nanoscale Res Lett 2011, 6:139. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Webber BT, Per MC, Drumm DW, Hollenberg LCL, Russo SP: Ab initio thermodynamics calculation of the relative concentration of NV and NV0 defects in diamond.

Conibeer G, PerezWurfl I, Hao X, Di D, Lin D: Si solidstate quantum dotbased materials for tandem solar cells.
Nanoscale Res Lett 2012, 7:193. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Dick R: Interdimensional effects in nanostructures.
Nanoscale Res Lett 2012, 7(1):581. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Budi A, Drumm DW, Per MC, Tregonning A, Russo SP, Hollenberg LCL: Electronic properties of multiple adjacent δdoped Si:P layers: the approach to monolayer confinement.

Sun HH, Guo FY, Li DY, Wang L, Wang DB, Zhao LC: Intersubband absorption properties of high Al content Al(x)Ga(1−x)N/GaN multiple quantum wells grown with different interlayers by metal organic chemical vapor deposition.
Nanoscale Res Lett 2012, 7:649. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

De Padova P, Ottaviani C, Ronci F, Colonna S, Olivieri B, Quaresima C, Cricenti A, Dávila ME, Hennies F, Pietzsch A, Shariati N, Le Lay G: Mnsilicide nanostructures aligned on massively parallel silicon nanoribbons.
J Phys: Condens Matter 2013, 25:014009. Publisher Full Text

Barnard AS, Russo SP, Snook IK: Ab initio modelling of B and N in C29 and C29H24 nanodiamond.
J Chem Phys 2003, 118:1072510728. Publisher Full Text

Erogbogbo F, Liu X, May JL, Narain A, Gladding P, Swihart MT, Prasad PN: Plasmonic gold and luminescent silicon nanoplatforms for multimode imaging of cancer cells.
Integr Biol 2013, 5:144. Publisher Full Text

Weber B, Mahapatra S, Ryu H, Lee S, Fuhrer A, Reusch TCG, Thompson DL, Lee WCT, Klimeck G, Hollenberg LCL, Simmons MY: Ohm’s law survives to the atomic scale.
Science 2012, 335:64. PubMed Abstract  Publisher Full Text

Tucker JR, Shen TC: Prospects for atomically ordered device structures based on STM lithography.
SolidState Electron 1998, 42:1061. Publisher Full Text

O’Brien JL, Schofield SR, Simmons MY, Clark RG, Dzurak AS, Curson NJ, Kane BE, McAlpine NS, Hawley ME, Brown GW: Towards the fabrication of phosphorus qubits for a silicon quantum computer.

Shen TC, Ji JY, Zudov MA, Du RR, Kline JS, Tucker JR: Ultradense phosphorus delta layers grown into silicon from PH3 molecular precursors.
Appl Phys Lett 2002, 80:1580. Publisher Full Text

Fuechsle M, Ruess FJ, Reusch TCG, Mitic M, Simmons MY: Surface gate and contact alignment for buried, atomically precise scanning tunneling microscopypatterned devices.
J Vac Sci Technol B 2007, 25:2562. Publisher Full Text

Pok W, Reusch TCG, Scappucci G, Ruess FJ, Hamilton AR, Simmons MY: Electrical characterization of ordered Si:P dopant arrays.

Ruess FJ, Goh KEJ, Butcher MJ, Reusch TCG, Oberbeck L, Weber B, Hamilton AR, Simmons MY: Narrow, highly Pdoped, planar wires in silicon created by scanning probe microscopy.
Nanotechnology 2007, 18:044023. Publisher Full Text

Ruess FJ, Pok W, Goh KEJ, Hamilton AR, Simmons MY: Electronic properties of atomically abrupt tunnel junctions in silicon.

Ruess FJ, Pok W, Reusch TCG, Butcher MJ, Goh KEJ, Oberbeck L, Scappucci G, Hamilton AR, Simmons MY: Realization of atomically controlled dopant devices in silicon.
Small 2007, 3:563. PubMed Abstract  Publisher Full Text

Fuhrer A, Füchsle M, Reusch TCG, Weber B, Simmons MY: Atomicscale, all epitaxial inplane gated donor quantum dot in silicon.
Nano Lett 2009, 9:707. PubMed Abstract  Publisher Full Text

Fuechsle M, Mahapatra S, Zwanenburg FA, Friesen M, Eriksson MA, Simmons MY: Spectroscopy of fewelectron singlecrystal silicon quantum dots.
Nature Nanotechnology 2010, 5:502. PubMed Abstract  Publisher Full Text

Wilson HF, Warschkow O, Marks NA, Schofield SR, Curson NJ, Smith PV, Radny MW, McKenzie DR, Simmons MY: Phosphine dissociation on the Si(001) surface.
Phys Rev Lett 2004, 93:226102. PubMed Abstract  Publisher Full Text

Koiller B, Hu X, Das Sarma S: Exchange in siliconbased quantum computer architecture.

Boykin TB, Klimeck G, Friesen M, Coppersmith SN, von Allmen P, Oyafuso F, Lee S: Valley splitting in lowdensity quantumconfined heterostructures studied using tightbinding models.

Qian G, Chang YC, Tucker JR: Theoretical study of phosphorus δdoped silicon for quantum computing.

Carter DJ, Warschkow O, Marks NA, McKenzie DR: Electronic structure models of phosphorus δdoped silicon.

Carter DJ, Marks NA, Warschkow O, McKenzie DR: Phosphorus δdoped silicon: mixedatom psuedopotentials and dopant disorder effects.
Nanotechnology 2011, 22:065701. PubMed Abstract  Publisher Full Text

Cartoixa X, Chang YC: Fermilevel oscillation in ntype δdoped Si: a selfconsistent tightbinding approach.

Lee S, Ryu H, Klimeck G, Jiang Z: Million atom electronic structure and device calculations on petascale computers. In Proc. of the 13th Int. Workshop on Computational Electronics. Tsinghua University, Beijing; Unknown Month 27.
vol 10. Piscataway: IEEE; 2009. doi:
Publisher Full Text 
Ryu H, Lee S, Weber B, Mahapatra S, Simmons MY, Hollenberg LCL, Klimeck G: Quantum transport in ultrascaled phosphorusdoped silicon nanowires. In Proceedings of the 2010 IEEE Silicon Nanoelectronics Workshop, Honolulu, USA, 1314 June 2010. Piscataway: IEEE; 2010.
doi:
Publisher Full Text 
Ryu H, Lee S, Klimeck G: A study of temperaturedependent properties of Ntype δdoped Si bandstructures in equilibrium. In Proc. of the 13th Int. Workshop on Computational Electronics. Tsinghua University, Beijing; Unknown Month 27.
vol 10. Piscataway: IEEE; 2009. doi:
Publisher Full Text 
Lee S, Ryu H, Campbell H, Hollenberg LCL, Simmons MY, Klimeck G: Electronic structure of realistically extended atomistically resolved disordered Si:P δdoped layers.

Scolfaro LMR, Beliaev D, Enderlein R, Leite JR: Electronic structure of ntype δdoping multiple layers and superlattices in silicon.
Phys Rev B 1994, 50:8699. Publisher Full Text

RodriguezVargas I, GaggeroSager LM: Subband and transport calculations in double ntype δdoped quantum wells in Si.
J Appl Phys 2006, 99:033702. Publisher Full Text

Drumm DW, Hollenberg LCL, Simmons MY, Friesen M: Effective mass theory of monolayer δ doping in the highdensity limit.

Delley B, Steigmeier EF: Quantum confinement in Si nanocrystals.

Delley B, Steigmeier EF: Size dependence of band gaps in silicon nanostructures.
Appl Phys Lett 1995, 67:2370. Publisher Full Text

Ramos LE, Teles LK, Scolfaro LMR, Castineira JLP, Rosa AL, Leite JR: Structural, electronic, and effectivemass properties of silicon and zincblende groupIII nitride semiconductor compounds.

Zhou ZY, Brus L, Friesner R: Electronic structure and luminescence of 1.1 and 1.4nm silicon nanocrystals: oxide shell versus hydrogen passivation.
Nano Lett 2003, 3:163. Publisher Full Text

Barnard AS, Russo SP, Snook IK: Ab initio modelling of band states in doped diamond.
Philos Mag 2003, 83:1163. Publisher Full Text

Kresse G, Joubert D: From ultrasoft pseudopotentials to the projector augmentedwave method.
Phys Rev B 1999, 59:1758. Publisher Full Text

Blöch PE: Projector augmentedwave method.
Phys Rev B 1994, 50:17953. Publisher Full Text

Artacho E, Anglada E, Dieguez O, Gale JD, Garcia A, Junquera J, Martin RM, Ordejon P, Pruneda JM, SanchezPortal D, Soler JM: The siesta method; developments and applicability.
J Phys Condens Matter 2008, 20:064208. PubMed Abstract  Publisher Full Text

Troullier N, Martins JL: Efficient pseudopotentials for planewave calculations.

Perdew JP, Burke K, Ernzerhof M: Generalized gradient approximation made simple.
Phys Rev Lett 1996, 77:3865. PubMed Abstract  Publisher Full Text

Monkhorst HJ, Pack JD: Special points for Brillouinzone integrations.
Phys Rev B 1976, 13:5188. Publisher Full Text

Blöchl PE, Jepsen O, Andersen OK: Improved tetrahedron method for Brillouinzone integrations.
Phys Rev B 1994, 49:16223. Publisher Full Text

Wilson HF, Warschkow O, Marks NA, Curson NJ, Schofield SR, Reusch TCG, Radny MW, Smith PV, McKenzie DR, Simmons MY: Thermal dissociation and desorption of PH3 on Si(001): a reinterpretation of spectroscopic data.

Bradley CJ, Cracknell JP: The Mathematical Theory of Symmetry in Solids: Representation Theory for Point Groups and Space Groups. Oxford: Clarendon Press; 1972.

Chelikowsky JR, Cohen ML: Electronic structure of silicon.
Phys Rev B 1974, 10:5095. Publisher Full Text