Vanishing Wilson ratio as the hallmark of quantum spin-liquid models

P. Prelovšek, 2 K. Morita, T. Tohyama, and J. Herbrych Jožef Stefan Institute, SI-1000 Ljubljana, Slovenia Faculty of Mathematics and Physics, University of Ljubljana, SI-1000 Ljubljana, Slovenia Department of Applied Physics, Tokyo University of Science, Tokyo 125-8585, Japan Department of Theoretical Physics, Faculty of Fundamental Problems of Technology, Wroclaw University of Science and Technology, 50-370 Wroclaw, Poland


I. INTRODUCTION
Various frustrated S = 1/2 Heisenberg models (HM) have been subject of intensive theoretical studies in last decades in connection with the possibility of spin-liquid (SL) ground state (g.s.). These efforts have been recently strengthened by the discovery of several classes of insulating materials revealing low-energy spin excitations behaving as quantum SL without any magnetic order down to low temperatures (for reviews see [1][2][3]). Among isotropic S = 1/2 two-dimensional (2D) models most numerical evidence for the SL g.s. accumulated for the antiferromagnetic (AFM) HM on the kagome lattice (KL) [4][5][6][7][8][9][10], but as well for J 1 -J 2 HM on the square lattice(SQL) [11][12][13][14][15][16][17][18], J 1 -J 2 HM on the triangular lattice (TL) [19][20][21][22][23][24] and the HM on TL with ring exchange [25,26]. While the character of the g.s. and its properties still offer several controversies and challenges, it is even less known about finite-temperature T > 0 behavior of several basic quantities. At least some of them have been already measured in experiments on SL materials and can thus serve as a test whether and to what extent actual materials can be accounted for by theoretical models.
Among measurable spin properties are thermodynamic quantities as the uniform magnetic susceptibility χ 0 (T ), magnetic (contribution to) specific heat C V (T ) and related spin entropy density s(T ). They are crucial to pinpoint the different characters and scenarios of SL behavior, in particular whether materials follow gapped or gapless SL. These quantities are mostly extracted from experiments on KL systems, the prominent example being herbertsmithite [27][28][29][30][31], but also related compounds in the same class [32][33][34][35][36][37]. Another example are organic compounds where the relevant lattice is triangular [38][39][40][41] as well as charge-density-wave system 1T-TaS 2 , recently established as SL with composite S = 1/2 spins on TL [42][43][44][45]. The basic spin exchange scale in most of these systems is modest and, as a consequence, the whole T range is experimentally accessible which allows for the test of the whole range of spin excitations. Nevertheless, it should be noted that the lowest T might be influenced by additional mechanisms such as Dzyaloshinski-Moriya interaction [46][47][48], interlayer couplings and random effects [49].
It has been rather well established with elaborate exact-diagonalization (ED) and series-expansion studies of the HM with nearest-neighbor (n.n.) exchange on KL [4-6, 8, 50, 51], that lowest excitations are singlets dominating over the triplet excitations, for which most ED studies reveal a finite spin (triplet) gap ∆ t > 0 although there are numerical indications also for gapless scenario [52,53]. It has been recently shown [54] that the same scenario can be traced via the temperature-dependent Wilson ratio R(T ) in a J 1 -J 2 HM on TL including the next-nearest-neighbor (n.n.n.) exchange J 2 > 0 in the regime where the SL g.s. is expected [19][20][21][22]. This is in contrast with the triplet (S = 1) magnon excitation being the lowest excitations in an ordered AFM. It is also qualitatively different from the scenario for the basic onedimensional (1D) HM with gapless spinon excitations.
In the following we will present numerical results for s(T ), χ 0 (T ) and R(T ), which reveal that the vanishing of R(T → 0) is quite generic property of a wide class of isotropic 2D Heisenberg models in their range of (presumable) SL parameter regimes. In this context we generalise previous numerical T > 0 studies of HM on KL [55,56] to include also the n.n.n. exchange J 2 = 0 and upgrade results for the J 1 -J 2 HM on TL [24], now studying also the HM on TL with the ring exchange, as well as another standard model of SL, i.e., frustrated J 1 -J 2 HM on SQL. Results in the SL regimes confirm the singlets as dominating low-energy excitations. For comparison we present results also for 1D J 1 -J 2 Heisenberg chain, which serve as the reference, depending on J 2 /J 1 , either for the gapless spinon Fermi-surface (SFS) and valencebond (VB) solid scenarios. Still, we show that results appear (as expected) qualitatively different from considered 2D models.
Considered models have their particular features and challenges, nevetheless our results on thermodynamic properties reveal quite universal properties in their (presumable) SL regimes which put also restrictions on the SL scenarios explaining their low-T behavior. In particular, very attractive scenario of gapless SL with SFS excitations requires finite g.s. Wilson ratio R 0 = R(T → 0) > 0. The latter is realized in 1D HM, but does not seem to be the case in planar models. Observed enhanced low-T entropy s(T ) and related vanishing R 0 = 0 demonstrate the dominant role of singlet excitations over the triplet ones [50,57], but still offer several possibilities. While it is hard to exclude the scenario of VB solid (crystal) with broken translational symmetry [6,50], it is more likely that the g.s. in the SL regime does not break the translation symmetry and all correlations are short-ranged, i.e. revealing a scenario of VB (or dimer) liquid. On the other hand, it is well possible that considered models might not be enough to represent the SL real materials, in particular not in their low-T regime.
The paper is organised as follows: In Sec. II we introduce T -dependent Wilson ratio R(T ) and comment different scenarios for its low-T behavior. In Sec. III we present numerical methods used to evaluate thermodynamic quantities, but also the lowest spin excitations in 1D and 2D models. As the test of methods as well as of concepts we present in Sec. IV results for 1D J 1 -J 2 Heisenberg chain. The central results for various 2D frustrated HM models are presented and analysed in Sec. V, and summarized in Sec. VI.

