Quantum phase transition of fracton topological orders

Fracton topological order (FTO) is a new classification of correlated phases in three spatial dimensions with topological ground state degeneracy (GSD) scaling up with system size, and fractional excitations which are immobile or have restricted mobility. With the topological origin of GSD, FTO is immune to local perturbations, whereas a strong enough global external perturbation is expected to break the order. The critical point of the topological transition is however very challenging to identify. In this work, we propose to characterize quantum phase transition of the type-I FTOs induced by external terms and develop a theory to study analytically the critical point of the transition. In particular, for the external perturbation term creating lineon-type excitations, we predict a generic formula for the critical point of the quantum phase transition, characterized by the breaking-down of GSD. This theory applies to a board class of FTOs, including X-cube model, and for more generic FTO models under perturbations creating two-dimensional (2D) or 3D excitations, we predict the upper and lower limits of the critical point. Our work makes a step in characterizing analytically the quantum phase transition of generic fracton orders.


I. INTRODUCTION
The discovery of topological quantum phases revolutionized the characterization and classification of fundamental states of quantum matter beyond the Landau symmetry-breaking paradigm. The extensive studies have brought about two broad categories of topological matter, the symmetry-protected topological phases [1][2][3][4] and topologically ordered phases [5][6][7], which exhibit short-and long-range quantum entanglement, respectively. Unlike the topological phases which necessitate symmetry-protection, the topological orders are characterized by finite ground state degeneracy (GSD) which is protected by the bulk gap and robust against any local perturbations, and host freely-moving fractional excitations [8][9][10][11]. The topological protection of the GSD can have potential application to fault tolerant quantum computation [12][13][14], particularly for the systems which host non-Abelian excitations [15][16][17][18][19].
Recently, a new basic class of topological phases, called fracton topological order (FTO) [20][21][22][23][24][25][26][27][28] were proposed and deeply broaden the understanding of the correlated topological states. In contrast to the conventional topological-ordered phases, the GSD of FTO scales up with the system size, and the fractional excitations of a FTO are immobile or have restricted mobility [29][30][31][32]. In particular, two types of FTOs are proposed and dis- * Corresponding author: xiongjunliu@pku.edu.cn tinguished by their excitations [24]. For type-I fracton orders, there are three categories of excitations, i.e. the 1D lineon excitations, 2D and 3D fractons, which are generated by 1D line, 2D membrane, and 3D body operators, respectively. The fractons are located at the ends of the line operator or corners of the membrane and body operators, and can move only along the line or in the way that the number of corners of the membrane/body operator remains the same. On the other hand, for the type-II fracton phases the excitations are located at the corners of a fractal operator and are completely immobile [20]. The restricted mobility of excitations is important for the stability of FTOs at finite temperature [33,34] compared with conventional topological orders, and may have advantageous in applying to topological quantum computation.
The exploration of the FTOs is still in the relatively early stage. Efforts have been made in relating the fracton orders to other novel physics [35][36][37], improving the basic concepts and characterizations of FTOs [38][39][40], and building new models for realization [28,38]. So far the studies are performed on the fine-tuned toy models, which are exactly solvable but hard to achieve in the real physical systems. On the other hand, with the protection by the topological bulk gap, the FTOs are expected to be stable against small external perturbations and may undergo quantum phase transition when the system is tuned far away from the fine-tuning point. To uncover the critical point of such transition is crucial to determine the phase diagram, while very challenging in general, and may provide insights into the possible realization of FTOs in the real physical systems. This therefore brings a challenging but highly nontrivial issue to the frontier of this research direction, namely, characterization of the quan-tum phase transition of FTOs, as we study in this work.
In this article, we aim to characterize quantum phase transition in type-I FTOs induced by external perturbation terms, with a generic formula to obtain the critical point for the phase transition. The GSD is unchanged when the strength of external term is less than the critical value, and breaks down when the phase transition occurs as the external term exceeds the critical value. The universal results of the critical points are obtained. The article is organized as follow. In section II we briefly review the basics of FTOs. Then in section III we define the quantum phase transition in FTOs and present the generic formalism. The analytic theory of the critical point of phase transition is presented in section IV, with a universal result being obtained by mapping the present framework to a random walk theory and applicable to a class of FTOs, including the X-cube model. In section V, we generalize the present study based on random walk theory to more generic FTOs, and show the upper limits of the critical points of the phase transition. A reasonable conjecture with numerical verification of the lower limit of the critical points is also presented. Finally, the conclusion is given in the last section.

II. REVIEW OF FTO MODELS
We start with a brief review of the basics of FTOs by introducing two typical models, the X-cube model [24] and Haláz-Hsieh-Balents (HHB) model [25]. After that we shall introduce the outstanding issue on the quantum phase transition of FTOs, as studied in this work.

A. X-Cube model
The X-Cube model is a spin-1/2 model defined on the edges of a cubic lattice [24], with the Hamiltonian where c runs through all cube centered at c and A c is a product of 12 σ x -operators defined on the edges of the cube, v runs through all vertexes and B The forms of the ground states shows clearly that the topological GSD originates from different choices of ∀v, i, B (i) v = 1 . The ground states thus generated can only be connected through flipping spins along straight lines that form noncontractible loops. In this model, the degeneracy is 2 6N −3 for the system with N 3 cubes.
There are two types of excitations [24]. (i) Lineons -generated by applying σ x on a ground state along an edge. The lineons are formed at the endpoints of the edge with energy w/2. If σ x is applied on consecutive edges that forms a straight line, two lineons will be formed on the head and on the tail. If the consecutive edges contain turning points, right at the turning points new excitations will be generated and cost extra energy. (ii) Fractons -generated by applying σ z on a ground state on the edges that forms a membrane perpendicular to the direction of the edges. The fractons then form on the corners of the membrane, and the energy of each group of four fractons is w.

B. Haláz-Hsieh-Balents (HHB) model
The HHB model is a spin-1/2 model defined on a bodycentered cubic (bcc) lattice [25], where each lattice site contains four spins, and the Hamiltonian Here r runs over different lattice sites, j = A, B, C or D for each flavor of spin and r A , r B , r C and r D are the shortest vectors connecting two sites in the direction {111}, {111}, {111} and {111} respectively. If λ → ∞, all spins in the same site are fully polarized. For λ 1, the first term of the Hamiltonian can be treated as a perturbation, and the effective Hamiltonian at the leading order can be written as whereÃ r = s=±jσ y r+srj s=±,j={B,C,D}σ and w ∼ J/λ 31 . The tilded Pauli matrices at r denote the corresponding Pauli matrices defined on the subspace of two ground states in each site r at λ = ∞. AllÃ r commutes and they have eigenvalues ±1 so that in the similarity of the deviation in the X-Cube model, the ground states can be written in the following form where e and o are distinct subsets of the lattice sites such that each forms a simple cubic lattice that are translated in the {111} direction. Similarly, the topological GSD originates from different choices of ∀o,Ã o = 1 , and is 2 12N −10 for 2N 3 lattice sites. The system also hosts two different types of excitations. (i) Applyingσ x on odd sites andσ y on even sites (or vice versa) along a straight line in the direction ±r j (j = A, B, C, D) yields a pair of lineon excitations. Unlike X-Cube model, the straight lines can be packed together to form 2D membrane or 3D body operators.
(ii) Applyingσ x on the sites along a line in ±(r A − r j =A ) directions also yields a pair of lineons. This operator can again be packed together to form 2D membrane or 3D body operators. Thus these types of excitations are classified as 3D fractons in our terminology.

