Introduction

Recent advances in scanning tunneling microscopy have enabled extensive control of single atoms placed on different surfaces1,2,3,4,5,6,7,8. These techniques paved the way for the creation and investigation of a new class of synthetic two-dimensional materials constituted by atomic structures allocated on top of different substrates, e.g., systems of tin (Sn) and lead (Pb) adatoms disposed periodically on silicon Si(111)9,10,11,12,13,14,15, germanium Ge(111)16,17,18, or SiC(0001)19 surfaces. The possibility of tuning the structure and the chemical composition in these two-dimensional systems allows for a direct modification of properties in a similar way to cold atoms or Moiré heterostructures, as proposed in ref. 20. As a consequence, they can be seen as a promising platform for simulating various quantum effects21,22,23.

At the band structure level, depositing a monolayer of group-IV atoms on a Si(111), Ge(111), or SiC(0001) surface leads to the formation of a half-filled narrow band that is well-separated from the bands of an insulating background. On the one hand, this situation could allow for an application of the most advanced theoretical many-body approaches developed to date for model single-band systems (see, e.g., refs. 24,25). On the other hand, these materials exhibit a number of non-trivial features that make the solution of the problem not straightforward. For instance, the wave function of single-particle states is very extended, which results in a strong and long-ranged Coulomb interaction26,27 that has to be taken into account. Another important aspect that has to be considered is the strong spin-orbit coupling (SOC) that emerges in the case of heavy adsorbents (Sn, Pb, etc.)28.

The interplay between collective excitations and structural effects has been extensively investigated using a combination of experimental and theoretical methods in two-dimensional materials29,30,31,32. However, until recently, the theoretical investigation of these surface nanostructures was mostly dedicated to the description of the metal-insulator transitions observed in scanning tunneling spectroscopy and photoemission spectroscopy experiments14,27,33. Much less attention has been paid to collective electronic effects and, in particular, to magnetic properties, and the obtained results were controversial34. First-principles simulations using density functional theory predicted an antiferromagnetic ground state for the Si(111):Sn material35. On the other hand, it has been shown that taking into account more distant hopping processes instead stabilizes a row-wise collinear order36. It should also be noted, that both these calculations were performed without considering the effect of SOC, which may substantially affect the magnetic state. Unfortunately, there is still no direct experimental confirmation of which magnetic ordering is actually realized in the material.

Theoretically, the phase diagram of Pb adatoms deposited on a Si(111) surface is one of the most poorly understood features in this class of compounds. Similarly to the Si(111):Sn material, the Bravais lattice of the Pb adatom system is rotated by 30∘ with respect to the substrate. A known peculiarity of the triangular lattice is a high degree of frustration that can lead to a non-trivial competition between different ordering phenomena. In addition, Si(111):Pb displays very strong on-site and spatial electron-electron interactions27, which makes the system an ideal candidate to study charge and spin fluctuations. Finally, one would also expect noticeable effects related to the SOC, since the Pb adatoms have a sufficiently large atomic number28. In particular, the SOC results in a splitting of the Fermi surface, which can be observed experimentally in the quasiparticle interference pattern. Additionally, the SOC gives rise to the magnetic Dzyaloshinskii–Moriya interaction, which, in turn, can lead to the formation of chiral spin textures with non-commensurate ordering vectors28. Previous calculations on the phases of this system made use of methods unable to properly account for the short and long-range correlations appearing in this system, such as DFT15, Hartree–Fock28 and cluster methods37, or approaches that do not include magnetic fluctuations26. Notably, DFT predicts a metallic behavior15, whereas more correlated methods converge to a Mott-insulating behavior26,28.

Experimentally, it was observed that the Si(111):Pb system indeed shows a non-trivial behavior related to the above-mentioned features. Several different arrangements of the atoms on the surface were identified, namely a \(\sqrt{3}\times \sqrt{3}\) phase with respect to the underlying Si surface, a 3 × 3 phase, and a \(\sqrt{7}\times \sqrt{3}\) phase. Recent findings indicate the presence of superconductivity in the \(\sqrt{7}\times \sqrt{3}\) at low temperatures38,39 as well as chiral superconductivity in the Si(111):Sn system40. Unknown superconductive phases could appear also in other surface reconstructions of Si(111):Pb, likely coexisting with magnetic phases or CDW. The lattice distortion would likely play a crucial role as well. Due to the similarities with Si(111):Sn, we would expect the chirality to be still present. In addition to that, the strong spin-orbit coupling could even lead to more exotic forms of superconductivity41.

There exists a compelling evidence that the system with 1/3 coverage exhibits a structural transition to a 3 × 3 charge density wave (CDW) phase at a temperature of 86 K42,43,44. It is still a matter of ongoing research to understand whether this transition has to be attributed to a Peierls-like mechanism, an intrinsic asymmetry induced by the interaction with the substrate, or to strong electronic correlations, as claimed in ref. 37. Remarkably, a similar transition takes place in Ge(111):Pb and Ge(111):Sn16,45, but not in the Si(111):Sn compound. In addition, by using scanning tunneling microscopy (STM) it has been found that the quasiparticle interference patterns are influenced by the strong value of the SOC giving rise to a chiral spin structure at low temperatures inside the CDW phase15,37.

