Valence-bond solids, vestigial order, and emergent SO(5) symmetry in a two-dimensional quantum magnet

We introduce a quantum spin-1/2 model with many-body correlated Heisenberg-type interactions on the 2D square lattice, designed to host a plaquette valence-bond solid (PVBS) ground state breaking $\mathbb{Z}_4$ symmetry. We carry out a detailed quantum Monte Carlo study of the quantum phase transition between the antiferromagnetic (AFM) and PVBS states. We find a first-order transition, in contrast to previously studied continuous 'deconfined' transitions between the AFM and columnar valence-bond solids. We show that the coexistence state at the AFM--PVBS transition is unexpectedly associated with SO(5) symmetry, which may indicate that the transition is connected to a deconfined critical point. We also discuss the first-order transition in the context of a recent proposal of spinons with fracton properties in the PVBS state, concluding that the fracton scenario is unlikely. Further, we discover a novel type of eight-fold degenerate VBS phase that breaks the remaining $\mathbb{Z}_2$ symmetry of the PVBS phase. The PVBS phase can then be regarded as an intermediate 'vestigial' phase, and we develop a general graph-theoretic approach to describe the two-stage discrete symmetry breaking. At finite temperature, we observe fluctuation-induced first-order transitions, which are hallmarks of vestigial phase transitions. We also mention possible connections to the SO(5) theory of high-T$_{\rm c}$ superconductivity.


I. INTRODUCTION
Competing interactions in a quantum antiferromagnet with spin-rotationally invariant interactions can lead to the destruction of the conventional Néel antiferromagnetic (AFM) O(3) symmetry-broken ground state in two or higher dimensions.The most well studied transition into a gapped quantum paramagnet is in dimerized spin-1/2 systems (two spins per unit cell), where the ground state is non-degenerate and can be approximated as a product of singlets on the dimers.This quantum phase transition has an experimental realization in TlCuCl 3 under pressure [1].For a uniform system with one spin per unit cell, a uniform non-degenerate paramagnetic ground state is possible only under special conditions [2,3].With a half-integer spin per unit cell (e.g., a single spin-1/2) a more complex state must necessarily obtain which has either non-local entanglement and topological order (a spin liquid [4]) or strong local entanglement leading to a degenerate symmetry-breaking pattern of singlets (a valence-bond solid, VBS, sometimes also called valencebond crystal) [5][6][7][8][9][10][11].
While experimental searches for spin liquids have been a dominant theme in quantum materials science for more than a decade [12], comparatively less efforts have been devoted to quantum magnets with VBS states.A dimer VBS in the quasi-one-dimensional material GeCuO 3 was exhaustively studied in the 1990s [13].More recently, signs of a VBS with singlets forming on four-spin plaquettes were detected in the quasi-two-dimensional system SrCu 2 (BO 3 ) 2 under high pressure [14].This material also hosts an adjacent AFM phase at still higher pressures [15], thus for the first time opening prospects for detailed experimental studies of a direct transition between AFM and VBS states in two dimensions.
The AFM state has an obvious classical counterpart, and its low-energy properties, including its quantum phase transition into the dimer paramagnet, can be understood within essentially classical field theories in one higher dimension (imaginary time, to account for quantum fluctuations) [16,17].In contrast, spin liquids and VBSs are exotic from the classical standpoint, and lowenergy descriptions of them and their quantum phase transitions require more drastic deviations from the conventional field theories of statistical physics [18,19].The mathematical complexity of the quantum-field theories used to described exotic states of quantum magnets is formidable, and direct computational studies of the lowenergy properties of lattice models are indispensable for testing and guiding analytical approaches, and can also open new research directions [20].
In this article we introduce a two-dimensional (2D) quantum spin model in which the AFM ground state successively transitions into two different types of VBS states; a four-fold degenerate plaquette VBS (PVBS) followed by an eight-fold degenerate state we name the alternating VBS (AVBS).These states are schematically depicted in Fig. 1, along with the more commonly studied columnar VBS (CVBS) state.In the idealized PVBS state illustrated in Fig. 1(b), a plaquette is a resonating pair of valence bonds, with equal amplitude for horizontal and vertical orientation.Such a singlet state also goes under other names, e.g., plaquette-singlet solid or valence-plaquette solid.Here we adopt the name VBS as a generic term for non-magnetic states breaking lattice symmetries, forming ordered unit cells with either static or resonating valence bond descriptions.In general, quantum and/or thermal fluctuations perturb the completely ideal singlet product states in Fig. 1, but the depicted patterns survive in actual ground states by having finite valued order parameters.The pattern can also be explicitly visualized by observing the singlet-density maps, as we will do later.
Using quantum Monte Carlo (QMC) simulations of the J-Q 6 model illustrated in Fig. 2, we find an unusual firstorder AFM-PVBS quantum phase transition with emergent SO(5) symmetry of the five-component combined AFM (three components) and PVBS (two components) order parameters.To analyze the second transition, a continuous PVBS-AVBS transition, we introduce a new unified generic graph-theoretic description of multi-stage discrete symmetry breaking transitions, which we use here to analyze the PVBS and AVBS order parameters and their possible symmetry breaking paths.
The AFM-PVBS transition is an interesting analogy to the AFM-superconducting transition within the SO (5) FIG. 2. The J-Q6 model in which PVBS and AVBS ground states (Fig. 1) are realized.The two types of interactions between the spins on the square lattice are defined using two-spin singlet projectors, here depicted with ovals.The standard Heisenberg exchange of strength J between nearest neighbor spins is equivalent to projectors on the horizontal and vertical links, while the interaction of strength Q introduced here is expressed with products of six projectors, arranged in four different ways.The site labels are used in the formal definition of the Hamiltonian.
theory of the cuprates [21,22].Overall our results connect to several of the major concepts currently under debate from the field-theory perspective; deconfined quantum critical (DQC) points [23], spinons with fracton properties [24], and vestigial phase transitions [25].In the remainder of this introductory section we provide further background on these notions and summarize our aims and findings.