C. Overview and FTO breaking down
All FTO models have topological GSD. The ground states can only be connected to each other through global operation, for example the spin-flipping along noncontractible loops. Therefore, any local perturbations cannot break the GSD. However, a strong enough global perturbation can break the GSD. For example, we consider a transverse magnetic term on X-cube model with where e runs through all edges of the cubic lattice. This external term can generate lineons on the endpoints of the edges. It is easy to see that for v → ∞, there is only one ground state with all spins being aligned along −x direction, so a quantum phase transition must occur at some finite critical v = v c , beyond which the FTO breaks down. This brings an outstanding question -How does the external term breaks the GSD? Answering the question is important to understand the topological phase region of FTO, and is crucial for the realization of FTO in real physical systems without fine tuning. Currently the FTOs are given in toy models which can be solved analytically. With an addition of external term, the system typically becomes not exactly solvable, and how to characterize the phase transition is a highly nontrivial but challenging issue. To approach this question, we consider the Hamiltonian in the following form where H 0 is unperturbed Hamiltonian of FTO with system size N 3 , the external perturbation V (v) is considered to the system with strength parameter v, andv is a dimensionless factor that will be taken as unity latter in the study. The perturbation may induce couplings in the ground states and thus break the GSD. Let ∆E(N, v) be the energy splitting of the ground state subspace. We shall see that, This implies v = v c is a critical value, below which GSD is unchanged, while beyond which the GSD shall break down. In particular, we shall show that when v < v c , for any two degenerate ground states |g and |g , we have lim N →∞ g |H eff | g = 0, and for v > v c , for any |g , there exists other ground states |g and Here H eff is an effective Hamiltonian defined on the subspace of the ground states that captures the transitions among them. Note that to avoid finite-size effect, the system size N must be taken as infinity.

B. A general framework
We consider that the external extra term is written in the following generic form of fracton operators whereÔ r generates or annihilates two lineons ( Fig. 1(a)), four 2D fractons ( Fig. 1(b)) or eight 3D fractons ( Fig.  1(c)) at its vicinity. In Fig. 1, the operatorÔ r is denoted by red line, square or cube, respectively, and the excitations created by it are represented by the red spheres.
To facilitate the study on how the external term V (v) changes GSD, we split the Hilbert space into two subspaces where ψ represents the full Hilbert space, ψ g denotes the ground state subspace, and ψḡ denotes the subspace of the states with excitations. With this notation, we project the whole Hamiltonian onto the ground state subspace, and obtain the effective Hamiltonian (Appendix) P i−1 is the set of integer partitions of i − 1, i.e. s 1 + s 2 + · · · + s m = i − 1. The subscripts gg (orḡḡ) means that the corresponding operator is projected onto the ground state subspace (or the excited state subspace). Similarly, the subscripts gḡ andḡg mean that the corresponding operators connect the ground state with excited state subspace. Finally, Q (s) 's are the factors of inverse energy scale, which denotes the energy difference between ground state and intermediate excited states of the FTO (see more details in the Appendix). Since only the operation on noncontractible loops can induce transition between two different ground states, the leading order correction to the ground state subspace must be the N -th order. In the leading order contribution, only the terms with Q (1) = (E 0 −H 0,ḡḡ ) −1 have non-zero matrix element between two different ground states. Therefore, the effective Hamiltonian by taking into account the leading order contribution reads Any lower order perturbation cannot couple different ground states, since a noncontractible loop must be formed by at least N times ofÔ r operations. The energy splitting in the ground states depends on two aspects. The first is the coupling between each two ground states (e.g. |g and |g ) up to all higher orders. In particular, the matrix element for the (N + r)th-order contribution is given by H (N +r) eff,gg = g|H (N +r) eff |g , with r = 0 being the leading-order contribution. The second is the number of ground states |g which are coupled to a certain fixed ground state |g by the external term, denoted by f (N +r) . The energy splitting from the asymptotic behavior should be given by Intuitively, the ground state subspace can be treated as a flat band of strongly correlated many-body states before phase transition. For the strength of external perturbation exceeding a critical value, the flat band is deformed to be a dispersive one, with ∆E being the bandwidth.

IV. CRITICAL POINT OF THE QUANTUM PHASE TRANSITION
We proceed to study the key issue, the critical point of the quantum phase transition in the FTOs. We shall show that, for the external perturbation creating only lineon type excitations, the critical point of the phase transition can be given by a universal formula and in some special cases the result given by the formula is exact, as presented in the following theorem.
Theorem 1: Let w be the energy of a pair of separated lineons and w i be the energy of two lineons along different directions crossing at an intersecting site. Then the quantum phase transition occur exactly at This theorem can be proved by solving the leading-order and all higher-order contributions in the thermodynamical limit. We first begin our discussion on the general properties of the contributions of different order in the subsections A and B. In the latter part of subsection B and subsection C, we assume that f = w and prove the theorem in this regime. In subsection D, we argue that the exact result can extend to the case f > w and employ a similar method to give a universal formula for f ≤ w under certain approximation.

A. The leading order contribution
We start with the leading order contribution to the matrix element of H eff between two ground states |g and |g . Let the two ground states be connected by a noncontractible straight line ofÔ r operators, i.e.
where r runs through all sites along a straight line. Denote by V γγ = γ |V | γ the matrix element of V (v) and F γ the number of pairs of generated lineons at the intermediate states |γ . The matrix element of the effective Hamiltonian reads The denominator implies that the energy of the intermediate state depends purely on the number of lineons in the system. The summation is over all permutation of {β 1 , · · · , β N −1 }, the last line in above equation is obtained by noticing that V γγ is either 0 or v, and the prime means that the summation is over all the intermediate states with V γγ = v. As an example, when N = 6, one of the terms in the summation is shown in Fig. 2(a,b), where the number of lineons F β1∼5 of the intermediates states |β 1∼5 are 1, 2, 2, 1 and 1 respectively. Exactly counting the summation in Eq. (15) is highly nontrivial when N is large. With the details presented in the Supplementary Material [41], we find that the exact value of the leading order perturbation takes a novel form The leading order contribution is deeply related to the so-called Catalan number [41], C N −1 2N −2 /N , as illustrated in Fig. 2(c). The N th Catalan number is defined as the number of paths along the edges of a grid with N × N square cells from the bottom-left corner (purple dot) to the top-right corner (green dot) with two restrictions. (i) Each step of the path has to be either rightward or upward. (ii) The whole path has to be below the diagonal line (painted red). Or equivalently, the path has to pass through N odd lattices below the diagonal, including the lattice rightward to the bottom-left corner (b i ), the lattice downward to the top-right corner (b f ) and N − 2 more odd lattices (circles on the vertexes of the grid). In particular, the N − 1 th Catalan's number C N −1 is related to the leading order perturbation by H To better understand the leading order contribution, we show below that it can be mapped to a random walk problem, which has the important merit to facilitate the further study of higher-order contributions and the more generic extra terms creating 2D/3D fractons. More specifically, we consider a random walk starting at L 1,1 in a network L i,j , where 1 ≤ i ≤ min(j, N −j) equals the number of pairs of lineons in the j-th step of the intermediate process, with 1 ≤ j ≤ N − 1, and moving N − 2 steps with probabilities in each step [41] For each step, we set the reward gained by the random walk as 1/i so that the mean of the product of rewards matches the matrix element H (N ) eff,gg after summing over all paths, except for the prefactor v (−v/w) N −1 . The exact result is given by with i 1 = i N −1 = 1. A simple example with N = 6 of the random walk network is shown in Fig. 3. The shaded path shows a case where i 2 = i 3 = 2 and i 4 = 1 and it corresponds to the sum of all possible intermediate states so that F β1∼5 are 1, 2, 2, 1 and 1 respectively. One example of such states are shown in Fig. 2(a). The sum of the rewards of the random walk equals the N times of the Catalan number, recovering the result in Eq. (16).