The experimental study of the low-temperature 3 × 3 phase of Si(111):Pb is technically non-trivial44. This phase is difficult to grow as an extended phase limiting the experimental probes that can be used. For this reason, the investigation of this system has so far been limited to STM experiments. In order to perform STM measurements, it is necessary to induce a finite conductance in the system. To this aim, slightly doped substrates have to be used15,37. Depending on the doping level of the substrate, adatoms can exert an attractive or repulsive force on the impurities in the bulk, that can strongly affect the doping level of the surface band46,47. As a consequence, the use of Si substrates with strong electron-doping15 or hole-doping37 can induce a significant doping on the surface states. Additional data with accurate experimental control over doping conditions would be crucial to shed a light on the observed phases, since the doping in the system may strongly affect collective electronic effects and related phases. The effect of doping could also explain a crucial difference between theoretical results and experiments. Calculations with correlated theories predict a Mott-insulating behavior, while the measured STM spectrum is metallic15,37. This apparent contradiction may be explained by noting that in a Mott insulator, an arbitrarily small doping level can induce a metallic behavior. For this reason, a careful investigation of the temperature vs doping phase diagram is absolutely necessary to explain the experimentally observed effects in Si(111):Pb.

In this work, we use advanced many-body techniques to analyse collective electronic effects in Si(111):Pb as a function of temperature and doping. We find a very rich phase diagram comprising charge- and spin-density wave phases characterized by different ordering vectors. By comparing results for the \(\sqrt{3}\times \sqrt{3}\) and 3 × 3 structures, we find that different CDW orderings can originate from either a structural transition due to an asymmetric interaction of adatoms with the substrate, or from strong electronic correlations depending on the doping level. Further, we observe that the spin ordering in the system also depends on the doping. These results illustrate that varying the doping level in the Si(111):Pb material represents an efficient way of switching between different CDW and magnetic phases. In addition, we argue that a simultaneous detection of the charge- and spin-density orderings in an experiment can help to understand in which part of the complex temperature vs doping phase diagram the measured system is located.

Results

Model

According to density functional theory (DFT) calculations, the Si(111):Pb system with 1/3 coverage in the high-temperature \(\sqrt{3}\times \sqrt{3}\) phase (Fig. 1) exhibits a narrow half-filled band at the Fermi level, well separated from the rest of the bands26,27,28. In the maximally localized Wannier basis, this band has a pz character, and the corresponding Wannier orbitals are centered at the Pb adatom sites. We thus employ the following single-band interacting electronic model derived from the first-principle DFT calculations:

$$\begin{array}{ll}\hat{H}=\mathop{\sum}\limits_{ij,\sigma {\sigma }^{{\prime} }}{c}_{i\sigma }^{{\dagger} }\left({t}_{ij}\,{\delta }_{\sigma {\sigma }^{{\prime} }}+i\,{{{{\boldsymbol{\gamma }}}}}_{ij}\cdot {{{{\boldsymbol{\sigma }}}}}_{\sigma {\sigma }^{{\prime} }}\right){c}_{j{\sigma }^{{\prime} }}+\mathop{\sum}\limits_{i}{\Delta }_{i}{n}_{i}+U\mathop{\sum}\limits_{i}\,{n}_{i\uparrow }{n}_{i\downarrow }\\\quad+\frac{1}{2}\mathop{\sum}\limits_{i\ne j,}{V}_{ij}\,{n}_{i}{n}_{j}+\frac{1}{2}\mathop{\sum}\limits_{i\ne j}{J}_{ij}\,{{{{\bf{S}}}}}_{i}\cdot {{{{\bf{S}}}}}_{j}.\end{array}$$
(1)

In this expression, \({c}_{i\sigma }^{({\dagger} )}\) corresponds to an annihilation (creation) operator for an electron on the lattice site i with the spin projection σ ∈ {↑, ↓}. tij corresponds to the hopping amplitude between i and j lattice sites, while Δi indicates the local on-site potentials. The considered Hamiltonian accounts for the SOC in the Rashba form48,49 of a spin-dependent imaginary hopping \({{{{\boldsymbol{\gamma }}}}}_{ij}={\gamma }_{| i-j| }\,({\hat{r}}_{ij}\times \hat{z})\). The Coulomb interaction between electronic densities ni = ∑σniσ, where \({n}_{i\sigma }={c}_{i\sigma }^{{\dagger} }{c}_{i\sigma }\), is explicitly divided into the local U and the non-local Vij parts. Jij represents the direct ferromagnetic exchange interaction between the magnetic densities \({{{{\bf{S}}}}}_{i}={\sum }_{\sigma {\sigma }^{{\prime} }}{c}_{i\sigma }^{{\dagger} }{{{{\boldsymbol{\sigma }}}}}_{\sigma {\sigma }^{{\prime} }}{c}_{i{\sigma }^{{\prime} }}\), where σ = {σx, σy, σz} is a vector of Pauli matrices.

