Thermodynamics of the spin-half square-kagome lattice antiferromagnet

Over the last decade, the interest in the spin-$1/2$ Heisenberg antiferromagnet (HAF) on the square-kagome (also called shuriken) lattice has been growing as a model system of quantum magnetism with a quantum paramagnetic ground state, flat-band physics near the saturation field, and quantum scars. Here, we present large-scale numerical investigations of the specific heat $C(T)$, the entropy $S(T)$ as well as the susceptibility $\chi(T)$ by means of the finite-temperature Lanczos method for system sizes of $N=18,24,30,36,42,48$, and $N=54$. We find that the specific heat exhibits a low-temperature shoulder below the major maximum which can be attributed to low-lying singlet excitations filling the singlet-triplet gap, which is significantly larger than the singlet-singlet gap. This observation is further supported by the behavior of the entropy $S(T)$, where a change in the curvature is present just at about $T/J=0.2$, the same temperature where the shoulder in $C$ sets in. For the susceptibility the low-lying singlet excitations are irrelevant, and the singlet-triplet gap leads to an exponentially activated low-temperature behavior. The maximum in $\chi(T)$ is found at a pretty low temperature $T_{\rm max}/J=0.146$ (for $N=42$) compared to $T_{\rm max}/J=0.935$ for the unfrustrated square-lattice HAF signaling the crucial role of frustration also for the susceptibility. We find a striking similarity of our square-kagome data with the corresponding ones for the kagome HAF down to very low $T$. The magnetization process featuring plateaus and jumps and the field dependence of the specific heat that exhibits characteristic peculiarities attributed to the existence of a flat one-magnon band are as well discussed.


