First-principles study of the atomic and electronic structure of the Si(111)-(5x2-Au surface reconstruction

We present a systematic study of the atomic and electronic structure of the Si(111)-(5x2)-Au reconstruction using first-principles electronic structure calculations based on the density functional theory. We analyze the structural models proposed by Marks and Plass [Phys. Rev. Lett.75, 2172 (1995)], those proposed recently by Erwin [Phys. Rev. Lett.91, 206101 (2003)], and a completely new structure that was found during our structural optimizations. We study in detail the energetics and the structural and electronic properties of the different models. For the two most stable models, we also calculate the change in the surface energy as a function of the content of silicon adatoms for a realistic range of concentrations. Our new model is the energetically most favorable in the range of low adatom concentrations, while Erwin's"5x2"model becomes favorable for larger adatom concentrations. The crossing between the surface energies of both structures is found close to 1/2 adatoms per 5x2 unit cell, i.e. near the maximum adatom coverage observed in the experiments. Both models, the new structure and Erwin's"5x2"model, seem to provide a good description of many of the available experimental data, particularly of the angle-resolved photoemission measurements.


I. INTRODUCTION
The low-energy electronic spectrum of a one-dimensional metal is dominated by collective spin and charge excitations [1,2,3]. This is in contrast with the behavior of typical metals, that can be understood in terms of independent particle-like excitations usually called quasiparticles. These predictions are clear and well established. However, the observation of the Luttinger liquid behavior in real systems has proven quite elusive. One of the reasons for this might be that one-dimensional metals are in principle unstable with respect to the Peierls distortion that drives them into an insulating ground state [4]. A possible route to avoid this limitation is the fabrication of metallic chains absorbed on surfaces. The hope is that the rigidity of the substrate will make the energy cost for the structural distortions too large and, therefore, the one-dimensional chains would remain metallic. Semiconductor surfaces are specially attractive for this purpose. The existence of an energy gap prevents the coupling of the electronic states of the chain in the vicinity of the Fermi level with the substrate. The one-dimensional character of these states is thus preserved. It is in this context, that the fabrication of monatomic wires of metal atoms on silicon substrates has attracted much attention in recent years.
Monatomic wires of gold atoms are spontaneously formed on flat and vicinal Si(111) surfaces after the deposition of gold in the sub-monolayer regime (see Ref. 5 and references therein). Of particular interest are the vicinal substrates, where gold wires run parallel to each step-edge and the coupling between the chains, and thus the one-dimensional character of the electronic states, can in principle be controlled by changing the miscut angle. Examples of these systems that have attracted much attention in the last few years are the Si(557)-Au [6,7,8,9,10,11,12] and the Si(553)-Au [5,13,14,15] reconstructions. The photoemission spectra of these surfaces close to the Fermi energy are dominated by bands with a strong one-dimensional character that exhibit several interesting phenomena including peculiar splittings associated with the spin-orbit interaction [9] (although first interpreted as signature of spin-charge separation [6]), fractional fillings [13] and metal-insulator transitions [11]. The analogous to these one-dimensional structures in the case of the flat Si (111) surface is the so-called 5×2 reconstruction.
Already the first structural models, based on LEED measurements, considered two atomic gold chains per 5×2 unit cell running in parallel [19,20]. This was later confirmed by HREM [32] and seems to be firmly established. The gold chains run along the [110] and equivalent directions (parallel to the ×2 periodicity of the unit cell). Therefore, three different domains are possible for the 5×2 reconstruction on the flat Si(111) surface. Singledomain surfaces, necessary for ARPES, can be fabricated using vicinal surfaces with a slight cut-off angle [5,23]. The presence of one-dimensional structures in this reconstruction has also been confirmed by the ARPES studies. Early studies found a strong anisotropic signal near the Fermi level [26,31], but no evidence of Fermi-level crossing for this band was found [31]. Later studies at low temperature found a one-dimensional band with a strong dispersion along the direction of the gold chains [27,28]. The top of this band appears near the 5×2 zone boundary and disperses downward, reaching its minimum close the 5×1 zone boundary. This band has been reported to change its dimensionality from strongly one-dimensional near the Fermi energy to two-dimensional at lower energies [28]. In these studies a gap of ∼0.3 eV was also identified for this band. The presence of this gap and its apparent closing with increasing temperature was related to a Peierls instability [27,28].
More recent ARPES results [29,30], both at low and room temperatures, have been able to identify some additional surface bands. However, the metallic or semiconducting character of the surface is still a matter of debate. In fact, it has been proposed that the metallic or semiconducting character can depend on the concentration of silicon adatoms [30,33], and even that semiconducting and metallic segments can alternate along the gold chains in the surface [34].
The STM images are characterized by the presence of bright, irregular protrusions [24,25], and "Y"-shaped features [23,35] with a well defined orientation respect to the underlying substrate. The protrusions have been established to be silicon adatoms [36], which are present on the surface with an optimum coverage close to 1/4 adatoms per 5×2 unit cell.
In spite of all these experimental studies, the structure of the Si(111)-(5×2)-Au reconstruction has not been completely established yet. Earlier structural models only considered the adsorption sites of the gold atoms. Many of them could be ruled out on the bases of more detailed STM studies [24] and the knowledge of the exact gold coverage [17]. A few more refined models exist [25,32]. They consider both the position of the gold atoms on the substrate and the rebonding of the silicon atoms in the surface layer. Probably the most detailed structural model proposed to date is the one by Marks and Plass (MP) [32]. The MP model is based on a combination off-zone HREM, transmission electron diffraction and heavy-atom holography data.
The first theoretical studies using first-principles electronic structure calculations appeared only quite recently. This is due to the complicated structure and the large unit cell of the Si(111)-(5×2)-Au reconstruction. Kang and Lee [37] studied the MP and the Hasegawa-Hosaka-Hosoki (HHH) [25] models using density functional theory. Their main conclusion is that both models fail to reproduce some of the key features of the STM images and the experimental band structures. Using a similar methodology, Erwin [33] proposed and studied new structures which are characterized by the presence of the so-called honeycombchain silicon structure [38]. One of these models (the so-called "5×2" model) seems to fulfill many of the constraints imposed by the empirical evidence. An interesting point raised by In this work, we present a comprehensive study of the atomic and electronic structure of different models of the Si(111)-(5×2)-Au reconstruction using electronic structure calculations based on the density functional theory. We have used two different methodologies, the SIESTA code [39,40,41] using a basis set of localized atomic orbitals and the VASP code [42,43] using a basis set of plane-waves. We analyze the MP model [32], the models proposed by Erwin [33], and a new model that we found during our structural optimizations. We study in detail the energetics and the structural and electronic properties of the different models. We also calculate the change in the surface energy as a function of the content of silicon adatoms for the two most stable models. In order to do so, we perform calculations for large supercells containing realistic concentrations of adatoms: 5×4, 5×6, and 5×8 supercells. Our new model is the most favorable in the range of low adatom concentrations, while Erwin's "5×2" model becomes favorable for larger adatom concentrations. The crossing between the surface energy of both structures occurs close to 1/2 adatoms per 5×2 unit cell, i.e. near the maximum adatom concentration observed in the experiments. Both models, our new structure and Erwin's "5×2" model, seem to provide a good description of most of the experimental data. Particularly, we find a general agreement between the calculated and measured band structures along the direction parallel to the gold chains.

II. CALCULATION METHOD
Most of our calculations have been performed using the SIESTA [39,40] code. We have used here the local approximation (LDA) to the density-functional theory [44,45,46] (the generalized gradient approximation (GGA) has been also used for a few test calculations), Troullier-Martins pseudopotentials [47], and a basis set of numerical atomic orbitals obtained from the solution of the atomic pseudopotential at slightly excited energies [39,40,48,49].
We have used an energy shift [40,49]  points for the 5×2 unit cell (and a consistent sampling for other cells [61]) and the fineness of the real-space grid used to compute the Hartree and exchange-correlation contributions to the total energy and Hamiltonian matrix elements was equivalent to a 100 Ry plane-wave cutoff. This guarantees the convergence of the total energy, for a given basis set, within ∼20 meV/Au (∼0.5 meV/Å 2 ).
We modeled the surface using a finite slab, similar to that depicted in Fig. 1 adatom reconstruction as a test. Our results using a DZ basis and a slab containing three bilayers reproduced very well the results for the geometry and energetics given by previous calculations [51]. In fact, using a more complete DZP basis set, which includes a shell of d orbitals, or increasing the thickness of the slab (up to five bilayers) did not change appreciably the results.
A scalar-relativistic pseudopotential similar to that utilized in Ref. 52 have been used for gold. The gold basis set included double and polarized 6s orbitals (i.e. two different radial shapes to describe 6s orbital plus a 6p shell) and a single 5d shell. We refer to this basis set as DZPs-SZd. This basis set was already used in calculations for several gold clusters, where it was shown to lead to results in very good agreement with more complete basis sets [53].
We have also tested that these basis set and pseudopotential yield to the correct bulk lattice parameter and band structure.
In this work we study the relaxed structures and the energetics of several models of the Si(111)-(5×2)-Au surface reconstruction. The energy differences between different mod-els are of key importance since we would like to determine the most plausible structures.
Whenever it is necessary to compare the energies of structures containing different numbers of silicon atoms, the silicon chemical potential is set to the total energy of bulk silicon at the equilibrium lattice parameter. This choice is justified by the fact that the surface should be in equilibrium with the bulk. A summary of our results can be found in Table I. One can see that the relative energies are quite small in some cases. However, they are larger than the estimated error bar for the total energy (see above). Furthermore, the relative energies usually exhibit a faster convergence than the total energy of a single structure. It is also necessary to check the convergence of the results as a function of the slab thickness and the completeness of the basis set. Table II shows the results of these tests for the most stable structures. In one case, the slab thickness was increased by one silicon bilayer while, in the other, a DZP basis set was used for the silicon atoms. In both cases the systems were relaxed. The results are quite stable against the change of the slab thickness. In particular, the energy order of the structures is not changed and the variation of the relative surface energies is smaller than ∼0.5 meV/Å 2 in all the cases. The variations with the size of the basis set are somewhat larger. From the results in Table II we can estimate an error bar smaller than 2 meV/Å 2 for the relative surface energies of the different structures calculated using SIESTA.
In order to check the accuracy of our predictions we decided to perform calculations for some of the systems with another electronic structure code that utilizes a different methodology. We used the VASP code [42,43] for this purpose. We used projected-augmented-wave potentials and a well converged plane-wave basis set with a cutoff of 312 eV. All structures were relaxed (the equilibrium lattice parameter of bulk silicon obtained with VASP is 5.41Å). In Table II we can see some of the results obtained with VASP. They are in good agreement with the SIESTA results, especially with those obtained with the more complete DZP basis set. The order between our more stable models is preserved, although the energy difference is somewhat decreased. In particular, the new structural model found in the present work (model N in Tables I and II) is confirmed to be the most stable surface reconstruction between those studied here. It is also interesting to note that the energy associated with the addition and removal of adatoms for a particular structural model seems to be quite independent of the details of the calculation.
The surface BZs of the studied systems are shown in Fig. 2 (a). For the 5×1 system the BZ is a stretched hexagon while, for the remaining periodicities, the hexagons are distorted.
We plot the electronic band structures of the different models along the Γ-ZB ×2 -ZB ×1 -ZB ′ ×2 -M-Γ line (see the dotted line in Fig.2 (b)). The Γ-M path runs parallel to the gold wires in the surface, crossing the 5×2 BZ through three different regions. The M-Γ line is perpendicular to gold wires. The surface/bulk and main atomic character of the different bands is identified by means of a Mulliken population analysis [54].
Although a DZ basis is usually sufficient to obtain a quite good description of the occupied electronic states and the relaxed geometries in silicon systems, the use of a more complete basis set is necessary to describe the unoccupied part of the band structure even at low energies. For this reason all the band structures shown in the paper are calculated using a DZP basis set and slabs containing three underlying silicon bilayers (even if the relaxed geometry is obtained from a calculation using a DZ basis and/or a thinner slab).
Finally, the Scanning Tunneling Microscopy (STM) images are simulated using the theory of Tersoff and Hamann [55].

III. RESULTS
In this section we present our results for the different models of the Si(111)-(5×2)-Au surface. We first focus on the energetics, relaxed geometries, and the electronic band structures. We then turn our attention to the effect of the different silicon adatom contents and the simulated STM images, which we only analyze in detail for the most stable structural models. A summary of the relative energies of the calculated configurations, accompanied with a brief description of each of them, can be found in Table I. Before starting with the description of the results, it is interesting to point out some brief comments about the concentration of silicon adatoms on the Si(111)-(5×2)-Au surface. A detailed study of the equilibrium situation has recently been performed by Kirakosian et al. [36,56] using STM. Their results indicate that, at equilibrium, only one adatom site is occupied out of every four possible sites, corresponding to a 5×8 adatom periodicity (if all the adatom sites were occupied we would recover a perfect 5×2 periodicity). The analysis of the adatom-adatom correlation functions obtained from the STM images reveals a strong suppression of those configurations with small adatom-adatom distances, a clear maximum corresponding 5×4 periodicities, and a long range oscillatory tail [56]. This was interpreted in terms of a short range repulsion between adatoms plus a long range interaction term. In Ref. 36, Kirakosian and collaborators showed that the density of adatoms can be increased by depositing additional amounts of silicon reaching an almost perfect 5×4 arrangement of the silicon adatoms. Further deposition of silicon does not create a stable 5×2 adatom structure. Instead the extra silicon atoms decorate the step-edge of the terraces on the surface. These observations seem to have at least two implications: (i) the optimal adatom concentration must be certainly lower than one adatom per 5×2 cell and, (ii) the structure of the reconstruction must be stable against relatively large changes of the content of adatoms [62] since the density of silicon adatoms can be increased by a factor of two without, at least apparently, dramatic structural changes [36].
A systematic study of the energetics of the surface as a function of the adatom concentration by means of first-principles electronic structure calculations is quite complicated. This is for two main reasons. First, the energies involved are rather small, which implies the need of very well converged calculations. A more serious limitation, however, is the necessity to use large supercells consistent with the low adatom densities. For this reason we have concentrated most efforts in the two limiting cases, involving respectively 0 and 1 adatoms per 5×2 cell. The intermediate concentrations usually require drastic approximations. For example, Erwin [33] assumed that the main effect of the adatoms in the Si(111)-(5×2)-Au surface is to dope the gold chains with electrons and studied the energetics of the system as a function of this doping. Here we go a step beyond and present explicit calculations for adatom contents down to 1/4 adatoms per 5×2 cell, consistent with a 5×8 periodicity, which indeed can be reached in experimental conditions [36]. Due to the large size of these systems we limit this study to our two most stable models, and only use the smaller DZ basis set.

A. Marks and Plass model
We start our investigation of the structure of the Si(111)-(5×2)-Au surface using the model proposed by Marks and Plass [32] from experimental data obtained with heavy-atom holography and high resolution electron microscopy. We use the label MP + for this structure (see Table I). The + superscript indicates that the structure contains silicon adatoms saturating some of the silicon dangling-bonds in the structure, what we call "conventional" silicon adatoms. A schematic view of this structure, as proposed in Ref. 32, can be found in Fig. 1. It has to be taken into account that, due to the limitations of the experimental techniques, there are several uncertainties and assumptions in this structure. Only the atomic coordinates within the surface plane are accurate. The heights of the atoms over the substrate are only approximately resolved. The experimental beam error in combination with the size and complexity of the structure also limits the sensitivity to possible subsurface relaxations. As a consequence, the experimentally proposed structure only considers the reconstruction of the outermost bilayer and contains limited information about the registry between this surface bilayer and the underlying material. It is necessary to eliminate these uncertainties before one can undertake any serious study of the electronic structure of the MP + model. In order to do this while preserving all the information originally present in the MP + proposal, we started our study by performing constrained relaxations of the structure. The structure in Fig only implies approximately equal z-coordinates ( the z-axis is taken here along the surface normal). For this reason, in a second step, the atoms were allowed to relax in the z-direction while keeping their coordinates within the xy-plane. The resulting geometry preserves the bonding pattern of the original MP + proposal, and provides a reasonable initial guess to start our search of the most stable models by performing full structural optimizations.
We now consider in detail some of the structural patterns appearing in the MP + model of the surface. For this analysis we find useful the comparison with the Si(557)-Au surface, a closely related reconstruction that has been quite well characterized during recent years [6,7,8,9,10,12]. The stepped Si(557)-Au is formed after the deposition of ∼0.2 monolayers of gold on vicinal (111), with the misorientation chosen along the [112] direction. The size and orientation of the terraces of the Si(557)-Au represent an analogous to the flat 5×2 unit cell but including a single silicon step [27]. With half the gold coverage than the Si(111)-(5×2)-Au surface, the terraces of the Si(557)-Au contain only one Au wire running parallel to the step edges. Gold atoms occupy silicon substitutional positions in the surface layer. This is supported both, by recent X-ray diffraction data [12], and density functional calculations performed using a methodology similar to the one utilized here [8], which provide a consistent structural model of the surface. In particular, the highest stability of the silicon substitutional sites for gold has been unambiguously demonstrated by the ab initio calculations. For example, the substitutional site was determined to be at least 1 eV/Au more stable that adatom-like positions, where gold sits on the surface saturating one of the silicon dangling-bonds, or even ∼0.5 eV/Au more favorable than the adsorption decorating the step edges [8]. It seems, therefore, that the Au atoms on the Si(557)-Au surface exhibit a strong tendency towards three-fold silicon coordination. Gold atoms adapt to this situation without much strain, with typical Si-Au distances only a few percents larger than the bulk silicon bond length.
In the light of these observations the bonding pattern of some of the gold atoms in the MP + model (Fig. 1) seems quite peculiar. In particular, the gold atoms in the chain situated at the left side of the "gold trench" (marked with an L in Fig. 1) present a fourfold coordination. They are connected to three silicon atoms within the surface layer and, additionally, to the silicon atom immediately below. Furthermore, the Si atoms neighboring to the mentioned gold atoms (see atoms a and a ′ in Fig. 1) present an unsaturated dangling bond which might be avoided with a slight structural change.
It is interesting to note that the tendency of the gold atoms to occupy silicon substitutional positions in the top most layer cannot help to completely rationalize the structure.
A three-fold bonding pattern of the gold atoms is inherently frustrated by the presence of a surface dislocation. In the MP + model this dislocation is located at the position of the right-hand gold wire (marked with R in Fig. 1). Due to the change of the bonding sequence there are not three unpaired silicon electrons available for each of these Au atoms, but rather two. Therefore, they do not occupy a normal three-fold position and are quite likely to be displaced from the initial symmetric positions after relaxation. It should be noted that, in principle, the surface dislocation can be moved to different locations. In fact, we will see below that this provides a simple route to generate alternative structural models of the surface.
The comparison between the structure of the Si(557)-Au reconstruction [8,12] and the MP + model of Fig.1 raises another interesting point. In the case of the Si(557)-Au surface the silicon atoms in the proximity of the step edge suffer a considerable rebonding. They form characteristic silicon structure which has been identified with the so-called "honeycomb chain" (HC) by several authors [5,8]. The presence of the silicon HC seems instrumental to understand the stability of the Si(557)-Au and related reconstructions (see, for example, Ref. 5). The silicon HC was initially proposed by Erwin and Weitering [38] as the main building-block of the (3×1) reconstruction induced on Si(111) by the deposition of metals like Ag, Li, Na, K, Mg or Ba. The HC structure represents a large disturbance from the usual bonding pattern of silicon. The stability of the silicon HC stems fundamentally from the formation of a double-bond between two three-fold coordinated silicon atoms on the surface [38]. In the case of the (3×1) reconstruction a further stabilization mechanism comes from the fact that, after the transfer of the valence electrons from the metal atoms, the HC structure becomes electronically closed-shell. [63] It seems somewhat surprising that the silicon HC structure, common to the (3×1) and Si (557) We now proceed further with the structural relaxations of the MP + system. It is instructive to focus first on a optimization were only the silicon degrees of freedom are taken into account. The gold atoms are constrained to remain at their initial coordinates. Due to the more directional bonding of silicon we can expect the structural changes to be simpler to analyze and somewhat less dependent on the particular choice of the initial guess in this case. Besides, as a stronger scatterer, we can assume that the gold positions to be better resolved in the experiment. The resulting geometry is plotted in Fig. 3(a). We observe two main effects. On the one hand, the HC configuration clearly emerges. One of the driving forces behind this rebonding is the movement forward of a and a ′ atoms in order to form an additional covalent bond with the silicon atoms in the underlying layer. The doublebonded "dimers" of the HC structure are formed by atoms b and c. This questions the location of the adatoms in the surface since, in principle, the dangling-bond associated with atom b ′ could disappear with the formation of a silicon double bond. We can observe the disturbing effect of the adatom on the HC structure. The appearance of the HC bonding pattern during the relaxation of the MP + structure confirms the results of recent density functional calculations by Kang and Lee [37], who also made a geometrical optimization of the MP + model. The electronic bands calculated for this structure (not shown here) are also in quite good agreement with those presented by these authors in Ref. 37. Fig.3 (a) also shows clearly what could be classified as a "stacking fault" in the structure (bonds of atoms d and e coincide with those in the underlying silicon layer). This stacking fault, which probably is energetically unfavorable, can be easily avoided by moving the position of the surface dislocation from the right-hand to the left-hand of the gold trench. Alternatively we can visualize this change (at least approximately) as a 180 o rotation of the surface layer with respect to the underlying silicon structure. This transformation gives one of the structures discussed in the next section, which incidentally is almost identical to the "5×1" structure proposed recently by Erwin in Ref. 33.
When the relaxation of the MP + system is continued without any constraints, the monatomic gold wires are strongly distorted as can be seen in Fig. 3 (b). This distortion was not observed in the density functional calculations of Kang and Lee [37]. The reasons for this discrepancy are not completely clear at the moment. The break of the monatomic gold wires seems to be related with the presence of adatoms. If they are eliminated from the structure the gold atoms remain in two well separated parallel wires. Additionally, the strain introduced by the adatoms in the structure, results in the weakening of some Si-Si bonds in the surface layer (see the increased distance between atoms f and g). In spite of these strong structural distortions, the presence of adatoms in the structure is still energetically favorable as can be seen in Table I In particular, it can be translated to the other gold wire, this can also be assimilated to a rotation of the surface bilayer with respect to the underlying bulk silicon. This eliminates the "stacking fault" commented in the previous section, and produces a new structural model. This structure is very similar to the "5×1" model recently proposed by Erwin [33], and we refer to it as E(5×1). In Table I we can see that the E(5×1) model without silicon adatoms is slightly more stable than the relaxed MP structure. Fig. 4 (a) shows the relaxed structure of the E(5×1) model. The left (L) gold wire, where the surface dislocation is located, suffers a considerable dimerization, which is much smaller for the right (R) wire. The alternating Au-Au distances as obtained with VASP are, 4.06Å and 3.59Å for the L wire, and 3.82Å and 3.83Å for the R wire. The geometries obtained with SIESTA are very similar, specially those obtained with the more complete DZP basis set. Hereafter we name "SiAu complex" the structure formed by the two gold wires and the central silicon atom connecting them.
The silicon structure in between two of such SiAu complexes is quite flat and resembles what could be described as a double honeycomb chain (DHC) silicon structure [33]. The band structure along the direction parallel to the gold wires is shown in Fig. 4 (b). It shows several surface bands and has a metallic character. Those surface bands mainly associated with the Si-Au complex has been highlighted using solid symbols. Most of these bands are occupied and appear in the gap region. The unoccupied surface bands appearing in the gap are mainly associated with the silicon DHC. The most prominent feature is a dispersive band associated with the weakly dimerized (right) gold wire and the central silicon atom in the SiAu complex. This band is, in principle, metallic and close to half occupied. Although small gaps are opened associated with the crossings with other bands and slight geometrical distortions, it can be easily followed in Fig. 4 (b) extending from ∼1.3 eV below to ∼2.3 eV above E F . A similar band, with a similar origin, also dominates the band structure of the Si(557)-Au surface [7,8]. This band comes mainly from the sp 3 lobes of the central silicon atom in the SiAu complex. There is also a strong hybridization with the 6p states of the gold atoms in the R wire. For this reason, they are better assigned to the Si-Au bonds connecting the central silicon with the R gold wire. Its large dispersion is due to the large overlap between these Si-Au bonds along the wire. The metallicity stems from the inability of gold (each gold atom only provides one valence electron) to saturate the bonds with all its silicon neighbors [7]. The other states in the Si-Au complex give rise to relatively flat surface bands associated either with weakly overlapping silicon states or with the gold dimers.
In Ref. 33 it was also proposed that, under certain conditions, it could be energetically favorable to remove some of the over coordinated silicon atoms in the neighborhood of the surface dislocation. Our relaxed structure for this model (hereafter E(5×2)) is shown in Fig. 5. In this case both gold wires present an appreciable dimerization with alternating Au-Au distances of 4.37Å and 3.35Å for the left gold wire, and 4.16Å and 3.49Å for the right wire. Our SIESTA calculations with the smaller DZ basis set predict the E(5×2) model to be more stable, by at least 3.4 meV/Å 2 , than both the E(5×1) model and the different variants of the MP model (see Table I). However, the difference between the E(5×1) and E(5×2) models is reduced with the use of more complete basis set. In particular, our planewave calculations predict both models to be degenerate within 0.1 meV/Å 2 (the E(5×1) slightly more stable). This agrees with the results of Ref. 33 where the E(5×1) model is predicted to be more stable than the E(5×2) variant by less than 1 meV/Å 2 , and only after the addition of silicon adatoms the E(5×2) structure becomes favorable.
The band structure of E(5×2) with zero adatom coverage is plotted in Fig. 6. The band structure along the wires is in good agreement with that reported in Ref. 33 for this structure. Again, the surface bands close to the Fermi energy come mainly from the SiAu complex. Like in the case of the E(5×1) model, the band structure is metallic. This is in disagreement with one of the latest and more detailed ARPES experiments which suggests that the Si(111)-Au-(5×2) surface is a semiconductor with a band gap of at least 0.2 eV [29].
However, the metallic versus semiconducting character of this surface is still a matter of controversy. For example, the recent ARPES study by Himpsel and collaborators finds several metallic bands [30]. In fact, this reference and the scanning tunneling spectroscopy (STS) data of Ref. 34 indicate that the surface could be composed of alternate metallic and semiconducting regions along the gold wires. Our calculated band structure for the E(5×2) model is very close to being semiconducting. Just by shifting the S 1 band to higher energies by a few tenths of eV we could obtain a semiconducting surface. This might indicate that the metallic behavior is simply related to the limitations inherent to the local density approximation used here and the very simplified assumption that the monoelectronic eigenvalues can be directly identified with the photoemission peaks. In spite of its metallicity, several characteristics of the photoemission spectra are recovered by the band structure in Fig. 6. The most prominent band observed experimentally starts at the boundary of the 5×2 zone (ZB ×2 ) dispersing downwards until it reaches a minimum at the boundary of the 5×1 zone (ZB ×1 ) [27,28,29,30]. This band appears at binding energies between ∼0.2 eV and ∼1.3 eV. Following Erwin [33], we can try to identify this band with our S 2 band, whose maximum appears close to E F in the neighborhood of ZB ×2 . However, it becomes difficult to follow the dispersion of this surface band as we move to higher binding energies for two reasons: i) the band enters the region of the projected bulk bands, becoming a surface resonance and, ii) other surface bands coming from the same region of the surface appear in the energy interval between -0.5 and -1.2 eV. This last point is widely consistent with the experimental data in Ref. 29, where three additional bands are identified for binding energies larger than 0.5 eV.
Losio and collaborators [28] reported an interesting effect, a continuous dimensionality transition of the main surface band. The character changes from strongly one-dimensional at the band maximum (i.e. only dispersing in the direction parallel to the gold wires) to twodimensional at its minimum (i.e. with a non-negligible dispersion also in the perpendicular direction). The strong one-dimensional character of the surface states close to E F has also been confirmed in the most recent ARPES measurements [29,30]. The dispersions in the direction perpendicular to the wires can be found in Fig. 6 (b) and (c). The band widths are rather small for most surface bands. An effect similar to the reported dimensionality transition can be seen in the case of the S 1 band. It is tempting to assign the experimentally observed effect to the S 2 band (see the different dispersion of bands S 2 in panel (b) and S ′ in panel (c)). However, as commented above it is not so simple to follow the S 2 band as it disperses downwards. In fact, we can locate what seems to be an avoided crossing between the S 2 and the S 3 bands half way along the ZB ×2 -ZB ×1 path in Fig. 6 (a). Therefore, we think that the S ′ band in panel (c) is rather related to the S 3 band than to the S 2 band, and the dimensionality change would be absent from our results. Also the energy position of the band S ′ (∼-0.5 eV) is quite far from the ∼-1.3 eV found experimentally for the band minimum. Therefore, in contrast to Erwin [33] we conclude that our calculated band structure for the E(5×2) model does not provide a direct explanation to the observation by We now explore the role of the silicon adatoms in these structures. We first studied the stability of the adatoms in the E(5×1) model when they are located over the silicon part of the surface reconstruction, i.e on sites equivalent to those occupied by the adatoms in the original MP proposal. It is interesting to note that the role of the silicon adatoms in such positions is indeed not very clear. The stability of the adatoms in typical silicon reconstructions stems from the fact that each adatom can saturate three dangling bonds in the surface at the expense of creating just an additional dangling bond. The energy gained in this process usually overcomes the strain energy caused by the addition of the adatoms. However, the E(5×1) model in Fig. 4 (a) does not have silicon dangling bonds.
The appearance of unsaturated dangling bonds is avoided by the formation of the doublebonded silicon dimers that characterize the HC configuration. In fact, the only metallic band in this model comes from the SiAu complex as explained above. In accordance with these observations, we found extremely difficult to reach a stable configuration, i.e. with all the components of the forces below our threshold, for the silicon adatoms over the silicon sites of the E(5×1) model. Finally, after several hundreds of optimization steps this model spontaneously relaxed into a new structure. This structure, labeled N + in Table. I, belongs to a new family of structural models for the Si(111)-Au-(5×2) surface found in this work for the first time and described in more detailed in the next subsection.
In the light of the previous comments, a more stable adsorption site for the silicon adatoms would be on top of the SiAu complex. This has been previously proposed by Erwin [33], and is confirmed by our calculations. Table I shows the changes in surface energy after the addition of one silicon adatom per 5×2 unit cell. The behavior is opposite for the E(5×1) and E(5×2) models, with the addition being energetically favorable for the later model.
The E(5×1) remains metallic after the addition of the adatom, and the dispersive band associated with the SiAu complex remains quite unchanged. The situation with the E(5×2) model is different. In agreement with the results in Ref. 33 we find that the band structure becomes semiconducting after the addition of the adatoms. The corresponding atomic and electronic structure can be found in Fig. 7 (a) and (b) respectively. The surfaces bands with a larger contribution from the atoms in the SiAu complex has been highlighted using solid symbols. It has been impossible to identify a band that can be solely assigned to the adatoms.
We can see that the band structure of the E(5×2) suffers major modifications after the addition of adatoms, at least for the large concentrations considered here. Besides the fact that the structure becomes semiconducting, the agreement with the detailed photoemission experiments of references 29 and 30 seems to be somewhat degraded.