Fig. 1: Crystal structure of Si(111):Pb.
figure 1

a Top view of the Si(111):Pb surface reconstruction. Shown is one 3 × 3 unit cell. b Side view of the slab geometry adopted in the DFT structural relaxations consisting of 1/3 monolayer of Pb adatoms on top of three Si bi-layers, bottom terminated by hydrogen capping. The uniform distance between the Pb adatoms and the substrate corresponds to a high-temperature \(\sqrt{3}\times \sqrt{3}\) phase. In the low-temperature 3 × 3 reconstruction the Pb adatoms form a “1-up-2-down” configuration with respect to the substrate.

Full size image

In momentum-space, one can write the Fourier transform of the hopping amplitudes as \({\varepsilon }_{{{{\bf{k}}}},l{l}^{{\prime} }}^{\sigma {\sigma }^{{\prime} }}={t}_{{{{\bf{k}}}},l{l}^{{\prime} }}{\delta }_{\sigma {\sigma }^{{\prime} }}+i\,{\overrightarrow{\gamma }}_{{{{\bf{k}}}},l{l}^{{\prime} }}\cdot {\overrightarrow{\sigma }}_{\sigma {\sigma }^{{\prime} }}\) with \({l}^{({\prime} )}\) labeling nonequivalent lattice sites within the unit cell. Further, we focus on the two distinct structures of the Si(111):Pb material. In the high-temperature \(\sqrt{3}\times \sqrt{3}\) structure the Pb adatoms form a triangular lattice with identical lattice sites, so we set \(l={l}^{{\prime} }\). Upon decreasing the temperature, the system undergoes a structural transition, which results in a 3 × 3 reconstruction of the adatoms. The resulting structure has the form of an effective triangular lattice, but the unit cell contains three Pb atoms. Lattice relaxations within the generalized gradient approximation (GGA) and experiments show that these three Pb atoms display a corrugated “1-up-2-down” configuration with respect to a flat surface42,43,44. We find that a local potential Δl with l ∈ {1, 2, 3} is sufficient to describe the position of nonequivalent sites within the unit cell. This potential is set to zero in the \(\sqrt{3}\times \sqrt{3}\) structure, while it is non-zero in the 3 × 3 structure because of the substrate-induced deformation, which corresponds to a static electron-phonon interaction18. In this regard, the high-temperature \(\sqrt{3}\times \sqrt{3}\) phase can be seen as a time-averaged 3 × 3 structure, due to dynamical fluctuations of the adatom height50,51. The values of all model parameters and details of the DFT calculations are given in the Methods section.

Detection of collective electronic instabilities

Instabilities related to collective electronic fluctuations in the charge (c) and spin (s) channels can be detected via the momentum-dependent static structure factor (see, e.g., refs. 52,53,54)

$$\begin{array}{r}{S}^{c/s}({{{\bf{q}}}})=\mathop{\sum}\limits_{l{l}^{{\prime} }}{e}^{-i{{{\bf{q}}}}\cdot ({{{{\bf{R}}}}}_{l}-{{{{\bf{R}}}}}_{{l}^{{\prime} }})}{X}_{l{l}^{{\prime} }}^{c/s}({{{\bf{q}}}},\omega =0),\end{array}$$
(2)

where the vector Rl depicts the position of the atom l within the unit cell. In the high-temperature \(\sqrt{3}\times \sqrt{3}\) phase, where \(l={l}^{{\prime} }\), the static structure factor coincides with the static susceptibility Xc/s(q, ω = 0) obtained at zero frequency ω. The divergence of the structure factor at momenta q = Q indicates a transition to a symmetry-broken ordered state associated with Bragg peaks at Q. Transitions without symmetry-breaking, such as the metal to Mott insulator phase transition, can be observed by inspecting the spectral function. In this work, the introduced many-body problem (1) is solved using the dual triply irreducible local expansion (D-TRILEX) method55,56,57. This method provides a consistent treatment of the local correlation effects and the non-local collective electronic fluctuations in the charge and spin channels58,59,60,61. Importantly, D-TRILEX is also able to account for the long-range Coulomb interaction59 and the SOC57, which are the two important aspects of the considered material. More details on the many-body calculations are provided in the Methods section.

Phase diagram for the \(\sqrt{3}\times \sqrt{3}\) structure

The phase diagram for the Si(111):Pb material in the \(\sqrt{3}\times \sqrt{3}\) structure is shown in Fig. 2 as a function of doping level δ and temperature T. In the considered system, the value of the local Coulomb interaction is approximately 3 times larger than the electronic bandwidth27,28. As a consequence, at high temperature the half-filled system lies deep in the Mott insulating phase (black line at δ = 0%). A small amount of hole- or electron-doping causes a phase transition to a Fermi liquid regime (gray area). For this reason, the electronic behavior in doped Si(111):Pb is a characteristic manifestation of the physics of a doped Mott insulator.

Fig. 2: Phase diagram for Si(111):Pb in the \(\sqrt{3}\times \sqrt{3}\) structure.
figure 2

