Introduction

The discovery of Weyl semimetals, which host spin-split massless quasiparticles in three-dimensional (3D) crystals, is particularly thrilling as it represents an experimental realization of Weyl fermions, a concept proposed long ago in the realm of particle physics1,2,3,4,5. The chiral nodal points and 2D Fermi arc surface states of 3D Weyl semimetals bring about exotic properties such as chiral anomaly, unusual optical conductivity and nonlocal transport6,7,8,9,10,11,12,13,14,15,16. While the Weyl equation was derived only for odd spatial dimensions, the generalization of a 3D Weyl fermion state in 2D leads to a distinct topological state of matter, labeled as 2D Weyl semimetals, that exhibit a spin-polarized analog of graphene. The dimension reduction gives rise to a multitude of unconventional physical properties, including parity anomaly in (2 + 1)-D (space-time) quantum field theory17,18,19,20, charge fractionalization with zero modes of charge e/221, spin-valley Hall effects22,23, giant Berry curvature dipole (BCD)24,25, and topological quantum criticality22. Among the various topological Dirac/Weyl solid-state materials (see Fig. 1a), 2D Weyl semimetals stand out as the final frontier that has yet to be extensively explored in experiments. A winding number of π can be obtained by integrating the Berry phase along a loop encircling each Weyl node26. The nonzero winding number can be regarded as the topological charge of 2D Weyl semimetals, which guarantees the existence of topologically protected edge states27. These topological edge states take the form of Fermi strings with one end attached to the projection of bulk Weyl nodes at the Fermi level, as schematically shown in Fig. 1a28. The Fermi string edge states serve as the 1D counterparts to the Fermi arc surface states observed in 3D Weyl semimetals. In this context, 2D Weyl semimetals present an unprecedented paradigm of bulk–boundary correspondence in topological materials. The highly unusual properties associated with 2D Weyl fermion states have inspired a myriad of theoretical and experimental works20,28,29,30,31,32,33,34,35,36,37. Nevertheless, the realization of an intrinsic 2D Weyl semimetal in experiments remains elusive to date.

Fig. 1: Topology, lattice property, and band structure of 2D Weyl semimetal, α-bismuthene grown on SnS.
figure 1

a Overview of Dirac/Weyl semimetals and their topological boundary states. b Side and top views of the lattice structure of bismuthene and SnS substrate. The red dashed squares indicate the unit cell of bismuthene and SnS. c Large-scale STM image of bismuthene grown on SnS substrate. d The height profile is taken along the red arrow in (c). e Zoom-in STM images of bismuthene (top) and the surface of SnS (bottom). The red dashed squares indicate the unit cell of SnS and bismuthene. f Calculated band structure of free-standing bismuthene. g Calculated band structure of bismuthene on SnSe.

In this work, we report the realization of a 2D Weyl semimetal in an intrinsic 2D crystal, epitaxial bismuthene (a single atomic layer of bismuth stabilized in phosphorene structure38,39). Our spin- and angle-resolved photoemission spectroscopy (spin-ARPES) results unambiguously demonstrate the spin-polarized Weyl fermion states in bismuthene grown on SnS(Se) substrates. Furthermore, our scanning tunneling spectroscopy (STS) measurements unveil a robustly enhanced local density of states at the edge of bismuthene, which aligns seamlessly with the calculated edge spectrum featuring Fermi string states. These experimental results establish epitaxial bismuthene on SnS(Se) as an ideal 2D Weyl semimetal.

Results and discussion

The material base for hosting 2D Weyl fermion states

