Quantum Monte Carlo simulations of highly frustrated magnets in a cluster basis The two-dimensional Shastry-Sutherland model

. Quantum Monte Carlo (QMC) simulations constitute nowadays one of the most powerful methods to study strongly correlated quantum systems, provided that no “sign problem” arises. However, many systems of interest, including highly frustrated magnets, suﬀer from an average sign that is close to zero in standard QMC simulations. Nevertheless, a possible sign problem depends on the simulation basis, and here we demonstrate how a suitable choice of cluster basis can be used to eliminate or at least reduce the sign problem in highly frustrated magnets that were so far inaccessible to eﬃcient QMC simulations. We focus in particular on the application of a two-spin (dimer)-based QMC method to the thermodynamics of the spin-1/2 Shastry-Sutherland model for SrCu 2 (BO 3 ) 2 .


Introduction
Quantum many-body systems constitute a challenge for a quantitatively accurate treatment. In this context, Quantum Monte Carlo (QMC) methods have become one of the most powerful tools [1][2][3]. However, in some cases, including geometrically frustrated magnets, their use is limited by a phenomenon known as the "sign problem", namely a cancellation of positive and negative weights of configurations [4][5][6]. In the present contribution, we illustrate how these problems can be overcome or at least alleviated by working in a suitable cluster basis [7][8][9][10]. For the purpose of simplicity, we will focus on the application of a dimer basis to a highly frustrated two-dimensional spin-1/2 model known as the "Shastry-Sutherland model" [11] that is realised in the compound SrCu 2 (BO 3 ) 2 [12,13].