II. TEMPERATURE-DEPENDENT WILSON RATIO
Besides thermodynamic quantities: uniform magnetic susceptibility χ 0 (T ) and the entropy density s(T ), together with related specific heat C V (T ) = T ds/dT , it is informative to extract also their quotient in the form of temperature-dependent Wilson ratio R(T ), defined as [54], being dimensionless assuming theoretical units k B = gµ B = 1. It should be reminded that the standard quantity is the (zero-temperature) Wilson ratio as R W = 4π 2 χ 0 0 /(3γ) where χ 0 0 = χ 0 (T = 0) and γ = lim T →0 [C V /T ]. R W has its usual application and meaning in the theory of Fermi liquids and metals, as well as in gapless spin systems [58]. We note that in normal Fermi-liquid-like systems where s = C V = γT the definition, Eq. (1), coincides at T → 0 with the standard R W . Although at low T (in most interesting cases) both s(T ) and C V (T ) have the same functional T dependence, it is more convenient to employ in Eq. (1) the entropy density s(T ) being monotonously increasing function.
It should be also pointed out that R(T ) is a direct measure of the ratio of the density of excitations with finite z component of total spin S z tot = 0 relative to density of all (spin) excitations, including S z tot = 0, as measured by s(T ). To make this point evident we note that χ 0 (T ) = (S tot z ) 2 /(N T ) where N is the number of lattice sites, so that From above expression it is also follows that R(T ) has a well defined high-T limit which is for isotropic S = 1/2 HM R(T → ∞) = π 2 /(3 ln 2) = 4.746. Moreover, R 0 ≡ R(T → 0) can differentiate between distinct scenarios: a) In the case of magnetic long-range order (LRO), e.g., for AFM in HM on SQL, at T → 0 one expects in 2D isotropic HM χ 0 (T → 0) = χ 0 0 > 0 (where the finite value can be interpreted as the contribution of spin fluctuations transverse to g.s. magnetic order) whereas the effective magnon excitations lead to s ∝ T 2 [59], so that R 0 ∝ 1/T → ∞, b) In a gapless SL with large SFS one would expect Fermi-liquid-like finite R 0 ∼ 1 [2, 41,44]. The evident case for such scenario, as for reference considered later on, is the simple Heisenberg chain where R 0 = 2 [60], in contrast to the value R 0 = 1 for noninteracting Fermi systems. c) Vanishing R 0 → 0, or more restricted from Eq. (1) R 0 ∝ T η with η ≥ 1, would indicate that low-energy singlet excitation dominate over the triplet ones [2, 50,51]. In the following we find numerical evidence that this appears to be the case in the SL parameter regime of considered 2D frustrated isotropic HM.
Within the last scenario one should still differentiate different possibilities with respect to gapless spin systems or systems with the gap. One option for SL is that both singlet and triplet excitations are gapped, but the effective triplet gap is larger ∆ t > ∆ s (in the limit of large systems N → ∞) which would lead (in a simplest approximation) to R 0 ∝ T η exp[−(∆ t − ∆ s )/T ] → 0. More delicate case could be when ∆ t = ∆ s = ∆. Then Eq. (1) offers several scenarios with, e.g., R(T < ∆) ∝ T η . Such situation appears, e.g., for 1D chain J 1 -J 2 model around the Mazumdar-Ghosh point J 2 /J 1 = 0.5. Since s(T ) measures both singlet and triplet excitations (as well as higher S tot > 1) possible case ∆ s > ∆ t should be similar to the previous scenario.
When classifying options for T → 0 we should also consider the possibility of VB solid (crystal), i.e., the g.s. with broken translational symmetry. In finite systems (with short-range spin correlations) the signature of VB solid should be the degenerate or (due to finite-size effects) nearly degenerate g.s. with degeneracy N d > 1. This should be reflected in a finite g.s. entropy for finite system with N sites, Such remanent s 0 > 0 does not contribute to C V (T ) and moreover vanishes in the limit N → ∞. A clear VB solid case is 1D J 1 -J 2 HM in the dimerized regime where N d = 2. It then makes sense to consider in the evaluation of R(T ), Eq. (1), besides full s(T ) also reduced ones = s − s 0 . Still, it is not always straightforward to fix proper N d in finite-size systems.

