Dynamics and Instabilities of the Shastry-Sutherland Model

TL;DR

This paper investigates the excitation spectrum and phase transitions in the Shastry-Sutherland model using an unbiased variational method, revealing multiple instabilities and order parameters in the system.

Contribution

It introduces a variational approach that works in the thermodynamic limit to identify susceptibilities and instabilities in the Shastry-Sutherland model.

Findings

Identified four different instabilities: antiferro spin nematic, plaquette spin nematic, stripe magnetic order, and plaquette order.

Revealed the nature of phase transitions in the model.

Confirmed some previously reported phases and discovered new potential orderings.

Abstract

We study the excitation spectrum in the dimer phase of the Shastry-Sutherland model by using an unbiased variational method that works in the thermodynamic limit. The method outputs dynamical correlation functions in all possible channels. This output is exploited to identify the order parameters with the highest susceptibility (single or multitriplon condensation in a specific channel) upon approaching a quantum phase transition in the magnetic field versus the $J'/J$ phase diagram. We find four different instabilities: antiferro spin nematic, plaquette spin nematic, stripe magnetic order, and plaquette order, two of which have been reported in previous studies.

Figures

Tables

Table 1. Table 1: Chatacter table of group C 4 ​ v subscript 𝐶 4 𝑣 C_{4v}   Elliott and Dawber ( 1979 ) .

	$E$	$C_2$	$2C_4$	$2{\sigma }_v$	$2{\sigma }_d$
$A_1\left[1\right]$	1	1	1	1	1
$A_2\left[xy\left(x^2-y^2\right)\right]$	1	1	1	-1	-1
$B_1\left[x^2-y^2\right]$	1	1	-1	1	-1
$B_2\left[xy\right]$	1	1	-1	-1	1
$E[x,y]$	2	-2	0	0	0

Equations

$$\mathcal{H}=J\sum _{⟨ij⟩}{\mathbf{S}}_i\cdot {\mathbf{S}}_j+J^{\prime }\sum _{⟨⟨ij⟩⟩}{\mathbf{S}}_i\cdot {\mathbf{S}}_j-h\sum _i{S}_i^z,$$

$$|{\phi }_i(\mathbf{k})⟩\equiv \frac{{\hat{P}}_{\mathbf{k}}|{\phi }_i⟩}{\sqrt{⟨{\phi }_i|{\hat{P}}_{\mathbf{k}}|{\phi }_i⟩}},$$

$$S^{-+}(\mathbf{k},\omega )=2\pi \sum _{\nu }|⟨\nu |S_{\mathbf{k}}^{+}|0⟩{|}^2\delta (\omega +E_0-E_{\nu }),$$

$${\chi }_{S2E}(\mathbf{k},\omega )=2\pi \sum _{\nu }|⟨\nu |A_{\mathbf{k}}^{S2E}|0⟩{|}^2\delta (\omega +E_0-E_{\nu }).$$

$$|{\Psi }_1(\mathbf{Q}){⟩}_{S2E}$$

$$\approx \alpha \left[|{\phi }_7^{(2)}(\mathbf{Q})⟩+|{\phi }_8^{(2)}(\mathbf{Q})⟩\right]+\beta \left[|{\phi }_2^{(2)}(\mathbf{Q})⟩-|{\phi }_1^{(2)}(\mathbf{Q})⟩\right],$$

$$|{\Psi }_2(\mathbf{Q}){⟩}_{S2E}$$

$$\approx \alpha \left[|{\phi }_6^{(2)}(\mathbf{Q})⟩-|{\phi }_5^{(2)}(\mathbf{Q})⟩\right]+\beta \left[|{\phi }_4^{(2)}(\mathbf{Q})⟩-|{\phi }_3^{(2)}(\mathbf{Q})⟩\right],$$

$$|{\phi }_i^{(S_z)}(\mathbf{Q})⟩=\frac{{\hat{P}}_{\mathbf{Q}}|{\phi }_i^{(S_z)}⟩}{\sqrt{⟨{\phi }_i^{(S_z)}|{\hat{P}}_{\mathbf{Q}}|{\phi }_i^{(S_z)}⟩}}.$$

$$A_{\mathbf{Q}}^{S2E}\equiv |{\Psi }_1(\mathbf{Q}){⟩}_{S2E}⟨0|,.$$

$$|\Psi (\mathbf{K}){⟩}_{S2A_2}$$

