Unidirectional Kondo scattering in layered NbS2

Crystalline defects can modify quantum interactions in solids, causing unintuitive, even favourable, properties such as quantum Hall effect or superconducting vortex pinning. Here we present another example of this notion - an unexpected unidirectional Kondo scattering in single crystals of 2H-NbS2. This manifests as a pronounced low-temperature enhancement in the out-of-plane resistivity and thermopower below 40 K, hidden for the in-plane charge transport. The anomaly can be suppressed by the c-axis-oriented magnetic field, but is unaffected by field applied along the planes. The magnetic moments originate from layers of 1T-NbS2, which inevitably form during the growth, undergoing a charge-density-wave reconstruction with each superlattice cell (David-star-shaped cluster of Nb atoms) hosting a localised spin. Our results demonstrate the unique and highly anisotropic response of a spontaneously formed Kondo lattice heterostructure, intercalated in a layered conductor.


Introduction
Layered van der Waals materials, such as transition metal dichalcogenides (TMDs), have attracted major interest thanks to their rich variety of ground states and the possibility of their exfoliation down to an atomically thin level, which remarkably modifies their electronic properties 1,2 . Recent observations of intriguing physics in artificiallyassembled heterostructures highlight the importance of interlayer interactions. Examples include the outstanding stability of interlayer excitons in semiconducting TMDs 3 , and strongly correlated states in twisted bilayer systems 4 .
Relevant aspects of the inter-plane coupling can be deduced by probing out-of-plane charge transport, even in bulk materials 5 . However, enforcing the current flow strictly along the c axis can be rather challenging due to the crystals' common flake-like appearance and their propensity for delamination. Such a pitfall can distort measurement results by orders of magnitude, as demonstrated in our recent study of microstructured samples of 1T-TaS2 with a well-defined current flow 6 . This observation motivates a careful re-examination of the out-of-plane charge transport properties in this class of materials by adopting the latest state-of-the-art for quantum matter microfabrication 7 .
Here we present data on the out-of-plane electrical resistivity of bulk monocrystalline 2H-NbS2. This material is one of the three known structural variants of layered NbS2. The two other polytypes are 3R and 1T, the latter occurring only in atomically thin form 8,9 . As illustrated in Figure 1, 1T-NbS2 consists of corner-sharing octahedral NbS6 cells. Layers of 2H-and 3R-NbS2 both contain NbS6 units of trigonal prismatic geometry, but exhibit different stacking configurations. The 1T polytype has been attracting interest recently as a candidate for realising a two-dimensional magnetic system 11,12 . 2H-NbS2 has been actively featured in the literature due to a superconductivity below 6 K, proposed to have a multiband character. It also does not show any charge density wave (CDW) order, which is uncommon for metallic TMDs [13][14][15][16][17][18] . Another distinguishing feature of 2H-NbS2 is its non-trivial synthesis procedure. This polytype is thermodynamically stable in a relatively narrow range of 3 temperatures and reactant stoichiometries [19][20][21] . Crystals formed during high-temperature growth must be rapidly quenched in order to capture 2H-NbS2 in a metastable room-temperature state. However, X-ray diffraction studies have shown that the resultant material has up to 18% of pairs of neighbouring layers stacked in a 3R-like manner 22,23 . Additionally, diffuse X-ray scattering experiments 23 revealed weak traces of the √13 × √13 CDW reconstruction, which appears as a triangular superlattice of David-star-shaped clusters defined by 13 Nb atoms 24 .
Such a reconstruction is not expected for pure 2H-NbS2 or 3R-NbS2. Earlier theoretical investigations have predicted 1T-NbS2 to be particularly prone to developing such a CDW order 11,12 . One can therefore conclude that single crystals of 2H-NbS2 contain rare, atomically thin inclusions of the 1T polytype.
Our study of 2H-NbS2 revealed a remarkably strong low-temperature anomaly in the compound's out-of-plane resistivity (ρc), manifesting as a minimum at around 40 K, followed by a pronounced upturn upon further cooling.
The feature is simultaneously invisible in the in-plane resistivity (ρab), and shows a highly anisotropic response to magnetic field. Neither 2H-NbSe2 nor 3R-NbS2 display such an anomaly, implying that the phenomenon is linked to the structural defects specific to 2H-NbS2. 1T-NbS2, layers of which are one of such defects, were predicted to form a lattice of unpaired localised spins located at the centre of each David-star CDW superlattice cluster 11,12 . We argue that planes of magnetic moments, hosted by the inclusions of 1T-NbS2, cause a Kondo effect observable only when the current flows across these planes.