III. METHODS
We calculate entropy density s(T ), uniform susceptibility χ 0 (T ) and via Eq. (1) the Wilson ratio R(T ), using the finite-temperature Lanczos method (FTLM) [61,62], previously used in numerous studies of T > 0 static and dynamical properties in various models of correlated electrons [63]. Since in the case of considered thermodynamic quantities only conserved quantities are involved, in particular the Hamiltonian H and S z tot , the memory and CPU time requirement for given system size N are essentially that of the Lanczos procedure for the g.s., provided that we scan over all (different) symmetry sectors S z tot and wave-vector q due to translational symmetry and periodic-boundary conditions (p.b.c.), in case of the code with translational symmetry. A modest additional sampling N s over initial wave-functions is then used. Limitations of the present method are given by the size of the many-body Hilbert space with N st basis states which can be handled efficiently within the FTLM, restricting in our study lattice sizes to N ≤ 36. In the following we use two FTLM codes for considered models: a) To calculate largest systems with N = 36 sites for the 2D TL, KL as well as SQL J 1 -J 2 HM with N st ∼ 10 10 , we develop a code that equips a technique to save the memory for the Hamiltonian by dividing the H into two subsystems. In addition, to improve the accuracy, we use replaced FTLM technique [64]. b) The code for more modest computers takes into account translational symmetry, thus able to reach N st < 10 7 and sizes N ≤ 30, was used for the 1D HM chain and the TL with ring exchange.
When discussing the accuracy of FTLM results we have to distinguish results for given system from finite-size effects due to restricted N . The central quantity evaluated is the grand-canonical sum [61,62], where E 0 is the g.s. energy. For reachable systems FTLM provides accurate results provided that we use modest random sampling N s ≤ 30 over (random) initial wavefunctions. This is in particular important to get correct low-T limit, i.e., Z(T = 0) = 1 in the case of nondegenerate g.s. The main restriction of FTLM results are, however, reachable N and related finite-size effects most pronounced at T → 0: a) In isotropic HM with T → 0 LRO (in dimension D ≥ 2), or long-range spin correlations in 1D, spin excitations are gapless in the thermodynamic limit. Such case is correlated with finite-size effects in evaluated quantities. One can expect that results reach the N → ∞ validity only for Z > Z * = Z(T f s ) 1. Since Z is intimately related to entropy, the criterion for T f s can be the smallest value for s. Actually, in reached systems N ∼ 36 we get estimate s(T f s ) ∼ 0.07 − 0.1 (see, e.g., the finite-size analysis in [56]). In such systems s(T ) and χ 0 (T ) results at T < T f s are dominated by finite-size effects and are not representative for N → ∞. In any case, due to frustration and consequently enhanced s(T J 1 ) in SL models FTLM generally allows to reach lower effective T f s . E.g., while for HM on a (unfrustrated) SQL (even at largest N = 36) T f s ∼ 0.4J 1 [62,63], SL models allow for considerably lower T f s ≤ 0.1 [24,56]. b) For systems with only short-range spin correlations one can reach situation where spin correlation length (even at T → 0) is shorter that the system length, ξ ≤ L. In such a case, FTLM has no obvious restrictions even at T → 0, so T f s ∼ 0. This can be the situation for gapped SL, including some examples discussed further on.
Besides thermodynamic quantities, it is also instructive to monitor directly lowest excited states and their character. For the largest 2D N = 36 systems excited states are obtained within ED (without translational symmetry) by eliminating Lanczos-ghost states while comparing results for different number of Lanczos steps. For TL with ring exchange we employ ED results of systems with N = 28 and evaluate the lowest (singlet and triplet) energies in different q sectors.
In 1D models we use also density matrix renormalization group (DMRG) method to investigate the J 1 -J 2 HM with open boundary conditions (o.b.c.). The method allows for accurate computation of the S z tot = 0 g.s., and in the same way also first excited triplet state with S z tot = 1. In order to get also excited (singlet) states within S z tot = 0 sector we evaluate the g.s. wave-function |ψ 0 and then construct effective Hamiltonian for the excited states 65], and then repeat the standard DMRG algorithm for H 1 . The requirement of orthogonality is, however, difficult to meet for excited states which are (due to o.b.c.) edge states, e.g., as within 1D dimerized regime.