$$=\alpha \left[|{\phi }_5^{(2)}(\mathbf{K})⟩-|{\phi }_6^{(2)}(\mathbf{K})⟩-|{\phi }_7^{(2)}(\mathbf{K})⟩+|{\phi }_8^{(2)}(\mathbf{K})⟩\right]$$

$$+\beta \left[-|{\phi }_1^{(2)}(\mathbf{K})⟩-|{\phi }_2^{(2)}(\mathbf{K})⟩+|{\phi }_3^{(2)}(\mathbf{K})⟩+|{\phi }_4^{(2)}(\mathbf{K})⟩\right].$$

$$|\Psi (\mathbf{K}){⟩}_{S0A_2}$$

$$\displaystyle=\alpha\Big{[}|\varphi_{1}^{(0)}(\bm{K})\rangle+|\varphi_{2}^{(0)}(\bm{K})\rangle+|\varphi_{3}^{(0)}(\bm{K})\rangle+|\varphi_{4}^{(0)}(\bm{K})\rangle$$

$$\displaystyle\quad-|\varphi_{5}^{(0)}(\bm{K})\rangle-|\varphi_{6}^{(0)}(\bm{K})\rangle-|\varphi_{7}^{(0)}(\bm{K})\rangle-|\varphi_{8}^{(0)}(\bm{K})\rangle\Big{]}$$

$$\displaystyle\quad-\beta\Big{[}|\varphi_{9}^{(0)}(\bm{K})\rangle+|\varphi_{10}^{(0)}(\bm{K})\rangle-|\varphi_{11}^{(0)}(\bm{K})\rangle-|\varphi_{12}^{(0)}(\bm{K})\rangle\Big{]}.$$

$${\Delta }_M(\mathbf{k})=\sqrt{\frac{1}{\Lambda }{\int }_0^{\Lambda }d\omega {\left(\frac{S^{-+}(\mathbf{k},\omega {)}_M-S^{-+}(\mathbf{k},\omega {)}_{M-1}}{S^{-+}(\mathbf{k},\omega {)}_M+S^{-+}(\mathbf{k},\omega {)}_{M-1}}\right)}^2},$$

Full text

The Shastry-Sutherland model (SSM) has become a paradigmatic Hamiltonian of frustrated quantum magnetism Shastry and Sutherland (1981); Miyahara and Ueda (2003) because it includes an exactly solvable ground state Shastry and Sutherland (1981), very heavy low-energy excitations Miyahara and Ueda (1999); Weihong et al. (1999); Fukumoto (2000); Knetter et al. (2000a, b); Totsuka et al. (2001), exotic phases obtained upon varying the ratio $J^{\prime}/J$ between two competing exchange constants Albrecht and Mila (1996); Weihong et al. (1999); Koga and Kawakami (2000); Koga (2000); Knetter et al. (2000b); Zheng et al. (2001); Läuchli et al. (2002); Munehisa and Munehisa (2003); Lou et al. ; Corboz and Mila (2013), and a series of magnetic field induced magnetization plateaux Miyahara and Ueda (1999); Müller-Hartmann et al. (2000); Momoi and Totsuka (2000a, b); Fukumoto and Oguchi (2000); Fukumoto (2001); Misguich et al. (2001); Miyahara et al. (2003); Dorier et al. (2008); Abendschein and Capponi (2008); Isaev et al. (2009); Takigawa et al. (2010); Nemec et al. (2012); Corboz and Mila (2014); Schneider et al. (2016); Takigawa and Mila (2011). Its realization in SrCu2(BO3)2 Kageyama et al. (1999); Miyahara and Ueda (1999); Ueda and Miyahara (1999) enabled various experimental studies, including magnetization Kageyama et al. (1999); Onizuka et al. (2000); Sebastian et al. (2008); Jaime et al. (2012); Takigawa et al. (2013); Matsuda et al. (2013); Haravifard et al. (2016), specific heat Tsujii et al. (2011), inelastic neutron scattering (NS) Kageyama et al. (2000); Cépas et al. (2001); Gaulin et al. (2004); Kakurai et al. (2005); Zayed et al. (2014); McClarty et al. (2017); Zayed et al. (2017), far-infrared Rõõm et al. (2000), electron spin resonance (ESR) Nojiri et al. (1999, 2003), Raman scattering Lemmens et al. (2000), and nuclear magnetic resonance (NMR) Kodama et al. (2002, 2005); Takigawa et al. (2013). These studies revealed that a finite Dzyaloshinskii-Moriya (DM) interaction Dzyaloshinsky (1958); Moriya (1960) must be added to the SSM in order to account for several properties of SrCu2(BO3)2 Nojiri et al. (1999); Cépas et al. (2001); Nojiri et al. (2003); Gaulin et al. (2004); Miyahara et al. (2004); Kodama et al. (2005); Kakurai et al. (2005); El Shawish et al. (2005); Cheng et al. (2007); Romhányi et al. (2011, 2015); McClarty et al. (2017).

