Thermodynamic properties of the Shastry-Sutherland model throughout the dimer-product phase

The thermodynamic properties of the Shastry-Sutherland model have posed one of the longest-lasting conundrums in frustrated quantum magnetism. Over a wide range on both sides of the quantum phase transition (QPT) from the dimer-product to the plaquette-based ground state, neither analytical nor any available numerical methods have come close to reproducing the physics of the excited states and thermal response. We solve this problem in the dimer-product phase by introducing two qualitative advances in computational physics. One is the use of thermal pure quantum (TPQ) states to augment dramatically the size of clusters amenable to exact diagonalization. The second is the use of tensor-network methods, in the form of infinite projected entangled pair states (iPEPS), for the calculation of finite-temperature quantities. We demonstrate convergence as a function of system size in TPQ calculations and of bond dimension in our iPEPS results, with complete mutual agreement even extremely close to the QPT. Our methods reveal a remarkably sharp and low-lying feature in the magnetic specific heat around the QPT, whose origin appears to lie in a proliferation of excitations composed of two-triplon bound states. The surprisingly low energy scale and apparently extended spatial nature of these states explain the failure of less refined numerical approaches to capture their physics. Both of our methods will have broad and immediate application in addressing the thermodynamic response of a wide range of highly frustrated magnetic models and materials.