B. Higher order contributions
The mapping to random walk problem can be generalized to the study of all higher-order contributions H (N +r) eff,gg , with r being even integers. Note that the odd r is forbidden since the effective Hamiltonian with an addition of odd number ofÔ r operations must scatter ground states to excited ones. It is convenient to consider the higher order terms in the two different case, namely, the higher-order contribution with finite r and the contribution with infinite r, respectively.
To get the intuition of the higher-order contributions, we consider first the lowest higher order term with r = 2. This type of higher order terms contain N + 2 number ofÔ r operators and some of them has to be repeated. SinceÔ 2 r = 1, the repetition is actually equivalent to a backward step in the random walk process. For r = 2, the network L is allowed to move back by one step at an intermediate step. As a simple example with N = 6 and r = 2 shown in Fig. 4(a), theÔ 4 operation is repeated at the third step. When represented by the random walk, this process is equivalent to moving back at the second step. Specifically, the network L is expanded to a new one L (1) , in which the random walk includes a new boxed subnetwork to account for the process of moving backward by two steps, as shown in Fig. 4(b). For the generic configuration, let P i be the probability that the walker moves back at i-th step. The average reward, and so the matrix element, can be approximated as Here R N,i is the average reward of the random walk process which includes the backward moving at the i-th step. It can be shown that whenever P i is finite, the asymptotic behavior of R N,i remains the same as H (N ) eff,gg for large N (see Supplementary Material [41]). This further determines the asymptotic behavior of H (N +2) eff,gg . For generic finite r, the above argument iterates r/2 times to form a new network L (r/2) and the similar result can be obtained. It then follows that that the matrix element H (N +r) eff,gg with finite r should be dominated by (4v/w) N +r and the sum of the perturbation series up to such higher orders is exponentially small if v < w/4.
The remaining question is how the perturbation grows as r tends to infinity, noting that the iteration cannot be done by infinite number of times. For this purpose, we first assume that the energy of 2 lineons coming from mutually different directions meeting at the same point is 2×w/2 = w, i.e. lineons coming from different directions are independent. Let r = ζN , where ζ is a real number. After some algebra we find for higher order perturbations with infinite r that [41] lim N →∞ whereλ = N/4 is the reward average over all possible configurations, i.e. all possible intermediate states connecting |g with |g . The above is an important novel result, and can be understood in the following way. For the higher order term with infinite r, the network L (r/2) allows to move backward infinite number of times so that the process are close to unrestrained random walk, namely, the random walk can perform in arbitrary directions. Thus the growth of the average of product of rewards is captured by the reward averaging over all possible configurations. With this one can show that the matrix element grows by 4v/w as r increases one order. [41].

C. Energy splitting
Summing over all-order contributions yields the total transition matrix element As the growth of H (N +r) eff,gg is 4v/w as r → ∞, the growth of prefactor C (r) is subexponential. On the other hand, the number f (N ) obeys power-law of the system size and is of the order N 2 . This is because for a reference ground state |g , the number of other ground states that can be coupled to it through NÔ r operators is ∼ N 2 (see Fig. 5). Further, the number f (N +r) does not grow until r is comparable to N , since an addition of finite number ofÔ r operations couples no new ground states. Thus the number f (N +r) also grows subexponentially versus N and r. With these results we conclude from Eq. (13) that the system, which we assume the lineons coming from different directions are independent, has a critical point given by For the perturbative term weaker than the above value the GSD of the FTO is unaffected. We provide an intuitive understanding of the above critical point by partially relating the perturbative study to the 1D transverse Ising model. Consider the case for two ground states |g and |g which can be connected through applyingÔ r operators on a single straight line, i.e. |g = r∈LÔ r |g . Here L is a straight line forming a non-contractible loop. We further consider the following HamiltonianH =H 0 +H 1 : whereξ i,i+1 is an operator with eigenvalue 1 when there is a lineon between site i and i + 1, and 0 when there is not.Ô i remove (create) when there is (not) a lineon in the vicinity of the site i. It then follows that With the above results we can verify that the matrix element ofH between |g and |g is identical to H eff,gg . Consider now the mapping thatξ i,i+1 → (1 + σ z i σ z i+1 )/2 andÔ i → σ x i , which keeps the above commutation relation (23). This leads toH → H Ising = w i , which is the 1D transverse Ising model with the critical point given by v c = w/4. From this result we can see that at the critical point the bulk gap of the FTO transition should close. It is noteworthy that this mapping is valid for part of the configurations, not exactly identical to the total 3D effective Hamiltonian (11), but gives the correct critical point if the lineons coming from different directions are independent. Also the mapping is valid only for the regime before the phase transition.

