Emergent fractons and algebraic quantum liquid from plaquette melting transitions

Paramagnetic spin systems with spontaneously broken spatial symmetries, such as valence bond solid (VBS) phases, can host topological defects carrying non-trivial quantum numbers, which enables the paradigm of deconfined quantum criticality. In this work, we study the properties of topological defects in valence plaquette solid (VPS) phases on square and cubic lattices. We show that the defects of the VPS order parameter, in addition to possessing non-trivial quantum numbers, have fracton mobility constraints deep in the VPS phase, which has been overlooked previously. The spinon inside a single vortex cannot move freely in any direction, while a dipolar pair of vortices with spinon pairs can only move perpendicular to its dipole moment. These mobility constraints, while they persist, can potentially inhibit the condensation of vortices and preclude a continuous transition from the VPS to the Néel antiferromagnet. Instead, the VPS melting transition can be driven by proliferation of spinon dipoles. For example, we argue that a 2d VPS can melt into a stable gapless phase in the form of an algebraic bond liquid with algebraic correlations and long range entanglement. Such a bond liquid phase yields a concrete example of the elusive 2d Bose metal with symmetry fractionalization. We also study 3d valence plaquette and valence cube ordered phase, and demonstrate that the topological defects therein also have fractonic dynamics. Possible nearby phases after melting the valence plaquettes or cubes are also discussed.

The search for exotic quantum phases and their transitions has gained a wide audience in recent years, due to a range of unusual properties which cannot be understood within the Landau paradigm of symmetry breaking. For example, topologically ordered phases, which do not possess a local order parameter, can exhibit topologically protected ground state degeneracies and deconfined fractionalized excitations. Furthermore, even for seemingly conventional symmetry-breaking phases, it is possible to have Landau-forbidden phase transitions with similar unusual phenomenology. The most well-known example of such a transition is between a Néel antiferromagnet and a valence bond solid (VBS). These phases break two different types of symmetries (spin rotation and lattice symmetries, respectively), so Landau theory would predict that a generic transition between these phases is either first order, or possesses an intermediate regime where the two orders coexist. In contrast, a more in-depth analysis reveals the existence of a generic second-order transition between these phases [1]. While the Néel and VBS phases have fairly simple phenomenology, characterized by Landau order parameters, the unusual critical point between them hosts deconfined fractionalized excitations, earning it the name of a "deconfined quantum critical point." This phenomenon of deconfined quantum criticality has now been intensely studied in a variety of physical systems [2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18].
Intuition for this unusual transition can be gained by arXiv:1908.08540v1 [cond-mat.str-el] 22 Aug 2019 studying the defects of each type of order, which carry quantum numbers associated with the other symmetry. For example, a vortex of VBS order naturally carries spin-1/2, as we review in Section II. Destruction of the VBS phase via condensation of these vortices then necessarily breaks spin-rotation symmetry, leading to an antiferromagnetic Néel phase. This point of view leads to a general picture for deconfined quantum criticality, with implications well beyond the Néel-VBS transition. While the VBS phase is the most commonly studied spatially ordered phase with spin quantum numbers associated with its topological defects, it is certainly not the only one. Another order of this type, often overlooked in discussions of deconfined quantum criticality, is a valence plaquette solid (VPS), in which spins are entangled in clusters of fours, instead of the pairwise entanglement associated with a VBS phase [19][20][21][22][23][24][25]. Both VPS and VBS order, as candidates of paramagnetic crystals, have been oberserved in a wide variety of materials [26,27]. Like a VBS defect, vortices of this plaquette order also carry spin-1/2, so we expect that destruction of the VPS phase via vortex proliferation will lead to a phase with broken spin-rotation symmetry, such as a Néel antiferromagnet. However, as we demonstrate in detail, this simple picture of vortex proliferation is complicated by the fact that vortices of VPS order behave like fracton quasiparticles, with their characteristic mobility restrictions, leading to important consequences for phase transitions out of the VPS phase.
The tensor gauge theory description allows us to make an immediate connection between plaquette order and the physics of fractons. As we review in Section II, a conventional valence bond solid can be mapped onto the confined phase of a vector gauge theory, with vortices behaving as linearly confined charges. In similar fashion, we show in Section III that a VPS phase can be mapped onto the confined phase of a symmetric tensor gauge theory, with vortices of the plaquette order acting as the fractons. Even though these fractons are confined within the VPS phase (corresponding to the large energy cost associated with vortices), their mobility restrictions still can have important consequences for phase transitions driven by vortex proliferation. For a valence bond solid, the spin-carrying vortices become gapless deconfined excitations at a quantum critical point, then subsequently condense to drive the system into a Néel antiferromagnet. For a VPS phase, on the other hand, the vortices are immobile fractons, for as long as the description in terms of plaquettes remains valid. (The mobility restrictions could break down in a regime where a plaquette can easily break down into a pair of dimers.) Assuming we remain in a regime well-described in terms of plaquettes, potential phase transitions and quantum critical points can be strongly impacted by the fractonic nature of the vortices. Even if the fractons become deconfined at a quantum critical point, their mobility restriction serves as an impediment to direct condensation. Fracton systems therefore have a tendency to first exhibit a condensation transition of mobile dipoles, which relaxes the mobility restrictions and allows a subsequent fracton condensation transition. This two-stage nature of fracton condensation transitions is dramatically realized in the analysis of two-dimensional quantum melting, which predicts that two-dimensional crystals must pass through a hexatic phase before fully melting to an isotropic liquid [34,37,38].
We therefore conclude that VPS order can generically host an intermediate phase in which only dipoles of vortices are condensed, while individual vortices remain gapped. The precise nature of this intermediate phase depends on the microscopic details governing the interaction of two vortices within a dipole pair. In certain cases, the intermediate phase may be a simple bondordered phase, such as a valence bond solid. As another illustrative example, we argue that a transition between two-dimensional VPS and Néel phases can feature a stable intermediate gapless phase in the form of an algebraic bond liquid [71][72][73][74] with quasi-long-range order between dipoles formed by spinon pairs. Such a featureless gapless liquid carries many features akin to the concept of the 'Bose metal,' including power law correlations and zero-energy nodal lines. In particular, its thermodynamic and entanglement properties have the hallmarks of a 2d Fermi liquid and hence can be regarded as the 'boson descendant' of a Fermi surface [75]. A gapless intermediate phase of this type is consistent with existing numerics on the VPS-Néel transition in 2d Heisenberg models [20]. In Section III, we review the stability of this algebraic bond liquid phase and study its transition with a VPS phase. We also describe some novel properties which may be used to detect an algebraic bond liquid in experiments or numerics, such as specific heat, structure factor, and entanglement entropy. The bond liquid exhibits novel characteristics, including T ln(T ) depen-dence for specific heat [72] and long-range entanglement, with entanglement entropy scaling as L ln(L) which exceeds the boundary law [76].
In Section IV, we extend these ideas to threedimensional cube-ordered and valence plaquette phases on a cubic lattice. In each case, the fundamental topological defects of the order behave as immobile fractons. In the case of cube order, even spinon dipoles are locked in place, with spinon quadrupoles behaving as onedimensional lineons. We also study different types of plaquette order in three dimensions, focusing on a strongly anisotropic VPS phase which can exhibit a continuous transition to a VBS phase. This phase transition features a reduced dimension allowing the physics of certain 3dimensional VPS melting transitions to be mapped onto the problem of a 2-dimensional VBS melting transition.