IV. ONE-DIMENSIONAL HEISENBERG MODEL
We consider first the 1D J 1 -J 2 HM, which can serve as the reference for further discussion of 2D HM results. The AFM isotropic S = 1/2 J 1 -J 2 HM is given by where we further on put J 1 = J = 1 as the unit of energy. We will investigate with FTLM only J 2 ≥ 0 case on systems of finite length N with p.b.c. Thermodynamic properties are well known and understood for simple J 2 = 0 Heisenberg chain [60], as well the g.s. and the triplet excited state for the frustrated chain with J 2 > 0 [66]. Beyond critical J 2 > J * 2 ∼ 0.241 the g.s. is dimerized (N d = 2) in the thermodynamic limit [66]. At the same time, lowest excited states are degenerate triplets and singlets with the gap ∆ t = ∆ s , consistent with the unbound spinons as elementary excitations.
Numerical results for s(T ), χ 0 (T ) and finally R(T ), as obtained on a system with N = 30 sites, are presented in Fig. 1 for different 0 ≤ J 2 ≤ 1: a) For the simple J 2 = 0 chain we get s(T ) ∼ γT in very broad range T < 0.6. Finite-size effects are most pronounced in this case, so that below T < T f s ∼ 0.2 we get s < 0.1 and finite-size effects prevent any further firm conclusions. Still, for T > T f s numerical results are consistent with analytical and previous numerical results, in particular with the known limit R 0 = 2 [60]. Moreover, it is remarkable that R(T ) is nearly constant in a wide range T < 0.6. b) The gap becomes pronounced for the Mazumdar-Ghosh point J 2 = 0.5 and even more for J 2 = 1.0 (where ∆ t ∼ 0.25 [66]). In the gapped case FTLM finite-size effects are less pronounced, and one can expect T f s → 0. In fact, for the J 2 = 0.5 and J 2 = 1.0 results appear sizeindependent for reached N = 30, apart from the dimerization degeneracy N d = 2 leading via Eq. (3) to s 0 > 0. The latter has influence on the R(T ∼ 0), so we present in Fig. 1 also the result taking into account subtracted s(T ). In both analyses the behavior is consistent with R 0 = 0. For the J 2 = 0.5 and J 2 = 1.0 modified results are still consistent with vanishing R(T < ∆) ∝ T η with η ≥ 1, but this behavior remains to be clarified. For the marginal case J 2 = 0.3 ∼ J * 2 , the behavior of all quantities is similar to J 2 = 0, except that we find γ larger and consequently also T f s smaller.
It is instructive to investigate in connection with finitesize effects also lowest triplet and singlet excitations in the model. While triplet excitations have been in detail studied using DMRG already in Ref. [66], to establish singlet excitations requires more care, see Sec. III. In Fig. 2 (a) we present the DMRG (with o.b.c.) N = 60 result for excitations: lowest triplet t and lowest singlet s vs. J 2 , together (as the inset) with their 1/N scaling in the gapless regime J 2 = 0.2 < J * 2 . Due to o.b.c. DMRG is unable to properly resolve the dimerized partner of g.s. since it represents in open chain excited edge states. Hence, we present in Fig. 2(a) the first singlet excited state only for J 2 ≤ 0.4. Still, DMRG results confirm that no other singlet is stable below the triplet for J 2 > J * 2 , unlike seen later on in 2D SL models. In Fig. 2(b) we display also DMRG results for g.s. bond spin correlations S z i S z j . It is also apparent that for J 2 > J * 2 the g.s. is dimerized (in n.n. bond correlations), whereby the particular case is J 2 = 0.5 with alternating n.n. correlations S z i S z j = −1/4 and 0. Stronger correlations remain AFM in the whole J 2 > J * 2 , while it is easy to recognize the change of character from AFM for J * c < J 2 < 0.5, to a ferromagnetic state for J 2 > 0.5.