Results
Optimisation of sample geometry with focused ion beam (FIB) micro-milling greatly improves charge transport study precision 6,7 . Using this approach, we shaped single crystals into samples with well-defined few micron thick and wide current channels, oriented along the two principal directions: normal and parallel to the atomic planes ( Figure 2a shows a sample of 2H-NbS2 produced this way). Such a design allowed simultaneous measurement of both ρab and ρc via the four-point technique. Probing ρc on two segments of different surface-to-volume resulted in mutually consistent values, allowing us to ensure that our results were not distorted by the presence of surfacerelated effects or macroscopic defects. 4 Figure 2b shows the plots of ρab and ρc of 2H-NbS2 against temperature (T), as well as their ratio in the inset. Note that in contrast to the earlier study which reports an anisotropy of the order of 1000 (Ref. 25), our measured value was as low as 10 at room temperature, monotonically increasing to 180 on cooling. As it was shown for the case of 1T-TaS2 using finite element simulations (Ref. 6), such an overestimate by the older study could be a result of a sub-optimal measurement geometry and an incorrect prior assumption that the anisotropy is very large. While ρab has a conventional metallic temperature dependence, ρc is also metallic, but shows a few noteworthy features.
First, the residual out-of-plane resistivity is very high, presumably due to a significant concentration of static defects. Second, ρc approaches saturation in the high-temperature region. This flattening of resistivity may be attributed to the mean free paths decreasing to the point of becoming comparable to the interlayer separation, a concept known as the Mott-Ioffe-Regel limit 26 . Third, at low temperatures, ρc displays a minimum at around 40 K, with a major upturn at lower temperatures. No corresponding feature exists in ρab (in agreement with previous results 27 ). All studied samples of 2H-NbS2 showed qualitatively identical behaviour, with slight differences in the absolute values of resistivity-related to slight impurity content variations-and temperatures of the minimum distributed in the 30 K -40 K range.
We compared ρc of 2H-NbS2 to that of the isostructural and isovalent compound 2H-NbSe2 (dashed line in Figure   2b). The latter material did not exhibit a similar low-temperature anomaly. Based on the nominal lattice parameters, density functional theory calculations predict that the two compounds will have nearly identical electronic band structures (see Supplementary Note 1 and Supplementary Figures 1,2). We therefore conclude that the upturn of ρc of 2H-NbS2 is not intrinsic to the nominal structure of the compound, but is caused by crystalline lattice defects.
The Seebeck coefficient (S) is a useful quantity for sensitively probing energy landscape variations near Fermi level. The open-circuit voltage generated by the thermal gradient is unaffected by the presence of static defects such as vacancies or non-magnetic stacking faults. On the other hand, the Seebeck coefficient is a function of the energy dependence of the conduction electron scattering rate. This is then strongly affected by the occurrence of resonance peaks in the density of states close to the chemical potential. Seebeck coefficient of 2H-NbS2 revealed a prominent peak at approximately 15 K, appearing only for the out-of-plane thermal gradient (Figure 2c). The 5 phonon drag phenomenon produces a similar feature in the temperature dependence. However, it manifests only in conductors with long phonon and charge carrier mean free paths-such as semimetals or extremely pure metals-where momentum-conserving scattering is dominant 28 . This is highly unlikely for 2H-NbS2, as its high content of static defects should clearly favour momentum-relaxing scattering. Furthermore, absence of the corresponding Sab peak rules out the phonon drag from the possible Sc anomaly origins. An alternative interpretation of the peak, the Kondo effect, will be discussed further below.
The out-of-plane resistivity anomaly of 2H-NbS2 demonstrated a particularly curious response to magnetic fields.
As can be seen in Figure 3a, the transverse and longitudinal out-of-plane magnetoresistances of the material are strikingly different. Applying the field along the c axis suppresses the resistivity upturn, shifting the minimum to lower temperature. Yet even at 63 T the anomaly is still present. Consequently, at 50 K and below, ρc decreases when magnetic field is increased, with signs of saturation appearing around 50 T (Figure 3b). In contrast, in transverse magnetic field, ρc behaves as a more typical orbital magnetoresistance, common to metals. It is positive, and about three times weaker in magnitude, than the longitudinal one ( Figure 3c) without significantly affecting the shape of the upturn in ρc(T).
In order to emphasise the observed phenomenon's highly anisotropic nature, we also report the in-plane magnetotransport of 2H-NbS2 up to 14 T, presented in Figure 4. The magnetoresistance is weak for all field directions, but similarly to ρc, ρab is also reduced by the field along the c axis (clearly depicted in Figure 4b), which could be a trace of the same anomaly. The in-plane field-dependence of ρab ( Figure 4c