II. REVIEW OF THE VBS-NÉEL TRANSITION
Let us first recall the prominent deconfined quantum critical point between VBS and Néel phases on a 2D square lattice spin- 1 2 system. For an interacting spin-1 2 model on a square lattice, the Lieb-Schultz-Mattis theorem forbids any featureless gapped paramagnet ground state [77][78][79]. Two common symmetry breaking states are the Néel order, which spontaneously breaks spin rotation and lattice translation symmetries, and the valence bond solid (VBS) order, which breaks the lattice translation and C 4 rotation. In Ref. [1][2][3], it was proposed that, despite their distinct symmetry breaking patterns, the VBS and Néel phases in a square lattice spin- 1 2 system can be connected by a continuous phase transition which is dubbed as a deconfined quantum critical point. The key ingredient for such a Landau-forbidden phase transition between two symmetry broken phases lies in the fact that the topological defect of the order parameter for one broken symmetry carries a nontrivial quantum number of the other symmetry group. In order to restore a broken symmetry, one can consider the condensation of topological defects of the order parameter. However, this condensate would drive the system into another broken symmetry state simultaneously due to the charge carried by the defect. In the VBS-Néel example, the VBS phase spontaneously breaks the C 4 lattice rotations, which has four degenerate ground state patterns related by the C 4 rotations, as depicted in Fig.2. The VBS order can be characterized as a complex scalar with four-fold anisotropy. To quantum disorder the VBS order, one should proliferate the Z 4 vortices. Since the VBS vortex carries an unpaired spin, as shown pictorially in Fig. 2, consequently the vortex condensation breaks spin rotation symmetry. To formulate the VBS-Néel transition, it is useful to note that the VBS phase can be mapped to the confined phase of a compact U (1) gauge theory. The dimer coverage on the edge of the square lattice can be mapped to an electric field, where D i (r) refers to the dimer coverage on the edge adjacent to site r along i-direction, which only takes values 0 or 1. The index i r = 0 or 1 for r ∈ A or B sublattice respectively. In the VBS phase, each spin can be paired with one and only one of its neighboring spins to form spin singlets. This translates to the local Gauss's law for the electric fields, where q(x) is the number of unpaired spinons at site r.
(It is worth noting that the 'spinon quantum number' corresponds to the gauge charge of the emergent U(1) gauge field carried by the fractionalized spinon, which should be distinguished from the S z charge.) A free spinon charge appears if no dimer is touching the site. Breaking a dimer can create a pair of spinons. Further, let us define the gauge connection, A i , which is the conjugate variable to E i , namely [A i (x), E j (y)] = i 2π δ ij δ xy . In the spin picture, A i operator breaks (creates) a dimer on link-i if the link does (not) have an existing dimer. The magnetic flux operator B = ∇×A flips local dimer orientation on a plaquette, as in Fig. 1. Thus, the low-energy physics of the VBS phase is characterized by a compact U (1) gauge theory with background charges whose Hamiltonian is given by the following, The pure (2+1)d compact U (1) gauge theory, due to the instanton events, is always in its confined phase at low energy, where the spinon excitations experience a linearly confining potential. The confined phase is mapped to the VBS phase where the VBS vortices have linear confinement due to the energy cost of the domain walls connecting the spinons.
Consider the phase transition out of the VBS state by proliferating the VBS vortices. The quantum critical point can be described by the following field theory, where z is a CP 1 field that captures the spinon degree of freedom inside the vortex core. A µ is the emergent gauge field described previously. Since the spinon is charged under the emergent gauge field as in Eq. 2, they are minimally coupled to the gauge fields. The hopping of a spinon thus requires the change of dimer configurations along the path. In this theory, we implicitly include the 4fold monopole creation and annihilation operators which correspond to the possible instanton events that respect the lattice C 4 rotation symmetry [1][2][3]. When the spinon is gapped, namely r > 0, the gauge theory is confined due to monopole proliferation and the system is in the VBS phase. At the critical point, r = 0, the spinon field becomes massless and there is evidence that the compact U (1) gauge field dynamically becomes deconfined at the fixed point [1], i.e. the 4-fold instanton events are irrelevant under renormalization group flow. As r decreases below 0, the VBS vortices/spinons condense which restores the C 4 rotation and spontaneously breaks the spin rotation symmetry.

III. PLAQUETTE PARAMAGNET IN 2D
Apart from columnar valence bond order, another widely observed paramagnetic crystalline phase is the VPS (valence plaquette solid) state which breaks C 4 symmetry and lattice translation T x , T y for both directions. Such a plaquette paramagnet has been found and fabricated in frustrated magnets and AMO systems [21,26,80]. There are various types of VPS wave functions which respect the same symmetries. For example, for a spin-1/2 system, each valence plaquette represents a symmetric combination of vertical and horizontal dimer pairs on the same plaquette. For an SU (4) spin system with a fundamental representation on each site, one can have plaquettes in an entangled state of four SU (4) spins. Regardless of the microscopic configuration inside the valence plaquette, each spin participates in only one of the four plaquette clusters adjacent to the site. The VPS order enlarges the unit cell into four plaquettes, so there are four distinct VPS patterns related by site-centered C 4 rotation, as shown in Fig. 3.