V. PLANAR FRUSTRATED HEISENBERG MODELS
A. J1-J2 Heisenberg model on kagome lattice HM on KL is the prototype model for the existence of SL in planar models. It has been the subject of numerous studies, devoted mostly to the g.s. using ED [4,5,8,51], series expansion [6], DMRG [7,10,67] and variational methods [9,52]. We consider here the extended model with p.b.c., involving also the n.n.n. exchange J 2 as shown in the inset of Fig. 3(c), whereby the role of J 2 > 0, as well as J 2 < 0, is to reestablish the magnetic LRO [68]. The basic HM on KL has been the clearest case for a dominant role of lowlying singlet excitations over the triplet ones [6,50,51]. The latter fact and related large entropy, persistent at low T 1, has been well captured within block-spin [4,5,69] and recently within related reduced-basis approach [54], whereby the singlet excitations can be attributed to chiral fluctuations, distinct from (higher-energy) triplet excitations.
Thermodynamic quantities for the basic J 2 = 0 HM on KL have been calculated via FTLM previously [56] up to the size N = 42. Here we extend the study, evaluating via FTLM also for J 2 = 0 for N = 36. Results in Fig. 3 reveal that increasing |J 2 | > 0 suppreses strongly s(T J 1 )) while leaving χ 0 (T ) less affected (at least for T > T f s ). The result for J 2 = ±0.2 indicates on divergent R 0 → ∞ consistent with the emergent magnetic LRO [54,68]. On the other hand, at J 2 = 0, 0.1 the behavior of χ 0 (T ) and s(T ) are consistent with finite triplet gap ∆ t ∼ 0.15 and smaller or even vanishing singlet gap ∆ s < ∆ t .
The transition from the singlet-dominated SL regime to the phases with magnetic LRO one can monitor also via low-lying levels in considered systems. In Fig. 4 we present the evolution of excitation energies for lowest lying triplet t as well as several low-lying excited singlets s,i (i = 1, 2, · · · , 6), as obtained via ED on N = 36 sites, and in part for N = 32 sites for the HM on TL with ring exchange. It should be pointed out that we monitor only nondegenerate excited states, whereby in general the degeneracy is present and depends on particular lattice and related p.b.c. The level evolution, plotted vs. J 2 (or J r discussed lateron) serves primarily as another test where to expect SL with macroscopic number (in the limit N → ∞) of singlet excitation below the triplet ones, but also to locate some possible g.s. first-order transitions. In Fig. 4(a) the level scheme for KL is consistent with the previous ED studies of (J 2 = 0) KL model [6,50,51] which reveal a massive density of singlet levels with s ∼ 0 below the lowest triplet one t . Introducing |J 2 | > 0 reduces the degeneracy and might lead to ∆ s > 0 even in the N → ∞ limit. Still, a large density of singlet levels appear below the triplet in a wide (SL) range J c1 2 < J 2 < J c2 2 where J c1 2 ∼ −0.1, J c2 c ∼ 0.1 from Fig. 4(a) and we define J c1,c2 2 with the crossing of (all) lowest s,1−6 < t . We note that marginal J c1,c2 2 are consistent with Fig. 3 where J 2 = ±0.2 already reveal magnetic LRO with R 0 → ∞.

