Quantum phases and spin liquid properties of 1T-TaS2

Quantum materials exhibiting magnetic frustration are connected to diverse phenomena, including high Tc superconductivity, topological order, and quantum spin liquids (QSLs). A QSL is a quantum phase (QP) related to a quantum-entangled fluid-like state of matter. Previous experiments on QSL candidate materials are usually interpreted in terms of a single QP, although theories indicate that many distinct QPs are closely competing in typical frustrated spin models. Here we report on combined temperature-dependent muon spin relaxation and specific heat measurements for the triangular-lattice QSL candidate material 1T-TaS2 that provide evidence for competing QPs. The measured properties are assigned to arrays of individual QSL layers within the layered charge density wave structure of 1T-TaS2 and their characteristic parameters can be interpreted as those of distinct Z2 QSL phases. The present results reveal that a QSL description can extend beyond the lowest temperatures, offering an additional perspective in the search for such materials.


INTRODUCTION
The idea of destabilizing magnetic order by quantum fluctuating resonating valence bonds (RVB) originated with Anderson in 1973 1 . Systems with frustrated interactions are particularly susceptible to these effects and much work on highly frustrated systems and their quantum spin liquid (QSL) phases has followed. In parallel with extensive theoretical work 2-4 , many candidate QSL materials have been proposed. For triangular lattices these include molecular [5][6][7][8] and inorganic [9][10][11] compounds. However, to date, little experimental effort has been devoted to testing for crossover between different topological phases in a QSL, or to exploring the nature of the fractionalized spinon excitations (bosons, fermions, or anyons) 12 . This report concerns 1T-TaS 2 , which is a layered compound with a series of charge-density-wave (CDW) instabilities terminating in a fully commensurate CDW phase (C-CDW) below 200 K. In the C-CDW state the distortion pattern forms a triangular-lattice of Star-of-David (SD) clusters, each cluster containing 13 Ta atoms and one unpaired spin-½ (Fig. 1a). It shows no evidence for magnetic ordering and was one of the triangular-lattice materials that originally inspired the RVB theory 1 .
The electronic and magnetic scenario offered by 1T-TaS 2 is far from being simple and is currently being revisited. Theoretical studies have shown that a single layer of 1T-TaS 2 is indeed a QSL 13 and experimental results in the selenium counterpart, 1T-TaSe 2 , have revealed that single layers are Mott insulators that behave as QSLs 14,15 . However, when considering bulk crystals of 1T-TaS 2 , where the layers are stacked in the out-of-plane direction, new challenges for the interpretation of its properties appear. Taking into account band theory, 1T-TaS 2 should be a metal, since there are an odd number of electrons per unit cell. As it has been experimentally observed to have an insulating state below 200 K, the Mott mechanism has been invoked to explain this apparent contradiction and, thus, the possibility of a QSL has emerged 13 . However, if the layers form dimers in the out-of-plane direction, the unit cell doubles and there are an even number of electrons, thus transforming the system to a conventional band insulator where the QSL scenario has no place 16,17 . Several recent studies have made further exploration of the evidence that 1T-TaS 2 is a QSL, both from experimental and theoretical perspectives 13,[18][19][20][21] , but some alternative scenarios have also been proposed, such as a Peierls mechanism or the formation of domain wall networks, among others [22][23][24] .
Although the QSL may seem directly linked to the presence or absence of dimerization of the layers in the out-of-plane direction, there is still another possibility, since the out-of-plane stacking order can be arranged in different configurations and exhibit a temperature dependence, as recently proposed [25][26][27] . Thus, below 200 K the layers can start off being uncoupled, allowing a QSL to exist, but then they can become coupled progressively while cooling down, moving towards a band insulator at low temperatures if all the layers become dimerized. However, at intermediate temperatures, there is a scenario where observed QSL properties may be the result of coexisting dimerized and undimerized layers. This relatively unexplored picture is examined more closely in the present study, where we report on muon spin relaxation (μSR) and specific heat as a function of the temperature on high-quality samples of 1T-TaS 2 . We observe distinct regions versus temperature for both of these properties. We discuss these different observed thermal regions in relation to possible changes in the out-of-plane stacking order of the layers and in relation to some models of closely competing QSL states that are currently found in the literature.