I. INTRODUCTION
The spin-1/2 Heisenberg antiferromagnet (HAF) exhibits Néel semi-classical magnetic long-range order on most of the two-dimensional (2D) lattices. Among the eleven 2D Archimedean tilings [1,2], only the celebrated kagome lattice and the so-called star lattice have nonmagnetic spin-liquid ground states [3][4][5][6][7][8][9]. Both lattices have low coordination number z, and they are highly frustrated. A non-Archimedean lattice without magnetic long-range order is the square-kagome (sometimes also called shuriken) lattice [10][11][12][13]. Similar to the kagome lattice it is a 2D tiling of corner-sharing triangles, i.e. the classical ground state of the square-kagome Heisenberg antiferromagnet (SKHAF) is highly degenerated. There are two non-equivalent sites A and B and nearestneighbor bonds J and J , see the inset in Fig. 1. In our paper we will focus on the balanced situation J = J = 1 which is most similar to the kagome HAF. The corresponding Hamiltonian augmented with a Zeeman term is given by where quantum mechanical operators are marked by a tilde. Over the last decade, the interest in the spin-1/2 SKHAF has been growing as a model system of quantum magnetism with a quantum paramagnetic ground state as well as flat-band physics near the saturation field and quantum scars [13][14][15][16][17][18][19][20][21][22][23][24][25]. It is expected that the low-energy physics of the SKHAF is just as fascinating as for the related kagome HAF. While it is not questioned that the strong frustration prevents groundstate magnetic ordering, the specific nature of the quantum ground state is under debate. Thus, a recent Schwinger-boson approach leads to a topological spinliquid ground state with weak nematicity and a small gap [22], whereas a variational Monte-Carlo study is in favor of a gapped pinwheel valence-bond crystal state [25]. Previous exact-diagonalization studies seem to also suggest a finite singlet-triplet gap (spin gap) [12,13], though this issue seems as yet unresolved. However and again similar to the kagome HAF [26][27][28], there are many low-lying singlet excitations filling the singlet-triplet spin gap [12]. In addition to the zero-field properties, the magnetization process is of particular interest. Promi-nent features of the magnetization curve M (B) are a wide plateau at 1/3 of the saturation magnetization M sat and a jump to saturation with a preceding plateau at (2/3)M sat [12,13]. A further motivation to study this model comes from recent experimental investigations. A gapless spin liquid was dicovered in the square-kagome magnet KCu 6 AlBiO 4 (SO 4 ) 5 Cl [29,30] which is, however, not a well-balanced (i.e. J = J ) SKHAF; rather in this compound three different exchange couplings are relevant. For Na 6 Cu 7 BiO 4 (PO 4 ) 4 [Cl,(OH)] 3 measurements of the magnetization as well as of the heat capacity indicate the absence of long-range order which might be a signal of spin liquid behavior [30,31].
While almost all previous studies of the SKHAF are focused on ground-state properties, the thermodynamics of the model is much less investigated. Only, in the early paper Ref. [32] the specific heat was calculated by a simple renormalization group approach. In our paper we want to fill this gap of missing finite-temperature studies of the spin-1/2 SKHAF by large-scale numerical simulations of finite lattices of up to N = 54 sites (see Fig. 11 in the Appendix) by means of full exact diagonalization (ED) [33] and of the finite-temperature Lanczos method (FTLM) [34][35][36][37][38][39][40][41][42]. We will use the well investigated kagome HAF as reference system when discussing the data of the SKHAF.
The very existence of an excitation gap is crucial for low-temperature thermodynamics at low temperatures T . Thus, a singlet-triplet gap leads to an exponentially activated low-temperature behavior of the susceptibility, whereas low-lying non-magnetic singlet excitations are relevant for the specific heat. As indicated by previous Lanczos data for finite lattices up to N = 36 sites [12] similarly to the kagome HAF we may expect a small singlet-triplet gap filled with a noticeable number of singlet states.
It is appropriate to mention here, that for smaller cluster sizes some related data have already been reported in previous works. That concerns (i) the zerotemperature magnetzation curves shown in Sec. III B, cf. Refs. [12,13]. However, here we add some new data for larger lattice sizes of N = 48 and N = 54 sites and we present the finite-size dependence of the widths of the magnetization plateaus. Furthermore, (ii) it concerns the specific heat in strong magnetic fields near saturation. In Ref. [12] one can find related specific-heat data for much smaller sizes of N = 18 and N = 24. In our paper we show specific-heat data up to N = 54.
The paper is organized as follow. In Section II we introduce our numerical scheme, in Section III we present our results for the SKHAF and compare them with corresponding data for kagome HAF. In the last Section IV we discuss and summarize our findings. For convenience we show the finite lattices considered here in an appendix.
II. CALCULATIONAL SCHEME The investigated spin system is modeled by the spin-1/2 Heisenberg Hamiltonian given in Eq. (1). The complete spectrum for the spin-half system can be calculated only for the smallest finite lattice of N = 18 spins. For larger systems we perform full diagonalization in high sectors of magnetization M (e.g., for N = 42 in subspaces with M = 16, . . . , 21). The obtained exact energy spectrum yields the contribution of these M sectors to the partition function Z(T, B).
For subspaces that are not accessible by full exact diagonalization we use the finite-temperature Lanczos method (FTLM) [34,35,43,44] which provides approximations of thermodynamic quantities with remarkable accuracy [36,37,42]. Within the FTLM scheme the sum over an orthonormal basis in the partition function is replaced by a much smaller sum over R random vectors: Here the | ν label random vectors for each symmetryrelated orthogonal subspace H(γ) of the Hilbert space, where γ denotes the respective symmetry. In Eq. (2) the exponential of the Hamiltonian is approximated by its spectral representation in a Krylov space spanned by the N L Lanczos vectors starting from the respective random vector | ν , where | n(ν) is the n-th eigenvector of H ∼ in this Krylov space. This method is known to provide accurate data for typical observables such as magnetization, uniform magnetic susceptibility and specific heat, see, e.g. [42].
Naturally, we take into account the commutation of the Hamiltonian H ∼ with the z-component of the total spin S ∼ z = i s ∼ z i to decompose the full Hilbert space into much smaller orthogonal subspaces which can be labeled by the magnetization M = S ∼ z . For a further decomposition of the Hilbert space by employing lattice symmetries we use Jörg Schulenburg's publicly available spinpack package [45,46].
For the number of random vectors R used to approximate the partition function we chose at least R = 20 which is sufficiently large to ensure very accurate FTLM data, cf. Refs. [39,40]. However, for N = 42, where the largest Hilbert-space dimension is N H = 7.34 × 10 10 , we used R = 5 in the subspaces of M = 0 (that contains the ground state and the lowest energy levels), R = 2 for M = 1 and M = 2, and then R = 10 for 2 < M < 16. This strategy was used already in Ref. [39], where also the accuracy of this approach was evaluated.

