Half-filled orbital and unconventional geometry of a common dopant in Si(001)

The determining factor of the bulk properties of doped Si is the column rather than the row in the periodic table from which the dopants are drawn. It is unknown whether the basic properties of dopants at surfaces and interfaces, steadily growing in importance as microelectronic devices shrink, are also solely governed by their column of origin. The common light impurity P replaces individual Si atoms and maintains the integrity of the dimer superstructure of the Si(001) surface, but loses its valence electrons to surface states. Here we report that isolated heavy dopants are entirely different: Bi atoms form pairs with Si vacancies, retain their electrons and have highly localized, half-filled orbitals.


I. INTRODUCTION
Semiconductors such as silicon (Si) derive remarkable utility from the ease with which they are converted from insulating to conducting behaviour via substitution of dopants from neighbouring columns in the periodic table. Isolated phosphorus (P) and bismuth (Bi) atoms, both with one more valence electron than Si, substitute for Si atoms, and behave approximately as hydrogenic atoms with radii and binding energies described by the Bohr theory where the electron mass and dielectric constant  are those provided by the diamond lattice of Si rather than the vacuum.
Dopants in the bulk are well understood, but their properties at surfaces are poorly understood; for instance, we have recently shown that heavy dopants at Si(111) can induce ambipolar charge states 1 . There is growing interest in heavier dopants of conventional semiconductors because of their more tightly bound outer electrons, which increases operating temperatures for applications ranging from THz impurity lasers to quantum computers exploiting their spin 2 . We have accordingly chosen the heaviest column V element for the current study. Here we report that isolated heavy dopants in the Si(001) surface are entirely different to light dopants: Bi atoms form pairs with Si vacancies, retain their electrons and have highly localized, half-filled orbitals. Thus, we finally have a dopant of silicon with a half-filled orbital at the surface, enabling both new applications in spintronics, as well as Si surface-based quantum many-body physics, hitherto confined to the more complex Si(111) surface.
Scanning tunnelling microscopes (STMs), where a metallic tip is positioned to sub-Angstrom precision using piezoelectric drives and the tunnelling current between tip and sample is directly measured and/or used to control the drives, are uniquely suited to determine the electronic structure of surfaces. Here we employ an Omicron low-temperature (LT)-STM, at 5 K in ultrahigh vacuum (UHV) with a base pressure of ∼ 4×10 −9 Pa. For the Si(001) surface, strong tip-surface interactions can induce a complicated, unstable structural change between arXiv:1307.1984v1 [cond-mat.mtrl-sci] 8 Jul 2013 c(4×2) and p(2×2) reconstruction at low temperature 4,5 . To avoid this, we ensure a relatively large tip-sample separation during measurements by using small tunnelling currents, typically 100 pA for topographical imaging.
The surface was prepared in UHV by a standard flashing (direct current heating) method at ∼1100 • C for 10 seconds. The pressure was kept below 1× 10 −8 Pa during the flashing to prepare a clean surface. When searching for features resulting from Bi dopants, we must first eliminate other features. The electronic structure of standard defects in a clean Si(001) surface have already been studied by STM [6][7][8] , and so such defects can be identified and eliminated from our study by bias-voltagedependent imaging. The native B dopant might also be observed; however, it has a much lower density than the ion-implanted Bi impurities, and bias-dependent STM images of a single B dopant in Si(001) surface have already been reported 9 , so also allowing B impurities to be identified, if present, and excluded from the current study.
The electronic structure calculations used density functional theory (DFT) within the PBE generalized gradient approximation (GGA) 12 . The system was modelled with a slab at least six layers deep with the last two layers fixed in bulk-like positions and the bottom layer terminated in hydrogen. The Brillouin zone was sampled with the same density for all cells at a level which gave meV convergence in total energies. We used two different basis sets and codes: plane waves and gaussians.
The plane-wave basis calculations as implemented in the VASP code 10 used a cutoff energy of 200 eV, which converged energy differences to meV accuracy. We used ultrasoft pseudopotentials 11 and spin-polarized density functional theory, with the PBE GGA functional 12 for the exchange-correlation terms. After performing spinpolarized calculations we also tested the effect of including spin-orbit interactions, which made only a small difference to the electronic structure. For calculations on a single dopant, we chose the cell sizes to vary the spacing between its periodic images. We used cells with between 1 and 3 dimer rows, giving lateral sizes of 0.77 nm, 1.54 nm and 2.30 nm, with four to eight dimers in the cell, giving lengths of 1.54 nm and 3.07 nm. The basic density of k-points was set by the 1.54 nm × 1.54 nm cell, which we sampled with a mesh of 4 × 4 × 1 k-points; other cells had the number of k-points varied to achieve the same density, and convergence checks were made on all cells. Structural relaxations were carried out until all components of all forces were less than 0.25 eV/nm. Exact exchange is relatively inefficient to implement in plane-wave DFT because the exchange kernel is nonlocal; we therefore use the CRSYTAL06 code 13 , which has a Gaussian basis set, using the same atomic coordinates as were found by the VASP calculations. The basis sets used were 88-31G for Si and Bi and 31G for the terminating H atoms on the reverse of the slab, with an all-electron treatment of the Si atoms and an effective core potential 14 for the Bi atoms. The k-point sampling was a (2×1×1) mesh, giving the same density as in the VASP calculation.