RESULTS AND DISCUSSION
μSR data μSR provides a good way to check a candidate QSL for magnetic ordering and several previous reports using zero-field μSR (ZF-μSR) have confirmed the absence of magnetic ordering in 1T-TaS 2 down to 20 mK [18][19][20] . The muon probe is fully spin polarized on implantation and the forward/backward asymmetry of the detected muon decay positrons, a(t), reflects the timedependent polarization of the muon spin. In the case of ZF (Fig.  2a), the relaxation of a(t) reflects the presence of multiple asymmetry components, corresponding to different local environments. Diamagnetic environments show characteristic nuclear dipolar relaxation, which is initially Gaussian, whereas paramagnetic environments typically show a Lorentzian (exponential) relaxation due to the dynamics of the unpaired electron. This electronic contribution is in a fast fluctuation regime where the relaxation rate λ is proportional to the electronic spin fluctuation time τ 28 . μSR can also measure critical parameters 6 that can be compared against QSL models. A QSL can support fractionalized spinon excitations and various experimental properties will depend on the characteristics of these spinons. Here we use μSR in longitudinal magnetic fields (LF-μSR) to quench both the diamagnetic relaxation and any slow local dynamics. The relaxation then becomes a slow exponential (Fig. 2a), reflecting the dynamics of spinons rapidly diffusing through the triangular lattice (Fig. 1a). Further details are given in the "Methods" section. The coupling of the muon to the spinons is determined by the stopping sites (Fig. 1b) and the associated hyperfine interactions at these sites, which we have computed using density functional theory (DFT), further details are given in the Supplementary Methods.
The μSR relaxation rate is relatively slow in this material (Fig. 2a), as found previously [18][19][20] . A detailed study of λ(B LF ) at fixed T provides information about the spin fluctuations 29 . This method has been used in μSR for studying the propagation of spinons in spin-½ Heisenberg antiferromagnetic chain systems 30 , where 1D diffusive motion leads to a weak power law dependence for λ(B LF ), and the method can be extended to systems of higher dimensionality 29 . Further details are given in the Supplementary Methods. Figure 2b shows example data at 120 K, fitted against spinon diffusion models on 1D, 2D, or 3D lattices. The 2D model clearly provides the best fit and global analysis of the LF data between 6 and 180 K shows consistency with the 2D model across the whole data set (Fig. 2c). Figure 3a shows the 2D diffusion rate D 2D obtained from fits to the LF-μSR scans. Several distinct thermal regions are observed and effective thermal power laws are used to describe the data in the different regions. The power law parameterizations will be useful for making comparison with QSL models, which make their predictions in terms of power laws. On cooling through T nC and entering the C-CDW phase, D 2D rises significantly, consistent with the presence of a QSL state that supports diffusing spinons. This higher temperature region of the C-CDW phase is labeled III and on cooling, D 2D continues rising until a maximum occurs at T 0~1 10 K. The region below T 0 labeled II shows a relatively strong power law D 2D / T nD with n D = 1.74 (14). Below T 1~2 5 K another region, labeled I, is found, in which n D = 0.47 (17). Figure 3b shows f 0, the amplitude parameter of the diffusive signal, which is the product of f diff , the fraction of spectral weight in the diffusion process, and f para , the fraction of muon sites coupled to unpaired electronic spins. The values for f 0 have been obtained using the muon hyperfine coupling calculated from DFT analysis Ā = 4.5(6) MHz (see Supplementary Methods  and Supplementary Table 2). It can be seen that the f 0 parameter also shows multi-region behavior, with region boundaries matching those of D 2D .