A. Defect structure
To restore the spatial symmetry, it is essential to proliferate the C 4 symmetry defect. For the plaquette crystalline phase, one can define the four distinct plaquette patterns as a Z 4 boson. The C 4 symmetry defect formed by the vortex configuration of the order parameter carries a spinon. Naively, one expects the vortex condensate will restore the crystalline symmetry and concurrently develops magnetic order. However, in order to drive the defect condensate, vortices need to be able to fluctuate in spacetime. As opposed to the VBS phase, where a spinon in the background of dimers can hop among sites by reconstructing local valence bond configuration, a spinon in the background of plaquette order is frozen -cannot move away from the original vortex center without breaking additional plaquettes, as depicted in Fig. 4. In contrast, a pair of spinons living on the link between adjacent sites can hop along the stripe perpendicular to that link without breaking additional plaquettes, as depicted in Fig. 5. Such a spinon pair, which we refer to as a spinon dipole, is a 1d subdimensional particle which moves transversely to the dipole's orientation.
Based on these observations, the topological defect of the plaquette order displays restricted motion which exactly resembles the behavior of fractons. In fracton phases of matter, the fundamental deconfined quasiparticle excitation is immobile due to a finite energy barrier associated with the creation of additional excitations. Meanwhile, a pair of fractons (dipole) has a certain degree of mobility, though it too is often restricted to motion only within a submanifold, such as a line, plane, or fractal.
To make the connection between VPS defects and fractons precise, we introduce a higher rank gauge theory description for the valence plaquette order on a square lattice. By analogy to the VBS order, the plaquette order can be mapped to a rank-2 symmetric tensor electric field defined at the center of each square as the following.
where P = 1 (0) corresponds to the valence plaquette occupancy (vacancy) on each square. The index i r is the same as defined before. As opposed to the VBS state, where dimers can have two orientations corresponding to E x and E y , the plaquette electric field is a singlecomponent field, effectively a scalar. We can also define a conjugate variable A xy , satisfying [A xy (r), E xy (r )] = i 2π δ r,r . The operator e ±iAxy creates/annihilates a valence plaquette. As each spin on the site is only entangled with one of the four adjacent plaquette clusters, one can define a Gauss's law for the rank-2 electric field as, where q(r) is the number of unpaired spinon at site r. As long as there is one plaquette adjacent to a site, there is no free spinon on that site. If plaquettes are absent from all four squares surrounding the site, then there exists a free spinon charge at the center. This Gauss's law is precisely the two-dimensional version of the Gauss's law seen in the fracton phase of matter described by a hollow rank-2 symmetric tensor gauge theory [44,65,66]. Due to the particular double derivative in Eq. 7, the spinon number is conserved on each row and column of the system, so the theory respects an emergent subsystem U (1) symmetry: A similar equation holds in the y-direction. Again, we emphasize that q is the spinon charge quantum number corresponding to the emergent U (1) symmetry, which should be distinguished from the total S z charge.
In particular, only the spinon charge is conserved on row/columns while the S z charge is conserved globally. Due to the emergent subsystem symmetry, single spinon motion is prohibited. However, a pair of spinons, which we refer to as a dipole, can hop only along the stripe perpendicular to its orientation. Based on the Gauss's law in Eq. 7, the low-energy sector is invariant under the following gauge transformation, for any function α with arbitrary spatial dependence.
Since there is only one component of the gauge field, there is no local gauge invariant operator representing a flux degree of freedom. The absence of a local flux operator indicates that there is no resonant process that can flip plaquette configurations locally into each other. However, one can define global flux operators, These global flux operators shift a row (or column) of the valence plaquette configuration by one unit cell, as in Fig. 6.

B. Bond ordered phases
Due to the restricted mobility of the spinon corresponding to the VPS defect, it is difficult for the spinon to condense directly. We therefore do not expect a continuous transition from the VPS phase to a simple Néel antiferromagnet. Based on this observation, our next goal is to figure out the possible phases near the plaquette melting transition, driven by condensation of dipolar pairs of vortices. Since these vortex pairs live along the links of the lattice, it is natural to expect that the resulting phase will be well-described in terms of valence bonds, instead of valence plaquettes. The precise nature of this intermediate phase will depend on the microscopic structure of the dipolar pair, which is dictated by the details of the underlying Hamiltonian.
In the simplest case, the spinon dipole may form an SU (2) singlet. In this situation, the result of dipole condensation will be a type of valence bond solid. Even in this case, there are several types of possible VBS states which can result, depending on how precisely the dipoles condense. For example, condensation of ydirected dipoles will naturally lead to a VBS state with all valence bonds aligned with the y-direction. This physics is borne out by studies on two-dimensional quantum dimer models on a square lattice, which can host a continuous transition between plaquette order and a staggered VBS phase, whose critical point is described by the quantum Lifshitz theory with z = 2 dynamical exponent [81].
Beyond such VBS phases, obtained by condensing singlet spinon dipole, another possibility is that the two spinons of a dipole tend to align their spins in a triplet state. In this case, since the dipoles are carrying a net spin, the resulting condensed phase breaks spin rotation symmetry. However, this state will not be the simple Néel phase, but rather will have antiferromagnetic order coexisting with a form of bond order. An example of a system with this type of order is depicted in Figure 7. Further condensation of defects of the bond order can then drive this phase into a simple Néel antiferromagnet.