III. THE SQUARE-KAGOME LATTICE ANTIFERROMAGNET
In what follows we discuss the Wilson ratio, the density of states, the specific heat, the uniform susceptibility, the entropy and the magnetization process.

A. Zero-field properties
First we study the modified Wilson ratio [47,48] where χ is the uniform magnetic susceptibility and S is the entropy. It measures the ratio of the density of magnetic excitations with M > 0 and the density of all excitations including singlet excitations with M = 0. As demonstrated for the kagome HAF [47,48], a vanishing P as temperature T → 0 is a hallmark of quantum spinliquid ground state with dominating singlet excitations at low T . On the other hand, for quantum spin models with semi-classical magnetic ground-state order, such as the square-lattice HAF the Wilson ratio diverges as P (T → 0) ∝ T η , η ≥ 1. . Note that the data for N = 36 (not shown) practically coincide with the data for N = 42. Note further that corresponding data for N = 36 for the kagome and the square-lattice HAF where previously presented in Ref. [47]. Inset: Sketch of the square-kagome lattice. Here A and B label the two non-equivalent sites and J and J label the two non-equivalent nearest-neighbor bonds. All calculations are preformed for the balanced case J = J.
In Fig. 1 we show the Wilson ratio for the N = 42 SKHAF in comparison with the N = 42 kagome HAF (both with non-magnetic spin liquid ground state) as well as for the N = 40 square-lattice HAF (with a magnetically ordered ground state). The striking accordance of the kagome and square-kagome data is obvious signaling the dominance of singlets as T → 0. This behavior is in agreement with the findings reported in Ref. [12] where a noticeable number of singlets was found below the first triplet excitation. Next we discuss the specific heat C(T ), the entropy S(T ) and the uniform susceptibility χ(T ). We use a logarithmic scale for T that makes the low-temperature features transparent, see Figs. 2(a), 4(a), and 5(a). In corresponding panels (b) we compare kagome and squarekagome results using a linear temperature scale. The position T max = 0.67 and height C max = 0.189N of the main maximum of C(T ), set by the exchange coupling J = 1, coincide for both models [49].
Below T ∼ 0.2 the curvature of C(T ) changes and a shoulder-like profile or an additional low-temperature maximum appears for 0.05 T 0.2. The very existence of such an unconventional low-temperature feature below the main maximum seems to be size-independent, where the shrinking height of the low-T maximum with as a function of the respective excitation energy E * : total density of states -black, total density of states for |M | = 1red, for |M | = 2 -green, and for |M | = 3 -blue (curves from left to right). For technical details, see [50].
growing N [marked by the red arrow in Fig. 2(a)] indicates that there will survive rather a shoulder than an additional maximum in the thermodynamic limit. At very low temperatures the singlet excitations dominate the temperature dependence of C. For N = 36 and N = 42 the singlet-singlet gap is the smallest and therefore C remains non-zero even below T ∼ 0.01.
In Fig. 2(b) we compare C(T ) of the SKHAF and of the kagome HAF for N = 42. There is an almost perfect coincidence down to T = 0.3. Below this temperature deviations between the curves are obvious, which can be attributed to subtle details in the low-energy singlet excitation spectrum (see also the discussion of the susceptibility given below, where singlet excitations do not play a role). However, the shoulder-like profile is present in both systems, whereas the sharp low-T peak in C(T ) below the shoulder observed for the kagome HAF is absent for the SKHAF. Note that very recently such a shoulder of the low-temperature specific heat has been found in the kagome quantum antiferro- To shed light on the relevance of low-lying excitations of different sectors of M for the low-temperature behavior of C(T ) we show the contributions of the different sectors of total magnetization M to the density of states n(E * ) as a function of the respective excitation energy E * in Fig. 3, where n(E * ) is calculated by histograming the Krylov space energy eigenvalues in combination with their respective weights. The dominance of the singlet excitations below E * ∼ 0.1 is obvious. This observation is further supported by the behavior of the entropy S(T ) as shown in Fig. 4(a), where a change in the curvature is present just at about T = 0.2. As for the specific heat finite-size effects become visible below about T = 0.2.  The comparison of the S(T ) profiles of the kagome HAF and SKHAF, see Fig. 4(b), confirms again the striking similarity of both models up to pretty low T .
Next we turn to the zero-field susceptibility χ displayed in Fig. 5. As shown in Fig. 3 the relevant singlettriplet gap (spin gap) is significantly larger than the singlet-singlet gap leading to an exponentially activated low-temperature behavior. Since the singlet-triplet gap shrinks with growing N , cf. Ref. [13], this feature sets in at lower T as increasing N . However, the question of a finite spin gap as N → ∞ seems to be not clarified so far. Since the non-magnetic singlet excitations are irrelevant for the susceptibility, the temperature profiles of χ of the SKHAF and the kagome HAF show an excellent agreement also below T = 0.3, where the specific heat starts to deviate for both models. As already discussed in Refs. [39,52] for the kagome HAF, the maximum in χ(T ) is at a pretty low temperature T max = 0.146 (for the N = 42 SKHAF) compared to T max = 0.935 for the square-lattice HAF [53,54], thus signaling the crucial role of frustration also for the susceptibility.  The magnetization process of strongly frustrated quantum magnets exhibits a number of interesting features, such as plateaus and jumps [55]. We start with a brief discussion of the ground-state magnetization curve M (B) of the SKHAF, see Fig. 6. In previous studies [12,13]  . The saturation field gµ B B sat = 3 is the same as for the kagome HAF. Magnetization plateaus exist at 1/3 and 2/3 of the saturation magnetization M sat for the infinite system at T = 0. The widths W of the plateaus are pretty large. From our data we estimate W 1/3 ≈ 0.368B sat and W 2/3 ≈ 0.123B sat (see the inset in Fig. 6), which is larger than corresponding plateau widths of the kagome HAF [56,57]. Moreover, there is the typical macroscopic jump to saturation due to the presence of independent localized multi-magnon eigenstates stemming from a flat onemagnon band [58][59][60][61][62]. Its existence is analytically proven and it does not exhibit finite-size effects [58,59]. The very existence of a flat one-magnon band is connected with destructive quantum interference which is related to the geometric structure of corner-sharing triangles. To get an imperession of the band-structure of the one-magnon excitations above saturation we refer the reader to Ref. [63], where the band structure is shown for a fermionic model on the square-kagome lattice which is closely related in terms of the one-particle spectrum.
Similar to the kagome HAF, see [56,57,64], the plateau states of the spin-half SKHAF are non-classical valence bond states. In the upper plateau it is the exactly known magnon-crystal state built of one-magnon states on the squares and z-aligned spins on intermediate A sites on the connecting triangles [12,17,58]. The lower 1/3 plateau state is not exactly known. The inspection of the expectation values s ∼ z i and s ∼ i · s ∼ j provides strong evidence that the squares carry an almost perfect We present the influence of the temperature on the magnetization curve in Fig. 7. It is obvious that for elevated temperatures the experimental detection of plateaus becomes difficult, particularly, if the plateau width is only of moderate size, see the upper part of the M (B, T ) curve. Moreover, we notice that the jump of the magnetization to saturation at T = 0 is washed out already for small temperatures.
To detect plateaus and jumps in experiments the first derivative dM/dB as a function of T is more suitable, cf., e.g., Ref. [65]. We show dM/dB for N = 42 in Fig. 8. We notice first that the oscillations present for the lowest temperature T = 0.03 (red curve) are finite-size effects. The plateaus at 1/3 and 2/3 show up as pronounced minima in dM/dB, however, the detection of such minima requires sufficiently low temperatures. The jump of the magnetization to saturation at T = 0 (washed out in the M (B) curve at T > 0) leads to a high peak in dM/dB at the saturation field. The smaller peak at the upper end of the 1/3 plateau is related to the enlarged step (twice as high as the normal finite-size step). Contrary, to the jump to saturation we may expect that it will disappear as N → ∞. The melting of the pronounced plateaus at 1/3 and 2/3 with growing temperature happens more rapidly for the upper 2/3 plateau. Furthermore, the 2/3 plateau melts asymmetricly, i.e. the minimum in dM/dB is below the midpoint of the plateau.
The influence of the magnetic field B on the specific heat C for N = 42 is depicted in Fig. 9(a) for selected B values as shown in the inset of Fig. 9(a). Below the 1/3 plateau and at very low temperatures the influence of B is determined by the shift of the low-lying magnetic excitations with M = 1, 2 and 3 towards and even beyond the zero-field singlet ground state. As a result, the low-T shoulder [present for B = 0, line 1 in Fig. 9(a)] may become a low-temperature peak in C(T ) (line 2), where its position and height depends on B. A similar low-T scenario is observed for B between the 1/3 and the 2/3 plateaus (line 4) and for B above the 2/3 plateau (line 6). A striking low-T feature for B-values inside the plateaus (lines 3 and 5), where the ground state is a valence-bond state (see above), is the well-pronounced pretty high extra peak. The extra peak for B within the 2/3 plateau is well-understood. It is a flat-band effect and is related to the huge manifold of low-lying localized multi-magnon states [19,39,40,61,62,66]. Apparently, it is size-independent (see Fig. 9(b)), i.e., it persists for N → ∞. It is worth mentioning the contrast to the kagome HAF which is related to the different structure of the respective localized-magnon states. While these states for the SKHAF are located on isolated trapping cells (squares), the trapping cells (hexagons) of the kagome HAF are connect leading to a repulsive interaction of the traps. As a result, this extra maximum of the specific heat of the kagome HAF increases with N and in the thermodynamic limit it becomes a true singularity indicating a low-temperature order-disorder transition into a magnon-crystal phase [40,61]. For the extra peak when B is inside the 1/3 plateau there is no straightforward explanation. Our data for N = 30, 36, 42 shown in Fig. 9(b) indicate that there is only a small finite-size effect and the height of the maximum is slightly growing with N . We may conclude that most likely it also persists for N → ∞.
Let us turn to the main maximum around T ∼ 1, see Fig. 9(a). At a first glance, the variation of this maximum around T ∼ J seems to be not very systematic. However, as discussed for the kagome HAF [39] the effect of B on the main maximum of C is influenced by the huge manifold of localized-magnon (flat-band) states. We compare the height C max and the position of the main maximum T max for the N = 42 SKHAF, the N = 42 kagome HAF and the N = 40 square-lattice HAF in Fig. 10. At low magnetic fields the value of T max is determined by J and the coordination number z and therefore the behavior of T max is very similar for all models. While C max remains almost constant until B ∼ 0.8B sat and then it increases smoothly for B > B sat the variation of T max as a function of B is more pronounced. It exhibits two maxima at B = 0 and B ≈ 1.1B sat and two minima at B ≈ 0.5B sat and B ≈ 1.4B sat as well as two regions with an approximately linear increase of T max . Again, the remarkable agreement of the SKHAF and the kagome HAF is obvious. On the other hand, the difference to the square-lattice HAF in C max and more pronounced in T max is striking. While the difference in C max at low T is related to the different low-energy physics which affects C at higher T according to the sum rule ∞ 0 (2), the different behavior of T max (B) beyond B ∼ 0.5B sat signals flat-band effects, i.e. the exponentially growing number of localized multi-magnon states [66] yields an increase of T max up to a noticeable maximum around B = B sat . It follows a linear increase of T max above B sat present in all models which is linked to the paramagnetic phase.