Specific heat data
An independent measurement of the behavior of the quantum phase (QP) regions I, II, and III is given by the power law electronic contribution to the specific heat C e1 p / T nC . Details of the measurements and their analysis are provided in the "Methods" section. These provide a set of values of n C for each region (Fig. 4). The overall behavior of the specific heat is summarized in Fig. 4a. The border between regions II and III is clearly signified by a maximum value of C p /T at T 0 , corresponding to the distinct peak found in the μSR parameters of Fig. 3. In contrast to the strong feature at T 0 , the border between regions I and II at T 1 is not so obvious because the lattice term is rapidly varying in this region. In the low T region the lattice contribution is taken to follow the standard asymptotic T 3 Debye form and the power law n C for the electronic contribution in region I can then be obtained from a fit, as shown in the associated scaling plot of Fig. 4b. The obtained value n C = 1.46 is significantly larger than expected for a spinon Fermi liquid, where the predicted value is 1, reducing to 2/3 when the system is close to a quantum critical point (QCP) 31 . For regions II and III the lattice contribution is modeled by the sum of an anisotropic Debye term and two Einstein terms, as detailed in the "Methods" section. The peak in C/T is then seen to be the result of n C reducing by around one on warming through T 0 , as shown in Fig. 4c. A consistent parameter set for the different regions is summarized in Table 1.

Analysis
We consider here possible interpretations of the different regions observed in these measurements. First, we look at how the experimentally obtained power laws compare with the critical parameters of appropriate QSL models. Next, we consider the nature of the transitions between the regions and a possible interpretation in terms of a delicate balance between competing QSL phases, as found in numerical studies of the triangular lattice Heisenberg antiferromagnet. Finally, we consider the role of the degree of dimerization of the layers defined by the interlayer stacking of the SD clusters.
Critical parameters and QSL models The many possible 2D QSL states for a triangular lattice can be broadly classified as SU(2), U(1), or Z 2 3,13 . Of these, the Z 2 states are expected to be most stable 4,13 and there are well-developed theories for describing the quantum critical (QC) properties of Z 2 QSL phases [32][33][34] . Thus, Z 2 models provide a natural starting point for interpreting the data. A key feature of QC phases is that they occupy extended regions of T, in contrast to classical critical regions, that generally require T to be very close to the transition point. For a Z 2 QSL the spinon transport in the QC regime 34 is determined by a momentum relaxation time that follows power law in T of 2/ν − 3, where ν is the correlation length exponent. Applying the Einstein diffusion relation, the power law for the diffusion rate is then larger by one, i.e. n D = 2/ν − 2. Besides this critical power law, there may also be present an additional noncritical contribution, related to an energy dependence for the spinon density of states (DOS). This can be taken into account by including an additional parameter q, i.e. n D = 2/ν − 2 + q. The value q = 1 describes a QSL with a linear DOS, such as would be found for a nodal excitation spectrum, whereas q = 0 for a QSL with constant DOS. The ν exponent in the QC model [32][33][34] is that of the O(4) criticality class, calculated to be 0.748(1) 35 . For the specific heat, the corresponding critical exponent α can be obtained from ν via the hyperscaling relation α = 2 -3ν, giving α = −0.244(3).
The power law for the specific heat is then taken to follow the form n C = 1 − α + q.
Experimental values for n D and n C are compared with the model values for q = 0 and q = 1 in Fig. 5a. The corresponding values for ν and α are shown in Fig. 5b. For region I the experimental values derived from both LF-μSR and specific heat are reasonably consistent with the Z 2 QSL model with q = 0. Data for region II requires q = 1, suggesting a Z 2 -linear QSL. In region III, it is found that q = 0 again gives the closest match to the n C value.
We note that the strongly negative n D value measured in region III (Fig. 3a) is not predicted within this model (it would require an unphysical negative value for ν), suggesting that the QC response is overwhelmed by another effect that is specific to LF-μSR. We have considered whether muon diffusion might be responsible for this behavior, but the energy barriers between adjacent sites in the ab planes and also from one side of the layer to the other are found to be too large to support diffusion in this temperature range (see Supplementary Methods, Supplementary Fig. 3 and  Supplementary Table 3). It is more likely that the fall reflects interplane spin correlations sensed by the muon (see Supplementary Note and Supplementary Fig. 6), as well as the influence of charge excitations, which would be expected to become significant in the region that borders on the nC-CDW phase. We note that a negative power might be the result of a U(1)-gapless QSL state with fermionic spinons 36 , but the predicted exponent n D = −1/3 is very much smaller than measured here (Fig. 3a) and the corresponding specific heat prediction 31 of α = 1/3 is not found (Fig. 5b).