C. Algebraic bond liquid phase
Beyond analytical exploration of transitions to bondordered phases, there are also various numerical simulations of frustrated spin models and hardcore boson models which appear to reveal a possible symmetric gapless phase nearby the VPS phase. This indicates that a plaquette ordered system may melt into a more exotic phase under certain conditions. In this section, we propose a stable algebraic bond liquid phase with power-law correlations which can potentially emerges near the VPS phase.
In our previous discussion, we have elucidated that the spinon pair along a link can hop and fluctuate along the 1d stripe perpendicular to the dipole orientation. During the melting of the VPS paramagnet, the dipole proliferates and fluctuates along the stripe with quasi-onedimensional dispersion. Because of the restricted dimension, the spinon dipole cannot establish long range order. Instead, the dipoles form a state with quasi-long range order. The system form an algebraic bond liquid.
To elucidate the behavior of the emergent spinons during the plaquette melting transition, we implement the Schwinger boson representation, with the constraint z † 1 z 1 + z † 2 z 2 = 1. The two component field (z 1 , z 2 ) T are the bosonic spinons. The total charge density is fixed in the original microscopic model. The relative density 1 2 (z † 1 z 1 − z † 2 z 2 ) corresponds to the S z quantum number. Both z 1 and z 2 couple to the same emergent U (1) gauge field, while each of them also carries ±1/2 charge for S z . The magnon excitation S + = z † 1 z 2 corresponds to the exciton pair of Schwinger bosons with opposite flavors. The subsystem charge conservation law in Eq. 7 indicates the gauge charge of the Schwinger boson is conserved on each line. We can approximately express the Schwinger boson operator in terms of a U (1) rotor field as z † a ∼ e iθa , writing the effective Hamiltonian of the bosons in the following form, where the u term is to implement the onsite constraint n 1 + n 2 = 1, which can be released when coarse graining in the continuum limit. Note that the Hamiltonian contains no usual kinetic term for the phase variables such as (∂ i θ a ) 2 , as a single spinon is immobile in any direction. Instead, the leading order kinetic term ∂ x ∂ y θ corresponds to an xdirected dipole hopping along the y direction, or vice versa. Such dipole hopping requires reconstruction of the valence plaquette pattern along the hopping path, so the dipole current minimally couples with the gauge field A xy . Since there is no local gauge flux for such a higher rank gauge field, the role of gauge fluctuation merely projects the Schwinger boson to the physical Hilbert space of n 1 + n 2 = 1.
In the easy-plane anisotropy limit where S z =0, we take n 1 = n 2 = 1 2 so each slave boson is at half-filling. This theory resembles the exciton bose liquid phase studied in Ref. [71][72][73][74], where the bosons interact via a ringexchange interaction which respects a subsystem U (1) symmetry. Such an exciton bose liquid constitutes a stable gapless phase with power-law correlation at fractional filling. When the u/t becomes small and the theory tends to have order in θ variables, we can implement a "spinwave" approximation, expanding the cosine terms, and obtain a Guassian theory, The action is invariant under the following transformation, which are the remnant of the subsystem symmetry of the VPS phase. We can decomposed the two branches as We here rescale space and time to set the two coefficients equal. There is then just one remaining dimensionless paraxmeter as K = t/u. Such a quadratic Lagrangian, in which all terms involve derivatives of the fields, describes a scale invariant phase at long length scales. In a sense it can be viewed as a "fixed point" Lagrangian. The legitimacy of this approximation and the relevance of compactness will be discussed in detail. The θ + mode and the gauge field A xy gap out each other through an analog of the Higgs mechanism. Subsequently, only the θ − branch is physical and that is the degree of freedom we will consider here and after. We first scrutinize the confinement energy between spinons and dipoles in the VPS phase. Due to the subsystem symmetry, the spinons must be excited in quadrupole pairs with two living on the same row/column of the lattice. Following the duality and bosonization argument introduced in Ref. [71,72], we can map the theory in Eq. 16 to its dual representation. As the θ − fields are compact with the identification θ − = θ − + 2πZ, particular types of topological defects will be allowed, which can be most conveniently addressed by passing to a dual representation N, φ defined on the plaquette centers [82].
N and φ are a pair of conjugate variables with φ being discrete-valued andN being compact with N ∈ [0, 2π].
The Gaussian part of the dual action is, In this dual picture, due to the discreteness of φ − , one can also add vertex operators such as cos(4π∂ i φ − ) in the effective theory. In the VPS phase, these vertex terms are relevant so the spinon and dipole are confined, on which we elaborate further below. In the presence of these vertex operators, the separation of four spinons costs a linear confinement energy proportional to the domain wall length.
We now demonstrate the stability of the algebraic bond liquid phase obtained via melting of the valence plaquette order. In the original representation in Eq. 16, the correlation between two charges cos(θ − (r)) cos(θ − (0)) vanishes at long-wave length due to the subsystem U (1) symmetry. The long-range correlation between Schwinger bosons implies the spinon cannot be excited in pair but instead are created in quartet forms. The leading order non-vanishing correlation functions are between two dipole operators living on bonds, Notice that the dipole-i correlation function is only nonzero when they are at the same row transverse to the dipole orientation. Thus, the dipoles effectively behave as a 1d Luttinger liquid with restricted motion and algebraic correlation on each stripe. This quasi-onedimensional behavior is crucial for the stability of the bond liquid phase. As the dipole displays 1d motion within the same stripe, the quantum fluctuation forbids any dipole condensation with off-diagonal long-range order as a consequence of the Mermin-Wagner theorem.
Instead, there appears a quasi-long range order between dipoles akin to the Kosterlitz-Thouless transition. The quasi-long range order between dipoles corresponds to a spinon-quadrupole excitation with four spins living on the corner of a thin stripe.
To explicate the stability of the bond liquid phase, we need to consider the compactness of the θ − term. In the dual representation in Eq. 18, the φ − is chosen to be half-integers so we consider the vertex operators V = cos(4π∂ i φ − ), V = cos(2π∂ 2 i φ − ). When K is small, these vertex operators turn out to be relevant, so the system is driven into a Mott phase which essentially breaks the crystalline symmetry. In Ref. [21,72], it was shown that there is a finite region for K > K c where all vertex operators are irrelevant, so the algebraic bond liquid is stable. Hence, there could potentially appear an algebra bond liquid phase near the VPS phase.
When K < K c , the vertex operator proliferates and thus creates a spin gap. According to the Lieb-Schultz-Mattis theorem, such a gapped spin-1/2 model does not support a featureless Mott phase, so the ground state must break lattice symmetry, which corresponds to the VBS order or VPS order depending on the microscopic symmetry of the vertex operator. For V i = cos(4π∂ i φ − ), the vertex operator creates a kink for the dipole-i along the transverse stripe. If both V x , V y proliferates, the system breaks translation on both directions and thus falls into the stripe order or a plaquette order, depending on the sign of V i [71]. When ∂ i φ − = 0, the system becomes prone to plaquette order as which exactly agrees with our VPS picture. Based on this observation, we conclude that the proliferation/suppression of V i drives the transition between a VPS phase and the algebraic bond liquid. When ∂ i φ − = π/2, the system favors a stripe order as Finally, we mention that, in Ref. [73,74], the authors present a microscopic boson model with various ringexchange terms which supports a phase transition between plaquette order and the bond liquid phase, with a possible competing order toward a charge density wave. This opens a search for the unconventional plaquette melting quantum phase transition in concrete spin models.