A. Valence-bond solid states and deconfined quantum criticality
VBSs are often regarded as less exotic and interesting than spin liquids, because they have local order parameters.However, the objects that form the long-range order in a VBS-the singlets-are not the elementary microscopic degrees of freedom of the model but emergent entangled composites with no direct counterparts in conventional classical spin models.It was realized in 1980s and 1990s that continuum field-theory descriptions of one-dimensional (1D) [26][27][28] and 2D [8,29,30] VBS states require careful consideration of topological defects and their quantum interference, which is still an active field of research [24,[31][32][33].While the early field-theory studies already indicated that 2D AFM-VBS transitions may be unusual, it was only after intriguing numerical results were found [34,35] that it was recognized that direct quantum phase transitions between AFM and VBS ground states may fall outside the standard Landau-Ginzburg-Wilson (LGW) paradigm [36,37] (as was already known to be the case with many transitions in 1D systems).Concurrently, some other exotic non-LGW 2D transitions related to VBSs were also proposed in quantum-dimer systems [38].
Within the LGW framework, direct transitions between two ordered phases with unrelated symmetries are generically first-order.However, the two order parameters at the AFM-VBS transition are not independently fluctuating but represent different manifestations of the confinement of more fundamental objects-spinonswhich deconfine at the phase transition (the DQC point).The expected first-order transition may then be supplanted by a generically continuous transition.
Although the DQC scenario predicts fascinating connections between the AFM and VBS states, it relies on theoretical assumptions such as the conservation of skyrmion numbers at the transition point [35][36][37].Some of the assumptions are hard to justify rigorously, and, thus, unbiased computational research plays a crucial role in this developing field.Numerical confirmations of the DQC phenomenon face several challenges, however.To begin with, it is difficult to distinguish a continuous transition from a weak first-order transition in general.In order to deal with this issue in a reliable manner, one needs to utilize a method which can reach large lattice sizes while still being free from systematic errors.QMC methods serve such purposes, but in many cases they suffer from the sign problem [39,40] and are not useful.
The conditions for DQC transitions are at first sight most naturally realized in geometrically frustrated SU(2) symmetric quantum spin systems, which is one of the model classes for which the QMC sign-problem is most severe [41,42].The lattices accessible with exact numerical diagonalization methods are too small to even allow definite characterization of the non-magnetic states in these systems [10,11,43].While in principle the densitymatrix renormalization (DMRG) method or variational states based on tensor networks can be used, and some results for potential DQC systems have been reported [44][45][46][47], these techniques have not yet reached the length scales and reliability of modern QMC techniques.
The invention of the J-Q class of Hamiltonians [48,49] has been a solution to the challenge of QMC simulations of the DQC transition in quantum spin systems.In these models the Heisenberg interaction J is supplemented by multi-spin correlated singlet projectors Q, which open possibilities of designing sign-problem free Hamiltonians with the desired phase diagrams containing AFM and VBS phases.A large number of QMC studies with different variants of 2D J-Q models have been reported [48][49][50][51][52][53][54][55][56][57][58][59][60][61][62][63][64][65].Concurrently, various discrete versions of the proposed DQC action have also been studied numerically [66][67][68], and three-dimensional classical loop [69,70] and dimer [71,72] models with arguably the same low-energy fixed points have also been studied extensively.Although the observed quantum phase transitions in many of these models appear to be continuous, there are puzzling scaling anomalies that have been interpreted either as signs of an eventual weakly first-order transition (with discontinuities presumably developing on larger lattices than currently reachable) [23,68,73,74], or as manifestations of new physics associated with the DQC scenarionovel finite-size and finite-temperature scaling behaviors related to the presence of two divergent length scales [59].
While the ultimate continuous or first-order nature of the AFM-VBS transition is a question of fundamental interest, it should be noted that the DQC scenario does not stand or fall with it-what matters is whether spinon deconfinement takes place on large length scales, and this is now beyond doubt.The critical or near-critical state is at the very least a good approximation to a particular case of a gapless Dirac spin liquid [60,63,64], as expected at the DQC transition.Experimental realizations of spinon deconfinement [75,76] and DQC transitions are plausible [14,15,77], though no conclusive observations have been reported as of yet.
The broader interest in the DQC point and related phenomena, such as emergent symmetries and anomalous scaling, lies in the many unresolved puzzles in quantum matter where "beyond-LGW" physics may be at play.For instance, some of the unusual metallic properties of the high-T c cuprates may arise from doping the DQC point [78,79].Being accessible to large-scale QMC simulations, the DQC phenomenon and other exotic aspects of VBS states and transitions offer unique opportunities to study important cases of beyond-LGW physics in detail.With a suitable sign-free "designer Hamiltonian" [20], definite numerical results can be obtained and compared with approximate and often speculative predictions of field theories, and unique insights can be gleaned from the simulations in their own right.

B. J-Q6 model and main findings
Stimulated by the success of designer Hamiltonians suitable for QMC studies of the DQC phase transitions and VBS physics, many other sign-free quantum spin models have recently been constructed for studies of a wide range of exotic quantum phases and quantum phase transitions, e.g., symmetry-enhanced first-order transitions [65], Z 2 spin liquids [80], and Haldane nematics [81].In this article we introduce a new type of J-Q Hamiltonian that exhibits a host of fascinating phenomena within a single phase diagram.Our variant of the J-Q 6 model (where the subscript denotes the number of singlet projectors) contains 12-spin interactions (see Fig. 2) and may appear contrived from the standpoint of experimental realizations.However, in the spirit of designer Hamiltonians [20], it opens access to studies of quantum states and quantum phase transitions which can likely be realizable also with other microscopic interactions (e.g., frustrated exchange interactions) and whose quantum field-theory descriptions are of great current interest.
Being sign free, the J-Q 6 model illustrated in Fig. 2 (and define in detail in Sec.III) is amenable to QMC simulations on lattices with thousands of spins.We will construct the ground-state phase diagram as a function of the ratio Q/J of the twelve-spin and Heisenberg couplings.On increasing Q/J, we first find a direct quantum phase transition between the standard Heisenberg AFM ground state and a four-fold degenerate PVBS state [Fig.1(b)].When further increasing the control parameter, the Z 4 symmetry breaking PVBS state undergoes a Z 2 -breaking continuous quantum phase transition into another state, the eight-fold degenerate AVBS state [Fig.1(c)].The initial aim of our study was to design a model hosting an AFM-PVBS transition as a possible realization of a DQC transition.The J-Q 6 model was the first and so far only successful sign-free Hamiltonian realizing a four-fold degenerate PVBS state, and in the course of studying it we also discovered the unexpected AVBS state.
Previous QMC studies of DQC transitions have focused on the CVBS state [1(a)], but the likewise fourfold degenerate PVBS is also a possible DQC candidate on equal footing with the CVBS state according to theory [18,36].However, the AFM-PVBS transition found here is clearly first-order though it exhibits the emergent SO(5) symmetry which was proposed in one variant of the DQC theory [82].Emergent symmetry would not in general be expected at a first-order transition [70,83], but recently other examples have been found [65,84,85] where the coexistence state appears to be described by a vector or pseudo-vector combining all the components of the two different order parameters and transforming under a spherical symmetry.In the case at hand here the combined order parameter comprises three AFM components and two VBS components.
The emergent SO (5) symmetry may indicate the proximity of a continuous DQC transition with this symmetry.We will also discuss possible relevance of the first-order transition to a recently proposed alternative fracton theory [24] of the PVBS state and AFM-PVBS transition.Furthermore, we demonstrate that the twostage breaking of the symmetries of the PVBS and AVBS phases fits into the framework of "vestigial" phase transitions [25,86], where our case is the first example where both phases break discrete symmetries-the previous cases involved the breaking of a discrete symmetry followed by the breaking of a continuous symmetry.We describe the double-discrete ground-state vestigial transition using a six-dimensional order parameter and also by a novel graph-theoretic approach.The first-order nature of the transitions at finite (non-zero) temperature from a paramagnet into either the PVBS or AVBS state lends further support to the vestigial phase scenario.

C. Article outline
The remainder of the article is structured as follows: In Sec.II we provide further background on VBS states and numerical studies of the AFM-VBS transition, setting the stage for our new developments related to PVBS ordering, emergent symmetries, and vestigial phase transitions.In Sec.III we define our J-Q 6 model in detail and briefly describe the QMC method we use to study it.In Sec.IV we discuss results for the AFM-PVBS groundstate transition and contrast it with the often studied AFM-CVBS transition.In Sec.V we present results for the second phase transition into the AVBS state.In Sec.VI we construct an order parameter that captures both the PVBS and AVBS phases, and also present our graph-theoretic approach (the "order graph") for classifying the two-stage symmetry breaking.We present results for finite temperature in Sec.VII and explain them based on the graph approach and the scenario of vestigial phase transitions.Finally, in Sec.VIII we discuss implications of our results to existing theories, including the fracton scenario and the analogy of the AFM-PVBS transition to the cuprate SO(5) theory [21].In Appendix A we provide a detailed discussion on the subtle issue of exactly what symmetries are broken in the CVBS and PVBS phases.In Appendix B we further discuss the motivations behind the concept of the order graph, and also describe its symmetry properties in more detail.

II. VALENCE-BOND SOLIDS AND EMERGENT SYMMETRIES
The square-lattice CVBS depicted in Fig. 1(a) is the most well-studied non-magnetic state in the context of the DQC transition.In addition to pointing to a continuous quantum phase transition between the CVBS and AFM states, numerical studies of J-Q and related models have revealed that the fluctuations among the four degenerate dimer patterns develop emergent U(1) symmetry as the critical point is approached [48-50, 70, 87].This emergent symmetry confirms an important aspect of the DQC scenario, where the ordered CVBS patterns can be assigned angles φ = nπ/2, n = 0, 1, 2, 3, and tunneling between these discrete angles corresponds to n traversing a range of continuous values (alternatively, the order parameter is a complex scalar) in an effective potential ∝ cos(4φ) [88].The emergent U(1) symmetry corresponds to a flat distribution of the coarse-grained CVBS angle φ at the critical point.
The proposal of emergent U(1) symmetry, which is also directly related to the conjectured absence of topological defects at the DQC transition [36,37], was partially motivated by an analogy between the VBS order parameter and the magnetization vector of a classical 3D XY model with a four-fold symmetric (q = 4) potential ∝ cos(qθ i ) for the microscopic spin angles θ i .In these clock models, when q ≥ 4 the cosine perturbation of the U(1) symmetric XY model is "dangerously irrelevant", i.e., it reduces the symmetry in the ordered state but not at the critical point (in the thermodynamic limit when the order parameter is coarse-grained over large regions), and the universality class of the phase transition remains to be that of the 3D XY model [89][90][91][92][93][94][95].In analogy, at the AFM-VBS transition the effective low-energy interactions responsible for locking the dimerization to one of the four static patterns may be dangerously irrelevant, so that the coarse-grained VBS order parameter can take any angle between 0 and 2π when the critical point is approached from the VBS side.
If there indeed is emergent U(1) symmetry, then the field theory describing the AFM-VBS critical point does not need to include the ingredients causing the four-fold degeneracy-quadrupled monopole instanton events in the path integral [8,30].In the originally proposed DQC field theory, the CP 1 model [36,37], this aspect can be taken into account by not including any topological defects in the U(1) gauge field corresponding to the continuously fluctuating VBS order parameter; hence the proposal that the transition is described by the non-compact CP 1 model.Here it should be noted that the analogy with the clock model is not precise, because the critical points are different, and it has not been demonstrated within the field theory that the quadrupled monopoles really are dangerously irrelevant.
Numerically, emergent U(1) symmetry has been established up to the largest length scales reachable in simulations of the J-Q model [48][49][50]87] and in related classical 3D loop models [69].Moreover, there are also indications of an SO(5) symmetry of the combined O(3) AFM and emergent U(1) VBS order parameters [60,69].While the original DQC scenario did not involve any such symmetry, only O(3)×U(1), and allowed for different exponents η A and η V of the critical AFM and VBS correlation functions, a later proposal was explicitly formulated with η A = η V .In this alternative (perhaps dual) theory the three components of the AFM order parameter and the two VBS components are treated on equal footing as a five-component vector transforming as SO(5) [82].This treatment is possible only with SU(2) spins, and the SU(N) generalizations of the CP N −1 theory do not appear to allow such a higher symmetry.Numerically, η A ≈ η V has been observed for SU (2) spins [48,55,87] and in the loop model [69], while generalizations of the J-Q model and other related models to SU(N) spins show clearly η A = η V , with the values of both exponents in remarkably good agreement with 1/N expansions of the SU(N ) theory for large N [48,55,96].Whether or not the observed SO(5) symmetry for SU(2) spins is exact or only approximate (i.e., breaking down on some large length scale) is still an open question, especially in light of the fact that the numerical values of the exponents η A and η V of the lattice models do not satisfy a bound obtained from conformal bootstrap calculations with SO(5) symmetry [97].
Although the DQC scenario has been probed in detail with models hosting CVBS ground states, the PVBS state depicted in Fig. 1(b) still remains to be studied thoroughly, mainly due to the lack of microscopic models realizing it without QMC sign problem.Models in which a PVBS has been proposed include the Heisenberg model with first-and third-neighbor AFM interactions [98], which has a severe QMC sign problem that has prohibited detailed studies of the putative AFM-PVBS transition.Resonating VBS states which spontaneously breaks the lattice symmetry similarly to the square-lattice PVBS state have also been studied on the honeycomb lattice, with frustrated Heisenberg Hamiltonians with sign problems [100,101] and J-Q models without [32,56] sign-problem.While it is apparently easier to construct sign-free J-Q type Hamiltonians with PVBS ground states on the honeycomb lattice, the symmetry broken in that case is Z 3 instead of Z 4 on the square lattice.The smaller number of degenerate states has an additional complication in whether it allows for emergent U(1) symmetry or not [61].We will only discuss the square lattice here.
The DQC scenario does not distinguish in any crucial manner between the CVBS and PVBS states, as they both are four-fold degenerate and subject to the same mechanism of emergent U(1) symmetry-essentially the difference is in the sign of the effective cosine potential experienced by the coarse-grained VBS order parameter.Contrary to the expectation of a DQC transition, there is a recent suggestion that spinons in the PVBS states will be immobile fractons [24], thus prohibiting the deconfinement mechanism underlying the emergent U(1) symmetry and the continuous quantum phase transition into the AFM state.In the simulations of the J-Q 6 model presented here, we indeed find a first-order AFM-PVBS transition.Nevertheless, we will argue that the absence of spinon deconfinement is not necessarily confirming the fracton picture, and that the fracton mechanism is not universal for PVBS states and may only be realized under very particular conditions.
Moving deeper into the PVBS phase, the J-Q 6 model exhibits another phase transition out of the PVBS phase into the eight-fold degenerate AVBS state illustrated in Fig. 1(c).This previously not anticipated transition breaks an additional Z 2 symmetry, and, thus the fourfold degenerate PVBS phase can be regarded as an intermediate phase that only partially breaks the D 4 symmetry of the square lattice.This two-stage symmetry breaking fits into the scheme of "vestigial" transitions [25,86], although previous examples of such multi-stage transitions have involved the breaking of a discrete (normally Z 2 ) symmetry followed by the breaking of a continuous symmetry, in contrast to both transitions breaking discrete symmetries in our case.We will introduce a general systematic way of describing such multi-stage discrete symmetry breaking using a graph-theoretic approach.
In addition to serving as important testing grounds for the theory of quantum magnetism and quantum phase transitions, PVBS states are important also considering potential realization of DQC transitions in real materials.So far, the most promising candidate is the quasi-2D quantum magnet SrCu 2 (BO 3 ) 2 under high pressure [14,15], which appears to realize a certain type of PVBS state (though there are also views opposing this notion [102]).The 2D magnetic interactions in SrCu 2 (BO 3 ) 2 realize the Shastry-Sutherland model [103], and the ratio of the inter-and intra-dimer coupling constants change with pressure in such a way that the three phases of the model are realized within the accessible pressure range; first a non-degenerate dimer singlet state, then a PVBS followed by an AFM state, with the latter further stabilized by inter-layer couplings [15].However, the Shastry-Sutherland PVBS state is only two-fold degenerate owing to the structure of the lattice, while the original DQC scenario applies to a four-fold degenerate PVBS such as the one illustrated in Fig. 1(b).Nevertheless, there are theoretical proposals suggesting that also a two-fold degenerate VBS state may give rise to a DQC point or DQC-like physics [77,104].
Beyond the sign-problematic (in the most interesting parameter regimes) Shastry-Sutherland model and extensions of it [99], sign-free 2D J-Q models with two-fold degenerate PVBS ground states can also be designed.In a "checker-board" J-Q (CBJQ) model, an unusual firstorder transition between the AFM state and a Z 2 breaking PVBS ground state was found [65].Surprisingly, an emergent O(4) symmetry of the combined O(3) AFM and scalar PVBS order parameters in the coexistence state of the system was observed, despite the clear presence of discontinuities at the transition.While until recently emergent symmetries had only been expected at continuous transitions [70,83], subsequently a first-order transition with enhanced symmetry was also found in a Z 2 deformed classical 3D loop model [84].It is currently not clear whether the enhanced symmetry is only manifested up to some large length scale or truly asymptotically exact, but a theory in a different context of boson-fermion supersymmetry does suggest that first-order transitions can be associated with exact emergent symmetry [85].Here we will demonstrate SO(5) symmetry at the firstorder AFM-PVBS transition of the J-Q 6 model, thus providing yet another example of the emergence of enhanced symmetries under unexpected conditions.

III. MODEL AND METHOD
The J-Q 6 model we consider here is defined with S = 1/2 spins on a 2D square lattice.As illustrated in Fig. 2, the model combines the standard Heisenberg exchange of strength J with a 12-body interaction of strength Q in the Hamiltonian were, P (i, j) denotes a bond operator; a singlet projector on the pair of sites i and j, Thus, a single J-term −JP (i, j) is equivalent to an AFM Heisenberg coupling, up to an additive constant.The sum over i, j runs over all nearest-neighbor bonds.The Q 6 terms act on twelve spins specified by the set H of groups of twelve sites i 1 , . . ., i 6 , j 1 , . . ., j 6 , with their relative positions and pairings into bonds with singlet pro-jectors illustrated in Fig. 2. In order to preserve all symmetries of the square lattice, H includes all possible translations of the four bond arrangements shown.Note that two of the four patterns correspond to a π/2 rotation of the other two, preserving the four-fold rotational symmetry of the lattice.We use L × L(= N ) square lattices with periodic boundary condition, yielding 2N J-terms and 4N Q-terms.In some cases we will also consider open boundary conditions, in which case we remove all operators acting on sites beyond the L × L edge.Previously studied J-Q models have Q-terms with two or three bond operators and were occasionally referred to as J-Q 2 or J-Q 3 models, respectively [48].Different J-Q models are further distinguished by the relative arrangement of the bond operators of the Q terms, with columnar and stair-case arrangement considered in the past-the former leading to DQC transitions, as discussed in the previous section, and the latter type of Q 3 interaction causing a strongly first-order transition between the AFM state and a staggered VBS [53] (and a stair-case Q 2 interaction does not destroy the AFM order).The combination of vertical and horizontal singlet projectors in the new Q 6 terms introduced here is the feature that promotes the formation of a PVBS instead of the CVBS obtaining if all projectors are stacked in columns.For comparison, we will also present some results for the case of a columnar Q 6 interaction, in which case we find a strongly first-order transition into a CVBS phase.We also note that we have not been able to realize a PVBS state with less than six bond operators.
For all calculations reported here we used the stochastic series expansion (SSE) method [105], a QMC method without discretization errors that incorporates efficient loop updates [106,107].In order to search the entire sign-free ground state phase diagram (positive J and Q), we vary Q ∈ [0, 1] in Eq. (1) while keeping J +Q = 1 fixed to serve as the unit of energy in the simulations.Unless otherwise noted, the inverse temperature is set to β = L, which is sufficient for studying ground state ordering and scaling properties of quantum phase transitions with dynamic critical exponent z = 1, which is the expected z at both the DQC point and a quantum Ising critical point (which corresponds to the symmetry-breaking at the PVBS-AVBS transition).The choice of β(L) is also suitable for detecting a first-order phase transition.

IV. DIRECT AFM-TO-PVBS TRANSITION
The J-Q 6 model exhibits a direct quantum phase transition from the standard AFM phase to the PVBS phase at a point Q c ≈ 0.2758, making this transition a candidate for the DQC mechanism.However, we here instead find a first-order transition with clearly finite coexisting order parameters at the transition point.This is in sharp contrast to previously studied variants of the J-Q model, which host apparently continuous AFM-CVBS transitions, or possibly very weak first-order transitions   with coexisting order parameters too small to be detected currently.Here we will first present our numerical simulation results in detail and leave discussion of implication of our result to the DQC and fracton theories mostly to Sec.VIII.

A. Order parameter definitions
Let us first define the two components Π 0 and Π 1 of the PVBS order parameter Π = (Π 0 , Π 1 ) as where a = 0, 1 is a sublattice label corresponding to a checker-board pattern of the lattice coordinates (x, y), i.e., a ≡ x + y mod 2. P z x,y is a projection operator to zero z-magnetization of the four spins on the plaquette with its low-left corner at (x, y), i.e.,
(4) Eq. ( 3) is just one out of many possible order parameters capable of detecting the PVBS state.Ideally, one might prefer a spin-rotation invariant definition, e.g., some direct measure of the singlet density.Evaluating such offdiagonal operations involving more than two spins is very time consuming, however [108].Though the condition for P z x,y = 1 is a necessary but not sufficient condition for the four sites to form a singlet, the order parameter Π still detects a modulation of the mean singlet density.
It is useful to compare the plaquette order parameter with the frequently used dimer order parameter D = (D x , D y ), where, to conform with the plaquette order parameter in Eq. (3), we also use a diagonal definition (instead of the rotationally invariant definition that can also be evaluated efficiently [108]): For a CVBS or PVBS ordered state, the plaquette order parameter (Π 0 , Π 1 ) essentially behaves as a π/2 rotated columnar order parameter (D x , D y ), and both squared order parameters are non-vanishing.This is a consequence of the fact that the CVBS and PVBS states both break the D 4 lattice symmetry into Z 2 .To be more precise, when we regard symmetry transformations of the model that do not change the order parameter, they will naturally form the eight-component dihedral group D 4 , which breaks into Z 2 in both the CVBS phase and the PVBS phase.The two remaining Z 2 symmetries are isomorphic via an automorphism of the D 4 group.This point is explained in detail in Appendix A. For the discussion in this section, it suffices to recognize that Π and D are both valid order parameters for detecting CVBS and PVBS order, and these orders can be distinguished by examining the 2D distribution of either of the order parameters; here we will analyze both of them.