B. J1-J2 Heisenberg model on triangular lattice
While numerical studies for the basic (J 2 = 0) HM on TL [70][71][72] confirm magnetic LRO with moments pointing into 120 0 -angle directions, modest additional frustration with J 2 > 0 allows for the possibility of SL g.s., with the evidence for either gapless [19] or gapped SL [20][21][22][23] in the intermediate regime J 2 ∼ 0.15. Beyond that, for J 2 > 0.2 again stripe AFM is expected. Thermodynamic (and some dynamic) quantities for the J 1 -J 2 HM, Eq. (7), on TL have been recently calculated using FTLM [24] up to N = 30 sites and employing the reduced-basis approach [54], whereby the similarity of Here we upgrade previous FTLM studies with the calculation of J 1 -J 2 HM on TL to N = 36 sites. Results, presented in Fig. 5, are qualitatively consistent with previous ones for N = 30 [24] but due to larger size and consequently smaller T f s more clearly reveal the small entropy s(T ) and related diverging R(T ) below T ∼ 0.2 for J 2 ∼ 0, where the g.s. possesses magnetic LRO. A similar behavior is expected for J 2 > 0.2 where the stripe AFM g.s. has been established [19]. In reachable system N = 36 the upturn of R(T ) is partly masked by finite-size s 0 > 0, Eq. (3), due to the degeneracy N d > 1 of striped magnetic LRO, evident in Fig. 5(a)   χ 0 (T → 0) (indicating a finite triplet gap ∆ t > 0) leads to vanishing R 0 = 0. In Fig. 4(b) we plot the corresponding evolution of excitations vs. J 2 , as obtained with ED on N = 36 lattice. The triplet gap apparently remains substantial, i.e. t > 0.38 for considered N in the whole range of J 2 < 0.3. Still, singlet excitations s,1−6 all cross t for small J 2 ∼ 0.1. This leads effectively to a g.s. level crossing s,1 = 0 at J 2 ∼ 0.17 exchanging the character of the g.s. into a striped AFM. But most important, in the intermediate range 0.1 < J 2 < 0.17, which should be the relevant SL regime, singlet-excitations collapse is consistent with the conclusions from thermodynamics in Fig. 5 and R 0 = 0. It should be, however, acknowledged that the singlet collapse is not as pronounced as for J 2 ∼ 0 KL in Fig. 4(a).