D. Signatures of the algebraic bond liquid phase
The algebraic bond liquid phase is engendered by the quasi-one-dimensional dipole fluctuations. Due to the restricted motion of spinons, a single spinon excitation is gapped and the spinon correlation e iθa(r) e iθa(r ) vanishes at long wavelength due to the subsystem U (1) symmetry. The bond correlator, denoting the dipole correlation of spinon pairs, has quasi-long-range order as, It is worth mentioning that the dipole correlation is anisotropic and only displays algebraic order along the transverse direction, as a consequence of the subdimensional behavior of the dipoles. In particular, the quasione-dimensional motion of the dipole, originating from the subsystem charge conservation law, is crucial for the stability of the bond liquid phase. While spontaneous U (1) symmetry breaking is generally expected in 2d quantum systems, a subsystem U (1) symmetry consists of independent symmetry operations acting on an extensively large set of 1-dimensional lines. The quantum fluctuations thereby suppress dipole long range order, so the Mermin-Wagner theorem still applies. As a result, the bond correlations decay as a power law due to the absence of spontaneously broken subsystem U (1) symmetry in 2d. However, as the spinon is the emergent fractionalized degree of freedom, its four point correlation cannot be measured directly. Instead, we can measure the magnon pair correlator, Which also renders an algebraic correlation. The above result is based on the one-loop correction where the magnon correlation is simply the product of two Schwinger boson correlators. Including higher loop corrections can potentially change the exponent of the algebraic correlation.
If we go back to the square lattice structure with unit length a, the spinon pair between links can only hop along the transverse direction with even lattice units 2a (we will take a = 1 henceforward). Thus the periodicity of the unit cell is doubled and the Brillouin zone sits in the region k i ∈ [−π/2, π/2] with dispersion as, The low-energy dispersion displays a zero-energy nodal line which qualitatively changes the IR behavior, including transport features and entanglement. Such an excitation spectrum with nodal lines can be measured in terms of the static structure factor for the S z correlator. In our slave boson representation, the S z number corresponds to the imbalanced charge density between the two Schwinger bosons as S z ∼ n 1 − n 2 , so the static structure factor for the S z correlator can be written in terms of the slave boson correlator as, S zz (k) ∼ (n 1 (k) − n 2 (k))(n 1 (−k) − n 2 (−k)) = (n 1 (k)n 1 (−k) + n 2 (k)n 2 (−k)) ∼ | sin(k x ) sin(k y )| This structure factor can be measured using conventional experimental techniques, such as inelastic neutron scattering and electron spin resonance.
As opposed to the usual bosonic superfluid theory, whose low energy mode condenses at zero momentum, the bond liquid phase contains two nodal lines along the two axes. At each fixed momentum k x , the dispersion resembles a relativistic 1d theory, E = v f k y . Such a bond liquid phase, with low-energy modes containing nodal lines, is termed as a "Bose surface", in analogy with the 2d Fermi liquid with a Fermi surface. For both Fermi surfaces and Bose surfaces, each patch with fixed transverse momentum carries a linearly dispersing 1d mode. Due to the existence of the Bose surface, the low energy transport behavior of the bond liquid phase is qualitatively different from the usual photon gas or weakly interacting bosons. For each fixed nonzero momentum k x , the low energy modes display linear dispersion with respect to k y , akin to 1d relativistic bosons. In particular, due to the nodal lines at k x , k y = 0 with a subextensive number of quasi-1d modes, the specific heat at low temperature scales as [72], which is similar to the marginal Fermi liquid theory in 2d.
The excitation in Eq. 24 corresponds to the dipole excitation containing two spinons between bonds. However, as the spinon is conserved on each line, the dipole excitations must be created in pairs. To be concise, dipoles are created in pairs on the same stripe as a consequence of subsystem symmetry so the spinon appears in quadrupole form. This phenomenon is in close analog to the spinon excitations in 1d spin chains where the magnon excitations can fractionalized into two spinon excitations. Consequently, the spin spectrum function covers a broad continuum whose upper and lower limit corresponds to the parallel and anti-parallel motion of the two spinons. Such continuum in the spectral function distinguishes the spinon excitations with the regular magnons with sharp dispersion.
To seek the collective mode of the dipoles, we calculate the spectral function for B † x B x . The B x = S − (r)S − (r + e x ) operator creates a pair of magnons between an x-link. For a non-fractionalized bond liquid [71,72], whose dipoles are composed of magnon pairs with integer S z charge on each row, the spectral function has contributions from the magnon pair excitations with a sharp dispersion relation. In our algebraic bond liquid state, the dipole excitations, created in pairs, carry two spinons between the link with half S z charge on each row/column. Such a collective excitation can be interpreted as two dipoles on the same stripe moving along the transverse direction with independent dynamics.
In our slave boson theory, the bond operator B x = S − (r)S − (r + e x ) = e −i∂x(θ1−θ2) creates the x-dipole excitation for both z 1 , z 2 slave particles, each carrying half S z charge. As these slave dipole pairs are deconfined excitations in the bond liquid phase, each propagates along the stripe with independent motion and the collective excitation corresponds to the combination of the two. To reach such a collective excitation among spinonpairs, we calculate the static structure factor for the B † x B x = e i∂x(θ1−θ2) e −i∂x(θ1−θ2) correlator, = dkdω e ik 2 G1(k,ω) e i(k+Q) 2 G2(k+Q,ω+Ω) To extract the collective excitation spectrum, we expand the correlator as, The poles in each Π m,n (Q, Ω) correspond to a collective mode. We start with the leading order expansion Π 1,1 (Q, Ω), This propagator renders a series of poles as Ω(Q) = E(Q − k) + E(k). The energy spectrum at fixed momentum Q covers a broader range depending on the choices of k. This continuum spectrum can be understood as the two-dipole excitation with dispersion E(Q − k) and E(k). The total momentum is fixed as Q with Q − 2k being relative momentum due to the independent motion of the dipole pair. In particular, if we fix a slice of the momentum space by taking Q y = 0, π, the energy spectrum has an upper and lower limit as, The upper or lower limit of the spectrum denotes the parallel or opposite motion for the two dipoles. This collective excitation fills the continuum region between these two limits, as depicted in Fig. 8, which resembles the spinon spectrum in 1d anti-ferromagnetic spin chains. Actually, when fixing Q y , we focus on the collective motions of y-dipoles along the x-direction. As the dipole's motion is restricted along the x-stripe, the collective mode is attributed by the dispersion of two dipoles along the stripe [83], which can be regarded as the descendent of the spinon in 1d spin chains.
Another signature for the algebraic bond liquid is the violation of the area law for entanglement entropy, arising as a consequence of the bose surface. In Ref. [84], it was demonstrated that if we equally bipartition the bond liquid state on a lattice along the x-axis, the entanglement entropy scales as L x ln(L y ), which resembles the entanglement entropy for the 2d Fermi surface [85][86][87]. Such a violation of the area law can be roughly understood by dividing the Bose surface into small patches over which the surface looks approximately flat, such that each patch can be regarded as a one-dimensional relativistic boson whose entanglement entropy scales as ln(L). Summing over the contributions from all patches, the total entanglement entropy should behave as L ln(L). Such longrange entanglement is smoking-gun evidence for the Bose surface which can be detected in numerical simulations.

