Quantum phases and spin liquid properties of 1T-TaS₂

Abstract

Quantum materials exhibiting magnetic frustration are connected to diverse phenomena, including high T_c 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-TaS₂ 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-TaS₂ and their characteristic parameters can be interpreted as those of distinct Z₂ 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¹. 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,3,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)¹². This report concerns 1T-TaS₂, 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¹.

Fig. 1: Excitations in 1T-TaS₂ and the muon probe site.

a The basic hopping step for 2D spinon diffusion in a QSL involves rearrangement of the valence bonds (VB). b A typical interstitial muon probe site obtained using DFT.

The electronic and magnetic scenario offered by 1T-TaS₂ is far from being simple and is currently being revisited. Theoretical studies have shown that a single layer of 1T-TaS₂ is indeed a QSL¹³ and experimental results in the selenium counterpart, 1T-TaSe₂, have revealed that single layers are Mott insulators that behave as QSLs^14,15. However, when considering bulk crystals of 1T-TaS₂, 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₂ 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¹³. 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₂ 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₂. 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₂ 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 time-dependent 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 τ²⁸. μSR can also measure critical parameters⁶ 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.

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 χ² 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 χ² 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 χ² 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.

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²⁹. This method has been used in μSR for studying the propagation of spinons in spin-½ Heisenberg antiferromagnetic chain systems³⁰, 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²⁹. 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₀ ~ 110 K. The region below T₀ labeled II shows a relatively strong power law $D_{\mathrm{2D}} \propto T^{{n}_{\mathrm{D}}}$ with n_D = 1.74(14). Below T₁ ~ 25 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₀ 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₀ parameter also shows multi-region behavior, with region boundaries matching those of D_2D.

Fig. 3: Spinon diffusion obtained from LF-μSR in the C-CDW region below T_nC.

a Temperature dependence of the 2D spinon diffusion rate D_2D. b The scaling parameter for the amplitude of the diffusive signal (taking Ā = 4.5 MHz). Three distinct regions are found in the measurements, labeled I, II, and III. In each case the error bars indicate standard errors of the fit parameters.

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_{\mathrm{p}}^{\mathrm{e1}}} \propto T^{{n}_{\mathrm{C}}}$. 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₀, corresponding to the distinct peak found in the μSR parameters of Fig. 3. In contrast to the strong feature at T₀, the border between regions I and II at T₁ 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 ³ 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)³¹. 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₀, as shown in Fig. 4c. A consistent parameter set for the different regions is summarized in Table 1.

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 ³ 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.

Table 1 Summary of the properties of the three QSL regions obtained in this study.

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₂^3,13. Of these, the Z₂ states are expected to be most stable^4,13 and there are well-developed theories for describing the quantum critical (QC) properties of Z₂ QSL phases^32,33,34. Thus, Z₂ 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₂ QSL the spinon transport in the QC regime³⁴ 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 non-critical 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)³⁵. 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₂ QSL model with q = 0. Data for region II requires q = 1, suggesting a Z₂-linear QSL. In region III, it is found that q = 0 again gives the closest match to the n_C value.

Fig. 5: Critical parameters in the different regions.