C. New structural model
In this section we present a novel structural model for the Si(111)-(5×2)-Au surface reconstruction that has been found during our investigation. Our slab spontaneously relaxed to this new structure while trying to optimize a modified version of the E(5×1) model commented in the previous section. The new structure can be found in Fig. 8, and will be referred here as model N. Table I shows that the energy of the new model compares favorably with those of the other structures proposed to date. In fact, within our calculational scheme it is the most favorable structure. The difference with the second most stable model without adatoms, the E(5×2), is 4.7 meV/Å 2 . This difference is reduced to 4.1 meV/Å 2 when using a more complete DZP basis set as shown in Table II. These energy differences are quite small, so further studies have been performed in order to drive more definitive conclusions. First, we have repeated our calculations using the PBE [57] GGA exchange and correlation functional instead of LDA. The new structure continues to be more stable by 3.6 and 5.1 meV/Å 2 using, respectively, a DZ and a DZP basis set. As a second step, the energy ordering between the N and the E(5×2) structures has been confirmed using VASP and slabs containing three and four silicon double-layers. The new model is more stable than the E(5×2) by at least 2.6 meV/Å 2 . These results convincingly establish, at least within the framework of density functional calculations, the larger stability of our new structural model compared to previous proposals in the limit of negligible adatom coverage Given the small energy differences between both models it may be interesting to estimate the effect of the vibrational degrees of freedom in the surface free energy γ(T ). The vibrational contribution can affect the energy ordering even at zero temperature due to the zero-point energy, and its importance grows with temperature. Unfortunately, an accurate estimation of the vibrational surface free energy γ vib (T ) is a formidable task that would require the detailed calculation of the dynamical properties (phonon band structure) of the different surface models. This is a computationally very demanding calculation that is beyond the scope of the present paper. We can obtain a rough estimation of the vibrational contribution to the difference of the surface free energies between the different structures ∆γ(T ) following Ref. 58. We have ∆γ(T ) = ∆E surf + ∆γ vib (T ), Here ∆E surf is given in Table I  deduced from the HREM studies of the surface. [32] Similarly to the E(5×1) structure, most of the surface of the N model is covered with a silicon double honeycomb chain structure [33]. One of the silicon atoms in the DHC appears at a higher position over the surface. This indicates that this atom has a charged danglingbond and, therefore, is trying to develop a sp 3 hybridization. This atom is expected to be more visible in the STM images and to provide a preferential site for adsorption on the surface, in particular for possible silicon adatoms. The boundaries between the DHC stripes are occupied by the SiAu complex, in which a central silicon atom appears bonded with three gold dimers.
The band structure of the new structure is plotted in Fig. 9. The general features are in good agreement with the most recent ARPES studies [29,30], although some of the details are different. The most dispersive and prominent surface bands are quite similar to those found for the E(5×2) model. The surface is predicted to be semiconducting, which agrees with the results of Ref. 29. The bands named S 1 and S 2 by Matsuda et al. [29] can be easily identified in our calculation, and we use the same notation. Other less dispersive surface bands are also observed in our calculated band structure. These can be tentatively identified with those labeled S 3 and S 4 by Matsuda et al.. However, the S 3 band appears shifted to lower binding energies by a few tenths of eV. We can relate this upward energy shift to the use of the LDA in our calculations, which is likely to be less suited to describe more localized (less dispersive) states. Besides this energy shift, the sole major discrepancy with the experimental band structure in Ref. 29 is the absence of the S ′ 3 band. However, this band is not so clearly resolved in the experiments as the others.
Different symbols are used in Fig. 9 according to the main atomic character of the bands. We find several unoccupied surface bands whose atomic character is difficult to determine.
One of these bands is located at energies very close to E F , particularly near the Γ point.
The metallic/semiconducting character of the surface is thus governed by the position of this band. This situation is very similar to that already observed for the E(5×2) model, although in this case the band reaches to lower energies and becomes partially occupied driving the system to metallic.
In agreement with experiment, most surface bands show a strong 1D character in our new structural model as can be seen in Fig. 9 (b) and (c). This is particularly clear in panel (b), where most states are located within the bulk gap. In the region displayed in panel (c) (at the zone boundary of the 5×1 Brillouin zone) the S 2 and S 1 bands merge with the bulk bands, becoming surface resonances. It is no longer possible to identify the S 1 and S 2 resonances with a single band of our finite slab and, as a consequence, it is difficult to follow the band dispersion of these spectral features in the direction perpendicular to the gold wires. However, from the data in panel (c) it is clear that the combined effect of the possible dispersion, plus the broadening of the resonances extends over a range of ∼0.2 eV, much larger than its dispersion for energies closer to E F . This is broadly consistent with the 1D to 2D transition reported in Ref. 28 for the most prominent photoemission feature as the binding energy increases.
We now explore the structure and energetics of the model N under the addition of one silicon adatom per 5×2 unit cell. We tried several different adsorption sites: directly on the SiAu wire following the proposal by Erwin [33] (referred as N ⋆ ), and bonded to the prominent dangling bond in the DHC structure occupying hollow H 3 (N + ) or top T 4 (N +′ ) sites [59].
As shown is Table I, this high coverage of adatoms is energetically unfavorable in all cases by at least 2.7 meV/Å 2 . This is in contrast with the situation for the E(5×2) model, where the addition of one silicon adatom per unit cell is slightly favorable. In the N ⋆ structure (not shown) the silicon adatoms tend to locate in a peculiar bridge position between two gold dimers along the [110] direction. The structure of the N + model is shown in Fig. 10 (a). The silicon atoms bonded to the adatom adopt a typical silicon configuration although, contrary to what is observed for the clean Si(111) surface, the hollow site is preferred over the top site [51].
The band structure of the N + surface is shown in Fig. 10 (b). It is very similar to that found for the model without adatoms. The S 1 and S 2 are largely unchanged, which clearly indicates its origin in the SiAu complex. The flat S 3 band disappears from the gap region as a consequence of the saturation of the dangling bond with the adatom. A new unoccupied band, associated with the adatoms, appears instead. This new band can be found around ∼0.6 eV above E F in Fig. 10 (b).