C. Heisenberg model with ring exchange on triangular lattice
While J 1 -J 2 HM on TL is conceptually simple, it is less obvious to justify in connection with experiments and with more basic models. The organic SL materials [38][39][40][41] and 1T-TaS 2 [42][43][44][45] are closer to the metal-insulator transition where simple S = 1/2 n.n. HM is presumably not enough. Assuming as the starting point the singleband Hubbard model on the insulator side of the Mott transition U > U c the lowest correction to the n.n. HM comes in the form of the ring exchange term [25,26,73,74], with where ijkl are taken over different four-cycles on TL, as shown in the inset of Fig. 6(c). H r , Eq. (9), has been confirmed as the leading correction in the numerical study of the half-filled Hubbard model [73] in the insulating regime where J r ∼ 80t 4 /U 3 ∼ 20t 2 /U 2 J < 0.2 [74], taking into account that the Mott insulator on TL requires U > U c ∼ 8t−10t and J ∼ 4t 2 /U . It should be also mentioned that in Eq. (9) we do not consider higher-order t/U corrections to J. It has been already proposed that modest ring exchange J r > 0 on TL destroys the magnetic LRO and induces SL g.s. [25,26], including the observation of several possible singlet excitations below the lowest triplet one. In Fig. 6 we present results for the HM on TL with J r > 0 , Eq. (9), as obtained via FTLM on N = 28 sites (smaller size due to more complex H). It is evident that J r > 0 steadily increases low-T entropy s(T ), while increasing χ 0 (T ). Resulting R(T ) looses magnetic LRO character already for J r ≥ 0.05 being followed by SL-like regime with vanishing R 0 → 0.
The same message follows from the consideration of lowest levels on N = 28 lattice, presented in Fig. 4(d). Analogous to Figs. 4(a) and 4(b), there is a clear collapse of singlet levels s,1−3 (here we employ a q-resolved code and cannot monitor all singlet excitations) below the triplet one t for J r > 0.1. In the latter regime t represents already a reasonable estimate of the limiting N → ∞ triplet gap ∆ t > 0 [25], whereas to establish a proper singlet gap (the lowest singlet in N → ∞ limit) ∆ s < ∆ t requires more detailed finite-size analysis.
D. J1-J2 Heisenberg model on square lattice Finally, we turn to the J 1 -J 2 HM, Eq. (7) on SQL. The latter has been one of first considered for the possible (plaquette) VB solid g.s. at intermediate J 2 ∼ 0.5 [11,12,16,75], but also for the SL g.s. [13-15, 17, 18]. Results for corresponding thermodynamic quantities presented in Fig. 7 are consistent with the diverging R 0 → ∞ indicating magnetic LRO outside quite narrow parameter regime, i.e. outside 0.5 ≤ J 2 ≤ 0.6. In the latter regime we again find substantial entropy s(T 1) and consequently R 0 → 0, whereby for J 2 ∼ 0.6 there are already some indications for possible degeneracy s 0 > 0 which could be in favor of broken translational symmetry, e.g., a plaquette VB solid [11,12,16,75].
Caveats for the SL interpretation emerge also when considering the excitation evolution vs. J 2 [see Fig. 4(c)], as obtained from ED results on N = 36 cluster. For given system size, the singlet levels reveal s,1−6 < t only in a very narrow regime 0.55 < J 2 < 0.62. Even then, higher singlets (apart from s,1 ) are not well below t . Consistent with previous works [11][12][13][14][15][16][17][18] level scheme indicates on a change of the g.s. character for J 2 > 0.6. As a consequence, the SL in the intermediate regime, and even more on the singlet-dominated regime is less conclusive, and other options [17,75] have to be also considered.