Discussion
The question of the origin of the anomalies of ρc in 2H-NbS2 will now be addressed. A number of phenomena could result in a finite-temperature resistivity minimum in a metal. Resistivity upturns can be caused by electronelectron interactions in the presence of static disorder [29][30][31] . However, the corresponding quantum mechanical correction is either weakly enhanced by a magnetic field, or is effectively field-independent. A closely linked phenomenon of weak localisation (WL) is also known to produce an additional contribution to resistivity at low 6 temperatures 32 . In this scenario, when a series of scattering events cause an electron to follow a closed path, quantum interference favours the net backward scattering over the forward one. Magnetic flux threading these scattering loops shifts the phases of the wavefunctions, diminishing the effect. In our case, when electrons are scattered between different planes, the closed paths should have comparable projections along the in-plane and out-of-plane directions. However, the upturn is only influenced by the c-axis-oriented field, contradicting the WLbased interpretation. A metal-insulator transition or conduction based on a thermally activated hopping between defects 33 would cause a divergence of ρc at the lowest temperatures, which was not the case. The possibility of quantum tunnelling playing a significant role is ruled out based on a linear relation between current and voltage (Supplementary Note 2 and Supplementary Figure 3). Resistivity upturns have also been observed in strongly doped cuprate superconductors 34,35 . In those materials, the effect is believed to be caused by scattering from magnetic droplets forming around non-magnetic impurities. This interpretation, however, relies on the existence of strong electronic correlations, and therefore does not apply to our system.
Finite-temperature resistivity minimum in a metal is also a well-known signature of the Kondo effect; a scattering of conduction electrons off dilute localised magnetic moments 36 . Besides the upturn, the characteristic features of the phenomenon, observable in charge transport, include a negative curvature of ρc(T) at the lowest temperatures and a suppression of the upturn by magnetic field, which causes spin-flip scattering to become inelastic 37,38 . The observed peak in the Seebeck coefficient is also characteristic to dilute and concentrated spin systems and, Kondo lattices 39,40 . It originates from the resonant scattering in the Kondo channel at the Fermi level. Above the Kondo temperature, TK, the resonance is smeared out, and depending on the specifics of a system, the peak in S appears at a temperature between 0.3TK and 0.9TK (Refs. 40,41). When temperature is low enough, the localized spins are We therefore argue that scattering off magnetic impurities is the most fitting explanation of our observations. The temperature dependence of ρc in 2H-NbS2 is consistent with the one expected from the numerical renormalisation group (NRG) theory calculations for Kondo effect 37,45 , as illustrated by the fit in Figure 3a. We modelled ρc with a sum of three contributions: a temperature-independent residual resistivity ρ0, an electron-phonon scattering term ρep (captured by the Bloch-Grüneisen formula), and the Kondo term ρK (the only magnetic-field-dependent term), for which we used the common empirical expression closely following the results of the NRG theory 37,45 : This explanation immediately raises a question regarding the nature of our system's magnetic impurities. The standard scenario where magnetic atoms are uniformly distributed clearly does not fit our picture. Doping 2H-NbS2 with Fe results in the upturn observable in ρab as well as the disappearance of superconductivity 46 .
Additionally, the undoped material does not display a corresponding signature in the heat capacity 47 . Lack of pronounced anomaly effects on ρab implies that the responsible defects take form of sparse planes, extending along the layers. When the current then flows along the layers, only a small fraction of the conduction electrons move in close proximity to these planes. But for the out-of-plane current flow, effectively all charge carriers have to pass through them, resulting in a particularly strong influence. Although we observed planar irregularities in the crystalline lattice via transmission electron microscopy, their atomic structure could not be determined due to a limited resolution (see Supplementary Note 5 and Supplementary Figures 7,8). Taking a closer look at the in-plane magnetoresistance for the c-axis-oriented field reveals that the difference between ρab at 0 and 14 T increases linearly with respect to ln(T) between 30 K and 15 K (Figure 4a inset), which further supports our hypothesis. This means that a minute contribution of Kondo scattering is present in ρab, but it is not strong enough to change the sign of the gradient of ρab(T).
The presence of 1T-NbS2 layers, evidenced by the characteristic CDW signatures 23 , offers a fascinating interpretation of our findings, illustrated in Figure 5. As we mentioned in the introduction, the √13 × √13 CDW order, associated with the 1T polytype, forms a triangular superlattice of David-star-shaped clusters defined by 13 Nb atoms 24 . The electronic structure of monolayer 1T-NbS2 as well as 1T-NbSe2 in such a configuration has been predicted to contain one very flat band around the Fermi level. This makes the materials susceptible to electronic instabilities like Mott localisation, with a concomitant magnetic order 11,48,49 . The referenced works found the ferromagnetic insulating state as the most stable, although others have proposed that such triangular lattices can host antiferromagnetic spin-liquid phases 50,51 . These magnetic planes play the role of scatterers in the Kondo effect.
The described scenario is conceptually similar to the Kondo effect occurring in artificially fabricated magnetic tunnel junctions 52,53 , yet in our case the phenomenon is observed in a spontaneously formed system. The same kind of Kondo interaction has been very recently observed in a 2H/1T-NbSe2 heterostructure, grown by molecular beam epitaxy 54 . One outstanding question is the anomaly's markedly different response to the two orientations of magnetic field. This difference is probably coming from the localised electron's highly anisotropic g-factor, causing a very small spin splitting (less than TK), but could also be related to the magnetic ordering. Sizeable anisotropy of the g-factor is expected for systems with strong spin-orbit coupling, such as TMDs 51 .
Since 2H-NbS2 is known to contain frequent 3R-like stacking faults, it is natural to ask whether the anomaly is somehow caused by the inclusions of 3R-NbS2. We measured the latter compound's out-of-plane resistivity, and while the corresponding temperature dependence was surprisingly found to be non-metallic, the extremely weak reaction of the interlayer conduction to the longitudinal magnetic field (Δρc/ρc ≈ 0.1% at 14 T) was incompatible with the behaviour observed in 2H-NbS2 (see Supplementary Note 6 as well as Supplementary Figure 9 for the relevant data on 3R-NbS2). The abundance of these stacking faults could explain the high residual component of ρc. The current understanding is that the poor conductivity of 3R-NbS2 is not intrinsic, but rather originates from the disorder due to self-intercalated Nb atoms 55 . However, in 2H-NbS2 the abundance of stacking faults results in high residual component of ρc and good in-plane metallicity.
In summary, we have demonstrated that a delicate alternation of the interlayer crystalline structure of 2H-NbS2 by introducing different polymorphs of the same atomic composition, dramatically affects the material's physical properties. In particular, the crystal's 2H stacking is occasionally disrupted by the 1T layers, which undergo a CDW instability. This then results in a triangular superlattice of David-star-shaped clusters, each hosting a lone 9 spin at the centre. Such a texture of localised magnetic moments can be seen as a two-dimensional Kondo lattice, After extraction, the lamella was glued to a sapphire substrate with a tiny amount of Araldite Rapid epoxy, keeping the external face exposed. Besides anchoring the lamella, the epoxy also formed a meniscus around it which smoothly connected the substrate's surface to the lamella's exposed face. The setup was then sputter coated with a 100 nm layer of gold. Next, the Ga FIB milling at 10 nA was used for defining the probing electrodes by selectively removing the sputtered gold layer, and for patterning the lamella in order to form the current channel and voltage probing points. The procedure was concluded with polishing the exposed side faces of the sample with a 1 nA ion beam in order to clean the surface of the re-deposited material and define the final dimensions of the device. Since the entire bottom face of the sample was rigidly attached to the substrate, differential thermal contraction and compressibility were expected to produce inhomogeneous stresses throughout the lamella. In our study, these stresses did not have a significant influence on the measured data. More detailed information about the FIB-assisted sample preparation can be found in the relevant review paper 7 and references therein.
Resistivity measurements. Resistivity was measured via the four-point technique with direct or alternating excitation currents in the 20-40 µA range. Temperature sweeps rate was limited to 1 K/min for the ambient pressure and of 0.5 K/min for the high-pressure measurements in order to reduce the thermal lag and gradients.
Resistivity at high pressure was measured using a piston cylinder cell produced by C&T Factory. Daphne oil 7474 was used as a pressure-transmitting medium. Pressure was determined from the changes in resistance and superconducting transition temperature of a sample of Pb located next to the 2H-NbS2 sample.
Measurements in high magnetic fields were conducted at the high magnetic field facilities in Grenoble (up to 34 T DC field) and Toulouse (up to 63 T pulsed field). Quantum Design PPMS was used for measurements in fields up to 14 T.
Seebeck coefficient measurements. Seebeck coefficient was measured using an in-house setup. For the in-plane Seebeck coefficient measurement, a thin and long sample was mounted on a ceramic bar. One end of the bar was connected to the thermal bath, while the other one had a resistive heater attached. A differential thermocouple was used to measure the temperature difference across the sample. The out-of-plane Seebeck coefficient measurement was performed using a setup displayed in Figure 2c and described in the corresponding caption.