D. Adatom coverage
So far we have only considered the limiting cases with zero or maximum adatom coverage, which correspond to a number x of silicon adatoms per 5×2 unit cell equal, respectively, to 0 and 1. However, the experimental evidence indicates that the equilibrium concentration is x ∼ 1/4, corresponding to one adatom per 5×8 supercell. Under silicon rich conditions the adatom coverage can be increased in the experiment only up to x ∼ 1/2, consistent with a 5×4 periodicity. We have performed explicit calculations for x = 1/2, x = 1/3, and x = 1/4 for our two most stable models of the reconstruction in order to simulate these situations that can be reached experimentally. Due to the very large supercells necessary for these calculations (up to 273 atoms), we have performed them with the SIESTA code and restricted to the use of a DZ basis set for silicon. The results of the energetics as a function of the adatom content can be found in Table III and in Fig. 11. The behavior is opposite for both models, N and E(5×2). It should be kept in mind that model N favors the adatoms in hollow sites over the silicon surface, while in the E(5×2) structure the adatoms sit on the gold chains are more favorable.
The surface energy monotonously decreases as a function of the number of adatoms for the E(5×2) model. We do not find any evidence of an energy minimum as a function of the adatom concentration. This is in contrast with the suggestion made by Erwin in Ref. 33. In that reference the addition of adatoms was studied using the following simplification: it was assumed that the sole effect of the adatoms is to dope the gold chains with electrons and the energy of the system was studied as a function of the doping. Erwin found a minimum of the total energy for 0.5 extra electrons per 5×2 unit cell. Since each adatom was found to donate two electrons to the surface, this would correspond to the observed adatom concentration at equilibrium of x ∼ 1/4. However, our simulations introducing explicitly the adatoms in the structure do not confirm this behavior. The surface energy of the E(5×2) structure always decreases as the adatom concentration is increased. However, the slope of the curve becomes very small for intermediate adatom concentrations, showing a weak dependence of the surface energy in that region. We cannot completely rule out the presence of a minimum for the surface energy at very low adatom concentrations. However, it seems quite unlikely looking at Fig 11. It could also be argued that the DZ basis set is not flexible enough to produce the correct behavior. This seems quite improbable looking at the data in Table II, which clearly show that the energy changes induced by the addition of adatoms are weakly dependent on the details of the calculation.
In the case of the new model N, the Fig. 11 shows that the surface energy systematically increases as a function of the adatom concentration. With the DZ basis set the N model is always more stable than the E(5×2) structure. Using a more complete basis set and a converged plane-wave calculation we find a crossing: the new model is always more stable at low adatom coverage, but becomes unstable compare with the Erwin "5×2" model at larger coverages. Scaling the data calculated with the DZ basis set to reproduce the VASP results at the end points (i.e. x = 0 and x = 1) we can estimate that the crossing occurs at x ∼1/2.
We can conclude then that the N model is, at least in the framework of density functional calculations, more favorable than the E(5×2) model for adatom concentrations below ∼1/2 adatoms per 5×2 cell.