Despite the great theoretical efforts devoted to the SSM, the problem is still far from being solved. Perturbative approaches are only applicable in narrow regimes and conventional numerical methods suffer from severe size effects. As a consequence, the nature of the quantum phase diagram has been debated for a long time Albrecht and Mila (1996); Weihong et al. (1999); Koga and Kawakami (2000); Koga (2000); Knetter et al. (2000b); Zheng et al. (2001); Läuchli et al. (2002); Munehisa and Munehisa (2003); Lou et al. ; Corboz and Mila (2013). It is thus desirable to develop and apply alternative approaches. The infinite projected entangled-pair states (iPEPS) is an example of an alternative approach that works in the thermodynamic limit Jordan et al. (2008); Corboz and Mila (2013, 2014). However, it relies heavily on the initial guess of the physical states and it is difficult to extract dynamical responses.

In this Letter, we introduce an unbiased variational method to calculate the excitation spectrum and dynamical responses (susceptibilities) of the SSM in the dimer phase 111A numerical method is said to be “variational” when it satisfies the variational principle. Other examples of unbiased variational methods include the density matrix renormalization group White (1992, 1993) and matrix product method Östlund and Rommer (1995); Rommer and Östlund (1997), etc.. The method works in the thermodynamic limit and it complements alternative approaches like iPEPS. The basic idea was originally introduced to compute the single-hole dispersion of the square lattice $t$-$J$ model Trugman (1988, 1990). The same method was applied to the Shastry-Sutherland lattice $t$-$J$ model El Shawish and Bonča (2006); Haravifard et al. (2006). Here we exploit this idea to compute dynamical correlators and dominant instabilities. By working in a reduced Hilbert space, which preserves all model symmetries, we obtain low energy excitations classified by quantum numbers. We then predict the character of the neighboring phases by detecting the order parameter with highest susceptibility. Besides confirming the previously reported plaquette order and antiferro spin-nematic phases, we find two new phases; namely, a plaquette spin-nematic phase and stripe magnetic ordering, induced by simultaneously increasing the magnetic field and $J^{\prime}/J$. In particular, the plaquette spin-nematic phase explains the nature of the two-triplon states (pinwheels) that crystallize at higher field values Corboz and Mila (2014).

We consider the spin-$\frac{1}{2}$ SSM under a magnetic field Shastry and Sutherland (1981):

[TABLE]

where $\langle ij\rangle$ and $\langle\langle ij\rangle\rangle$ denote intradimer and interdimer neighbors. The unit cell has 4 sites (see Fig. 1). The exact ground state for small enough $J^{\prime}/J$ and $h$ is a direct product of singlet states on all dimers Shastry and Sutherland (1981). The elementary excitations of this “dimer phase” are singlet-triplet excitations known as triplons. Triplons are dressed by quantum fluctuations with a magnetic correlation length $\xi$ that increases with $J^{\prime}/J$. Methods that can account for the spatial extent of these quantum fluctuations should provide a good description of the low-energy excitation spectrum of the dimer phase.

We start the process by creating local excited states $|\varphi_{i}\rangle$ (e.g., single and two triplons). We then project these representative states into subspaces with fixed momentum $\bm{k}$,

[TABLE]

where the projector is defined as $\hat{P}_{\bm{k}}\equiv\frac{1}{N}\sum_{\bm{r}}e^{i\bm{k}\cdot\bm{r}}\hat{T}(\bm{r})$. $N\rightarrow\infty$ (thermodynamic limit) is the total number of unit cells, and $\hat{T}(\bm{r})$ is the translation operator. Application of $\mathcal{H}$ to $|\varphi_{i}(\bm{k})\rangle$ generates new states that dress the corresponding quasiparticle excitation. This procedure can be applied iteratively to systematically improve the variational space. After $M$ iterations, we obtain a basis $\{\varphi_{i}(\bm{k})\}$ with good quantum numbers $\bm{k}$ and $S_{\text{tot}}^{z}$. The number of iterations determines the spatial range $l$ of the fluctuations that dress the quasiparticle, so the method is then expected to produce accurate results for $l\gtrsim\xi$.