Different phases as a function of doping δ and temperature T are highlighted by colors (color code in the legend). Calculations have been performed by fixing the temperature and conducting a scan over doping levels on a finite grid, which defines the error bars. Positive (negative) values of δ correspond to electron (hole) doping. The horizontal dashed black line depicts the temperature T = 86 K at which the material exhibits a structural phase transition according to refs. 42,43. The vertical line that divides the magnetic phases below the transition points is only meant as a guide to the eye, since we are not able to distinguish between the two different phases in symmetry-broken regime.

Full size image

Upon solving the many-body problem (1) we identify several different spin-density wave (SDW) and CDW orderings at different values of doping, as illustrated in Fig. 2. Since these phases are realized for a non-integer filling of electrons, they are likely metallic. However, we cannot confirm this in our actual calculations because our method does not allow us to perform calculations inside phases induced by dynamic symmetry breaking. Specifically, around half-filling, we observe a CDW ordering (orange area around δ = 0%) characterized by the divergence of the static charge structure factor at the Q = K point of the Brillouin zone (BZ). This ordering is analogous to the 120∘-Néel phase of the Heisenberg model on a triangular lattice with three inequivalent sites in the unit cell (see, e.g., ref. 62). For this reason, hereinafter we call this type of ordering a “tripartite CDW”. Importantly, we find that this instability does not appear if instead of the full long-range Coulomb potential Vij one considers the interaction only between nearest-neighbor lattice sites. In the presence of only local interactions, the Mott phase and a CDW would be mutually exclusive. Here, we note that the effective long-range interaction is enhanced by correlations as the temperature is reduced, while the local interaction barely depends on temperature. We would also like to note that competing tripartite CDW and Mott phases have been observed experimentally in the other adatom system Ge(111):Sn63.

Additionally, we identify two other CDW phase transitions at dopings around δ = ± 10%. These instabilities appear to be weakly temperature-dependent and approximately symmetric with respect to half-filling. At hole doping, the CDW ordering vector remains Q = K (orange area), as in the half-filled case. However, in the electron-doped regime the divergence of the static charge structure factor occurs at the Q = M point of the BZ, which can be associated with a “row-wise CDW” ordering (red area). One can speculate, that this ordering might be related to the isoelectronic mosaic phase observed in Si(111):Pb64 or to the intermediate stripe-like order in the alkali-doped Si(111):Sn surface65. However, a direct observation of the row-wise CDW phase in Si(111):Pb has not been performed yet. The momentum-resolved static charge structure factors obtained close to both these CDW instabilities are shown in Fig. 3, where the Bragg peaks clearly indicate the corresponding ordering vectors.

Fig. 3: The static charge structure factor Sc(q).
figure 3

The result is obtained close to the CDW transition points δ = − 7%, T = 25 K (a) and δ = 10%, T = 67 K (b). In the hole-doped case, the Bragg peaks in the structure factor appear at the Q = K points of the BZ indicating the tripartite CDW ordering. In the electron-doped case, the ordering vector Q = M corresponds to the row-wise CDW instability.

Full size image

In addition to the CDW instabilities, we also observe magnetic structures with different ordering vectors depending on the doping level (cyan and blue areas in Fig. 2). Around half-filling, we observe a SDW characterized by Bragg peaks in the static spin structure factor that lie at an incommensurate point \({{{\bf{Q}}}}\simeq \frac{2}{3}\,{{{\rm{M}}}}\) of the BZ (Fig. 4a). At δ ≳ 2% of electron-doping the SDW ordering vector changes, and the peaks shift to another incommensurate position \({{{\bf{Q}}}}\simeq \frac{3}{4}\,{{{\rm{K}}}}\) (Fig. 4c). The appearance of the Bragg peaks at incommensurate points of the BZ signals the formation of a chiral magnetic order that can be viewed as a superposition of spin spirals. According to the position of the Bragg peaks, we call these magnetic structures “chiral-M” (cyan area) and “chiral-K” (blue area) SDW, respectively. The presence of the chiral magnetic orderings in Si(111):Pb suggests that this material might be a suitable candidate for the realization of skyrmionic phases that can possibly be stabilized under an external magnetic field28.

Fig. 4: The static spin structure factor Ss(q).
figure 4

The results are obtained with (left column) and without (right column) SOC, respectively. The upper row corresponds to the half-filling δ = 0%, the bottom row to δ = 7.4% electron doping. The chosen temperature, T = 67 K, is close to the SDW transition. Without SOC, the Bragg peaks in Ss(q) indicate the row-wise (b) and the Néel (d) magnetic structures. Taking into account the SOC, the Bragg peaks in both cases shift to incommensurate positions with chiral-M (a) and chiral-K (c) SDW orderings.

Full size image

Remarkably, the obtained chiral SDW structures partially coexist with the CDW orderings. In the considered Si(111):Pb material, such coexistence was recently observed by means of STM measurements15, but an estimate of the doping level in the system was not provided, presumably due to difficulties in the determination of the effective doping. Remarkably, we find that the chiral-M SDW structure coexists only with the tripartite CDW ordering, which appears around half-filling. Instead, the row-wise CDW ordering coexists only with the chiral-K SDW at a relatively large electron doping. This observation suggests a simple way for a qualitative estimation of the doping level in the experimentally measured material, which is difficult to probe directly (see refs. 46,47 and related supplemental materials for discussion).