B. First-order transition
Examples of SSE results used in order to locate the AFM-PVBS transition point are shown in Fig. 3. Fig. 3(a) shows the squared order parameters for three system sizes versus Q.It is apparent that the point where the AFM order dies out is also where the PVBS order emerges, implying a direct transition.To precisely analyze the transition, we show in Fig. 3(b) the Binder cumulants U m and U Π of the order parameters, defined such that U X → 1 with increasing system size if there is long-range order of the X kind (X = m or X = Π) and U X → 0 otherwise.The Binder cumulant of the AFM order parameter is given by where the factors are chosen for a single component of the three-component AFM order parameter, For the two-component plaquette order we have where Π = (Π 0 , Π 1 ), with the components defined in Eq. ( 3).The crossing value and U X (Q, 2L) as the control parameter Q is varied (for any order parameter X), is known to have much smaller finite-size drifts compared to other quantities defining finite-size transition points.L-2L crossing values Q * (L) obtained from data such as those in Fig. 3(b) are graphed versus the inverse system size in Fig. 3(c).The estimated transition point Q c = Q * (L → ∞) in the thermodynamic limit agrees well between different order parameters.Another quantity, which serves even better as the Q c estimator for this system, is the crossing between the two different Binder cumulants computed with the same lattice size, thus providing a larger number of L points with the available data sets.We also show this estimate in Fig. 3(c).It is well described by a single power law correction for L ≥ 16, and we use it for the most precise extrapolation of the transition point in the thermodynamic limit; Q c = 0.2758 (5).
With the estimated transition point at hand, the nature of the AFM-PVBS transition can be further studied.We find that both order parameters approach finite values with increasing system size at the estimated Q c value, as shown in Fig. 4.Here we also show another quantity characterizing the AFM state; the spin stiffness ρ S .In SSE calculations it is obtained using the simple formula [105] where n ± x denotes the number of S ± x,y S ∓ x+1,y operators in the SSE operator string.The stiffness also clearly extrapolates to a finite value.
Finite values of both the AFM and PVBS order parameters at a single point implies phase coexistence and a first-order transition.However, even though the extrapolated order parameters are substantial, the Binder cumulants of these quantities [Fig.3(b)] did not show the negative peaks that are typically present at first-order transitions.A negative peak in a cumulant originates from the Order parameters characterizing the AFM and PVBS states graphed vs the inverse system size at the estimated transition point, Qc = 0.27581.The curves show polynomial fits based on which infinite-size extrapolated values are obtained (the yellow squares shown at 1/L = 0).When not visible, error bars are smaller than the graph symbols.
double-peaked distribution of the squared order parameter; one peak at 0 reflecting the disordered phase and a second peak at a non-zero value reflecting the ordered phase.In the conventional classical case [109,110], the volume increase in free-energy barriers between the two degenerate phases implies that the negative peak grows linearly with the system volume, and in the quantum case a corresponding behavior is expected on account of the increasing tunneling barrier with the system size (and this has been observed at the phase transition in a staggered J-Q 3 model [53]).The tunneling barrier is absent if the order parameters form an enlarged spherical symmetry at the transition point, so that moving from one phase to the other corresponds to rotating the order parameter without energy cost.Such an unexpected mechanism at play at a first-order quantum phase transition was recently proposed to explain results for the transition between the AFM state and a two-fold degenerate PVBS in the CBJQ model [65], and subsequently other potential cases were also identified [84,85].

C. Emergent symmetry
Anticipating an emergent symmetry also in the present case, we proceed to examine the symmetry properties of the order parameters at the AFM-VBS transition.Fig. 5 shows the probability distribution of the PVBS order parameter P (Π 0 , Π 1 ) near the transition point and deeper inside the PVBS phase.Since we here focus on qualitative aspects of the angular dependence of the distributions, the (linear) color scales in the histograms are unimportant and not shown for simplicity.The fact that the histogram has four peaks on the horizontal and vertical axes once the system is sufficiently far inside the VBS phase, seen clearly in Fig. 5(d), implies that this phase is indeed a PVBS and not a CVBS.Furthermore, an emergent U(1) symmetry in the histogram at the critical point is clearly visible in Fig. 5(b).In Fig. 5(c) the distribution takes an approximate ring shape, with maximum weight away from the center, indicating a large magnitude of the order parameter.The distribution is still nearly U(1) symmetric.
Figure 6 shows the probability distribution of the dimer order parameter P (D x , D y ).In the VBS phase the four peaks of D are located at the diagonal angles, also indicating PVBS order (which implies coexisting D x and D y dimer order).The near-U(1) symmetry of the distribution P (D x , D y ) is similar to that seen in P (Π 0 , Π 1 ).
Emergent U(1) symmetry is a characteristic feature of the DQC phenomenon [49,88], and there is a divergent (on approach to the critical point) length-scale associated with the break-down of the symmetry inside the VBS phase.This length scale is reflected in an increasing sharpness of the four peaks of the order-parameter distribution when the system size is increased inside the VBS phase, which can be used to extract the exponent governing the relevant length scale [48] (in a way which has recently been further refined in the context of classical clock models [95]).However, as we showed in Fig. 4, in the present case the transition is clearly first-order, and the associated L dependence of the histogram features inside the PVBS phase should then be dictated by an exponent analogous to the dimensionality of a classical system [65].We have not studied the L dependence inside the PVBS systematically in the present case, but qualitatively the four peaks become sharper as L increases Distribution of the dimer order parameter P (Dx, Dy), with the components defined according to Eq. ( 5), accumulated in the same simulations as the plaquette order distributions in Figs.5(c,d); (c') and (d') are for Q = 0.28 and Q = 0.3, respectively.
(approaching δ-functions as L → ∞).We will next show that the emergent U(1) symmetry at the transition point actually reflects an even higher SO(5) symmetry, thus establishing the AFM-PVBS transition in the J-Q 6 model as another example of a symmetry-enhanced first-order transition.The spherical symmetry allows the AFM and VBS order parameters to continuously rotate into each other at fixed energy, which explains the absence of negative peaks in the Binder cumulants.
In the case of the CBJQ model, the O(3) symmetric AFM order and scalar (two-fold degenerate) PVBS order combine into an O(4) vector or SO(4) pseudo-vector.The methods used could not distinguish between the two possible cases [65] as physical reflections are not accessible in the simulations (while the meandering of the order parameter vector on the sphere is easily accessible).In the J-Q 6 model, we similarly expect O(5) or SO(5) symmetry of the components (Π 0 , Π 1 , m x , m y , m z ).This kind of emergent SO(5) symmetry was previously proposed within the the DQC scenario [70,82].
To demonstrate the unexpected enlarged symmetry at the first-order transition, we show the joint probability distribution of m z and Π 0 at and near the AFM-PVBS transition point in Fig. 7. Inside the AFM phase, Fig. 7(a), long-range order with O(3) symmetry implies a line segment when the five-dimensional distribution is projected down to the plane (m z , Π 0 ).Since we are near the transition, where the finite-size fluctuations of the amplitude of the order parameter are significant, the line is broadened.In the VBS phase, Figs.7(c,d), the fluctuation of Π 0 among two non-zero values (one positive and one negative) results in two blobs on the vertical axis.The central probability maximum visible in Fig. 7(d) corresponds to Π 0 being small when |Π 1 | is large, and the significant probability density remaining between the maximas show that tunneling between the four PVBS patterns is still rather prominent at this Q value for this system size.At the transition point, Fig. 7(b), we observe an U(1) symmetric histogram between m z and Π 0 after simple rescaling.Along with the exact SO(3) symmetry of the AFM order parameter and the numerically observed O(2) symmetry of the PVBS order parameter, this implies O(5) or SO(5) symmetry of the combined order parameter.
Note that, in constructing Fig. 7, we had to account for the fact that the definitions of the two different order parameters are associated with essentially arbitrary factors.The spherical symmetry at the transition point is apparent only if one of the order parameters is suitably rescaled by the ratio of the standard deviations of the two order parameters.Away from the transition point, no such rescaling results in a rotationally symmetric histogram.The features away from the transition point are most clearly visible if the scale factor is held fixed at its value calculated at the transition point, which is what we have done in Fig. 7.
We have carried out numerous tests of the angular uniformity of the histograms near the transition point, including those discussed in Ref. [65], for different system sizes.Here we present a different method, analyzing conditional probabilities to test the expected form of the radial distribution and to quantify some aspects of the angular distribution.We use the accumulated data points (m z , Π 0 ) for which m z = 0 to collect the conditional probability P (Π 0 |m z = 0), and similarly accumulate P (m z |Π 0 = 0) and P ( Π 2 0 + m 2 z |m z = Π 0 ).Here we point out that the computed equal-time values of m z and Π 0 in a given SSE sampled configuration are integerbased, and the above conditional probabilities are unambiguously defined.Related to this issue, note that the m z = Π 0 case does not correspond exactly to a diagonal cut in histograms such as Fig. 7(b), because of the different factors needed to scale from an elliptical to a circular-symmetric histogram.The required scale factor is not too far from unity, and m z = Π 0 without rescal-ing (which we consider for this analysis) corresponds to a radial line cut at an angle of about ≈ 70 • in Fig. 7(b).
The first-order transition with higher symmetry corresponds to the order parameter vector living on the surface on an O(5) sphere, in contrast to the case of a spherical symmetry at a critical point, where the radial distribution reflects critical fluctuations and are centrally peaked.The finite radius of the sphere reflects the nonzero magnitude of the coexisting order parameters in the thermodynamic limit.The radius of the sphere observed in simulations (defined using the sum of squared order parameters, including the scale factor discussed above) will not be constant because of finite-size effects, but will exhibit relative fluctuations that vanish when L → ∞.
Without fluctuations, the conditional probabilities defined above would be semi-circles if there is SO(5) symmetry.To take into account the radial fluctuations, we average the semi-circular distributions over Gaussian distributions of their radia.In practice, this is a oneparameter fit, because we can fix the resulting standard deviation of the distribution to that of the SSE computed histogram.We carry out this matching procedure only for one out of the three computed distributions, and rescale to the other two by a factor fixed by the standard deviations of those distributions.As shown in Fig. 8, the distributions obtained in this way for our largest system size, L = 48, agrees essentially perfectly with all the SSE computed conditional probabilities.As a contrast, we also show a Gaussian distribution, which has significantly fatter tails than the broadened semi-circle.
The conclusion drawn based on the results presented in this section is that the AFM-PVBS transition in the J-Q 6 model is first-order and associated with emergent O(5) or SO(5) symmetry, thus apparently similar to the first-order transition with emergent O(4) or SO(4) symmetry observed in the CBJQ model [65].Numerically, if no deviations from the symmetric distribution (or if the measures of symmetry continue to improve with increasing system size, as is expected for an emergent symmetry [65]) it is not possible to prove positively that the symmetry survives asymptotically as the system size increases, only that any perturbation that may eventually break the symmetry is weak and does not affect the system up to the largest sizes studied so far.

