\[sec:zfold\] to provide a clear contrast to the properties deriving from the $\delta$-doping of the silicon discussed in the paper.

Methodology {#sec:methods}
===========

Density functional theory calculations have been carried out using both plane-wave and LCAO basis sets. For the plane-wave (PW) basis set, the Vienna *ab initio* simulation package (<span style="font-variant:small-caps;">vasp</span>) [@Kresse99] software was used with projector augmented wave (PAW) [@Blochl94a; @Kresse99] pseudopotentials for Si and P. Due to the nature of the plane-wave (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.

Localized basis set calculations were performed using the Spanish Initiative for Electronic Simulations with Thousands of Atoms (<span style="font-variant:small-caps;">siesta</span>) [@Artacho08] software. In this case, the P and Si ionic cores were represented by norm-conserving Troullier-Martins pseudopotentials.[@Troullier91] The Kohn-Sham orbitals were expanded in the default single-$\zeta$ polarized (SZP) or double-$\zeta$ 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 $d$-type functions. These calculations were found to be converged for a mesh grid energy cutoff 