The eigenvalues and eigenvectors of the Hamiltonian restricted to the variational space are obtained by applying the implicitly restarted Arnoldi method Lanczos (1950); Lehoucq et al. (1998). The eigenvectors are classified by the Little Group of $C_{4v}$ for each momentum $\bm{k}$. A continuous phase transition manifests via a vanishing gap (condensation) that signals a phase transition into a broken symmetry state. The symmetry of the new state is determined by the irreducible representation (IREP) of the eigenstate that becomes gapless. To keep the method unbiased, the initial basis must not break the point group symmetry of $\mathcal{H}$ 222This means that the basis generates an invariant subspace of the group. In addition, a fully symmetric initial basis leads to better convergence (as a function of $M$) in comparison with other choices Ref. El Shawish and Bonča (2006). .

For illustration, we first focus on the $S_{\text{tot}}^{z}=0$ sector relevant to $h=0$. We include $\mathcal{D}=14$ $S_{\text{tot}}^{z}=0$ initial states to start the iteration sup and then apply Eq. (1) onto this basis to systematically increase the variational space 333The computational cost is exponential in the linear size of the quasiparticle. In the present case, the dimension of the truncated Hilbert space roughly scales as $\mathcal{D}\propto 6^{M}\propto 6^{l}$. For the $S_{\text{tot}}^{z}=0$ sector, the dimension of the variational space is $\mathcal{D}=74$ for $M=1$, $\mathcal{D}=396$ for $M=2$, $\ldots$, $\mathcal{D}=17\,730\,330$ for $M=8$.. After obtaining the lowest energy eigenstates, we use the eigenfunction to calculate $S_{\text{tot}}$ and its IREP 444In principle, one can switch the order of “diagonalization” and “point group symmetry analysis,” to reach larger $M$ and get better convergence. Similar tricks can be found in Refs. Sinitsyn et al. (2007); Läuchli (2011)..

In contrast to the result obtained with perturbative continuous unitary transformations (CUTs) Knetter et al. (2000b), we find that the first instability as a function of $J^{\prime}/J$ (for $h=0$) takes place in the $S_{\text{tot}}=0$ channel with IREP $A_{2}$ 555The $S=1$ instability found in Ref. Knetter et al. (2000b) is checked by our method both in $S_{\text{tot}}^{z}=0$ and $S_{\text{tot}}^{z}=1$ sectors. The energy difference between the two sectors is negligible.. Figure 2 shows the evolution of the gap as a function of $M$. Convergence is reached beyond $M=3$ for $J^{\prime}/J\lesssim 0.5$, but the increase of $\xi$ slows down the convergence for larger $J^{\prime}/J$. Although Eq. (1) does not conserve the triplon number, the state that condenses is adiabatically connected with the corresponding $S_{\text{tot}}=0$ IREP $A_{2}$ pure two-triplon state in the $J^{\prime}/J\to 0$ limit (see Fig. 2).

We can read out the critical value of $J^{\prime}/J$ when this state condenses. Figure 3(a) shows the evolution of the critical value as we increase $M$ (circles). At $M=8$, $(J_{c}^{\prime}/J)^{(M=8)}\approx 0.702$. Previous tensor network based calculations Lou et al. ; Corboz and Mila (2013) showed that the transition is actually of first order and the transition point is at $J_{c}^{\prime}/J=0.675$ Corboz and Mila (2013). A susceptibility analysis, like the one presented here, is in general inadequate to detect first order transitions. However, it can still be used to detect the nature of the order parameter if the system still transitions into the broken symmetry state with highest susceptibility 666A typical example is the case of attractive interaction between the modes (particles) that become soft. Because of the attraction (negative quartic term in a Ginzburg-Landau expansion) the particle density changes discontinuously (first order transition) before the single-particle excitation becomes gapless. . Given that the first-order transition takes place when this susceptibility is still finite, $J_{c}^{\prime}/J$ turns out to be smaller than the value at which the susceptibility becomes divergent. This observation explains the difference between the values of $J_{c}^{\prime}/J$ obtained with both approaches. In addition, it illustrates their complementary nature. The unbiased susceptibility analysis can be used to detect candidates for broken symmetry states. These candidates can then be tested with biased variational techniques, such as iPEPS, which can produce more accurate values of the transition point.