D. Correction due to intersecting lineons
More generally, the result of critical value for system with lineon excitation depends on the energy w i of two lineons along different directions crossing at intersecting site. (i) w i = w. The global line operators can be treated as independent so the critical point is Then the matrix elements between two arbitrary ground states may be finite only at some value v > w/4. However, the critical point remains at v = v c1 = w/4 since the leading order contribution between |g and |g that are connected by a single global line operator is finite  (20) so that H eff,gg = r C (r) 4vf (r) /w N +r . An estimation of the correction factor can be given below (details are found in the Supplementary Material [41]. We start the estimation by noting that the correction factor can be written as where i denotes the i th intermediate step, · · · β denotes the weighted average over all configurations β and f β,i captures the ratio of energy changes due to intersecting lineons. For example, if there are 4 non-intersecting lineons and 2 intersecting in a given configuration and the energy of each intersecting lineon is w/3, then f βi = (4 × w/2 + 2 × w/3)/(6 × w/2) = 9/8. We start the estimation of f (r) by assuming f β,i to depend on F βi only. This assumption is reasonable since the probability of lineons to be intersecting is proportional to the density of lineons. With this assumption, the correction factor becomes where M βi is the number of configurations for a given F βi .
The correction factor f where w i ≤ w is the energy of two intersecting lineons and the critical point of this more generic case can be approximated by This result is universal for the type of extra terms that create only lineon excitations. In particular, the quantum phase transition of the transverse field X-cube model, as given in Eq. (6), can directly obtained from the above result, and v X-Cube ≈ (0.87 ± 0.04) × w/4. The origins of the error are discussed in the supplementary material [41]. The value agrees with the numerical result v pCUT X-Cube ≈ 0.9196 × w/4 by the method of perturbative continuous unitary transformations [42].

V. CRITICAL POINTS FOR THE GENERIC EXTRA TERMS
The random walk theory developed based on the effective Hamiltonian approach can be further applied to study the more generic cases beyond lineon type perturbations. For those cases, the external perturbation can induce also the 2D and/or 3D fractons, e.g. for transverse field coupling toσ x,y in the HHB model, namely, the termÔ r generates or annihilates four 2D fractons (as shown in Fig. 1(b)) or eight 3D fractons (as shown in Fig. 1(c)) nearby the site r. We expect that when the strength v of the extra term V (v) exceeds certain critical value, the phase transition occurs and the FTO breaks down. Similar to the study for lineon type perturbations, the exact critical points can be obtained in three key steps. (i) Identify how the leading order contribution grows as N → ∞. This encodes also how the higher-order contributions grow with finite r. (ii) With fixed N , identify how the infinite-r higher-order contributions grow. (iii) Summing over all order contributions to obtain the energy splitting. Similarly, for a fixed ground state |g , the number f (N +r) of |g coupled to it through leading-order perturbation must be power-law function in N . This value does not grow until r is comparable to or over N (for 2D and 3D fractons). Thus for N → ∞ the growth of f (N +r) is again subexponential.
Upper limits.-The exact result of the critical point is not analytically available in the current stage. The challenge arises from working out the leading-order contribution, as studied later. However, the higher-order contributions with infinite r can be exactly obtained by the random walk theory. For the case with 2D or 3D fractons, the sole consideration of (ii) gives an upper limit v In particular, we show for the 2D fracton case, the growth of H (N +r) eff,gg is 16v/3w as r = ζN increases by one order when |g and |g is connected by a 2D membrane of operatorsÔ r . This can be seen by noticing that at any point r 0 (on which a 2D fracton is defined) that are surrounded by four operatorsÔ r1 ,Ô r2 ,Ô r3 andÔ r4 (listed in clockwise direction, as shown in Fig. 6(a)), the number of fractons F 2,r0 at r 0 can be written as where s 5 = s 1 and s i = (1 +Ô ri )/2. As shown through random walk theory in section IV B, the growth of high order contributions with infinite r is proportional to the average numberF 2,r0 of excitations over all possible configuration. Note that each term of s i contributes a factor of 1/2. It then follows that thatF 2,r0 = 3/4. Therefore, the growth required is 4 × 4v/3w = 16v/3w, giving the upper limit of the critical value v c2 ≤ v (u) c2 = 3w/16. The upper limit of the critical point can be obtained for the case with 3D fractons in a similar way. At any point r 0 that is surrounded by eightÔ ri operators (as shown in Fig. 6(b)), the number of fractons F 3,r0 at r 0 is given by Here a i1···i8 are coefficients whose values can be found in [41]. The average number of fractons is then obtained thatF 3,r0 = 29/32, and the growth of H (N +ζN ) eff,gg for infinite r is 256v/29w as r increases by one. Thus the upper critical point reads v c2,c3 drives the FTO into a trivial phase. Lower limits.-There remains an interesting issue on the lower limits of the critical point of the phase transition when the extra term can create higher dimensional (2D or 3D) fractons. The lower limit lies in identifying the growth of the leading-order contribution, whose exact solution is currently not available. However, we provide below a reasonable conjecture of the lower limits.
Conjecture: The critical points satisfy v c2 ≥ w/16 and v c3 satisfies v c3 ≥ w/64 when the extra term can create 2D and 3D fractons, respectively. The lower limits can be approached by estimating the maximum value of the transition matrix element between the ground states. The conjecture is based on the observation that for the 2D fractons, the number of configurations connecting two ground states in the same order perturbation is more than the case for 1D fractons due to the higher dimension in the operations. This is already seen in the exact study of the upper limits. The maximum number of the configurations could be the square of that for the case with 1D fractons, namely the square of the Calatan number as shown in the case with lineons. In this case, the lower limit of the critical point is given by v In the similar way, we conjecture that the lower limit for the case of 3D fractons is given by v (l) c2 = w/64. The conjecture is clearly supported by numerical results, as given in the Supplementary Material [41]. To analytically confirm the lower limits and further the exact result of the phase transition critical point in the generic case deserves future great efforts.

VI. CONCLUSION
We have presented an analytic theory of the quantum phase transition in FTO induced by external terms which create fracton excitations. We developed an effective Hamiltonian approach which is model independent to study the critical point, beyond which the ground state degeneracy (GSD) breaks down. This approach can be mapped to the random walk problem, allowing to study the critical point of quantum phase transition induced by fractons of different dimensions. For the lineonic external term, we provide a universal result of the critical point, × 0.26 for w i < w, with w being the energy of a pair of lineons and w i being the energy of two lineons meeting at the same point. This result can be directly applied to X-cube model with transverse Zeeman field [24]. For the external terms creating 2D or 3D fractons, e.g. the transverse field in the HHB model [25], we show an upper limit and conjectured a reasonable lower limit with numerical verification for the critical phase transition point. Beyond the upper limit (Below the lower limit) the FTO breaks down (survives) in general. The analytic theory presented here shows insights into the understanding of the quantum phase transition in FTOs, and can be applied to predicting the full phase diagram of generic FTOs, which may facilitate the physical realization of such exotic orders. In this Appendix, we present the rigorous calculation of arbitrary n th order degenerate perturbation on the manybody ground states, whose degeneracy is assumed to be unchanged up to the n − 1 th order. Let H = H 0 +vV , wherev 1, H 0 denotes the unperturbed Hamiltonian, and V is the external perturbation of the system. Write down that H 0 = H 0,gg +H 0,ḡḡ and the external term V = V gg +Vḡḡ+V gḡ +Vḡ g . The subscripts gg (orḡḡ) means that the corresponding operator is projected onto the ground state subspace (or excited state subspace). Similarly, the subscripts gḡ andḡg mean that the operators connect the ground and excited state subspaces.
Denote by P the projection operator onto the ground state subspace whose energy is E. For the total Hamil-tonian H, the secular equation is Then the effective Hamiltonian projected on the ground state subspace is given by: With the power series E = j E jv j+1 we obtain that where E 0 is the energy for the unperturbed ground state, P i is the set of integer partitions of the integer i, and the partition with different orders are counted different. For example, 1 + 2 = 2 + 1 = 3 are two different partitions of i = 3. In the summation s j is the j th term in the integer partitions s and |s| is the number of integers partitioned. For example, if s = 3 + 1 + 5 + 4 ∈ P 13 , one has that s 3 = 5, s 4 = 4, and |s| = 4). For convenience, we further , with Q (i) in the expansion given by For example, the 4 th order term Q (5) With these results the effective Hamiltonian projected onto the ground state subspace can be further simplified in the following way  Rules 1 to 5 are conventions to denote a specific perturbation term. The 6 th rule is to guarantee that the transition must be from ground state to ground state, and that all intermediate states are excited states. If a horizontal line can be drawn, at that point the state transforms back to the ground state, which is prohibited (c.f. Appendix). The 7 th rule is to forbid a simultaneous creation and annihilation of a lineon. With this seven rules, it is possible to consider the terms in the formv n V gḡ Q −1 0 Vḡḡ n−2 Q −1 0 Vḡ g . In particular, Fig. S1 shows an example. (a) shows a unit cell of a simple cubic lattice with size N 3 with periodic boundary condition. Suppose N = 4, then (b) may express terms likev 4Ôŷ Q −1 0Ô 0 Q −1 0Ô 2ŷ Q −1 0Ô 3ŷ . Note that the leading order perturbations must be drawn with a single loop and describe a non-contractible loop in FTO. For diagrams with more than one loop, for example Figure (c), it may express terms likev 4Ô 2ŷ Q −1 0Ôx Q −1 0Ô 2ŷ Q −1 0Ôx , so that the lineons at the neighborhood of 2ŷ andx just annihilate with their partners nearby without mutual communication. In addition, it is worth noting that in the diagrammatic language if the smallest non-contractible loop has N sites, there must be at least a loop of order at least N (i.e. with at least N vertexes) to connect different ground states. Such loops are called main loops. Therefore, among two diagrams, diagram (b) represents off-diagonal terms in the effective Hamiltonian. More generally, for the terms of the formv n V gḡ Π r−1 i=1 Q (αi) Vḡḡ Q (αr) Vḡ g , where some of the α i are larger than 1 and (α 1 · · · α r ) ∈ P n−1 , the last two rules (viii) and (ix) are needed. The transitional lines and subdiagrams are to record the places at which Q appear and the corresponding E s terms as in the Appendix. For example, a sixth order term is shown pictorially in Fig. S2. The diagram represents a termv 6  decomposition principle. This physical properties can be transformed into diagrammatic language as follows: given a diagram D with n > 1 loops, draw a family of diagrams D in the following step. (i) Add arbitrary transitional lines and assign that number on the external vertexes of at least one loop. (ii) If any subdiagram contains more than one loop, stretch the loops arbitrarily so that the internal vertices can be of different order by having different relative height. (iii) All diagrammatic rules are still enforced. For example, the 6 th rule -when stretching the loops of the subdiagrams with more than one loop, there should not be a horizontal line that can pass through a subdiagram without intersecting any solid lines. Then, we claim that Claim: Given a diagram with n > 1 loops, the sum of the contribution of a family of diagrams containing such diagram is zero. Label the proposition for the case with n loops and each loop have r i + 1 internal vertexes to be P n (r 1 , · · · , r n ). Take the proposition of P 2 (1, 1) as an example (Fig. S3). Assume at the moment that the left loop produce excitations with energy of E a and the right loop E b . Label the sum of the family of diagrams as S 2 (1, 1). Then (S1) Since a diagram with more than 2 loops can be decomposed into 2 groups, it is actually suffix to prove P 2 (n, m). So the proof can be done by mathematical induction. Note here that P 2 (1, 1) is explicitly shown. Suppose now P 2 (n ≥ 1, 1) is correct. Denote a i to be the energy under the i th internal vertex for the order n loop and b is the energy of the order 2 loop. The sum of the family of diagrams S 2 (n, 1) is S 2 (n, 1) = 1 b 2 1 (a 1 + b) · · · (a n + b) + n i=1 1 b 1 (a 1 + b) · · · (a i + b)a i · · · a n + n i=1 1 b 1 a 1 · · · a i (a i + b) · · · (a n + b) + n i2=i1 n i1 1 a 1 · · · a i1 (a i1 + b) · · · (a i2 + b)a i2 · · · a n − 1 a 1 · · · a n n i 1 Denote the i th term of S 2 (n, 1) be q i,n , then S 2 (n + 1, 1) = q 1,n 1 a n+1 + b + q 2,n 1 a n+1 Therefore, by mathematical induction, the proposition P 2 (n, 1) is true for all n. A similar calculation gives induction on m as well so that P 2 (n, m) is true for all n and m. The above sections tell how an arbitrary order of valid perturbation looks like. In particular, the only leading order term that survives look like where there are N 's V in the equation, since all other terms that are of the same order have fewer V so that they together cannot form a non-contractible loop. On the other hand, the appearance of Q (s>2) in the perturbation series makes the diagram contain subdiagrams so that it is canceled by all possible combinations of such kind of diagram (as shown in the last section). Therefore, the higher order term can also be written in a similar form where there are N + r's V in the equation.