Each Bismuth atom typically forms three covalent bonds with its neighbors. Two allotropic structural phases of bismuthene exist in the 2D limit, namely, the orthorhombic phosphorene-like phase40,41 and the hexagonal honeycomb-like phase42,43. Here we grow the phosphorene-like bismuthene (bismuthene for short in the following) by molecular beam epitaxy (MBE). We selected SnS(Se) as the substrate due to its van der Waals semiconductor properties, and the (001) surface of SnS(Se) shares a lattice structure similar to that of bismuthene. Figure 1b schematically shows the lattice structure of the sample. Bi atoms form a single-layer phosphorene structure on the (001) surface of SnS(Se). The Bi atoms exhibit a pronounced sp3-hybridization, resulting in the three Bi–Bi bonds adopting a nearly tetrahedral configuration. This configuration leads to two atomic sublayers, which are distinguished by red and orange colors in Fig. 1b. The unit cell has four atoms labeled as AU, AL, BU, and BL (see the top view in Fig. 1b), where “U” and “L” indicate the upper and lower sublayer, respectively. The (001) surface of SnS(Se) shares the same lattice structure with bismuthene, with S(Se) occupying A sites and Sn occupying B sites. We performed scanning tunneling microscope (STM) measurements to map the surface topography of the Bi sample grown on SnS(Se). The STM result in Fig. 1c indicates the high structural quality of bismuthene. The apparent height of bismuthene (including the thickness of bimuthene and the spacing between bismuthene and SnS(001) surface) is 8.0 Å, which can be seen from the height profile in Fig. 1d. The zoom-in STM images in Fig. 1e demonstrate the surface unit cell of bismuthene and SnS(Se). We found the in-plane lattice constants (a, b) are (4.5, 4.8 Å) for bismuthene, (4.1, 4.5 Å) for SnS surface, and (4.3, 4.6 Å) for SnSe surface, respectively.

Band structure of the 2D Weyl semimetal

Here we present the first-principles band structures and angle-resolved photoemission spectroscopy (ARPES) spectra to verify the existence of the 2D Weyl fermion states in the epitaxial bismuthene films. The calculated band structure of free-standing bismuthene is shown in Fig. 1f. The only prominent band feature near the Fermi level is a gapped Dirac cone located at a generic k point between \(\overline{\Gamma }\) and \({\overline{{{{{{{{\rm{X}}}}}}}}}}_{1}\). The gapped Dirac bands originated from pz orbitals of Bi atoms. The effective k ⋅ p model around the Dirac point can be written as \({{{{{{{{\mathcal{H}}}}}}}}}_{0}^{{{{{{{{\rm{Dirac}}}}}}}}}({{{{{{{\boldsymbol{k}}}}}}}})={\tau }_{\pm }({v}_{x}{k}_{x}{\sigma }_{x}+\Delta {k}_{x})+{v}_{y}{k}_{y}{\sigma }_{y}+{\tau }_{\pm }{\lambda }_{{{{{{{{\rm{SOC}}}}}}}}}{\sigma }_{z}{s}_{z}\), where (kx, ky) are measured from the Dirac nodes at ( ± k0, 0), σi (i = x,y,z) are Pauli matrices with respect to the basis of \(\{\left\vert {{{{{{{\rm{A}}}}}}}},{p}_{z}\right\rangle,\left\vert {{{{{{{\rm{B}}}}}}}},{p}_{z}\right\rangle \}\) (the pz orbitals at the two sublattices A and B), si (i = x,y,z) are the spin matrices, Δ describes tilting of the Weyl cone in the kx direction and is crucial for generating a non-zero Berry curvature dipole (BCD)24, λSOC is the effective spin–orbit coupling, τ± = ± 1 represents the chirality of the Dirac nodes located at ( ± k0, 0), and vx,y are Fermi velocity along kx and ky directions, respectively. vx = 3.17 × 105 m/s, Δ = 0.19 × 105 m/s, vy = 4.23 × 105 m/s, λSOC = 55 meV according to the first-principles results. The energy gap induced by spin-orbit coupling is Egap = 2λSOC = 0.11 eV. Every band is doubly degenerate to the spin degree of freedom since the lattice of bismuthene is centrosymmetric.