Data availability
The data that support the findings of this study are available from the authors (E.M. and K.S.) upon reasonable request.  (c) 20 out-of-plane (Sc) and in-plane (Sab) directions as functions of temperature, measured on bulk single crystals (note that Sab is negative). The setup for measuring Sc is illustrated schematically. The crystal was approximately 1 mm long in the c axis direction (indicated in the drawing), and 2-3 mm long laterally. The value of Sc is the ratio of the voltage across the sample (Vs) and the thermal gradient across it, determined from the differential thermocouple voltage (VDTC). The sample sat between two copper plates, which homogenised temperature at its two faces and was electrically decoupled from the heatsink by a thin sapphire plate.  The results are shown in Supplementary Figure 1 and are in agreement with various parts of the electronic structure of the compound presented in earlier reports [6][7][8] . Electronically, 2H-NbS2 is nearly identical to 2H-NbSe2 9,10 . Fermi surfaces of both compounds include a three-dimensional hole pocket at the Г point of the Brillouin zone, bearing the chalcogen pz orbital character, and two kinds of quasi-tubular Fermi sheets originating from Nb d orbitals.

Figures
Based on these elements, particularly the hole pocket, one would expect the charge transport properties of 2H-NbS2 and 2H-NbSe2 to be nearly identical. is approximately the same (in the normal states), 2H-NbS2 has a much higher residual resistivity, shifting its ρc to higher values over the entire temperature range. This implies a significantly larger content of static defects in the c axis direction, compared to 2H-NbSe2.