I. INTRODUCTION
Frustrated quantum spin systems, particularly those forbidding magnetically ordered ground states, provide one of the most important avenues in condensed matter physics for the realization and investigation of phenomena ranging from fractionalization to many-particle bound states, from massive degeneracy to topological order, and from quantum entanglement to quantum criticality [1].Significant progress has been made in describing the ground states of many such quantum magnets in two and three dimensions (2D and 3D), including different types of valence-bond crystal and both gapped and gapless quantum spin liquids [2].However, a full understanding of the excitation spectrum and thermodynamic response of frustrated quantum spin models, which is the key to experimental interpretation, has remained elusive in all but a small number of special cases.This is because no unbiased analytical or numerical methods exist by which to study these properties for general systems at low but finite temperatures and on lattices whose sizes approximate the thermodynamic limit.* awietek@flatironinstitute.org A case in point is the Shastry-Sutherland model [3], which was constructed explicitly to have as its ground state a product state of dimer singlets.Also referred to as the orthogonal dimer model [4], it is defined by the Hamiltonian which is represented schematically in Fig. 1.J D is an intra-dimer coupling (denoted by ij ) and J a mutually frustrating inter-dimer coupling (denoted ij ).The ground state is clearly a dimer-product phase at small J/J D and a square-lattice antiferromagnet (AF) at large coupling ratios.As intimated above, the groundstate phase diagram is in fact quite well known, with the most quantitatively accurate results provided by the ansatz of infinite projected entangled-pair states (iPEPS, which we will extend here to finite temperatures) [5].
As shown in Fig. 1, the singlet dimer-product state is found up to J/J D = 0.675 (2) and the AF state beyond J/J D = 0.765 (15), with a gapped, plaquette-based state occurring in the intermediate regime.The lower quantum phase transition (QPT) is of first order; the upper one has also been found numerically to be first-order, but only weakly so [5], and this has also been questioned [6].
It is the physics in the vicinity of the dimer-to-plaquette dimer plaquette Néel FIG. 1. Geometry and coupling constants, JD and J (center panel), of the Shastry-Sutherland model [Eq.(1)].The three panels show schematic representations of the singlet dimerproduct state, which is exact, the plaquette state, and the Néel phase.Shaded ellipses and squares denote respectively dimer-and plaquette-singlet states.Phase boundaries are taken from Ref. [5].
QPT that has posed a long-standing challenge.Somewhat remarkably, the Shastry-Sutherland model has a faithful realization in the material SrCu 2 (BO 3 ) 2 , which has exactly the right geometry and only weak non-Heisenberg (primarily Dzyaloshinskii-Moriya) interactions.Discovered in 1991 [7] and first studied for its magnetic properties 20 years ago this year [8], SrCu 2 (BO 3 ) 2 lies in the dimer-product phase [9] and shows a very flat band of elementary triplet ("triplon") excitations [10,11].Excellent magnetic susceptibility [8,12] and specific-heat [13] measurements made at this time were used to estimate the coupling ratio as J/J D 0.635, thereby placing SrCu 2 (BO 3 ) 2 rather close to the dimer-plaquette QPT.Detailed experiments performed in the intervening two decades have revealed a spectacular series of magnetization plateaus [8,[14][15][16][17][18][19][20][21], as well as a curious redistribution of spectral weight at temperatures very low on the scale of the triplon gap [11,22].Of most interest to our current study is the result [23] that an applied pressure makes it possible to increase the ratio J/J D to the extent that, at approximately 1.9 GPa, the material is pushed through the QPT into a plaquette phase.Data for the specific heat under pressure have only appeared [24] concurrently with the present study, and indicate that the low-temperature peak moves to a lower temperature at 1.1 GPa before evidence of an ordered phase appears at 1.8 GPa.
In contrast to the ground state, the excited states and thermodynamics of the Shastry-Sutherland model are not at all well known around the QPT, despite the attention focused on this regime due to SrCu 2 (BO 3 ) 2 .Full exact diagonalization (ED) has been limited to systems of up to 20 sites and, as we will demonstrate in detail here, suffer from significant finite-size limitations at J/J D > 0.6.Quantum Monte Carlo (QMC) methods suffer from extreme sign problems, and despite a recent study by some of us [25] that expands the boundaries of what is possible by QMC, the low-temperature behavior in the regime J/J D > 0.6 remains entirely inaccessible.iPEPS methods work in the thermodynamic limit and are immune to frustration problems, but to date have been limited to ground-state properties.As a consequence, two decades after their measurement, the best interpretation of the thermodynamic properties of SrCu 2 (BO 3 ) 2 [12,13], and the origin of the estimate J/J D = 0.635, remain based on ED of 16-and 20-site clusters [4,26].
Here we introduce two qualitative technical advances in numerical methods for computing thermodynamic properties and benchmark them by application to the Shastry-Sutherland model in the dimer-product phase up to the QPT.For a quantum many-body system, a thermal equilibrium state can be represented by a typical pure state, which is known as a thermal pure quantum (TPQ) state [27,28].This result can be understood as a consequence of a more general phenomenon, known as quantum typicality [29,30].All statistical-mechanical quantities can be obtained from averaging over a few TPQ states.For practical purposes, obtaining TPQ states by ED allows access to thermodynamic quantities without having to perform a full diagonalization, and hence admits a qualitative increase in accessible system sizes [31,32].We will exploit TPQ methods to reach system sizes of N = 40 in the Shastry-Sutherland problem.
iPEPS are a tensor-network ansatz to represent a quantum state on an infinite lattice [33][34][35], here in 2D, whose accuracy is controlled by the bond dimension, D, of the tensor.In fact this ansatz provides not only a compact representation for the ground states of gapped local Hamiltonians, but also for representing a purification of the thermal density operator, and hence the thermal states of local Hamiltonians.Here we exploit newly developed algorithms [36] allowing the efficient calculation of imaginary-time evolution processes, which in general require additional truncation to retain the value of D at each time step, to generate thermal states of the Shastry-Sutherland Hamiltonian and hence to compute its thermodynamic properties.
Our TPQ and finite-T iPEPS methods enable us to capture the physics of the Shastry-Sutherland model, and thus of SrCu 2 (BO 3 ) 2 , in a way not hitherto possible near the QPT.To illustrate this clearly, in Fig. 2 we show sample results from both methods for the magnetic specific heats and susceptibilities of Shastry-Sutherland models with coupling ratios J/J D = 0.60, 0.63, 0.65, and 0.66.Over this narrow coupling regime near the QPT, the specific heat develops a pronounced two-peak form whose lower peak becomes extremely sharp and moves to a remarkably low energy as the QPT is approached.The corresponding susceptibilities rise increasingly abruptly before an increasingly "square" turnover at their maximum.We draw attention to the fact that the results from both our methods are in excellent quantitative agreement, apart from the specific-heat peak heights extremely close to the QPT.These characteristic shapes, and their systematic evolution, have not been obtained before and shed new light both on the thermodynamic quantities themselves and on the underlying excitation spectrum responsible for their form.
The manuscript is organized as follows.In Sec.II we introduce TPQ states and their application to the calculation of thermodynamic quantities.In Sec.III we provide a brief introduction to iPEPS and discuss the generation of thermal states by their imaginary-time evolution.Section IV presents our complete results from both methods for the magnetic specific heat and susceptibility of the Shastry-Sutherland model over the coupling range 0.60 ≤ J/J D ≤ 0.66, which are summarized in Fig. 2 and which it has not been possible to obtain with any reliability in previous studies.In Sec.V we discuss the physics of the sharp features we observe in the specific heat near the QPT, the limits of our numerical methods in the context of the Shastry-Sutherland model, and the application of our results to the material SrCu 2 (BO 3 ) 2 .Section VI contains a brief summary and perspective concerning the impact of our results in frustrated quantum magnetism.

II. THERMODYNAMICS FROM THERMAL PURE QUANTUM STATES
To evaluate the thermal average of a quantum mechanical observable, A, in the canonical ensemble, one computes where H denotes the Hamiltonian, β = 1/(k B T ) the inverse temperature, k B the Boltzmann constant, and Z = Tr(e −βH ) the canonical partition function.Thermal pure quantum (TPQ) states provide an alternative approach to the evaluation of thermodynamic properties [27,28].The trace of any operator, A, can be evaluated by taking the average of its expectation values over a set of random vectors, where |r denotes a normalized vector with random normal-distributed components and d denotes the dimension of the Hilbert space.Hence, any thermal average can be estimated by evaluating where denotes the "canonical TPQ" state [28] and . . .denotes the average over the random vectors |r .Importantly, the TPQ state, |β , is still random.The variance of the estimate in Eq. ( 4) decreases as 1/ √ d, where d = 2 N for a S = 1/2 lattice model on N sites, and hence becomes exponentially small as a function of the system size under only mild assumptions about the form of the operator H; a precise statement may be found in Ref. [28].For larger system sizes, it is therefore necessary to generate only relatively few random TPQ states, |β , in order to achieve small statistical errors.Because the samples in the denominator and the numerator of Eq. ( 4) are correlated when evaluated from the same set of random vectors, |r , we perform a jackknife analysis [37] to estimate the error bars.
The key advantage of using TPQ states, when compared with evaluating statistical averages using the standard formalism of Eq. ( 2), is that the state |β can be approximated numerically to arbitrary precision without resort to a full diagonalization of the Hamiltonian, H. Variants of the TPQ method have been employed in several numerical studies of static and dynamical observables [31,32,38,39].Whereas Ref. [28] proposed a Taylor-expansion technique, here we improve the speed of convergence by using the Lanczos basis [40] of the Krylov space of H to evaluate the states |β .
The nth orthonormal basis vector in the Lanczos basis, |v n , is given recursively [41] by where |v 1 = |r , b 1 = 0, and By defining the matrix of Lanczos vectors the nth Lanczos approximation of H can be written as where T n is a tridiagonal matrix given by We use this approximation to evaluate where e 1 = (1, 0, . . ., 0) T .In the case that A is a power of the Hamiltonian, A = H k , we have We stress that the only quantities to be computed are the sequences a n and b n .Further, results at all temperatures can be obtained from a single evaluation of the Lanczos basis.The exponentiated tridiagonal matrices, e − β 2 Tn , are computed by the eigendecomposition of the matrices T n and exponentiation of the diagonal matrix of eigenvalues.The eigendecomposition of T n is computed using an implicit QR method [42][43][44] by calling the LA-PACK routine dsteqr [45,46].The dimension, n, of the Krylov space is increased until the desired precision is achieved.The computations in this manuscript required a dimension n < 150 to achieve a precision superior to the statistical errors.
In our computations we use a block-diagonal form of the Hamiltonian, where the blocks H ρ correspond to the Hamiltonian in different lattice-symmetry and total-S z sectors.Statistical averages are evaluated in each individual block and then summed, We employ the techniques proposed in Ref. [47] to perform the computations in the individual symmetry sectors.
It has recently become possible to compute Lanczos approximations in particular symmetry sectors for systems of up to N = 50 S = 1/2 sites for the Heisenberg model on the square lattice and N = 48 S = 1/2 on the kagome lattice [47].It is worth noting here that the latter geometry and model are particularly favorable for TPQ [28] because the density of states is high at very low energies [48].However, deriving thermodynamic quantities from the TPQ method requires that these Lanczos approximations be computed multiple times in each symmetry sector for a number, R, of random "seeds", i.e. independent random realizations of the vector |r in Eq. ( 5), and it is R that determines the statistical error.Thus restrictions on CPU time currently limit the practical application of the TPQ approach to an upper limit of 40-42 S = 1/2 spins.More specifically, the calculations in this manuscript were performed with R = 80 seeds for all coupling ratios J/J D on clusters of N = 32 spins and R = 20 for the N = 36 cluster, except at J/J D = 0.63, where R = 40 seeds were evaluated.For the N = 40 cluster we used R = 5.In total, the TPQ calculations in this manuscript required approximately 5.5M CPU hours.
We comment that our implementation of the TPQ method is reminiscent of the finite-temperature Lanczos method [49], which has also recently been pushed to a cluster of N = 42 S = 1/2 spins for the kagome lattice [50].Although the precise relation between these two different methods remains to be understood, currently we believe that the TPQ approach has a more systematic theoretical foundation and as a consequence that the errors within this approach are better controlled.
For system sizes N ≤ 20, full diagonalization of the Hamiltonian yields numerically exact results with no statistical errors.An intermediate method is to compute a large number of low-lying eigenstates, for example by using the variant of the Lanczos method outlined in Ref. [51], and to approximate the low-temperature thermodynamics using these.Despite the significant advantage of also providing results that are entirely free of statistical errors, this method remains restricted by the fact that the required numbers of eigenvalues are large even at the surprisingly low temperatures where the dominant physical processes occur in the present case (Sec.IV).Thus we present benchmarking Lanczos ED results for cluster sizes up to N = 28, although we comment that N = 32 would be accessible were we to invest CPU times comparable to those used in our TPQ calculations.
As an illustration of the manner in which TPQ system sizes allow access to the thermodynamic properties of the Shastry-Sutherland model, in Fig. 3(a) we show the specific heat per dimer, where E ≡ H is the energy, obtained from our TPQ calculations for coupling ratio J/J D = 0.63 with N = 32, 36, and 40. Figure 3(b) shows the corresponding magnetic susceptibility per dimer, where h is a magnetic field applied along the z-axis and m ≡ M , with denotes the total magnetization.In both panels of Fig. 3, we compare these TPQ results with ED calculations for systems of sizes N = 16 and 20, and at low temperatures N = 24 and 28.The details and context of these results will be discussed in full in Secs.IV and V; for the purposes of this introduction, we stress that both C(T ) and χ(T ) have sharp physical features, namely the low-temperature peak in the former and shoulder in the latter, whose shape is manifestly not well reproduced at small N but becomes very much clearer as N is increased.From Fig. 2, this type of access to larger N values becomes increasingly important as the QPT is approached.

III. THERMODYNAMICS FROM IPEPS
An infinite projected entangled-pair state (iPEPS), as a tensor-network ansatz representing a 2D quantum state in the thermodynamic limit [33][34][35], can be seen as a natural generalization of infinite matrix-product states (MPS) to two, or indeed higher, dimensions.While the rank of the tensors used in the representation is determined by the physical problem, the amount of information they contain is determined by their bond dimension, D, which is used to control the accuracy of the ansatz.Although iPEPS were introduced originally for representing the ground states of local Hamiltonians, more recently several iPEPS methods have been developed for the representation of thermal states [36,[52][53][54][55][56][57][58][59][60][61][62].Here we focus on the approaches discussed in Ref. [36], where an iPEPS is used to represent a purification of the 2D thermal density operator.We comment that MPS-based methods have also been applied recently [63][64][65] to compute the thermal response of finite 2D systems in cylindrical geometries; these approaches provide accurate results for cylinders up to a circumference W , which is limited by an exponential increase of the required bond dimension with W .
An iPEPS consists of a unit cell of tensors which is repeated periodically on an infinite lattice.For a ground state (pure state), each tensor has one physical index, p, describing the local Hilbert space of the dimer and four auxiliary indices, r, t, l, and b [Fig.4(a)], each with bond dimension D, which connect the tensor to its four nearest neighbors on a square lattice.For the Shastry-Sutherland model [Fig.1] we require one tensor per dimer, arranged in a unit cell with two tensors, A and B, in a checkerboard pattern, as shown in Fig. 4(b).
We represent a thermal density operator, ρ(β) = e −β Ĥ , at inverse temperature β by its "purification," |Ψ(β) , which is a pure state in an enlarged Hilbert space where each physical site, j (here a dimer with local dimension d = 4), is accompanied by an ancilla site with the same local dimension, described by a local basis |p, a j , where p and a denote respectively the local physical and ancilla basis states.The representation of the thermal state is through its thermal density operator, ρ(β) = Tr a |Ψ(β) Ψ(β)|, where the trace is taken over M 2 m 6 I N w V t 9 e Z 2 0 q x X v u l J t 1 s r 1 W h 5 H A c 7 h A q 7 A g x u o w z 0 0 o A U M E J 7 h F d 6 c R + f F e X c + l q 0 b T j 5 z B n / g f P 4 A w f m M 3 A = = < / l a t e x i t > t < l a t e x i t s h a 1 _ b a s e 6 4 = " R w v g B W 0 9 7 n X t t m U g 7 x z J J Z all ancilla degrees of freedom, a. Physically, the ancilla states act as a perfect heat bath, and taking the trace gives exact thermodynamic averages.In tensor-network notation, each ancilla site is incorporated through the additional index a [blue lines in Fig. 4(c)], leading to the diagrammatic representation of |Ψ(β) shown in Fig. 4(d) and of ρ(β) in Fig. 4(e) [66].
Each application of Ûij increases the bond dimension between the tensors on sites i and j to a larger value, D > D, which for reasons of efficiency must be truncated back to D at every step.As in the imaginary-time evolution of ground states, there exist several approaches for the truncation of a bond index, the two most common being the simple- [68] and full-update methods [34,69].In the simple-update approach, the truncation is performed using a local approximation of the state, which is computationally cheap but does not yield the optimal truncation.In the full-update approach, a bond is truncated in a manner optimal under the method by which the entire state is taken into account; however, this is considerably more computationally expensive and in general it is not possible to reach bond dimensions as large as by the simple-update technique.Below we compare results obtained by these two approaches.Further details concerning the iPEPS ansatz and algorithms for thermal states may be found in Ref. [36].
For contraction of the thermal-state tensor network, which is required to compute all physical expectation values, we use a variant [70] of the corner-transfer-matrix (CTM) renormalization-group method [71,72] in which the accuracy is controlled by the boundary bond dimension, χ.To improve the efficiency of the calculations we exploit the global U(1) symmetry [73,74] of the model of Eq. ( 1).Further details of the iPEPS methods we employ here may be found in Refs.[5,69,75].
Turning to the evaluation of the thermodynamic properties of the Shastry-Sutherland model with iPEPS, we obtain the magnetic specific heat [Eq.( 16)] from the numerical derivative of the energy as a function of the temperature and the magnetic susceptibility [Eq.(17)] from the numerical derivative of the magnetization with respect to a small external magnetic field, h.To illustrate the dependence of our results on the parameters τ , D, and χ intrinsic to the iPEPS procedure, in Figs. 5 and 6 we show example data for the coupling ratio J/J D = 0.63.
Our primary results are obtained using a second-order Trotter-Suzuki decomposition with a time step τ = 0.04.To demonstrate that this value is sufficiently small, meaning that the Trotter discretization error is negligible compared to the finite-D errors (below), in Fig. 5 we show the comparison with results obtained using a time step τ = 0.005.Further, it is easy to show that a relatively small boundary bond dimension, such as χ = 50 for D = 14, is large enough for the accurate computation of expectation values.As both panels of Fig. 5 make clear, results obtained using χ = 80 are almost identical to χ = 50.To compute the magnetic susceptibility we use a field h/J D = 0.01, which is sufficiently small that the discretization error is not significant; a comparison with h/J D = 0.005 is shown in Fig. 5(b).
To investigate the dependence of our results on the finite value of D, in Fig. 6 we present iPEPS data for J/J D = 0.63, obtained by the simple-update method for a range of different bond dimensions.We find in Fig. 6(a) that the position of the specific-heat peak shifts slightly (4%) downwards in temperature with increasing D, while its height is also reduced (2%).However, these changes take place largely from D = 12 to D = 16, whereas changes between the D = 16 and D = 18 data are scarcely discernible on the scale of the figure.In the susceptibility [Fig.6(b)], finite-D effects are clearly very small for T /J D > 0.2; for 0.05 < T /J D < 0.2 a small increase can be observed with increasing D. Finally, to investigate the difference between results by using the simple-and full-update approaches, we return to Fig. 5.At fixed D = 14, this difference is relatively small in the specific heat, amounting to changes of less than 2% in the position and height of the peak [inset, Fig. 5(a)].For the susceptibility, differences between the two approaches are minimal [Fig.5(b)].Because this difference is small in comparison with the finite-D effects shown in Fig. 6, we focus henceforth on the simple-update approach, which allows us to reach bond dimensions as large as D = 20 in the present problem, rather than the full-update approach, where D = 14 represents a practical upper limit.

IV. THERMODYNAMICS OF THE SHASTRY-SUTHERLAND MODEL WITH
0.60 ≤ J/JD ≤ 0.66 We begin the systematic presentation and comparison of our thermodynamic results by making contact with the previous state of the art, as shown in Ref. [25].The analysis of these authors, which included ED, QMC, and high-temperature series expansions (HTSEs), terminated at the coupling ratio J/J D = 0.60, where all three methods were beyond their limits.However, it was shown that QMC remains a highly accurate technique, with only minimal sign problems and thus applicable with large system sizes, for J/J D < 0.526.Thus we first show, in Fig. 7, the magnetic specific heat and susceptibility for J/J D = 0.50 computed by TPQ for system sizes N = 32 and 36 and by iPEPS with D = 16 and 18.By comparison with the large-system QMC benchmark (N = 200 sites), we observe that our iPEPS data, arguably a product of the newest and least well-characterized technique, agree precisely for both bond dimensions shown.
For ED-based methods, it was known already [25] that clusters of size N = 20 give an adequate account of the thermodynamics at this coupling ratio, with the exception of the very peak in C(T ) [Fig.7(a)].Nevertheless, our TPQ results contain two surprises, in the form of the rather large errors and the significant difference between results for the two cluster sizes.Regarding how accurately TPQ reproduces both C(T ) and χ(T ) at temperatures around the C(T ) peak, we note that our N = 32 results do indeed lie close to the QMC benchmark for this system size, with the discrepancy being well within the error tube.Quite generally, the origin of these somewhat large errors lies in the very low density of low-lying states in the spectrum, as pointed out in Ref. [28], and this is certainly the case in the Shastry-Sutherland model when the gap is as large as its value at J/J D = 0.5 (more details may be found in Sec.VA).We will show that our TPQ results become progressively more accurate as the system approaches the QPT, which reduces the gap, until coupling ratios very close to the transition.With regard to the effects of the shape or geometry of the cluster, regrettably the method allows little control beyond the absolute value of N .
In Fig. 8 we move into the previously intractable parameter regime by showing the specific heat and susceptibility at J/J D = 0.60 computed by TPQ for system sizes N = 32 and 36 and by iPEPS with D = 14, 16, and 18. Focusing first on the specific heat [Fig.8(a)], we still observe a single peak at low temperatures (T /J D 0.12), but followed by a very flat plateau extending to T /J D 0.45.It is clear that our TPQ results, augmented by an ED calculation with N = 28 reaching temperatures just above the peak, exhibit finite-size effects in this region.Both the position and the sharpness of the peak show a significant evolution as the system size is augmented, although from N = 32 to 36 the difference is within the error bars of the method.As expected from Sec. III, the three iPEPS bond dimensions show good convergence, in fact to a peak shape exactly between the two TPQ estimates.Outside the peak region, all methods and sizes agree extremely well, and also converge to the values obtained by QMC, a method that returns results of sufficiently low statistical uncertainty at high temperatures (T /J D 0.3) for this coupling ratio [25].
The corresponding susceptibility [Fig.8(b)] is less sensitive to differences in either system size or bond dimension.Indeed, only an ED calculation with N = 20 is not capable of capturing the maximum of χ(T ) for this coupling ratio.However, the two TPQ results do differ concerning the exact location in temperature of the rapid rise in χ(T ), and the fully convergent iPEPS results appear to offer the benchmark.All of our calculations agree closely on the location of the broad maximum.We comment that the characteristic temperature of χ(T ) need not match that of C(T ), given that the former has contributions only from magnetic states but the latter includes the contributions of singlets.However, at this coupling ratio, where the one-triplon excitations of the system lie below all two-triplon bound states [25], the peak in C(T ) and the rise of χ(T ) do indeed coincide.We discuss the excitation spectrum of the model in detail in Sec.VA.
As noted in Sec.I, the coupling ratio J/J D = 0.63 has special significance in the Shastry-Sutherland model because of its proposed connection to the material SrCu 2 (BO 3 ) 2 , and for this reason we have already used it for benchmarking purposes in Secs.II and III.In Fig. 9 we compare our TPQ and iPEPS results for C(T ) and χ(T ).Independent of either method, we observe that the specific-heat peak has become sharper, although not any taller, and that the drop on its high-T side has turned into a true minimum, before C(T ) recovers to a broad maximum around T /J D = 0.45.The position of the sharp peak has moved down to T /J D 0.10, which is a remarkably large shift for such a small change in coupling ratio.Although the corresponding features are far less pronounced in χ(T ), it it clear that onset is considerably steeper than at J/J D = 0.60 and that the temperature scales both for it and for the flat maximum are lower.Focusing on the details of C(T ) [Fig.9(a)], once again our TPQ results for N = 32 and 36, particularly in combination with N = 28 ED, show nontrivial changes around the peak, both in its position and its height.However, the N = 40 data appear to offer a systematic interpolation of both quantities.In our iPEPS results, D = 14 no longer appears representative of the large-D limit at this coupling ratio, but finite-D effects seem to become very small at higher D. Qualitatively, our optimal iPEPS peak appears to be very close to the form at which one might expect the TPQ results to converge with system size, and we will quantify this statement below.In χ(T ) [Fig.9(b)], the two methods achieve full convergence for all N and D, despite the lowering effective temperature scales, and we draw attention to the fact that previous ED and QMC studies at this coupling ratio could not accurately capture this low-temperature behavior; we revisit this point in the context of SrCu 2 (BO 3 ) 2 in Sec.VC.
Proceeding yet closer to the QPT, Fig. 10 shows the specific heat and susceptibility at the coupling ratio J/J D = 0.65.Clearly the "peak-dip-hump" shape of C(T ) has become significantly more pronounced in all respects than at J/J D = 0.63; the peak is sharper, the dip is deeper, and the separation of the peak and the broad hump has increased, with the peak moving below T /J D = 0.09 and the high-T maximum centered around T /J D = 0.5.Once again the progression from our ED data at N = 28 to TPQ with N = 32 and then 36 shows an increasing ability to capture the tip of the peak, but it cannot be said that convergence has been reached [inset, Fig. 10(a)].As at J/J D = 0.63, our iPEPS data for D = 16 and 18 are in essence indistinguishable and do suggest the peak to which the TPQ results are converging.At all other parts of the C(T ) curve, convergence is complete by all methods.In the corresponding susceptibility, again the ED and TPQ size sequence provides a systematic convergence towards the "squareness" of the turnover from rapid rise to broad maximum, while the higher iPEPS bond dimensions have converged to the shape of this feature.We take the coupling ratio J/J D = 0.66, shown in Fig. 11, as the upper limit to the validity of comparing our methods.Physically, it is obvious that the peakdip-hump shape of C(T ) becomes yet more pronounced, that the dip is deeper, and that the peak is even lowerlying and sharper than at J/J D = 0.65 [Fig.11(a)].Numerically, it appears that both our methods show slower convergence, not only around the peak but also in the dip.In more detail, our N = 32 and 36 TPQ results, and our D = 16 and 18 iPEPS results, do both offer internally consistent peak shapes.However, it is no longer clear that TPQ and iPEPS might extrapolate to the same limit, although this could be because the effective location of the QPT differs for the two methods due to finite-D and -N corrections, implying a closer proximity to a competing phase in one than in the other.On the overall shape of C(T ), we comment that the hump can be under-stood as the energy scale of local spin-flipping processes on the scale of J D , while the low-T peak is the focus of our discussions in Sec.VA.Concluding with the susceptibility at J/J D = 0.66 [Fig.11(b)], again our iPEPS results offer a consistent picture of the very steep onset and square maximum, while our TPQ results suggest better agreement with iPEPS at N = 32 than at N = 36.
The results in Figs. 8 to 11 capture the thermodynamic properties of the Shastry-Sutherland model with a degree of precision that was entirely unattainable by all previous techniques.It is safe to say that no prior numerical studies had even identified the peak-dip-hump form of the purely magnetic specific heat with any reliability, let alone given an indication of the sharpness of the low-T peak in the regime near the QPT.While the square shape of χ(T ) has been known since the early measurements on SrCu 2 (BO 3 ) 2 [8,12], again no numerical approaches had reproduced this form.We discuss both the physics underlying these features and the experimental comparison in Sec.V.
The thermodynamic properties of the Shastry-Sutherland model nevertheless set an extremely challenging problem in the regime close to the QPT, even in the dimer-product phase where the ground state is exact.It is clear in Figs. 8 to 11 that neither our TPQ nor our iPEPS calculations can be judged to be fully convergent, and thus accurate, as J/J D approaches the presumed critical value of 0.675.Because both are limited by their control parameters, N for TPQ and D for iPEPS, an optimal reproduction of the exact thermodynamics would be obtained by extrapolating both sets of results to infinite N or D.
By inspection of the full peak shape in C(T ), our TPQ results can only be said to have reached convergence at N ≤ 40 for J/J D ≤ 0.63.Similarly, it is not clear that our iPEPS results for the largest available D values (we have performed calculations with D = 18 at all temperatures and with D = 20 for temperatures around the peak) can be said to have converged to the shape of the peak at J/J D = 0.66.For extrapolation purposes, we first extract the primary characteristic features of this low-T peak, namely its position and height.In Fig. 12 we show the convergence of both quantities as functions respectively of 1/N in our ED and TPQ calculations and of 1/D in our iPEPS calculations, for each of the coupling ratios J/J D = 0.60, 0.63, 0.65, and 0.66.
Focusing first on our iPEPS results, it is evident that the values we have obtained for both the position and the height of the C(T ) peak show an excellent linear convergence with 1/D for all coupling ratios.Both quantities decrease with increasing D and in fact their slopes are remarkably insensitive to J/J D .Thus we estimate the putative infinite-D limit of the peak characteristics by linear extrapolation in 1/D (shown by a filled triangle) and we take the difference of this limit from the result at the largest computed D value as an estimate of the error (shown by a one-sided error bar).On the one hand, all iPEPS estimates decrease monotonically with increasing D. On the other hand, one ultimately expects exponential convergence in D rather than linear in 1/D.Consequently, we expect the true D = ∞ limit to lie somewhere inside the interval that is indicated.
Turning to our ED and TPQ results in Fig. 12, the peak heights at all coupling ratios show a quite systematic convergence, from below and with changing slopes, towards values fully consistent with our iPEPS estimates.By contrast, the finite-size evolution of the peak positions is less clear, but it is true to state for every coupling ratio that their values for all system sizes are either clustered closely around the extrapolated iPEPS values or are trending towards these.
Thus we may conclude that our TPQ and iPEPS results are consistent and convergent for even the most challenging features of the most challenging region of the Shastry-Sutherland phase diagram.As such they may be used as a reliable statement of the excitation spectrum as a function of the coupling ratio, which is the topic of Sec.VA.We discuss the implications of the convergence details in Fig. 12 for future numerical development in Sec.VB.

A. Spectrum of the Shastry-Sutherland Model
In Fig. 2(a) we showed our collected results for the specific heat of the Shastry-Sutherland model across the range of coupling ratios 0.60 ≤ J/J D ≤ 0.66.As the system approaches the QPT from the dimer-product to the plaquette state, there is a systematic emergence of a sharp low-energy peak, followed at higher temperatures by a deepening dip before a second broad maximum on the scale of J D /2.The energy scale of the concentration of low-lying states contributing to the peak falls steeply as the QPT is approached.From the corresponding susceptibilities, shown in Fig. 2(b), it is clear that the decreasing energy scale is largely similar for triplet (and higher spinful) excitations, and hence at minimum that it is not a special property of singlets; the physics of the approach to the QPT seems to affect all excitations in the same general manner.
Before discussing the nature of the low-lying spectrum, we characterize the height and temperature scale of the specific-heat peak, and the temperature scale of the susceptibility onset, as the coupling ratio approaches the QPT.The peak position, T C max , shown in Fig. 13(a), clearly falls faster than linearly in J/J D , but does not approach zero at the QPT.In general form, the peak position as a function of J/J D is somewhat reminiscent of the behavior of the singlet gap of the system, with a constant of proportionality of approximately 1/4 [Fig.13(a)].The fact that these falling energy scales do not vanish at the QPT is consistent with the first-order nature of this transition [5].The behavior of the peak height [Fig.13(b)] is less evident, in that it falls on approach to the QPT in our TPQ calculations but increases in our iPEPS ones, the rate of change accelerating with J/J D in both cases.Given the strong finite-size effects visible in our ED and TPQ calculations for J/J D > 0.63 [Fig.12], here it appears that our iPEPS results are more representative of the large-N limit.The final piece of information one may wish to extract about the peak is its width.Because of the dip and hump on the high-T side, we take a characteristic width from the quantity T C max − T C 1/2 determined on its low side [Fig.13(a)].We observe that the half-height and full-height positions scale to good approximation in the same way with J/J D , separated only by a constant factor, and hence that all the peak shapes are actually self-similar on the low-temperature side, with only one characteristic J/J D -dependent energy scale.
Because the maximum of χ(T ) is broad and flat, the determination of T χ max is delicate, and for an accurate characterization of the susceptibility we rather take the temperature, T χ 1/2 [76], at which χ(T ) reaches half its maximum value during its rapid onset.Figure 13(c) shows again QMC, ED, TPQ, and extrapolated iPEPS results for T χ max and T χ 1/2 .We note that finite-D corrections become comparatively important around T χ 1/2 and 13. (a) Position, T C max , (b) height, Cmax, of the peak in the specific heat, and (c) position T χ max of the peak of the susceptibility of Shastry-Sutherland model in the regime of coupling ratios 0.50 ≤ J/JD ≤ 0.66, as determined from QMC, ED, TPQ, and iPEPS calculations.In the case of iPEPS, triangles designate 1/D-extrapolations and one-sided error bars the difference to the largest available D. Also shown in panels (a) and (c) are the half-heights, T C 1/2 and T χ 1/2 on the low side of the peaks, from which one may deduce an effective peak width (shaded gray region connecting the extrapolated iPEPS data).Shown for reference in panels (a) and (c) are the functions ∆S/4 and ∆T /4, respectively, where ∆S and ∆T are respectively the gaps of the lowest singlet and triplet excitations obtained for a system of size N = 36.
as J/J D approaches 0.66 such that the corresponding extrapolated iPEPS results are subject to relatively large error bars.Within these error bars, Low-lying energy spectrum of the Shastry-Sutherland model in the S = 0 and S = 1 sectors, computed by ED using clusters of sizes N = 28 (to Ẽ = 2.6JD), N = 24 (to 3.75JD), and N = 20 (to 4.2JD) and shown over the full range of coupling ratios, J/JD, in the dimer-product phase.The excitation energy is denoted by Ẽ = E − EGS, where EGS is the energy of the ground state.The shading represents regions with discrete but dense singlet (blue) and triplet excitations (pink).The lowest pink line shows the one-triplon sector, whose lower edge is ∆T (J/JD); the yellow line shows the two-triplon scattering threshold, 2∆T (J/JD).The lower edge of the lower blue sector, which consists of singlet bound states of two triplons, is ∆S(J/JD).We draw attention to how few states in the three-triplon sector fall below the green line (∆S + ∆T ) and in the four-triplon sector below the black line (2∆S).0.66.Again within error bars, T χ 1/2 /T χ max is also indistinguishable from a constant such that the shape of the low-temperature side of the peak of the magnetic susceptibility χ is even less sensitive to J/J D in the region 0.5 ≤ J/J D ≤ 0.66 than the peak of the specific heat C. In this case the behavior of T χ 1/2 as a function of J/J D is correlated most closely with the triplet gap ∆ T , again with a constant of proportionality close to 1/4, as shown by the solid orange line in Fig. 13(c).
Figure 13 is a very specific diagnostic for the nature of the spectrum of low-lying energy levels, reflecting in particular a large number of states that become more nearly degenerate and are characterized by an ever-smaller (but always finite) gap as the system approaches the QPT.In Ref. [25] we showed a schematic illustration of the energy spectrum of the Shastry-Sutherland model as a function of the coupling ratio J/J D , based on Lanczos ED calculations to obtain the lowest levels of a cluster of size N = 36 sites.For more specific insight into the nature of the states revealed by the Lanczos procedure, here we have computed a much larger number of low-energy states, but for smaller clusters (of sizes up to N = 28).In Fig. 14 we gather these results to show the low-lying spectrum of the system in the S = 0 and S = 1 sectors over the full range of coupling ratios.We stress that, given the large number of states in these calculations, it is not easy to observe the discrete structure of the spectrum even in its lowest-lying regions, and as in Ref. [25] we use shading to indicate the support of these excitations.
In Fig. 14, the thin line at the bottom of the triplet spectrum is the elementary "triplon" excitation, which has energy Ẽ = J D at J = 0 and remains extremely narrow (non-dispersive) throughout the dimer-product phase due to the perfect frustration (Fig. 1).It is clear that significant numbers of two-and other multi-triplon states, in both the spin sectors shown, move to very low energies as the system approaches the QPT.Because the lowest states in the singlet sector must be of twotriplon origin, and these cross below the one-triplon energy around J/J D = 0.60 (the exact value depends on N ), it is equally clear that the physics of the Shastry-Sutherland model in the regime near the QPT is dominated by states one may classify as "strongly bound."To specify the meaning of this term, the yellow line shows the quantity 2∆ T (J/J D ), i.e. twice the singleparticle gap, which corresponds to the two-triplon scattering threshold, and near the QPT many states lie significantly below this energy.We comment that the separate "families" of excitations that can be traced back to independent two-triplon states at energy Ẽ = 2J D at J = 0, and similarly for higher triplon numbers [Fig.14], exhibit a much larger dispersion across the Brillouin zone than do the one-triplon states, due to the smaller effects of geometric frustration on the propagation of bound states [4,77,78].Although some authors have performed perturbative studies of the lowest-lying bound states of the Shastry-Sutherland model [77], we are not aware of previous attempts to count them systematically for comparison with thermodynamic properties, as we do next.
Beyond illustrating significant binding-energy effects in the two-and higher-triplon sectors, Fig. 14 does not provide further information about the structure of the spectrum.As a step in this direction, Fig. 15 shows the numbers of singlet and triplet states lying below the threshold given by the yellow line in Fig. 14.The number of triplet states excludes the N/2 states of the onetriplon band.For comparison, the expected number of two-triplon scattering states, 1  4 N ( 1 2 N − 1), is shown by the green line.Clearly the numbers of "bound" states in both sectors are significant fractions of the number of scattering states throughout the region J/J D 0.55.Further, the number of S = 1 states is comparable to that in the S = 0 sector, despite the fact that the energy gain is greater in the latter.As the QPT is approached, specifically once J/J D > 0.62 for N = 28 [Fig.15], the numbers of low-lying states in both sectors increase well beyond the number of independent two-triplon states.
Before discussing the nature of these states, for further information about the low-lying spectra we show in Fig. 16 the total numbers of singlet (red) and of triplet states (blue) lying below a given energy, Ẽ, in ED calculations performed with four different system sizes for three fixed values of J/J D in the interval near the QPT.Here the N/2 one-triplon excitations are included in the triplet count, appearing as the plateau region in all the blue lines directly above their onset.We note that all excitations with higher spin, S ≥ 2, lie above the two-triplon scattering threshold for all J/J D values shown; while binding effects have also been observed in the S = 2 sector [79], these are extremely small and are not detectable in clusters up to N = 28 due to finite-size effects.
The evolution of this spectrum as a function of the coupling ratio contains two new messages beyond the ever-increasing number of very low-lying states shown in Fig. 15.First, the steadily decreasing gap in the onetriplon sector can be expected to control the response to ∆S = 0 processes over most of regime 0.60 ≤ J/J D ≤ 0.66, but very near the QPT it is overwhelmed by the high density of other triplet states that finally fall to this gap.Second, the dominant feature of the spectrum in this regime is the large number of low-lying singlet states, which are governed by an energy scale that in turn falls significantly more quickly with J/J D than the triplet gap.
To discuss the nature of these states, and hence the consequences of our calculated spectra for the thermodynamic response functions computed in Sec.IV, it is helpful to compare and contrast the Shastry-Sutherland model with the fully frustrated two-leg spin ladder.In this 1D model, some of us [76,80] also found a striking thermodynamic response at low temperatures, which could be traced to a proliferation of low-energy multimagnon bound states upon approaching a first-order QPT.While it is clear that the proliferation of low-lying states is the origin of the low-temperature properties of the Shastry-Sutherland model, there are some important differences between the 1D and 2D systems, most notably concerning the mechanism behind the formation of these states.
First, in the fully frustrated two-leg ladder, it was clear that these states were truly bound, meaning that they were intrinsic combinations of n contributing triplons, with n ≥ 2 becoming very large near the QPT.By contrast, in the Shastry-Sutherland system one may identify essentially all of the additional low-lying states, whose origin clearly lies in n ≥ 3 triplons, with the expected scattering states of two singlet two-triplon bound states or of a singlet bound state and an elementary triplon.In fact it is clear that when the energy of a singlet bound state of two triplons falls beyond the one-triplon gap, the lowest-lying scattering states of a pair of singlet bound states must fall below the two-triplon scattering threshold in the S = 0 sector, as shown by the black line in Fig. 14, and so does the scattering state of a singlet bound state and an elementary triplon in the S = 1 sector, shown by the green line in the same figure.Thus the proliferation of low-lying excitations in the regime J/J D > 0.6 can, due to the differences in state-counting in 2D, be driven (almost) completely by the strength of the two-triplon binding effect in the singlet sector.It is for this reason that we do not refer to all of the low-lying states we observe here as "bound states."Concerning the possibility of true multi-triplon binding effects, as noted above we may comment only that, if present, these are sufficiently weak that we are unable to detect them.
Second, at the QPT in the fully frustrated ladder, the bound states cluster strongly at a single energy value of approximately 85% of the one-triplon gap [76,80], but in the Shastry-Sutherland model we observe no such sharp structures in the density of states.Specifically, although we find for J/J D = 0.66 that low-lying singlet states far outnumber the triplets at any energy Ẽ < 0.5J D [Fig.16(c)], their number continues to rise with Ẽ, at an increasing rate, and shows no sign of a discrete set of degenerate levels that could be related directly to the specific-heat peak.While these states are distributed over a relatively broad range of energy, their distinguishing feature, at least in the singlet sector, is that their binding effects are energetically stronger than in the 1D system.Although finite-size effects prevent us from reaching the QPT (J/J D = 0.675), at J/J D = 0.66 the singlet two-triplon bound states have already descended to only 54% of the one-triplon gap (on the N = 28 cluster).Thus the low-temperature peak in C(T ) observed in Sec.IV emerges at a relatively lower energy scale than in the fully frustrated ladder, and we suggest that it is this relative location that contributes to the anomalous sharpness of the peak.We note that the importance of only one energy scale is consistent with our discovery in Fig. 13(a) that the peak shapes are identical.Regarding the susceptibility, the decrease of the one-triplon energy scale in Fig. 16 appears to provide a reasonable account of the decrease in the "onset" temperature indicator T χ 1/2 [Fig.13(c)] over the entire range of J/J D , while the steep fall in energy of the other low-lying triplet states may be responsible for the increasingly abrupt ("square") turnover at the maximum in χ(T ).
In summary, our detailed ED investigation of the lowenergy spectrum [Figs.14, 15, and 16] reveals the rapid descent in energy of a multitude of singlet and triplet states as the QPT is approached.Despite the narrow nature of the C(T ) peak, these states do not converge to a single characteristic energy in either sector.In contrast to the n-triplon bound states of the fully frustrated ladder near its QPT, the energetics of the low-lying states of the Shastry-Sutherland model are dominated by the binding energy of two triplons into a net singlet.Because this local bound state involves two neighboring orthogonal dimers [77,78], it has little or no relation with the plaquette singlets that control the physics on the other side of the QPT (Fig. 1), and hence one may conclude that the first-order transition involves a complete rearrangement of spin correlations in the system.

B. Numerical Methods
We comment only briefly on the limitations and prospects for improvement of our numerical techniques.In a sense the Shastry-Sutherland model near the QPT presents an ideal test-bed for any method, in that the low-T peak in the specific heat becomes progressively narrower, and thus more challenging to reproduce, as the transition is approached.In our ED and TPQ calculations, it is both clear and completely unsurprising that larger cluster sizes improve the ability of both methods to capture this peak (and the corresponding susceptibility shoulder).Whereas the thermodynamic response at coupling ratio J/J D = 0.60 could not be reproduced with complete accuracy by N = 20 ED, or indeed N = 28 Lanczos ED (Fig. 8), we may assert that J/J D = 0.63 can be treated accurately by TPQ with N ≤ 40 (Fig. 9).
However, the limits to TPQ appear to be reached at J/J D > 0.65, where our calculations no longer reproduce the extremely sharp C(T ) peak that forms in this regime (Fig. 11).With reference to our discussion of Sec.VA, one may only speculate that the physics of this peak resides in scattering states of the singlet two-triplon bound states whose fully developed real-space extent begins to exceed this cluster size very close to the QPT.As noted in Sec.II, reaching larger cluster sizes with acceptable statistical errors is not possible with present computing power.Similarly TPQ, like ED, is not limited by algorithmic complexity, and indeed our present calculations already exploit all available spin and spatial (cluster) symmetries, implying that more complex Hamiltonians would be subject to stricter limits on N .
By contrast, our iPEPS calculations have relatively modest CPU and memory requirements.To date they have been run only on a single node, and thus larger-scale computations would certainly become feasible through the future development of a highly parallelized code.In comparison with more mature numerical methods, tensor-network calculations remain in their relative infancy, and as a result retain significant scope for algorithmic improvement.In current implementations, both CPU and memory requirements increase as high powers of the bond dimension, D; although most developments focus on reducing these exponents, it is not proven beyond doubt that large D is the only route to an accurate tensor-network description.As one example, we comment that higher accuracies can be achieved for the same D value within our present calculations by the use of disentangling gates acting on the ancillas [36,81].
Nevertheless, as demonstrated in Sec.III, the influence of the other variables in the iPEPS method (Trotter step, CTM boundary bond dimension, type of update) is small in comparison with the effects of D. Physically, there is at present no known means of relating D as a figure of merit in finite-T (i.e.thermal-state purification) iPEPS implementations with D in a ground-state iPEPS calculation.As with the problem of which absolute value of D is able to capture the physics of any given Hamiltonian, the most reliable approach to date is the empirical one of increasing D and monitoring changes and convergence.As shown in Secs.III and IV (Fig. 12), for the present problem we have indeed reached a well-defined limiting regime of D, even at J/J D = 0.66.Thus one may conclude that the Shastry-Sutherland model does not present a physical system where calculations are limited by the ability to capture long-ranged entanglement, as may perhaps be expected when the ground state is a product state, all the excitations are gapped, and the nearby QPT is first-order.

C. Experiment
In the light of the data and insight obtained from applying our two advanced numerical methods to the thermodynamics of the Shastry-Sutherland model, we revisit the experimental situation in SrCu 2 (BO 3 ) 2 .We note that low-lying singlet excitations in SrCu 2 (BO 3 ) 2 were documented by Raman scattering soon after the discovery of the material [82], and later that the same method was used to study the low-lying spectrum in an applied magnetic field [83].The triplet excitations were studied at first by electron spin resonance (ESR) [84] and inelastic neutron scattering [10].Subsequent high-field ESR measurements were used to demonstrate the presence of two-triplon bound states [85], while a recent application of modern neutron spectroscopy instrumentation studied the lowest bound triplet mode, and its response in a magnetic field, in unprecedented detail [86].The evolution of the triplet spectrum has been studied under pressure [23], which as noted in Sec.I acts to increase the ratio J/J D through the first-order QPT out of the dimerproduct phase.Although an analogous study of the thermodynamic properties has appeared in parallel with the completion of the present study [24], only one dataset for C(T )/T is shown at a pressure between ambient and the QPT; the low-T peak moves down by approximately 25% at this pressure and undergoes a definite increase in sharpness.Around the QPT, some extremely sharp C(T )/T peaks are found in some, but not all, samples.We observe that the numerical modeling that accompa-nies these measurements does not advance the state of the art beyond that of Ref. [26].
We comment first that a rather accurate understanding of the high-temperature, and indeed high-field, properties of SrCu 2 (BO 3 ) 2 can be obtained using the Shastry-Sutherland model with a coupling ratio close to J/J D = 0.63, and it is the remaining mystery surrounding the low-T response that we address here.We use our D = 18 iPEPS results as a consistent indicator of the evolution of the specific heat and susceptibility with the coupling ratio, and in Fig. 17 we show both quantities for all J/J D = 0.60, 0.61, . . ., 0.66.For the specific heat [Fig.17(a)] we compare our results with the data of Ref. [13] for SrCu 2 (BO 3 ) 2 , from which we have subtracted a phonon contribution C ph = γ T 3 with γ = 0.5 mJ/(mol K 4 ).For the susceptibility [Fig.17(b)] we compare our results with the data of Ref. [12], using a g-factor of g c = 2.28 [84].The optimal values of J D required to reproduce the position and shape of the low-T peak in C(T ) and to reproduce the rapid upturn and maximum height of χ(T ) are determined separately; in both cases they turn out to be strong functions of the coupling ratio, and hence can be used for accurate estimation of this energy scale.
In Fig. 17(a) it is safe to say that the best account of the low-T peak in C(T ) is given by J/J D = 0.62, with J D = 77 K. Somewhat surprisingly, it is not necessary to exploit to its limits our abilities to resolve a very narrow peak in C(T ), as would be the case at higher coupling ratios.Because J D is constrained by many other aspects of the experimental data for SrCu 2 (BO 3 ) 2 , its strong dependence on J/J D serves as a strict constraint on the coupling ratio.Turning to Fig. 17(b), the best account of χ(T ) is given by J/J D = 0.64, with J D = 91 K.Here it is safe to say that the fits in both panels of Fig. 17 are difficult to reconcile, certainly at the experimentally determined value of g c , and as a result we can conclude that we have found the limits of the Shastry-Sutherland model when applied to the SrCu 2 (BO 3 ) 2 problem.It is known that the appropriate spin Hamiltonian for this material contains additional terms, most notably an interlayer coupling, whose magnitude has been estimated at 9% of J D [4], and Dzyaloshinskii-Moriya interactions both within and normal to the plane of the system, whose magnitudes have been estimated [85] for both components at 3% of J D .Because it is not the purpose of our present analysis to address these details, for further comparison in the present context we restrict our attention to establishing the optimal Shastry-Sutherland fit for SrCu 2 (BO 3 ) 2 , to which end we adopt the compromise coupling ratio J/J D = 0.63 and find that the corresponding optimal J D is 87 K.
In Fig. 18 we show a quantitative comparison between the results of Sec.IV, specifically Fig. 9, and the C(T ) data of Ref. [13].We observe that our best TPQ results (N = 40) and our iPEPS results separately provide quantitatively excellent fits to the measured data.Our conclusion that such a fit is optimized with the coupling states and infinite projected entangled-pair states (iPEPS), for computing the thermodynamic properties of quantum magnets.Both approaches are unaffected by frustration and we apply them to compute the magnetic specific heat and susceptibility of the Shastry-Sutherland model near the quantum phase transition (QPT).This challenging region of coupling ratios has remained entirely impervious to meaningful analysis by any previous numerical techniques, and now we are able to reveal why.The specific heat develops a sharp low-temperature peak, which is separated from the broad hump due to local processes by an emerging dip.This peak is caused by a proliferation of low-lying energy levels whose origin lies in scattering states of multiple two-triplon bound states, and its anomalously low energy scale moves ever lower (albeit not to zero) as the QPT is approached.The analogous feature in the magnetic susceptibility is an increasingly rapid rise at an ever-lower temperature, followed by an abrupt maximum.This physics is a consequence of both high frustration and the proximity to a QPT, and its manifestations have been discussed previously in highly frustrated 1D systems, but now we have captured its realization in 2D.
Our methods and benchmarks offer very wide scope for immediate application.The list of open problems in frustrated quantum magnetism is long, and we await with interest definitive results, meaning for extended or spatially infinite systems, for the thermodynamic response of the triangular, kagome, and J 1 -J 2 square lattices.The ability to fingerprint the energy spectrum at the level of thermodynamic properties will provide a new dimension of physical understanding, and in some cases may even yield new insight into the ground state.Away from the classic problems of geometrical frustration in Heisenberg models, there is a pressing need for thermodynamic information to characterize many other types of frustrated system, including spin-ice materials, candidate Kitaev materials, coupled anisotropic (Ising and XY) chain and planar materials, candidate chiral spin liquids, materials with Dzyaloshinskii-Moriya interactions, topological magnets, and certain skyrmion systems.
Technically, neither of our methods presents any barrier to the application of a magnetic field, which could be used to control the relative singlet and triplet gaps and create qualitative changes in the thermodynamic response in many of the above types of system.For the SrCu 2 (BO 3 ) 2 problem, the next step enabled by the capabilities developed here is clearly to follow the evolution of the system to and through the pressure-induced QPT into the plaquette phase [23,24].This process presents a clear need for the investigation of thermodynamic properties, even if the physics of the material beyond the dimer-product phase is critically dependent on additional Hamiltonian terms beyond the Shastry-Sutherland model.

65 J
FIG. 2. (a)Magnetic specific heat, C(T ), and (b) magnetic susceptibility, χ(T ), of the Shastry-Sutherland model in the regime of coupling ratios 0.60 ≤ J/JD ≤ 0.66.The number of sites in the TPQ calculations is N = 36 and the shaded regions show the standard error of the TPQ method for the respective coupling ratios (Sec.IV).The iPEPS bond dimension is D = 18.

b
< l a t e x i t s h a 1 _ b a s e 6 4 = " n K 3 d i u m 9 c O 0 v N v c m N D J y r 3 1 p 8 0 E = " > A A A B 6 H i c b V B N S 8 N A E J 3 4 W e t X 1 a O X x S J 4 K k k t 6 L H g x W M L 9 g P a U D b b S b t 2 s w m 7 G 6 G E / g I v H h T x 6 k / y 5 r 9 x 2 + a g r Q 8 G H u / N M D M v S A T X x n W / n Y 3 N r e 2 d 3 c J e c f / g 8 4 a 0 A I G C M / w C m / O o / P i v D s f y 9 Y N J 5 8 5 g z 9 w P n 8 A 2 j m M 7 A = = < / l a t e x i t > l < l a t e x i t s h a 1 _ b a s e 6 4 = " P h 8 u s h d n e 6 m b + 7 I b o o 4 a 0 A I G C M / w C m / O o / P i v D s f y 9 Y N J 5 8 5 g z 9 w P n 8 A 0 S G M 5 g = = < / l a t e x i t > p < l a t e x i t s h a 1 _ b a s e 6 4 = " T Q 6 S p B 9 Y K W e C J G b 5 m U h H b b 6 4 a 0 A I G C M / w C m / O o / P i v D s f y 9 Y N J 5 8 5 g z 9 w P n 8 A 2 j m M 7 A = = < / l a t e x i t > l < l a t e x i t s h a 1 _ b a s e 6 4 = " P h 8 u s h d n e 6 m b + 7 I b o o 4 a 0 A I G C M / w C m / O o / P i v D s f y 9 Y N J 5 8 5 g z 9 w P n 8 A w H W M 2 w = = < / l a t e x i t > t e x i t s h a 1 _ b a s e 6 4 = " R w v g B W 0 9 7 n X t t m U g 7 x z J J Z K j M n Q = " > A A A B 6 H i c b V B N S 8 N A E J 3 4 W e t X 1 a O X x S J 4 K k k t 6 L H g x W M L 9 g P a U D b b T b t 2 s w m 7 E 6 G E / g I v H h T x 6 k / y 5 r 9 x 2 + a g r Q 8 G H u / N M D M v S K Q w 6 L r f z s b m 1 v b O b m G v u H 9 w e H R c O j l t m z j V j L d Y L G P d D a j h U i j e Q o G S d x P N a R R I 3 g k m d 3 O / 8 8 S 1 E b F 6 w G n C / Y i O l A g F o 2 i l J g 5 K Z b f i L k D W i Z e T M u R o D E p f / W H M 0 o g r Z J I a 0 / P c B P 2 M a h R M 8 l m x n x q e U D a h I 9 6 z V N G I G z 9 b H D o j l 1 Y Z k j D W t h S S h f p 7 r l S b t X K 9 l s d R g H O 4 g C v w 4 A b q c A 8 N a A E D D s / w C m / O o / P i v D s f y 9 Y N J 5 8 5 g z 9 w P n 8 A 3 U G M 7 g = = < / l a t e x i t > b < l a t e x i t s h a 1 _ b a s e 6 4 = " n K 3 d i u m 9 c O 0 v N v c m N D J y r 3 1 p 8 0 E = " > A A A B 6 H i c b V B N S 8 N A E J 3 4 W e t X 1 a O X x S J 4 K k k t 6 L H g x W M L 9 g P a U D b b S b t 2 s w m 7 G 6 G E / g I v H h T x 6 k / y 5 r 9 x 2 + a g r Q 8 G H u / N M D M v S A T X x n W / n Y 3 N r e 2 d 3 c J e c f / g 8 FIG. 4. (a) Graphical representation of an iPEPS tensor, A rtlb p , in which each leg corresponds to an index of the tensor.The physical index, p (green leg), represents the local Hilbert space of a lattice site (here a dimer) and the other four, r, t, l, and b are auxiliary indices, each with bond dimension D, which connect neighboring tensors.(b) A standard iPEPS representing a pure state, shown here for a unit cell containing two different tensors, A and B, arranged in a checkerboard pattern on the infinite lattice.A connection between two tensors implies taking the sum over the corresponding index.(c) iPEPS tensor with an extra index, a (blue leg), representing the local Hilbert space of an ancilla site.(d) iPEPS ansatz for the purification, |Ψ(β) , of the thermal density operator, ρ(β).(e) Representation of ρ(β) = Tra |Ψ(β) Ψ(β)|, where the trace is taken over all ancilla degrees of freedom (shown by the connected blue lines).

FIG. 5 .
FIG. 5. Comparison of D = 14 iPEPS data for (a) the magnetic specific heat and (b) the magnetic susceptibility for the Shastry-Sutherland model with J/JD = 0.63.The reference line corresponds to a D = 14 simple-update calculation with time step τ = 0.04, boundary bond dimension χ = 50, and in panel (b) a field h/JD = 0.01.For each of the other lines, one parameter has been varied in the calculation with respect to the reference line, as indicated in the legend.The insets magnify the maxima and the minimum.

FIG. 6 .
FIG.6.iPEPS results for (a) the magnetic specific heat and (b) the magnetic susceptibility of the Shastry-Sutherland model with J/JD = 0.63, obtained by the simple-update method for different values of the bond dimension, D. The insets magnify the maxima and the minimum.

= 16 iPEPS D = 18 QMC N = 32 QMC N = 200 FIG. 7 .
FIG. 7. (a) Magnetic specific heat of the Shastry-Sutherland model at J/JD = 0.50, computed by TPQ with N = 32 and 36 and by iPEPS with D = 16 and 18. Shown for comparison are ED results with N = 20 and QMC results with N = 32 and 200.(b) Corresponding magnetic susceptibility.

= 32 FIG. 8 .
FIG. 8. (a) Magnetic specific heat of the Shastry-Sutherland model at J/JD = 0.60, computed by TPQ with N = 32 and 36 and by iPEPS with D = 14, 16, and 18. Shown for comparison are ED results with N = 20 and with N = 28 at temperatures T /JD ≤ 0.14, as well as QMC results with N = 32 at T /JD 0.3.(b) Corresponding magnetic susceptibility.

= 32 FIG. 9 .
FIG. 9. (a) Magnetic specific heat of the Shastry-Sutherland model at J/JD = 0.63, computed by TPQ with N = 32, 36, and 40 and by iPEPS with D = 14, 16, and 18. Shown for comparison are ED results with N = 20 and with N = 28 at temperatures T /JD ≤ 0.125, as well as QMC results with N = 32 at T /JD 0.35.(b) Corresponding magnetic susceptibility.

= 32 FIG. 10 .
FIG. 10.(a) Magnetic specific heat of the Shastry-Sutherland model at J/JD = 0.65, computed by TPQ with N = 32 and 36 and by iPEPS with D = 14, 16, and 18. Shown for comparison are ED results with N = 20 and with N = 28 at temperatures T /JD ≤ 0.11, as well as QMC results with N = 32 at T /JD 0.4.(b) Corresponding magnetic susceptibility.

= 18 FIG. 11 .
FIG. 11.(a) Magnetic specific heat of the Shastry-Sutherland model at J/JD = 0.66, computed by TPQ with N = 32 and 36 and by iPEPS with D = 14, 16, and 18. Shown for comparison are ED results with N = 20 and with N = 28 at temperatures T /JD ≤ 0.1.(b) Corresponding magnetic susceptibility.

FIG. 12 .
FIG. 12. Characteristic position and height of the low-T peak in the magnetic specific heat of the Shastry-Sutherland model, shown for our ED and TPQ results as functions of 1/N and for our iPEPS results as functions of 1/D.(a,b) J/JD = 0.60.(c,d) J/JD = 0.63.(e,f) J/JD = 0.65.(g,h) J/JD = 0.66.
FIG. 14.Low-lying energy spectrum of the Shastry-Sutherland model in the S = 0 and S = 1 sectors, computed by ED using clusters of sizes N = 28 (to Ẽ = 2.6JD), N = 24 (to 3.75JD), and N = 20 (to 4.2JD) and shown over the full range of coupling ratios, J/JD, in the dimer-product phase.The excitation energy is denoted by Ẽ = E − EGS, where EGS is the energy of the ground state.The shading represents regions with discrete but dense singlet (blue) and triplet excitations (pink).The lowest pink line shows the one-triplon sector, whose lower edge is ∆T (J/JD); the yellow line shows the two-triplon scattering threshold, 2∆T (J/JD).The lower edge of the lower blue sector, which consists of singlet bound states of two triplons, is ∆S(J/JD).We draw attention to how few states in the three-triplon sector fall below the green line (∆S + ∆T ) and in the four-triplon sector below the black line (2∆S).

1 FIG. 15 .
FIG.15.Numbers of low-lying singlet and triplet states on an N = 28 lattice, shown as functions of J/JD in the dimerproduct phase.Symbols show ED data while connecting lines are guides to the eye.The horizontal green line shows the number of independent two-triplon scattering states.

S = 1 ,FIG. 16 .
FIG.16.Total number of excited states lying below the energy Ẽ at coupling ratios (a) J/JD = 0.60, (b) J/JD = 0.63, and (c) J/JD = 0.66.Red lines indicate S = 0 excitations, blue lines S = 1.Increasing line widths indicate increasing system sizes, N , on which ED calculations were performed.Vertical green lines denote the two-triplon scattering thresholds for the corresponding system sizes.
FIG. 19.(a)Comparison of experimental measurements of the magnetic susceptibility, taken from Ref.[12], to numerical calculations performed by ED, TPQ, and D = 18 iPEPS for J/JD = 0.63.We use g = 2.28 and JD = 87 K. (b) Detail of the comparison at low temperatures.