Since the two-triplon bound state has $S_{\text{tot}}=0$, the new ground state (condensate) must be nonmagnetic. Furthermore, since the $A_{2}$ state is odd (even) under reflection (rotation) sup , the new ordered state should only break reflection but not rotation symmetries. These features are consistent with the previously reported plaquette ordering Koga and Kawakami (2000); Läuchli et al. (2002); Lou et al. ; Corboz and Mila (2013). Figure 4(b) shows a schematic plot of the corresponding bond ordering. As expected, $\langle\bm{S}_{i}\cdot\bm{S}_{j}\rangle$ becomes different on different plaquettes and there is no magnetic order. In other words, the plaquette order parameter can be defined as $\langle\bm{S}_{i}\cdot\bm{S}_{j}-\bm{S}_{i}\cdot\bm{S}_{j^{\prime}}\rangle$, where $ij$ and $ij^{\prime}$ are two bonds related by a mirror reflection [see Fig. 4(b)].

We consider now the case of nonzero magnetic field. The energy of excited states with finite $S_{\text{tot}}^{z}$ decreases linearly in $h$. Figure 3(b) shows that the dominant instability for $h/J=0.2$ corresponds to condensation of a state with $S_{\text{tot}}=1$, $\bm{k}=(0,0)$, and IREP $E$, leading to the stripe magnetic order depicted in Fig. 4(c). The IREP E is a two-dimensional representation corresponding to the two possible directions of the stripes (along $\bm{a}_{1}$ or $\bm{a}_{2}$). We note that the two same-color spins in the same unit cell are not identical (the corresponding mirror symmetry is broken).

The stripe state has the highest susceptibility over a narrow range $0.66\lesssim J^{\prime}/J\lesssim 0.70$ for $M=8$ iterations [see Fig. 4(a)]. Because of the frustrated exchange interactions, the energies of a few other states are not much higher than the stripe state sup . Among them, the lowest one is a two-triplon state with $\bm{k}=(0,0)$ and IREP $B_{1}$, corresponding to vector chiral order Chubukov and Starykh (2013); Parker and Balents (2017); Wang et al. (2017). Although their energies are slightly higher than the stripe magnetic instability within $M\leq 8$, the situation may change in the $M\rightarrow\infty$ limit, or if small perturbations are added to the original model.

The $S_{\text{tot}}=2$ excited states take over for higher magnetic field values. Figure 3(c) shows that for $h/J=0.6$ the lowest excited state is the $S_{\text{tot}}=2$ two-triplon bound state with momentum $\bm{k}=(0,0)$ and IREP $A_{2}$. The fact that this state and the $S_{\text{tot}}=0$ state that condenses at zero field belong to the same point group IREP $A_{2}$ indicates that the condensation of the $S_{\text{tot}}=2$ $A_{2}$ state also leads to “plaquette” ordering [shown in Fig. 4(d)]; the difference being that the $S_{\text{tot}}=2$ condensate also breaks the U(1) symmetry group of global spin rotations along the field direction, leading to spin-nematic ordering. In other words, the local bond order parameter is $\langle S_{i}^{+}S_{j}^{+}-S_{i}^{+}S_{j^{\prime}}^{+}\rangle$ instead of $\langle\bm{S}_{i}\cdot\bm{S}_{j}-\bm{S}_{i}\cdot\bm{S}_{j^{\prime}}\rangle$ ($ij$ and $ij^{\prime}$ denote two bonds connected by a mirror reflection $\sigma_{d}$, see Fig. 4).

As indicated in Fig. 4(a), the “plaquette spin-nematic” state covers a wide range $0.40\lesssim J^{\prime}/J\lesssim 0.66$. It has been shown in Ref. Corboz and Mila (2014) that the $\frac{1}{8}$ plateau at slightly higher magnetic field values and $J^{\prime}/J=0.63$ is induced by crystallization of $S_{\text{tot}}^{z}=2$ bound states. A closer scrutiny of the “pinwheel” structure of these bound states shows that they locally preserve rotational symmetries, while breaking reflection symmetries Corboz and Mila (2014); i.e., they are the same two-triplon bound states that we are finding in the dilute limit.