A. Atomic Structure
The Si(001) surface shows characteristic 110 rows of dimerized Si atoms 15 . Dimers along the row are separated by a 0 / √ 2 = 0.38nm, where a 0 is the cubic lattice constant of Si. Tunnelling images of the Bi-implanted samples after in-situ surface preparation contained four easily distinguished features attributable to Bi donors. The bismuth-vacancy feature we describe in this work, referred to as Bi-V, is located at the surface; we have found other features which are clearly sub-surface features and will be discussed in a future publication; these sub-surface donors are considerably less common than the surface Bi-V pairs. Of our large inventory of 500 observed features, ∼84 % are Bi-V pairs and are located at the surface. Figs. 1a and 1b display empty-and filled-state topographs of a single Bi-V defect and how it fits within the dimer rows. For filled states ( Fig. 1b; a sample-bias-voltage V S of -2 V) the Bi-V feature is similar to the surrounding surface, though a little brighter. For empty states ( Fig. 1a; V S = +2 V), it appears as a bright protrusion neighbouring a depression surrounded by the characteristic dimer π -bonding states; this is a very unusual structure in the Si(001) surface. The apparent flatness of the surface surrounding the Bi-V feature indicates that the Bi is not charged-strong brightening is seen around charged dopants in Si(001) 16 , due to band bending. Considering the topographic structure, it is extremely likely that the Bi atom substitutes for a Si atom in the surface, as illustrated in Fig. 1c.
As a further check on the identity of the new features we investigated the stability of the Bi-V structure and its population density by repeated sample flashing. The density is estimated from topographic images (Fig. 2a) to be ∼ 2 × 10 11 cm −2 after repetitive surface flashing (Fig. 2b). If we assume that Bi atoms in the first two layers are visible in STM, and have the same density as produced by the initial ion implantation/annealing protocols, then we expect to see (2.2 × 10 19 cm −3 ) × (0.14 nm) ∼ 3 × 10 11 cm −2 (as Si layer spacing along (001) is 0.14 nm). This is in excellent agreement with the defect density observed by STM, strengthening our identification of these features as Bi donors. Moreover, the Bi-V structure is always found, with a population density within a factor of two of that associated with the bulk dopant density. Even after the sample was exposed to air intentionally and heated again, the same defects are seen, although with a barely significant trend towards lower density with repeated surface treatments.
To identify the structure of the feature and further characterize it, we have performed total energy spinpolarized density functional theory (DFT) calculations. As the structure clearly involves a vacancy and the prepared surface contains a certain number of vacancies, we have considered structures involving a vacancy; the energetics of structures in the absence of vacancies have also been considered and are given below. Various candidate structures were considered: a missing dimer vacancy (1DV) and a mixed Bi-Si dimer with different charge states, by analogy with the P-Si hetero-dimer system [16][17][18] ; a Bi adatom on the Si(001) surface with a 1DV; a Bi atom symmetrically replacing an entire Si dimer in the surface; and a Bi-vacancy pair replacing a Si dimer in the surface, again with different charge states. Kinetic effects during sample preparation produce an abundance of non-equilibrium structures 19 . Accordingly, the main criterion used to judge the fitness of the structure was the agreement of the appearance between simulated and real STM images, rather than total energy from DFT; the stability of different features is discussed further below. The only structure which accounts for the STM images is the neutral surface Bi-vacancy pair (Bi-V). Fig.  1c is the associated schematic, while Fig. 3a allows comparison of STM topographic images and DFT calculations. In this model, one atom of a surface Si dimer is replaced by a two-fold coordinated Bi atom, while the other Si atom is removed, and the second-layer Si atoms rebond across the vacancy as indicated by a red line in Fig. 1c. The resulting Bi atom has a neutral configuration with two 6s and three 6p electrons, unhybridized, in contrast to that in bulk Si where the donor ionizes leaving four electrons in sp 3 hybrid orbitals. This structure, and other unusual Bi structures at the Si(001) surface such as Bi nanolines 20 , can be explained by considering the large energy difference between the 6s and 6p Bi valence electrons 21 : the hybridization energy required to form sp 3 orbitals, giving the tetrahedral bonding char-acteristic of Si, is very high, and structures allowing p 3 bonding (as in native Bi bulk) will be more stable. Density of states (DOS) plots for the Bi-V structure show that the Bi 6s orbital can almost be considered as a part of the core, with only the 6p electrons forming valence bonds. Two of the three 6p orbitals form bonds with second-layer Si atoms, and the final 6p orbital remains unbonded and half-filled, aligned in the plane of the surface and perpendicular to the Si-Si dimer row direction as shown in Fig. 1c. The DOS of the Bi thus originates mainly from 6p orbitals, and the contributions from s and d orbitals are small in the energy range studied here.
Of the structures which we have modelled covering the possible atomic arrangements, only the Bi-V structure matches the bias-dependent STM images and spectroscopy. Of all structures involving vacancies, DFT finds the symmetric Bi structure to be most stable, with the Bi-V structure the next most stable; however, a hybrid DFT calculation finds that this ordering is reversed, with the Bi-V structure most stable (see below). A possible formation mechanism for the Bi-V structures at the surface is the diffusion of Bi atoms through the bulk in association with Si vacancies 24 . Further DFT calculations, also detailed below, indeed show that the Bi and Si vacancy will remain associated at the surface, once paired.
We give the total energies of various structures in Table I, from both DFT GGA and hybrid DFT calculations, as DFT will underestimate the energy of localised states thus also underestimating the stability of features such as the Bi-V. The structures modelled are illustrated in Fig. 4. The energies are all calculated for unit cells with the same number of atoms (hence we have included a surface vacancy in the calculations where Bi is not associated with a vacancy).
It is important to note various points about the energies. First, once a vacancy is present at the surface, the Bi-V or Symmetric Bi structures are most stable, as shown in Table I. Second, given the non-equilibrium preparation technique, structures which are energetically unfavourable can form because kinetic pathwaysderived from the highly non-equilibrium ion implantation process and subsequent incomplete thermal annealing steps-favour their formation. Third, while the symmetric Bi feature has a lower energy, it does not match the feature observed in STM. Fourth, calculations with hybrid functionals indicate that, for the B3LYP functional, the Bi-V structure is more stable than the symmetric Bi 38 . Finally, significant numbers of energetically unfavourable structures (e.g. defects, steps) are found in clean, well-prepared Si(001) surfaces in all experiments.