The calculated band structure of bismuthene on SnSe is plotted in Fig. 1g. The presence of SnSe substrate breaks space-inversion symmetry and causes spin splitting in the bands of bismuthene. Specifically, the electronegativity difference between Sn and Se induces an in-plane dipole electric field, leading to an average energy shift of A sites relative to B sites of bismuthene. In addition, the surface potential of SnSe generates an electric field perpendicular to the surface, creating a potential difference between the two Bi sublayers. The substrate perturbations can be effectively expressed as \({\lambda }_{{{{{{{{\rm{Dip}}}}}}}}}{\sigma }_{z}+{\lambda }_{V}{\sigma }_{y}{s}_{x}+{\tau }_{\pm }{\lambda }_{V}^{{\prime} }{\sigma }_{x}{s}_{y}\), where λDip describes the in-plane dipole field on the SnSe surface and \({\lambda }_{V}({\lambda }_{V}^{{\prime} })\) is the Rashba coupling caused by the electric field perpendicular to the substrate surface (see Supplementary Information for details). The spin splitting eliminates the SOC-induced energy gap, leading to the formation of a linear band crossing around the Fermi level. The 2D Weyl states generated through this method represent a fine-tuned critical state at the transition between two topologically distinct gapped phases. A detailed discussion on the formation mechanism of 2D Weyl states can be found in Supplementary Information-Section II. The states originating from the SnSe substrate are exclusively positioned either 0.9 eV below or 0.4 eV above the Fermi level. This unique arrangement can be attributed to the semiconductor nature of SnSe, characterized by a band gap of 1.3 eV. As a consequence, the 2D Weyl fermion states near the Fermi level primarily arise from Bi orbitals, confining them spatially within the Bi monolayer. This characteristic defines the Weyl fermion states as inherently 2D in nature.