A. Calculation of the leading order perturbation
To break the ground state degeneracy, the effective Hamiltonian must have non-zero off-diagonal elements. In particular, given the external term as V = v rÔ r , the matrix element of leading order is g H (N ) eff g , where |g and |g differs only by a global line operator. Denote α i to be the place whereÔ i is defined so that iÔ i connects two ground states. Write V γγ = γ |V | γ , then g H where V βiβj is the matrix element β i |V | β j and is equal to either equal to v or 0, F β is the number of pairs of lineons at the state β and the prime in the summation sign means the summation over intermediate states are limited to those with V γγ = v. Let v ps be the result of the primed summation in (S1) . Then Claim: If N is the system size, then v ps (N ) = If the claim is true, then as N → ∞, the asymptotic behavior of the matrix element is This value is zero whenever 4v/w < 1 and thus v c = w/4 to this order. Suppose there are N lattices α 1 , α 2 , · · · , α N so that the transition occurs when N i=1Ô i is applied on |g . For leading order perturbation, the sum is actually over permutations group of α 1 , α 2 , · · · , α N . The permutation corresponds to the order of howÔ i are applied on the ground state |g , as mentioned in the main text (or described in Fig. S4

Exact calculation of the leading order perturbation
where p = (p 1 , p 2 , · · · , p N −1 , p N ) ∈ S N and F pi 's are the number of pairs of lineons after applyingÔ on the lattices α p1 , α p2 , · · · , α pi . Thus for every permutation p, we can define the following angle notation p <> = a 1 , a 2 , · · · , a i , · · · , a N −1 . The notation with N − 1 numbers refers to a permutation p ∈ S N . The numbers a i 's mean that there are a i pairs of lineons after the operation i j=1Ô pj . For example, for the examples shown in the Fig.  S4(a,c,e), the angle notations of the permutations are 1, 1, 1 , 1, 2, 1 and 1, 2, 2, 1, 1 , respectively. All permutations with the same angle notations are called equivalent permutation. Below are the statements and proofs of two claims, which count the number of equivalent permutations in two different cases.
Counting Equivalent Permutations (I)-.Note that |a i+1 − a i | ≤ 1, since an application ofÔ on a lattice either creates (1) or destroys (−1) two lineons or moves one lineon to another place (0). Here, we shall first consider the case that the angle notation satisfies a i+1 − a i = 0. In this case, there are N a i equivalent permutation. To prove this statement, suppose the first element of any permutation is 1 so that at the end of calculation, the equivalent permutation should be multiplied by N . Now, there are three different cases: (i) Suppose p <> is of the form 1, 2, · · · , n − 1, n, n − 1, · · · , 1 . Then the first n terms in p ∈ S N =2n must be odd so that there are in total (n − 1)! ways to paint the first n edges and there are also n! ways to paint the last n edges. Therefore the total number of equivalent permutation is 2n × n!(n − 1)! = N i a i .
(ii) Suppose p <> is of the form 1, · · · , n, n − 1 · · · , m, m + 1, · · · , p, p − 1, · · · , 1 (i.e. a i 's increase from 1 to n, decrease from n to m, increase from m to p, and finally decrease from p to 1). Here, whenever a i − a i−1 = ±1, the p i must be odd (even), otherwise there will be some j that satisfy a j − a j−1 = 0. The counting will be done in four steps: (I) Similar to the first n terms in the last case, and it contributes (n − 1)!.
(iii) Suppose p <> is in the most general case such that a i −a i−1 = ±1. Suppose the numbers a i increase and decrease so that they peak at n 1 , · · · , n r+1 and trough at m 1 , · · · , m r successively. Here N = 2 r (n r − m r ). Following the same spirit of the last case, we split the counting into several steps, where some steps will be iterated. (I, II) are similar and they contribute (n 1 − 1)! × P n1−m1 n1 . Then the following two steps are iterated. (a) The next n 2 − m 1 action are to increase the number of segments by adding odd segments in the list. It is similar to the first step when we treat the old m 1 segments as m 1 odd numbers with a particular ordering. Then total number of such action is (n 2 − 1)!/(m 1 − 1)! (The division is present since the order is chosen in the first step). (Make the replacements n 2 → n j+1 , m 1 → m j and m 2 → m j+1 for j th iteration process.) (b) Similar to second step, there are P n2−m2 n2 to connect the n 2 odd segment to m 2 segments. (Make the replacements n 2 → n j+1 and m 2 → m j+1 for j th iteration process.) (III, IV) are again similar and they contribute the factors H nr+1−mr mr × n r+1 ! × N . So This completes the counting of equivalent permutations so that whenever a i+1 − a i = 0 for all i the equivalent permutations is N i a i .
Counting Equivalent Permutations (II)-.A generalization to the above case is that a i+1 − a i = 0 for some i. For this, we let q to be the number of such i. Similar to the step (a) of case (iii) in the last proof, we imagine there are a i segments, and they behave likes a i numbers with specific ordering. When a i+1 − a i = 0, the i + 1 th action is to lengthen an existing segment. Since each action of lengthening gives a factor of 2a i . (Each segment has two ways to lengthen -on the left and on the right.) Therefore, there are extra factors of 2a i /a i = 2 for every i such that a i+1 − a i = 0. So, Knowing the number of equivalent permutations in the most general case, we can count v ps by summing over all families of equivalent permutations. Then we arrive the following result.
Summing over Equivalent Permutations-.For every permutation p ∈ S N with angle notation p <> , define f (p) = is equivalent to (S2). Now we try to prove this formula in order to complete the calculation of v ps . Let S N to be the set of families of equivalent permutations in S N . Its element p ∈ S N means p ∈ p, so p 1 = p 2 implies p 1 and p 2 are equivalent. Then the summand simplifies to where a i 's are the numbers in the angle notation p <> and q( p) is just the number of i such that a i+1 − a i = 0 for any p in the family of equivalent permutation p. Thus q( p) can also be represented by any element in p, q(p). Define q r (p) to be the number of i ≤ r such that a i+1 − a i = 0. Then for all λ, where 1 ≤ λ ≤ N − 1, the following equation holds: where the summation sign with a i means to sum over all possible values of i th term in p <> and b ρ,λ = a1···a λ−1 2 q λ (p) a λ =ρ such that the λ th term of the angle notation p <> is ρ. The first equal sign is true because the summation over family of equivalent permutations is equivalent to the summation over all possible angle notations. The second equal sign is true because the exponent 2 q(p) can be factored into two parts, with each part depends only the corresponding a i 's. Now, the value of the summation is just b 1,N −1 since the last element of the angle notation must be 1. Then from the definition of b ρ,λ , the following recursion relation holds where the third equal sign is true because ρ = ρ ± 1 or ρ = ρ, and for the last case, the summand contributes an extra factor of 2. The boundary condition are given by b 0,λ = b λ>ρ,ρ = 0 and b 1,1 = 1. This problem is equivalent to a counting problem in mathematics: suppose there is a grid of N by N square. Count the number of way to walk from the bottom-left corner to the top-right corner on the edges of small square with the conditions: (i) Going leftward or downward at any step is prohibited; (ii) The first step must be a step towards right; (iii) Stepping over the diagonal is prohibited. For example N = 6, as shown in Fig. S5(a). The number of steps to get to the lattice b i,j from the bottom-left corner satisfy the above relation (S11). The counting result is well known (S12)

B. Higher order perturbation: random walk theory
We proceed to study the higher order contributions by mapping the present problem to a random walk process. We shall see that this mapping has key merit in studying not only the higher order contributions, but also the phase transition induced by extra terms creating 2D and 3D fractons. To map calculation to a random walk theory, we first turn it into a probability problem. In the most general case, consider a function f : A → R and suppose a summation S over A is wanted: The problem can be turned into a probability question given that the size of A can be calculated. Suppose A is a set of events that occur with uniform probability and the reward of a particular event of a ∈ A is f (a) × |A|, the expected reward S on the set of events A is If a non-zero normalized function g (i.e. a∈A g(a) = 1) is defined on A, a further generalization can be made: Each particular event a has g(a) probability to appear with reward f (a)/g(a), then In the same spirit, the calculation v ps can be turned into a probability question: where |p <> | / |S N | can be interpreted as the probability distribution and |S N | / i a i the reward. For higher order perturbation, a similar summation is needed, and a similar process can be made, so that the reward function is simple to find out, whereas for the probability distribution, an estimation can be made so that the size of the ultimate sum can be approximated. Note that the calculation of the equivalent permutation splits into different stages, and each stages are independent of the other stages (c.f. the proofs recorded in last section). Therefore, the calculation of probability |p <> | / |S N | can be split into products of independent probabilities, each describing how the system to jump from a i to a i+1 with different probability to change to number of lineons. Let L i,j be the places that the random walker passes through, where the subscript i represents the number of pairs of lineons at j th step. Then quantitatively, where P N | Li,j →L i+δ,j+1 is the probability of the random walk at j th step go from L i,j to L i+δ,j+1 (δ = ±1, 0).
FIG. S6: (a) An example of equivalent random walk for N = 6. Li,j represents the j th number in the angle notation is i. The walker must go rightward throughout the process. (b) An example of equivalent random walk for r = 2, N = 6. The walker is allowed to move backward once. The step that move from primed to unprimed lattice is the step when the walker move backward. The steps encircled by the black dashed rectangle are the inserted steps as described in the main passage. For both graphs, the probability for each step can be read from its color -Cyan: 3/5, Orange: 2/5, Purple: 1/6, Blue: 2/3, Red: 1/3 and Black: 1. (c) Generic form of network for r = 2. M denotes the master node, from which the walker has been assigned to go to Network i randomly, which is a network describing the walker must step backward after i th step. D is the dummy designation as described in the main text.
Consider N = 6 (see Fig. S6(a)) for an example. Let α i , where i = 1, 2, · · · 6, be the lattices on which the operator O defined. For any permutation with angle notation a 1 , a 2 , a 3 , a 4 , a 5 , reinterpret it as a random walk with four steps, where probability of path taken in the j th step represents how often the system (i) creates a pair of lineons (for a j+1 − a j = 1), (ii) move a lineon (for a j+1 − a j = 0), or (iii) destroy a pair (for a j+1 − a j = −1) of lineons in the j + 1 th action. Then the average reward is which is just C N −1 2N −2 for N = 6. This current reinterpretation helps generalize the calculation to the higher order perturbation, although only estimation can be made. For the higher order perturbation, by the cluster decomposition principle, only terms with one loop are needed to be considered. They are the terms in the form λ r Vḡ g (Q −1 0 Vḡḡ) N +r−2 Q −1 0 V gḡ (c.f. (A4)), for the N + r th order perturbation. Let S (r) N be the set of the permutations of N + r numbers taken from {1, 2, · · · , N }, where each number appears odd number of times. In particular, the set S N is a special case of this definition for r = 0, i.e. S where (p 1 , p 2 , · · · p N +r−1 ) ∈ S ps .

High order terms with finite r
If r is a finite number, the sum v (r) ps can be modeled by the same random walk problem, with the exception that the walker can walk backward. Take r = 2 and N = 6 as an example. The network L is allowed to move backward once. Suppose the walker move backward at step 2 (as in Fig. S6(b)). Then it can be reinterpreted as a new network, where everywhere are the same except that an extra network (dashed boxed) is inserted. More generally, pick an element p from S (2) N , the numbers p 1 , p 2 · · · p N +2 are all unique except that three numbers are same. Let them be p β = p i+1 = p γ , where β < i + 1 < γ. So for the first i th actions, the effect of i + 1 th and i + 2 thÔ r operators are to shift an appliedÔ r to other places. These two actions are captured by the boxed region as shown (b).
Suppose the connections and the probabilities in the boxed region has been already given exactly. The only thing left is to calculate the probability of the walker to walk backward after i th step. Now the number of permutations p ∈ S (2) N such that p i+1 is the same as exactly one of p 1,2,··· ,i is (i + 1)(N − i)(N − 1)! , therefore the probability P i is Thus, the full network can be represented as in Fig. S6(c). The step going from the master node to the i th network, constructed in the above method and has the probability P i . At the end of the network, the walker moves to the dummy designation (which do nothing to the reward) to finish the journey. In this way, the value v (2) ps can be exactly calculated. Although the probabilities in the boxed region can be calculated exactly, they are unimportant as N → ∞. This can be seen in the following argument. (i) In this limit, nearly all walkers walk to the designation through network number i, such that both i and N − i are comparable to N . (Other paths have negligible probabilities to walk on.) (ii) Within these network, nearly all random walkers pass through L j,i such that j is comparable to N , since the average number of segments at i is of the order N and the standard deviation √ N . Now, since the boxed region contains 2 steps, the reward multiplier j within the boxed region must lies between j ± 2, which is unimportant compare to N . Therefore, the average reward R N of the random walker R N is where P (q) is the probability density of the random walker to go backward at qN step, and λ(q) is the average number of pairs of lineons, in the original system. Since the integral is slowly varying in N , the ratio R N /R N −1 is still the same as v ps (N )/v ps (N − 1). Now for r to be any finite even number, the arguments are still the same, except the above step are iterated finite number of times, and so the growth H eff,gg remains the same after including higher finite order terms. The argument can be further quantified by finding the mean and standard deviation of the number of segments. They are of order N and √ N , respectively, when both i and N − i are comparable to N . Assume that the system is on open boundary condition, which is unimportant to the counting as N → ∞. After qN steps, suppose that there are i pairs of lineons. Each pairs of lineons correspond to a line segment and the number of the configuration of the length of the segments is S i qN , where S n r means the number of ways to put r identical items into n categories in which there are at least one item. The number of ways to distribute the remaining sites between the segments is S i+1 N −qN +2 . The number 2 is present since the head and the tail can contain no unchosen segments. Since the combinatorial operator S can be calculated as S r n = C r−1 n−1 , For the last approximate sign, Stirling's approximation is applied and terms with order smaller than 1/N are omitted. The integrand is sharp at its maximum and can be approximated by a Gaussian integral. The maximum can be obtained by finding the root of the derivatives of the logarithm of the integrand. Mathematically, let x = i max /N , where i max is the place where λ(q) is maximum: where the last equal sign is correct since the third term dominates as N → ∞. Keeping the second order term, the integrand can be written as where A = 2 log 2π − N + (1 + N ) log N + 2qN log q + (2 + 2N (1 − q)) log(1 − q) and B = N/(2(1 − q) 2 q 2 ). Therefore, the final result of λ can then be acquired The variance σ(q) can be similarly calculated In addition, the average number of pair of lineons λ and average steps i over all configuration can also be calculated. and where the number after the semicolon denotes the step of the increment.

High order terms with infinite r
Given a random walk network for a specific N → ∞. For infinite r, the growth of the reward as r increases are dominated by the middle process that moves forward and backward a larger number of times. The result of the growth of the average reward can be proven to be 1/λ = 4/N . To capture the most essential elements, consider an one-dimensional network of random walk with sites s i , i = 0, 1, 2, · · · , N and is placed from left to right. The probability of going right is (N − i)/N and left is i/N for each site i. The probability p i of finding the random walker at i after a long time satisfy which can be solved by p i = C i N /2 N . Note that the probability is sharply peaked at i max = N/2. Then if the reward are written in following expansion centered at i max -r i = 1 + j r (j) (i − i max ) j /N j , the expected value of the geometric mean of the rewards after walking a long time is r pi i = 1 + r (2) /(4N ) + O(N −2 ). Indeed, the problem posted is related to a 2D random walk problem similar to Fig. S6(a), except that the random walker can move both forward or backward. The probabilities are replaced by (leading order in N ) The probability p i,j of finding the random walker at L i,j can be solved and is sharply peaked at i max = N/4 and j max = N/2, and the reward expanded about N/4 is Therefore, the geometric mean of the reward after a long time is 4/N × (1 + O(N −1 )), where the coefficient of N −1 can be found by carefully considering the probabilities and the rewards, but is unimportant in this discussion. From here, one can see clearly that as r increases by 1, the growth of the average reward is multiplied by 4/N . Now that the size of the set S (r=ζN ) N can be calculated as follows which is evident when assuming the parity of the occurance of each site is independent except that the total parity matches N (1+ζ). Putting everything together, the growth of the matrix element at order r = ζN is N ×4/N ×v/w = 4v/w, so that when v < v c = w/4, the matrix element is asymptotically zero. Therefore, the total perturbation can be written in the following form where the coefficient C r increases subexponentially so that as v < v c = w/4, the matrix element is zero.

C. Energy splitting of H eff
Let H eff be the effective Hamiltonian of the system in the subspace of ground states, which can be written in the following form where H (r) is a matrix independent of v, and the norm of its elements satisfy H (r) eff,ij 1, and N is the size of the system so that the total lattice sites are of the order N 3 . Then the ground state degeneracy is of the order 2 N and since all ground states are on equal footing, the Hamiltonian can be written in the following form where h (r,q) G 1, A q are some permutation matrices, and the total number of q is of the order N 2r . Therefore the maximum eigenvalue λ max where α r are constants of order 1. This shows that λ max H G = 0 if v < v c and N → ∞.

D. Correction due to intersections of lineons
The asymptotic behavior of the contribution of g |H eff | g , where |g and |g are separated by multiple global line operators depends on the energy w i of two lineons from different directions meeting at a single point. Here we specifically discuss the case that w i < w, where the critical point should be corrected. To model the correction, consider the case that |g and |g are separated by a full sheet intersecting line operators. Let g H βi , where β is a configuration of successive intermediate states to connect two ground states, F ( ) βi are the corresponding energy in the unit w after i th action of the operators and the primed F means that the energy are calculated with the assumption that lineons in different directions are independent, and f β,i = F βi /F βi . The average correction factor at order r writes As mentioned in the main text, the factor f β,i is assumed to depend on F βi only. Note her that for sufficiently small F βi , the average number of the F βi endpoints that meet perpendicularly is proportional to F 2 βi , so the correction factor satisfies f βi ∼ F βi / F βi − kF 2 βi /N 2 ∼ 1 + kF βi /N 2 for some constant k. Therefore, only when F βi is comparable to N 2 , the correction factor f (r) F β i starts to affect the critical point. So we consider the high order terms with infinite r, which have largest proportion of β i to have F βi ∼ N 2 . To get a quantitative result, we first calculate F βi,dom that dominates the contribution of g H (r=∞) eff g , and then get the corresponding average correction factors We start by considering the matrix M that specify the random walk problem is the probability of the system with i segments to increase (remains unchange/decrease), and δ = 1/(2N 2 ). Therefore, the asymptotic behavior of the sum of higher order terms are captured by the matrix , where denotes the element-wise multiplication since a successive operations of the matrix S is equivalent to the calculation of the average reward of a random walk characterized by the probability matrix M and the reward matrix 1 √ ij ij . For sufficiently large N , the square root of the matrix S is The eigenvector of S 1/2 with largest eigenvalue describes the F βi that dominates the contribution of the summation of high order terms with infinite r. Under the similarity transformation defined by a diagonal matrix The eigenvector problem can be translated into a eigenfunction problem with a differential equation ]. An approximation is made by assuming j = i and it is good except at x ≈ 0, at which j − i is comparable to i, and the true denominator of the first term in the region x ≈ 0 should be smaller. The solution of the above differential equation can be approximated by f (x) ≈ Ai 2 2/3 N 4/3 x + x 0 , where Ai(x) is the Airy function, which is the solution of the differential equation f = xf and x 0 is the first zero of the Airy function. The approximation comes from the assumption that 1/ √ 1 − x is roughly linear except at large x. Therefore, the eigenfunction of S 1/2 is where the proportional factor comes from the backward similarity transformation. As N → ∞, the root of g (x) can be found by comparing coefficient of the series expansion of g(x) at N → ∞, which ultimately simplifies to log 1 The equation has a solution at x ≈ 0.26367, so the dominating F βi in the high order contribution with infinite r is The average correction can then be calculated, by assuming that at any cross, each endpoint has x dom = F dom /N 2 probability to appear. Then the correction factor of the dominating F βi is where w i is the energy of two intersecting lineons. Therefore, the random walk calculation is enhanced by a factor f dom for each order, and so the critical point is reduced by a factor of f −1 dom , i.e. v = v c1 = f −1 dom × w 4 . For X-cube model, f −1 dom ≈ 0.87. We estimated the errors for the X-cube model below -There are three types of errors (i) the average correction factor are calculated using arithmetic mean of F βi among all β i that have the same F βi , whereas the correct mean should be harmonic. The discrepancy can be checked for small N , where for N = 4, the discrepancy is about −0.02 of f −1 dom . (ii) The approximations employed in developing and solving the different equation tend to overestimate F dom . This gives a discrepancy of about 0.02. (iii) Since the correction factor depends on f βi , the dominating F βi is shifted when including the correction factor. This gives another 0.02. Therefore, an estimation can be given as below v c = (0.87 ± 0.04) × w/4.

S-III. THE UPPER LIMIT OF 2D AND 3D FRACTONS
The mapping to random walk problem highly facilitate the calculation of the perturbation series as shown in section VI. A similar mapping can also apply in the case of 2D or 3D fractons, where there are two important steps to identify the critical point. (i) Identify how the leading order perturbation grows as N → ∞. This encodes also how the higher order perturbation with finite r grows. (ii) With fixed N , identify how the higher order term with infinite r grows. Interestingly, the sole consideration of (ii) provides a upper limit in the case of 2D and 3D fractons. It can be seen as follows: suppose the leading order term can be written in the following form where N l is the minimum size of a membrane operator that connects two ground states |g and |g , α is the growth of the leading order term as N → ∞ and β is a function of N l that is subexponential. If the growth of higher order term with infinite r is known, the matrix element between the grounds state is where β , β are subexponential function of r, α is the growth of the higher order term of infinite r. From here, it can be seen that as long as r is sufficiently large comparing to N , the last term in the summation must dominate so that whenever v > w/α = v (u) c , the matrix element cannot be zero, and v (u) c is the upper limit of the critical point. As shown in section VI, the growth of the higher order term with infinite r is N/F d × v/w, whereF d is the average number of 2D (d = 2) or 3D (d = 3) fractons over all configuration. Therefore the upper limit for the two cases is where d is the dimension of the fractons andF d,r0 is the average number of dd fractons at some particular point r 0 . The factor 2 d is present since the energy of a group of 2 d dD fractons is defined to be w. For the case with 2D fractons, pick any point r 0 on which a fracton can define. LetÔ r1 ,Ô r2 ,Ô r3 andÔ r4 to be the four surrouding operators listed in clockwise direction. The number of fractons at r 0 is where s 5 = s 1 and s i = (1 +Ô ri )/2. Therefore, the average number of fractons over all possible configuration is F 2 = 3/4 so that the upper limit of the critical point is 3w/16. For the case with 3D fracton, pick any point r 0 on which a 3D fracton can define. Let r i = r 0 + l× ((−1) i1 , (−1) i2 , (−1) i3 ), where the binary representation of i − 1 is i 1 i 2 i 3 and l is smallest possible length such thatÔ ri is well defined. The number of fractons at r 0 is where the binary representation of j is j 1 j 2 j 3 j 4 j 5 j 6 j 7 j 8 , and the values of a j 's are recorded in Table S1. Therefore the average number of fractons isF and the upper limit of the critical point is 256v/29w.

S-IV. LOWER LIMIT OF THE CRITICAL POINTS
The main idea is to calculate the leading order perturbation H (N ) eff , where N is the size of the system and then estimate how the matrix element grows as N → ∞ to get the critical point. For lineon excitation, the leading order perturbation can be written as where S N is the permutation group of order N , a p is the number of pairs of lineons multiplied through the process characterized by p. (e.g. N = 6, p = (1, 3, 5, 2, 4, 6), then the numbers of pairs of lineons in each process are 1, 2, 3, 2 and 1, so that a p = 12.). As shown above, this value can be calculated exactly: so that the critical point calculated in this way is v c = w/4. For fracton excitation, consider L 1 numbers of lines of length L 2 packed together. We conjecture that the growth of the matrix element is bounded by the following function where the notation means that the inequality is satisfied for large N . Write where a p is the multiplication of the number of pair of lineons by treating all horizontal lines as independent, a p is the corresponding number by treating all vertical lines as indepedent, and b p is the number of tetrad of fractons, then the last term is almost smaller than n! except a few terms. In addition, we observe that a p is only weakly related to a p so that v v w n−1 −1 Numerical simulations of small L 1 and L 2 is shown in Fig. S7. The ratio R of F 2,L1L2 (w L1L2−1 /v L1L2 ) to f 2 (L 1 , L 2 ) is exponentially decaying as L 1 or L 2 increases. If the conjectured inequality is correct, when L 1 = L 2 = N , f 2 (L 1 , L 2 ) ∼ 16 L1L2 = 16 N 2 is its asymptotic behavior. Therefore, the lower limit of the critical point for the case with 2D fractons is v l c2 = w/16. In the same spirit, we conjectured a similar inequality for 3D fractons: where L 1 , L 2 , L 3 are the width, length and height of the body operator. The asymptotic behavior of F 3,L1L2L3 is bounded by f 3 ∼ 64 L1L2L3 = 64 N 3 as L 1 = L 2 = L 3 = N , so the lower limit of the critical point for the case with 3D fractons is v l c3 = w/64. R is the ratio of F2,L 1 L 2 (w L 1 L 2 −1 /v L 1 L 2 ) to f (L1, L2)