Moving away from the plaquette spin-nematic phase towards the $J^{\prime}/J\ll 1$ limit, it is already known from an early perturbative calculation that $S_{\text{tot}}=2$ two-triplon bound states with $\bm{k}=(\pi,\pi)$ give the highest susceptibility Momoi and Totsuka (2000b). This is confirmed by our variational method [see Fig. 3(d)]. Since the two-triplon bound state has momentum $\bm{k}=(\pi,\pi)$, the corresponding ordered state also breaks translational symmetry. As shown in Fig. 4(e), $\langle S_{i}^{+}S_{j}^{+}\rangle$ changes sign going from one unit cell to its nearest neighbors. Similar to the case of the stripe ordering (which also comes from condensation of IREP E states), there are two choices for aligning the bonds. We note that the breaking of the $C_{4}$ lattice rotational symmetry leads to a modulation of $\langle S_{i}^{z}\rangle$ that can be detected with NS experiments.

It is worth mentioning that the two spin-nematic states found in this Letter are different from nematic phases discussed in various other contexts Fradkin et al. (2010); Kamiya et al. (2011); Fernandes et al. (2012, 2014). The so-called “Ising-nematic” ordering corresponds to (discrete) lattice rotation symmetry breaking. In contrast, “spin-nematic” ordering corresponds to broken spin rotational symmetry. The spin-nematic orderings discussed in this Letter break both the point group symmetry and spin rotation symmetry Penc and Läuchli (2011). In other words, they are simultaneously Ising nematic and spin nematic.

The frustrated nature of the SSM makes the calculation of dynamical response a difficult task. To date, the only calculation including multitriplon contributions is the perturbative CUTs, which breaks down for $J^{\prime}/J\gtrsim 0.63$ Knetter and Uhrig (2004). The variational Hilbert space generated by our method thus provides a more reliable access to dynamical responses via the continued fraction method Gagliano and Balseiro (1987).

Near the phase boundaries, we expect the susceptibilities of corresponding order parameters to diverge at $\omega=0$. Magnetic orderings, such as the stripe phase, are detected by computing the dynamic structure factor (DSF) flu ; sup :

[TABLE]

which is measured with inelastic NS. As shown Fig. 5(b), the lowest peak of $S^{-+}(\bm{k},\omega)$ approaches $\omega=0$ near the phase boundary indicating condensation of an $S_{\text{tot}}=1$ state.

The divergent susceptibilities of the other phases are revealed by computing two-point dynamical correlation functions of the corresponding order parameters. These order parameters are the operators that create a state that has finite overlap with the one that is condensing. For $J^{\prime}/J\ll 1$, the lowest energy $S_{\text{tot}}=2$ eigenstates are known to be a linear combination of triplons located on nearest (and next-nearest) neighbors Momoi and Totsuka (2000b); sup . Denoting the order parameter as $A_{\bm{k}}^{S2E}$, the corresponding susceptibility is

[TABLE]

Similarly, using the approximate wave functions of the $S_{\text{tot}}=2$ $A_{2}$ state and the $S_{\text{tot}}=0$ $A_{2}$ state sup , we can also construct the order parameters and compute the corresponding susceptibilities $\chi_{S2A_{2}}(\bm{k},\omega)$ and $\chi_{S0A_{2}}(\bm{k},\omega)$. Figure 5 shows the nearly divergent susceptibilities in each channel by picking appropriate Hamiltonian parameters near the phase boundaries.

While the tendency toward stripe ordering can be detected with inelastic NS, the experimental detection of the other phases is nontrivial. Lattice distortions induced by the order parameter through magnetostriction can provide indirect evidence if they are large enough to be detected Chubukov and Starykh (2013); Parker and Balents (2017); Wang et al. (2017). Experimental knobs, such as pressure, doping, and magnetic field can drive the material into different instabilities Haravifard et al. (2016); Zayed et al. (2017). Thus, the method presented in this Letter provides valuable insight for revealing the nature of the new phases in such experiments. However, the model relevant to SrCu2(BO3)2 also includes DM interactions that modify the single triplon dispersion and can potentially change the phase diagram reported here. In addition, DM interactions reduce the spin rotational symmetry of the model, implying that they can change the nature of the order parameters.

We emphasize that the applicability of the method is not restricted to the SSM considered here. The same method can be used to detect the instabilities of other quantum paramagnets Zapf et al. (2014). Especially, it is very difficult to enumerate all the possible instabilities for highly frustrated systems. The low energy spectrum produced by our method provides a valuable educated guess for biased numerical approaches. Given that the method works in the thermodynamic limit, it can also detect incommensurate instabilities, which cannot be handled by most numerical methods.