The calculated 2D Weyl band structure is confirmed by our ARPES measurements. The ARPES results taken from the bismuthene/SnSe(S) sample are plotted in Fig. 2. The Fermi surface (Fig. 2a) contains two circular electron pockets in the direction of \({\overline{{{{{{{{\rm{X}}}}}}}}}}_{1}-\overline{\Gamma }-{\overline{{{{{{{{\rm{X}}}}}}}}}}_{1}\). (We note that a similar pair of electron pockets show up in the direction of \({\overline{{{{{{{{\rm{X}}}}}}}}}}_{2}-\overline{\Gamma }-{\overline{{{{{{{{\rm{X}}}}}}}}}}_{2}\) but with much lower intensity. This extra pair of pockets is due to the existence of Bi domains rotated by 90° in the MBE sample.) The ARPES spectra of bismuthene on SnSe along the lines of “cut1” and “cut2” are plotted in Fig. 2b, e. In the ARPES spectra, we found the band dispersion along “cut1” ("cut2”) as well as that from a rotated Bi domain. This can be better seen in Fig. 2c, f with overlays of calculated bands on top of the ARPES spectra. The magenta lines are bands along “cut1” ("cut2”), while the green lines are bands along a direction perpendicular to “cut1” ("cut2”). The sole presence of linear Weyl cones near the Fermi level indicates that the transport and optical properties of this system are exclusively determined by the low-energy Weyl fermion states. No apparent gap was found at the nodal point, as evidenced by the second derivative spectrum (Fig. 2h) and the map of energy distribution curves (Fig. 2i). We notice that one linear branch of the Weyl cone is much dimmer than the other in the spectrum along “cut1” due to the photoemission matrix element effects. Figure 2d and g show the calculated spectra with the inclusion of the photoemission matrix elements, which agrees well with the ARPES result. The bismuthene/SnSe sample is electron-doped as the nodal point lies 0.1 eV below the Fermi level. We also performed ARPES measurements on bismuthene grown on SnS, and the results are plotted in Fig. 2j–l. Compared with bismuthene/SnSe, the Fermi level of bimuthene/SnS lies right at the Weyl node. The shift of the Fermi level can be attributed to the electronegativity difference and the doping deference between SnSe and SnS substrates. So, bismuthene/SnS is a perfect 2D Weyl semimetal with charge neutrality. Considering the different surface conditions of SnSe and SnS, the observation of gapless Weyl cones in both sample configurations indicates the robustness of the 2D Weyl fermion states against weak perturbations.

Fig. 2: ARPES and first-principles band structure of Weyl fermion states in epitaxial bismuthene.
figure 2

a ARPES Fermi surface taken from bismuthene on SnSe. b ARPES spectra of bismuthene taken along the line of “cut1” marked in (a). c Overlay of calculated band structure on the ARPES spectrum along “cut1''. The magenta lines are bands along the direction of \(\overline{\Gamma }\)-\({\overline{{{{{{{{\rm{X}}}}}}}}}}_{1}\) while the green lines are bands in the direction perpendicular to \(\overline{\Gamma }\)-\({\overline{{{{{{{{\rm{X}}}}}}}}}}_{2}\). d Calculated band spectra with the inclusion of photoemission matrix elements. e–g Same as b–d but for bands of bismuthene on SnSe along the line of “cut2'', which passes across the Weyl node. h Second derivative of the ARPES spectrum in the red box in (e). i Energy distribution curves (EDC) from the ARPES spectrum inside the red box in (e). The blue dotted lines mark the maximum of each EDC. The red solid line plots the EDC taken at the momentum of the Weyl point. j–l Same as b–d, but for bands of bismuthene on SnS along “cut2''.

Spin texture of the 2D Weyl cone

The defining characteristic of Weyl fermion states is the spin polarization of the linear band dispersion. The dipole term (∝λDip) together with the SOC term cause the spin to be polarized in the z direction while the vertical-field terms (∝λV and \({\lambda }_{{{{{{{{\rm{V}}}}}}}}}^{{\prime} }\)) give rise to an in-plane spin polarization (see Supplementary Information). As a result, a canted spin texture is expected for the Weyl cone. The calculated spin-resolved band spectrum along “cut1” is plotted in Fig. 3b. To verify the spin texture of the Weyl fermion states, we performed spin-resolved ARPES measurements on bismuthene/SnSe. Spin-resolved momentum distribution curves (MDC) along “cut1” at E = −0.1 eV are shown in Fig. 3c. The blue and red lines are photoemission intensity recorded in the “spin-down” and “spin-up” channels for the corresponding spin component, respectively. The spin polarization can be calculated from the spin-ARPES spectrum according to the formula, \(P=\frac{1}{{S}_{{{{{{{{\rm{eff}}}}}}}}}}\frac{{I}_{+}-{I}_{-}}{{I}_{+}+{I}_{-}}\), where the effective Sherman function Seff = 0.275 for our spin detectors44,45,46. The spin polarization extracted from the results in Fig. 3c is shown in Fig. 3d. The observed spin polarization of 〈sy〉 is in good agreement with the theoretical result. The sx component is absent along the line of “cut1”, because the \(\overline{\Gamma }-{\overline{{{{{{{{\rm{X}}}}}}}}}}_{1}\) direction corresponds to a glide mirror of the lattice. We observed that the two valleys of Weyl fermion states exhibit opposite spin polarizations, a consequence of their partnership under time-reversal symmetry. The spin-integrated and spin-resolved ARPES spectra along “cut2” are plotted in Fig. 3e and f. The red and blue dots in Fig. 3f represent the “spin-up” and “spin-down” signals, respectively, recorded in the spin detector. This 2D “snapshot” of the sy component of the Weyl cone aligns with the calculated result presented in Fig. 3g, where the two linear branches of the Weyl cone display the same sign in the sy component. This characteristic is also evident in the spin-resolved MDC curves taken at E = −0.25 eV in Fig. 3h and the extracted spin polarization in Fig. 3i. In contrast, the two linear branches along “cut2” exhibit opposite signs in the sx component, as shown in the theoretical spin spectrum (see Fig. 3j). This spin behavior is confirmed by the experimental results of spin-resolved MDCs (Fig. 3k) and spin polarization (Fig. 3l). The in-plane spin texture of the Weyl fermion states is schematically summarized in Fig. 3m. We also measured the out-of-plane spin component. The magnitude of the observed spin polarization 〈sz〉 is suppressed due to the existence of rotated domains in the MBE samples. Nonetheless, the observed 〈sz〉 exhibits a similar characteristic to the 〈sy〉 component, consistent with the first-principles results (see the Supplementary Information for details). Our spin-ARPES results demonstrate the unique spin texture of Weyl fermion states. This distinctive spin texture stands as a defining feature that sets apart 2D Weyl semimetals from Dirac semimetals like graphene.

Fig. 3: Spin texture of 2D Weyl fermion states.
figure 3

a ARPES spectrum along “cut1” marked in Fig. 2a. b Calculated band structure along “cut1” with the inclusion of photoemission matrix elements. The bands are colored according to the calculated expectation value 〈sy〉 of each state. c Spin-resolved momentum distribution curves (MDC) taken at E = −0.1 eV along the line marked by the pink dashed arrows in (a) and (b). The blue and red curves are the photoelectron intensity recorded in the “spin-down” and “spin-up” channels of the spin detector, respectively. d The spin polarization extracted from the MDCs in (c). The shaded area indicates net spin polarization of 〈sy〉. The formula for error bars can be found in the Supplementary Information. e ARPES spectrum along “cut2” marked in Fig. 2a. f 2D spin-resolved ARPES map of 〈sy〉 taken along “cut2''. The blue and red dots represent the photoelectron signals recorded in the “spin-down” and “spin-up” channels of the spin detector, respectively. g Calculated electron band structure along “cut2'', which is colored according to the expectation value 〈sy〉 of each state. h Spin-resolved MDCs of 〈sy〉 taken at E = −0.25 eV along the line marked by the green dashed arrows in (e) and (g). i Spin polarization of 〈sy〉 extracted from the MDCs in (h). j–l Same as g–i, but for the spin component 〈sx〉. m In-plane spin texture of iso-energy contours at E = −0.25 eV. The red and blue circles schematically depict the iso-energy contours of two valleys at E = −0.25 eV. The black arrows indicate the in-plane spin orientation of the Weyl fermion states.

Bulk–boundary correspondence

To show the unique bulk–boundary correspondence of 2D Weyl semimetals, we calculated the bands of a semi-infinite bismuthene/SnSe heterostructure with an open boundary in the (010) direction. The result is plotted in Fig. 4a. The Fermi string edge band emanates from the bulk Weyl nodes as required by the band topology of the 2D Weyl cone. More precisely, the Fermi string band comprises two chiral edge bands, each corresponding to the critical phase of a specific spin sector. To illustrate this property, we utilized a tight-binding model incorporating sz conservation to simulate the edge electronic structure of the 2D Weyl semimetal. The edge spectrum is shown in Fig. 4b. In this case, the Weyl bulk bands comprise two independent spin sectors with \({s}_{z}=+ \frac{1}{2}\) and \(-\frac{1}{2}\), given that sz is a good quantum number. Each spin sector (colored in red and blue, respectively) corresponds to the critical phase of Haldane’s quantum anomalous Hall model47. Consequently, each sector possesses a chiral edge band that disperses from the bulk nodal point at one valley to the bulk band edge at the other valley. Interestingly, the two edges demonstrate distinct connection patterns between the edge states and the bulk band. The Fermi strings connect to the valence band at the top edge and to the conduction band at the bottom edge, in accordance with the requirement of charge conservation. In real materials with strong spin–orbit coupling like bismuthene, sz is not a good quantum number, and thus the two spin sectors are allowed to hybridize with each other. On the other hand, the connection of Fermi string edge bands to the Weyl nodes remains topologically protected, because the winding number of linear bands is unaffected by spin–orbit coupling. (It is worth noting that, in addition to the Fermi string bands, there are additional edge state bands located within the bulk band gap (in Fig. 4a). The existence of those additional edge state bands can be attributed to the fact that the 2D Weyl semimetal is at a critical point in connection to two topologically distinct insulator phases (see Supplementary Information-Section XI for a detailed discussion).) The Fermi string edge bands give rise to an enhanced local density of states (LDOS) at the edge, especially, in a narrow energy window around the energy of bulk Weyl nodes. The LDOS can be directly probed by the differential conductivity, dI/dV spectrum in STS experiments. An STM topography of bismuthene on SnSe is shown in Fig. 4c. The surface of SnSe, the interior of bismuthene, and the edge of bismuthene are marked by black, green, and magenta dots, respectively. We measured the averaged dI/dV spectrum at 4.6 K from the three regions (black, green, and magenta), and the result is plotted in Fig. 4d. A large gap of ~1.3 eV is observed in the dI/dV curve from the SnSe surface, indicating the Weyl fermion states are entirely confined within the Bi overlayer. Remarkably, the edge dI/dV spectrum (magenta) shows a notably enhanced LDOS compared to the bulk spectrum (green) in a narrow energy window around EW (the inset in Fig. 4d). This enhanced LDOS can also be seen in the dI/dV maps taken around a step edge of bismuthene (Fig. 4e). The edge appears brighter than both the SnSe surface and the interior of bismuthene only at Vbias = −12 and 17 meV, which correspond to energies close to EW. A similar edge dI/dV spectrum was also observed in bismuthene/SnS samples at different temperatures (see Supplementary Information). The results consistently affirm the presence of Fermi string edge states near the energy of Weyl nodes.

Fig. 4: Bulk–boundary correspondence in 2D Weyl semimetals.
figure 4

a The edge bands and projected bulk bands of a semi-infinite bismuthene film on SnSe with an open boundary in the (010) direction. The bands are weighted with the charge density near the edge. d0 = 3.89 Å is the relaxed interlayer spacing between bismuthene and SnSe substrate. b The connection of Fermi string edge bands to bulk Weyl nodes. The band structure was calculated by using a tight-binding model in which sz is conserved. c STM topography of bismuthene on SnSe. d Differential conductivity dI/dV spectra with the bias voltage aligned with the energy of Weyl nodes. The green curve is the averaged spectrum taken at the green grid inside the bismuthene patches shown in (c). The black curve is the averaged spectrum from the black grid on the surface of SnSe. The magenta curve is the averaged spectrum from the magenta points at the edges of the bismuthene patches. The inset shows the zoom-in dI/dV curves around the energy of Weyl nodes. e dI/dV map over the area marked by the blue box in c at bias voltage V= −305, −12, 17, and 1000 meV. f High-resolution dI/dV map at bias voltage V = +10 meV. g The corresponding Fourier transform of the dI/dV map in f, showing the quasiparticle interference (QPI) pattern. h The calculated quasiparticle interference pattern with spin-dependent scattering probability is considered.

The high-resolution dI/dV map unveils a unique plane wave-like pattern near the energy of Weyl nodes, as illustrated in Fig. 4f. The Fourier transform of the dI/dV map produces the quasiparticle interference pattern (Fig. 4g), featuring an oval contour at the center of the q space along with two smaller satellite contours. The separation between the satellite and the center of the oval, Δq = 0.42 Å−1, precisely matches the separation between the two Weyl nodes measured from ARPES results (Fig. 2a). Hence, the central oval contour in the QPI arises from intravalley scatterings, while the two satellites are induced by intervalley scatterings. The experimental QPI is consistent with the calculated result shown in Fig. 4h. The plane wave-like QPI, characterized by a single wavevector, reflects the nodal Fermi surface of the 2D Weyl semimetal.

The 2D Weyl semimetal establishes a solid-state realization of Weyl fermions in 2D space. The two spin-polarized valleys in the epitaxial 2D Weyl semimetal correspond to two Weyl fermions with opposite helicities. Our demonstration of the 2D Weyl semimetal and the Fermi string edge states paves the way for the exploration of fascinating topological quantum properties of Weyl fermions in reduced dimensionality.

Methods

Growth of bismuthene on SnS and SnSe

Bi was deposited on the cleaved surface of SnS and SnSe crystals in an MBE-ARPES-STM ultrahigh vacuum (UHV) system. The SnS and SnSe crystals are n-type doped with Br. The base pressure was lower than 2 × 10−10 mbar. High-purity Bi was evaporated from a standard Knudsen cell with a flux of 0.3 Å/min. The temperature of the substrate was kept at 50 °C during the growth. The substrate temperature is critical for growing a smooth Bi monolayer in the phosphorene structure.

Scanning tunneling microscopy measurement

An in-situ Aarhus-150 STM was used to characterize the surface topography and the lattice parameters of the α-Bi films. The topography was measured under room temperature, and the base pressure was lower than 2 × 10−10 mbar. The bias voltage and the tunneling current were set to be 1.5 V and 0.01 nA for the surface topography measurement, 5 mV and 0.15 nA for the zoom-in atom-resolved STM measurement. The dI/dV spectrum, dI/dV mapping, and the STM QPI were produced by an Omicron LT-Nanoprobe system at 4.6 K. The sample was transferred through an ultra-high vacuum (UHV) suitcase with pressure below 1 × 10−9 mbar. The tunneling current is set to be 200 pA during the dI/dV measurement.

Spin- and angle-resolved photoemission spectroscopy measurements

Spin and angle-resolved photoemission spectroscopy measurements were performed in a lab-based system coupled to the molecular beam epitaxy system, using a Scienta DA30L hemispherical analyzer with a base pressure of P < 5 × 10−11 mbar and a base temperature of T ~ 8 K. Samples were illuminated with linearly polarized light using an Oxide hν = 11 eV laser system. The light polarization was set perpendicular to the sample and the slits of the detector. For electronic dispersion measurements, a pass energy of 2 eV and 0.3 mm slit was used for a total energy resolution ~2.5 meV and momentum resolution ~0.01 Å−1. Dual VLEED ferrums that utilize exchange scattering are coupled to the electron analyzer and used to determine the spin sx, sy, and sz polarizations of the measured electrons. For spin-resolved measurements, a pass energy of 10 eV and a 1 mm × 2 mm spin aperture was used, yielding a total energy resolution ~50 meV and momentum resolution ~0.033 Å−1.

First-principles calculation

First-principles calculations with density function theory (DFT) were performed by using the Vienna ab Initio Simulation Package (VASP) package48. The Perdew–Burke–Ernzerhof (PBE)49 exchange-correlation functional was used. The experimental lattice parameter was applied for bismuthene; the lattice parameter of SnSe(S) was modified to match the covered Bi. The spin–orbit coupling (SOC) was included self-consistently in the calculations of electronic structures with a Monkhorst–Pack 11 × 11 × 1 k-point mesh. The vacuum thickness was >20 Å to ensure the separation of the slabs. Atomic relaxation was used until the residual forces were <0.01 eV/Å.

We constructed a tight-binding Hamiltonian for Bi/SnSe(S), where the tight-binding model matrix elements were calculated by projecting onto the Wannier orbitals50,51,52, which used the VASP2WANNIER90 interface53. The Bi p orbitals, Sn p orbitals, and Se p orbitals were used to construct the Wannier functions without performing the maximizing localization. The edge state electronic structure was calculated by Green’s function technique, which computes the spectral weight near the edge of a semi-infinite system. To simulate the photoemission matrix element effects in the ARPES spectra, we consider a non-trivial structure factor in ab initio calculations54,55. To simulate this effect, we construct a unitary matrix U(k):

$$U({{{{{{{\boldsymbol{k}}}}}}}})=\left(\begin{array}{ccc}{{{e}}}^{i{{{{{{{\boldsymbol{k}}}}}}}}\cdot {{{{{{{{\boldsymbol{r}}}}}}}}}_{1}}&\cdots \,&0\\ \vdots &\ddots &\vdots \\ 0&\cdots \,&{{{e}}}^{i{{{{{{{\boldsymbol{k}}}}}}}}\cdot {{{{{{{{\boldsymbol{r}}}}}}}}}_{n}}\end{array}\right),$$

where k is the momentum vector, ri is the real space coordinates of the ith atom in the original Bi/SnSe(S) unit cell. We simulate the unfolded band structures of Bi/SnSe by applying this unitary matrix to the tight-binding Hamiltonian, U(k)H(k)U(k)†. The details of the unfolding procedure are described in the work by W. Ku et al.56.

The quasiparticle interference pattern was calculated based on Green’s function method by using the spin-dependent scattering probability (SSP) method57,58, which can be written as

$${J}_{{s}}(q)=\frac{1}{2}\sum\limits_{k}\sum\limits_{i=0,1,2,3}{\rho }_{i}(k){\rho }_{i}(k+q),$$

where \({\rho }_{0}(k)={{{{{{{\rm{Tr}}}}}}}}[G(k)]\) is the total spectral density and \({\rho }_{i}(k)={{{{{{{\rm{Tr}}}}}}}}[{\sigma }_{i}G(k)]\) is the spin density, in which G(k) = [ω+iη−H(k)]−1 is Green’s function of the system and σi=1,2,3 are the Pauli matrices for spin.