We have made a very crude estimation of the doping level by calculating the area of the Fermi surface that can be deduced from the STM map shown in ref. 15. The obtained result is compatible with up to ≃11% electron-doping, which coincides with the region of coexisting chiral-K SDW and row-wise CDW orderings. This result appears to be consistent with the use of an electron-doped substrate15.

Effect of the SOC

We observe that the large SOC, which is an intrinsic feature of Si(111):Pb, manifests itself in the magnetic properties of the material. In particular, the effect of the SOC can be seen in the spin structure factors shown in Fig. 4. As we have shown above, the SOC results in the formation of the chiral-M (a) and chiral-K (c) SDW orderings in the system. Instead, if the SOC is not taken into account, the Bragg peaks in the static spin structure factor calculated close to the SDW phase transitions appear at the Q = M (b) and Q = K (d) points of the BZ. These instabilities correspond to commensurate row-wise and Néel magnetic structures, respectively. Remarkably, despite the shift of the peaks in the BZ and the consequent change of the ordering of the system, we find that the position of the phase boundaries is not affected by the SOC (up to the error bars of our calculations), similarly to what has been found in ref. 66 for a square lattice. Based on this result, one can argue that the phase boundaries in the considered system can be obtained correctly without taking into account the SOC. However, considering the SOC is absolutely necessary for an accurate determination of the ordering vectors.

Effective Heisenberg model

The observed changes in the spin structure factor as a function of doping level can be explained by analyzing the exchange interactions67,68,69. These quantities are accessible in D-TRILEX calculations57. To this aim, we consider the following effective Heisenberg-like classical spin Hamiltonian with bilinear magnetic exchange interactions:

$$\begin{array}{r}H=J\mathop{\sum}\limits_{\langle ij\rangle }{{{{\bf{S}}}}}_{i}\cdot {{{{\bf{S}}}}}_{j}+{J}^{{\prime} }\mathop{\sum}\limits_{\langle \langle ij\rangle \rangle }{{{{\bf{S}}}}}_{i}\cdot {{{{\bf{S}}}}}_{j}+D\mathop{\sum}\limits_{\langle ij\rangle }{{{\bf{z}}}}\cdot \left({{{{\bf{S}}}}}_{i}\times {{{{\bf{S}}}}}_{j}\right).\end{array}$$
(3)

In this expression, J and \({J}^{{\prime} }\) are the nearest-neighbor 〈ij〉 and the next-nearest-neighbor 〈〈ij〉〉 exchange interactions, respectively. D is the nearest-neighbor Dzyaloshinskii–Moriya interaction (DMI), which appears due to the SOC. We have also calculated the symmetric anisotropy, but we omit it for simplicity as it hardly affects the following considerations. The value of its only non-zero component is Γyy ≈ 0.5D in the whole range of δ considered here.

Figure 5 shows the evolution of \({J}^{{\prime} }\) and D, normalized by the value of J, as a function of doping. Remarkably, we find that the magnitude of D in Si(111):Pb is of the order of the nearest-neighbor exchange interaction J, which is very unusual for magnetic systems. Moreover, D and J even become equal in the electron-doped case. At half-filling the value of D/J coincides with the one obtained in ref. 28 using the strong-coupling approximation. This fact confirms that the half-filled Si(111):Pb material lies in the strong-coupling regime. Further, we observe that the ratio D/J has an approximately linear dependence on doping with different slopes in the hole- and electron-doped regimes. In the hole-doped case, D/J substantially decreases upon increasing the doping. Instead, in the electron-doped regime, D/J slowly increases with increasing δ. This behavior explains the formation of the chiral SDW orderings in the regime of doping levels δ ≳ −7%, where DMI is strong enough (D/J ≳ 0.4) to be able to shift the Bragg peaks from a commensurate to an incommensurate position, as shown in Fig. 4.

Fig. 5: Magnetic exchange interactions as a function of doping.
figure 5

The orange line depicts the value of the nearest-neighbor Dzyaloshinskii–Moriya interaction D/J. The blue line corresponds to the next-nearest-neighbor exchange interaction \({J}^{{\prime} }/J\). Both quantities are normalized by the value of the nearest-neighbor exchange J. The results are obtained at T = 50 K. The black hexagon and the black dot represent the values of D/J and \({J}^{{\prime} }/J\) obtained at half-filling in ref. 28 using the strong-coupling approximation. The vertical dashed black line at δ = 1.8% indicates the transition from the chiral-M to the chiral-K phases according to our calculations. The horizontal dashed line at \({J}^{{\prime} }/J=0.12\) represents the prediction for the M to K transition in the J-\({J}^{{\prime} }\) Heisenberg model obtained from Monte Carlo calculations in ref. 70.

Full size image