E. SU (3) Plaquette defect on triangle lattice
Our previous discussion based on plaquette order on square lattices can be straightforwardly applied to other bipartite lattices, such as the honeycomb lattice. In particular, the spinon living in the C 3 plaquette defect obeys a conservation law on each θ = 2π/3 line, so the spinon is also a fracton. A spinon pair, regarded as a dipole, can hop on transverse zigzag stripes as a 1d subdimensional particle.
Finally, we would like to bring attention to the SU (3) plaquette order on a triangular lattice with possible fractal dynamics. A typical plaquette order contains SU (3) singlets between three SU (3) spins living on the three sites of left-oriented triangles, as illustrated in Fig. 9. As the SU (3) valence plaquette only lives on the leftoriented triangles, the system breaks C 3 rotation symmetry and the ground state contains three VPS configurations related by C 3 rotation.
A typical topological defect renders the breaking of SU (3) singlet on one triangle which creates three dangling SU (3) spins. It is obvious that one can not move or separate a single SU (3) spin out of the defect due to fractal conservation laws. In particular, if we define

IV. FRACTONS IN 3D CUBE ORDERED AND VALENCE PLAQUETTE PHASES
In this section, we extend the scope of our analysis to a 3d cubic lattice. We consider properties of the topological defects in 3d valence cube solid (VCS) and valence plaquette solid (VPS) phases. We show that the topological defects of these orders are generically fractons. Then we consider possible outcome of the melting transitions of these orders indicated by the fracton dynamics.

A. 3D valence cube solid (VCS) order
A natural generalization of the 2d plaquette order in 3d is the cube order. The cube order for an SU (2) spin-1/2 model on a cubic lattice is a state which has resonating clusters of 8 spins on every other cube in the lattice. One can also imagine similar state for models with larger spin symmetry such as SU (4) and SU (8). We argue that topological defects of these cube ordered phases are emergent fractons just as in the 2d plaquette order. We also propose possible neighboring phases assuming the fractonic dynamics persists across the phase transition. In the following, we restrict ourselves to the case of spin-1/2 systems on cubic lattices.
The fractonic nature of the topological defects becomes clear after mapping the cube order to a rank-3 tensor gauge theory. Since cubic lattice is a bi-partite lattice, on each cube, one can define an "electric" field E xyz (r) = (−1) ir C(r), where i r =odd/even if r ∈ A/B sublattice, and C(r) = 1 (or 0) denotes that there is (or not) a resonating cluster on cube r [88]. E xyz furnishes a rank-3 hollow gauge theory in 3d. In the low energy Hilbert space, the electric field satisfies the following Gauss's law around a site on the lattice, where q = 0/1 denotes the number of free spinon at site r, and all the derivatives should be treated as lattice derivatives. The q = 1 state corresponds to a topological defect of the cube order, which maps to a point charge of the rank-3 gauge field. One can introduce the conjugate field A xyz (r) for E xyz (r), namely [A xyz (r), E xyz (r )] = i 2π δ r,r . The e ±iAxyz(r) is the creation/annihilation operator for E xyz (r) on cube r. The Gauss law in the low energy subspace in Eq. 32 implies the following gauge transformation for A xyz , With this gauge transformation, we can locally remove the gauge field. Therefore, there is no local "magnetic" flux in this rank-3 gauge theory. Physically, this means there is no local resonating process for the cube order. Namely any local adjustments of the cube order parameter inevitably break the Gauss's law constraint. Of course, there are still global flux operators, which can adjust the cube order either on a whole plane or along a straight line, analogous to the global flux operators appearing in the 2d case in Eq. 11. Since the topological defect, which traps a single spinon, appears as the matter field that couples to the rank-3 gauge field, it is a fracton that cannot move along any direction. Specifically, the Gauss's law of this theory implies conservation of all components of dipole and quadrupole moments, along with certain components of the octupole moment. The mobility of other point defects of the cube order are also easy to determine. From the physical picture of cube order, one can see that the spinon dipoles are also immobile. While the spinon quadrupole on a plane is movable along the normal direction, hence it is a linenon. Now we consider possible melting transitions if this fractonic constraint is kept all the way through. Since the spinon monopoles and dipoles have no mobility at all, it is hard to consider their condensation. The most probable way to drive the system out of the cube order is to proliferate the planar spinon quadrupole, which is movable along different lines. However, 1-dimension cannot host true long range order of continuous symmetry. Therefore, the resultant phase may be an algebraic spin liquid similar as the 2d case.
Let us write down the low energy field theory which encodes the coupling between the rank-3 gauge field and spinon matter field, where we have again used the CP 1 map as in Eq. 12 to fractionalize the spin at the topological defects. As before, the n 1 and n 2 are the number operators of the two bosonic spinons. We adopt a roton approximation for the spinons. Correspondingly, the θ 1 and θ 2 are the phases of the two bosons. This theory is invariant under the following symmetry transformation, Consider an easy plane limit, namely adding a term u a=1,2 (n a − 1 2 ) 2 to favorn 1 =n 2 = 1 2 . With a large t, the system tends to fall into an ordered state of θ 1 and θ 2 . In a spin-wave approximation, we can expand the cosine term and take the gaussian theory. In the resultant theory, the gauge fields and the ∂ x ∂ y ∂ z (θ 1 + θ 2 ) mode gap out each other through an analog of the Higgs mechanism. The only physical mode left is ∂ x ∂ y ∂ z (θ 1 − θ 2 ). Let us label θ = θ 1 − θ 2 , and n = n 1 − n 2 . The continuum theory now reads where we have rescaled spatial coordinates to simplify the Hamiltonian to this form with K ∼ t/u . Because of the symmetry in Eq. 34, there is no cosine terms for θ. The gaussian theory describes a 3d generalization of the 2d algebraic liquid phase, in which the planar quadrupole operators acquire algebraic correlations. For example, the quadrupole operators Q z = cos(∂ x ∂ y θ) have zero equal time correlation along x and y direction due to the symmetry in Eq. 34, however, power-law correlation along z-direction, Planar quadrupole operators along different planes have similar power law correlations. The directional powerlaw correlations of the planar quadrupoles are the remnant of their fracton dynamics in VCS phase.
In the gaussian theory, we have ignored the compactness of the θ variable. To justify the stability of the gaussian theory, one has to consider the relevancy of the vertex operators. To do this, it is most convenient to go to the dual description. We can have a duality map where φ and N are conjugate variables defined on the dual lattice sites. φ should take values in Z/2, while N ∈ [0, 2π). The dual hamiltonian reads In this continuous field theory, we have regarded φ fields as real number. However, since the φ variable actually takes values in Z/2, we have to include cosine terms, namely the vertex operators, to reinstall the integral constraint of the φ fields. There are various vertex operators that can appear in the theory, for example, V = cos(4πφ), V x = cos(4π∂ x φ), V xy = cos(4π∂ x ∂ y φ) and so on. The task is to determine the scaling dimensions of these vertex operators and to inspect if they are relevant at the gaussian fixed point. The gaussian theory has an emergent symmetry, which restricts the correlations of the vertex operators. The correlations between V at different points are zero, similarly for V x , due to the emergent symmetry. The V xy operator can have non-zero correlation functions only along z direction. The most relevant vertex operator has correlation function as the following, with η > 0 depending on the UV definition of the vertex operator. We can see that for large enough coupling K, the vertex operators in the theory can be all irrelevant, hence, the algebraic spin liquid is stable.