In summary, we have used an unbiased variational method to study the excitation spectrum of the SSM in the dimer phase. Several instabilities are found next to the dimer phase corresponding to condensations of single-triplon or two-triplon bound states. Two of the instabilities (antiferro spin nematic and plaquette) are known from previous studies and the others (plaquette spin nematic and stripe) are newly discovered in this Letter. The same method can be used to compute relevant dynamical correlation functions.

We thank S. Haravifard, B. Shastry, G. Ortiz, S. Zhang, and H. Suwa for helpful discussions. Z.W. and C.D.B. are supported by funding from the Lincoln Chair of Excellence in Physics. This work used the Extreme Science and Engineering Discovery Environment (XSEDE) Towns et al. (2014) through allocation TG-DMR170029, which is supported by NSF Grant No. ACI-1548562. This research used resources of the Compute and Data Environment for Science (CADES) at the Oak Ridge National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC05-00OR22725.

Appendix A Character table of $C_{4v}$ group

The point group of Shastry-Sutherland lattice is $C_{4v}$, with 5 classes of symmetry operations: identity $E$, two rotations $\{C_{4},C_{2}\}$ and two reflections $\{\sigma_{v},\sigma_{d}\}$ (see Table. 1).

Appendix B Choice of initial basis

The initial basis is always chosen such that it does not break the point group symmetries explicitly. The choices for different values of $S_{\text{tot}}^{z}$ are listed in Figs. S1, S2 and S3.

Appendix C Construction of order parameters

In the $J^{\prime}/J\ll 1$ limit, the $S_{\text{tot}}=2$, $\bm{k}=(\pi,\pi)$, IREP $E$ states are known from Ref. Momoi and Totsuka (2000b). Besides the $S_{\text{tot}}^{z}=2$ initial states shown in Fig. S3, there is another set of relevant states (see Fig. S4) generated at $M=2$. Here we write down the diagonalized wavefunctions for $\bm{Q}\equiv(\pi,\pi)$:

[TABLE]

where

[TABLE]

The above wavefunctions have $S_{\text{tot}}=2$, as long as the normalization condition $2|\alpha|^{2}+2|\beta|^{2}=1$ is satisfied. The values of $\alpha,\beta$ depend on $J^{\prime}/J$ Momoi and Totsuka (2000b). To construct the nematic order parameter, we only need a bound state wavefunction that has significant overlap with the exact one. We then choose $\alpha=0.7$, $\beta=0.1$ to construct the order parameter

[TABLE]

The momentum dependence of the corresponding susceptibility is obtained by computing the two-point correlator after replacing $\bm{Q}$ with $\bm{k}$.

We have demonstrated in the main text that the lowest bound state wave function has $S_{\text{tot}}=2$, $\bm{k}=(0,0)$, IREP $A_{2}$ for larger values of $J^{\prime}/J$. Similarly, with the eight states in Figs. S3-S4, we can write down the approximate $\bm{K}\equiv(0,0)$ wavefunction :

[TABLE]

Once again, $S_{\text{tot}}=2$ is guaranteed by the normalization condition $4|\alpha|^{2}+4|\beta|^{2}=1$. To construct the corresponding order parameter $A^{S2A_{2}}_{\bm{k}}$, we use $\alpha=0.47$, $\beta=0.17$ in the main text.

We note that the exact bound state in this region (larger $J^{\prime}/J$), has a larger size than the one in the $J^{\prime}/J\ll 1$ limit. Consequently, $A^{S2A_{2}}_{\bm{k}}$ produces more spectral weight at high energies in comparison to $A^{S2E}_{\bm{k}}$ [see Fig. 5(c)].

For even larger $J^{\prime}/J$ with $h=0$, the lowest states have $S_{\text{tot}}=0$, $\bm{k}=(0,0)$, IREP $A_{2}$. we can similarly construct the approximate wavefunction at $\bm{K}\equiv(0,0)$:

[TABLE]

In this case, the normalization condition $8|\alpha|^{2}+4|\beta|^{2}=1$ is different from the $S_{\text{tot}}=0$ condition, which requires $\alpha=\beta$. Combining both conditions we have $\alpha=\beta=\frac{1}{2\sqrt{3}}$.

Once again, since this approximate wavefunction has a reduced overlap with the exact bound state, we see a significant amount of spectral weight at high energies in Fig. 5(a) in the main text.

Appendix D Dominant instability at zero field