The nature of the phases and their transitions
One question concerns the nature of the transitions between the regions, in particular whether they are crossovers or phase transitions. The absence of any sharp features in the specific heat in the whole C-CDW/QSL region rules out first-order phase transitions. A conceptual framework that could be used for interpreting the transitions in terms of underlying quantum phase transitions (QPTs) 37,38 is illustrated in Fig. 6. The phase diagram has a magnetic (AF) region and a non-magnetic QSL region that is separated into a high T QC phase and low T quantum disordered (QD) phase. We assign the three experimentally observed phases to the QSL-QC region. The QCP separating the spin liquid from the  Spin diffusion model c Fig. 2 Muon spin relaxation. a LF relaxation for two different fields at 55 K, with the ZF relaxation shown for comparison. Error bars indicate standard errors defined by the counting statistics. b LF dependence of the relaxation rate at 120 K with fits to 1D, 2D, and 3D spinon diffusion models. Error bars indicate standard errors of the fit parameters. c Quality of fit, as described by the χ 2 parameter for global analysis of the set of λ LF values measured at 13 temperatures between 6 and 180 K, comparing different diffusion models. The expected value of χ 2 for a good model fit is the number of degrees of freedom N D = 65, defined as the difference between the number of data points (91) and the number of free parameters in the fit (two for each temperature). N D is shown here as the solid line, with dotted lines indicating the expected χ 2 distribution width. From this analysis, we conclude that the 2D model is fully consistent with the data over the full temperature range, whereas the 1D and 3D models are not.
AF ground state lies at a specific value of control parameter g 1 (Fig. 6a). A possible association of g 1 could be made with ringexchange, which is known to destabilize the magnetic phase of triangular-lattice Mott-Hubbard insulators 39 . The energies of the competing QSL states are assumed to depend on a second control parameter g 2 and at T = 0 continuous second-order QPTs would be expected between pairs of QSL-QD states (Fig. 6b).
Moving to the finite T case, appropriate to our measurements, several additional factors come into play. Firstly, since the excitation spectrum is different between the adjacent QSL phases, as indicated by their difference in q, there will in general be a difference in entropy between the phases at the transition temperature, which would tend to turn the continuous secondorder transition into a discontinuous first-order transition. A second factor is that at finite T we expect the QPTs to evolve into broader QC regions separating the QD phases, with the transitions becoming crossovers. In the current system it appears that the second mechanism dominates over the first. In this scenario, g 2 is taken to have a T dependence that produces smooth crossover transitions between the QSL states. The finite temperature transitions at T 0 and T 1 then would directly correspond to multicritical QCPs in a suggested underlying (g 1 , g 2 ) quantum phase diagram for T = 0, as shown schematically in Fig. 6b.
Comparison with a QSL phase diagram Spinon pairing instabilities leading to Z 2 spin liquid states have been discussed previously in several studies [40][41][42][43] . In general, there are 64 possible symmetric Z 2 spin liquid states for the triangular lattice [44][45][46] . Some of the more stable states have been revealed via variational studies and it is useful to discuss the multiple phases revealed by our measurements in relation to previously reported phase diagrams from these studies. The most comprehensive study to date is by Mishmash et al. 40 . This uses a model Hamiltonian that includes terms representing the nearestneighbor exchange J, the next-nearest-neighbor exchange J 2 , nC-CDW IC III II I Fig. 4 Specific heat. a Measured specific heat over T plotted across a wide T range. The border between regions II and III is clearly signified by a maximum in C p /T. Within the C-CDW region, the data are fitted by the sum of contributions from the crystal lattice and the electrons. b For the low-temperature region I the lattice contribution follows its asymptotic T 3 form and the power law exponent for the electronic contribution n C is straightforwardly obtained from a fit to the data below 10 K. c In regions II and III the lattice contribution is large, but the T dependence of the electronic term can still be modeled using two effective values for n C . It is notable that region II has a significantly larger n C value than both regions I and III. Error bars indicate the estimated standard error of the measurement. and the ring-exchange K. The parameters J 2 and K are expressed in units of J. There is general agreement from different studies that for this type of model the AF state gives way to the U(1) spin liquid state with a spinon Fermi surface when K exceeds values of order 0.2-0. 25 (refs. 21,39,40,43 ). At lower K a wider range of states have been suggested. Figure 6c shows a representation of the part of the phase diagram obtained by Mishmash et al. 40 that is most relevant to a Mott-Hubbard insulator, where the exchange parameters are determined by the U/t ratio. In this study a BCSlike d-wave paired Z 2 spin liquid was found over an extended region of the phase diagram (Fig. 6c) 40 . The properties of this nodal state were discussed by Grover et al. 43 . Below a characteristic pairing temperature it shows a linear (q = 1) DOS. In their study Mishmash et al. 40 considered six possible pairing modes and found a second paired state, labeled d + id, that is very close in energy to the nodal state and the most stable state for a small region of parameter space indicated in Fig. 6c. Having such strongly competing states corresponds well with the picture put forward in Fig. 6b and a clear mapping of the phenomenological g 1 and g 2 control parameters to the K-J 2 space is suggested by the structure of the phase diagram (Fig. 6c, dashed lines). The d + id state has quadratic dispersion around the central Γ point 40 which gives a constant DOS state with q = 0. A possible cooling path and assignment of the two lowest temperature thermal phases is indicated by the purple arrowed line in Fig. 6c. Within this picture phase II could be assigned to the q = 1 nodal d-wave state and phase I could be assigned to the q = 0 quadratic d + id state.
The scenario of Fig. 6c would place K in the region of 0.12, corresponding to a U/t ratio of order 11. This can be compared to the U/t estimates of 7 to 8 for the well-studied triangular-lattice QSL system κ-(ET) 2 Cu 2 (CN) 3 , with corresponding K values in the region 0.5-0.3. Although model phase diagrams such as Fig. 6c provide a good starting point for discussion, further distinct properties of 1T-TaS 2 need also to be taken into account for a more complete picture, as discussed in the next section.