Method and the Shastry-Sutherland model
One flexible QMC approach is based on the representation of the partition function Z as a high-temperature series [14,15], leading to the "stochastic series expansion" (SSE). The SSE is based on a splitting of the Hamiltonian H into bond terms H b . These terms are split in turn into different types, such that the action of each bond term H (b,t) on a given basis state |α either vanishes or yields exactly one other basis state |α . Note that one of the types t is diagonal while all others are off-diagonal. Using H = (b,t) H (b,t) , one finds in the given basis where β = 1/T is the inverse temperature (we set the Boltzmann constant k B = 1). Important developments concern the representation and efficient Monte Carlo sampling of the sum in Eq. (1), as well as measurement of physical quantities, and we refer to Refs. [10,[16][17][18] for detailed discussions of these issues.
Here, our main concern will be a different one: the weights W c of a configuration c in Eq. (1) are given by products of matrix elements α i+1 | −H (b i ,t i ) |α i and there is in general no guarantee that they are non-negative, i.e., that the expression Eq. (1) does indeed have a probabilistic interpretation. For diagonal terms of the Hamiltonian, one can add a constant shift to render −H (b i ,t i ) + non-negative. When off-diagonal terms are negative, one uses the absolute value of the weights |W c | for sampling such that the average of an observable A is given as the ratio of the averages of the observable and of the sign [1,4], In order to be specific, we will now focus on the Shastry-Sutherland model, which is represented schematically in Fig. 1. The Hamiltonian in the singlesite basis is where S i are spin-1/2 operators for site i. J D is the intra-dimer coupling (denoted by i, j ) and the interdimer coupling ( i, j ) J corresponds to a square lattice, as shown in Fig. 1. This model was constructed [11] to ensure that the ground state is an exact product of singlets formed on the dimer bonds for small and intermediate coupling ratios J/J D . More recently, this model was also found to be realised in the quantum spin system SrCu 2 (BO 3 ) 2 [12,13], which increased the interest in calculating its thermodynamic properties. As a consequence of the triangular arrangement of the J and J D interactions, this problem is geometrically frustrated. Consequently, QMC simulations in the single-spin basis (3) suffer from a severe sign problem, as we will show in section 3. However, one may pass to the total-spin and spin-difference operators, respectively , for the two spins S d,1 and S d,2 on a given dimer d [19]. The coupling of the two spins in dimer d can then be expressed through T 2 d and the inter-dimer coupling for the two dimers d and d indicated in Fig. 1 takes the form The first term corresponds just to a composite-spin model on a square lattice, which is unfrustrated and thus poses no QMC sign problem. The sign of the coefficient of the second term in Eq. (4) depends on the assignment of the labels 1 and 2 in dimer d and varies for the different inter-dimer coupling terms throughout the lattice. There are also signs involved in the construction of the operator D d such that it is not immediately clear to which extent this second term may still give rise to a sign problem. Nevertheless, the operator D d yields dimer triplet states when applied to the singlet state of the dimer d, Because the ground state is an exact product of such dimer singlets for J J D , one may expect the difference operators D d not to contribute to the leading orders at low temperatures, such that a potential sign problem remains manageable. We will show in the next section that indeed this turns out to be true in the regime where J is not too large.

Results
We now present some QMC results for the Shastry-Sutherland model where we will focus on the two representative cases J/J D = 0.5 and 0.55. Figure 2 presents the average sign sign |·| . In the single-spin basis (open symbols), the geometric frustration leads to the expected sign problem (see section 2), which manifests itself as the rapid drop of the average sign at low temperatures. If a sign problem arises, the average sign generically becomes exponentially small in N/T [1,4]. This is verified by the single-site results, where the average sign not only drops at low temperatures but also decreases with increasing system size N . Working in the dimer basis improves the sign problem very significantly, as is illustrated by the full symbols in Fig. 2. There is still a sign problem, as evidenced by the regime in which sign |·| < 1. However, for N = 32 and 64 spins, the average sign is always bigger in the dimer basis than in the single-site basis. Most remarkably, at least in the case J/J D = 0.5 ( Fig. 2(a)) the average sign does not vanish as temperature goes to zero, but instead recovers after going through a minimum around T /J D ≈ 0.25. This minimum still drops with increasing N , but one can go to systems as large as N = 200 and still ensure that sign |·| 0.06, such that the resulting loss of accuracy can be compensated by improving the statistics.
At J/J D = 0.55 ( Fig. 2(b)), the average sign no longer recovers for T → 0, even in the dimer basis. This difference between the cases J/J D = 0.5 and 0.55 can be traced to an algorithmic phase transition at J/J D = 0.526(1) [19], where the effective model defined by the absolute value of the weights in Eq. (2) undergoes a transition out of the dimer phase. Nevertheless, sign |·| in the dimer basis for N = 200 spins retains a larger value than that in the single-site basis for the much smaller system N = 32 at temperatures T /J D 0.7, as shown in Fig. 2(b).
Next we present results for thermodynamic observables at J/J D = 0.5 and 0.55 which are shown respectively in Figs. 3 and 4. As before, QMC results in the single-spin basis are shown by open symbols and results obtained in the dimer basis by filled symbols. For reference, we have included exact diagonalisation (ED) results for N = 20 spins [19] and infinite projected entangled-pair states (iPEPS) data for a bond dimension D = 16 [20]. As thermodynamic observables we have selected the magnetic susceptibility χ and specific heat C; the temperature window 0 ≤ T ≤ J D is chosen to include the maximum of these two quantities.
The single-site basis is barely able to resolve the maximum of χ at J/J D = 0.5 ( Fig. 3(a)) or get close to it for J/J D = 0.55 (Fig. 4(a)), provided that one restricts a system size as small as N = 36. However, single-site QMC generally loses accuracy already at significantly higher temperatures, such that the maxima are out of reach and only the high-temperature regime is accessible, as was to be expected according to the average sign shown in Fig. 2. On the other hand, for J/J D = 0.5, the dimer basis allows one to obtain accurate QMC data for systems with up to N = 200 spins and at all temperatures, as shown in Fig. 3. For J/J D = 0.55, the sign problem of Fig. 2(b) prevents one from reaching the zero-temperature limit, but the maximum of the magnetic susceptibility can still be well resolved in the dimer basis for N ≤ 128 spins ( Fig. 4(a)). In order to resolve the maximum of C (Fig. 4(b)) by QMC in the dimer basis, one needs to compromise on system size, i.e., restrict to N = 36, and on accuracy, but even N = 128 still works significantly better in the dimer basis than N = 36 in the single-spin basis, again as expected from the results for the average sign shown in Fig. 2.
The iPEPS data for D = 16 can be considered as representative of the thermodynamic limit in the present parameter regime and at the level of resolution in Figs. 3 and 4 [20]. Indeed, all of our QMC results agree with iPEPS within the corresponding statistical error bars. On the other hand, small (albeit distinct) finite-size effects can be seen in the N = 20 ED data, in particular around the maximum of C and at temperatures above it. These finite-size effects become increasingly pronounced at the larger value J/J D = 0.55 ( Fig. 4(b)). This highlights the importance of having access to larger systems in order to ensure quantitatively accurate results.
Concerning the physics of these thermodynamic quantities, we observe the emergence of lowtemperature maxima in the magnetic susceptibility and the specific heat that move to lower temperature and sharpen as J/J D increases [19,20]. In the present results, this is seen most clearly for the specific heat, shown in Figs. 3(b) and 4(b).

Conclusion
We have shown how a dimer basis can be used to perform accurate QMC simulations for the thermodynamic properties of a highly frustrated two-dimensional spin-1/2 model, the Shastry-Sutherland model, that would be out of reach of conventional approaches. The compound SrCu 2 (BO 3 ) 2 is believed to correspond to the parameter ratio J/J D ≈ 0.63 [13,20]. This ratio is still in the dimer phase of the Shastry-Sutherland model, but rather far beyond the algorithmic phase transition at J/J D = 0.526(1) [19]. Thus the parameter regime relevant to SrCu 2 (BO 3 ) 2 remains unfortunately inaccessible to accurate QMC simulations even in the dimer basis, and one needs to resort to other methods such as iPEPS [20,21].
Still, there are other cases to which the ideas illustrated here can be applied. For example, in the case of a frustrated ladder, only boundary terms contribute to the sign problem, such that it remains sufficiently mild to allow for a detailed study of the entire phase diagram [8,22]. Another case is a fully frustrated bilayer model, where the counterpart of the T d · D d term in Eq. (4) is absent and thus the dimer basis eliminates the sign problem completely [9,23,24]. A detailed QMC study of the phase diagram then becomes possible, allowing one to trace a first-order phase transition from zero temperature to finite temperature until it terminates at an Ising critical point [21,24].
The cluster-basis approach is is not restricted to dimers and we have demonstrated recently that it also works with a trimer (three-spin) basis [10]. This permits one to generalise the investigation of the fully frustrated bilayer model to a trilayer analogue, and study for example the influence of the ground-state entropy that arises in the latter model on the slope of the first-order line at finite temperatures [10].