Giant valley-Zeeman coupling in the surface layer of an intercalated transition metal dichalcogenide

Spin–valley locking is ubiquitous among transition metal dichalcogenides with local or global inversion asymmetry, in turn stabilizing properties such as Ising superconductivity, and opening routes towards ‘valleytronics’. The underlying valley–spin splitting is set by spin–orbit coupling but can be tuned via the application of external magnetic fields or through proximity coupling. However, only modest changes have been realized to date. Here, we investigate the electronic structure of the V-intercalated transition metal dichalcogenide V1/3NbS2 using microscopic-area spatially resolved and angle-resolved photoemission spectroscopy. Our measurements and corresponding density functional theory calculations reveal that the bulk magnetic order induces a giant valley-selective Ising coupling exceeding 50 meV in the surface NbS2 layer, equivalent to application of a ~250 T magnetic field. This energy scale is of comparable magnitude to the intrinsic spin–orbit splittings, and indicates how coupling of local magnetic moments to itinerant states of a transition metal dichalcogenide monolayer provides a powerful route to controlling their valley–spin splittings. The authors study the electronic structure of the intercalated transition metal dichalcogenide V1/3NbS2, showing that its bulk magnetism can lead to a strong tunability of spin–valley locked states at its surface.

Spin-valley locking is ubiquitous among transition metal dichalcogenides with local or global inversion asymmetry, in turn stabilizing properties such as Ising superconductivity, and opening routes towards 'valleytronics'. The underlying valley-spin splitting is set by spin-orbit coupling but can be tuned via the application of external magnetic fields or through proximity coupling. However, only modest changes have been realized to date. Here, we investigate the electronic structure of the V-intercalated transition metal dichalcogenide V 1/3 NbS 2 using microscopic-area spatially resolved and angle-resolved photoemission spectroscopy. Our measurements and corresponding density functional theory calculations reveal that the bulk magnetic order induces a giant valley-selective Ising coupling exceeding 50 meV in the surface NbS 2 layer, equivalent to application of a ~250 T magnetic field. This energy scale is of comparable magnitude to the intrinsic spin-orbit splittings, and indicates how coupling of local magnetic moments to itinerant states of a transition metal dichalcogenide monolayer provides a powerful route to controlling their valley-spin splittings. NbS 2 is one of the most intriguing members of the transition metal dichalcogenide (TMD) family. It is the only TMD to exhibit superconductivity that does not develop from a charge density wave state 1,2 , although it is thought to lie close to not only charge-ordering instabilities but also magnetic instabilities 3 . Intercalating transition metal ions (for example, M = V, Cr, Co, Fe) into the van der Waals gap between neighbouring NbS 2 layers can lead to the development of long-range magnetic orders that host a variety of non-trivial spin textures [4][5][6][7][8] . Here, we consider the critical composition of M 1/3 NbS 2 with M = V. The intercalated V atoms form an ordered structure, occupying distinct sites above and below each NbS 2 layer, V 1 and V 2 , respectively (Fig. 1a), making the crystal non-centrosymmetric. Below a critical temperature, T N ≈ 50 K, the V spins order. Detailed neutron diffraction measurements 9 have revealed two coexisting components to the magnetism: an A-type antiferromagnetic order, confined within the ab plane and consisting of ferromagnetic layers of V atoms coupled antiferromagnetically between layers; and a second component with moments aligned along the c axis and consisting of a one-up two-down longitudinal density wave-like modulation of the spins. Overall this forms a canted antiferromagnetic state (Fig. 1b), which is reported to lead to a net uncompensated ferromagnetic moment 9,10 . Article https://doi.org/10.1038/s41563-022-01459-z origins of the two-dimensional electron gas at the SrTiO 3 /LaAlO 3 interface 14 and for extreme self-doping at surfaces of, for example, some cuprates 15 and delafossite oxides 16,17 . We investigate this experimentally in Fig. 1c-h.
Our samples were cleaved in situ before commencing photoemission measurements (Methods). Given the covalent nature of the NbS 2 block, the NbS 2 -layer and V-layer surface terminations discussed above and shown in Fig. 1a represent the natural cleavage planes of the crystal. Both terminations can be expected to be formed with equal probability and distributed randomly across the sample surface (with mixed and/or disordered terminations also possible 12 ). We investigate this with spatially resolved core-level spectroscopy in Fig. 1c-e. With the photon energies used here, our measurements are extremely surface sensitive, with a probing depth limited to the top few atomic layers 18 . The relative ratio of the integrated spectral weight of, for example, Nb 4p and V 3p core levels therefore provides a good measure of the topmost atomic species at the surface. As shown in Fig. 1c, this ratio shows a strong spatial variation on a length scale of ~20-30 μm. We assign regions of higher (lower) intensity to NbS 2 (V) surface terminations.
Core-level spectra extracted from representative regions of our spatial mapping data (rectangles in Fig. 1c) are shown in Fig. 1d,e. As well as their relative intensity variations, our measurements As well as inducing magnetic order, the incorporation of nominally V 3+ ions can be expected to lead to substantial charge transfer into the NbS 2 layer as compared to pristine NbS 2 (refs. [11][12][13]. Indeed, in a simple ionic picture, the V-intercalated compound can be viewed as an alternating stack of positive [V 1/3 ] + and negative [NbS 2 ] − charged layers, with the NbS 2 layer thus electron doped relative to pristine NbS 2 . Consistent with this, our density functional theory (DFT) calculations of the charge density difference between V 1/3 NbS 2 and NbS 2 ( Fig. 1a; also Methods and Supplementary Fig. 1) indicate the expected charge transfer from the V to the NbS 2 layers, with the V 3+ cations acting as an ionic linkage between the otherwise weakly bonded NbS 2 layers.
Further modifications to the charge transfer can also be seen at the material surface. Charge transfer into a NbS 2 layer in the bulk from a V layer above can no longer occur at a NbS 2 -terminated surface. The surface NbS 2 layer would thus be expected to become hole doped as compared to its bulk counterpart. Conversely, a V-terminated surface would be expected to become electron doped as compared with the bulk. This simple picture is supported by our calculations of the surface electronic structure, shown in Supplementary Fig. 2, where the formation of new surface bands, which result from compensating the polar surface charge, can be viewed as a form of electronic reconstruction. This is similar to that discussed for the so-called 'polar catastrophe'   a, Crystal structure of V 1/3 NbS 2 , displaying the bulk and distinct surfacetermination-dependent charge transfers as compared to pristine NbS 2 , calculated using DFT. b, Bulk magnetic structure from the literature 9 , displaying the relative orientation of local magnetic moments on the V sites in two side-byside magnetic unit cells. c, X-ray photoelectron spectroscopy (XPS) spatial map displaying the ratio of the integrated spectral weight of Nb 4p to V 3p core levels. High (yellow) and low (blue) intensity regions correspond to the NbS 2 and V surface terminations, respectively. Max, maximum; Min, minimum. Article https://doi.org/10.1038/s41563-022-01459-z indicate a shift of the core-level peaks to a lower binding energy for the NbS 2 -terminated surface with respect to the V-terminated surface. This is a direct experimental signature of the additional surface charge transfer discussed above, with hole doping at the NbS 2 -terminated surface and electron doping at the V-terminated surface, as discussed further in Supplementary Note 1 and Supplementary Fig. 3. Interpretation of the S 2p core levels is more challenging due to the presence of well-defined surface core-level components. Nonetheless, mapping of a pronounced component that develops at low binding energies for the NbS 2 -terminated surface yields a spatial distribution that is in excellent agreement with that obtained by considering the Nb 4p/V 3p core-level ratios in Fig. 1c (Supplementary Fig. 4). Together, our core-level measurements thus allow ready identification of distinct surface terminations and point to a dominant electronic reconstruction as the route to compensate the corresponding polar surface charge. The resulting surface charge transfer leads, in turn, to marked differences in the surface electronic structure. Figure 1g,h shows the measured electronic structure integrated over the same spatial regions as for the core-level spectra discussed above. For the V termination, we observe somewhat diffuse but dispersive features, which we attribute to NbS 2 -derived bulk states, consistent with our calculations of such dispersive states in the bulk electronic structure of V 1/3 NbS 2 , shown in Supplementary Fig. 1. Additionally, we find a rather non-dispersive spectral weight, which likely results from the more localized V-derived states ( Supplementary Fig. 2).
The NbS 2 termination, on the other hand, hosts a much richer electronic structure. In particular, we find several sharp dispersive states in the vicinity of the Fermi level. These appear to be hole-doped analogues of the dispersive states observed for the V termination. From photon-energy-dependent measurements (Supplementary Figs. 5 and 6), we find the new states arising on the NbS 2 termination are strictly two-dimensional. Integrating the spectral weight over these dispersive states (dashed region shown in Fig. 1g,h), we find a spatial intensity distribution (Fig. 1f) that closely resembles that derived from the core-level spectra (Fig. 1c). We thus attribute these new dispersive states as surface states of the NbS 2 -terminated surface of V 1/3 NbS 2 .
We show in Fig. 2a,b higher-resolution measurements of the low-energy electronic structure of the NbS 2 -terminated surface, measured below the bulk magnetic ordering temperature. The surface bands are split off from the underlying bulk V 1/3 NbS 2 bands due to the self-doping at this polar surface, becoming located in projected bandgaps of the bulk electronic structure (also Supplementary Fig. 2). They are thus well localized at the surface and, in a simple picture, can be viewed as deriving from a monolayer-like surface NbS 2 layer. For the broken inversion symmetry inherent to such a monolayer, spin-orbit Article https://doi.org/10.1038/s41563-022-01459-z coupling lifts the spin degeneracy, leading to a locking of quasiparticle spin to a valley pseudospin as is evident in Fig. 2d [19][20][21][22][23][24] .
Signatures of the electronic structure that would be expected for such a spin-valley locked NbS 2 monolayer (Fig. 2c,d) can be identified in our experimental measurements. In particular, we observe a large hexagonal pocket around the Brillouin zone centre (labelled β in Fig. 2a,c) and two split-off states (which we denote α 1,2 ), which disperse upward towards the K point of the NbS 2 Brillouin zone, forming a pair of hole-like barrels centered at this zone corner point. The Fermi pockets observed experimentally (Fig. 2b) are qualitatively similar to those in the calculations of monolayer NbS 2 (Fig. 2d), although they are substantially smaller. This size mismatch can largely be understood as a consequence of the charge transfer discussed above: while the [NbS 2 ] − layer in bulk V 1/3 NbS 2 is formally expected to have 2 electrons filled in the Nb d-orbital manifold (a d 2 charge count), an additional nominal hole doping of 0.5 holes per formula unit at the NbS 2 -terminated surface should render the surface layer in a d 1.5 configuration, electron doped by 0.5 electrons per formula unit as compared to a pristine NbS 2 monolayer. Indeed, from a Luttinger analysis of the experimentally measured Fermi surfaces, we find an electron count of the NbS 2 surface layer of 1.48 ± 0.02, in agreement within experimental error of that predicted from our simple charge transfer arguments.
There are, however, a number of important differences in the measured electronic structure with respect to a simple doped NbS 2 monolayer. The first is an evident doubling of all of the states, with a replica of the large hexagonal pocket also found centered at each Brillouin zone corner (β bf ), and additional smaller pockets found at the Brillouin zone centre ( bf 1,2 ). This is a signature of a new periodicity induced by the subsurface V atoms. As shown in Fig. 2e, a V atom is located directly below one in every three surface Nb atoms, leading to a new periodic potential with a √3 × √3R30 ∘ periodicity as compared to a pristine NbS 2 monolayer. 4,9 This is evident in our low-energy electron diffraction measurements in Fig. 2f, generating a reduced and rotated Brillouin zone shown by the dashed line in Fig. 2b.
The electronic states from the larger parent NbS 2 zone will become back-folded about the new Brillouin zone boundary, leading to the additional pockets discussed above. The opening of hybridization gaps can be expected at the new Brillouin zone boundary M points, and these are evident in our measured dispersions in Fig. 2a (also Supplementary  Fig. 8). Even considering this back-folding, however, there remains an intriguing discrepancy with the expected electronic structure for a NbS 2 monolayer. Specifically, the α 1,2 bands along the Γ-K line exhibit substantial splittings, while along Γ-K′, the band splitting between these states appears to be absent (additional data in Supplementary Fig. 9).  The colours indicate temperature as in e and f. These values indicate a band splitting that is constant for both directions at high temperature, but that diverges below the magnetic ordering temperature, as evident by comparison to the field-cooled d.c. magnetic susceptibility (χ) shown in i, measured in an applied field of H = 330 Oe (H ⊥ c axis). The functional form of the susceptibility data here is a result of a net uncompensated ferromagnetic moment for the bulk magnetic structure of V 1/3 NbS 2 (ref. 9 ). Error bars in h reflect an approximate estimate of the uncertainty in extracting the underlying peak positions from the experimental measurements, incorporating statistical errors in peak fitting, systematic errors and experimental resolution. The greyscale colour bands indicate the measured ARPES intensity.
From the spatial mapping data shown in Fig. 3a-c, we find that there are distinct regions for which the splitting of the Nb-derived conduction band (α 1,2 ) appears on the positive side (Fig. 3b, orange colouring in the inset of Fig. 3a) or negative side (Fig. 3c, purple colouring in Fig. 3a) of our measured momentum. This suggests that the momentum-dependent modulation of the band splittings observed here may be of magnetic origin, with the signs of the splitting reversed for different magnetic domains of the bulk antiferromagnetic order. To confirm this interpretation, we show in Fig. 3d-g measurements of the evolution of the band splitting within a single domain as a function of temperature. At low temperature (Fig. 3g), the clear asymmetry in the splitting of the α 1,2 bands for positive and negative momentum is evident in the Fermi level momentum-distribution curves (MDCs) shown in Fig. 3e,f: a clear two-peak structure is observed near the K point, but a single sharp (albeit slightly asymmetric) peak is observed near the K′ point. By contrast, the high-temperature Fermi level MDCs extracted near K and K′ are mirror-symmetric copies of each other. The observed peaks are broad and rather flat topped, indicating the presence of two components and hence a remnant band splitting, but one that is the same for both positive and negative momenta (also Supplementary Fig. 7).
The MDC peak splitting evolves smoothly between the low-and high-temperature limits up to temperatures of ~50 K, and remains rather static above this. This is confirmed from the Fermi level momentum splitting, Δk F , extracted from fits to our MDCs (Fig. 3h): the high-temperature band splitting is found to be identical within experimental error for the positive and negative momentum sides, while it diverges below a critical temperature of ~53 K, growing on the positive momentum side and collapsing to near zero on the negative momentum side. Magnetic susceptibility data from our samples (Fig. 3i) indicate that the onset of the change in band splitting occurs precisely at the magnetic transition temperature.
Our temperature-dependent measurements thus indicate that the unusual momentum-dependent band splittings observed at low temperature (Fig. 2a,b) result from exchange coupling of the itinerant NbS 2 -derived surface states to local magnetic moments on the underlying V sites. To confirm this, we extract the corresponding energetic shifts of the NbS 2 surface band from our temperature-dependent angle-resolved photoemission spectroscopy (ARPES) measurements (Supplementary Figs. 10 and 11). The resulting temperature-dependent changes in band position reproduce the measured magnetic susceptibility well, indicating that they reflect an order parameter of the magnetization, and thus can be assigned as an exchange splitting.
To explore the origin of this exchange coupling further, we show in Fig. 4a calculations of the real-space projected spin density isosurfaces for bulk V 1/3 NbS 2 , performed for a simple antiferromagnetic configuration. We find a spatially alternating spin density, predominantly localized around the V ions but with partial fragmentation within the intermediate NbS 2 layers. The shape of the spin density distribution around the Nb sites is strikingly similar to that of a d z 2 orbital, which, as shown in Supplementary Fig. 1, is the primary contributor to the only band crossing the Fermi level in bulk V 1/3 NbS 2 . Moreover, the spin density around the Nb sites carry a sign opposite to that of their adjacent V sites, together directly indicating an exchange interaction between V magnetization and the spins of the itinerant Nb electrons.
To further confirm this picture, we have performed bulk-sensitive resonant photoemission measurements at the V L 2,3 edge, allowing us to selectively enhance the spectral weight derived from the V d orbitals. Our measurements ( Supplementary Fig. 12) indicate that the V states contribute a well-defined peak in the density of states at a binding energy of ~1 eV, with a negligible weight persisting to the Fermi level. This is unlike the case for the Cr-based and Co-based sister Article https://doi.org/10.1038/s41563-022-01459-z compounds 13,25,26 , where more pronounced band hybridization in the vicinity of the Fermi level has been proposed (in Supplementary Fig.  13 we compare surface-termination-dependent ARPES data of V 1/3 NbS 2 and Cr 1/3 NbS 2 ). Instead, in V 1/3 NbS 2 , our measurements demonstrate that there is minimal band hybridization between the V-derived and Nb-derived states at the Fermi level. While further studies are required to definitively identify the magnetic exchange pathways, these findings, together with the large V-V real-space distance, are suggestive of a carrier-mediated Ruderman-Kittel-Kasuya-Yosida-like exchange interaction 27,28 underpinning the magnetic ordering of V 1/3 NbS 2 . Crucially, and irrespective of the precise hierarchy of exchange couplings here, our observations of pronounced temperature-dependent shifts in the surface band structure can unambiguously be identified as the result of an induced moment in the Nb-derived itinerant states. To understand the critical role of this coupling on the electronic structure, we model the NbS 2 surface layer as a NbS 2 monolayer with time-reversal symmetry breaking from proximity coupling to a magnetic layer beneath. Specifically, we consider hopping terms for the NbS 2 monolayer, and add the coupling to V as a splitting on the spin subspace, to construct a tight-binding Hamiltonian for the Nb d xy , d x 2 −y 2 and d z 2 orbitals and S p x , p y and p z orbitals with the form where ℋ 0 is the Hamiltonian of an unperturbed NbS 2 monolayer and ℋ ex describes the spin splitting that results from the magnetic exchange where J i is the effective exchange coupling between the spin operator S i on the NbS 2 site i and the spin operator S V on the adjacent V site. This can be decomposed into two contributions, ℋ ∥ ex and ℋ ⟂ ex , arising from two distinct coupling processes J ∥ and J ⊥ between the V d orbitals and the in-plane {d xy , d x 2 −y 2 } and out-of-plane d z 2 NbS 2 orbitals, as shown in Fig. 4b,c, respectively.
Considering an exchange field directed along [001] (Supplementary Figs. 14 and 15 for a discussion of additional field configurations and a comparison of first-principles and experimentally derived exchange parameters), we examine the effect each coupling mechanism has on the spin splitting as a perturbation from the case with only spin-orbit coupling included. From the orbital projections of the calculation with only spin-orbit coupling shown in Fig. 2c, one can anticipate that at Γ, the exchange splitting will be dominated by the ℋ ⟂ ex term, and at K/K′, it will be dominated by ℋ ∥ ex . Consistent with this, we find in Fig. 4d,e that the ℋ ∥ ex term leads to marked modifications in the electronic structure close to the K and K′ points as compared to the calculations with only spin-orbit coupling, while the electronic structure is almost unchanged close to Γ. By contrast, for the ℋ ⟂ ex term, exchange splittings are observed close to the Brillouin zone centre, while the valley-spin splitting around the K and K′ points is little changed from that shown in Fig. 2.
We include both of these exchange contributions in our model, as well as the band folding arising from the subsurface V-induced superlattice potential, which is visible in our ARPES measurements, as discussed above, but with reduced spectral intensity. Through this, we find a Fermi surface that is in excellent agreement with the one we observe experimentally (Fig. 4f). This includes both the valleydependent spin splitting at the K and K′ points, and also the azimuthalangle-dependent splittings of the bf 1,2 Fermi pockets back-folded to the Brillouin zone centre.
From the Fermi surface shown in Fig. 4f, it is clear that at and below the Fermi level, the changes in electronic structure we observe are dominated by exchange couplings between V d orbitals and the NbS 2 planar orbitals (Fig. 4b,d), with the spin splittings induced by exchange and spin-orbit coupling acting to enhance each other at the K point, while opposing each other at K′ (Fig. 4d). This is because the exchange interaction is valley independent, while the spin-orbit coupling splitting is of equal magnitude but opposite sign for each valley. Therefore, the total valley splitting here can be captured by the simple expression where τ = ±1 is the valley index and Δ SO is the spin-orbit splitting.
We thus attribute the change in surface electronic structure observed here upon cooling through the magnetic ordering temperature as being due to a giant valley-selective Ising coupling (a so-called 'valley-Zeeman coupling' [29][30][31][32]. From fits of the band splittings from our measured ARPES dispersions ( Supplementary Fig. 11), we estimate a total exchange splitting, Δ ex = 52 ± 7 meV. This is of comparable magnitude to the full intrinsic spin-orbit splittings observed in the normal state, which we estimate as Δ SO = 59 ± 4 meV. As a result, we find that a near spin degeneracy is recovered for one valley at low temperature, while a total spin splitting in excess of 110 meV is obtained for the other. Typical g-factors for TMDs lead to a valley-Zeeman splitting in an externally applied field of only ~0.2 meV T -1 (ref. 29 ). The splittings realized here would therefore require application of a magnetic field exceeding 250 T. Even if it were possible to apply such a large external field, that in itself would act to decouple the electron spin from its orbital degree of freedom, while the giant orbitally driven Ising exchange coupling observed here remains accessible through the magnetic proximity effect.
Our measurements of V 1/3 NbS 2 thus indicate how exchange coupling between local moments and itinerant states provides a critical route to enhanced control over valley-spin splittings in TMDs. These can be readily tuned by modest changes in temperature and provide substantially larger responses compared to other proximity-coupling schemes used to date [33][34][35] . Key to this giant effect is the strong coupling between the magnetic and itinerant layers (Fig. 4a). A large and rich family of intercalated TMDs exists, hosting distinct TMD layers, magnetic orders and critical temperatures. Furthermore, evidence exists of more substantial hybridization between the intercalated metal and TMD states for some such compounds than was observed here for V 1/3 NbS 2 (refs. 4,25,26 ). This raises the tantalizing prospect of gaining additional control over spin splittings in the TMD layer, and coupling the effects observed here with other collective states.

Online content
Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41563-022-01459-z.

Sample preparation and characterization
Single-crystal samples of V 1/3 NbS 2 were grown using the chemical vapour transport method with iodine as a transport agent. An evacuated quartz ampoule was slowly heated until a temperature gradient of 950 °C to 850 °C was applied across it, with polycrystalline V 1/3 NbS 2 situated in the 950 °C region of the ampoule acting as a seed for the growth. After 500 hours, the ampoule was slowly cooled. Temperature-dependent d.c. magnetization measurements were performed on single-crystal samples mounted in General Electric varnish on a non-magnetic straw in a Quantum Design Magnetic Property Measurement System (MPMS). The single crystals were aligned so that the applied field was perpendicular to the c axis of the crystals. The crystal was measured over a temperature range of 5 to 80 K in a field of 330 Oe.

Microscopic-area spatially resolved ARPES
Spatially resolved ARPES and XPS measurements were performed at the Bloch beamline of the MAX IV Laboratory using linearly polarized light with photon energies between 25 and 200 eV, and with a probing spot size of 10 × 15 μm 2 . The samples were mounted on a conventional six-axis manipulator that allowed cooling to 20 K. The samples were cleaved in situ at the base temperature and measured using a Scienta DA30 electron analyser.

Calculations
Electronic structure calculations were performed within DFT using the Perdew-Burke-Ernzerhof exchange-correlation functional 36 , as implemented in the Vienna Ab-initio Simulation Package programme 37,38 . Relativistic effects, including spin-orbit coupling, were fully included. For the bulk calculations, we considered a √3 × √3 × 1 supercell containing six formula units of NbS 2 and two intercalated V atoms antiferromagnetically ordered along the [001] axis, as shown in Fig. 4a. The corresponding Brillouin zone was sampled by a 15 × 15 × 10 k mesh. An additional on-site Hubbard term (U) with an effective value of 2 eV was added to the V 3d orbitals to reproduce the experimentally observed alignment of V bands with respect to the valence continuum of NbS 2 .
To compute the excess and depletion charge distribution due to the V intercalation, we first calculated the charge density for a five-layer slab of V 1/3 NbS 2 stacked along the crystalline c axis with a vacuum thickness of 15 Å and two different terminations, shown in Fig. 1a. The Brillouin zone sampling was done using a 20 × 20 × 1 k mesh. We next calculated the charge densities from the NbS 2 and V centres individually by removing V and NbS 2 layers in the same slab, respectively. Finally, we subtracted these individual contributions from the V 1/3 NbS 2 charge density, treating the remaining as excess (wherever it was positive) or depletion (wherever it was negative) charge density.
To model valley-Zeeman coupling, we carried out a separate DFT calculation for a monolayer of NbS 2 within the same level of theory as that of the slab calculations. From this, we constructed a 22-band tight-binding model using Wannier functions 39 with Nb d and S p orbitals as the projection centres. This model was then reduced to a four-band model in the basis of {d xy , d x 2 −y 2 } and d z 2. An Ising-like exchange coupling term ℋ ex was further included to account for the magnetic interactions, as discussed further in the main text. The exchange coupling constants were first deduced from the nearest neighbour Nb-V hopping parameters interpolated from a first-principles calculation performed on bulk V 1/3 NbS 2 using maximally localized Wannier functions. They were then slightly modified empirically to best match the experimental data. This led to only small quantitative changes in the electronic structure, as shown in Supplementary Fig. 15. The best agreement with the experiment was for J ∥ = 50 meV and J ⊥ = 60 meV.