Interlayer stacking
The mode of interlayer stacking of the SD clusters has an important bearing on the local electronic properties and, in particular, the emergence of insulating behavior. However, this aspect of 1T-TaS 2 remains poorly understood. Recent experiments address possible origins for the insulating state, from the conventional Mott scenario 19 to layer dimerization 47 , unit-cell doubling 24 , or a Peierls transition 22 . Synchrotron X-ray studies reveal significant disorder in the interlayer stacking of the SD clusters at 100 K from the presence of diffuse scattering 47 , suggesting a transition from a band insulator at low temperatures to a Mott insulator at higher ones 26 . Besides a broad diffuse background, diffuse peaks were observed at positions suggesting correlations in the stacking sequence on length scales of three and five layers. In contrast, data taken for the nC-CDW phase at 300 K show a simple three-fold periodicity in the interlayer stacking 47 . A scanning tunneling microscopy (STM) study 24 found two types of cleaved surface at 77 K: large-gap surfaces that were assigned to spin-paired bilayers and narrow-gap surfaces assigned  to unpaired Mott-gapped layers, that could support the QSL state. The large-gap surfaces were assigned to surface bilayers (A-mode stacking with SDs aligned between layers). The weaker-gap surfaces were assigned to a surface layer stacking denoted as a C-mode, for which the central Ta site of the SD of the lower layer is located below the Ta site in one of the tips of the SD of the upper layer. A regular stacking sequence of the form ACACAC is one possibility (Fig. 7a). If the electronic stacking along the c-axis is more complex, e.g. taking the form ACACCAC (or, in general, … ACA⊥⊥AC…, where ⊥ is a layer stacking that differs from A), then a fraction of the cleavage planes will be unpaired monolayers. The experimental ratio of three to one found in the STM study 24 would then be consistent with an average density of electronic stacking defects around this one in seven level (Fig. 7b). For muons implanting in such a sequence, two out of seven interlayer sites will see the unpaired layer and couple to the paramagnetic electronic relaxation of the QSL layers and the remaining five sites will see the diamagnetic spin-paired bilayers, showing only nuclear relaxation, thus f para = 2/7 in this case. Figure 7c, d shows examples of local configurations with increasing QSL layer concentrations and correspondingly larger values for f para .