In the $S_{\text{tot}}=0$ sector at $h=0$, we can roughly estimate $J_{c}^{\prime}/J\approx 0.69$ in the $M\rightarrow\infty$ limit (see Fig. 3 in the main text). At this point $J^{\prime}/J=0.69$, we show In Fig. S5 that the gaps of the other $S_{\text{tot}}$ sectors remain finite in the $M\rightarrow\infty$ limit, even with a linear extrapolation. Note that the calculation converges for large enough $M$, implying that the slopes of the curves in Fig. S5 should finally become zero. In other words, the true value of the gaps should be larger than that obtained from linear extrapolation. Since all other gaps remain finite while the $S_{\text{tot}}=0$ gap is closing, we can safely conclude that the order parameter with highest susceptibility has $S_{\text{tot}}=0$ for $h=0$.

Appendix E Convergence of dynamic structure factors

Generally, the convergence speed of the DSFs is determined by $\xi$. Larger $\xi$ values require larger number of iterations ($M$) to obtain qualitatively correct results. For instance, $\chi_{S2E}(\bm{k},\omega)$ is almost fully converged at $M=3$ and $J^{\prime}/J=0.2$ [see Fig. 5(d) in the main text] because $\xi$ is a few lattice spaces.

To have a more quantative sense of the convergence as a function of $M$, we first benchmark the DSF $S^{-+}(\bm{k},\omega)$ with $J^{\prime}/J=0.63$, $h=0$, which is relevant for the compound SrCu2(BO3)2 Corboz and Mila (2014). In Fig. S6, we show the evolution of the calculated $S^{-+}(\bm{k},\omega)$ as a function of $M$. The lowest single and two-triplon bands are already converged at $M=5$. We note that the method produces not only the bound states, but also the multitriplon continuum, which is qualitatively caputured at $M=5$ and about to converge at $M\gtrsim 8$ for the parameters used in Fig. S6.

To measure the relative change between different iterations, we can define

[TABLE]

where $\Lambda$ is the energy cutoff. For the fully converged result, it is expected that $\Delta_{M}(\bm{k})\rightarrow 0$.

In Fig. S7, we show the linecuts of Fig. S6 at two momenta: $\bm{k}-(4\pi,0)=(0,0)\text{ and }(\pi,\pi)$. Correspondingly, the relative errors $\Delta_{M}(\bm{k})$ are shown in the bottom panels of the same figure. The relative errors decrease rapidly as we increase $M$, reaching around $15\%$ at $M=8$. We note that the value of $\Delta_{M}(\bm{k})$ depends on the energy cutoff $\Lambda$: the lowest relative errors in Fig. S7 are found when $\Lambda/J=1$, consistent with the fact that low-energy excitations converge faster than the high-energy ones. Note also that $\Delta_{M}(\bm{k})$ depends heavily on the choice of Lorentzian broadening factor used in the DSF calculation: larger broadening factors lead to smaller relative errors.

It is also worth benchmarking with Ref. Knetter and Uhrig (2004), which carried out the DSF calculation for SSM using the perturbative CUTs. Figure S8 shows a good qualitative agreement. However, there are two small differences: the band near $5$meV produced by our method is slightly higher than that obtained from CUT, which could arise from incomplete convergence for $M\lesssim 8$ in our method or missing higher order corrections in CUT. We also note that our method captures more accurately the states in the continuum.

Appendix F Competing chiral state

The main text describes the $S_{\text{tot}}=1$ instability that leads to the stripe magnetic order. As shown in Fig. 5(b), the lowest instability from the DSF arises from condensation of the $\bm{k}=(0,0)$ mode. Figure S9 examines this point in more detail. Besides reproducing Fig. 5(b) in higher resolution, we also denote the symmetries of the $\bm{k}=(4\pi,0)$ peaks [$\bm{k}=(0,0)\equiv(4\pi,0)$].

As we mentioned in the main text, the state that produces the lowest peak at $\bm{k}=(4\pi,0)$ belongs to the IREP $E$. In Fig. S9, we see that the next peak at $\bm{k}=(4\pi,0)$ corresponds to the IREP $B_{1}$ [with negligible weight at $\bm{k}=(4\pi,0)$]. Up to $M=8$, we find that this $B_{1}$ peak is always higher in energy (lower susceptibility) than the one corresponding to the IREP $E$. It is worth mentioning that condensation of the $B_{1}$ state leads to the vector chiral ordering depicted in Fig. S10.