B. Tensor gauge theory and fractons in 3D VPS
Consider an SU (2) spin model on a cubic lattice. In analogy with the previously discussed 2d valence plaquette order, the 3d plaquette order contains 12 distinct plaquette order patterns, as there are 12 plaquettes adjacent to a single site. Such plaquette order breaks cubic symmetries and translations while maintaining the SU (2) spin rotation symmetry. We now show the valence plaquette order on the cubic lattice, similar to the 2d case, can be mapped to a hollow rank-2 gauge theory [21,62,66].
We denote the plaquette order on each i-j square as a tensor electric field, where the binary variable P ij = 1 (0) corresponds to the valence plaquette occupancy (vacancy) on each square. i r is again the sublattice index for site r. As the plaquette lives on three i-j planes, there are three components of the tensor electric field E xy , E zy , E xz , forming a symmetric rank-2 hollow (i.e. purely off-diagonal) gauge theory [21,62,66]. The gauge field A ij , as the conjugate variable of E ij satisfying the usual commutation relation [A ij , E ij ] = i 2π , is the operator which creates/annihilates a valence plaquette on each square and thus enables plaquette fluctuation. This is different from the 2d case as in 2d the plaquette order cannot fluctuate locally. Since each spin on a site is only entangled with one of the 12 adjacent plaquette clusters, one can write down the analogy of the Gauss's Law for the rank-2 gauge field as the following, where ∂ i should be regarded as the lattice derivative on the cubic lattice and the Gauss's law respects the cubic symmetry. This Gauss's law exactly resembles the hollow rank-2 symmetric gauge theory in 3D as the U (1) generalization of the X-cube fracton model [21,29,30,62,66]. Due to the particular double derivative in Eq. 42, the spinon is conserved on each i-j plane, so the theory respects a subsystem planar U (1) symmetry. Consequently, a fundamental topological defect, which carries a single spinon, is a fracton that is restricted from moving in any direction. In addition, the topological defects which host a pair of spinons along a link can hop within the 2d plane which is perpendicular to its dipole orientation. As opposed to the 2d VPS order, the plaquette configuration on the cubic lattice can fluctuate and resonate locally on each cube. These local fluctuations defined on each cube can be mapped to three types of magnetic flux operators in the gauge theory language, which are invariant under the following gauge transformation, Different from the vector U (1) gauge theory in 3d, the flux operators here are point-like excitations obeying the identity a B a = 0. In particular, the flux excitations are also fractons, known as lineons [21,30], which only move along straight lines. Such flux operators flip the plaquette configuration on one cube, as in Fig. 10, which generates resonant states between different flippable plaquette configurations. A typical hamiltonian for the pure compact rank-2 gauge theory can be written as the following, the electric fields are subject to the Gauss's law on every site. Due to the proliferation of topological defects, namely the 2π instanton tunneling which turns out to be relevant, the pure rank-2 gauge theory is generically in a confined phase with crystalline orders [21]. Next we consider matter fields and possible Higgs-like transitions in this rank-2 gauge theory. The matter field that couples to the gauge theory are the topological defects of the VPS order. (46) where again we have use the CP 1 formalism to represent the spinon trapped in the topological defects. The single spinon are fractons that cannot move, while a pair of spinon can move in the plane that is perpendicular to its dipole moment. Let us use this hierarchy of matter field mobilities to infer the possible nearby phases. Due to the restricted motion of the spinon, a direct condensation of spinons is inhibited. The leading ordering instability should be dipole condensation, where the spinon pair between links acquires coherence along a 2d plane. Such condensation restores the mobility of the spinon along the dipole direction and thus breaks subsystem symmetry. Depending on the microscopic Hamiltonian, the spinon pair condensate could engender a valence bond solid or liquid phase. In the following, we will consider a case with strong anisotropy where the melting transition leads the system to a valence bond solid state.
C. An anisotropic VPS We now consider a special limit of the valence plaquette solid with strong anisotropy along the z-direction, namely the valence plaquettes energetically favor to be on the xz, yz squares. In this case, one can imagine a columnar plaquette ordered ground state along xz or yz direction, which spontaneously breaks the C 4 rotation along z direction. In this anisotropic case, the low energy Hilbert space has E xy = 0. Thus the Gauss's law is reduced to, ∂ x ∂ z E xz (r) + ∂ y ∂ z E yz (r) = (−1) ir (1 − q s (r)). (47) Notice there are only two electric operators, corresponding to the plaquettes on xz and yz planes. A single spinon is still immobile, while a pair of spinons along a z-link can hop on the xy-plane as a 2d subdimensional particle. A spinon pair along an x or y link is conserved in each z-stripe, so they are restricted to fluctuate along z.
The valence plaquette configuration can fluctuate and resonate between the side faces of the cube. Such a local plaquette configuration resonance corresponds to the flux operator, which is invariant under the following gauge transformation, The effective field theory description for this anisotropic limit can be written as which is derived from Eq. 45 by enforcing E xy = 0 and confining the A xy gauge field component. As the plaquette order fluctuates, there can appear vortex configurations where the four plaquette patterns related by C z 4 symmetry meet at a point, forming a 2d vortex, as shown in Fig. 11. This C z 4 vortex defect, car- 11. A pair of spinon can hop on the 2D plane perpendicular to its polarization rying a spinon pair along a z-link, can fluctuate on the xy-plane. Similarly, there are spinon pairs along x and ylink. However, they can only fluctuate along z-direction due to the fracton constraint and anisotropy. Now let us consider a zero temperature quantum melting transition driven by condensation of these spinon pairs. We start from a particular VPS state, say all the plaquette is along xz plane. The melting transition will try to break the plaquettes. We can imagine a situation where the plaquette can only break into two z-direction dimers in the hilbert space. In this limit, we can imagine a VPS to VBS phase transition. Due to the fractonic nature of the dimers, two things need to be accomplished during the transition. a) In each doublelayer, the z-dimer proliferates and destroys the plaquette order within a layer. b) The layers establish coherence. If a) happens before b), then we have a transition that features dimensional decoupling. If b) happens first, we have a transition that has reduced effective dimensions.
Let us first analyze the phase transition within a double layer. Due to the anisotropy, the resonating plaquettes are formed between the two layers. If we view the system from the top, the VPS order actually is projected to a 2d VBS pattern. In this top view, the interlayer singlet pair is the vortex core of the VBS pattern. Therefore, for the double layer system, the VPS to VBS transition is mapped to a VBS melting transition. The key difference between this transition and the 2d DQCP is that the VBS vortices here do not have spinon degree of freedom. The low energy field theory for such a transition is the following which is similar to the theory of DQCP as in Eq. 5 except the matter field φ here is a single scalar. We still need to include the 4-fold instantons in the theory, which is allowed by translation symmetry. In the particle-vortex dual picture, the transition can also be described as a 3d XY transition with four-fold anisotropy. It is known that the 4-fold anisotropy is irrelevant at the 3d XY transition. Thus, the VPS to VBS transition in the double layer system is in the 3d XY universality class. Now we consider the coupling between these 2dimensional critical theories along the z direction. Notice that for each bilayer there is an independent emergent gauge field, which constrains the possible coupling between layers. The lowest order gauge invariant coupling of the matter field between layers is where i, j are layer indices. This coupling is irrelevant under RG at the 3d XY fixed point. However, we also need to include the monopole tunneling between layers.
This monopole tunneling term is gauge invariant. More importantly, this term is highly relevant in the 2d critical theory. This term will lock the gauge field between different layers. In the end, there is only one gauge field, which has 2d characters, across the 3d system. This suggest that the layers should first establish coherence and then go through the VPS-VBS transition all together. Physically, this means that the plaquette patterns of different layers along z directions synchronize. The transition are driven by proliferation of straight vortex lines along z-directions. This suggests that the VPS-VBS transition looks like a 2d transition despite the system is 3-dimensional. We have a transition that features a reduced effective dimension.