Alternative scenarios for the higher T regions
The f 0 parameter shown in Fig. 3b can provide valuable information about evolution of the electronic state with T. In region I, it is essentially constant, indicating a stable electronic state and allowing the QC model to be applied with confidence. On the other hand, the significant changes of f 0 seen in regions II and III provide a greater challenge for the overall interpretation.
The largest value of f 0 is found at T 0 , where it approaches its maximum possible value of one, indicating that both f para and f diff approach one here. The reason for the drop in f 0 in region III above T 0 is not entirely clear, but besides reflecting AF interlayer correlations at the muon site, it may also reflect the increasing influence of charge excitations, which are not included in the basic spin-only model. The origin of the drop in f 0 on the lower side of T 0 when cooling through region II is also a key question. In one scenario, we could assume that f diff is one at all T, so that the drop in f 0 is entirely due to a reduction in f para at low T. This could be consistent with the STM study 24 , as discussed in the previous section. In this scenario one could propose that the q = 1 value obtained for the specific heat in region II would not reflect a QSL state with a linear DOS, but rather a linear increase in the concentration of the QSL layers above T 1 , reaching a maximum at T 0 . This scenario would not however be consistent with the observed large change in n D at T 1 , since the obtained 2D diffusion rate is a property of an individual QSL layer and it would not be expected to depend directly on f para . The opposite scenario for the drop in f 0 within region II is that f para remains at or close to one, so that the changes in f 0 observed in region II are primarily due to changes in f diff . In this case the q = 1 value in region II would reflect an underlying QSL having a linear DOS, as discussed earlier. This latter scenario would then be consistent with both the diffusion rate and specific heat data obtained in the present study.

Summary
The main result of our study is that the combined local-probe LF-μSR and bulk specific heat measurements in 1T-TaS 2 within the C-  Fig. 6 Phenomenological interpretation and comparison with theory. a The QC spin liquid model has a QCP between an AF phase, in this case the 120°non-collinear state (thick line) and a gapped quantum disordered QSL-QD Z 2 spin liquid phase. The QPT is tuned by a control parameter g 1 and at the QCP, the spin stiffness ρ s and the QSL gap Δ both reduce to zero. The experiments at finite temperature T indicate QSL phases that are gapless and thus close to being directly above the QCP in the QSL-QC region, suggesting that the actual g 1 value closely matches the position of the QCP and does not have a strong T dependence. The fan-shaped QC region is bounded by three crossovers (dashed lines). On the low side, these reflect ρ s and Δ. The exchange coupling J would normally provide the high limit, but T nC actually occurs first here. b A second parameter g 2 is introduced to reflect the stabilization of different QSL phases. c A reported triangular lattice phase diagram from variational studies 40 . The control parameters are the second-nearest-neighbor exchange J 2 and the four-site ring exchange K, both normalized to J. Close to the AF state, the d + id state and the nodal d-wave state strongly compete, leading to a pocket of stability for the d + id state. Dashed lines show possible J 2 -K mapping of the parameters g 1 and g 2 . The purple arrow shows a possible parameter path for 1T-TaS 2 on cooling through regions II and I.
CDW regime support the existence of QSL layers within the layered structure. The overall thermal behavior is however found to be complex, with three distinct thermal regions observed in the data. For region I, the most stable region that occurs at low temperature, it is found that both μSR and specific heat are consistent with the QC regime of a Z 2 QSL model. The specific heat in region III is also consistent with this model. In contrast, region II shows an additional T-linear contribution to the measured properties, which most likely reflects the presence of a competing QSL state with a nodal excitation spectrum. Such competing nodal phases have previously been found in various theoretical studies of the S = 1/2 triangular lattice Heisenberg antiferromagnet. Compared to the usual triangular lattice QSL systems with fixed lattice structure for the spin sites, the 1T-TaS 2 system has a specific additional degree of freedom related to weak interlayer correlations affecting the mode of stacking between the SD layers. Any tendency within the stacking to form dimerized bilayers will lead to a non-magnetic contribution from these bilayers. From our data, we conclude that in our sample a large fraction of the layers are undimerized at the boundary between regions II and III and still only a minor fraction become dimerized at low temperature.