B. Electronic Structure
To understand the electronic structure of the Bi-V feature, we performed STS measurements and compared them to non-spin-polarized and spin-polarized DFT DOS calculations. Fig. 5a shows the differential conductance   Fig. 5e and f. Note that the non-spin-polarized LDOS are very different to the experiment, while the spin-polarized agree well (discussed in more detail below). Our key new result is a doublepeaked Bi-induced resonance in the Si dimer gap, visible in both the spin-polarized DFT and STM. We compare experiment with theory first for clean Si. The dI/dV spectrum taken at the Si dimer shows characteristic peaks at ∼ -0.8, +0.6 V (Fig. 5a), which are in agreement with previously reported results 4,5 . Similar features representing π (-0.4 V) and π (+0.2 V) states of the Si-Si dimer are seen in the simulated LDOS plot (Fig. 5d and f). It is well-known that DFT significantly underestimates band gaps, though it will give reasonable state densities and charge density distributions once appropriate energy shifts are made 22 . This will contribute strongly to the energy differences between experiment and theory: ∼ 0.4 V in both empty and filled states. Another factor is likely to be tip-induced band bending 23 .
The experimental dI/dV spectrum of Bi in Fig. 5a shows a sharp peak at ∼ -0.4 V and a small peak structure at ∼ +0.7 V. The -0.4 V peak can be attributed to the sharp peak at -0.15 V in the simulated spin-polarized DOS (Fig. 5e), associated with the Bi p state. Similarly, the peak at +0.7 V in the dI/dV spectrum can be matched to the +0.15 V peak in the calculated spinpolarized DOS. The extra peaks seen in the calculated DOS are caused by slight mixing with neighbouring Si atoms and are due to the method of projecting the DOS onto individual atoms. Fig. 5b summarizes the comparison of peak bias-voltages found for dI/dV and calculated DOS. Fig. 3b shows dI/dV maps around the Bi-V complex at two characteristic peak energies, -0.4 V and +0.7 V. The filled state map at -0.4 V reveals the spatial structure of the Bi donor level. The DOS around Bi becomes larger reflecting the sharp Bi p state peak, and a substantial hybridization with neighbouring Si atoms can also be observed. In the empty states at +0.7 V, the Bi site appears as a localized protrusion. To compare the experimental results with DFT calculations, including the energy shifts due to DFT errors and band bending, we have shifted the simulated DOS map in accord with the shifts of the key spectral features in Fig. 5a. Thus, the simulated DOS map at -0.15 V (corresponding to -0.4 V in dI/dV map) can account for the spatial structure of the Bi donor level observed in the dI/dV map. Also, the localized bright spot in the dI/dV map at +0.7 V is clearly reproduced at +0.15 V in the DFT calculation. The electronic structure of a Bi-Si heterodimer is shown in Appendix A, and confirms that the Bi-V feature is not a heterodimer.