While DMI is responsible for the formation of chiral spin structures, the change in the ratio \({J}^{{\prime} }/J\) with doping explains the transformation of the magnetic ordering from the M- to the K-type, as observed in our calculations. The magnitude of \({J}^{{\prime} }\) is rather small compared to J and D, but it is not negligible. In addition, we find that the actual value of the more distant, next-nearest-neighbor exchange interaction \({J}^{{\prime} }\) is substantially larger than the one predicted by a strong-coupling estimate28. An important feature is that the ratio \({J}^{{\prime} }/J\) is nearly constant in the hole-doped regime, while in the electron-doped case it substantially decreases and even changes sign. We attribute this variation of \({J}^{{\prime} }/J\) to the shift of the Bragg peaks in the spin structural factor from M to K, which is consistent with Monte Carlo calculations for the J-\({J}^{{\prime} }\) Heisenberg model on a triangular lattice performed in ref. 70. It has been shown there, that the transition from a row-wise (Q = M) to a Néel (Q = K) magnetic order occurs for \({J}^{{\prime} }/J\simeq 0.12\). As shown in Fig. 5, this result coincides with our estimate for the transition point between the chiral-M to chiral-K SDW orderings. In this figure, the horizontal dashed black line depicts the \({J}^{{\prime} }/J=0.12\) value, and the vertical dashed black line marks the mean-point between the closest doing levels that correspond to chiral-M and chiral-K SDW orderings.

Phase diagram for the 3 × 3 reconstruction

At low temperature, Si(111):Pb undergoes a structural phase transition from \(\sqrt{3}\times \sqrt{3}\) to 3 × 3 periodicity. The 3 × 3 reconstruction exhibits a 1-up-2-down configuration of Pb adatoms, as confirmed in experiments42,43 and by DFT calculations15,33. In order to account for the effect of the structural phase transition, we also perform many-body calculations for the 3 × 3 reconstruction of adatoms. The 1-up-2-down configuration requires to consider a unit cell with three Pb atoms, which significantly increases the cost of the numerical calculations. As previously discussed, the inclusion of the SOC does not affect the position of the phase boundaries in the considered material. In order to make numerical calculations in the 3 × 3 phase feasible, we neglect the Rashba term in the model Hamiltonian (1).

Figure 6 shows the resulting phase diagram for the 3 × 3 reconstruction, which qualitatively agrees with the one obtained for the \(\sqrt{3}\times \sqrt{3}\) structure. Indeed, the phase diagram for the 1-up-2-down configuration of Pb atoms also contains row-wise and tripartite CDW phases that are nearly temperature-independent and appear at values of the hole- and electron-doping comparable to the \(\sqrt{3}\times \sqrt{3}\) case. We note that these dynamical CDW instabilities emerge on top of the structural phase transition, which affects the ordering vector of the row-wise CDW structure. Indeed, Fig. 7b shows that the Bragg peaks in the charge structure factor are now found at incommensurate positions in the vicinity of the M point of the BZ. This result can be explained by the observation that the divergence of the corresponding charge susceptibility \({X}_{l{l}^{{\prime} }}^{c}({{{\bf{q}}}},\omega =0)\), which enters the expression (2) for the structure factor, also appears at incommensurate positions in the vicinity of the M point of the reduced BZ. A wave-vector at the M point would mean row-wise ordering, as in the single-site case. However, here we have two overlapping orderings: a row-wise order induced by correlations and the 3 × 3-underlying broken symmetry due to the lattice distortion. The reason for this pattern is that a perfect row-wise arrangement would not be commensurate with the underlying 1-up-2-down structure. It means that the spontaneous symmetry breaking leading to the row-wise CDW ordering occurs between different unit cells on the lattice, but not within the unit cell of three Pb atoms. On the contrary, we find that the ordering vector Q = K of the tripartite CDW instability remains unchanged upon the structural transition (top left panel of Fig. 7). The tripartite CDW corresponds to the ordering, where all three Pb atoms in the unit cell are inequivalent. The fact that upon the tripartite CDW phase transition the charge susceptibility diverges at the Γ point of the reducible BZ confirms the statement that, in this case, the spontaneous symmetry breaking occurs within the unit cell. Consequently, the 1-up-2-down structure of Pb atoms in the unit cell transforms to a tripartite structure, and the Bragg peaks in the structural factor appear at the K point of the BZ as usual.

Fig. 6: Phase diagram for Si(111):Pb in the 3 × 3 reconstruction.
figure 6

Different phases that appear in the system as a function of the doping δ and temperature T are highlighted in color. The color code can be found in the legend. Calculations have been performed by fixing T and conducting a scan over doping levels on a finite grid, which defines the error bars. The vertical line dividing the magnetic phases below the transition points are only meant as a guide to the eye.

Full size image
Fig. 7: Charge (top row) and spin (bottom row) static structure factors.
figure 7

The depicted results are for the 3 × 3 reconstruction in the vicinity of the tripartite CDW (a), row-wise CDW (b), M SDW (c), and K SDW (d) phase transitions. The corresponding doping levels are specified for each panel. The temperature is chosen to be close to the phase boundaries, T ≈ 25 K in panel (a), T ≈ 35 K in panel (c), and T ≈ 67 K in (b) and (d).

