Finite-temperature crossover phenomenon in the $S=1/2$ antiferromagnetic Heisenberg model on the kagome lattice

Thermal properties of the $S=1/2$ kagome Heisenberg antiferromagnet at low temperatures are investigated by means of the Hams-de Raedt method for clusters of up to 36 sites possessing a full symmetry of the lattice. The specific heat exhibits, in addition to the double peaks, the third and the forth peaks at lower temperatures. With decreasing the temperature, the type of the magnetic short-range order (SRO) changes around the third-peak temperature from the $\sqrt{3} \times \sqrt{3}$ to the $q$=0 states, suggesting that the third peak of the specific heat is associated with a crossover phenomenon between the spin-liquid states with distinct magnetic SRO. Experimental implications are discussed.

Geometrically frustrated magnets have attracted special interest due to its unique and novel ordering properties. Among them, kagome antiferromagnets have long been studied extensively. Especially, much recent interest has been paid to the quantum spin-1/2 nearest-neighbor (n.n.) antiferromagnetic (AF) Heisenberg model on the kagome lattice because of the possible realization of a quantum spin-liquid (QSL) state having no magnetic long-range order. A large number of theoretical studies performed to understand the nature of its ground state have lead to various competing scenarios on the nature of its ground state, including the Z 2 spin liquid, 1-7 the algebraic U(1) spin liquid, [8][9][10][11][12] the chiral spin liquid, 13 the valence bond crystal, [14][15][16] etc. The true situation, however, still remains unclear.
Along with such intensive studies on the ground state, thermal properties of the S = 1/2 kagome antiferromagnetic Heisenberg (KAH) model at finite temperatures have also attracted much attention. A highlight issue might be the exotic temperature (T ) dependence of the specific heat, which exhibits multiple peaks. Namely, earlier numerical studies based on an exact diagonalization (ED) method 17,18 or a decoupledcell Monte Carlo (MC) simulation 18 indicated that, in addition to the broad peak at a higher T , the specific heat exhibited the second peak at a lower T . Whether this second peak identified for small-size systems really survives in the thermodynamic limit had been examined by various calculations: Mentioning some of them, the ED method up to 24 spins, 17-21 a high-T expansion, 19, 22 a high-T entropy method, 20 an approximate effective-Hamiltonian method, 23, 24 a transfer-matrix MC method 25 and a linked-cluster algorithm. 26,27 Recently, Sugiura and Shimizu have succeeded in computing the specific heat of the model up to the sizes N = 27 and 30, 28 by using the imaginary-time version of the equationof-motion method (the Hams-de Raedt algorithm), a powerful numerical technique of computing thermal properties of * t.shimokaw@gmail.com the quantum model at finite T developed some time ago 29,30 (Sugiura and Shimizu called the method the canonical thermal pure quantum state method). It was then observed that, on increasing the system size up to N = 30, the second peak was appreciably suppressed, with only a shoulder-like structure remaining. 28 Similar behavior was reported also by a finite-T Lanczos method 31 applied to the S = 1/2 KAH model of N ≤ 30. 32 In the present Letter, we wish to investigate the finite-T properties of the S = 1/2 KAH model by means of the Hamsde Raedt method, paying attention not only to the multiplepeak problem of the specific heat, but also to the type of the magnetic SRO as mentioned below. We extend the cluster size up to 36 spins possessing a full symmetry of the lattice, exceeding the previous works. While the method could provide us exact information for frustrated quantum systems, special care might be taken to specific circumstances peculiar to the KAH model. For example, possible significance of the √ 3 × √ 3 [see Fig. 1(a)] SRO calls for finite-size clusters of multiples of nine, while the possible singlet ground state calls for even-N clusters. In fact, the maximum sizes treated in recent exact finite-T calculations, N = 24 − 30, do not satisfy these requirements, and the results might be subject to stronger finite-size effects. [17][18][19][20]28 The maximum size treated in the present work N = 36 meets these criteria and is fa-vorable in that respect. Indeed, the 36-spin cluster possesses the full symmetry of infinite kagome lattice under periodic boundary conditions. We find that the second peak of the specific heat persists in the 36-spin cluster, even a bit more enhanced than that in the 27-and 30-spin clusters, suggesting that the second peak (or the shoulder) persists in the continuum limit. In addition, we observe in the 36-spin cluster the third and the fourth peaks at lower T . 21 Interestingly, the third peak turns out to be associated with a finite-temperature crossover phenomenon between the two distinct magnetic SRO states, i.e., the ones with the √ 3 × √ 3 SRO at higher T and the q = 0 [see Fig. 1 Our model is the S =1/2 AF Heisenberg model on the kagome lattice, whose Hamiltonian is given by is a spin-1/2 operator at the i-th site on the lattice, and J = 1 is the nearest-neighbor AF coupling. We treat several finite-size kagome clusters up to 36 spins with periodic boundary conditions. In computing the T dependence of various physical quantities, we employ the Hams-de Raedt method. [28][29][30] This method allows us to compute physical quantities by treating a small number of quantum states instead of taking an ensemble average over a full spectrum of the Hilbert space. While in the most direct ED method the memory limitation prevents us from treating more than 20 spins even for S =1/2, the method enables us to treat systems containing about 40 spins.
We briefly describe the computational method. A set of pure states for the inverse temperature β = 1/T and the system size N, |β, N , is constructed by operating an operator on a set of initial random vectors as |β, N = exp[−βH/2]|ψ 0 , where the initial vectors are given by |ψ 0 = 2 N i=1 c i |i with {c i } being random complex numbers satisfying the normalization condition i |c i | 2 = 1 and with {|i } being an arbitrary orthonormal basis set of the Hilbert space of H. Hams and de Raedt proved that the standard thermal average of a physical quantityÂ was given by where the overline denotes the average over the initial random vectors. 30 When this random average is performed over finite number I of realizations of initial vectors, as is necessarily the case in real numerical calculations, the deviation from the true value decays as ∼ 1/ √ ID where D represents the dimension of the entire Hilbert space of the model, D = 2 N in the present case. 30 This relation means that, for larger system size N, fairly accurate value can be obtained even from smaller-I calculations. In our computation, the average over initial random vectors is taken over 200 (N=18), 40 (N=24), 20 (N=27), 10 (N=30) and 3 (N=36) realizations.
The computed specific heat is displayed in Fig. 2 Our results for N = 18 and 24 agree with the ED results of Refs., [19][20][21] and those for N=27 and 30 with the results of Ref.. 28 A broad first peak arises around T ∼ 0.7 (not indicated in Figs.2), while the second peak arises at T ≃ 0.05 − 0.1, whose location gradually moves to lower T on increasing N. Interestingly, on increasing N for 18 ≤ N ≤ 30, the sharpness of the second peak tends to be gradually suppressed, while it is a bit more enhanced for N = 36 than for N = 30. This recovery of the second peak might be related to the fact that the N = 36 cluster retains a full symmetry of the infinite kagome lattice allowing for the √ 3 × √ 3 structure accommodated. As shown below, we find that the magnetic SRO around the second-peak temperature is indeed the √ 3 × √ 3 structure. In addition to the first and the second peaks, the third peak appears at a lower T around 0.01 T 0.02, although its sharpness is largely size-dependent. Furthermore, even the fourth peak appears at around the lowest T studied T ≃ 0.005 for N = 36, consistently with an earlier result by an approximate method. 24 For N = 36, the singlet gap was estimated to be ∼ 0.010, 34 while the triplet gap to be much larger ∼ 0.164. 34 Thus, the observed low-T structure of the specific heat is borne by singlet excitations. Our observation then sug- gests that at least the singlet gap would be quite small in the bulk, 0.01, or even to be gapless.
In the inset of the Fig. 2 in the T → 0 limit observing the third law, while, for N = 30, an additional specific-heat peak seems to be required at a still lower temperature of T/J < 0.005. In order to get information about the nature of the spin SRO, we compute the static spin structure factor S (q, β) defined by The computed S (q, β) for N = 36 are shown in Figs. 3 as an intensity plot in the wavevector q=(q x ,q y ) plane for (a) T =0, (b) 0.08, and (c) 0.01. At T = 0, S (q, β) exhibits a broad ridge-like structure along the zone boundary of the extended Brillouin zone (BZ), in which weak SRO peaks are observed at the wavevector points corresponding to the q = 0 state, i.e., at q = (0, ±2π/ √ 3) and (±π, ±π/ √ 3), the length unit here taken to be the n.n. distance of the original kagome lattice. The result is consistent with the earlier ED result by Läuchli. 35 Note that, although the ED calculation on smaller size clusters of N =18 and 27 favors the √ 3 × √ 3 structure than the q = 0 structure in contrast to the present result on N = 36 (refer also to Fig. 4 below), the recent DMRG calculation for larger systems up to N = 108 also reported the SRO peaks appearing at the points corresponding to the q = 0 state. 2 These observations lend support to the expectation that quantum fluctuations favors the q = 0 SRO than the √ 3 × √ 3 SRO in the ground state of the S = 1/2 model. If one recalls the fact that the semi-classical or large-S calculations suggest the preference of the √ 3 × √ 3 state, 36, 37 a natural expectation would be that the √ 3× √ 3 SRO is favored at higher T even in the S = 1/2 model. Then, a crossover associated with the change of the dominant type of SRO might occur at a certain finite T between the q = 0 SRO at lower T and the √ 3 × √ 3 one at higher T . Fig. 3(b) exhibits S (q, β) at a temperature T =0.08 close to the second peak of the specific heat. The dominant SRO peaks now appear at the points corresponding to the √ 3 × √ 3 order, i.e., q = (±4π/3, 0) and (±2π/3, ±2π/ √ 3), moving from the ones corresponding to the q = 0 state at T = 0 of Fig.3(a). Hence, even in the S = 1/2 system, thermal fluctuations select the √ 3× √ 3 SRO as in the classical case. [38][39][40][41][42][43] In order to get further detailed information about the magnetic SRO, we compute the temperature dependence of the S (q, β) intensity at the two representative q-points corresponding to the q = 0 and the √ 3 × √ 3 orders for the sizes N = 18, 27 and 36, each multiple of nine, and the result is shown in Fig.4. On decreasing T , finite-size effects get enhanced indicating the development of the magnetic SRO. At higher T , the √ 3 × √ 3 intensity exceeds the q = 0 one irrespective of the size N. The behavior at lower T , however, turns out to differ significantly between N = 18, 27 and N = 36. For N = 36, on decreasing T , the q = 0 intensity grows while the √ 3 × √ 3 intensity is suppressed, the former exceeding the latter at around T ≃ 0.01 close to the third-peak temperature of the specific heat as mentioned above. Thus, a finitetemperature crossover phenomenon between the two distinct types of magnetic SRO is likely to occur close to the thirdpeak temperature, T ∼ 0.01.
Of course, whether this finite-temperature crossover phenomenon survives or not in the thermodynamic limit is a nontrivial question. Yet, if one notices that the stabilization of the √ 3 × √ 3 SRO at higher T is strongly supported both by our data of Fig.4 and by the result of the semi-classical calculations, while the stabilization of the q = 0 SRO at lower T is supported both by our data for the N = 36 cluster and by the DMRG data of much larger sizes of N ≤ 108, 2, 44, 45 its existence is quite plausible.
The form of S (q, β) around the crossover temperature T cross ∼ 0.01 where the two types of SRO compete would be of special interest. In Fig.3(c), we show S (q, β) at T = 0.01. As can be seen from the figure, the intensity here forms an almost flat ring-like ridge along the BZ boundary. It resembles the intensity of the "ring liquid" proposed in the frustrated honeycomb AF in Ref. 46 and the frustrated square ferromagnet in Ref.. 47 We finally discuss possible implications of our present results to experiments. Recent low-energy inelastic neutron scattering measurements on S = 1/2 kagome AF herbertsmithite revealed the broad spots corresponding to the q = 0 SRO at a low T of T = 2K corresponding to T ∼ 0.01. 48,49 Although the authors of Refs. 48,49 invoked that the impurity effects coming from the Cu 2+ impurities on adjacent triangular (Zn) interlayers as an origin of these spots, our present results indicate that the emergence of the q = 0 SRO alone is understandable even without invoking such impurity effects. Of course, the impurity (or the quenched randomness) effect originated from the adjacent triangular layer could be important in understanding the observed spin-liquid-like behavior as recently emphasized in Refs.. 50 Furthermore, the Dzyaloshinskii-Moriya interaction might be playing a role in stabilizing the q = 0 SRO in real materials. 51 If the possible effects of the randomness and the DM interaction would be negligible in herbertsmithite (might not be the case !), the SRO peaks might move to the distinct q-points corresponding to the √ 3 × √ 3 order as the temperature is further raised. In summary, we studied the effects of thermal fluctuations on the ordering of the S =1/2 kagome Heisenberg antiferromagnet by means of the Hams-de Raedt method up to the N = 36 cluster retaining a full symmetry of the lattice, to find that the second peak of the specific heat persists, while the third and the fourth peaks appear at lower T . In particular, we observed a finite-temperature crossover phenomenon occurring at T ∼ 0.01 close to the third-peak temperature, which was associated with the changeover of the type of the magnetic SRO between the q = 0 (lower-T ) and the √ 3 × √ 3 (higher-T ) states.