Crystal growth
High-quality 1T-TaS 2 bulk crystals ( Supplementary Fig. 1a) were synthesized by chemical vapor transport and characterized using elemental analysis, ICP mass spectrometry, and XRD. The crystal growth conditions were as follows: first, powders of Ta (99.99%, Alfa-Aesar) and sulfur (99.8%, Sigma-Aldrich) were mixed in a stoichiometric ratio, sealed in an evacuated quartz ampoule (length: 250 cm, internal diameter: 14 mm; P: 2 × 10 −5 mbar) and heated from room temperature to 910°C in 3 h. The temperature was then kept constant for ten days and finally allowed to cool down naturally. In order to grow larger crystals, 1 g of the obtained material was mixed with iodine (99.999%, Sigma-Aldrich; [I 2 ] = 1 mg/cm 3 ), and sealed in an evacuated quartz ampoule (length: 500 mm, internal diameter: 14 mm; P: 2 × 10 −5 mbar). The quartz tube was placed inside a three-zone furnace with the material in the leftmost zone. The other two zones were heated up in 3 h from room temperature to 950°C and kept at this temperature for 2 days. After this, the leftmost side was heated up to 925°C in 3 h and there was established in the three-zone furnace a gradient of 950/900/925°C. This temperature profile was kept constant for 20 days and the final stage was a rapid quench into water.
The weight concentration of elements obtained by inductively coupled plasma optical emission spectroscopy (ICP-OES) was Ta: (74.6 ± 1.4)% and S: (25.6 ± 0.5)%. These values are in good agreement with the expected ones (Ta: 73.8% and S: 26.1%). Phase purity was confirmed by the refinement of the X-ray pattern (measured with a PANalytical Empyrean Xray platform with Cu radiation source) to the previously reported structure of 1T-TaS 2 (ICSD 85323) using the X´Pert Highscore Plus program. The hexagonal crystal system with a P-3m1 space group was obtained and the unit cell was determined to be α = β = 90°, γ = 120°, a = b = 3.3662(4) Å and c = 5.8975(7) Å. The obtained results are in accordance with the ones reported in the literature 48 . The density of structural stacking faults was estimated to be below 2% from the X-ray data (see Supplementary Methods and Supplementary Figs. 4 and 5).

Magnetic susceptibility
Magnetic susceptibility was measured using a Quantum Design MPMS SQUID system. The measured susceptibility is shown in Supplementary Fig.  1c. As a basic characterization to allow comparison with previous reports, the data below T nC were fitted to the form χ = χ 0 + C/(T + θ). Parameters for the sample in this study are compared with those of some previous studies in Supplementary Table 1. It can be seen that the sample used in the present study shows the largest Pauli-like term and the smallest impurity term, corresponding to a weakly interacting spin ½ concentration of 0.02%. A close inspection of supplementary Fig. 1c shows an additional weak feature around 50 K where there is a small drop in χ. There are no comparable features in any of our other measurements and we assign this drop to a modification of the impurity term. This may be related to a change in the magnetic environment of the impurity, e.g. we can speculate that it reflects an impurity in a host diamagnetic layer that becomes a paramagnetic layer above 50 K.

Electronic transport measurements
In order to confirm the hysteretic behavior related to the commensuration of the CDW and compare it with previous results in the literature, large crystals of 1T-TaS 2 were contacted in an in-plane (ab plane) 4-probe configuration using platinum wires. The transport measurements were carried out in a Quantum Design PPMS-9 with a cooling/warming rate of 1 K/min. The resistivity values as well as the thermal dependence ( Supplementary Fig. 1d) are in agreement with the previous results in  Fig. 7 Layer stacking and the muon relaxation signal. a Regular ACACAC stacking of the triangular-lattice layers of Star-of-David (SD) clusters produces completely spin-paired bilayers that do not support the QSL state. Cleavage at C planes (dashed lines) in such a structure results in bilayer surfaces. Muon sites in this case all see diamagnetic environments (shown as D here). b A defect in the SD stacking sequence at the level of one in seven layers on average gives a sequence such as ACA⊥⊥AC for which one in four cleaved surfaces are unpaired monolayers that can support the QSL state. This ratio matches that seen in a recent STM study 24 . The muon stopping sites in this case split between diamagnetic (D) and paramagnetic (P) environments, with f para , the P to D site ratio, being 2/7 for this concentration of stacking defects. c, d Example sequences with greater concentrations of the electronic stacking defects giving more unpaired QSL layers and larger values of f para . the literature 18,49 , showing the formation of the commensurate charge density wave superlattice at around 200 K.