Full size image

The structural transition also affects the phase boundaries of the temperature-dependent instabilities. All of them, namely the CDW around half-filling and both SDW instabilities, are pushed down to lower temperatures. This can be related to the appearance of an effective local potential Δl upon the structural transition to the 1-up-2-down structure. This potential acts as an on-site doping that differs from site to site and thus suppresses collective charge and spin fluctuations. Interestingly, the CDW ordering found around half-filling in the 3 × 3 reconstruction has a row-wise structure instead of the tripartite one observed in the \(\sqrt{3}\times \sqrt{3}\) case. As discussed above, the row-wise ordering does not break the 1-up-2-down structure of Pb adatoms in the unit cell. Probably for this reason, the formation of the row-wise CDW is more favorable in the 3 × 3 phase. Finally, we note that apart from decreasing the critical temperature for the SDW instabilities, the structural transition does not affect the magnetic ordering in the system. As in the \(\sqrt{3}\times \sqrt{3}\) case we find the M SDW ordering around half-filling and the K SDW ordering at δ ≳ 2% of electron-doping. In our calculations, the Bragg peaks in the corresponding spin structure factors appear at commensurate Q = M (top left panel of Fig. 7) and Q = K (top right panel of Fig. 7) positions. We expect that the inclusion of the SOC would shift the peaks to incommensurate positions and lead to the formation of the chiral magnetic structures also in the 3 × 3 case.

Discussion

We performed many-body calculations for a system of Pb adatoms on a Si(111) substrate, including the SOC and long-range Coulomb interactions. By investigating spatial collective electronic fluctuations in both, charge and spin channels, we observe a rich variety of different symmetry-broken charge- and spin-density wave phases in the low-temperature regime by varying the doping level. Regarding the Mott physics, our results show a picture similar to that of Sn on Si(111): the system is a Mott insulator at half-filling, but immediately turns into a metal as soon as some small doping is introduced in the system47. We find that the strong SOC in this material results in a very large Dzyaloshinskii–Moriya interaction comparable to the usual Heisenberg exchange interaction. This leads to the formation of chiral-M and chiral-K SDW phases, a signature of which have recently been observed in STM measurements15. These chiral spin structures are compatible with magnetic skyrmion textures, as highlighted in previous theoretical calculations28. Tuning the doping level allows one to switch between the two chiral SDW phases and thus realize different kinds of spin structures with potential topological structure in one material. We note that a similar change of the magnetic ordering was proposed for a Si(111):Sn system by means of varying the local Coulomb interaction36.

We also find that two different CDW orderings can appear in Si(111):Pb, and that their geometry is strongly affected by the doping level. The values of doping, at which the transition takes place, appear to be consistent with the intrinsic doping levels observed in this kind of systems46. There is an ongoing debate whether the 3 × 3 pattern of charge densities observed in experiments emerges in Si(111):Pb due to a dynamical symmetry breaking associated with strong electronic correlations37, or by means of a structural transition18. We argue that the corresponding 1-up-2-down structure of Pb adatoms can be realized in the system upon either the structural transition from the \(\sqrt{3}\times \sqrt{3}\) to the 3 × 3 phase, or the dynamical symmetry breaking towards the row-wise CDW phase, depending on the doping level and temperature. In addition, we find another CDW ordering in the system associated with the formation of a tripartite structure.

In order to realize these theoretically predicted phases in the experiment, it is necessary to use a probe sensitive to collective excitations, as well as to be able to give an accurate estimation of the doping level. Since the precise occupation of the isolated band is experimentally challenging to access, we propose an alternative way to identify the doping level. Using a probe sensitive to the underlying magnetic structure, such as spin-polarized STM71, could prove a valid alternative to the measurements of the doping, since the magnetic textures appearing at different doping levels exhibit different geometry and also coexist with different types of CDW ordering.

A recent study on a similar adatom system of Sn adatoms on germanium indicated the presence of strong electron-phonon coupling (EPC)72. This system has a different composition, so it is not known if a similar effect holds also for Pb on Si(111). We argue that the EPC scales as \(1/\sqrt{M}\) with the atomic mass M, so the contribution to the effective electron-electron interaction scales as 1/M and it is much smaller on the Pb surface than in the case of Sn. Additionally, in order to strongly affect the properties of the system, EPC would need to overcome the very strong Coulomb interaction present in this system. As this is very unlikely to occur, we conclude that we do not expect this contribution to be crucial to determine the phases of this system. However, it could modify the position of the phase boundaries, hence in the future, it would be desirable to devise a way to deal with EPC in D-TRILEX calculations. Further studies are also required in order to investigate superconductivity in the low-temperature regime.

Methods

Ab-initio DFT calculations