V. SUMMARY
In this work, we have shown how fracton physics gives important insight into melting transitions of valence plaquette solids. The topological defects of this type of spatial order are characterized by fractonic mobility constraints, as elucidated by a mapping onto a symmetric tensor gauge theory. Specifically, the individual vortices are completely immobile, while dipoles of vortices exhibit 1-dimensional behavior. This restricted mobility can prevent a direct condensation of single vortices, precluding a continuous transition from the VPS phase to a simple Néel state. Rather, a continuous melting transition of a VPS tends to involve condensation of vortex dipoles, giving rise to an intermediate phase between the VPS and Néel state. This intermediate phase can take different forms depending on the microscopic details of the dipoles, such as various types of bond order. A particularly interesting possibility in two dimensions is an algebraic bond liquid, which serves as a stable intermediate gapless phase, in agreement with numerics on certain 2d Heisenberg models [20]. We have discussed several signatures of this algebraic bond liquid, such as its structure factor, specific heat, and entanglement properties. We have also discussed the extension of these ideas to three-dimensional valence plaquette and valence cube solid phases. We show that, a particular anisotropic type of plaquette order can undergo a continuous transition to a bond-ordered phase via a quantum critical point that features a reduced effective dimension.

ACKNOWLEDGMENTS
We thank Cenke Xu, Meng Cheng, Chong Wang and T. Senthil for helpful discussion. YY is supported by PCTS fellowship at Princeton University. ZB is supported by the Pappalardo fellowship at MIT. This work was performed in part at Aspen Center for Physics, which is supported by National Science Foundation grant PHY-1607611.

APPENDIX: PARTON VIEW OF THE SPINON BOND LIQUID
In Ref. [73,74], the authors present a parton construction for the bond liquid phase, which captures several salient features of this phase including subdimensional motion and long-range entanglement. Here we apply their parton perspective to illustrate the emergence of subdimensional particle.
Writing the boson operator as b † a = v † a h † a . The Schwinger boson operator is further fractionalized into the vertical and horizontal bosons v † a , h † a . At the meanfield level, v † i v i+ex = h † i h i+ey = 0, as each vertical or horizontal parton only hops along the x or y direction, as dictated by the following Hamiltonian, H v,h = |(∂ x + a x )v| 2 + |(∂ y − a y )h| 2 + cos(∇ × a) + e 2 g (54) Each parton v † , h † is analogous to a 1d relativistic boson coupled with an emergent gauge field a. The strong fluctuation of the gauge field induces strong interaction between the two partons and projects the parton state to the physical Hilbert space. The vertical parton carries gauge charge a, so its algebraic correlation along x should be modified by the gauge fluctuation. Besides, a pair of vertical parton is charge neutral with power-law correlation v † i v i+ey v i+x v † i+x+ey = 1/|x| η . At the mean field level, this parton pair correlation exactly corresponds to the bond correlation, This parton construction provides a pictorial understanding of the bond liquid phase, which has low-energy behavior, including entanglement entropy, very similar to that of 1d relativistic bosons. In the mean field level, each parton forms a 1d spin chain along horizontal/vertical with elementary spinon excitation. After parton projection, the vertical/horizontal spinon motions are bound together as a spinon pair between a link hopping along the transverse direction.