C. Spin and Magnetic Structure
Fig. 5c-f present non-spin-polarized and polarized densities of states (DOS) for Bi and Si atoms. Without polarization, which is required unless the Bi is ionised, the calculated DOS of Bi is inconsistent with the measured dI/dV spectra, since a single peak at E F is predicted (Fig. 5c). Taking spin into account, the majority band associated with the unbonded Bi p-state is full and the minority band empty (Fig. 5e, left), implying full polarization, with a moment of 1 µ B /Bi, and the theory bears much closer resemblance to the data in Fig. 5a. The associated spin density is dominated by the Bi atom and one of the up Si atoms neighbouring the Si vacancy on the opposite side of the dimer row, but shows some delocalization along the row. There is no significant difference between polarized and non-spin-polarized results ( Fig. 5d and f) for the DOS at the Si sites far from the defect. Thus, the peaks seen in the dI/dV spectra for Bi-V must be due to a spin-polarized Bi state. The splitting between the occupied and unoccupied states at a single spin-polarized defect is approximately 0.2 eV, a large value supported by our experiments.
Though we are unable at present to show experimental data in support of this, we now turn to calculations of the energetics of magnetic ordering between pairs of Bi-V features. We compare the ferromagnetic and brokensymmetry antiferromagnetic total energies in order to obtain the couplings of the Bi-V pairs. Because of the p 3 bonding and the size of the atoms, the Bi atoms will always be found substituting for the 'up' position of a buckled dimer, with the vacancy at the 'down' site. We have computed only pairs where the Bi atoms are located on the same side of the dimer row; because of the 'up-down' tilt sequence found at low temperatures in the Si(001)c(4×2) structure, this implies a spacing of an even number of dimer separations (and hence an odd number of dimers separating the two defect sites).
DFT relies on a finite supercell, with two rows of eight dimers each, to obtain the DOS. It is natural to ask to what extent the splitting is derived from the ferromagnetism implied by the periodic boundary conditions. Accordingly, we have also considered the electronic structure of arrays of pairs of the Bi-V feature in the same supercell with different separations, shown in Fig. 6. For these pairs we can estimate the exchange interaction between neighbouring Bi atoms by calculating the energy difference between ferromagnetic (FM) and antiferromagnetic (AFM) arrangements of the spins localized on the Bi atoms. Since it is known that conventional DFT tends to delocalize electronic states artificially 26 , and hence over-estimate such magnetic interactions, we have also made this comparison using hybrid DFT.
The calculations considered interactions when Bi atoms are in the same row, using two rows of sixteen dimers each, with the size is constrained by the limits of DFT; all results are shown in Fig. 6. We consider Bi atoms on the same side of the row, with spacings of 4, 6 or 8 dimers, or 1.54, 2.30 and 3.07 nm. The FM/AFM energy differences derived from hybrid DFT are: -0.03eV for 1.54 nm (4 dimers); they also favour ferromagnetism but have magnitudes below 0.001 eV for 2.30 and 3.07 nm (6 and 8 dimers). This difference is significant for the 4-dimer spacing, but falls below the accuracy of our calculations for 6 and 8 dimers. The calculations indicate that the interactions are ferromagnetic, in contrast to the antiferromagnetic interaction of around +0.06 eV estimated for Bi at the same separation in bulk Si 27 , with considerable overlap of the associated Rydberg atoms. The difference is a result of the bonding of the Bi and hence its electronic structure.
There are two significant difficulties in analysing such magnetic couplings between defect pairs using DFT: first, conventional DFT tends to delocalize artificially localized spin densities, and hence over-estimate magnetic couplings 26 . Second, even in pairs on the insulating side of a Mott transition, the Kohn-Sham orbitals form delocalized molecular states over both members of the pair and the bonding state is doubly occupied in the case of anti-aligned spins (anti-ferromagnetic pair), whereas for aligned spins both bonding and anti-bonding states are singly occupied. This artificially raises the energy of AFM relative to FM structures. We solve the first problem by using the B3LYP hybrid functional 28 , which adds an empirical mixture of exact exchange into gradientcorrected DFT; this is found to give a better balance between localization and delocalization and to give magnetic couplings in better agreement with experiment in a range of structures 26,28 . We deal with the second problem by using a broken-symmetry approach 29 , in which the Kohn-Sham states are no longer required to transform according to irreducible representations of the structural symmetry group but are rotated so as to localize them on the magnetic centres and enable a fair comparison of the FM and AFM configurations.
The resulting spin densities are shown in Fig. 6. The strong two-lobed feature corresponding to the unbonded Bi p-state is clearly visible in all cases, as is the delocalization of the spin density onto the dangling bond for a neighbouring Si up-atom. In most cases this is on the opposite side of the dimer row, but for the smallest separation considered (4 dimer spacings) the strain in the surface has caused the neighbouring dimer of the righthand defect shown to reverse its tilt, and accordingly the spin density spills onto the neighbouring Si on the same side of the dimer row instead. The origin of the magnetic coupling in the delocalization of the spin density along the dimer row is clearly visible in the 4-dimer spacing case, and can also be seen in lower-density isosurfaces (not shown) for the larger separations. The localisation of the Bi state on the surface is significantly greater than in bulk, but it can couple to the delocalised π state on the surface along dimer rows; in the limit of no dimer tilt, where the surface bands are half-filled and therefore metallic, the resulting RKKY interaction will give alternating ferromagnetic and antiferromagnetic interactions along a row.