a Values for n_D and n_C derived directly from free fitting of the data (solid circles). Values for n_C in regions II and III obtained by a fit constraining n_C in region II to be consistent with n_D in region II (open circles). The open square for region III shows the n_D value that would be consistent with the n_C obtained in the same region. The measured n_D in this region at −4(1) deviates strongly from this value. Dashed and dotted lines show O(4) model values for q = 0 and q = 1. b The self-consistent set of derived critical exponents compared with the O(4) model values, shown as dashed lines. Open and closed symbols indicate the q values found for each region. Error bars indicate standard errors of the fitted parameters.

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 inter-plane 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³⁶, but the predicted exponent n_D = −1/3 is very much smaller than measured here (Fig. 3a) and the corresponding specific heat prediction³¹ 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 AF ground state lies at a specific value of control parameter g₁ (Fig. 6a). A possible association of g₁ could be made with ring-exchange, which is known to destabilize the magnetic phase of triangular-lattice Mott-Hubbard insulators³⁹. The energies of the competing QSL states are assumed to depend on a second control parameter g₂ and at T = 0 continuous second-order QPTs would be expected between pairs of QSL-QD states (Fig. 6b).

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₂ spin liquid phase. The QPT is tuned by a control parameter g₁ 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₁ 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₂ is introduced to reflect the stabilization of different QSL phases. c A reported triangular lattice phase diagram from variational studies⁴⁰. The control parameters are the second-nearest-neighbor exchange J₂ 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₂–K mapping of the parameters g₁ and g₂. The purple arrow shows a possible parameter path for 1T-TaS₂ on cooling through regions II and I.

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 second-order 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₂ is taken to have a T dependence that produces smooth crossover transitions between the QSL states. The finite temperature transitions at T₀ and T₁ then would directly correspond to multicritical QCPs in a suggested underlying (g₁, g₂) quantum phase diagram for T = 0, as shown schematically in Fig. 6b.

Comparison with a QSL phase diagram

Spinon pairing instabilities leading to Z₂ spin liquid states have been discussed previously in several studies^40,41,42,43. In general, there are 64 possible symmetric Z₂ 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.⁴⁰. This uses a model Hamiltonian that includes terms representing the nearest-neighbor exchange J, the next-nearest-neighbor exchange J₂, and the ring-exchange K. The parameters J₂ 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.⁴⁰ that is most relevant to a Mott–Hubbard insulator, where the exchange parameters are determined by the U/t ratio. In this study a BCS-like d-wave paired Z₂ spin liquid was found over an extended region of the phase diagram (Fig. 6c)⁴⁰. The properties of this nodal state were discussed by Grover et al.⁴³. Below a characteristic pairing temperature it shows a linear (q = 1) DOS. In their study Mishmash et al.⁴⁰ 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₁ and g₂ control parameters to the K–J₂ 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⁴⁰ 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)₂Cu₂(CN)₃, 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₂ 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₂ remains poorly understood. Recent experiments address possible origins for the insulating state, from the conventional Mott scenario¹⁹ to layer dimerization⁴⁷, unit-cell doubling²⁴, or a Peierls transition²². Synchrotron X-ray studies reveal significant disorder in the interlayer stacking of the SD clusters at 100 K from the presence of diffuse scattering⁴⁷, suggesting a transition from a band insulator at low temperatures to a Mott insulator at higher ones²⁶. 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⁴⁷. A scanning tunneling microscopy (STM) study²⁴ 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²⁴ 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.

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²⁴. 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.

Alternative scenarios for the higher T regions

The f₀ 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₀ seen in regions II and III provide a greater challenge for the overall interpretation. The largest value of f₀ is found at T₀, 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₀ in region III above T₀ 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₀ on the lower side of T₀ 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₀ is entirely due to a reduction in f_para at low T. This could be consistent with the STM study²⁴, 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₁, reaching a maximum at T₀. This scenario would not however be consistent with the observed large change in n_D at T₁, 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₀ within region II is that f_para remains at or close to one, so that the changes in f₀ 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₂ within the C-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₂ 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₂ 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.

Methods

Crystal growth

High-quality 1T-TaS₂ 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⁻⁵ 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₂] = 1 mg/cm³), and sealed in an evacuated quartz ampoule (length: 500 mm, internal diameter: 14 mm; P: 2 × 10⁻⁵ 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 X-ray platform with Cu radiation source) to the previously reported structure of 1T-TaS₂ (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⁴⁸. 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 χ = χ₀ + 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₂ 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 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₂ (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⁷ 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⁵⁰. 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⁻², 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 ³ 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⁵¹ 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⁵¹ 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⁵² with the generalized gradient approximation and the PBE functional⁵³. 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⁵⁴ 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.