Direct comparisons between ρc(T) and ρab(T) measured in 2H-
We therefore conclude that the low-temperature out-of-plane resistivity anomaly of 2H-NbS2 is not an intrinsic property originating from the nominal parameters of the crystal, but rather a consequence of either frequent stacking faults (only creating disorder along the c axis) or inclusions of alien phases, which extend along the layers and are rare enough to not produce signatures in the in-plane charge transport. b, Out-of-plane resistivities, same as the ones plotted in Figure 2 of the main text, but without any multiplicative scaling. The large difference between resistivity values is evident, while both curves show a very similar temperature dependence above 50 K. The CDW transition in 2H-NbSe2 has an extremely weak effect on the out-of-plane resistivity.
We considered a possibility that the out-of-plane resistivity anomaly in 2H-NbS2 is caused by the conduction partially occurring via a quantum tunnelling across non-conducting defects. Such a conduction results in a nonlinear IV curve, however out measurements of it showed no detectible deviations from linearity (Supplementary Supplementary note 4. Fitting the out-of-plane resistivity of 2H-NbS2 with a Kondo effect model.
Here we provide a detailed explanation of the approach used for fitting the temperature dependences of the outof-plane resistivity at different magnetic fields with the theoretical model for Kondo effect.
The overall temperature dependence of the out-of-plane resistivity is determined by multiple contributions. These can be separated into the conventional part, with a typical metallic resistivity, monotonically growing with temperature, and the anomalous component, resulting in the low-temperature upturn. The analysis is simplified if the anomalous contribution can be easily separated from the conventional one, by using the high-temperature section of the data (unaffected by the anomaly) for constraining the conventional part. In the case of 2H-NbS2 this approach does not work, since the effect of the anomaly is noticeable up to about 90 K, but just above this temperature the resistivity saturation already starts coming into effect.
We therefore focused on the section of the data below 90 K and assumed that resistivity in the selected region was dominated by three terms: a temperature independent residual resistivity ρ0, a contribution due to electron-phonon scattering ρep, and the anomalous component ρK.
We modelled the electron-phonon scattering term with a Bloch-Grüneisen formula 15,16 : The parameter Θ is the characteristic temperature, which is usually well approximated by Debye temperature. The exponent n is often equal to 5, but may take different values depending on specifics of the scattering 17 . C is a temperature-independent multiplier.
For capturing the anomalous contribution, ρK, presumably arising due to the Kondo effect, we used the common empirical expression, which precisely captures the temperature dependence predicted by the numerical renormalisation group theory 18,19 : The parameter α is related to the spin of magnetic impurity, ρK0 is the zero temperature limit of ρK, and TK is the characteristic temperature (Kondo temperature), such that ( ) = 0 2 ⁄ .
The model described above has 6 parameters that have to be fitted: ρ0, C, n, α, ρK0 and TK. Such a large number may lead to an underconstrained fit. To alleviate this problem, the least squares fit was conducted globally, for the data at different magnetic fields values (0 T, 10 T, 20 T, 30 T, 45 T, and 63 T), with only ρK0 being allowed to vary with field. For Θ, we used the Debye temperature value of 259 K, obtained from the low-temperature heat capacity data of 2H-NbS2 20 .
The resultant best fit to the data is plotted in Supplementary Figure 5a. While the model does a good job at describing the measured data, if one considers the individual components ρ0, ρep, and ρK separately, then the contribution (ρ0 + ρep) appears too small compared to the level expected from the visual extrapolation of the data to higher fields. If n is set to the more common value of 5, the fit becomes more realistic, however a substantial discrepancy between the measured and fitted data is visible (Supplementary Figure 5b). Setting the value of n to 3 (the Bloch-Wilson formula 17 ) gave a noticeably better fit than for the other choices, with only a minor deviation from the measured data at 50 K for the 10 T and 20 T data (Supplementary Figure 5c, this fit is displayed in Figure   3a of the main text). The corresponding values of the fitting parameters are provided in Supplementary Table 1.
The resultant value of A was vanishingly small, meaning that the term proportional to T 2 is not needed for adequately describing ρc of 2H-NbS2 The value of 1.71 obtained for α is an order of magnitude larger than the typical numbers expected for spin-1/2, 1, Our proposed idea, that the Kondo effect in 2H-NbS2 is caused by planar defects, indeed bears a resemblance to the case of a magnetic tunnel junction.
The fitting was performed with an assumption that the only effect of magnetic field was the reduction of ρK0 and TK. However, one could expect additional longitudinal magnetoresistance, associated with the conventional metallic resistivity. To estimate the magnitude of this component, we refered to 2H-NbSe2. Supplementary Figure   33 6 shows how the out-of-plane resistivity of 2H-NbSe2 is affected by the longitudinal magnetic field.
Magnetoresistance is positive and increases on cooling down, reaching nearly 10% at 10 K and 14 T. When comparing magnetotransport in different systems, one should consider the dependence not simply on B, but rather on ωcτ (where ωc is the cyclotron frequency and τ is the mean scattering time). Since at low temperature ρc of 2H-NbSe2 is around an order of magnitude lower than that of 2H-NbS2, we would expect the magnetoresistance of the latter material to vary an order of magnitude more slowly as a function of B. Such a contribution to magnetoresistance would not be significant compared to the observed suppression of ρK.
To summarise, the analysis presented here shows that based on the measured temperature and magnetic field dependences of the out-of-plane resistivity of 2H-NbS2, Kondo effect is a plausible interpretation of the observed resistivity upturn.  Supplementary Note 5. Evidence for the c-axis disorder from transmission electron microscopy.
Transmission electron microscopy (TEM) was used in an attempt to directly detect and characterise structural defects in our crystals. The results obtained confirm their presence but allow only a limited qualitative evaluation and more work is needed to fully understand their structure and topology. For our experiments a lamella of 2H-NbS2, with its plane perpendicular to a 〈112 ̅ 0〉 direction was extracted from a single crystal and polished down to sub 50 nm thickness with focused ion beam 23 . First, by analysing a transmission electron diffraction image, we concluded that majority of reflections belonged to both 2H and 3R-types of stacking, but some spots could not be accounted for (Supplementary Figure 7).
This result is in agreement with the earlier descriptions of this material 24,25 . In the 2H and 3R structures the layers are identical, but their stacking sequences are different. Due to the lower symmetry of 3R-NbS2, two different orientations of it needed to be considered (green and purple dashed lines). The pattern had to be vertically stretched by about 9% in order to get the best match with the expected spot positions.
High-angle annular dark-field imaging (HAADF) scanning TEM images revealed extended homogeneous regions separated by two kinds of planar defects, comprising a single or three dark lines running along the layers (Supplementary Figure 8a). In the former case, the intensity profile on two sides of the defect was shifted by half a period, while no such change occurred for the latter one (Supplementary Figure 8b). The phase shift indicates that the single-line defect represents a layer of increased thickness. It was not possible to achieve atomic resolution and deduce the metal-chalcogenide coordination. This would have been necessary in order to identify the 1T-type layers, where metal is octahedrally coordinated by suphur, as opposed to the corner-sharing trigonal prismatic coordination for 2H and 3R polytypes.
Dark appearance of the defects could be caused by a variety of mechanisms, such as stoichiometry variation, diffraction contrast or even local magnetic order. The fact that these defects take form of isolated and laterally extended planes makes Nb deficiency an unlikely explanation, as the vacancies are expected to be distributed more randomly. A more fitting explanation would be the atomic disorder/displacement, which can originate from a 1T-NbS2 layer with the √13 × √13 CDW order. Traces of this particular lattice reconstruction have been observed 39 in 2H-NbS2 using diffusive X-ray diffraction 25 . Increased thickness of the defective layer also agrees with the scenario proposed above, as the associated distortion of the lattice is known to increase the interlayer spacing 26 .
The role of local magnetism in generating the contrast cannot be excluded. This hypothesis could be tested by further studies by TEM with differential phase contrast.
The larger than expected interlayer separation observed in the studied lamella remains an open question, and additional studies are necessary in order to make conclusions about the topology of stacking faults.

Supplementary note 6. Out-of-plane resistivity of 3R-NbS2
Since 2H-NbS2 is known to contain frequent stacking faults, locally appearing as 3R-NbS2, we also probed resistivity of the latter material along the two directions. The available single crystals of the 3R polytype had 1.5% of S substituted by isovalent Se. An isovalent substitution in TMDs can strongly tune the stability of often competing CDW and superconducting states, but generally has little effect on the normal state charge transport properties. Since neither 3R nor 2H-NbS2 develop any CDW order, the aforementioned doping is of minor relevance. The measured out-of-plane and in-plane resistivities, as well as their ratio, are plotted against temperature in Supplementary Figure 9a. The shape of the in-plane resistivity is identical to the ones reported for undoped 3R-NbS2 [27][28][29] Low residual in-plane resistivity ratio of 1.3 implies a high content of static defects, yet similarly low values have been found in the earlier studies. Below 30 K, the in-plane resistivity exhibits a very weak upturn, commonly observed in this compound and explained by electron-electron interactions in presence of static disorder 27 . Surprisingly, the out-of-plane resistivity has a negative and nearly constant gradient throughout the entire explored temperature range. Resistivity anisotropy of 3R-NbS2 goes from 6 to 9 upon cooling. The