VI. CONCLUSIONS
Thermodynamic quantities: entropy density s(T ) (together with directly related specific heat C V (T ) = T ds/dT , not presented in this paper), uniform susceptibility χ 0 (T ), and consequently T -dependent Wilson ratio R(T ), offer another view on properties of frustrated spin models. We considered here prototype 2D isotropic S = 1/2 HM, which are at least in some parameter regimes best candidates for the SL g.s. For comparison, we investigated in the same manner also simplest 1D HM which can serve as reference for some concepts and scenarios.
R(T ), in particular its low-T variation, is the quantity which differentiates between different scenarios. Whereas 2D systems with magnetic LRO can be monitored via R 0 → ∞, we are more interested in the SL regimes with g.s. without magnetic LRO and even without any broken translational symmetry which could be classified as VB solid (crystal). As prototype case we present results for 1D J 1 -J 2 HM which does not have magnetic LRO, but offers already two firm scenarios: a) the gapless regime for J 2 < J * 2 with spinons (or 1D SFS) as elementary excitations, and consequently finite R 0 = R(T → 0) ∼ 2 (for J 2 ∼ 0), b) a gapped regime for J 2 > J * 2 with dimerized g.s. (being the simplest 1D form of VB solid) apparently also with R 0 = 0, although not yet fully resolved variation R(T → 0). SL regimes in considered 2D frustrated isotropic S = 1/2 HM are in our study located via enhanced low-T entropy s(T ) and gapped (or at least reduced) χ 0 (T ), resulting in vanishing R 0 = 0. Similar information and criterion (although less well defined) emerges from the excitation spectra, when differentiating singlet and triplet (or even higher S tot > 1) excitations over the S tot = 0 g.s. Most evident cases for such VB (dimer) liquid scenario appears within the KL around J 2 ∼ 0. Analogous, although somewhat less pronounced, case is obtained within HM on TL with ring exchange J r > 0.1 and for the J 1 -J 2 HM on TL in the intermediate regime 0.1 < J 2 < 0.17. For such systems the level evolution as well as R(T ) reveal massive density of singlet states below the lowest triplet excitation. On the other hand, the situation in the HM on SQL in the narrow regime J 2 ∼ 0.6 is less clear-cut in this respect, since singlets are not well below the lowest triplet.
Vanishing R 0 does not support the scenario of SL with large (or even Dirac-cone) spinon Fermi surface, which would require finite R 0 > 0 (as in 1D HM), although our finite-size studies should be interpreted with care and cannot give a final answer to this problem. Still, emergent scenario of VB (dimer) liquid should be critically faced with the possibility of VB solid. The difference should be that in the case of VB solid g.s. should be (due to broken translational symmetry) degenerate with finite N d > 1 (in the thermodynamic limit N → ∞). Except for the SQL at J 2 ∼ 0.6, we do not find much evidence for that in the presumable SL regimes, since results mostly indicate (besides finite triplet gap ∆ t > 0) also either on finite singlet gap ∆ s > 0 or vanishing ∆ s ∼ 0 for the pure KL, but evidently ∆ s < ∆ t . To establish (or exclude) possible N d > 1 and to determine ∆ s > 0 beyond doubt still requires further studies.
Finally, it should be stressed that evaluated thermodynamic quantities are (at least in principle) measurable in experimental realizations of SL materials. s(T ) is accesible via measured magnetic specific heat C V (T ) and uniform susceptibility χ 0 (T ) via macroscopic d.c. or/and Knight-shift measurement. Since known SL materials are characterized by modest exchange J, properties can be measured in the wide range T J. This offers the possibility of critical comparison with model results, whereby considered isotropic HM might still miss some ingredients relevant for the low-T behavior, in particular the Dzyaloshiniskii-Moriya interaction, the disorder influence, and the inter-layer coupling. ACKNOWLEDGMENTS P.P. is supported by the program P1-0044 and project N1-0088 of the Slovenian Research Agency. K. M. and T. T. are supported by MEXT, Japan, as a social and scientific priority issue (creation of new functional devices and high-performance materials to support next-generation industries) to be tackled by using a post-K computer. T.T. is also supported by the JSPS KAKENHI (No. JP19H05825). The numerical calculation was partly carried out at the facilities of the Supercomputer Center, the Institute for Solid State Physics, the University of Tokyo and at the Yukawa Institute Computer Facility, Kyoto University. J. H. acknowledges grant support by the Polish National Agency of Academic Exchange (NAWA) under contract PPN/PPO/2018/1/00035.