IV. CONCLUSION
We have discovered an unusual structure formed by the heavy dopant Bi at Si(001) surfaces following ionimplantation and surface cleaning. Fig. 7, where the band structures for bulk Si, Si(001) and the Bi-V feature, with and without spin polarization, are plotted, shows the essential physics. The reduction of the bandgap with creation of the surface is clear, along with the introduction of defect levels in the surface state gap. The on-site exchange, due to the Coulomb interaction, splits majority from minority spins and opens up a large gap for the Bi levels. Bi thus behaves quite unlike a bulk dopant, or even a conventional light dopant (e.g., P) at the surface. We understand the formation and nature of this structure, and its absence for lighter dopants, as the natural consequence of the much less significant bonding role of the outer s orbital for a heavy atom in column V of the periodic table. Beyond the new atomic arrangement associated with the incorporation of Bi in the surface, the electronic structure is unprecedented for a dopant at the surface of Si: the Bi atom has a half-filled, localized orbital, which carries an unpaired spin. Given the interest in ferromagnets based on semiconductors, it is especially encouraging that the DFT results show that superlattices of Bi-V defects favour ferromagnetic order, and indicate exchange couplings as large as 0.03 eV between surface Bi atoms. For the future, the calculations suggest that arrays of closely-spaced Bi-V features would be fully polarized at accessible temperatures. This would represent progress over a variety of heterogeneous thin film systems 30,31 including Si as well as dilute magnetic semiconductors such as Mn in GaAs 32,33 and Ge 34 . The half-filled orbitals which we have discovered provide an entrée to Si(100) surface-based quantum many body physics, thus complementing previous work 35 concerning adatom layers on Si(111), and could also find use in surface spin-based quantum information technology 36,37 .