V. THE ALTERNATING VBS PHASE
When Q is further increased beyond the AFM-PVBS transition point, there is a second phase transition at Q Ac ≈ 0.934 into the eight-fold degenerate AVBS phase illustrated in Fig. 1(c).This phase transition can be described as a continuous freezing of the resonating plaquette singlets into either horizontal or vertical valence bond pairs in a checker-board pattern.The formation of the checker-board pattern involves the breaking of a Z 2 symmetry.In the AVBS phase, the plaquette order parameter Π does not change appreciably from its value in Fitting functions in black dotted lines are computed under the assumption of a uniform distribution over a 5-dimensional spherical surface with Gaussian fluctuation of its radius.By fixing the second moment of the resulting function to match the SSE data, we conducted a one-parameter fitting with the relative width of the Gaussian to match P ( Π 2 0 + m 2 z |mz = Π0).The resulting distribution was rescaled by the standard deviations of the other two histograms P (Π0|mz = 0) and P (mz|Π0 = 0), without any other adjustments.A simple Gaussian distribution with the same second moment as P (mz|Π0 = 0) is drawn as the light blue dashed curve.The inset shows the same distributions on a log scale.
the PVBS phase prior to the transition, and we are unable to detect this Z 2 symmetry breaking in Π (though potentially there could be a higher-order singularity).
To analyze the PVBS-AFBS transition we define the following four-component "alternation order parameter" where the notation (x, y) ≡ (a, b) stands for coordinates for which x ≡ a mod 2 and y ≡ b mod 2 for one out of the four possible combinations of the subscripts a ∈ {0, 1} and b ∈ {0, 1}.These combinations correspond to the four different plaquette patterns in the PVBS phase.
The "plaquette orientation" operator Ψ x,y in Eq. ( 10) is defined as and, thus, it can detect the dominant direction of the valence-bond pairs in the plaquette with lower-left corner at (x, y).Note that, while the order parameter A can detect the Z 2 breaking associated with the AVBS phase (where one out of the four components A a,b takes a nonzero value when the symmetry is broken), it remains zero through the AFM-PVBS phase transition.
Squared AVBS order parameter on lattices with open boundaries close to the PVBS-AVBS transition.A0,0 is the component in Eq. ( 10) that is favored by the boundaries.The values have been rescaled to test the expected critical scaling form A 2 0,0 = L −2β/ν f (δL 1/ν ), where δ = Q/QAc − 1 and β and ν are 3D Ising exponents (β ≈ 0.607, ν = 0.630 [113]).The scaling function f is exhibited by the data collapse for sufficiently large L. The inset shows the corresponding Binder cumulant U (1) 2 )/2, from which the critical point QAc ≈ 0.934 is obtained as the asymptotic location of the crossing points of data for different system sizes.
Since the system is already in a PVBS phase before the second phase transition takes place, in a large system, where the Z 4 symmetry in practice is broken in QMC simulations (i.e., the tunneling time between the different pattern is much longer than feasible simulation times [111]), it is possible to only consider the single component of A corresponding to the specific plaquette pattern realized.We will later, in Sec.VI, discuss the order parameter and symmetry breaking in more detail, but for now focus on the simplest way to detect and characterize the second transition.
The best way we found to study the PVBS-AVBS transition is to use open boundary conditions on L×L lattices with even length L, in which case the PVBS plaquette pattern is non-degenerate, i.e., the boundaries themselves induce a specific plaquette pattern that is stable in the limit L → ∞ if Q is above the critical value Q Ac for the AFM-PVBS transition (while for Q < Q Ac the pattern only remains at the boundaries and decays to zero in the bulk with increasing L).This setup allows us to monitor only the corresponding single component of the four-component operator A a,b in Eq. (10).The choice of boundary condition does not change the universality class or any other essential physics.We can also observe the transition in systems with periodic boundary conditions, though with considerably larger fluctuations that make it harder to study the critical behavior.We only discuss open-boundary SSE results in this section.
Since this phase transition involves an additional Z 2 symmetry breaking in the (2 + 1)-dimensional quantum FIG. 10.An order graph depicting the relationships between the eight degenerate AVBS states, represented as the yellow circles.The thick black bonds connect states that can most easily fluctuate between each other, while the thinner grey bonds correspond to the second-strongest fluctuations, with all being equal due to symmetries.States not connected by any bond have the smallest probability of fluctuating into each other.The PVBS phase forms when the tunneling barriers between states connected by the thinner bonds vanish and the system locks into one of the four groups of two strongly-connected states.Subsequently AVBS order forms when also the thick bonds vanish and the system locks into one of the two states in the previously chosen PVBS pair.The vectors displayed close to the circles show an example of a faithful Euclidean embedding of the graph in six dimension (the minimal dimensionality of the combined order parameter).
system, we expect the criticality to be that of the 3D Ising universality class.Fig. 9 shows the squared order parameter rescaled with the expected powers of the system size, with the estimated critical point Q Ac = 0.934 obtained using the Binder crossing method (shown in the inset of Fig. 9).Except for the smallest system size, the data points collapse well onto a single curve (the scaling function), thus supporting the expected Ising nature of the PVBS-AVBS transition.Thus, we have demonstrated that the PVBS state, which is a well established ground state of a quantum magnet even though it has previously not been found in sign-free models, can be unstable to a continuous freezing of the resonating plaquettes into static dimers with a checker-board pattern on the already formed plaquette pattern.The D 4 symmetry of the square lattice is then fully broken in two stages, first with Z 4 broken in the PVBS phase and then another Z 2 breaking in the AVBS phase.
A state similar to our AVBS state was previously discussed in the context of a quantum dimer model with multiple potential and kinetic terms in the Hamiltonian [112].Perhaps because the dimers lack the intrinsic singlet nature and are restricted to connecting nearestneighbor sites (while there are no analogous constraints in a spin model, where the dimers of a VBS are emergent objects), this dimer model only exhibited first-order transitions of the AVBS-like state, and there was no discussion of multi-stage symmetry breaking of the kind found here.We next discuss how to formally describe the twostage breaking of the lattice D 4 symmetry.

VI. SYMMETRY BREAKING AND UNIFIED VBS ORDER PARAMETER
We here construct a unified framework for describing the two stages of discrete symmetry breaking and, within that scheme, interpret the PVBS phase as an intermedi-ate phase.The first step is to construct an adequate order parameter for detecting both of the phase transitions in some way, unlike the previously discussed order parameters D, Π, and A, which capture only one of the transitions.In order to do that in a systematic and compact way, we introduce the concept of an order graph, which we will outline in the following and explain in more detail in Appendix B.

A. The order graph
The eight-fold degenerate ordered states in the AVBS phase have relative fluctuations depicted in Fig. 10, which is what we will call the order graph.In the thermodynamic limit, such fluctuations do not occur inside the AVBS phase, but they are manifested on any finite system and become more prominent as a quantum phase transition is approached.The thick black bonds represent fluctuations among states which are separated by the smallest tunneling barriers, and two states which have no bonds between them have the smallest probability of fluctuations among each other.
The states connected by thin bonds all have equivalent relations although they may not have the exact same symmetry transformations connecting them.For example, starting from the state marked by A in Fig. 10, transforming to state B or B corresponds to lattice translations in the +x and −x direction, respectively.These transformations are different, but since the parity symmetry between ±x directions is not broken, these pairs of states should have equal tunneling barriers separating them.The fluctuations between primed and unprimed states with the same letter are the easiest, because they do not involve changing the plaquettes on which the bond pairs form, only the orientation of the bond pairs within the plaquettes.The tunneling barriers between these state pairs vanish at the AVBS-PVBS transition.
In the PVBS phase the accessible subspace of the Hilbert space corresponding to one out of the four possible plaquette patterns is much larger than just the combined space of the two states represented by a pair of strongly connected yellow circles in Fig. 10, as the bond pairs on different plaquettes fluctuate essentially independently of each other (with some correlations that vanish with increasing distance between plaquettes) and connecting the two states by tunneling involves the creation of domain walls.The physical interpretation of the thick bond thus changes from representing the largest tunneling probability between the eight states in the AVBS phase on a finite lattice (which vanishes when the system size is taken to infinity) to an enlarged "basin" in the Hilbert space in the PVBS phase, roughly comprising all the states on what was previously the tunneling path.In the AFM phase the fraction of the Hilbert space involved in the unique ground state is further enlarged to encompass equivalently all the points in the graph-the AFM phase of course has completely different symmetry structure and fluctuations that are not described by the order graph for the PVBS and AVBS states.Some remnants of the graph structure may still be manifested in the short-distance fluctuations and plaquette and dimer correlation functions, at least close to the AFM-PVBS transition.
Any symmetry transformation that preserves the Hamiltonian will correspond to an automorphism of the order graph.Therefore, the order parameter is required to have the same transformation properties as that automorphism group.One way to construct such an order parameter is to embed the order graph into an Euclidean space in a faithful way, i.e., where all pair of vertices have equal distance if and only if they have the same type of bond in the order graph.This is possible in six dimension (or more) in our case, which is exemplified in Fig. 10 with possible coordinates shown as vectors.From the embedding, a useful six-dimensional order parameter can be constructed naturally, which we will explain in detail in Appendix B. The idea is to construct an order parameter that has eight peaks at those points embedded in six-dimensional space.The projection of the six-dimensional order parameter into lower dimensions corresponds to the 2D PVBS order parameter Π (and/or alternatively D) and the four-dimensional AVBS order parameter A.

B. Classification of broken symmetries
The order graph also provides us with a convenient way to correctly classify what symmetry is broken at a given transition.As we discussed in Sec.IV and further in Appendix A, the CVBS and PVBS phases both can be regarded as realizing a Z 4 symmetry breaking of the square lattice D 4 symmetry [88].The usual prescription for such symmetry breakings is to consider the symmetry group G the system originally transforms under to-gether with the remaining symmetry group H that it still transforms with after a symmetry-breaking transition has taken place.Then the quotient G/H corresponds to relevant symmetry group of the order parameter.This can be translated into the automorphism group of the order graph G , and the relevant subgroup H of that graph, which is the automorphism of the order graph when we distinguish some vertices.The distinction of vertices into two kinds corresponds to the state(s) the system gets stuck into and the remaining ones.This point of view is equivalent to classifying symmetry transformations into equivalent classes under their operations on the order parameters, as we discuss further in Appendix B.
The automorphism group of the order graph in Fig. 10 is D 4 Z 4  2 , corresponding to the rotations and reflections of the over-all graph (D 4 ) and 4 independent swaps (e.g.A ↔ A ) possible at each corner (Z 4  2 ).The automorphism group becomes Z 5  2 when two pairs of vertices (connected with thick black bonds) in the order graph are chosen in the PVBS phase.More precisely, the individual swaps in the corners remain, and D 4 breaks down to Z 2 because we are only left with one reflection that leaves the chosen pair invariant.Therefore, the broken symmetry is expressed by the coset D 4 /Z 2 .This is exactly the same as the broken symmetry of the four-state clock model, which is frequently, somewhat inaccurately (in a way explained in Appendix A), referred to as Z 4 symmetry breaking [36,88].When the system is in the AVBS phase, it corresponds to a single vertex being selected in the order graph, and now the remaining symmetry is Z 4  2 .Therefore, the symmetry breaking at the PVBS-AVBS transition is Z 5 2 /Z 4 2 = Z 2 (as also naively expected).In the case of a direct transition into the AVBS phase, the order graph symmetry reduces from D 4 Z 4 2 to Z 2 × Z 3  2 , and the broken symmetry becomes D 4 .This could also be intuitively understood in the following way.In the PVBS phase, all rotations around a lattice site does not conserve the macroscopic state anymore, but there is always a remaining spacial reflection with an axis going through a lattice site diagonally (whether ±45 • depends on which of the four PVBS states).This is the remaining Z 2 symmetry of the original square lattice point group symmetry D 4 , but also breaks when the system becomes AVBS, thus the entire From the above symmetry considerations, it is clear that, as long as the individual AFM-PVBS and PVBS-AVBS phase transitions are considered, a 2D clock order parameter and a scalar order parameter, respectively, are sufficient.Note that in the latter case, we assume that the system is already completely locked into one of the four degenerate PVBS patterns, as is the case in an infinite or, in practice, in a very large system in QMC simulations [111], or when using open boundary conditions as we did in Sec.V. Otherwise, a four-dimensional order parameter, e.g., as defined in Eq. (10), would be necessary for studying a transition between the PVBS and AVBS phases.To study a direct transition from a disordered state into AFVBS, a six-dimensional order parameter is preferable; we will consider that in Sec.VII C in the context of a finite-temperature paramagnetic-AVBS transition.

VII. FINITE-TEMPERATURE TRANSITION AND TRANSITION PATHS
Since the PVBS and AVBS phases break discrete symmetries they extend to finite temperature T > 0. We here discuss the nature of the finite-temperature transition and also summarize all the possible ways in which the symmetries can be broken simultaneously or in multiple stages leading to the AVBS phase, not just in the case of the J-Q 6 model but generically based on the order graph.
A. Transitions at T > 0 Figure 11 shows SSE results for the PVBS order parameter and the corresponding Binder cumulant versus Q at several different temperatures T .Here we can see clearer signs of first-order transitions as T is increased.Not only does the order parameter grow more acutely for higher T , but the Binder cumulant develops an increasingly prominent negative dip in the neighborhood of the transition point at higher T -such a dip is a sign of conventional phase coexistence, i.e., with no emergent higher spherical symmetry.
In this case, we indeed observe histograms of the PVBS order parameter with five peaks (see Fig. 12); four of them corresponding to the ordered PVBS patterns and the one at the origin corresponding to the absence of PVBS order in the paramagnetic state.We can con- clude that no emergent symmetry of the PVBS order parameter is present at T > 0, and the transition is firstorder with a conventional coexistence state with freeenergy barriers.This is in contrast to the CVBS case, where the T > 0 transition stays continuous, with exponents varying according to one of the critical branches of the Ashkin-Teller model, becoming increasingly similar to a Berizinsky-Kosterliz-Thouless transition as T c → 0 when the quantum-critical point is approached (e.g., the correlation-length exponent ν grows as T c → 0 + and is different from the DQC exponent exactly at T = 0) [114].
We provide the full T -Q phase diagram of the model in Fig. 13, which is based on several transition points (2) (21) FIG. 14.All 16 possible order graphs that satisfy the symmetry of the AVBS phase.The numbers 1, 2, and 3 correspond, respectively, to the strongest (black), second strongest (grey), and absent fluctuation paths in Fig. 10.Those fluctuations that survive in the thermodynamic limit are coded according to their relative strengths by the order of the numbers within (), with strongest to weakest fluctuations from left to right.In the corresponding order graphs, the fluctuation paths are drawn in black (strongest), gray (second-strongest), and dotted gray (weakest).In order to avoid cluttering, only a few of the weakest bonds are shown in the cases where they are present (in the six outermost order graphs).Symmetry breaking corresponds to the order graph becoming disconnected because of vanishing fluctuations, i.e., the graph breaks into equivalent clusters of connected yellow vertices (states), and the system is trapped in one of the equivalent clusters.These are shown in blue (one out of two clusters), red (one out of four clusters), or green (one out of eight vertices).The thick bonds connect order graphs that only differ by a single swap of adjacent fluctuation strength or by the disappearance of the weakest fluctuation.
on the paramagnetic-PVBS phase boundary, extracted using fixed-T scans such as those shown in Fig. 11.To determine points on the PVBS-AVBS boundary we used open-lattice calculations, as described in Sec.V.An important aspect of the phase diagram is that the AVBS phase is completely inside the PVBS phase.A direct paramagnetic-AVBS transition takes place upon lowering the temperature for Q very close to 1.
The fact that the T > 0 transitions become more strongly first-order as the transition temperatures increase may seem counterintuitive at first sight, but such behavior has been observed and analyzed in the context of fluctuation-induced first-order transitions of vestigial phases [25,86,115].We next provide another perspective on this phenomenon, utilizing the order graph introduced in the previous section.