All model parameters used in the model Hamiltonian (1) have been obtained from ab-initio calculations. For the \(\sqrt{3}\times \sqrt{3}\) structure of adatoms, we adapted the parameters from ref. 28, where a Wannier projection on localized orbitals was performed to obtain the nearest-neighbor t01 = 41.3 meV and the next-nearest-neighbor t02 = − 19.2 meV hopping amplitudes. The Rashba parameters γ01 = 16.7 meV and γ02 = 2.1 meV are taken from the same work as the hopping amplitudes. The value of the local Coulomb interaction U = 0.9 eV is the one obtained from cRPA calculations26,27,28. The long-range Coulomb interaction with a realistic 1/r tail is parametrized by the nearest-neighbor interaction V01 = 0.5 eV as suggested in refs. 26,27,73. The direct exchange interaction between neighboring sites that enters Eq. (1) is rather small and reads J〈ij〉 = 1.67 meV28.

For the 3 × 3 reconstruction, we simulated the surface by a slab geometry consisting of 1/3 monolayer of Pb adatoms on top of three Si bi-layers, as established in previous works15,26,28,35,74. The Pb adatoms occupy the T4 positions. The dangling bonds of the bottom Si bi-layer are compensated by hydrogen capping, and 19 Å of vacuum are included in the simulation. For structural relaxations, we employ the WIEN2k75,76 program package, a full-potential linearized-augmented plane-wave code. We start with the relaxation of the \(\sqrt{3}\times \sqrt{3}\) structure, which contains one Pb per unit cell. We then construct the 3 × 3 supercell containing 3 Pb atoms (66 atoms in total, thereof 54 Si). To relax the 3 × 3 structure, which in the experiment is found in a 1-up-2-down configuration, we displace one of the three Pb adatoms by 0.4 Å perpendicularly to the surface in the first DFT self-consistent-field iteration. We then let the internal coordinates of all atoms in the supercell relax freely until convergence. We employed a multisecant approach76, as implemented in WIEN2k75,76. A k-grid with 6 × 6 × 1 k-points in the reducible Brillouin zone was used and internal coordinates were relaxed until forces were less than 2 mRy/bohr. We employed the generalized gradient approximation (PBE), spin-orbit coupling was neglected. In agreement with the experiment, we find the stabilization of a 3 × 3 reconstruction, where one Pb adatom is vertically displaced by 0.22 Å compared to the other two Pb adatoms in the supercell. The energy gain of this 1-up-2-down reconstruction is found to be 9.5 meV with respect to a flat adatom layer. These findings are in good agreement with previous ab-initio calculations15,42. We find that the computed band structure for the 3 × 3 reconstruction can be well interpolated with a three-band dispersion using the same parameters taken from ref. 28 by simply adding a local potential Δl to the inequivalent Pb atoms l ∈ {1, 2, 3} in the model Hamiltonian. We choose this approach to ensure better comparability between the calculations. The obtained values for the potential are Δ1 = Δ2 = 31.5 meV and Δ3 = − 55.4 meV. The effect of the substrate-induced deformation, which corresponds to a static electron-phonon interaction, can be a crucial ingredient for the formation of the 3 × 3 structure18. We stress that this effect of phonons is taken into account in our calculations of the 3 × 3 structure by keeping the lattice distortion appearing at the DFT level in the interacting problem (1).

Many-body D-TRILEX calculations

The interacting electronic problem (1) is solved using the finite temperature D-TRILEX method55,56,57. To this aim, we first perform converged dynamical mean-field theory (DMFT) calculations77 with the w2dynamics package78 in order to take into account local correlation effects in a numerically exact way. Furthermore, the effect of the non-local collective electronic fluctuations and of the SOC is taken into account diagrammatically as described in ref. 57. The spin susceptibility \({X}_{l{l}^{{\prime} }}^{s}({{{\bf{q}}}},\omega )\) required for the calculation of the structure factor (2) is defined as the maximum eigenvalue of the matrix

$$\begin{array}{r}{X}_{l{l}^{{\prime} }}^{s{s}^{{\prime} }}({{{\bf{q}}}},\omega )=\left\langle {S}_{l{{{\bf{q}}}}\omega }^{s}{S}_{{l}^{{\prime} },-{{{\bf{q}}}},-\omega }^{{s}^{{\prime} }}\right\rangle \end{array}$$
(4)

in the space of spin channel indices \({s}^{({\prime} )}\in \{{s}_{x},{s}_{y},{s}_{z}\}\). The charge susceptibility is defined as:

$$\begin{array}{r}{X}_{l{l}^{{\prime} }}^{c}({{{\bf{q}}}},\omega )=\langle {n}_{l{{{\bf{q}}}}\omega }{n}_{{l}^{{\prime} }-{{{\bf{q}}}},-\omega }\rangle .\end{array}$$
(5)

Note that in this work, the susceptibility is computed non-self-consistently, as, e.g., in ref. 59. This means that the susceptibility is calculated on the basis of the electronic Green’s functions dressed only by the local DMFT self-energy, which resembles the way the susceptibility is computed in DMFT56. This procedure allows one to treat collective electronic instabilities in the charge and spin channels independently without mutually affecting each other.

The magnetic exchange interactions used to construct the effective Heisenberg model (3) are also computed within the D-TRILEX scheme, as explained in ref. 57.