Signatures of a Deconfined Phase Transition on the Shastry-Sutherland Lattice: Applications to Quantum Critical SrCu$_2$(BO$_3$)$_2$

TL;DR

This paper investigates a deconfined quantum phase transition in a Shastry-Sutherland lattice model, providing numerical and theoretical evidence for an emergent O(4) symmetry and predicting experimental signatures in SrCu$_2$(BO$_3$)$_2$.

Contribution

It demonstrates the existence of a deconfined quantum critical point with emergent O(4) symmetry in a realistic lattice model, supported by iDMRG analysis and correlation spectrum analysis.

Findings

Evidence for a deconfined quantum critical point (DQCP) between N{\'e}el and pVBS phases.

Identification of emergent O(4) symmetry at the transition.

Predicted sharp signatures of DQCP in phonon and magnon spectra.

Abstract

We study a possible deconfined quantum phase transition in a realistic model of a two-dimensional Shastry-Sutherland quantum magnet, using both numerical and field theoretic techniques. Using the infinite density matrix renormalization group (iDMRG) method, we verify the existence of an intermediate plaquette valence bond solid (pVBS) order, with two fold degeneracy, between the dimer and N\'eel ordered phases. We argue that the quantum phase transition between the N\'eel and pVBS orders may be described by a deconfined quantum critical point (DQCP) with an emergent O(4) symmetry. By analyzing the correlation length spectrum obtained from iDMRG, we provide evidence for the DQCP and emergent O(4) symmetry in the lattice model. Such a phase transition has been reported in the recent pressure tuned experiments in the Shastry-Sutherland lattice material $\mathrm{SrCu_2 (BO_3)_2}$. The non-symmorphic lattice structure of the Shastry-Sutherland compound leads to extinction points in the scattering, where we predict sharp signatures of a DQCP in both the phonon and magnon spectra associated with the spinon continuum. The effect of weak interlayer couplings present in the three dimensional material is also discussed. Our results should help guide the experimental study of DQCP in quantum magnets.

Figures

Tables

Table 1. Table 1: ( J 1 / J 2 ) c ​ 1 subscript subscript 𝐽 1 subscript 𝐽 2 𝑐 1 (J_{1}/J_{2})_{c1} ( ( J 1 / J 2 ) c ​ 2 subscript subscript 𝐽 1 subscript 𝐽 2 𝑐 2 (J_{1}/J_{2})_{c2} ) is the transition point between the dVBS and pVBS (pVBS and Néel) phases. Transition points are extracted from the peak of the energy derivative at χ = 4000 𝜒 4000 \chi=4000 . At L = 12 𝐿 12 L=12 , the DMRG simulation do not converge for different initial states at χ = 4000 𝜒 4000 \chi=4000 , resulting in different transition points. Thus, the critical point is determined as the midpoint between two transition points obtained from the VBS-like and Néel-like initial states.

	$L=6$	$L=8$	$L=10$	$L=12$
${\left(J_1/J_2\right)}_{c1}$	0.682	0.677	0.675	0.675
${\left(J_1/J_2\right)}_{c2}$	0.693	0.728	0.762	0.77

Table 2. Table 2: Symmetry transformation of momentum ( k x , k y ) subscript 𝑘 𝑥 subscript 𝑘 𝑦 (k_{x},k_{y}) , Néel order 𝒏 𝒏 {\bm{n}} , monopole operator ℳ † superscript ℳ † \mathcal{M}^{\dagger} , spin operator 𝑺 i subscript 𝑺 𝑖 {\bm{S}}_{i} and fermionic spinon f i subscript 𝑓 𝑖 f_{i} (in the sense of PSG).

	$(k_x,k_y)$	$𝒏$	${ℳ}^{†}$	${𝑺}_i$	$f_i$
$G_x$	$(k_x,-k_y)$	$-𝒏$	$-{ℳ}^{†}$	${𝑺}_{G_x\left(i\right)}$	${(-)}^y{f}_{G_x\left(i\right)}$
$G_y$	$(-k_x,k_y)$	$-𝒏$	$-{ℳ}^{†}$	${𝑺}_{G_y\left(i\right)}$	$f_{G_y\left(i\right)}$
${\sigma }_{xy}$	$(k_y,k_x)$	$𝒏$	$ℳ$	${𝑺}_{{\sigma }_{xy}\left(i\right)}$	${(-)}^{xy}{f}_{{\sigma }_{xy}\left(i\right)}$
${\sigma }_{x\overline{y}}$	$(-k_y,-k_x)$	$𝒏$	$ℳ$	${𝑺}_{\sigma x\overline{y}\left(i\right)}$	${(-)}^{x\left(y+1\right)}{f}_{{\sigma }_{x\overline{y}}\left(i\right)}$
$𝒯$	$(-k_x,-k_y)$	$-𝒏$	$ℳ$	$-{𝑺}_i$	$𝒦i{\sigma }^2{f}_i$
$𝒯G_x$	$(-k_x,k_y)$	$𝒏$	$-ℳ$	$-{𝑺}_{G_x\left(i\right)}$	${(-)}^y𝒦i{\sigma }^2{f}_{G_x\left(i\right)}$
$𝒯G_y$	$(k_x,-k_y)$	$𝒏$	$-ℳ$	$-{𝑺}_{G_y\left(i\right)}$	$𝒦i{\sigma }^2{f}_{G_y\left(i\right)}$

Table 3. Table 3: Singlet configuration energy around different plaquette with different pairing symmetry, estimated from J 1 / J 2 ∼ 0.7 similar-to subscript 𝐽 1 subscript 𝐽 2 0.7 J_{1}/J_{2}\sim 0.7 for the pVBS phase.

	$s$-wave	$d$-wave
square	$-2J_1\sim -1.4J_2$	0
diamond	$-2J_1+0.25J_2\sim -1.15J_2$	$-0.75J_2$

Table 4. Table 4: Symmetry-allowed most-relevant monopole terms (apart from r ​ ℳ † ​ ℳ 𝑟 superscript ℳ † ℳ r\mathcal{M}^{\dagger}\mathcal{M} ) and the corresponding DQCP emergent symmetries on different lattices

Global Symmetry	${ℒ}_{ℳ}$	Emergent Symmetry
Square ($p4m$)	${\lambda }_{4}Re{ℳ}^4$	$\mathrm{SO}\left(3\right)\times {\mathbb{Z}}_4\to \mathrm{SO}\left(5\right)$
Easy-plane Square	${\lambda }_{4}Re{ℳ}^4$	$U\left(1\right)\times {\mathbb{Z}}_4\to O\left(4\right)$
Rectangular ($pmm$)	${\lambda }_2^{\prime }Im{ℳ}^2$	$\mathrm{SO}\left(3\right)\times {\mathbb{Z}}_2\to O\left(4\right)$
Shastry-Surtherland ($p4g$)	${\lambda }_{2}Re{ℳ}^2$	$\mathrm{SO}\left(3\right)\times {\mathbb{Z}}_2\to O\left(4\right)$

Table 5. Table 5: Transformation of the monopole operator ℳ ℳ {\cal M} and Néel vector 𝒏 𝒏 {\bm{n}} in the field theory of nonlinear sigma model in the Néel order phase.

Symmetries	Transformations
$T_x$	${ℳ}^{†}\mapsto -iℳ$  $𝒏\mapsto -𝒏$
$T_y$	${ℳ}^{†}\mapsto iℳ$  $𝒏\mapsto -𝒏$
$R_{\pi /4}^{\text{site}}$	${ℳ}^{†}\mapsto i{ℳ}^{†}$  $𝒏\mapsto 𝒏$
$R_{\pi /4}^{\text{plaq}}$	${ℳ}^{†}\mapsto -ℳ$  $𝒏\mapsto -𝒏$
${\sigma }_x$	${ℳ}^{†}\mapsto iℳ$  $𝒏\mapsto 𝒏$
${\sigma }_y$	${ℳ}^{†}\mapsto -iℳ$  $𝒏\mapsto 𝒏$
$𝒯$	${ℳ}^{†}\mapsto ℳ$  $𝒏\mapsto -𝒏$

Table 6. Table 6: The table shows the correspondence between the symmetries of the Shastry-Sutherland lattice ( p ​ 4 ​ g 𝑝 4 𝑔 p4g ) and square lattice ( p ​ 4 ​ m 𝑝 4 𝑚 p4m ). Bold symbols are for symmetries of the Shastry-Sutherland lattice.

$G_{p4g}$	$G_{p4m}$	Action
${𝑻}_x$	$T_x^2$	${ℳ}^{†}\mapsto {ℳ}^{†}$  $𝒏\mapsto 𝒏$
${𝑻}_y$	$T_y^2$	${ℳ}^{†}\mapsto {ℳ}^{†}$  $𝒏\mapsto 𝒏$
${𝝈}_{xy}$	$R_{\pi /2}{\sigma }_x$	${ℳ}^{†}\mapsto ℳ$  $𝒏\mapsto 𝒏$
${𝝈}_{x\overline{y}}$	$T_x{T}_y{R}_{\pi /2}{\sigma }_y$	${ℳ}^{†}\mapsto ℳ$  $𝒏\mapsto 𝒏$
${𝒈}_x$	$T_x{\sigma }_x$	${ℳ}^{†}\mapsto -{ℳ}^{†}$  $𝒏\mapsto -𝒏$
${𝒈}_y$	$T_y{\sigma }_y$	${ℳ}^{†}\mapsto -{ℳ}^{†}$  $𝒏\mapsto -𝒏$
${𝑹}_{\pi /4}$	$R_{\pi /4}^{\text{plaq}}$	${ℳ}^{†}\mapsto -ℳ$  $𝒏\mapsto -𝒏$
$𝒯$	$𝒯$	${ℳ}^{†}\mapsto ℳ$  $𝒏\mapsto -𝒏$

Table 7. Table 7: Momentum and symmetry transformations of the VBS order parameter 𝒗 𝒗 {\bm{v}} and finite-momentum phonon mode 𝒖 𝒖 {\bm{u}} .

	$v_x$	$v_y$	$u_x$	$u_y$
$𝒒$	$(\pi ,0)$	$(0,\pi )$	$(\pi ,0)$	$(0,\pi )$
$G_x$	$-v_x$	$-v_y$	$-u_x$	$-u_y$
$G_y$	$-v_x$	$-v_y$	$-u_x$	$-u_y$
${\sigma }_{xy}$	$v_y$	$v_x$	$u_y$	$u_x$
${\sigma }_{x\overline{y}}$	$v_y$	$v_x$	$u_y$	$u_x$

Equations

$$H=J_1\sum _{ij\in \text{n.n.}}{\mathbf{S}}_i\cdot {\mathbf{S}}_j+J_2\sum _{ij\in \text{dimer}}{\mathbf{S}}_i\cdot {\mathbf{S}}_j,$$

$$\text{imaginary}⟨\mathcal{M}⟩\sim \sum _i(-{)}^x⟨{\mathbf{S}}_i\cdot {\mathbf{S}}_{i+\hat{x}}⟩-(-{)}^y⟨{\mathbf{S}}_i\cdot {\mathbf{S}}_{i+\hat{y}}⟩.$$

$$\text{expectationvalue}\mathbf{n}(r)\cdot \mathbf{n}(0)\sim \text{expectationvalue}\mathbf{v}(r)\cdot \mathbf{v}(0)\sim 1/r^{1+\eta },$$

$${\mathcal{L}}_{\text{CP1}}=\text{absolutevalue}{(\partial -ia)z}^2+\kappa \text{quantity}(\nabla \times a{)}^2+\dots ,$$

$$\begin{array}{rl}& \mathbf{n}\sim (-{)}^{x+y}{\mathbf{S}}_i,\\ & v_x\sim (-{)}^x(\frac{1}{4}-{\mathbf{S}}_i\cdot {\mathbf{S}}_{i+\hat{x}}),\\ & v_y\sim (-{)}^y(\frac{1}{4}-{\mathbf{S}}_i\cdot {\mathbf{S}}_{i+\hat{y}}),\end{array}$$

$$\mathcal{M}^{\dagger}=\frac{1}{\sqrt{2}}\big{(}(v_{x}+v_{y})+\mathrm{i}(v_{x}-v_{y})\big{)},$$

$${\mathcal{L}}_{\mathcal{M}}=r{\mathcal{M}}^{†}\mathcal{M}+{\lambda }_2\Re {\mathcal{M}}^2+\cdots .$$

$${\mathcal{L}}_{\mathcal{M}}=\tilde{r}(\text{imaginary}\mathcal{M}{)}^2+{\tilde{\lambda }}_2(\Re \mathcal{M}{)}^2+\cdots ,$$

$$S(Q_x,Q_y)=e^{\mathrm{i}{Q}_x}S(Q_x,-Q_y)=e^{\mathrm{i}{Q}_y}S(-Q_x,Q_y),$$

$$\begin{array}{rl}{H}_{\text{MF}}=& -\sum _{⟨ij⟩}{t}_{ij}(f_i^{†}{f}_j+\text{h.c.})\\ & +g\sum _{[ijkl]}(f_i^{†}{f}_j+\text{h.c.})(f_k^{†}{f}_k-f_l^{†}{f}_l)+\cdots ,\end{array}$$

$${\omega }_{\text{min}}(\mathbf{q})=2t\sqrt{{\sin }^2{q}_x+{\sin }^2{q}_y},$$

$$S_{\text{magnon}}(\omega ,\mathbf{q}=X)\propto \omega .$$

$${\mathcal{L}}_{\text{VBS-phonon}}\sim v_x{u}_x+v_y{u}_y,$$

$$S_{\text{phonon}}(\omega ,\mathbf{q}=X)\propto {\omega }^{-2+\eta .}$$

$${\mathcal{S}}_{\text{DQCP}}^{(l)}=\int d^{3}x\frac{1}{2g}({\partial }_{\mu }{\Phi }^{(l)}{)}^2+2\pi \mathrm{i}{\Gamma }_{\text{WZW}}[{\Phi }^{(l)}]$$

$$\mathcal{S}=\sum _l{\mathcal{S}}_{\text{DQCP}}^{(l)}-\sum _l\int d^{3}x{g}^{ab}{\Phi }_a^{(l)}{\Phi }_b^{(l+1)},$$

$$\mathcal{L}=\bar{\psi }({\gamma }^{\mu }{D}_{\mu }+m+\mathrm{i}{\gamma }^5{\Phi }_a{\Gamma }^a)\psi .$$

$$S=\frac{1}{2g}\int d\tau d^{2}r\text{quantity}[\text{quantity}(\frac{\partial \mathbf{n}}{\partial \tau }{)}^2+c^2({\nabla }_r\mathbf{n}{)}^2]+S_B$$

$$g:{\mathcal{M}}_r^{†}\mapsto {\alpha }_g\cdot g[{\mathcal{M}}^{†}{]}_{g(r)},g[{\mathcal{M}}^{†}]={\mathcal{M}}^{†}\text{ or }\mathcal{M},$$

$$\frac{\eta (R_{\pi /2}^{\text{site}}r)\beta }{\eta (r)\beta }=i⟹{\mathcal{M}}^{†}\mapsto i{\mathcal{M}}^{†}.$$

$$\frac{{\eta }^{\ast }(r+\hat{x}){\beta }^{\ast }}{\eta (r)\beta }=-i⟹{\mathcal{M}}^{†}\mapsto -i\mathcal{M}.$$

$$H=J_1\sum _{\text{expectationvalue}i,j\in \text{blue}}{\mathbf{S}}_i\cdot {\mathbf{S}}_j+J_1^{\prime }\sum _{\text{expectationvalue}i,j\in \text{red}}{\mathbf{S}}_i\cdot {\mathbf{S}}_j$$

$$H=J_1\sum _{\text{expectationvalue}i,j}{\mathbf{S}}_i\cdot {\mathbf{S}}_j+J_2\sum _{\text{expectationvalue}\text{expectationvalue}i,j}{\mathbf{S}}_i\cdot {\mathbf{S}}_j$$

$$\begin{array}{rl}{\mathbf{r}}_A& =(1+\delta ,1+\delta )/2,\\ {\mathbf{r}}_B& =(1-\delta ,3+\delta )/2,\\ {\mathbf{r}}_C& =(3+\delta ,1-\delta )/2,\\ {\mathbf{r}}_D& =(3-\delta ,3-\delta )/2,\end{array}$$

$$\begin{array}{rl}\overline{AD}& =\frac{1+\delta }{2}\approx \frac{2.905\mathring{\mathrm{A}}}{8.995\mathring{\mathrm{A}}},\\ \overline{AB}& =\frac{\sqrt{1+{\delta }^2}}{2}\approx \frac{5.132\mathring{\mathrm{A}}}{8.995\mathring{\mathrm{A}}}.\end{array}$$

$${\mathbf{r}}_i\to {\mathbf{r}}_i^{\prime }={\mathbf{r}}_i+{\mathbf{u}}_i(i=A,B,C,D).$$

$$\begin{split}E_{\text{bond}}[{\bm{u}}_{i}]=&\frac{k_{1}}{2}\sum_{ij\in\text{n.n.}}\big{(}({\bm{r}}^{\prime}_{i}-{\bm{r}}^{\prime}_{j})^{2}-({\bm{r}}_{i}-{\bm{r}}_{j})^{2}\big{)}^{2}\\ +&\frac{k_{2}}{2}\sum_{ij\in\text{dimer}}\big{(}({\bm{r}}^{\prime}_{i}-{\bm{r}}^{\prime}_{j})^{2}-({\bm{r}}_{i}-{\bm{r}}_{j})^{2}\big{)}^{2},\end{split}$$

$$\begin{array}{rl}E[{\mathbf{u}}_i]& =E_{\text{bond}}[{\mathbf{u}}_i]+E_{\text{VBS}}[{\mathbf{u}}_i],\\ E_{\text{VBS}}[{\mathbf{u}}_i]& =\text{imaginary}\mathcal{M}\sum _p(-{)}^p\sum _{i\in p}({\mathbf{r}}_i^{\prime }-{\mathbf{R}}_p{)}^2,\end{array}$$

$${\bm{R}}_{p}=\left\{\begin{array}[]{ll}(1,0)&p\in X,\\ (0,1)&p\in Y.\end{array}\right.$$

$$\begin{array}{rl}{\mathbf{u}}_A& =\frac{\text{imaginary}\mathcal{M}}{4k_1}(-1,1),\\ {\mathbf{u}}_B& =\frac{\text{imaginary}\mathcal{M}}{4k_1}(-1,-1),\\ {\mathbf{u}}_C& =\frac{\text{imaginary}\mathcal{M}}{4k_1}(1,1),\\ {\mathbf{u}}_D& =\frac{\text{imaginary}\mathcal{M}}{4k_1}(1,-1).\end{array}$$

Full text

1. I Introduction

2. II Numerical Study

3. II.1 Detection of the pVBS Phase

4. II.2 Correlation Length Spectra, Monopole Fluctuations, and Emergent $O(4)$ Symmetry

5. III Symmetry Analysis

6. IV Dangerously Irrelevant Scaling and Its Absence

7. V Spectral Signatures of DQCP

8. VI Effects of Interlayer Coupling

9. VII Predictions for Experiment

10. VIII Conclusions

11. A Monopole Transformation

12. B Comparison with a DMRG simulation result of the $J_{1}$-$J_{2}$ model in the square lattice

13. C VBS-Phonon Coupling and Phonon Spectrum

14. D Detailed Numerical Data

I Introduction

Quantum magnets can host some of the most exotic phenomena in condensed matter physics, due to the strong quantum fluctuations of the microscopic spin degrees of freedom. Notably, deconfined gauge fluctuations and fractionalized spinon excitations can emerge in quantum magnets, which bear no analog in classical spin systems. Such behavior can exist either in a quantum spin liquid, which is a stable phase of matter with topological order Saminadayar et al. (1997); de Picciotto et al. (1997); Martin et al. (2004); or by tuning a single parameter to a critical point known as the deconfined quantum critical point (DQCP) Senthil et al. (2004a, b). The DQCP describes the possible continuous phase transition between two distinct symmetry breaking phases, which is beyond the conventional Landau-Ginzburg paradigm. While the search for quantum spin liquids is still an ongoing research effort in condensed matter physics Savary and Balents (2017), the possibility of observing the DQCP in materials could provide us with an alternative opportunity to study the properties of deconfined spinons and emergent gauge fields in quantum magnets, as well as in interacting fermion systems that realizes the DQCP Assaad and Grover (2016a); Sato et al. (2017); You et al. (2018); Zhang et al. (2018a); Ippoliti et al. (2018); Liu et al. (2019).

In a recent experiment Zayed et al. (2017), a phase transition between Néel antiferromagnet and plaquette valence bond solid (crystal) Mambrini et al. (2006) was observed in a single crystal of $\mathrm{SrCu_{2}(BO_{3})_{2}}$ under pressure. The material is a layered quantum magnet. Within each two-dimensional layer, the copper ions carry the spin-1/2 degrees of freedom and are arranged on a Shastry-Sutherland lattice as shown in Fig. 1(a). The spin system was proposed to be effectively described by the Shastry-Sutherland model Shastry and Sutherland (1981); Miyahara and Ueda (1999)

[TABLE]

where the $J_{1}$ and $J_{2}$ bonds are specified according to Fig. 1(a). The ratio $J_{1}/J_{2}$ between the coupling constants is tunable by pressure in experiments within certain range. In the large $J_{1}$ (or large $J_{2}$) limit, the model reduces to the square lattice Heisenberg model (or the decoupled dimerized model), which stabilizes the Néel phase (or the dimer valence bond solid (dVBS) phase). Between these two limits, numerical Koga and Kawakami (2000); Läuchli et al. (2002); Mambrini et al. (2006); Lou et al. (2012); Corboz and Mila (2013) and theoretical Chung et al. (2001) analysis of the model have revealed an intermediate plaquette valence bond solid (pVBS) phase, as illustrated in Fig. 1(c). Remarkably, the experiment in Ref.Zayed et al., 2017 seems to confirm this phase diagram. Since the pVBS and Néel phases separately break two distinct symmetries, the lattice and the spin rotation symmetry, a direct second-order transition between them would necessarily go beyond the Landau-Ginzburg paradigm and point to the possibility of the DQCP. Although the nature of the pVBS-Néel transition remains unresolved by experiments, there are promising signs for the exciting opportunity that $\mathrm{SrCu_{2}(BO_{3})_{2}}$ might provide the first experimental platform to realize DQCP.

Recent studies on different models with the same symmetry class showed that the transition between pVBS and Néel phases could be first-orderZhao et al. (2018); Serna and Nahum (2018). However, despite being first-order, the transition is accompanied with an extended region of quantum-critical-like scaling and an emergent $\mathrm{O}(4)$ symmetry, implying that the transition could be close to a DQCP (possibly as an avoided criticality). Thus the DQCP is still the best theory to account for these anomalous features in the critical region, even though it may eventually break down at longest scales. Note that the $J$-$Q$ model or loop model studied in Monte Carlo simulationsZhao et al. (2018); Serna and Nahum (2018) are designed differently from the original Shastry-Sutherland model to avoid the sign problem. Given that the first- or second-order nature of the transition can be tuned by model parameters Sen and Sandvik (2010); Qin et al. (2017); Zhang et al. (2018b); Jian et al. (2017a) and is therefore a model-dependent property, the fate of the pVBS-Néel transition in the Shastry-Sutherland model remains to be fully resolved yet.

The goal of this work is to investigate the pVBS-Néel transition in the Shastry-Sutherland model Eq. (1) in more detail using both field theory and the density matrix renormalization group (DMRG) approach, and to identify the unique signatures of DQCP that can be probed by inelastic neutron scattering (INS) or resonant inelastic X-ray scattering (RIXS) experiments. We use the infinite DMRG technique to overcome the sign problem. Our numerical simulation indicates (i) that the transition between pVBS and Néel phases appears continuous up to the largest available system size (infinite cylinder with the circumference of 10 lattice sites), although we can not rule out the possibility of a weakly first-order transition due to our limited system size. (ii) We also observe the asymptotic degeneracy between spin-triplet and spin-singlet excitations over a large length scale, demonstrating an approximate emergent $\mathrm{O}(4)$ symmetry which rotates among the Néel and pVBS order parameters. Our theoretical analysis further suggests that (iii) in the Shastry-Sutherland lattice, in contrast to previous realizations of DQCP, a dangerously irrelevant operator is absent which has consequences for numerics and that (iv) critical spinon continua appear at the extinction points of lattice diffraction peaks (c.f. Fig. 1(b)) in both the magnon and phonon channels at low-temperature around the DQCP. The universal critical behaviors of these continua are examined as well, which could guide the experimental study of the candidate DQCP in the $\mathrm{SrCu_{2}(BO_{3})_{2}}$ material.

The rest of the paper is organized as follows. In Sec. II, we perform an infinite DMRG simulation on the Shastry-Sutherland spin model and discuss the nature of the phase transition between Néel and VBS phases based on a correlation length spectra. In Sec. III, we analyze symmetry quantum numbers of a monopole operator whose proliferation induces the transition to the VBS phase. By investigating the transformation property of the monopole, we show that the single-monopole term is suppressed while the double-monopole term can appear in the action describing the Néel order in the Shastry-Sutherland lattice. We compare the differences among various microscopic models – easy-plane, rectangular, and Shastry-Sutherland – whose possible emergent symmetry is all $\mathrm{O}(4)$. One distinct feature of the Shastry-Sutherland lattice is the presence of the relevant anisotropy operator that breaks the four-fold lattice rotation symmetry, which stabilizes the VBS order and gives rise to a fast-growing spectral gap in the spin-0 channel as the system enters the VBS phase. In Sec. V, we propose the spectral signatures of DQCP in the Shastry-Sutherland model, including the $\mathrm{SO}(4)$ conserved current fluctuation in the magnon spectrum and the VBS fluctuation in the phonon spectrum. Both of these features appear at the extinction point of the Shastry-Sutherland lattice, which is detectable at low energy without overwhelmed by the elastic scattering signals. We conclude our discussion in Sec. VII.

II Numerical Study

In this section, we study the model in Eq. (1) using the infinite density matrix renormalization group (iDMRG) method White (1992); McCulloch (2008). We wrap the 2D system onto a cylinder which is infinite along $x$-direction but compact along $y$-direction, with a finite circumference $L$. In the simulation, spin-1/2s are mapped into the hard-core bosons with density $\expectationvalue{n}=1/2$ per site; an antiferromagnetic spin interaction would then be translated into hopping and density-density interaction terms for bosons. As the boson number is conserved in the simulation, we have an explicit $\mathrm{U}(1)_{z}$ symmetry which allows us to extract correlation functions for an operator with a specific $\mathrm{U}(1)_{z}$ quantum number. During the simulation, we fixed the value of $J_{2}=1$ and tuned the value of $J_{1}$ across the phase transition between the pVBS and Néel order phases. Since the unit cell of Shastry-Sutherland lattice contains $2\times 2$ square unit cells, the iDMRG unit cell becomes $2\times(2m)$ slice of the infinite cylinder, where the circumference of the cylinder $L=2m$.

In the iDMRG simulation, we have two limiting factors to describe the exact two-dimensional state: the circumference length $L$ and the bond dimension $\chi$. Due to limited computational capacity, it is difficult to conclusively identify whether the phase transition is continuous or weakly first order in the DMRG simulation. Still, we can extract useful information of the ground state by simulating the model with increasing values of $L$ and $\chi$. For the discussionon bond dimension scaling, see Appendix. D. In the DMRG simulation, we measured (i) energy, (ii) plaquette order parameter, and (iii) correlation length spectra.

II.1 Detection of the pVBS Phase

Although the matrix product state description of the state is exact in the infinite bond dimension limit, for a finite bond dimension, the cylindrical geometry of the iDMRG simulation provides some bias to the preferred entanglement structure for the ground state. As a result, in the pVBS phase, one of two symmetry broken phases would be automatically chosen and the order parameter would not vanish in the iDMRG simulation. To characterize the pVBS phase, we define the plaquette order parameter as follows (See Eq. (5) and Eq. (6) in Sec. III for a more rigors definition),

[TABLE]

In the pVBS phase, the dimer strengths $\langle{\bm{S}}_{i}\cdot{\bm{S}}_{j}\rangle$ on each bond $\langle ij\rangle$ can be visualized in Fig. 2, which clearly demonstrates the pattern of the pVBS ordering around empty (square) plaquettes. Indeed, we measured that $\imaginary\langle\mathcal{M}\rangle\neq 0$ and $\real\langle\mathcal{M}\rangle=0$ in the paramagnetic regime as expected for the pVBS phase. (See Fig. 4)

Although the plaquette order parameter $\imaginary\langle\mathcal{M}\rangle$ is a useful indicator, it is not precise because (i) the system size does not reach the thermodynamic limit and (ii) the geometry of the iDMRG simulation provides some bias toward a specific entanglement structure for the ground state. Albeit small, it has a non-vanishing value in the Néel phase, see Fig. 4. Thus, we use a discontinuity in the second derivative of the energy $\partial^{2}E/\partial J_{1}^{2}$ to locate the pVBS-Néel transition point. Up to $L=10$, the first order derivative of energy is continuous across the phase transition, implying that the transition is either a weakly first order or second order transition. On the other hand, from the energy plot in Fig. 3, the dVBS-pVBS transition point can be easily extracted because the dVBS state is an exact ground state of Eq. 1 with energy per site $E_{\textrm{site}}=-0.375J_{2}$. As we can see, the first order derivative of the energy is discontinuous here, signaling a clear first-order phase transition between the two spin singlet dVBS and pVBS phases. For $L<6$, we do not observe the pVBS phase. Transition points for different system sizes are summarized in Tab. 1.

II.2 Correlation Length Spectra, Monopole Fluctuations, and Emergent $O(4)$ Symmetry

Here, we present a signature of the DQCP in the correlation length spectrum data obtained from the iDMRG simulation. A continuous phase transition is characterized by the divergence of a correlation length, and an equal-time correlation function exhibits the power-law decaying behavior $\langle\mathcal{O}(r)\mathcal{O}(0)\rangle\propto r^{-2\Delta_{\cal O}}$ where $\Delta_{\cal O}$ is a scaling dimension of an operator ${\cal O}$. In particular, if there exists an emergent symmetry, operators unified under the emergent symmetry should share the same scaling dimension $\Delta_{\cal O}$ Nahum et al. (2015a). For example, at the conventional DQCP with the emergent $\mathrm{SO}(5)$ symmetry, the Néel ${\bm{n}}=(n_{x},n_{y},n_{z})$ and VBS ${\bm{v}}=(v_{x},v_{y})$ order parameter should have the same power-law behavior:

[TABLE]

where $\eta=2\Delta_{\cal O}-1$ is the anomalous exponent defined relative to the engineering exponent in (2+1)D which is $1$. However, in the iDMRG simulation, the finite circumference $L$ of the cylinder and finite bond dimension $\chi$ introduce a cutoff length scale Schollwöck (2011), which prevents us from observing the power-law behavior in the correlation function. Therefore, at long distances along the cylinder, the DMRG correlation function of an operator ${\cal O}$ will always decay exponentially with certain finite correlation length $\xi_{\cal O}$ in the simulation. Nevertheless, instead of trying to compare the scaling dimension $\Delta_{\mathcal{O}}$, we can determine the emergent symmetry by comparing the correlation lengths $\xi_{\mathcal{O}}$ between Néel (spin-1) and VBS (spin-0) excitations. In Zauner et al. Zauner et al. (2015), it has been found that the correlation length spectrum is inversely proportional to the energy of the excitations that mediate this correlation behavior. Thus, the correlation length spectra give access to the individual dynamics of different types of excitations.

Using the DMRG transfer matrix technique, one can readily obtain the correlation lengths $\xi$ along the cylinder. Moreover, since our DMRG simulation has an explicit $\mathrm{U}(1)_{z}$ symmetry, extracted correlation lengths are labeled by $S_{z}$ quantum numbers. Since the microscopic model has full $\mathrm{SO}(3)$ symmetry, there will be an exact degeneracy among correlation lengths with different quantum numbers. For example, a three-fold degeneracy among $S_{z}=0,\pm 1$ would imply that this correlation length corresponds to the excitation carrying the quantum number $S=1$ of the $\mathrm{SO}(3)$ symmetry. In this way, we can identify the $\mathrm{SO}(3)$ spin quantum number of each operator appeared in the correlation length spectrum.

In Fig. 4, we plot the inverse of the largest correlation lengths for spin-singlet and triplet operators. The spin-singlet operator corresponds to the pVBS order parameter (or the monopole operator $\mathcal{M}$) at low-energy. The monopole operator is gapped in both the Néel and pVBS phases and only become gapless at the critical point. Therefore, the divergence of the spin-singlet correlation length can signal the onset of DQCP in the thermodynamic limit. Indeed, we identified that the peak of the singlet correlation length exactly coincides with the critical point extracted from the singularity of the energy derivative (Fig. 3). Furthermore, the correlation length $\xi_{S=0}$ at the critical point increases as the bond dimension $\chi$ increases. Although here we present only the correlation length spectrum at $L=10$ and $\chi=4000$, we also performed numerical simulations for different system sizes and obtained the result that the critical point summarized in Tab. 1 coincides with the peak of the spin singlet correlation length. Therefore, we can infer that the transition is induced by the proliferation of monopoles, consistent with DQCP physics.

Moreover, to contrast our simulation result in the Shastry-Sutherland lattice with the conventional $\mathrm{O}(3)$ Wilson-Fisher transition between an explicitly dimerized VBS and Néel order phases Matsumoto et al. (2001), we performed the iDMRG simulation for a 2D $J_{1}$-$J^{\prime}_{1}$ model with the antiferromagnetic Heisenberg coupling $J_{1}$ on nearest neighbor bonds together with the coupling $J^{\prime}_{1}$ on a fixed set of dimer covering bonds (See Appendix. B), and hence a unique VBS pattern is pinned by $J^{\prime}_{1}$. Indeed, for this model, one can observe that the correlation length of the spin-triplet sector is always larger than the correlation length of the spin-singlet sector across the phase transition. This agrees with the picture discussed in Ref. Sachdev and Bhatt (1990); Kotov et al. (2009), where the transition is triggered by the condensation of spin-triplet excitations (triplon) from the VBS phase. Within the mean-field theory framework, it was shown that the energy of the spin-triplet excitation $E_{\textrm{triplon}}$ is smaller than the energy of the spin-singlet excitation $E_{\textrm{singlon}}$ throughout the whole transition. Under the further assumption that the singlon and triplon have the similar characteristic velocity $v$ Zauner et al. (2015), one expect the aforementioned ordering of correlation lengths $\xi\sim(v/E)$. Thus, the inversion of the magnitude of correlation lengths for spin singlet and triplet operators at the transition signifies the unconventional feature of the DQCP in the Shastry-Sutherland lattice. For a detailed analysis regarding DQCP physics, see Sec. IV.

Finally, we remark that the spin-singlet correlation length $\xi_{S=0}$ and the spin-tiplet correlation length $\xi_{S=1}$ approaches each other at the (finite size) critical point as we increase the system size. At $L=6$ and $\chi=4000$, the ratio \xi_{S=1}/\xi_{S=0}\big{|}_{\textrm{crit}}=0.33, but at $L=10$ and $\chi=4000$, \xi_{S=1}/\xi_{S=0}\big{|}_{\textrm{crit}}=0.95. From this trend, we expect to have $\xi_{S=1}/\xi_{S=0}=1$ in the thermodynamic limit, which indicates that the spin-singlet and spin-triplet excitations will become degenerate at the critical point, forming the four-component vector representation of a larger $\mathrm{O}(4)$ symmetry group. Put differently, we can observe that the crossing point of $\xi_{S=0}=\xi_{S=1}$ in Fig. 4 approaches to the critical point as we increase the system size, consistent with the emergent $\mathrm{O}(4)$ symmetry relating Néel and VBS order parameters.

III Symmetry Analysis

The field theory of DQCP, the so called “non-compact” $\text{CP}^{1}$ theory, has the following form Motrunich and Vishwanath (2004); Senthil et al. (2004a)

[TABLE]

where a two-component complex spinon ${z}=(z_{1},z_{2})^{T}$ is coupled to $\mathrm{U}(1)$ gauge field $a$. On top of this critical theory, one can have additional terms depending on symmetry of the system. The Shastry-Sutherland lattice has a $p4g$ space group symmetry, as shown in Fig. 1(a). The lattice respects two glide reflection $G_{x},G_{y}$ and two diagonal reflection $\sigma_{xy},\sigma_{x\bar{y}}$ symmetries as illustrated in Fig. 1(a). The glide reflections and the (spinful) time-reversal $\mathcal{T}$ symmetries can combine into composite symmetries $\mathcal{T}G_{x},\mathcal{T}G_{y}$, dubbed as the time-reversal glide symmetries. Note that glide-reflection symmetry is also broken in the Néel phase, while the time-reversal glide is not. Therefore, relative to the Néel phase, it is proper to think about the pVBS phases as to break the time-reversal glide symmetries.

To define the symmetry transformations more conveniently, we consider an ideal version of the Shastry-Sutherland lattice on a regular square lattice without distortion, as shown in Fig. 5. This does not change the symmetry group but allows us to label every site by the Cartesian coordinate $(x,y)$ conveniently (where $x,y\in\mathbb{Z}$). The length of the nearest Cu-Cu bond is set to 1, such that the unit cell is of the size $2\times 2$ (and hence the lattice constant is 2 here). With this, we can define the glide reflections $G_{x}:(x,y)\to(x+1,-y)$ and $G_{y}:(x,y)\to(-x,y+1)$, the diagonal reflections $\sigma_{xy}:(x,y)\to(y,x)$ and $\sigma_{x\bar{y}}:(x,y)\to(-y+1,-x+1)$, as well as the translations $T_{x}:(x,y)\to(x+2,y)$ and $T_{y}:(x,y)\to(x,y+2)$. Together, they generate the $p4g$ space group. The $p4g$ space group also contains a 90∘ rotation symmetry $C_{4}:(x,y)\to(-y+2,x-1)$ with respect to the center of the plaquette without a diagonal bond.

On each site $i=(x,y)$, we define the spin operator ${\bm{S}}_{i}=(S_{i}^{x},S_{i}^{y},S_{i}^{z})$, whose symmetry transformations are listed in Tab. 2 (assuming the spin rotation is not locked to the spatial rotation in lack of the spin-orbit coupling). The time-reversal symmetry $\mathcal{T}$ is also included, which can flip all spin components. The Néel order parameter ${\bm{n}}=(n_{x},n_{y},n_{z})$ and VBS order parameters ${\bm{v}}=(v_{x},v_{y})$ are defined as

[TABLE]

where each prefactor translates into the momentum $(\pi,\pi)$, $(\pi,0)$ and $(0,\pi)$ respectively in $k$-space, as marked out in Fig. 1(b), and the displacement vectors are defined to be $\hat{x}=(1,0)$ and $\hat{y}=(0,1)$. The VBS order parameters $v_{x},v_{y}$ can be combined to form the following operator given a certain gauge choice (Appendix.A)

[TABLE]

known as the monopole operator in the literature Senthil et al. (2004a); Levin and Senthil (2004), which corresponds to a hegehog-monopole of the Néel order parameter ${\bm{n}}$ in the spacetime. The monopole event changes the skyrmion number of the ${\bm{n}}$ field configuration by $+1$, and it is equivalent to inserting $2\pi$-flux of $\mathrm{U}(1)$ gauge field $a$ in the CP1 theory Eq. (4). The symmetry properties of ${\bm{n}}$ and $\mathcal{M}^{\dagger}$ follows from those of the spin operator ${\bm{S}}_{i}$ and is summarized in Tab. 2. For a detailed derivation, see Appendix A. Apart from these discrete symmetries, there is also an $\mathrm{SO}(3)$ spin rotation symmetry, under which ${\bm{n}}$ transforms as a $\mathrm{SO}(3)$ vector.

The expectation value of the monopole operator $\langle\mathcal{M}\rangle$ defined in Eq. (6) serves as a unified order parameter for various types of VBS orders. Depending on the phase angle of $\langle\mathcal{M}\rangle$, the columnar VBS (cVBS) is described by $\langle\mathcal{M}\rangle\sim\pm e^{\pm\mathrm{i}\pi/4}$ and the plaquette VBS (pVBS) is described by $\langle\mathcal{M}\rangle\sim\pm 1$ (diamond plaquette) or $\langle\mathcal{M}\rangle\sim\pm\mathrm{i}$ (square plaquette) as illustrated in Fig. 6. On the other hand, $\langle\mathcal{M}\rangle\neq 0$ could also be interpreted as the condensation of monopoles. So the DQCP, as a transition from the Néel phase into the VBS phase, can be thought as driven by the VBS ordering or equivalently by the monopole condensation (starting from the Néel phase), which can be tuned by a monopole chemical potential $r$ in the Lagrangian as $r\mathcal{M}^{\dagger}\mathcal{M}$. The transition happens as $r$ changes sign. The condensation of monopole establishes the VBS order on the one hand and simultaneously destroys the Néel order on the other hand, due to a nontrivial topological term among the Néel and VBS order parameters, which was analyzed in details in Ref. Senthil et al., 2004a, b. This scenario provides a plausible description of a direct continuous transition between the Néel and VBS phases.

However, apart from the apparent tuning parameter term $r\mathcal{M}^{\dagger}\mathcal{M}$, we must also include other symmetry-allowed (multi-)monopole terms in the Lagrangian, which could crucially influence the properties of the DQCP. On the Shastry-Sutherland lattice, to the leading order, they take the form of

[TABLE]

Here we adopt the short-hand notations $\real\mathcal{O}=(\mathcal{O}+\mathcal{O}^{\dagger})/2$ and $\imaginary\mathcal{O}=(\mathcal{O}-\mathcal{O}^{\dagger})/(2\mathrm{i})$ for generic operator $\mathcal{O}$. Given the symmetry properties in Tab. 2, one can see that the single monopole term, no matter $\real\mathcal{M}$ or $\imaginary\mathcal{M}$, is forbidden by the glide reflection symmetry $G_{x}$ or $G_{y}$. Furthermore, the imaginary part of the double-monopole term $\imaginary\mathcal{M}^{2}$ is forbidden by the diagonal reflection symmetry $\sigma_{xy}$ or $\sigma_{x\bar{y}}$. Note that these symmetries exist in the critical theory as the spontaneous symmetry breaking has not yet occurred and the microscopic model has the symmetries.

The higher-order monopole terms ($\mathcal{M}^{4},\mathcal{M}^{6},...$) are expected to be less relevant and are therefore not included in Eq. (7) explicitly. Therefore the double-monopole term $\lambda_{2}\real\mathcal{M}^{2}$ is the most relevant monopole perturbation allowed on the Shastry-Sutherland lattice. Depending on its sign, the system will favor a square plaquette (or diamond plaquette) VBS order in the VBS phase, describe by the order parameter $\imaginary\mathcal{M}$ (or $\real\mathcal{M}$), if $\lambda_{2}>0$ (or $\lambda_{2}<0$), as demonstrated in Fig. 6. The square and diamond pVBS orders have distinct symmetry properties. Under the reflection symmetries $\sigma_{xy}$ and $\sigma_{x\bar{y}}$, $\imaginary\mathcal{M}\mapsto-\imaginary\mathcal{M}$ while $\real\mathcal{M}$ stays invariant (see Tab. 2), so the square pVBS spontaneously breaks the reflection symmetries while the diamond pVBS does not. Additionally, square-plaquette-centered $C_{4}$ rotation symmetry is spontaneously broken in the diamond pVBS, while it is not in the square pVBS. Therefore the two different pVBS orders will lead to different lattice distortions that are symmetry-wise distinguishable in the experiments in the X-ray/neutron diffraction or NMR Waki et al. (2007); Takigawa et al. (2010); Zayed et al. (2017).

Previous studies Koga (2000); Corboz and Mila (2013) as well as our iDMRG data show that the pVBS phase has a square plaquette order. Thus $\lambda_{2}>0$ should be relevant to our discussion of the pVBS-Néel transition in the Shastry-Sutherland model Eq. (1). $\lambda_{2}>0$ can be also argued based on the microscopic Hamiltonian in Eq. (1). In the analysis of the pVBS phases, there are two types of singlet plaquette configurations: $s$-wave and $d$-wave types Mambrini et al. (2006) represented by the following singlet pairing configurations in a plaquette:

[TABLE]

For simplicity, consider a single plaquette with four spin-1/2s on the corners. For both square and diamond plaquettes, there is a AFM coupling $J_{1}$ along plaquette sides; for the diamond plaquette, there is an additional AFM $J_{2}$ coupling across one diagonal. Then, the $s$-wave and $d$-wave singlet configuration has the energy $-2J_{1}+1/4J_{2}$ and $-3/4J_{2}$ respectively for the diamond plaquette. As the pVBS phase exists at a parameter regime $J_{1}/J_{2}\sim 0.7$, an estimation of the singlet configuration energy as listed in Tab. 3 indicates that the $s$-wave pairing in the square plaquette has the lowest energy. In the iDMRG simulation, the wavefunction indeed exhibits the $s$-wave pairing symmetry in the pVBS phase. A recent exact-diagonalization study Nakano and Sakai (2018) on the small cluster of spins ($N_{s}=40$) also reported that the phase next to the Néel order phase hosts a spin-spin correlation which contradicts to the $d$-wave singlet in a diamond plaquette. Therefore it is natural to have $\lambda_{2}>0$ in the phenomenological field theory.

At the first glance, it seems that the double-monopole term $\lambda_{2}$ in Eq. (7) is relevant and may destroy the DQCP. However, it is realized in Ref. Metlitski and Thorngren, 2017 that $r$ and $\lambda_{2}$ actually recombine into a new tuning parameter $\tilde{r}$ and a new relevant perturbation $\tilde{\lambda}_{2}$. In the case of $\lambda_{2}>0$, the Lagrangian $\mathcal{L}_{\mathcal{M}}$ in Eq. (7) can be written as

[TABLE]

with $\tilde{r}=r-\lambda_{2}$ and $\tilde{\lambda}_{2}=r+\lambda_{2}$. The parameter $\tilde{r}$ still drives a transition at $\tilde{r}=0$ (or equivalently $r=\lambda_{2}$), as shown in Fig.7 with a modified emergent symmetry. At the transition point, the relevant perturbation $\tilde{\lambda}_{2}=2\lambda_{2}>0$ simply gaps out the diamond plaquette pVBS fluctuation $\real\mathcal{M}$ from the low-energy sector, leaving the square plaquette pVBS fluctuation $\imaginary\mathcal{M}$ quantum critical. It is further argued that the pVBS fluctuation $\imaginary\mathcal{M}$ will become degenerate with the Néel fluctuation ${\bm{n}}$ at the critical pointMetlitski and Thorngren (2017), because the perturbations ${\bm{n}}^{4}$, ${\bm{n}}^{2}(\imaginary\mathcal{M})^{2},(\imaginary\mathcal{M})^{4},\cdots$ that can break the symmetry that rotates Néel and pVBS are all rank-four operators, which are expected to be irrelevant at the critical point. Therefore the Néel and VBS order parameters can combine into a $\mathrm{O}(4)$ vector $({\bm{n}},\imaginary\mathcal{M})$, manifesting an emergent $\mathrm{O}(4)$ symmetry. The remaining topological $\mathrm{O}(4)$ $\Theta$-term still ensures that the development of the pVBS order $\imaginary\mathcal{M}$ will simultaneously destroy the Néel order ${\bm{n}}$, establishing a direct pVBS-Néel transition with emergent $\mathrm{O}(4)$ symmetry.

IV Dangerously Irrelevant Scaling and Its Absence

In this section, we discuss a peculiarity of the DQCP in the Shastry-Sutherland lattice compared to the other DQCP scenarios that have been extensively discussed.Senthil et al. (2004b); Motrunich and Vishwanath (2004); Senthil et al. (2004a); You and You (2016); Metlitski and Thorngren (2017); Wang et al. (2017); Jian et al. (2017a, b) In the presence of the $\tilde{\lambda}_{2}$ term in Eq. (9), the $\mathrm{U}(1)$ symmetry of the monopole operator (which acts as $\mathcal{M}\to e^{\mathrm{i}\theta}\mathcal{M}$) is explicitly broken down to $\mathbb{Z}_{2}$ (at the lattice level). This $\mathbb{Z}_{2}$ symmetry can be identified as the glide reflection symmetry ($G_{x}$ or $G_{y}$) on the Shastry-Sutherland lattice, which can be further broken spontaneously in the pVBS phase by its order parameter $\expectationvalue{\imaginary\mathcal{M}}$. At the DQCP, this $\mathbb{Z}_{2}$ symmetry will be restored and combined with the $\mathrm{SO}(3)$ spin rotation symmetry to form the larger emergent $\mathrm{O}(4)$ symmetry, denoted as $\mathrm{SO}(3)\times\mathbb{Z}_{2}\to\mathrm{O}(4)$. Although the $\mathbb{Z}_{2}$ symmetry is restored at the DQCP, it is never further enlarged to the $\mathrm{U}(1)$ symmetry of monopole conservation, because the explicit symmetry breaking term $\tilde{\lambda}_{2}$ is relevant. The presence of the relevant coupling $\tilde{\lambda}_{2}$ leads to an important difference between the DQCP on the Shastry-Sutherland lattice with the more conventional DQCP on the square lattice.

To expose the differences and connections, let us briefly mention the other two lattices: the square lattice and the rectangular lattice. Due to the different lattice symmetries, the allowed leading monopole terms will be different, as summarized in Tab. 4. They will lead to different VBS orders and different properties of the DQCP. For example, on a rectangular lattice, the other double-monopole term $\lambda^{\prime}_{2}\imaginary\mathcal{M}^{2}$ is allowed but $\lambda_{2}\real\mathcal{M}^{2}$ is forbidden, which favors the horizontal/vertical cVBS order depending on $\lambda^{\prime}_{2}>0$ (or $\lambda^{\prime}_{2}<0$). The DQCP on the rectangular lattice has a similar emergent $\mathrm{O}(4)$ symmetry, which was carefully analyzed in Ref. Metlitski and Thorngren, 2017. However, on a square lattice, the four-fold rotational symmetry forbids all the double-monopole terms, leaving the quadruple-monopole term $\lambda_{4}\real\mathcal{M}^{4}$ most relevant, which favors cVBS (or pVBS) if $\lambda_{4}>0$ (or $\lambda_{4}<0$). In the absence of the double-monopole term, the DQCP on the square lattice has an even larger emergent symmetry $\mathrm{SO}(3)\times\mathbb{Z}_{4}\to\mathrm{SO}(5)$. In the easy-plane model, the lattice symmetry would be that of the square lattice, but the spin-rotation symmetry is reduced from $\mathrm{SO}(3)$ to $\mathrm{U}(1)$. Here, the symmetry enhancement would be $\mathrm{U}(1)\times\mathbb{Z}_{4}\to\mathrm{O}(4)$ Qin et al. (2017); Wang et al. (2017).

Although the DQCP on the square lattice with the easy-plane deformation has the same $\mathrm{O}(4)$ emergent symmetry as the DQCP on the rectangular or Shastry-Sutherland lattice, there is a crucial difference between them. On the square lattice, the $\mathrm{U}(1)\to\mathbb{Z}_{4}$ symmetry breaking term $\lambda_{4}\real\mathcal{M}^{4}$ is dangerously irrelevant, which enhances $\mathbb{Z}_{4}$ to $\mathrm{U}(1)$ at the DQCP. As we move away from the DQCP toward the VBS phase, the system will exhibit two different length scales: the spin correlation length $\xi_{\text{spin}}$ and the VBS domain wall width $\xi_{\text{VBS}}$ Senthil et al. (2004a); Shao et al. (2016). The system is critical at the length scale below $\xi_{\text{spin}}$, meaning that the spin correlation decays in a power-law. Beyond this length scale, however, the system is still not fully in the VBS phase. Because the $\real\mathcal{M}^{4}$ operator is irrelevant at the critical point, its coupling coefficient $\lambda_{4}$ has decreased under the renormalization group (RG) flow initially, as in Fig. 8(a). However, once RG flow goes beyond the length scale of $\xi_{\text{spin}}$, $\lambda_{4}$ starts to increase with the RG flow until it becomes strong enough to break the $\mathrm{U}(1)$ symmetry down to $\mathbb{Z}_{4}$ at the length scale $\xi_{\text{VBS}}$ Lou et al. (2007). A careful analysis Senthil et al. (2004a) showed that $\xi_{\text{VBS}}\sim\xi_{\text{spin}}^{(\Delta-1)/2}$, where $\Delta>3$ (irrelevant) is the scaling dimension of $\real\mathcal{M}^{4}$ 111Beyond the length scale $\xi_{\text{spin}}$, the system behaves like XY model with order parameter $\mathcal{M}$. From the inspection of the XY model, the domain wall length scale is given by $\xi_{\text{VBS}}\sim\lambda_{4}^{-1/2}$, where $\lambda_{4}$ is a renormalized coefficient for $\real\mathcal{M}^{4}$ at that length scale. From the generic scaling argument, we have $\xi_{\text{VBS}}\sim Lf(\lambda_{4}L^{3-\Delta})$, where $L$ can be any length scale within the scaling regime. At the length scale $L=\xi_{\text{spin}}$ where the system behaves like XY model, we already know that $\xi_{\text{VBS}}\sim\lambda_{4}^{-1/2}$. Thus, we can deduce $f(x)\sim x^{-1/2}$, and derive that $\xi_{\text{VBS}}\sim\xi_{\text{spin}}^{(\Delta-1)/2}$. Thus, $\xi_{\text{VBS}}$ grows more rapidly than $\xi_{\text{spin}}$ as we approach the critical point.

However, on the Shastry-Sutherland lattice, there is no symmetry enhancement from $\mathbb{Z}_{2}$ to $\mathrm{U}(1)$ at the DQCP because $\mathrm{U}(1)\to\mathbb{Z}_{2}$ symmetry breaking term $\tilde{\lambda}_{2}(\real\mathcal{M})^{2}$ is relevant. As a result, there is no dangerously irrelevant RG flow around the DQCP. As soon as we move away from the critical point, $\mathbb{Z}_{2}$ anisotropy is there to break the emergent $\mathrm{O}(4)$ symmetry down to the microscopic $\mathrm{SO}(3)\times\mathbb{Z}_{2}$ symmetry. Thus, there does not exist a separation of length scales nor scenario expected in Fig. 8. The singlet gap (the gap of pVBS fluctuation) should open up immediately as we tune $\tilde{r}$ away from the critical point in the thermodynamic limit.

What is the consequence of all these observations? In the study of finite size systems, the RG flow should stop at a certain point beyond the system size. This means that the dangerously irrelevant scaling can significantly affect the correlation behavior or excitation spectrum in small system sizes, as illustrated in Fig. 8(b). Due to the dangerously irrelevant scaling, near the DQCP towards the VBS side, the correlation length of the VBS order parameter fluctuation ($S=0$) should be larger than the correlation length of the Néel order parameter fluctuation ($S=1$). However, in the thermodynamic limit, such a behavior is not guaranteed as the eventual fate under the RG flow is often difficult to understand. For example, in the mean-field treatment of the VBS phase, it has been argued that the spin-triplet excitation (triplon) has a lower energy than the singlet excitation (singlon) Kotov et al. (2009). On the other hand, in the finite size systems, the dangerously irrelevant scaling enforces the region with $\xi_{S=0}>\xi_{S=1}$ in Fig. 8(b) to appear regardless of the eventual RG behavior in the VBS phase. Moreover, since we consider a quasi two-dimensional system in the iDMRG simulation, the finite circumference size can also affect the behavior in Fig. 8. In principle, the Mermin-Wagner theorem prevents the spontaneous symmetry breaking of a continuous symmetry in dimensions $D\geq 2$ Mermin and Wagner (1966). In other words, in a quasi two-dimensional system, the strong fluctuation of the continuous order parameter (e.g. SO(3) spin-rotation) would gap out the system and reduce the correlation length for the associated Goldstone bosons (e.g. magnons). As the disordering effect would be stronger for the physical SO(3) spin-rotation symmetry than for the emergent $\mathrm{U}(1)$ symmetry enhanced from $\mathbb{Z}_{4}$, we would again expect the parameter region of $\xi_{S=0}>\xi_{S=1}$ to appear near the DQCP.

By contrast, in the Shastry-Sutherland model, the relevant $\mathbb{Z}_{2}$ perturbation can always gap out the monopole fluctuation away from the critical point. To confirm this, we performed the iDMRG simulation on the spin-1/2 $J_{1}$-$J_{2}$ model at the same circumference size $L=10$ and compared the correlation length spectra between two models. In Fig. 9(a), we observe that $\xi_{S=0}$, the correlation length of a $S=0$ local excitation, is always larger than $\xi_{S=1}$ in the entire VBS phase on the square lattice. However, In Fig. 9(b), $\xi_{S=0}$ is larger than $\xi_{S=1}$ only at the transition point and immediately becomes smaller than $\xi_{S=1}$ as we tune the system away from the critical point towards the VBS phase. This is one of the non-trivial prediction from the presence of relevant anisotropy operator at a finite-system size simulation. In the Néel ordered phase, these two models exhibit very similar correlation length spectra. For more details, see Appendix. B.

Our numerical result for the square $J_{1}$-$J_{2}$ model aligns with the recent work Wang and Sandvik (2017) on the finite-DMRG simulation which calculated first several excited states with different spin quantum numbers. If we replace the excitation energy in Ref. Wang and Sandvik (2017) with $\xi^{-1}$, we obtain the same crossing behavior. This can be justified by the fact that when a local excitation has a mass gap $m$, the correlation function mediated by the local excitation decays as $\sim e^{-mr}$, thus $m\propto\xi^{-1}$. Therefore, our theoretical scenario explains why a local 222In the VBS phase, there are four degenerate ground states, whose degeneracy is lifted by the finite system size. These nearly degenerate excited states should be distinguished from the other excitations that can truly be ‘local’. spin-singlet excitation in the Ref. Wang and Sandvik (2017) has lower energy than a spin-triplet excitation around the critical point. It also elucidates the reason why the previous DMRG results of the $J_{1}$-$J_{2}$ model on the square lattice Gong et al. (2014); Wang et al. (2013); Mambrini et al. (2006); Capriotti and Sorella (2000); Takano et al. (2003); Murg et al. (2009) were unable to identify the nature of the VBS phase without applying a pinning field. Because all previous numerics were also performed in the similar system size, they were in the regime where the VBS order parameter fluctuations were severe, disallowing one to confirm whether it is a plaquette or columnar VBS. Therefore, we conclude that the absence or presence of an irrelevant operator is essential to understand physical observables in any finite size system.

In summary, we note that previously the two promising scenarios to realize DQCP with spin 1/2 was either the SO(3) symmetric or easy plane deformations, both with four-fold degenerate VBS orders. Here, however, we have the two-fold degenerate VBS order. Nevertheless, as discussed in Ref. Metlitski and Thorngren (2017), the two-fold monopole with SO(3) symmetry is equivalent to the easy plane deformation, if an enlarged SO(5) symmetry is assumed in the absence of these perturbations. Furthermore, we remark that the scaling behavior of dangerously irrelevant operator enables us to understand multiple observations made in previous numerical simulations which is absent in the Shastry-Sutherland model here.

V Spectral Signatures of DQCP

A hallmark of the DQCP is the emergence of deconfined spinons at the critical point, which entails distinct features in both the magnon and phonon excitation spectra that can be probed in INS or RIXS experiments. To better appreciate the predicted spectral features at low-temperature around the pVBS-Néel transition, we need to first understand the background elastic scattering signal of $\mathrm{SrCu_{2}(BO_{3})_{2}}$ in its high-temperature paramagnetic phase without any symmetry breaking. We will focus on the scattering of neutrons or photons off of the copper sites. As shown in Fig. 1(a), there are four copper sites in each unit cell, coordinated at ${\bm{r}}_{A}=(1+\delta,1+\delta)/2$, ${\bm{r}}_{B}=(1-\delta,3+\delta)/2$, ${\bm{r}}_{C}=(3+\delta,1-\delta)/2$, ${\bm{r}}_{D}=(3-\delta,3-\delta)/2$ respectively, with the distortion parameter given by $\delta=0.544$ according to Ref. Miyahara and Ueda, 1999. In the high-temperature paramagnetic phase, an elastic scattering experiment will reveal lattice diffraction peaks at a set of momenta ${\bm{Q}}=\pi(H,K)$ ($H,K\in\mathbb{Z}$, note that the lattice constant is 2 in our convention) with the amplitude given by $S({\bm{Q}})=\int\differential^{2}{\bm{r}}\rho({\bm{r}})e^{\mathrm{i}{\bm{Q}}\cdot{\bm{r}}}\simeq\sum_{a=A,B,C,D}e^{\mathrm{i}{\bm{Q}}\cdot{\bm{r}}_{a}}$, where $\rho({\bm{r}})$ can represent either the electron density from Cu orbitals (which scatters X-ray photons) or the density of Cu nuclear (which scatters neutrons). The corresponding intensity $|S({\bm{Q}})|^{2}$ is plotted in Fig. 1(b). Notably, there are extinction points in the diffraction pattern, protected by the glide reflection symmetry. Note that the system at the DQCP has the glide-reflection symmetries $G_{x}$ and $G_{y}$ although both the Néel and pVBS phases do not. 333The Néel ordered phase has the time-reversal glide ${\cal T}G_{x,y}$, while the pVBS phase has the time-reversal ${\cal T}$. At the DQCP, we have both symmetries, thus the composition ${\cal T}\circ{\cal T}G_{x,y}=G_{x,y}$ should be there. The glide reflections $G_{x}$ and $G_{y}$ act as lattice translations by a half lattice constant followed by the reflections about the translation directions, as illustrated in Fig. 1(a). The fact that the density distribution $\rho({\bm{r}})$ at equilibrium respects all lattice symmetries (including $G_{x}$ and $G_{y}$) implies that $\rho(x,y)=\rho(x+1,-y)=\rho(-x,y+1)$. They impose the following constraints on the scattering amplitude

[TABLE]

which implies the extinction of diffraction peaks at ${\bm{Q}}\in\pi(2\mathbb{Z}+1,0)$ or $\pi(0,2\mathbb{Z}+1)$, as marked out in Fig. 1(b). In general, $\rho({\bm{r}})$ can describe the spatial pattern of any scatterers that interact with the probing particle (e.g. X-ray or neutrons). Eq. (10) holds under the assumption that the scatterer field $\rho({\bm{r}})$ is even (symmetric) under glide reflection symmetry. This constraint can be generalized to other type of scatterers including magnetic fluctuations at finite frequency. For example, the destruction of the scattering amplitude at these extinction points could extend to finite-frequency inelastic scattering, as long as no scatterer at that energy scale breaks the glide reflection symmetry. However, as we lower the temperature and approach the pVBS-Néel transition, certain glide-reflection-breaking excitations (meaning that the scatterer is odd under the glide reflection) may emerge at low energy as part of the quantum critical fluctuation. Indeed, we will show that the emergent $\mathrm{SO}(4)$ conserved current fluctuation and the pVBS order fluctuation are examples of glide-reflection-breaking scatterers, which will become critical at the DQCP and “light up” the extinction points. They will provide unique signatures of the DQCP that are also relatively easy to resolve in experiments, as there are no background scattering signals at the extinction points.

The proposed spectral signatures at the pVBS-Néel transition is most convenient to analyze using a fermionic spinon theory for the DQCP, which has been shown to be equivalent (dual) to the conventional CP1 theoryYou et al. (2017); Assaad and Grover (2016b); Ma et al. (2018); Xu et al. (2018). In the fermionic spinon theory, the spin operator ${\bm{S}}_{i}=\frac{1}{2}f_{i}^{\dagger}{\bm{\sigma}}f_{i}$ is fractionalized into fermionic spinons $f_{i}=(f_{i\uparrow},f_{i\downarrow})^{\intercal}$, which are then placed in the $\pi$-flux stateRantner and Wen (2002); Hermele et al. (2005); Senthil and Fisher (2006); Ran and Wen (2006) described by the following mean-field Hamiltonian

[TABLE]

where $\langle ij\rangle$ runs over all the nearest neighbor bonds and $[ijkl]$ runs over all the shaded diamond plaquette depicted in Fig. 10(a). The spinon hopping amplitude $t_{ij}$ takes $t_{ij}=-t$ on the bonds highlighted in red in Fig. 10(a) and $t_{ij}=t$ on the rest of the bonds, such that the spinon sees $\pi$-flux threading through each plaquette. The phenomenological parameter $t\sim J_{1}$ sets the energy scale of the spinon, which is expected to be of the same order as the spin interaction strength $J_{1}$. The $\pi$-flux state model Eq. (11) was originally proposed as an example of algebraic spin liquidsRantner and Wen (2002); Hermele et al. (2005); Senthil and Fisher (2006); Ran and Wen (2006) (including the DQCP as a special case). It is recently confirmed via quantum Monte Carlo (QMC) simulationsAssaad and Grover (2016b); Ma et al. (2018); Xu et al. (2018) that the $\pi$-flux state actually provides a pretty good description of the spin excitation spectrum at the DQCP.

The $\pi$-flux state mean-field ansatz determines a projective symmetry group (PSG)Wen (2002) that describes how the spinon should transform under the space group symmetry, as concluded in the last column of Tab. 2. It is found that the spinon hopping along the dimer bonds are forbidden by the $\sigma_{xy},\sigma_{x\bar{y}}$ symmetry, which is not broken at the DQCP. Therefore the effect of $J_{2}$ can only enter the Hamiltonian as a four-fermion interaction to the leading order, given by the $g$ term in Eq. (11). The $g$ term describes the spinon interaction around each diamond plaquette $[ijkl]$ where the site indices $i,j,k,l$ are arranged according to Fig. 10(b). It turns out that the interaction $g$ is directly related to the $\lambda_{2}\real\mathcal{M}^{2}$ term given that the monopole operator $\mathcal{M}^{\dagger}\sim(v_{x}+v_{y})+\mathrm{i}(v_{x}-x_{y})$ can be written in terms of VBS order parameters $v_{x},v_{y}$, which are further related to fermionic spinons via $v_{x}\sim(-)^{x}f_{i+\hat{x}}^{\dagger}f_{i}+\text{h.c.}$ and $v_{y}\sim(-)^{x+y}f_{i+\hat{y}}^{\dagger}f_{i}+\text{h.c.}$, such that the interaction can be written as $gv_{x}v_{y}\sim g\real\mathcal{M}^{2}$. Therefore $g>0$ would correspond to $\lambda_{2}>0$, which favors the square plaquette VBS order $\imaginary\mathcal{M}$.

Let us ignore the interaction $g$ for a moment. By diagonalizing the spinon hopping Hamiltonian, we find the spinon dispersion $\epsilon_{{\bm{k}}}=\pm 2t\sqrt{\cos^{2}k_{x}+\cos^{2}k_{y}}$, which gives rise to 4 Dirac fermions (or equivalently 8 Majorana fermions) at momentum $(\pi/2,\pi/2)$. Note that the distortion parameter $\delta$ does not affect the spinon band structure. Naively, the spinon mean-field theory has an emergent $\mathrm{SO}(8)$ symmetry rotating among the 8 low-energy Majorana fermions. However, a $\mathrm{SU}(2)$ gauge structure must be introduced with the fractionalization, which reduces the emergent symmetry to $\mathrm{SO}(5)$. The interaction term $g$ plays an important role to further break the $\mathrm{SO}(5)$ symmetry explicitly to $\mathrm{O}(4)$, matching the emergent symmetry observed in the DMRG simulation.

Based on the fermionic spinon mean-field ground state, the spin excitation spectrum $S(\omega,{\bm{q}})$ has been calculated in Ref. Assaad and Grover, 2016b; Ma et al., 2018; Xu et al., 2018 and is reproduced in Fig. 11(a) for illustration, where the high symmetry points $\Gamma,X,M$ are defined in Fig. 1(b). The lower edge of the spinon continuum is given by $\omega_{\text{min}}({\bm{q}})=\min_{{\bm{k}}}|\epsilon_{{\bm{k}}+{\bm{q}}}-\epsilon_{{\bm{k}}}|$, which reads

[TABLE]

as shown in Fig. 11(a). This provides us a way to estimate the mean-field parameter $t$ from the experimentally measured spin excitation spectrum, by fitting this lower edge. This is a rather robust result as its shape is unaffected by the distortion parameter $\delta$.

Remarkably, a gapless continuum will appear on top of the extinction point $X$ (as well as $Y$ by $\sigma_{xy}$ symmetry) as seen in Fig. 11(a). As previously analyzed in Eq. (10), the spin excitation at the $X$ point must be odd under the glide reflection symmetry $G_{x}$. This symmetry constraint enforces that the gapless continuum should correspond to the emergent $\mathrm{SO}(4)$ conserved current fluctuation ${\bm{J}}_{y}={\bm{n}}\partial_{y}\imaginary\mathcal{M}-\imaginary\mathcal{M}\partial_{y}{\bm{n}}$ which involves both the Néel ${\bm{n}}$ and pVBS $\imaginary\mathcal{M}$ order parameters and has been thoroughly studied in Ref. Hermele et al., 2005; Ma et al., 2018, 2018. Since ${\bm{n}}$, $\partial_{y}$ and $\imaginary\mathcal{M}$ are all odd under the glide reflection $G_{x}$ (that preserves the extinction point $X$), the current operator ${\bm{J}}_{y}$ is also odd under $G_{x}$. Therefore Eq. (10) does not hold anymore, such that the fluctuation of ${\bm{J}}_{y}$ can appear at the extinction point $X$ as a spinon continuum in the magnon channel (as ${\bm{J}}_{y}$ carries spin-1). On the other hand, fluctuations like ${\bm{n}}\imaginary\mathcal{M}$ is not allowed to appear at the $X$ point in the spin excitation spectrum because the bound state ${\bm{n}}\imaginary\mathcal{M}$ is even under glide reflection (cf. Eq. (10)). The ${\bm{J}}_{y}$ continuum will only become gapless if both the Néel and pVBS fluctuations are gapless, which only happens at the DQCP. Protected by the emergent $\mathrm{O}(4)$, the conserved current operator must have a scaling dimension exactly pinned at 2, which indicates that the magnon spectral weight at the extinction point $X$ must increase with the frequency linearly (at below the energy scale of $J_{1}$)

[TABLE]

as shown in Fig. 11(b). We propose this as a hallmark feature of the spin fluctuation at the pVBS-Néel transition in $\mathrm{SrCu_{2}(BO_{3})_{2}}$. Confirmation of this linear frequency growth of the spectral weight will provide direct evidence for the emergent $\mathrm{O}(4)$ symmetry at the DQCP.

Apart from the features in the spin excitation (magnon) spectrum, the DQCP also introduces new features to the phonon spectrum, due to the pVBS-phonon coupling. The pVBS order has a linear coupling to the lattice displacement as its representation under the lattice symmetry group matches with that of strain fields

[TABLE]

where $u_{x}$ (or $u_{y}$) is the lattice displacement in the $x$ (or $y$) direction with a momentum $(\pi,0)$ (or $(0,\pi)$). It is crucial that $u_{x}$ and $u_{y}$ here are not acoustic phonon modes around momentum $(0,0)$ in a continuum theory, otherwise they can only enter the field theory in the form of derivatives $\partial_{i}u_{j}$. The coupling is allowed by lattice symmetry and is evident from the pVBS induced lattice distortion as demonstrated in Fig. 6(a). This leads to a hybridization between phonon and pVBS fluctuations. As the pVBS fluctuation becomes critical (gapless) at the DQCP, the quantum critical fluctuation will also appear in the phonon spectrum due to the hybridization effect, as illustrated in Fig. 11(c) (see Appendix C for detailed analysis). As the pVBS order carries the momentum $(\pi,0)$ and $(0,\pi)$, the phonon continuum will also get softened at these momenta, which happen to be the extinction points $X$ and $Y$ of the lattice diffraction pattern. Note that the pVBS fluctuation is odd under glide reflection, therefore new phonon continuum is allowed to emerge at the extinction points with the spectral weight diverging at low-frequency following a power-law,

[TABLE]

The anomalous dimension $\eta$ should match that of the pVBS order parameter at the DQCP, which, by the emergent $\mathrm{O}(4)$ symmetry, is also the same $\eta$ of the Néel order parameter. Based on the QMC simulations in Ref. Qin et al., 2017; Zhang et al., 2017, $\eta$ has been estimated to be $\eta=0.13\sim 0.3$. Observation of such critical phonon fluctuations at the extinction points with a power-law divergent spectral weight as shown in Fig. 11(d) will be another direct evidence of DQCP.

In conclusion, extinction points are protected by the glide reflection symmetry, but both the conserved current and the pVBS order parameter breaks the glide reflection. Their critical fluctuations are therefore allowed to appear at the extinction point. This is rather a nice property that there will be no background signals form lattice diffraction, which makes these spectral signatures of DQCP more easier to resolve in experiments. Our analysis indicates that the conserved current fluctuation (which carries spin-1) should appear in the magnon spectrum and the pVBS fluctuation (which carries spin-0) should appear in the phonon spectrum. Observation of these critical fluctuations in the scattering experiment will be strong evidence for the potential realization of DQCP in $\mathrm{SrCu_{2}(BO_{3})_{2}}$.

VI Effects of Interlayer Coupling

So far, we have discussed the DQCP physics assuming that the system is two-dimensional. However, the material realization of Eq. (1), SrCu2(BO3)2, has a three-dimensional structure which is a stack of the Shastry-Sutherland lattices with a relative shift given by the lattice vector (1,1) between neighboring layers. Each layer is separated by the layer of oxygen, but due to the super-exchange term mediated by the oxygen, there is a small interlayer anti-ferromagnetic interaction $J_{3}\sim 0.1J_{2}$ Miyahara and Ueda (2003) between the spin-1/2s located at the crossed dimers in Fig. 12. Since the stacking structure preserves $G_{x,y}$ and $\sigma_{xy,x\bar{y}}$ symmetries, previous monopole analysis still holds for each two-dimensional layer. Here, we would like to better understand the effect of the three-dimensional interlayer coupling to the DQCP physics.

To analyze, we first consider coupling layers of spin systems at the conventional DQCP point with the emergent $\mathrm{SO}(5)$ symmetry (other than the $\mathrm{O}(4)$ case in the previous discussion). Here, each layer is described by the following nonlinear sigma (NLSM) model in Euclidean spacetime Senthil and Fisher (2006); Tanaka and Hu (2005):

[TABLE]

where $\Phi^{(l)}=({\bm{n}}_{x},{\bm{n}}_{y},{\bm{n}}_{z},\imaginary\mathcal{M},\real\mathcal{M})$ is the order parameter of the $l$-th layer and $\Gamma_{\textrm{WZW}}[\Phi^{(l)}]$ is a Wess-Zumino-Witten (WZW) term at level $k=1$. With the interlayer coupling, the total action is given by the following general form Bi et al. (2015)

[TABLE]

where the interlayer coupling coefficients $g^{ab}$ can be arranged into a matrix $g$. Given that the scaling dimension of $\Phi$ is estimated to be around $\Delta_{\Phi}\simeq 0.6$ in the literatureLou et al. (2009); Harada et al. (2013); Block et al. (2013); Pujari et al. (2013); Sreejith and Powell (2014, 2015); Pujari et al. (2015); Dyer et al. (2015); Nahum et al. (2015b); Nakayama and Ohtsuki (2016), the interlayer coupling is expected to be relevant (as $2\Delta_{\Phi}<3$), locking the order parameters across different layers. Then, the coefficient matrix $g$ becomes important since the sign of its determinant $\det g$ crucially affects how the topological (WZW) term from each layer will be added up (cf. Bi et al. (2015)). More precisely, if the sign of $\det g$ is positive (negative), the WZW terms are added up in a uniform (staggered) manner. The reason is that we can always redefine $\Phi^{(l)}$ in alternate layers to make $g$ positive definite, but the price to pay will be that the WZW term will change sign alternatively if the original $\det g$ is negative. Note that the procedure requires the topological term to be invariant under the change of origin (translation symmetric), i.e. $T_{x,y}:\Gamma_{\textrm{WZW}}[\Phi]\mapsto\Gamma_{\textrm{WZW}}[\Phi]$. Otherwise, the addition or subtraction of topological terms across layers would not be well-defined due to the arbitrariness of the choice of origin across the layers for the locked order parameter.

For example, consider coupling the two-dimensional square lattice $J$-$Q$ models Sandvik (2007) into a 3D cubic lattice with AFM spin-spin interaction along vertical bonds. Each layer putatively realizes the DQCP physics with $\Phi=({\bm{n}},v_{x},v_{y})$ describing Néel and cVBS order parameters. Under the AFM interlayer coupling, the coupling matrix $g$ has the sign structure of $g^{aa}=(-,-,-,+,+)$. This is because the vertical AFM coupling prefers ${\bm{n}}^{(l)}=-{\bm{n}}^{(l+1)}$ while the vertical plaquette ring exchange 444The AFM spin-spin interaction along vertical bonds induces the effective vertical plaquette ring exchange for the low energy subspace. favors $v_{x,y}^{(l)}=v_{x,y}^{(l+1)}$. In this case, $\det g<0$, so the WZW terms in neighboring layers tend to cancel each other. However the cancellation will not be exact, as the $\Phi$ field still admits (smooth) fluctuation over layers, this results in a residual topological term, namely a topological $\Theta$-term. Staggering a WZW-term at level $k$ would give a $\Theta$-term at $\Theta=\pi k$. Now the problem of the coupled DQCPs boils down to understand the fate of the $\mathrm{SO}(5)$ NLSM with $\Theta=\pi$ in (3+1)D. There are some hints from the fermionic parton analysis. One can consider fractionalizing the $\Phi$ vector to the bilinear from of a fermionic parton field $\psi$ following $\Phi_{a}\sim\bar{\psi}\mathrm{i}\gamma^{5}\Gamma^{a}\psi$, such that the $\psi$ fermion is in the $\mathrm{SO}(5)$ spinor representation (or the $\mathrm{Sp}(2)$ fundamental representation). The emergent gauge structure will be $\mathrm{SU}(2)$, which points to the $\mathrm{SU}(2)$ quantum chromodynamics (QCD) model in (3+1)D with $\mathrm{Sp}(2)$ flavor symmetry,

[TABLE]

Integrating out the fermion and gauge fluctuation is expected to reproduce the $\mathrm{SO}(5)$ NLSM with $\Theta=\pi(1+\operatorname{sgn}m)$, such that $\Theta=\pi$ is realized at $m=0$. However, the number of Dirac fermion flavors ($N_{f}=2$) is not enough to avoid a chiral symmetry breaking in 3D. Therefore, under the interlayer coupling, it is likely that the $\mathrm{SO}(5)$ DQCP flows to a discontinuous transition point induced by the quantum fluctuation. Considering that the $\mathrm{O}(4)$ DQCP from easy-plane anisotropy or rectangular deformation is descended from the $\mathrm{SO}(5)$ DQCP Metlitski and Thorngren (2017), breaking $\mathrm{SO}(5)$ down to $\mathrm{O}(4)$ makes the situation worse.

In a similar way, one can analyze the three-dimensional stacking of the Shastry-Sutherland lattice. If we allow the possibility of the diamond pVBS phases, the order parameter is written as $\Phi=({\bm{n}},\real\mathcal{M},\imaginary\mathcal{M})$ where $\real\mathcal{M}$ ($\imaginary\mathcal{M}$) represents the square (diamond) pVBS order parameter. Note that now each layer is shifted by $(1,1)$ vector (see Fig. 5) relative to the layer below, and the interlayer AFM coupling is given by Fig. 12(a) instead of the direct vertical coupling. In Fig. 12(b), we show two identical layers of Néel ordered pattern relatively shifted by $(1,1)$. If we focus on the four spins surrounding the diagonal bond crossing, we found that the two spins from the upper layer and the two spins from the lower layer are aligned oppositely, which is favored by the AFM interlayer spin exchange $J_{3}$. So the interlayer coupling prefers to lock the Néel order parameter in the same direction across the layer as ${\bm{n}}^{(l)}={\bm{n}}^{(l+1)}$. In Fig. 12(c), we show two identical diamond pVBS patterns displaced from each other. This configuration can gain the effective interlayer ring exchange energy induced by the $J_{3}$ coupling, which resonates the nearby dimers across the layer. Thus we conclude that $\real\mathcal{M}^{(l)}=\real\mathcal{M}^{(l+1)}$ is more favorable. In Fig. 12(d), we show two opposite square pVBS patterns displaced from each other. This configuration also gains interlayer ring exchange energy by resonating the dimers lying on top of each other. But this would require the square pVBS order parameter to be opposite between neighboring layers as $\imaginary\mathcal{M}^{(l)}=-\imaginary\mathcal{M}^{(l+1)}$. In conclusion, the interlayer coupling prefers $({\bm{n}}^{(l)},\real\mathcal{M}^{(l)},\imaginary\mathcal{M}^{(l)})=({\bm{n}}^{(l+1)},\real\mathcal{M}^{(l+1)},-\imaginary\mathcal{M}^{(l+1)})$. Here, instead of using $(v_{x},v_{y})$, we use $(\real\mathcal{M},\imaginary\mathcal{M})$ to parameterize the VBS order parameter, which is a basis change with a positive determinant. In this basis, the interlayer coupling matrix $g$ takes the sign structure of $g^{aa}=(+,+,+,+,-)$. As a result, we again have $\det g<0$, which implies that $\mathrm{SO}(5)$ DQCP would flow to the first order transition point in 3D. Since our $\mathrm{O}(4)$ scenario is considered to be a perturbed $\mathrm{SO}(5)$ DQCP (see Fig. 7), it is likely that the $\mathrm{O}(4)$ DQCP would also flow to the first order transition point.

If $\det g$ happens to be positive, the leading order effect is that the WZW terms would add up together, as the interlayer coupling $g$ tends to lock the order parameters $\Phi^{(l)}$ across the layers. Admittedly, the locking effect will become weaker at longer distance (along the perpendicular direction of layers), but we can still analyze the problem by first grouping the neighboring layers to a renormalized layer and then considering the residual coupling between renormalized layers. Across neighboring layers, the order parameters are expected to bind together, such that the renormalized model can be viewed as a NLSM with WZW term at large level $k$. The intuition from (0+1)D $\mathrm{O}(3)$ WZW term is that the large level limit corresponds to the large spin limit, where the quantum fluctuation of the order parameter is suppressed. Coupling spin-1/2 into a spin chain ferromagnetically in $S_{x,y,z}$-channels results in the ferromagnetic ground state with giant spin and classical spin wave excitations. In this limit, the low-energy physics can be captured within Landau-Ginzberg (LG) theory. Further adding different easy-axis anisotropies to the ferromagnetic spin chain will drive first-order transitions between different Ising ordered phases according to the LG theory. We conjecture that in higher dimension, similar effect will render each renormalized layer into a classical magnet which should be described within Landau-Ginzberg paradigm, such that the DQCP is not available. So in the presence of interlayer coupling, the Néel and VBS phases will likely be separated by a first order transition or intermediate coexisting phases. Thus, our analysis shows that in both $\det g>0$ and $\det g<0$ cases, the interlayer coupling would ultimately destabilize the DQCP, and drive it, for example, to a first order transition. However, our analysis also indicates that the $\det g<0$ case, corresponding to the real material, will have stronger quantum fluctuations, potentially leading to a weaker first order transition or a smaller region of coexisting order parameters than the $\det g>0$ case.

One additional remark is that while the interlayer Néel order coupling enters directly from the $J_{3}$ term, the interlayer VBS coupling arises from the higher-order perturbation (resonance) of the $J_{3}$ term. As a result, the critical point would be shifted to expand the Néel order phase. This is consistent with the phase diagram studied in Koga (2000), where the interlayer coupling drives a system into the AFM order and shrink down the VBS region.

VII Predictions for Experiment

Before discussing experimental consequences of the DQCP, we will make a few remarks on the nature of the plaquette VBS phase. In the Sec. III, we discussed two different possibilities for the pVBS phases (see Fig. 6): the square and diamond pVBS. While the square pVBS breaks the reflection symmetries $\sigma_{xy}$ and $\sigma_{x\bar{y}}$, the diamond pVBS breaks the empty-plaquette-centered $C_{4}$ rotation symmetry, see Fig. 5. As a result, when the system is at the pVBS phase, magnetic excitations initially degenerate under the $p4g$ symmetry would split differently depending on whether the plaquette is formed at a square or diamond. While our iDMRG simulation of the Shastry-Sutherland model points to the square pVBS phase, in the recent experiments on SrCu2(BO3)2 using INSZayed et al. (2017) and NMRWaki et al. (2007); Takigawa et al. (2010), the magnetic excitations in the pVBS phase seems to break the $C_{4}$ rotation symmetry, indicating the diamond pVBS phase. This discrepancy implies that the effective spin model for the real material could deviate from the Shastry-Sutherland model studied here. For example, three-dimensional interlayer coupling may induce some effective further-neighbor couplings beyond Shastry-Sutherland model. Therefore the type of pVBS phase to be stabilized at low energy is model dependent. Nevertheless, this does not affect to the DQCP scenario and the emergent $\mathrm{O}(4)$ symmetry because it only corresponds to a different sign of $\lambda_{2}$ in Eq. (7).

As discussed earlier, the DQCP naturally realizes a quantum spin liquid, a long sought after state of quantum magnets. Furthermore it realizes a particularly exotic variety - a critical spin liquid - with algebraically decaying correlations arising from the gapless emergent degrees of freedom. Moreover, an experimental realization of the DQCP would be a crucial manifestation of the many-body Berry phase effect that intertwines different order parameters. A dramatic experimental consequence of the DQCP is the emergent symmetry and resultant spectroscopic signatures expected from INS or RIXS. In particular, the model for SrCu2(BO3)2 studied here exhibits the $\mathrm{O}(4)$ emergent symmetry with two promising spectroscopic signatures at $X$-point in the Brillouin zone, summarized as the following:

• •

Magnon ($S=1$) Channel: This gives the information about the critical fluctuation of the emergent $\mathrm{O}(4)$ conserved current. As a result, the spectral intensity increases linearly with the frequency, $S\sim\omega$. If the emergent symmetry did not exist, there should not exist low energy spectral weight at this momentum. The deviation from the linear relation would give us a measure of how accurate the emergent symmetry is.

• •

Phonon ($S=0$) Channel: This gives the information about the pVBS order parameter fluctuation. For a given anomalous dimension $\eta_{\textrm{VBS}}$ of the pVBS order parameter, the spectral intensity diverges with the frequency, $S\sim\frac{1}{\omega^{2-\eta_{\textrm{VBS}}}}$. The DQCP scenario implies that the VBS order parameter is fractionalized, resulting in a non-zero $\eta_{\textrm{VBS}}$.

Moreover, the emergent symmetry implies that $\eta_{\textrm{VBS}}$ is equal to the anomalous dimension of the Néel order parameter fluctuation, $\eta_{\textrm{N\'{e}el}}$. However, the Néel order parameter fluctuation is located at the $M$-point, which has a pronounced Bragg peak in addition. In principle, the Bragg peak corresponds to spin-0 channel and it must be possible to extract $\eta_{\textrm{N\'{e}el}}$ and compare the $\eta_{\textrm{N\'{e}el}}$ with $\eta_{\textrm{VBS}}$ to tell whether the emergent $\mathrm{O}(4)$ symmetry exits. In the earlier work Miyahara and Ueda (2003), it is estimated from the low-$T$ magnetic susceptibility and heat capacity data that $J_{1}\approx 4.7\textrm{ meV}$, $J_{2}\approx 7.3\textrm{ meV}$, and $J_{3}\approx 0.7\textrm{ meV}$. Moreover, the Debye frequency of the acoustic phonon branch has been measured to be $\omega_{D}\approx 10\textrm{ meV}$ in Ref. Haravifard et al. (2012). Therefore, for experiments to confirm the theoretical predictions, it is required to have the energy resolution smaller than the milli-electron volt.

According to the present numerical simulation, the phase transition between the pVBS and Néel order can be a second-order or weakly-first order transition. The result does not contradict recent numerical work with a similar phase diagram where a first-order transition behavior was observed, because these models Serna and Nahum (2018); Zhao et al. (2018) are different from the microscopic Hamiltonian in Eq. (1), which is more likely to capture the couplings in the real material. Since the nature of the phase transition may be tunable, it is possible that the experiments on SrCu2(BO3)2 may realize the transition that is either a second order or a weakly first order with large correlation length.

Would all these predictions become meaningless if the transition is actually a weakly first order? In fact, even if the two-dimensional system hosts the DQCP, we argued in the previous section that the interlayer coupling might drive the system into a first-order transition point in three-dimension. Indeed, at the weakly first-order transition point, the system would have a finite excitation gap, and experimental spectroscopic data at zero temperature would deviate from Fig. 11 due to the absence of a gapless critical fluctuation. However, if we examine the system within the length scale smaller than the (large) correlation length $\xi$, the system would still exhibit the DQCP physics. In other words, if we only examine the system above the energy scale set by the correlation length $\omega>\omega_{\text{gap}}\sim 1/\xi$, we would observe the predicted spectral intensity trends in Fig. 11. It is our hope that future experiments will be able to use these results to clarify the physics behind the interesting properties of SrCu2(BO3)2.

VIII Conclusions

We studied a two dimensional $S=1/2$ model which captures key features of the Shastry-Sutherland material SrCu2(BO3)2. We obtained the phase diagram using numerical iDMRG simulations and observed a potentially continuous transition between a plaquette VBS state with two-fold degeneracy and a Néel ordered phase. The transition, studied using both numerical and field theoretical techniques, is proposed to be a deconfined quantum critical point, and we discussed its special features including the lack of a dangerously irrelevant scaling and an emergent O(4) symmetry. Concrete predictions are made for future experiments in SrCu2(BO3)2, where a pressure tuned transition between Néel order and a putative plaquette VBS state has already been reported Zayed et al. (2017). The predicted experimental signatures include the form of spectral intensity of spin singlet and spin triplet excitations at extinction points, which should be accessible in future resonant X-ray and neutron scattering experiments. These can provide a smoking gun signature of the deconfined criticality and emergent O(4) symmetry. Complications arising from the coupling between layers in the third dimension of the bulk material are briefly discussed, although further work in this direction is needed. We hope this study will trigger future experimental investigation of this quantum critical point in an interesting material, and more generally provide a road map for the experimental study of deconfined quantum criticality.

Acknowledgements.

We thank Y.-C. He, C. Wang, M. Takigawa, A. Sandvik, B. Zhao, L. Wang, Z.-Y. Meng, S.E. Han, and M. Metlitski for stimulating discussions. The DMRG numerics were performed using code developed in collaboration with Roger Mong and the TenPy collaboration. The computations in this paper were run on the Odyssey cluster supported by the FAS Division of Science, Research Computing Group at Harvard University. J.Y. Lee and A. Vishwanath were supported by a Simons Investigator Fellowship. S. Sachdev was supported by the National Science Foundation under Grant No. DMR-1664842.

Appendix A Monopole Transformation

In this section, we calculate how the monopole operator transforms under the symmetry of the Shastry-Sutherland lattice. To do this, we think of the Shastry-Sutherland lattice as the lattice being deformed from the square lattice with spin-1/2 per site. Starting from the $2+1$D antiferromagnetic ordered phase (Néel) of the square lattice, one can derive the action in terms of local Néel order parameter ${\bm{n}}(r_{n})$ in a path integral formalism,

[TABLE]

where ${\bm{n}}\sim\epsilon_{r}{\bm{S}}$ is a Néel order parameter, and $\epsilon_{r}=(-1)^{r_{x}+r_{y}}$ is a factor required for an alternating spin direction. In addition to the continuum action of a classical $O(3)$ nonlinear $\sigma$-model, there exists a Berry phase contribution due to the quantum nature of spin dynamics, which manifestly has the lattice origin. Haldane (1988); Read and Sachdev (1990) Let $\omega(n_{r})$ be a solid angle swept by a local Néel vector located at $r$ throughout the imaginary time evolution from 0 to $\beta$, with respect to the reference direction $n_{0}$ (See Fig. 13(b)). Then, $S_{B}=iS\sum_{r}\eta_{r}\omega({\bm{n}}_{r})$ where $S$ is the spin of an each site and $\eta_{r}=(-1)^{r_{x}+r_{y}}$ is an alternating phase factor coming from the antiferromagnetic nature of the magnetic order.

For any spatial slice, one can define the Skyrmion number Q(\tau)=\frac{1}{4\pi}\int\hat{{\bm{n}}}\cdot(\partial_{x}\hat{{\bm{n}}}\times\partial_{y}\hat{{\bm{n}}})\Big{|}_{\tau}, which is a topological invariant. Then, a monopole creation operator is defined as a topological defect at the spacetime point which changes $Q(\tau)$ by $+1$ across it. By the further duality mapping, this monopole in the NL$\sigma$M will be mapped into the monopole of the CP1 theory. As the center of the monopole cannot have a finite spin direction, the monopole is located at a dual lattice. When we have a monopole event in the spacetime, it must give a branch-cut structure on the image of $\omega(r)$ because $\omega(r)$ must change by $4\pi$ around the monopole location in the space. More intuitively, when monopole residing on the dual lattice encircles a single site, the imaginary time trajectories of all spins except the one encircled by the monopole would oscillate just back and forth, while the trajectory for the encircled one would entirely wind its $\omega$ by $4\pi$ upon the completion of encircling as in Fig. 13(b). Thus, one can view this problem as a monopole hopping around the dual 2D lattice in which each site in the original lattice gives a $4\pi S\eta_{r}$ flux through the plaquettes of the dual lattice. (Monopole is a charged object under this “flux” emanated from a spin-$S$.) The associated phase factor is independent of the exact imaginary time location of the monopole event. Thus for $S=1/2$ case, one can fix the system into a certain gauge and view this as if $\pm\pi/2$ Aharonov-Bohm phase factor gets accumulated for each hopping process for monopoles.

The monopole transformation rule is summarized in Tab. 5. $R^{\textrm{site}}_{\pi/2}$ is a site (spin) centered rotation. Note that it is important to fix the convention for the rotation center because translation symmetry is projective. Here, we choose a rotation to be defined with respect to the black spin in Fig. 13(a). $R^{\textrm{plaq}}_{\pi/2}$ is a plaquette centered rotation, which is defined with respect to the center of four spins in Fig. 13(a). Under the unit translations $T_{x,y}$, time reversal ${\cal T}$, and plaquette centered rotation $R^{\textrm{plaq}}_{\pi/2}$, the Néel order changes its sign. It means that the flux pattern is reversed under such transformations, thus we need to transform a monopole into an anti-monopole to compensate for the change. For reflections $\sigma_{x,y}$, although ${\bm{n}}$ does not flip, the definition of the Skyrmion number tells that we need to change the sign for the number of monopoles. Thus, a monopole transforms into an anti-monopole again. After figuring out whether a monopole transforms into a monopole or anti-monopole, we need to multiply it by an additional phase factor $\alpha_{g}$ to account for the Berry phase effect.

Assume the topological term $S_{B}$ is absent momentarily, which is how CP1 theory in Eq. (4) is derived. This is a usual practice because unlike the first term inside the parenthesis in Fig. 19, the second term, $S_{B}$, cannot be straightforwardly extended to the continuum field theory description. Under the absence of $S_{B}$, the monopole insertion operator $\mathcal{M}^{\dagger}$ just makes a global adjustment of the Néel order configuration to increase the Skyrmion number by one, without any additional phase factor.

However, we know from the existence of $S_{B}$ in the lattice description that the Berry phase effect is important. In order to take into account the Berry phase effect, we need to examine how the monopole transforms under each symmetry action. Under the active transformation where the coordinate system remains the same, we have

[TABLE]

where the action of $g$ on ${\cal M}^{\dagger}$ is determined by the previous argument on whether the monopole transforms into the monopole or anti-monopole upon the symmetry transformation.

To determine $\alpha_{g}$, we need to fix a gauge first. Fixing a gauge is important because a monopole is always created in a pair with an anti-monopole. Thus, the monopole event always has a reference point (anti-monopole) connected by the branch-cut. The Berry phase factor is independent of the way we draw the branch-cut because going around a single spin-lattice site gives a $2\pi$ phase. Let’s fix the reference gauge such that the monopole created at $(\bar{0},\bar{0})$ gives a Berry phase $\beta$. Then, for a generic coordinate $r$, inserting a monopole would give the Berry phase factor $\eta(r)\beta$ where $\eta(r)=1,i,-1,\textrm{ or }-i$ depending on whether the coordinates $(r_{x},r_{y})$ are (even, even), (odd, even), (odd, odd), (even, odd) Haldane (1988). This is shown explicitly in Fig. 13(c) as moving the monopole along the arrow gives an additional phase factor $i$. Inserting an anti-monopole at $r$ gives a phase factor $\eta^{*}(r)\beta^{*}$ since it would give an exactly opposite contribution to the Berry phase term by having $\omega\mapsto-\omega$. Now, by fixing $\beta=1$, the insertion of the monopole and anti-monopole at each dual lattice site simply gives a phase $\eta(r)$ and $\eta^{*}(r)$.

To illustrate the further procedure, let us consider two examples, $R^{\textrm{site}}_{\pi/2}$ and $T_{x}$. Under $R^{\textrm{site}}_{\pi/2}$, a monopole remains monopole, but its location changes as $r\mapsto R^{\textrm{site}}_{\pi/2}r$. By calculating the relative phase factor between the monopole created at $r$ and $R^{\textrm{site}}_{\pi/2}r$, we obtain its transformation rule in the continuum description:

[TABLE]

In the case of $T_{x}$, a monopole transforms into an anti-monopole under the symmetry action. At the same time, its location changes as $r\mapsto r+\hat{x}$. Following the similar procedure, we obtain the following rule:

[TABLE]

Following the similar analysis, we can obtain $\alpha_{g}$ for all symmetry transformations summarized in Tab. 5. In the case of time-reversal symmetry, we do not need an extra phase factor because time-reversal already complex-conjugates the phase factor of a monopole operator, which matches with the phase factor given by the anti-monopole at the same site. Moreover, any monopole condensation amplitude would be time-reversal symmetric because ${\cal T}:\expectationvalue{\mathcal{M}}\mapsto\expectationvalue{\mathcal{M}^{\dagger}}^{*}=\expectationvalue{\mathcal{M}}$.

So far, we assumed $\beta=1$. However, a different choice of $\beta$ can be made, which would affect to the transformation rule for the symmetries that filps monopole to anti-monopole. In fact, we can show that this is related to the identification rule between monopole ${\cal M}^{\dagger}$ and the VBS order parameters $v_{x}$ and $v_{y}$. For $\beta=1$, we get the relation Eq. (6); for example, $T_{x}:\,\mathcal{M}^{\dagger}\mapsto-i\mathcal{M}$ implies that the $T_{x}$-invariant monopole condensation corresponds to the condition that $\real\mathcal{M}+\imaginary\mathcal{M}=0$. However, if $\beta=\beta=e^{i\pi/4}$, we would get $T_{x}:\,\mathcal{M}^{\dagger}\mapsto-\mathcal{M}$, implying that the $T_{x}$-invariant condition is $\real\mathcal{M}=0$. In such a case, we can deduce that $\mathcal{M}^{\dagger}\sim v_{x}+iv_{y}$.

Once we fix the gauge choice and determine how the monopole transforms under the square lattice symmetries ($p4m$), how the monopole transforms under the Shastry-Sutherland lattice symmetries ($p4g$) can be deduced easily. This is because the Shastry-Sutherland lattice symmetries can be expressed in terms of the square lattice symmetries if we take the Shastry-Sutherland lattice in Fig. 5 is deformed from the smaller square lattice. The result is summarized in Tab. 6.

Appendix B Comparison with a DMRG simulation result of the $J_{1}$-$J_{2}$ model in the square lattice

In this section, we present our analysis on the square lattice with spin-1/2. Using the iDMRG simulation, we studied these models on an infinite cylinder with a circumference size up to $L=10$ lattice sites. Here, we focused on the correlation length spectra. The simulation is explicitly $U(1)_{z}$ symmetric, and we can plot the correlation spectra for each $U(1)_{z}$ quantum number, $S_{z}$. Since following models are $\mathrm{SO}(3)$ symmetric in the microscopic Hamiltonian, there must exist some degeneracy between different $S_{z}$ sectors, which can be interpreted as the spectrum for a higher spin.

First, let us consider the case where the square lattice symmetry is broken. The following model realizes the phase transition between Néel order and dimerized phase:

[TABLE]

where red and blue bonds are shown in Fig. 14(a). Here, the dimerized phase does not break any symmetry because the square lattice symmetry is already broken in the model. Therefore, the transition should be described by the Landau-Ginzburg theory. Indeed, it is known from the quantum Monte Carlo simulation Matsumoto et al. (2001) that the system realizes $\mathrm{O}(3)$ Wilson-Fisher critical point at $J_{1}/J_{1}^{\prime}=0.523$. In the iDMRG simulation, we also observed that the Néel order parameter develops at $J_{1}/J_{1}^{\prime}\sim 0.52$. At the transition, a single monopole event is not suppressed because different configurations for single monopole event cannot exactly cancel each other due to the absence of symmetries. As a result, the gauge fluctuation (VBS order parameter fluctuation) becomes confining, and the CP1 theory is no longer valid. Instead, the critical theory is described by the classical NLsM with $\mathrm{O}(3)$ Néel vector. The correlation spectra in Fig. 14(c) shows that the spin-triplet correlation length is the largest across the phase transition while the spin-singlet correlation length is much smaller than that. This behavior is consistent with the critical theory described by the classical NLsM.

Next, we study the $J_{1}$–$J_{2}$ Heisenberg model with a square lattice symmetry. The model is defined by the following Hamiltonian:

[TABLE]

where ${\bm{S}}_{i}$ is a spin-1/2 operator, $J_{1}$ is the nearest-neighbor AFM coupling, and $J_{2}$ is the next nearest-neighbor AFM coupling, see Fig. 14(b). When $J_{2}=0$ ($J_{1}=0$), the model is known to realize Néel ordered (conventional AFM stripe) phase. For the intermediate value of $J_{1}/J_{2}$, the system is frustrated and known to realize the disordered phase.

In accordance with the recent IPEPS study Haghshenas and Sheng (2018), we obtained the Néel, columnar VBS (cVBS), and conventional stripe phases as we increase $J_{2}/J_{1}$. However, in order to obtain the VBS order, we had to apply some bias (pinning field). Under the absence of the bias, the system looks totally symmetric, implying the existence of the symmetric superposition of symmetry broken states, namely a Cat state. The Cat state can be preferred over the symmetry broken phase if the circumference size is comparable to the length scale associated with the fluctuations, which is the size of the monopole in this case.

Before getting into the discussion of the correlation length spectra, we want to elaborate on some simulation details. In our iDMRG simulation, the iDMRG unit cell consists of two columns of lattices along the circumference, as otherwise translational symmetry broken phases (AFM order or VBS order) cannot develop. The price to pay is that $k_{x}=0$ and $k_{x}=\pi$ momenta become indistinguishable. (In our simulation, $k_{y}$ cannot be measured) However, at the critical point where the explicit symmetry breaking order has not developed yet, we can use a single column to distinguish $k_{x}=0$ and $k_{x}=\pi$ momenta. Indeed, in the iDMRG simulation of the single-column unit cell, we observe that the largest correlation length for $S=1$ spectra carries momentum $k_{x}=\pi$ in the single-column simulation, which is consistent with the momentum $(\pi,\pi)$ of the gapless magnon in Néel order. For $S=0$ case, we obtain that the lowest one carries $k_{x}=0$ while the second lowest one carries $k_{x}=\pi$, which runs almost parallel to the lowest one. They correspond to the $\mathbb{Z}_{4}$ VBS order parameter fluctuations (spin-singlet) at $(0,\pi)$ and $(\pi,0)$, but the degeneracy is lifted due to the iDMRG geometry which breaks the $C_{4}$-rotation lattice symmetry.

Surprisingly, in this model, we obtained the correlation length spectra that exactly agrees with the level crossing behavior of the excitation spectrum in the finite DMRG algorithm (Fig.2 in Ref. Wang and Sandvik, 2017). Our result is consistent with Ref. Zauner et al., 2015, which discussed the agreement between correlation length spectrum and local excitation spectrum in the DMRG simulation. Moreover, we want to comment on the argument in Ref. Wang and Sandvik, 2017. In this previous work, it was argued that the small region where $\xi_{S=1}>\xi_{S=0}>\xi_{S=2}$ corresponds to the gapless spin liquid phase Wang and Sandvik (2017). However, this reasoning is inconsistent with the iDMRG simulation result of the Shastry-Sutherland lattice because this region is clearly a symmetry broken phase with a non-zero Néel order parameter from the numerics. Thus, we can conjecture that such a region would shrink into a critical ‘point’ rather than remain as an extended phase of Dirac spin liquid both in the $J_{1}$-$J_{2}$ square lattice model and Shastry-Sutherland model. Indeed, if we perform a single-column iDMRG simulation for the $J_{1}$-$J_{2}$ model, the simulation does not converge well for the $J_{1}/J_{2}>2.0$, which means that the hypothesized gapless spin liquid phase is, in fact, more like the AFM phase where the double-column iDMRG unit cell is required.

Finally, we remark on the evidence that supports our discussion in Sec. IV. In the previous finite DMRG works on this model Gong et al. (2014), the plaquette VBS appeared instead of the columnar VBS. In fact, it was found in Ref. Haghshenas and Sheng (2018) that these two states have almost the same energy ($\Delta E/E<0.1\%$). This again implies that the dangerously irrelevant operator ${\cal M}^{4}$, which is responsible for the VBS ordering, has not flowed large enough to condense monopole to a certain direction. This can be supported by Fig. 9(a), where the correlation length for the spin-singlet operator is smaller than the correlation length of the spin-triplet operator throughout the whole intermediate regime between the Néel and conventional AFM stripe order.

Appendix C VBS-Phonon Coupling and Phonon Spectrum

The square pVBS order breaks the glide reflection symmetries $G_{x},G_{y}$ and the diagonal reflection symmetries $\sigma_{xy},\sigma_{x\bar{y}}$. Due to the lattice symmetry breaking, the pVBS order should induce lattice distortion as shown in Fig. 6. Therefore the pVBS fluctuation must couple to the lattice vibration mode, i.e. the phonon mode. Here we would like to determine the specific form of the pVBS-phonon coupling.

We will focus on the copper lattice in the following discussion. Although the lattice also contains other atoms and the phonon spectrum can be complicated, we choose to work on the symmetry level to demonstrate the universal consequences of pVBS fluctuation on the phonon spectrum without diving into the details. For this purpose, we first specifies the coordinate of copper sites in the each unit cell. As shown in Fig. 15(a), there are four copper sites in each unit cell. At equilibrium, they locate at

[TABLE]

where $0\leq\delta<1$ parameterize the deformation of the Shastry-Sutherland lattice from the square lattice. According to Ref. Miyahara and Ueda, 1999, the lattice constant is 8.995Å, the shortest Cu-Cu bond is 2.905Åand the second shortest Cu-Cu bond is 5.132Å. This implies that, in unit of the lattice constant, we have

[TABLE]

The optimal solution is $\delta\approx 0.544$. Using this deformation parameter, we can write down equilibrium positions of copper atoms through out the lattice.

Our strategy to figure out the pVBS-phonon coupling is to first investigate the pattern of lattice distortion induced by the pVBS order. Because the pVBS order does not enlarge the four-copper unit cell, its induced lattice distortion will also be identical among unit cells. Therefore the distortion can be described by four displacement vectors ${\bm{u}}_{A},{\bm{u}}_{B},{\bm{u}}_{C},{\bm{u}}_{D}$, translating each sublattice separately as

[TABLE]

The energy cost associated with the distortion can be modeled by summing up the bond energies

[TABLE]

where ${\bm{u}}_{i}$ dependence is implicit in ${\bm{r}}^{\prime}_{i}={\bm{r}}_{i}+{\bm{u}}_{i}$. The energy will increase whenever a bond is stretched or compressed. The shape of the potential in Eq. (28) captures this physics when the distortion ${\bm{u}}_{i}$ is small. The two stiffness coefficients $k_{1}$ and $k_{2}$ are expected to be different in general. Of course, in the realistic material, Sr, B, O atoms will all be involved and the bond energy model will be more complicated. However the toy model Eq. (28) respects all the symmetry property and provides the stiffness to the copper lattice, which can be used to analyze the pVBS induced lattice distortion. Finally we note that the energy model Eq. (28) written with respect to a single unit cell with periodic boundary condition (i.e. on a torus geometry), so for those bond across the unit cell, their bond length must be correctly treated by considering the periodic boundary condition.

Upon introducing the pVBS order, we will add an additional term to the energy model,

[TABLE]

where $p$ denotes the square plaquettes and $i\in p$ denotes the four corner sites around the plaquette $p$. ${\bm{R}}_{p}$ coordinates the plaquette center,

[TABLE]

The staggering factor $(-)^{p}$ is $+1$ for $X$-type plaquette (yellow) and $-1$ for $Y$-type plaquette (green) as shown in Fig. 15(a). $\imaginary\mathcal{M}=v_{y}-v_{x}$ denotes the square pVBS order parameter. The physical meaning of $E_{\text{VBS}}$ is that the pVBS order will contract one type of the square plaquette and expand the other type, exerting forces on copper atoms that points towards or away from the plaquette center.

Given the full energy model in Eq. (29), we can expand $E[{\bm{u}}_{i}]$ to the quadratic order of ${\bm{u}}_{i}$. The linear term will be proportional the pVBS order parameter $\imaginary\mathcal{M}$, as $\imaginary\mathcal{M}$ is the force that distort the lattice. The quadratic term determines how the lattice responses to the distortion force in the linear response regime. We found that independent of the choice of $\delta$ and $k_{2}$, the response is always given by

[TABLE]

Under the Fourier transformation to the momentum space

[TABLE]

the solution in Eq. (31) corresponds to

[TABLE]

This calculation indicates that the square pVBS order will lead to a lattice distortion that corresponds to a the simultaneous condensation of phonon modes $u_{x}$ at momentum $(\pi,0)$ and $u_{y}$ at momentum $(0,\pi)$. Given that $\imaginary\mathcal{M}=v_{y}-v_{x}$, we conclude that there must be a linear coupling between lattice displacement and the VBS order parameter in the form of

[TABLE]

in order to produce the linear response in Eq. (33). The coupling in Eq. (34) can be further justified by symmetry arguments. Tab. 7 shows the momentum quantum number and the symmetry transformation of the VBS order parameter ${\bm{v}}$ and finite-momentum phonon mode ${\bm{u}}$. One can see ${\bm{v}}$ and ${\bm{u}}$ have identical symmetry properties and hence a linear coupling as in Eq. (34) is allowed.

Given the VBS-phonon coupling, we can investigate the effect of low-energy VBS fluctuation on the phonon spectrum near the DQCP. We first write down the field theory action describing both degrees of freedom,

[TABLE]

where $q=(\omega,{\bm{q}})$ represents the energy-momentum vector. $\Omega_{{\bm{q}}}$ describes the phonon dispersion relation. $G_{\text{VBS}}$ is the correlation function of the VBS critical fluctuation, whose low-energy behavior is given by

[TABLE]

where $\eta$ is the anomalous exponent of the $\mathrm{O}(4)$ vector at the $\mathrm{O}(4)$ DQCP. Based on the previous numerical measurementsQin et al. (2017); Zhang et al. (2017), $\eta$ is estimated to be $\eta=0.13\sim 0.3$. ${\bm{Q}}=(\pi,0)$ or $(0,\pi)$ denotes the momentum point where the VBS fluctuation gets soften. The VBS-phonon coupling $\kappa_{\bm{q}}$ is expected to be momentum dependent, because the VBS order parameter only couples to the high-energy phonon around $X$ and $Y$ points but not the acoustic phonon around $\Gamma$ point. By the acoustic phonon around $\Gamma$ point, we mean the low energy part of the phonon, i.e. the ‘segment’ of the acoustic branch around the gapless point, which usually appears in the field theory description of phonons. Given these setup, we can integrate out the VBS fluctuation and obtain the dressed propagator of phonon,

[TABLE]

Then the phonon spectral function can be obtained from

[TABLE]

The phonon dispersion $\Omega_{\bm{q}}$ is unknown to us, as we did not have the full model of the lattice vibration. For demonstration purpose, we can use the following toy model

[TABLE]

which captures the gapless acoustic phonon at the $\Gamma$ point and gapped phonons at $X$ and $Y$ points (see Fig. 1(b)). We also take the anomalous exponent $\eta=0.13$ and use $\kappa_{\bm{q}}=\kappa_{0}\Omega_{\bm{q}}$ with $\kappa_{0}=0.05$ so that

[TABLE]

With these, we show the phonon spectrum calculated from Eq. (38) in Fig. 15(b). The prominent feature is a V-shape continuum at the $X$ point (as well as the $Y$ point) in the Brillouin zone. This continuum in the phonon spectrum represents the critical fluctuation of the VBS order parameter at the DQCP. Although the spectral weight is expected to be weak, since the $X$ point is an extinction point, it is still feasible to collect spectral signals of this continuum. In particular, the frequency dependence of the spectral weight at the $X$ point is predicted to follow

[TABLE]

which can be checked experimentally. It will be meaningful to compare the measured $\eta$ with large-scale quantum Monte Carlo simulation result.

Appendix D Detailed Numerical Data

In this appendix, we discuss the evolution of the iDMRG simulations results as we increase the bond dimension $\chi$ for system sizes $L=8$ and $L=10$. Since the accuracy of the iDMRG simulation is determined by the bond dimension, a reliable analysis requires one to examine the results as a function of the bond dimension. Here, the iDMRG simulation results of the Shastry-Sutherland lattice model in Eq.1) at $L=10$ for a range of bond dimensions are presented in Fig. 16. Although the the truncation error $\epsilon_{\textrm{trun}}$ is very large ($>10^{-4}$) at the low bond dimensions, as we increase the bond dimensison upto $\chi=4000$, $\epsilon_{\textrm{trun}}$ goes below $10^{-5}$ and the iDMRG results becomes sensible.

Fig. 16(b) shows the first derivative of energy, whose change of the slope correspsonds to the transition between the pVBS and Néel ordererd phases. As the bond dimension increases, the transition point shifts leftward, implying that the parameter regime for the pVBS phase shrinks. The behavior aligns with the intuition that the gapped pVBS phase would be favored over the the gapless Néel orderd phase for a low entanglement MPS state. Accordingly, the peak of the spin-singlet correlation length which coincides with the phase transition point also shifts leftward, presented in Fig. 16(d). At the same time, the peak of the spin-singlet correlation length becomes larger and more pronounced as the bond dimension increases, signaling the continuous or weakly first order phase transition.

Although the peak location of the spin-singlet correlation length changes with the bond dimension, a further indication of the phases boundary can be obtained from the order parameter plotted versus the truncation error Hubig et al. (2018); Szasz et al. (2018). In principle, the ground state is fully symmetric and local order parameter cannot be non-zero. However, in the iDMRG simulation, the numerical process favors a minimally entangled state, giving rise to a non-zero local order parameter in a spontaneous symmetry breaking phase at finite bond dimensions. This can even happen when the system is outside but close to the spontaneous symmetry breaking phase, meaning that one needs to plot the order parameter as a function of the truncation error in order to see whether a non-zero order parameter is truly physical. In Fig. 16(f), we plotted the pVBS order parameter defined in Eq. 5 as a function of the truncation error, and extrapolated them. We observe that the extrapolated order parameter disappears at $J_{1}=0.76$, which agrees well with the peak location $J_{1}=0.762$ of the spin-singlet correlation length at $\chi=4000$ with $L=10$. As mentioned earlier, this extrapolation method would be benefited a lot if a wider range of bond dimensions is available. However, at $\chi=4000$, each data point already takes about 50 hours of simulations times for 12 multi-threads with 80GB RAM, and both the time and RAM scale roughly as $\chi^{2}$. Therefore, a significantly higher bond dimension is currently inaccessible at our capacity.

Finally, we remark the scaling of the result for different system sizes. In Fig. 17, we plotted the iDMRG simulation results at $L=8$ for a range of MPS bond dimension as in Fig. 16. Note that the truncation error is an order of magnitude smaller than the results at $L=10$. Moreover, the value of $J_{1}$ where the pVBS order parameter disappears is much smaller for the smaller system size $L$. However, this does not imply that the Néel ordered phase develops for $J_{1}>0.726$, as we can see from the correlation length plot in Fig. 18. In Fig. 18(a), while the peak of the spin-singlet correlation length implies that the peak corresponds to the phase transition point for the pVBS phase, the spin-triplet correlation length remains almost constant, which implies that the Néel order does not develop, since the Néel order would give rise to the increase of the spin-triplet correlation length originated from the gapless magnon excitation. On the other hand, at $L=10$, we observe that the peak of the spin-singlet correlation length is immediately followed by the rise of the spin-triplet correlation length which signals the onset of the Néel ordered phase. $L=6$ behavior is similar to that of $L=8$, and the signature of the Néel ordered phase, such as the staggering magnetization or the rise of the spin-triplet correlation lengh, is not observed near the pVBS critical point. This is related to the discussion in the main text. For a quasi one-dimensional system with a finite-size circumference size, the spontaneous symmetry breaking of the continuous group would be suppressed due to the disordering effect. As a result, the disappearance of the pVBS ordering is not immediately followed by the Néel ordering for a finite circumference system.