B. Multi-stage symmetry breaking
We now discuss how transition paths between different states can be described systematically, which will provide generic insights into multi-stage discrete symmetry breaking also beyond the specific J-Q 6 model and D 4 symmetry considered here.
Figure 14 shows all possible order graphs which respect the final symmetry of the AVBS state.This figure shows connections (thick lines) between order graphs if and only if continuously varying the fluctuation strengths in the system can cause the order graph to change from one to another.The fluctuations can change by varying parameters in the Hamiltonian or by changing the temperature.We require these changes to not affect the symmetries, i.e., the three sets of fluctuation bonds introduced in Fig. 10  This kind of map, although conceptually simple, can restrict the topology of the phase diagram.The two assumptions (a) and (b) may appear to imply that we are considering only continuous transitions, where the evolution of the non-vanishing fluctuation strengths upon varying a parameter is smooth.Clearly, there can be first-order transitions where this assumption is not valid, but if a transition is weakly first-order, such that the correlation length of an order parameter is rather large even before the transition, and the transition itself is related to these fluctuations, we also expect the transition to respect the connectivity of Fig. 14.A transition that does not respect this connectivity, jumping directly between order graphs without direct connection, should be expected to be strongly first-order.
The sequence of symmetry breakings involved in the J-Q 6 model when it transitions first into the PVBS phase (either from the AFM phase at T = 0 or from the paramagnet at T > 0) and then into the AVBS corresponds to starting in Fig. 14 from the box marked (12) to the one marked (1), and then finally to (∅) at the center.If we add some (unknown) term in the Hamiltonian that in the notation of Fig. 10 enhances the fluctuations A ↔ B, A ↔ B , etc., the order graph should eventually change from (12) in Fig. 14 to (21).When the order graph is altered in this way, all eight states are still connected even when the second strongest fluctuation bonds have died out in the case (2), and, therefore, a phase transition will take place directly into the AVBS if the fluctuations represented by the strongest bonds are brought to zero and the graph (∅) is reached.This transition breaks the D 4 symmetry in a single step.
It is generally considered that breaking a symmetry which has a reducible representation, such as D 4 = Z 4 Z 2 , should only happen by way of a strongly firstorder transition [116].As we have seen, our arguments based on the natural connectivity of the order graphs suggest that the direct D 4 -breaking transition can also be continuous in principle and not just by fine-tuning.However, it is physically not a very likely scenario, because of the requirement that all the fluctuation paths in the graph (2) in Fig. 14 must be of equal strength, with the other ones vanishing.If we imagine embedding the order graph in a higher dimensional space such that the eight points are equidistant, a minimum of six dimensions is required (as discussed further in VII C).In naturally occurring systems we would not expect this A t p J l P q E y J 0 Q = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " 6 K + A m s 7 A a o E L a w 9 n y c g N x 6 B 4 k w A t p J l P q E y J 0 Q = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " 6 K + A m s 7 A a o E L a w 9 n y c g N x 6 B 4 k w A t p J l P q E y J 0 Q = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " 6 K + A m s 7 A a o E L a w 9 n y c g N x 6 B 4 k w

T < l a t e x i t s h a 1 _ b a s e 6 4 = " w i q d U 9 M Y H g J Y I w 2 X 5 P D d h o e l T 7 k = " > A A A C Z H i c h V H L S s N
c 7 + Z l 0 4 r m 6 Z e a 9 h i + 2 q s m / q e 7 q m e E x l 8 j u x B C X J j 5 m f Q A 5 A A k G k r d g N t r A L C x p q q E L A h M f Y g A K X W x k y C D Z z 2 2 g y 5 z D S / X 2 B I 0 R Y W + M s w R k K s x U e 9 3 l V D l i T 1 5 2 a r q / W + B S D u 8 P K G c z R I 9 1 S m x 7 o j p 7 p / d d a T b 9 G x 0 u D Z 7 W r F f Z O 9 H g y 9 / a v q s q z h 4 N P 1 Z + e P e x h x f e q s 3 f b Z z q 3 0 L r 6 + u F 5 O 7 e a n W v O 0 x W 9 s P 9 L a t E 9

e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " w i q d U 9 M Y H g J Y I w 2 X 5 P D d h o e l T 7 k = " > A A A C Z H i c h V H L S s N
c 7 + Z l 0 4 r m 6 Z e a 9 h i + 2 q s m / q e 7 q m e E x l 8 j u x B C X J j 5 m f Q A 5 A A k G k r d g N t r A L C x p q q E L A h M f Y g A K X W x k y C D Z z 2 2 g y 5 z D S / X 2 B I 0 R Y W + M s w R k K s x U e 9 3 l V D l i T 1 5 2 a r q / W + B S D u 8 P K G c z R I 9 1 S m x 7 o j p 7 p / d d a T b 9 G x 0 u D Z 7 W r F f Z O 9 H g y 9 / a v q s q z h 4 N P 1 Z + e P e x h x f e q s 3 f b Z z q 3 0 L r 6 + u F 5 O 7 e a n W v O 0 x W 9 s P 9 L a t E 9

e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " w i q d U 9 M Y H g J Y I w 2 X 5 P D d h o e l T 7 k = " > A A A C Z H i c h V H L S s N
c 7 + Z l 0 4 r m 6 Z e a 9 h i + 2 q s m / q e 7 q m e E x l 8 j u x B C X J j 5 m f Q A 5 A A k G k r d g N t r A L C x p q q E L A h M f Y g A K X W x k y C D Z z 2 2 g y 5 z D S / X 2 B I 0 R Y W + M s w R k K s x U e 9 3 l V D l i T 1 5 2 a r q / W + B S D u 8 P K G c z R I 9 1 S m x 7 o j p 7 p / d d a T b 9 G x 0 u D Z 7 W r F f Z O 9 H g y 9 / a v q s q z h 4 N P 1 Z + e P e x h x f e q s 3 f b Z z q 3 0 L r 6 + u F 5 O 7 e a n W v O 0 x W 9 s P 9 L a t E 9

e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " w i q d U 9 M Y H g J Y I w 2 X 5 P D d h o e l T 7 k = " > A A A C Z H i c h V H L S s N
c 7 + Z l 0 4 r m 6 Z e a 9 h i + 2 q s m / q e 7 q m e E x l 8 j u x B C X J j 5 m f Q A 5 A A k G k r d g N t r A L C x p q q E L A h M f Y g A K X W x k y C D Z z 2 2 g y 5 z D S / X 2 B I 0 R Y W + M s w R k K s x U e 9 3 l V D l i T 1 5 2 a r q / W + B S D u 8 P K G c z R I 9 1 S m x 7 o j p 7 p / d d a T b 9 G x 0 u D Z 7 W r F f Z O 9 H g y 9 / a v q s q z h 4 N P 1 Z + e P e x h x f e q s 3 f b Z z q 3 0 L r 6 + u F 5 O 7 e a n W v O 0 x W 9 s P 9 L a t E 9 A t p J l P q E y J 0 Q = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " 6 K + A m s 7 A a o E L a w 9 n y c g N x 6 B 4 k A t p J l P q E y J 0 Q = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " 6 K + A m s 7 A a o E L a w 9 n y c g N x 6 B 4 k A t p J l P q E y J 0 Q = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " 6 K + A m s 7 A a o E L a w 9 n y c g N x 6 B 4 k T < l a t e x i t s h a 1 _ b a s e 6 4 = " w i q d U 9 M Y H g J Y I w 2