IV. CONCLUSIONS
In our study we performed large-scale calculations of thermodynamic quantities such as the magnetization M (T ), the specific heat C(T ), the entropy S(T ) and the susceptibiliy χ(T ) for the highly frustrated spin-half square-kagome Heisenberg antiferromagnet (SKHAF). For that we used the finite-temperature Lanczos method (FTLM) applied to seven different finite lattices including the 'large' lattices of N = 42, 48 and 54 sites, where for N = 48 and 54 only the thermodynamics near the saturation field was considered. At zero magnetic field, we find a remarkable accordance of the thermodynamic properties of the SKHAF with those of the paradigmatic kagome HAF. Our results for N = 42 and smaller sizes indicate that the specific heat very likely has got a low-temperature shoulder instead of an additional lowtemperature maximum. There is also a considerable influence of frustration on the susceptibility and on the entropy. Thus we find a noticeable shift of the maximum of χ(T ) to low T compared to two-dimensional unfrustrated HAFs, and, the entropy per site S(T )/N acquires about 40% of its maximum value ln 2 already at T /J ∼ 0.1.
Other interesting properties of the SKHAF are related to the magnetization process in an applied magnetic field B. There are two well pronounced plateaus at 1/3 and 2/3 of the saturation magnetization and a jump from the 2/3 plateau directly to saturation caused by the flat onemagnon band. The melting of the plateaus with growing temperature is faster for the upper 2/3 plateau and it melts asymmetricly. For magnetic fields inside the plateaus the specific heat shows a well-pronounced extra peak at low temperatures, which exhibits only small finite-size effects. Furthermore, we find that the influence of strong frustration is not only visible at low temperatures T J, it is also noticeable at moderate (and high) temperatures T ∼ J. Especially, the presence of a flat one-magnon band leading to a huge manifold of low-lying flat-band states yields pronounced effects in the magnetization process and the temperature dependence of the specific heat at magnetic fields above B ∼ 0.5B sat .
Though our investigation of the spin-1/2 SKHAF as a highly frustrated quantum spin system is of interest in its own right, it is also motivated by the recent discovery of a spin liquid in the square-kagome magnet KCu 6 AlBiO 4 (SO 4 ) 5 Cl [29], which exhibits, however, three different exchange couplings. Moreover, the large variety of magnetic insulators [67,68] as well as the progress in synthesizing new magnetic molecules and compounds with predefined spin lattices may open the window to get access to the observation of the discussed phenomena.
Bearing in mind the numerous studies of the lowenergy physics of the related kagome HAF we argue that our work may also stimulate other studies using alternative techniques, such as tensor network methods, DMRG, numerical linked cluster expansion or Green's function techniques [69][70][71][72][73].
ACKNOWLEDGMENT This work was supported by the Deutsche Forschungsgemeinschaft (DFG RI 615/25-1 and SCHN 615/28-1). Computing time at the Leibniz Center in Garching is gratefully acknowledged.
Appendix A: Finite square-kagome lattices used for the exact diagonalization and the finite-temperature Lanczos method Here, we provide the employed lattice structures in Fig. 11. [1] An Archimedean lattice is built by uniform polygons. Importantly, it has only one set of topologically equivalent sites.