E. Simulated STM images
The STM images of the Si(111)-(5×2)-Au surface are characterized by the presence of bright "protrusions" and "Y"-shaped features with a definite orientation respect to the underlying lattice [23,34,35]. It seems quite well established that the protrusion correspond to silicon adatoms [36,56,60]. However, the origin of the "Y"-shaped features is less clear.

IV. CONCLUSIONS
We have performed a systematic study of different models of the Si(111)-(5×2)-Au surface reconstruction by means of first-principles density-functional calculations using the SIESTA [39,40] and the VASP [42,43] codes. We start our investigation with the structural model proposed by Marks and Plass [32] (MP). This is the most detailed model of this surface reconstruction solely based on experimental information to date. Therefore, it provides a logical starting point for our study. We have also considered different variants of the relaxed MP model, including the structures recently proposed by Erwin [33], and a new structure found during our simulations. Within the computational schemes used here this new structure is the most favorable energetically, at least in the regime of low concentration of silicon adatoms. In general, we find a reasonable agreement between our results and those of the two existing theoretical studies of the surface [33,37]. The energy differences between different models are quite small, with most structures lying in a narrow range of surface energies of less than 10 meV/Å 2 (the estimated error bar for our energies is of the order of 1-2 meV/Å 2 ). This, together with the uncertainties arising from the use of the local approximation to the density functional theory, make difficult to draw definitive conclusions solely based on the energetics. The comparison of the calculated band structures and local density of states, respectively, with the available ARPES data [28,29,30] and the STM images [23,24,34] becomes then instrumental in order to identify the most plausible candidates for the equilibrium structure. In the following we summarized some our main conclusions: i) Like in the case of the reconstructions formed by the deposition of gold on stepped silicon surfaces [5,7,8,12], the silicon honeycomb chain (HC) [38] structure emerges as a fundamental building block of the reconstruction. In agreement with the result of Hang and Lee [37], the silicon HC is formed spontaneously during the relaxation of the MP model. The HC is also present in the optimized geometries of all the other structural models considered in our work.  SiAu complex. We arrive in this way to a new structure, the N model. According to our calculations this new structure is more stable, at least for low coverages of silicon adatoms, than any of the models proposed to date. The distance between the gold wires in this model is ∼3Å, which seems somewhat small compared to the ∼3.9Å deduced from the HREM measurements [32].
v) The calculated band structures of the E(5×2) and N models without adatoms are quite similar and appear to be in reasonable agreement with the available ARPES data [28,29,30]. The other models fail to reproduce the main features observed experimentally.
The agreement seems to be particularly good in the case of the N model. According to our analysis the most prominent and dispersive surface bands, named S 1 and S 2 in Ref. 29, come from the atoms in the SiAu complex. In the case of the N model the silicon adatoms tend to adsorb on the silicon part of the surface, i.e. bonded to three silicon atoms in the surface layer. As a consequence, the topology and the energy position of these bands are quite insensitive to the coverage of silicon adatoms. This contrast with the situation found for the E(5×2) model. Here the silicon adatoms tend to adsorb directly on the SiAu complex, thus causing a notable modification of the surface bands that worsens the agreement with the experimental ARPES spectra. unfavorable. As a consequence of this opposite behavior, the E(5×2) structure becomes more stable than the N structure in the limit of relatively large adatom concentrations, for x > ∼ 1/2 (notice that x = 1/2 corresponds to a 5×4 periodicity). According to this picture the exact content of adatoms is instrumental to determine the equilibrium structure of the reconstruction within the range of experimentally realizable adatom coverages. This introduces a new degree of complexity that should be taken into account when analyzing the experimental information. In particular, this might be behind the observed phase separation into 5×4 and 5×2 patches [30,56].
vii) The simulated STM images of the most stable models, N and E(5×2), are in broad agreement with the experimental images. The silicon atoms produce bright spots which are located in the middle of the underlying row structures for the E(5×2) and in a somewhat more lateral position for the N model. In both cases "Y"-shaped features similar to those observed in the experiment can be found. However, they are more clear in the case of the E(5×2) model [33] where the structure surrounding the gold chains is less symmetric.   [59]. The presence of adatoms located on top of the Au wires is indicated by a ⋆ superscript. The data in this table have been calculated using the SIESTA code with a DZ basis for silicon and DZPs-SZd basis for gold. The slabs contained two silicon bilayers below the surface layer (see Fig. 1 (a)). All energies are referred to that of the structure recently proposed by Erwin in Ref. 33.     Table I).  The results obtained with the DZP basis set (diamonds) and with plane-wave VASP calculations (triangles) for the two limiting cases are also shown for comparison. All energies are referred to those of the E(5×2) ⋆ model.