X 5 P D d h o e l T 7 k = " >
c 7 + Z l 0 4 r m 6 Z e a 9 h i + 2 q s m / q e 7 q m e E x l 8 j u x B C X J j 5 m f Q x U e 9 3 l V D l i T 1 5 2 a r q / W + B S D u 8 P K G c z R I 9 1 S m x 7 o j p 7 p / d d a T b 9 G x 0 u D Z 7 W r F f Z O 9 H g y 9 / a v q s q z h 4 N P 1 Z + e P e x h x f e q s 3 f b Z z q 3 0 L r 6 + u F 5 O 7 e a n W v O 0 x W 9 s P 9 L a t E 9 c 7 + Z l 0 4 r m 6 Z e a 9 h i + 2 q s m / q e 7 q m e E x l 8 j u x B C X J j 5 m f Q x U e 9 3 l V D l i T 1 5 2 a r q / W + B S D u 8 P K G c z R I 9 1 S m x 7 o j p 7 p / d d a T b 9 G x 0 u D Z 7 W r F f Z O 9 H g y 9 / a v q s q z h 4 N P 1 Z + e P e x h x f e q s 3 f b Z z q 3 0 L r 6 + u F 5 O 7 e a n W v O 0 x W 9 s P 9 L a t E 9 c 7 + Z l 0 4 r m 6 Z e a 9 h i + 2 q s m / q e 7 q m e E x l 8 j u x B C X J j 5 m f Q x U e 9 3 l V D l i T 1 5 2 a r q / W + B S D u 8 P K G c z R I 9 1 S m x 7 o j p 7 p / d d a T b 9 G x 0 u D Z 7 W r F f Z O 9 H g y 9 / a v q s q z h 4 N P 1 Z + e P e x h x f e q s 3 f b Z z q 3 0 L r 6 + u F 5 O 7 e a n W v O 0 x W 9 s P 9 L a t E 9 c 7 + Z l 0 4 r m 6 Z e a 9 h i + 2 q s m / q e 7 q m e E x l 8 j u x B C X J j 5 m f Q x U e 9 3 l V D l i T 1 5 2 a r q / W + B S D u 8 P K G c z R I 9 1 S m x 7 o j p 7 p / d d a T b 9 G x 0 u D Z 7 W r F f Z O 9 H g y 9 / a v q s q z h 4 N P 1 Z + e P e x h x f e q s 3 f b Z z q 3 0 L r 6 + u F 5 O 7 e a n W v O 0 x W 9 s P 9 L a t E 9 15. Two phase diagrams containing AFM, PVBS, and AVBS phases, together with a possible Z2 symmetry broken phase.By analyzing the order graphs in Fig. 14, we can infer that case (a) is impossible without fine tuning or altering the symmetry, while (b) is generically possible.
fluctuation map to be easily realized, and a more likely scenario (which is realized in the J-Q 6 model), is a firstorder transition corresponding to a direct jump from the graph ( 12) to (2) in Fig. 14.
In contrast to the path in Fig. 14 where D 4 = Z 4 Z 2 is broken in a single step, the previously considered ( 12) → (1) → (∅) path (non-VBS to PVBS and then to AVBS) corresponds to a sequential breaking of the D 4 symmetry; first Z 4 and then Z 2 .At first sight, it may seem that the Z 4 symmetry breaking line and the Z 2 symmetry breaking line could cross each other, forming a purely Z 2 symmetry breaking phase as exemplified in Fig. 15(a), similar to a nematic phase.In fact, this is what happens in the vicinity of a tetracritical point [117] where four critical phase boundaries meet at one point.Also, the order graphs ( 13), (31), and (3) in Fig. 14 all correspond to a purely Z 2 broken phase.However, this does not imply that such phase diagrams with a tetracritical point are possible.To the contrary, the order graph analysis tells us that, in general, it is impossible to have such a tetracritical point when the AVBS phase is set as the final phase.This is because it is necessary that the order graph relations form a minimal loop that contains order graphs corresponding to all of the different phases surrounding a multicritical point, in this case paramagnetic, PVBS, AVBS and the Z 2 symmetry breaking phase.None of the minimal loops in Fig. 14, square or pentagon, has all of those four phases.Nevertheless, we have minimal loops that contain three phases, thus tricritical points such as those in Fig. 15(b It is of course possible in general to have a tetracritical point as in Fig. 15(a) with Z 2 and Z 4 symmetry breaking if we consider systems with less than the full D 4 symmetry considered so far.More specifically, if we allow order graphs such as those drawn in Fig. 16, a minimal pentagon loop contains four phases and we can have a phase diagram like Fig. 15(a).However, in order to have these order graph relations, we would need to specify eight of the type 2 bonds to be stronger than the rest, explicitly breaking the symmetry the system had when we were considering the AVBS phase and original D 4 symmetry.Note that all the eight states of the less symmetric 8-fold degenerate phase are still equivalent (i.e.vertex transi-)

FIG. 16.
Relationships between order graphs such as in Fig. 14, but for an eight-fold degenerate final phase arising in a model with lower symmetry than the J-Q6 model.Only those order graphs are shown that are needed to demonstrate the possibility of a tetracritical point such as in Fig. 15(a).
tive, as we will explain in Appendix B), but they need not correspond to the AVBS patterns that we have considered explicitly in the J-Q 6 model.
To demonstrate explicitly that we cannot achieve the phase diagram in Fig. 15(a) with the current symmetry of the AVBS phase, let us look back again on Fig. 10.Consider, e.g., the spatial reflection along the x axis that conserves the states of A, A , B, and B but inverts C and C (and also D with D ).This is an automorphism of the original order graph Fig. 10, but not for the order graph shown in Fig. 16 (12).Thus, while the order graph analysis rules out the possibility of such tetracritical points when considering the AVBS phase with current symmetry, in models with different symmetry, they are possible in general.A trivial example would be a ferromagnetic model with both Ising degrees of freedom {σ i } and 4state clock degrees of freedom {θ i }, with a Hamiltonian such as (12) In this case, the Z 2 and Z 4 symmetry breakings are decoupled and by tuning J Ising /J clock , we can trivially obtain a phase diagram with a tetracritical point.Note that now the type 2 fluctuations in Fig. 16 connect a vertex to only two neighbors, namely states with the same Ising order but with ±π/2 rotated clock order.The x axis reflection we considered above would correspond to some simultaneous transformation in the spins {σ i } and {θ i } for this model, which would not conserve the total energy.This shows that the naive effective Hamiltonian Eq. ( 12) with both Z 4 and Z 2 symmetry breakings does not explain the physics of the J-Q 6 model, and we would need an effective Hamiltonian that correctly captures the symmetry of the order graphs in Fig. 14 instead.
Even with strongly first-order transitions that do not have to obey the above rules for the adjacency of order graphs, we can still rule out generic phase diagrams such as Fig. 15(a), since it will require fine-tuning so that all four first-order transition lines come across at one single point.We can also predict that first-order transi-tions connecting order graphs that are farther away are usually stronger.This is because in general they would correspond to farther minima of the free energy.
Consider now again the J-Q 6 model.Ruling out a tetracritical point purely from the connectivity of order graphs is not the only aspect of the Q-T phase diagram that can be deduced with this tool before empirically running simulations.From Fig. 14, we see that if we want a purely Z 2 symmetry broken phase as in Fig. 15(b), then we would need to enhance the fluctuation of type 3, the weakest one in Fig. 10.However, even when we consider finite temperature, there is no reason to expect such fluctuations to become stronger.On the other hand, if we consider having large enough Q and gradually lowering the temperature, the fluctuations of type 1 should become relatively weaker, because they are quantum fluctuations of singlets induced by the J-term (remember that we set J + Q = 1 now).Therefore, if this relative weakening of the type 1 bonds continues enough so that the fluctuation strength ranking actually is reversed to (21), we would have a direct transition to the AVBS phase as explained earlier.As we argued before, the order graph (2) is unlikely to occur naturally in the J-Q 6 model, and it is more likely that there would be a direct transition from ( 21) to (∅) in such case, which should be strongly first-order.
Even if the effect of small J is not strong enough to cause such reversals in the fluctuation strength, it will definitely drag the symmetry breaking path to the left side of Fig. 14 rather than the right side where the purely Z 2 broken phase is.If there is a direct transition from (12) to (∅), this also must be strongly first-order since they are not connected.We will see in the following section that this is indeed the case, by observing the sixdimensional order parameter we construct from the order graph.
In either cases, it is expected that larger Q with thermal fluctuations will result in a direct paramagnetic-AVBS transition.Since the region where there is a direct paramagnetic-AVBS phase transition must be strongly first-order (because the transition is between nonconnected order graphs in Fig. 14), the nature of the paramagnet-PVBS phase transition should also become more and more strongly first-order as we approach closer to the paramagnetic-AVBS region, i.e., when we increase temperature.This would be the order graph perspective of understanding the so-called fluctuation induced firstorder transitions, where it has been argued previously that increasing temperature makes such transitions more strongly first-order [25,86,115].
We should note here that we have also assumed that intermediate phases always preserve symmetries that are remaining in the final phase i.e., no reentrant transitions are considered.We argue that this assumption, along with other requirements discussed above, are natural for analyzing the J-Q 6 model, and other models satisfying the above criteria may have an even richer set of phases and symmetry breaking paths according to Fig. 14.The order graph provides us with a compact way of enumerating all the possibilities, restricting the possible topology of the phase diagram with multicritical points.Furthermore, we can know that phase transitions connecting faraway order graphs will always be strongly first-order, and vice-versa: it provides good reasons to expect weak first-order transitions to respect the paths discussed here.
C. Unified order parameter for direct D4 symmetry breaking and other transitions Here, we further analyze the phase diagram of the J-Q 6 model by looking into the unified order parameter for PVBS and AVBS phases.As we mentioned in Sec.VI A and further explain in Appendix B, the order graph can be used to construct a unified order parameter that captures both PVBS and AVBS order.The unified order parameter, which we name V , is a six-dimensional vector and takes nonzero values in both the PVBS and AVBS phases, but in a different way.In the two phases, the histogram of the order parameter will have 4 and 8 peaks, according to the degeneracy of the phases.The unified order parameter is essentially the direct sum of the plaquette order parameter Π and the alternating order parameter A, defined in Eq. ( 3) and (10), respectively; Here the factors 8/7 and 2 are chosen to normalize all the components to become 1 in an ideal AVBS state that is perfectly dimerized.
For visualization, we project V into two dimensions in a way that clearly shows the eight-fold degeneracy.The following projection results in a histogram similar to the depiction of the order graph in Fig. 10: Fig. 17 shows the distribution P (X, Y ), where in panels (a) and (b) we can observe the coexistence of thermal paramagnetic and VBS states and distinguish the type of VBS the system transitions into from the number of peaks.The fact that Fig. 17(a) has five peaks, similar to Fig. 12(a), shows that the phase coexistence at Q = 0.9 is indeed between the paramagnetic phase (center peak) and the PVBS phase (four surrounding peaks).In contrast, in Fig. 17(b) there are nine peaks in total at Q = 0.95, implying a phase coexistence between the paramagnetic phase and the eight-fold degenerate AVBS phase.We can also see that the quantum fluctuation induced by decreasing Q is connecting the domains into four pairs in the order parameter space.Because of this, the four peaks of Fig. 17  Fig.6(c') but with the distribution projected in a different way.In Fig. 17(d) we can see that the strongest fluctuation path inside the AVBS phase are still within the four groups of two states.By observing histograms of the unified order parameter we can also confirm the corresponding order graphs at each point in the phase diagram.For example, the phase coexistence in Fig. 17(a) suggests a transition between the order graphs ( 12) and (1) of Fig. 14.As we argued in the previous section, in principle it is possible that the transition between these two order graphs is continuous, but it turns out that it is actually first-order in this case.Fig. 17(b) corresponds to a transition between ( 12) and (∅), and we have argued in the previous section that this must be a strong first-order transition since the two order graphs are not directly connected in Fig. 14.
Finally, in Fig. 18, we can see that the Binder cumulant of the unified order parameter, U V = 4 − 3 V 4 / V 2 2 , behaves far better than the cumulant of just the AVBS order parameter, U , when investigating the paramagnetic-AVBS transition with periodic boundary condition.These cumulants are again normalized so that they must necessarily converge to 1 as L → ∞ in the ordered phase (since both order parameters have a nonzero value with a fixed absolute value).However, even for L = 32 we have a completely negative U (4) A in the ordered phase, while we see a more rapid con- , 18. Temperature dependance of the Binder cumulants of the unified order parameter, UV , and the alternating order parameter, U A , at Q = 0.95 for different system sizes.
vergence of the final value of U V respect to the system size (though still not reaching very close to 1).We close this section by noting that the construction of the unified order parameter and the analysis of the symmetry breaking paths with order graphs should be in principle possible with standard Landau theory.However our approach with the order graph based on strengths of fluctuation paths (combined with symmetries) provides a more intuitive and compact way of understanding multistage discrete symmetry breaking.We expect it to be a useful approach in many other cases beyond the eightfold degenerate AVBS paths analyzed here.

VIII. DISCUSSION
In summary, the J-Q 6 quantum spin Hamiltonian introduced here exhibits several interesting physical phenomena.To begin with, it is, to our knowledge, the first sign-problem free 2D model hosting a four-fold degenerate PVBS ground state, thus enabling detailed QMC studies of the AFM-PVBS transition-a DQC candidate alongside the often studied AFM-CVBS transition.Instead of a continuous transition, we here found a clearly first-order transition, the coexistence state of which hosts an enlarged symmetry of the order parameter, with the combined three-component AFM order parameter and the two-component PVBS order parameter forming an SO(5) symmetric psudo-vector.We further discovered that the PVBS phase can be considered as a vestigial phase; an intermediate phase on the way to an eight-fold degenerate AVBS state which breaks a Z 2 symmetry in the already Z 4 -symmetry broken PVBS phase.The twostage breaking of two discrete symmetries stimulated us to introduce the order graph as an intuitive tool for analyzing such phase transitions.
We here further discuss some of the most intriguing aspects of our findings.In Sec.VIII A we discuss the nature of spinons in the PVBS state.At variance with the fracton scenario [24], we argue that the spinons are mobile and similar to those in a CVBS state.In Sec.VIII B we further discuss possible reasons for the first-order nature of the AFM-PVBS transition, using a different model with a first-order AFM-CVBS transition for comparison; the columnar J-Q 6 model.In Sec.VIII C we discuss the intriguing analogy between the symmetry-enhanced AFM-PVBS transition and the superconducting transition in the SO(5) theory of the high-T c cuprates.We discuss future prospects in Sec.VIII D.

A. Domain-wall structure and spinons
As we discuss in detail in Appendix A, in terms of the effective symmetry being broken the CVBS and PVBS states are dual and equivalent.There is still a possibility that the domain wall structure and types of spinons emerging as a defects in these two states may be different, and this could be the reason for the first-order transition found here.Recently, it was argued that the AFM-PVBS transition should in general be first-order and not a continuous DQC transition [24], as a consequence of the spinons in the PVBS phase being immobile fractons.At first sight, this conjecture appears to be con- firmed by our results.However, we do not believe that the spinons in the PVBS considered here are fractons, and we will argue that the fracton scenario in general is very unlikely in quantum magnets.
In Ref. [24] the singlet plaquettes of the PVBS state were treated as rigid objects, and it was pointed out that a spinon (a site not belonging to any plaquette) forming at the nexus of four domain walls in such a state cannot move because the constraints prohibit local fluctuations that gradually shift the spinon together with some of the plaquettes.Our primary objection to this scenario is that the singlet plaquettes in an actual PVBS state in a quantum spin system are not rigid objects, except perhaps in some extreme case that is not likely to be realized in practice with a naturally-arising spin Hamiltonian.The PVBS state hosts various quantum fluctuations, and when viewed in the valence-bond basis the singlet plaquettes are resonating pairs of horizontal and vertical valence bonds on 2 × 2 sites on the square lattice.There are also longer valence bonds, though they likely are less important for capturing the essential physics.In a minimal effective description, the PVBS state should be viewed as a mixture of rigid plaquettes and dimers, and there must be processes converting a plaquette into a pair of dimers and vice versa.There are then two possible types of domain walls, and, therefore, two kinds of spinons.We illustrate these in Fig. 19.As in the DQC scenario [36,88], in the same way as a domain wall between CVBS phases can be regarded as consisting of plaquette singlets in a PVBS pattern, a domain wall between PVBS states should be primarily made up out of dimer singlets in CVBS pattern.Whichever type of domain wall is realized in a given system should depend on the details of the Hamiltonian.
In either type of domain wall, in the presence of dimers the spinon can also move locally inside the domain wall through processes where a dimer and a spinon are shifted together as depicted in Fig. 20.In both cases, the first step only involves the spinon itself moving as a consequence of the action of a single J-term (singlet projector) P (i, j).Note that, with a bipartite Hamiltonian, the spinon is constrained to move only within a given sublattice.The second step exemplifies how the domain walls around the unpaired spin can adapt in order to form the configuration shown in Fig. 19 again with the new spinon position.All the changes in this second step are plaquette-dimer conversions, or a pair of dimers resonating to point in the other direction.The dimer resonance is easily induced by a single J-term, and plaquette-dimer conversions happen even more easily because of the large overlap between the two states.Considering these processes, the domain wall itself can fluctuate, and, thus, the collective object (the spinon) consisting of the meeting point of the domain walls and the unpaired spin can move through the lattice and should not be regarded as a fracton.

B. First-order transitions and emergent symmetry
The arguments in the preceding section make the fracton scenario unlikely as an explanation for the first-order AFM-PVBS transition, and we therefore discuss other possible reasons here.First, we note that the DQC scenario itself does not guarantee that the AFM-VBS transition is necessarily continuous (with the VBS being either a CVBS or a PVBS); the claim is that the transition is generic and not a fine-tuned multi-critical point, but first-order transitions for some Hamiltonians can never be excluded [36].Second, the presence of an emergent SO (5) symmetry observed here at the transition point may suggest that the system is close to a DQC point with such a higher symmetry [23], so that expected perturbations destroying that symmetry is weak and only observable on large lattices.Another possibility is that the DQC phenomenon is connected to first-order transitions with asymptotically exact emergent symmetries.Both of these scenarios were discussed in the context of first-order transitions with apparent O(4) symmetry in systems with AFM-PVBS transitions where the PVBS state is only two-fold degenerate [65,84].While an ex-act emergent symmetry at a first-order transition would be outside previous expectations, this scenario has some support also in a recent work claiming exact emergent supersymmetry at certain phase transitions with fermionic and bosonic degrees of freedom [85].Here we do not provide any definite conclusions on the root causes of these unusual first-order transitions, but we address a related issue using illuminating computational results.
The first issue concerns the possibility of first-order transitions into the columnar CVBS state.So far, firstorder transitions have been observed in spin-anisotropic J-Q models in which the AFM order parameter is O(2) symmetric, where the discontinuities weaken as the anisotropy is decreased and there may be a change to a continuous transition at some critical value of the anisotropy [62,64].In spin-isotropic systems such as the J-Q 2 and J-Q 3 models with a columnar arrangement of the singlet projectors, no clear-cut signs of firstorder transitions have been observed, though there are unusual scaling corrections [52,59] that have also been interpreted as a consequence of a weak first-order transition [23,50,57,73,74].By introducing other types of Q interactions in a J-Q 3 model with columnar Q 3 interactions, it was recently found that the DQC point can evolve to a clearly first-order transition [122].It is also interesting to study columnar J-Q n models with Q terms containing more than n = 3 singlet projectors.Here we consider the n = 6 case, i.e., the same number of projectors as in the model studied previously in this paper but with a different spatial arrangement of the projectors.
As shown in Fig. 21, the columnar J-Q 6 model exhibits CVBS order for large Q, but, unlike the columnar J-Q 2 and J-Q 3 models, the AFM-CVBS transition here is strongly first-order.Evidence for this type of transition is seen in the CVBS Binder cumulant, which in Fig. 21 develops increasingly negative peaks as the system size increases.We also observe (not shown here) coexisting finite AFM and VBS order parameters similar to Fig. 4. We believe that the first-order transition here is caused by the large number of coupled spins in the Q 6 terms, which may cause the system to nucleate PVBS order locally on large enough "droplets" to cause an instability of the AFM state.
Having identified a first-order AFM-CVBS transition, we can now address the issue of the generality of emergent symmetry at first-order transitions at AFM-VBS transitions.Since the Binder cumulants in Fig. 21 exhibit the tell-tale signs of conventional phase coexistence, we do not expect any emergent spherical symmetry.Indeed, there are no signs of emergent higher symmetries in order-parameter distributions.Fig. 22 shows the distribution of the CVBS dimer order parameter P (D x , D y ) at four Q values close to the transition point.Here we see how the distribution evolves from a central peak in the AFM phase (where the CVBS order parameter is peaked at 0) to a four-peak distribution reflecting the four-fold degeneracy in the CVBS phase.In a narrow window close to the transition we observe [Fig.22(b,c)] five peaks re- flecting coexistence of the AFM and CVBS phases.
On the one hand, these results for the columnar J-Q 6 model support a standard nucleation mechanism at play in this case, as no emergent symmetry is present.On the other hand, in the plaquette J-Q 6 model we have the same number of interacting spins but do observe emergent symmetry, which speaks against the standard nucleation scenario.The continuous nature of the PVBS-AVBS transition we observe also shows that nucleation does not necessarily follow from a large number of interacting spins.Nevertheless, in this case, a plaquette consisting of four spins becomes one degree of freedom in the effective Ising model, so the Q-term may be considered as an effectively three-body term.
In any case, the emergent symmetry is clearly not uni-versal at general first-order AFM-VBS transitions (see also Ref. [81]), and the possibility also remains that the emergent symmetry at the AFM-PVBS transition in the present case only holds up to some length scale larger than the systems studied here.If so, the transitions in the two different J-Q 6 models may be qualitatively the same but differ in the length scale at which conventional coexistence is apparent-in which case it still is remarkable and unexpected to have SO(5) symmetry up to such large length scales in the plaquette case.It is interesting to note that the two J-Q models where emergent symmetry of the coexistence state has been observed both have PVBS states; the two-fold degenerate one in Ref. [65] and the four-fold degenerate case studied in the present paper.This may be an indication of a symmetry-breaking perturbation that exists (or is strong) at first-order AFM-CVBS but vanishes (or is very weak) at AFM-PVBS transitions.

C. SO(5) theory of high-Tc superconductivity
The possible emergence of SO( 5) symmetry in condensed matter systems has received significant attention due to the fact that phases with O(3) AFM order often exist adjacent to superconducting phases, which break U(1) symmetry.One may then speculate that the two types of orders share a common origin in a unified degree of freedom that collectively can rotate between the two phases [22].The SO(5) scenario for high-T c superconductivity [21] postulates that doping away from the half-filled-band, where the cuprate materials are AFM insulators, eventually leads to a "flop" on an SO(5) sphere from the a direction spanned by the three AFM components into the plane spanned by the superconducting order parameter.This mechanism is very similar to what we have discussed here for the transition of the AFM into the PVBS state.
In the case of the cuprates, to study the SO(5) scenario with numerical simulations, the underlying Hubbard or t-J model first has to be projected down to an effective model (because the electronic models are too difficult to study on sufficiently large length scales), which is bosonic and can be simulated with QMC methods.Such studies were carried out with the SSE method in Ref. [118].Though a first-order transition was identified at a critical doping fraction, the coexistence state did not exhibit SO(5) symmetry, but instead conventional phase coexistence was found.It was argued that long-range Coulomb interactions might eventually act against phase separation and lead to a quantum-critical point.This scenario thus differs from the first-order coexistence state with SO(5) symmetry of the J-Q 6 model at the AFM-PVBS transition, where there is neither phase separation nor conventional criticality.
Experimentally, the existence of an excitation mode at 41 meV, detected in inelastic neutron scattering experiments, was taken as support of the SO(5) scenario [21,119], but arguments to the contrary have also been voiced [120].Since the J-Q 6 model has a different, more exotic kind of coexistence state than what was found in the projected SO(5) model, it would be interesting to study also the dynamical spectral functions of the J-Q 6 model (which can be done with SSE simulations in combination with numerical analytic continuation methods [76]) in and close to the coexistence state.It is tempting to speculate that the SO(5) predictions for the cuprates would come out differently with the exotic coexistence state, and further studies of the J-Q 6 model may serve as an analogy where reliable results can be obtained.

D. Future prospects
Our work presented here illustrates the power of the J-Q designer Hamiltonian approach in engineering sign-free Hamiltonians with exotic quantum states and quantum phase transitions.The results and remaining open issues prompt several possible follow-up studies, some of which we summarize here.
It would clearly be interesting to design a J-Q model with a continuous AFM-PVBS transition, which was our initial goal.This may not be so easy however, as we have so far not even succeeded in creating a PVBS state with less than six singlet projectors (with Q 4 terms similar to our Q 6 terms in Fig. 2 leading to CVBS states).If nucleation due to a large number of coupled spins is indeed the root cause of the first-order transition (which is not clear, as discussed in Sec.VIII B), an even larger number of singlet projectors will likely take us even further away from the DQC scenario.It would still be worth trying, e.g., Q 8 interactions defined on 4 × 4 lattice sites.We here also note that the PVBS state can be regarded as a "superposition of superposition states" (locally resonating valence bonds), and this specific aspect of the state seems to be what makes it more difficult to realize a PVBS than a CVBS-at least within the J-Q approach.The broader family of models by Kaul [121], expanding on the J-Q approach, should also be explored more extensively, though our initial attempts have not been successful in producing PVBS states with less than twelve interacting spins.
Another interesting prospect is to mix different kinds of Q terms and investigate the interplay between interactions favoring PVBS and CVBS states (in a way similar to what was recently done with competing interactions favoring CVBS and staggered VBS states [122]).It may be possible in this way to design a Hamiltonian which can be fine-tuned to have an eight-fold degenerate VBS state, with both CVBS and PVBS order.In histograms such as those in Fig. 5 this kind of state would give rise to eight equidistant peaks, and the phase would break Z 8 symmetry.A potentially continuous transition between the AFM state and such a VBS would be very interesting, considering that clock-like Z q perturbations should be increasingly irrelevant as the number of states q is increased [90,95].It may then be possible to observe DQC physics with less anomalies and scaling corrections than with the spin models studied so far [59].Recently, a fermionic model argued to have a DQC-type transition without any discrete perturbation was studied [123], but in this case the system sizes are very small because of the unfavorable scaling of the fermion determinant QMC algorithm.
The finding here of a first-order transition also in the columnar J-Q 6 model calls for more systematic studies of columnar J-Q n models.Previous studies for n = 2 [49][50][51][52][56][57][58][59][60] and n = 3 [48] point to continuous transitions, and it would be useful to systematically increase n and find the smallest n for which the transition is clearly firstorder.Such a study might be helpful in determining whether the transitions for small n are truly continuous or only very weakly first-order.
The order graph approach we have introduced here to analyze multi-stage discrete symmetry breaking should have wide applicability to both classical and quantum phase transitions.The order graph provides a natural set of possible intermediate phases, and specifies in what ways they can be adjacent in the phase diagram.It would be interesting to explore possible phase diagrams in other challenging models, e.g., the classical frustrated spin models discussed in Refs.[124,125], which exhibit (possibly multiple-stage) discrete symmetry breaking with large numbers of states.
Here, translation of (s, t) in the (x, y) direction is denoted as T (s,t) and we categorize these according to the parity of s and t (mod 2).M x is reflection with respect to the x axis, and R θ is θ-rotation around the site at the origin.Any symmetry transformation of the square lattice can be written in the form R θ M {0,1} x T (s,t) , and thus belongs to one and only one of the above equivalence classes.
The above eight equivalence classes form a group isomorphic to D 4 , which exactly corresponds to the automorphism group of the order graph for four CVBS states-the order graph is described in Sec.VI the main text and in further detail in Appendix B. This is also equivalent to saying that, within each equivalence class, the transformations have exactly the same effect on the order parameter.For example, all transformations in a correspond to a π/2 rotation in the order parameter space of D and Π, and b results in a reflection along a diagonal axis.When the system is in the CVBS phase, the remaining symmetry is {I, ab} or {I, a 3 b} depending on which of the four states the system is in.
An important point is that neither {I, ab} nor {I, a 3 b} is a normal subgroup of D 4 .Therefore, the symmetry being broken here is technically D 4 /Z 2 , which is just a coset and not a quotient group.In the most strict sense, the symmetry breaking is neither a Z 4 nor Z 2 × Z 2 symmetry breaking.Nevertheless, we can still express the coset D 4 /Z 2 in a way that looks like Z 4 : e.g. the left coset is D 4 / L {I, ab} = (A1) {I, ab}, {a, a 2 b}, {a 2 , a 3 b}, {a 3 , b} , and we can choose the representatives as {I, a, a 2 , a 3 }.This corresponds to the fact that we can express the four degenerate CVBS states as and R 3π/2 ∈ a 3 .However, this is just because we chose an expression that is based on the Z 4 subgroup of D 4 for the left coset D 4 / L Z 2 .We can also choose the representatives to be {I, a 2 b, a 2 , b} which now corresponds to expressing the four degenerate states as which has a Z 2 × Z 2 structure.This ambiguity directly originates from the fact that D 4 /Z 2 is a coset, and either representations are equally valid.Thus, the terms "Z 4 symmetry breaking" and "Z 2 × Z 2 symmetry breaking" are both valid in the sense that subgroups of the original symmetry group G of the lattice that are isomorphic to Z 4 or Z 2 ×Z 2 are completely broken in the ordered phase.
The situation is almost exactly the same for the PVBS phase, including the way the symmetry transformations are divided into equivalence classes.The only difference is that the remaining symmetry in the PVBS phase is either {I, b} or {I, a 2 b}, depending on the symmetry broken state.Similarly to the case of CVBS phase, neither of the two are normal subgroups of D 4 and allow multiple representations of D 4 / L Z 2 .In this case, the left coset is D 4 / L {I, b} = {I, b}, {a, ab}, {a 2 , a 2 b}, {a 3 , a 3 b} , (A2) and the Z 4 and Z 2 ×Z 2 representations are {I, a, a 2 , a 3 } and {I, ab, a 2 , a 3 b}, respectively.We can represent the four different PVBS states as either Z 4 -like using the representation {I, a, a 2 , a 3 }: or Z 2 × Z 2 -like using the representation {I, ab, a 2 , a 3 b}: analogously to the CVBS case.
We conclude that the symmetry breaking of the CVBS and PVBS phases have exactly the same structure: D 4 /Z 2 .We have shown that this becomes apparent when we group all the symmetry transformations into equivalence classes according to how they act on the order parameters.This approach is practically very simple if one starts from the automorphism group of the order graph.The difference between the two phases amounts to whether the remaining symmetry is {I, b}&{I, a 2 b} or {I, ab}&{I, a 3 b}.The "duality" between CVBS and PVBS is most clearly understood when we observe that the remaining symmetries of Z 2 for either of the phases are actually isomorphic via automorphisms, e.g., a = a, b = ab of the group D 4 itself.To be more precise, this renaming will not affect the relations of the elements in the group (e.g.b a = a 3 b just as ba = a 3 b), and simply corresponds to changing which axis to choose for the reflection represented by b.In this sense, we can say that the effective symmetry breaking in the CVBS and PVBS phases are exactly the same, which can be expressed with Z 4 or Z 2 × Z 2 .
Note that the ferromagnetic clock model with four states exactly follows the above classification of symmetries, and can actually be seen as both Z 4 or Z 2 × Z 2 symmetry breaking.This is consistent with the fact that when we consider a "hard" clock model, where the spins can only point in discrete direction (instead of using a cosine potential), the model reduces to two decoupled Ising models, which obviously should be able to be seen as a Z 2 × Z 2 symmetry breaking.This ambiguity essentially comes from the peculiarity of the group D 4 that does not appear in D q with q > 4. Thus, for clock models with the number of states q larger than 4, the broken symmetry can only be classified as Z q .
Appendix B: Rigorous construction of the order graph and its embedding in Euclidean space Here, we precisely define the construction of the order graph introduced in Sec.VI and also explain exactly what we mean by a faithful embedding of it in Euclidean space.The topic of embedding graphs into Euclidean spaces of broad interest in fields ranging from machine learning [126] to pure graph theory [127], but the embedding we consider here has a constraint which has not been considered in depth before.
Let us consider a discrete symmetry breaking with possibly multiple steps.By assuming that it is always a discrete symmetry that is being broken, we will have only finitely many ordered (or disordered) states at any point.We consider a situation where we lower the (classical or quantum) fluctuation, starting from a disordered phase with no spontaneous symmetry breaking.Let us assume that, after one or several symmetry breakings there are M degenerate ordered states in the final phase, and that we already know those M states.

Definition: Order Graph
An order graph G = (V, E) for a multiple discrete symmetry breaking consists of the vertex set V where each of the vertices corresponds to one of the M ordered states in the final phase, and an edge function E. The edge function E : V 2 → {1, 2, . . ., K} is a symmetric function that tells what type of bond a pair of vertices have between them.The edge function will be defined according to the strength of fluctuation between the two ordered states.
Whereas the usual definition of a graph involves the edge set E ⊂ V 2 , our definition is a slight extension of it since now E is a function that specifies what type of relation two ordered states have out of K possibilities.Usual graphs can be considered as the case where K = 2, corresponding to the two possible relations either having or not having an edge.The strength of fluctuation between two states v 1 , v 2 in V can be quantified in several different ways.One is to compute the domain wall (free) energy between the two states.Another way is to evaluate the transition probability from one state to another under some local dynamics, e.g.local Monte Carlo updates.In a quantum state, tunneling amplitudes can be calculated in principle and are reflected in the probability distribution between peaks in histograms such as those in Figs. 5 and 6.
All of these methods should classify all M (M − 1)/2 pairs of states into groups with the same "distances".Then E(v 1 , v 2 ) will be set to 1 if they have the strongest possible fluctuation, and 2 if it is the next strongest, and so on.Exactly how the strength of the fluctuations are quantified may in some cases be further refined, but the essential point here is that we can in principle construct a well-defined function E, which orders and classifies the pairs of states according to their physical fluctuation strengths.For the model we analyze in this paper, if we consider the two-stage discrete symmetry breaking in the sequence AFM-PVBS-AVBS (or paramagnetic-PVBS-AVBS), then M = 8 and K = 3.We only draw edges corresponding to 1 and 2 in Fig. 10 for visualization, and all pairs of vertices without a drawn edge corresponds to the third type.
An automorphism on the order graph G is a map from V to itself f : V → V which satisfies Since this graph is representing the relations between the ordered states, which should have exactly the same (free) energy, this graph must be vertex transitive, which means that all vertices are in a way equivalent, i.e., The graph also must be edge transitive as well, since all pairs with the same value of E should be equivalent in the same way, i.e., Any symmetric transformation to a state that preserves energy corresponds to an automorphism of the order graph.To demonstrate this, let us take the model we have studied in the main text as an example.The order graph is as shown in Fig. 10, and we have M = 8 states in the final AVBS phase labeled from A to D .If the system locks into the pattern labeled A, a shift in the +x direction in unit distance transforms the state into pattern B. The state will become B if the shift is in the −x direction instead.Since the system is not breaking any symmetry between the +x direction and −x direction, the physical fluctuation between (A ↔ B) and (A ↔ B ) should be equivalent.This explains E(A, B) = E(A, B ).
The shift +x also maps other states as well, resulting in the automorphism of (A, B, A , B )(C, D, C , D ).Note that this permutation satisfies Eq. (B1).As another example, a π/2 clock-wise spatial rotation in the depiction will correspond to the automorphism (A, A )(B, D, B , D ).An ideal order parameter will have the same set of symmetries as the Hamiltonian does, which means it should respect the symmetry of the automorphism group of the order graph, as we explain next.When we think of an order parameter which is an d dimensional vector, we can regard it as a point in the d dimensional Euclidean space E d .If there are M different ordered states, they should correspond to M different points x 1 , x 2 , . . ., x M in E d , and these points should have equal distance from the origin.We set ∀i, |x i | = 1 for simplicity.Furthermore, let us think of two ordered states corresponding to x i and x j , where they are "adjacent", meaning that the fluctuation between those two ordered states are the strongest.This means that the corresponding vertices in the order graph v i and v j have a connecting edge E(v i , v j ) = 1.Ideally, all states (x k for example) which have the same relation to x i as x j does, should have the same relation, which requires |x i − x j | = |x i − x k |.Now, if we think of a symmetric transformation of the Hamiltonian, it transforms one ordered state to another one in general.Note that a symmetric transformation which does not change any ordered state to another corresponds to a symmetry that does not become broken even after a (partial) symmetry breaking has taken place.As we have exemplified in the previous paragraph, we can express the transformation by a permutation σ among M elements.Such transformations translate to isometric transformations of the M points in E d , corresponding to some rotation and/or reflection F .This is because the relative relation between states should not change under these transformations, and thus |x i − x j | = |x σ F (i) − x σ F (j) | holds for any F and i, j.If we want to construct an order parameter that reflects the broken symmetry, all symmetries of the Hamiltonian (note that they are always expressible by automorphisms of the order graph) should have a corresponding isometric transformation in the order parameter space E d .
From the above argument, we arrive to a general way for constructing order parameters for discrete symmetry breaking, by a faithful embedding of the order graph defined as the following.

Definition: Faithful Embedding
A faithful embedding of a graph G = (V, E) to a Euclidean space E d is a mapping Γ : V → E d of vertices to points, such that for all automorphism f of G, there is an isometry F in E d that satisfies ∀v ∈ V, Γ(f (v)) = F (Γ(v)).
With this embedding, order parameters could be defined to be the sum of the vectors Γ(v i ) = x i corresponding to each microscopic degrees of freedom in a given configuration.This is possible because global symmetry breaking always have long range order, meaning that looking at local degrees of freedom is enough for determining which state the system is locally in.
Let's consider the AVBS phase as an example.In this phase, it is sufficient to look at four spins in a plaquette to determine which of the eight-fold degenerate states the system is in.More precisely, if we observe two adjacent and parallel dimer singlets, depending on which direction (vertical or horizontal) they are pointing and which of the four plaquette patterns they are in, one of the eight states is uniquely determined to be associated with them.This corresponds to the local degrees of freedom determining the state, as explained above.Once we obtain the ordered state i that corresponds to a particular local configuration, then the vector Γ(v i ) = x i is the order parameter value that corresponds.Thus, a sum over all the sites of such vectors will be the most natural order parameter.In practice, since we use the z basis representation in our SSE simulation, it is more practical to define the correspondence between a local z basis configuration and the eight AVBS states.This essentially makes the resulting six dimensional order parameter to be equivalent to the direct sum of the plaquette order parameter Π and the alternating order parameter A defined in Eq. ( 3) and (10), respectively.
It should be noted that we only deal with discrete symmetry breaking here.Continuous symmetries are naturally associated with a corresponding perturbation in the physical system, and has a more direct physical meaning.For example, two ordered states which are connected with an infinitesimally small continuous symmetric transformation (Lie group) also are related with a perturbation (Lie algebra).Two spatially separated regions with different ordered states can be smoothly connected with a continuously varying order parameter, which is exactly the manifestation of the continuous symmetry.This corresponds to some shear in the order parameter space, such as spin stiffness.Discrete symmetry breaking on the other hand, are less connected to specific physical properties, since two ordered states which are connected with a symmetry transformation do not have a smooth connection between them characterized by a parameterized transformation of the symmetry in general.Discrete symmetry groups can be regarded as a subgroup of a continuous symmetry group, and we can expect that they become a very good approximations to the continuous symmetry when the subgroup is "large enough".This is indeed the case of q-state clock models with q large enough, which the discrete Z q symmetry becomes a good approximation of an O(2) symmetry which corresponds to an XY model.This results in those clock models with large q to have emergent XY criticality.Our method provides a systematic and intuitive way to obtain the possible emergent symmetry by analyzing the relation of the ordered states with a graph.
The embedding method discussed here is very close to that in Ref. [128], which considered only the case of M = 2 and furthermore imposed an additional constraint that edges are not allowed to cross each other in the embedding.The faithful embedding defined for our purpose is also interesting in it's own right, revealing interesting properties of vertex transitive graphs.For example, the minimum embedding dimension of a Petersen graph in the way we defined here is 5, which is nontrivial and also coincides with the embedding dimension defined in a more algebraically way [129], suggesting possible connection.

FIG. 1 .
FIG. 1. Depiction of the different types of VBS states discussed in this work; (a) four-fold degenerate CVBS, (b) fourfold degenerate PVBS, and (c) eight-fold degenerate AVBS.Long ovals in (a) and (c) represent spontaneously "frozen" singlets and the ovals forming squares in (b) represent plaquette singlets, which can be expressed as resonating pairs of horizontal and vertical singlet pairs.

FIG. 7 .
FIG. 7. Joint probability distribution P (mz, Π0) of the AFM and VBS order parameters for system size L = 48.The Q values are the same as in Fig. 5; (a) slightly inside the AFM phase at Q = 0.273, (b) at the transition point Qc = 0.2758, (c) slightly inside the PVBS phase at Q = 0.28, and (d) well inside the PVBS phase at Q = 0.3.

FIG. 8 .
FIG.8.Three conditional probability distributions of the PVBS and AFM order parameters; P (Π0|mz = 0) (blue), P ( Π 2 0 + m 2 z |mz = Π0) (purple), and P (mz|Π0 = 0) (red).The latter two histograms are shifted to the left for clarity.Fitting functions in black dotted lines are computed under the assumption of a uniform distribution over a 5-dimensional spherical surface with Gaussian fluctuation of its radius.By fixing the second moment of the resulting function to match the SSE data, we conducted a one-parameter fitting with the relative width of the Gaussian to match P ( Π 2 0 + m 2 z |mz = Π0).The resulting distribution was rescaled by the standard deviations of the other two histograms P (Π0|mz = 0) and P (mz|Π0 = 0), without any other adjustments.A simple Gaussian distribution with the same second moment as P (mz|Π0 = 0) is drawn as the light blue dashed curve.The inset shows the same distributions on a log scale.

2 〉FIG. 11 .
FIG. 11.Examples of constant-T , Q-dependent PVBS order parameters (inset) and the corresponding Binder cumulants (main graph) for L = 32 systems.The negative peaks in the cumulants indicate first-order transitions.

FIG. 12 .
FIG. 12.Probability distributions of the plaquette order parameter P (Π0, Π1) in (a) and the dimer order parameter P (Dx, Dy) in (b) at the transition value of Q (defined here as the position of the negative peak in the Binder cumulant-see Fig.11).The system size is L = 32.
are maintained.The only changes that are allowed in general in the order graph are either (a) a swap of ranking of two adjacent fluctuation strengths, or (b) removal of the weakest fluctuation or emergence of a new weakest fluctuation.
t e x i t s h a 1 _ b a s e 6 4 = " 6 K + A m s 7 A a o E L a w 9 n y c g N x 6 B 4 k w t e x i t s h a 1 _ b a s e 6 4 = " 6 K + A m s 7 A a o E L a w 9 n y c g N x 6 B 4 k w a j t g T V 6 3 a 7 q + W u N T D O 4 O K + N I 0 i P d U o s e 6 I 5 e 6 O P X W k 2 / R t t L g 2 e 1 o x X 2 X v R 4 K v f + r 6 r K s 4 e D L 9 W f n j 2 U s O x 7 1 d m 7 7 T P t W 2 g d f f 3 w v J V b y S a b c 3 R F r + z / k p 7 o n m 9 g 1 t + 0 6 4 z I a j t g T V 6 3 a 7 q + W u N T D O 4 O K + N I 0 i P d U o s e 6 I 5 e 6 O P X W k 2 / R t t L g 2 e 1 o x X 2 X v R 4 K v f + r 6 r K s 4 e D L 9 W f n j 2 U s O x 7 1 d m 7 7 T P t W 2 g d f f 3 w v J V b y S a b c 3 R F r + z / k p 7 o n m 9 g 1 t + 0 6 4 z I a j t g T V 6 3 a 7 q + W u N T D O 4 O K + N I 0 i P d U o s e 6 I 5 e 6 O P X W k 2 / R t t L g 2 e 1 o x X 2 X v R 4 K v f + r 6 r K s 4 e D L 9 W f n j 2 U s O x 7 1 d m 7 7 T P t W 2 g d f f 3 w v J V b y S a b c 3 R F r + z / k p 7 o n m 9 g 1 t + 0 6 4 z I a j t g T V 6 3 a 7 q + W u N T D O 4 O K + N I 0 i P d U o s e 6 I 5 e 6 O P X W k 2 / R t t L g 2 e 1 o x X 2 X v R 4 K v f + r 6 r K s 4 e D L 9 W f n j 2 U s O x 7 1 d m 7 7 T P t W 2 g d f f 3 w v J V b y S a b c 3 R F r + z / k p 7 o n m 9 g 1 t + 0 6 4 z I Fig.17shows the distribution P (X, Y ), where in panels (a) and (b) we can observe the coexistence of thermal paramagnetic and VBS states and distinguish the type of VBS the system transitions into from the number of peaks.The fact that Fig.17(a) has five peaks, similar to Fig.12(a), shows that the phase coexistence at Q = 0.9 is indeed between the paramagnetic phase (center peak) and the PVBS phase (four surrounding peaks).In contrast, in Fig.17(b) there are nine peaks in total at Q = 0.95, implying a phase coexistence between the paramagnetic phase and the eight-fold degenerate AVBS phase.We can also see that the quantum fluctuation induced by decreasing Q is connecting the domains into four pairs in the order parameter space.Because of this, the four peaks of Fig.17(a) have a rather extended shape in the radial direction as opposed to Fig. 12(a).Fig. 17(c) shows emergent U(1) symmetry at very low temperature (essentially in the ground state; T = 1/L) in the vicinity of the AFM-PVBS transition, just as in Fig. 5(c) and Fig.17shows the distribution P (X, Y ), where in panels (a) and (b) we can observe the coexistence of thermal paramagnetic and VBS states and distinguish the type of VBS the system transitions into from the number of peaks.The fact that Fig.17(a) has five peaks, similar to Fig.12(a), shows that the phase coexistence at Q = 0.9 is indeed between the paramagnetic phase (center peak) and the PVBS phase (four surrounding peaks).In contrast, in Fig.17(b) there are nine peaks in total at Q = 0.95, implying a phase coexistence between the paramagnetic phase and the eight-fold degenerate AVBS phase.We can also see that the quantum fluctuation induced by decreasing Q is connecting the domains into four pairs in the order parameter space.Because of this, the four peaks of Fig.17(a) have a rather extended shape in the radial direction as opposed to Fig. 12(a).Fig. 17(c) shows emergent U(1) symmetry at very low temperature (essentially in the ground state; T = 1/L) in the vicinity of the AFM-PVBS transition, just as in Fig. 5(c) and

FIG. 19 .
FIG. 19.Depictions of spinons forming at the nexus of domain walls in a PVBS state in two possible ways; with dimers perpendicular and parallel to the walls in (a) and (b), respectively.Extreme cases of the thinnest domain walls are shown to the left, and to the right the walls are widened by additional dimers.Dimer pairs can dynamically convert into resonating plaquettes, and vice versa, thus allowing the domain wall to fluctuate.We have here drawn plaquettes closest to the unpaired spin in order to clearly show that all four phases meet at the spinon.Fluctuations leading to changes in the spinon location are illustrated in Fig. 20.In both (a) and (b), the relative relations of the domain colors are the same, e.g., red and blue differ by a y direction shift in the PVBS pattern.

FIG. 20 .
FIG. 20.Illustration of quantum fluctuations allowing a spinon to move in two different cases, corresponding to (a) and (b) in Fig.19.Colors were changed to reflect the final domain wall configuration, while the grayed singlets reflect the initial domain wall position.

FIG. 21 .
FIG. 21.Binder cumulants of the AFM and and CVBS order parameters of the the columnar J-Q6 model.The inset shows one of the Q6-terms in the model; all translations of this operator and its π/2 rotated version are included in the Hamiltonian, defined in analogy with Eq. (1).