Muon measurements
The μSR experiments were carried out on the HiFi instrument at the ISIS Neutron and Muon Source. This muon instrument is optimized for LF measurements ranging from the mT region up to 5 T. The measurements reported here focus on the range 2 to 125 mT. A mosaic sample was made from several large crystals of 1T-TaS 2 (0.7272 g in total) all oriented with the ab plane parallel to the surface (Supplementary Fig. 1b). For the measurements, the mosaic was mounted on a silver sample plate and covered with a thin silver foil. This was placed in a closed cycle refrigerator and a beam of positive surface muons was implanted into the sample. The muon probe is fully spin polarized on implantation and the forward/ backward asymmetry of the detected muon decay positrons a(t) reflects the time-dependent polarization of the muon spin. In the LF measurements, the magnetic field was applied perpendicular to the ab plane of the crystals. Detailed measurements of the field and temperature dependence of the μSR were made in this study, counting 5 × 10 7 muon decay events in each μSR spectrum to get sufficient accuracy in the fitted parameters to estimate critical exponents that can be compared with those of QSL models. Data analysis was carried out using the WiMDA program 50 . The key advantage of the LF measurements over ZF studies is that both the diamagnetic contribution and the contribution from localized fluctuations are quenched with field as B LF −2 , leaving a relaxation at higher fields that is entirely due to the spin diffusion, which has a weaker field dependence. The LF relaxation can then be fitted to a single slow relaxing component. Further details of the relaxation model used for the analysis are given in the Supplementary Methods.

Specific heat
Specific heat was measured using a Quantum Design PPMS. The results obtained over a wide temperature range are shown in Fig. 4a. For the data analysis, a good representation of the lattice contribution is required. For the low-temperature region I a simple Debye T 3 dependence is appropriate where the only parameter is the Debye temperature θ D . At higher temperatures, the data become affected by the anisotropy of the acoustic phonons and the additional excitation of optical phonons and a more elaborate model is required. In regions II and III the phonon term is modeled by the sum of an anisotropic Debye term 51 representing the acoustic phonons and two Einstein terms to represent optical phonons. In this case C p lat /R = 3 F D (θ D , r) + 4 F E (θ E,1 ) + 2 F E (θ E,2 ), where F D is the anisotropic version of the Debye function 51 for anisotropy r and F E is a standard Einstein specific heat function. The three additional phonon parameters r, θ E,1 , and θ E,2 were determined by fitting the data in regions II and III between 25 and 180 K. Without imposing some additional constraint on the fitting, it is not possible to determine simultaneously the absolute values of n C in regions II and III, although a difference of order 1 between them is found for fits across a range of assumed values for n C in region II. In order to establish a consistent set of parameters between the LF-μSR and specific heat measurements, we therefore constrained the value of n C in region II to be the value corresponding to the experimental n D value found in region II, i.e. n C = 2.19 (15). The overall errors of the fitted parameters take into account the additional uncertainty associated with this constraint, as well as the statistical error of the fitting itself. The procedure works well in describing the data (Fig. 4c) and allows the clear drop in n C on going from region II to region III to be quantified. Parameters obtained for the phonon terms are listed in Table 1, together with an overall comparison of the experimental parameters obtained in this study.

Muon sites and hyperfine coupling
In order to identify the muon sites and their hyperfine coupling to the unpaired spin, calculations were made with plane-wave DFT using the CASTEP program 52 with the generalized gradient approximation and the PBE functional 53 . Calculations were carried out on a √13 × √13 × 2 supercell comprising two C-CDW unit cells stacked along the c-axis (the doubling of the c-axis is necessary to ensure sufficient isolation of the implanted muon from its periodic images). We use a plane-wave basis set cut off energy of 500 eV and a 5 × 5 × 5 Monkhorst-Pack 54 grid for k-point sampling, resulting in total energies that converge to an accuracy of 0.01 eV per cell. Muons, represented by a hydrogen pseudopotential, were initialized in 59 random positions within this supercell, generated by requiring initial positions being at least 0.25 Å away each other and from any of the atoms in the cell. Furthers details are given in the Supplementary Methods and Supplementary Fig. 2.