Evolution of Dynamical Signature in the X-cube Fracton Topological Order

As an unconventional realization of topological orders with an exotic interplay of topology and geometry, fracton (topological) orders feature subextensive topological ground state degeneracy and subdimensional excitations that are movable only within a certain subspace. It has been known in the exactly solvable three-dimensional X-cube model that universally represents the type-I fracton orders, that mobility constraints on subdimensional excitations originate from the absence of spatially deformable string-like operators. To unveil the interplay of topology and geometry, in this paper, we study the dynamical signature in the X-cube model in the presence of external Zeeman fields via large-scale quantum Monte Carlo simulation and stochastic analytic continuation. We compute both real-space correlation functions and dynamic structure factors of subdimensional excitations (i.e., fractons, lineons, and planons) in the fracton phase and their evolution into the trivial paramagnetic phase by increasing external fields. We find in the fracton phase, that the correlation functions and the spectral functions show clear anisotropy exactly caused by the underlying mobility constraints. On the other hand, the external fields successfully induce quantum fluctuations and offer mobility to excitations along the subspace allowed by mobility constraints. These numerical results provide the evolution of a dynamical signature of subdimensional particles in fracton orders, indicating that the mobility constraints on local dynamical properties of subdimensional excitations are deeply related to the existence of fracton topological order. The results will also be helpful in potential experimental identifications in spectroscopy measurements such as neutron scattering and nuclear magnetic resonance.


I. INTRODUCTION
The exploration of novel phases of matter beyond the Landau symmetry breaking paradigm is one of the major themes in modern condensed matter physics. Within numerous important findings, topological orders that feature long-range entanglement and their experimental detection in correlated materials and numerical simulation in quantum lattice models, have attracted a lot of attention [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15]. Topological excitations are one of the defining features of topological orders, such as anyons in the celebrated toric code model [3] and its frustrated magnet realizations [8,11,[16][17][18][19][20][21][22][23]. While a single topologically trivial excitation can be created by a local operator, point like topological excitations, must be created in pairs, which are located at the two endpoints of string-like operators. As endpoints are exactly the boundaries of open strings, the topological excitations are inherently related to string operators. Since the string operators can be arbitrarily deformed, the associated topological excitations can move freely in the whole space.
Recently, fracton orders as an unconventional realization of topological orders were proposed and have stimulated intensive activities in many areas ranging from condensed matter physics and quantum information science, to high energy physics . Exactly solvable models that support fracton orders (phases) are usually stabilizer codes defined in three and higher dimensions. In sharp contrast to the conventional topological orders, locally-indistinguishable ground state degeneracy (GSD) in fracton system grows with the system size, such interesting phenomena can be traced back to the foliation structure built in the lattice models of three [45] and higher dimensions [46,47]. Furthermore, while topological excitations in higher-dimensional conventional topological orders are characterized by topological properties, e.g., topological field theoretical description and topological braiding data among particles and loop-like excitations , topological excitations in fracton orders are unexpectedly exotic due to mobility constraints, which can be seen from the fact that spatially deformable string operators do not exist [27] in exactly solvable fracton models, e.g., the X-cube model [26] and the Haah's code [27] in three-dimensional (3D) lattice. In the literature, topological excitations in fracton orders are dubbed fractons and subdimensional particles according to the degree of their mobility [25], with fracton denoting excitations that are completely immobile (0-dimensional), and subdimensional particle denoting excitations that are still mobile within a subspace including 1-dimensional lineons moving along straight lines (e.g., straight lines formed by nearest links along the x direction) and 2-dimensional planons moving within flat planes (e.g., all flat planes formed by nearest plaquettes that are parallel to xy planes). In this paper, for the sake of convenience, we also include fractons as subdimensional particles. Interestingly, the spirit of this nomenclature has been generalized to higher dimensions, where spatially extended excitations are subject to complicated constraints on both mobility and deformability [46,47]. In addition, as a common feature of stabilizer codes, all these exactly solvable fracton models have exactly zero correlation length and triv-transformation (pCUT) method, the ground state phase transitions of the X-cube model under the single field perturbation have been confirmed to be first-order [55,56], which is tightly connected to immobility nature of fracton excitations. The nature of the transitions is also systematically studied by generalizing the Z 2 X-cube model to its Z N cousins in terms of the tensor network state representations [52]. Moreover, dispersion relations of single quasiparticle excitations are perturbatively computed [56].
Here, with the QMC+SAC scheme, we systematically compute real-space density-density correlation functions and frequency-momentum-dependence of spectral functions (dynamical structure factors) of subdimensional particles. From the QMC numerical results, we find that mobility constraints on subdimensional particles deeply influence real-space correlation functions and spectral functions in many aspects. Starting from the zero field limit where all subdimensional particles are dispersonless, our numerical results indicate that, in the presence of quantum fluctuations generated by the Zeeman fields, the subdimensional particles acquire dispersive dynamics strictly under mobility constraints. After the quantum phase transition where the fracton phase is turned into the trivial paramagnetic phase, the anisotropic features under mobility constraints disappear. In addition to QMC simulation, we also perform a mean-field+RPA (random phase approximation) analysis to fit the QMC results of spectral functions, which shows a qualitative consistency with the QMC results.
Our QMC results provide a set of dynamical signature of subdimensional particles in the type-I fracton orders, which clearly demonstrates an exotic interplay of topology and symmetry in fracton physics and shows sharp distinction from the conventional topological orders. As spectroscopy measurements, e.g., inelastic neutron scattering and nuclear magnetic resonance, play important roles in identifying novel phases of matter in strongly correlated materials, we expect the QMC results in the present work will be helpful in experimentally identifying fracton orders in real materials. In addition, the QMC results also present various features in correlation functions and spectral functions, which are yet to be fully understood, so we expect our numerical study will stimulate further studies in the field of fracton physics from both theoretical and numerical aspects. This paper is organized as follows. In Sec. II A, we introduce the theoretical background of the X-cube model and the density operator of three different subdimensional excitations: fracton, lineon and planon. Theoretical expectation on correlation functions and dynamical structure factors is given in Sec. II B. In Sec. III, QMC+SAC results are presented in details, where the comparison between QMC and meanfield+RPA analysis is also provided. We summarize the paper in Sec. IV. Appendix A contains several technical details of our numerical method. The Hamiltonian of the X-cube model with external Zeeman-like magnetic fields is expressed as [26,[54][55][56]58] FIG.
1. An illustration of the stabilizer operators in the X-cube model. A c,i denotes the product of all σ x in the each cube (the red shadowed cube), while B v,i denotes the product of σ z in the each cross (the red shadowed planes).
As shown in Fig. 1, spin-1 2 s are defined on links, and A c,i (B v,i ) denotes the product of all the σ x (σ z ) in the each cube (cross), respectively [26,29,54,58,112,113]. Here the summation index i in general refers to a spatial location, while its specific meaning depends on the summed operator, for example, in A c,i (B v,i , σ z i ) i refers to the coordinates of the center of a cube(vertex, link). As in this paper we mainly consider the distance between operators of the same kind, this ambiguity should not lead to any misunderstanding. There are three different kinds of subdimensional excitations in the X-cube model, which are lineon, fracton, and planon. A lineon corresponds to the flipped eigenvalues of two different B v,i terms at the same site, which combinatorially leads to three different species of itself. And a fracton corresponds to the flipped eigenvalue of the cube operator A c,i , which are completely immobile. A pair of fractons along the x (y, z) direction correspond to a planon which can move in the 2D plane perpendicular to x (y, z) direction.
We can obtain mobility constraints by considering the creation operators of excitations. First, fractons are generated by membrane operators of the form W(M) = i∈M σ z i , where M is a square membrane composed of parallel links, i ∈ M means that link i within M (see Fig.2(e)). Such a membrane operator generates 4 fractons respectively located at the four corners.
Therefore, moving such a fracton with another membrane operator (i.e. apply a membrane operator to annihilate the original fracton, and create another one at another site) will generate 2 extra fractons, that costs additional energy. Meanwhile, a planon, as a bound state of two fractons, is mobile along a two-dimensional manifold (see Fig. 2(f)). As for lineons, that are created by string operators of the form W(S ) = i∈S σ x i . Here, S is a straight string composed of links (see Fig.2(d)). Similar to fractons, moving a lineon at an endpoint of S with another string operator perpendicular to S will create an extra lineon at the intersection point of the two strings. Equivalently, we can say that for a W(S ) = i∈S σ x i operator, where S is an arbitrary string, lineons are not only generated at its endpoints, but also its turning points. In brief, any movement perpendicular to the direction of S is energetically unfavorable.
To investigate all these excitations in the QMC+SAC scheme [8,9,103,108,114,115], we compute the real-space correlation and their spectra (obtained from the SAC of dynamical correlation functions) in the X-cube model in both σ x and σ z basis with linear system size both L = 6 and L = 10 (the total size is N = L 3 ) and β = 2L. For the dynamic correlation, the imaginary-time correlation functions in the lineon channel ( Fig. 2(a)) Fig. 2(d)); the fracton channel ( Fig. 2(b)) with n f,i = 1 2 (A c,i − 1), where fractons are created by W(M) = i∈M σ z i (Fig. 2(e)), and the planon channel ( Fig. 2(c)) with n x,i = 1 4 (A c,i − 1)(A c,i+x − 1), where planons are created by W(M p ) = i∈M p σ z i ( Fig. 2(f)),, are measured in the QMC simulation. Here, operators O x , n f and n x denote the number operator of lineon, fracton and planon [30], respectively. Since the ground state by definition should be vacua of all excitations, one can verify that the expectation values of n f , O x and n x indeed vanish in the X-cube model's ground state. However, under the perturbation of Zeeman couplings h z and h x , their expectation values would become non-zero caused by noncommutation between A c,i (or B v,i ) and Zeeman fields, which corresponds to the emergence of subdimensional excitations (see also Appendix. A 3).

B. Correlation and spectral functions influenced by mobility constraints
Theoretically, we expect correlation functions to be an important way to characterize the fracton order of the X-cube model with perturbative external field. To see this point, we can compare the X-cube model with 3D toric code model (TCM) [3,116,117]. As a model with pure topological order, in 3D TCM, a pair of point-like charge excitations is generated at the 2 endpoints of a string operator, and a loop-like flux excitation is generated on the boundary of an open membrane operator. That is to say, the relation between creation operators and excitations in 3D TCM is completely topological: excitations are always generated at the boundaries of creation operators. While as a fracton order model, in the X-cube model, as discussed in Sec. II A, there is no such completely topological relation between excitations and creation operators. For example, in the X-cube model, fractons are generated at the corners of a membrane, and lineon are not only generated at the endpoints, but also at the turning points of a string. As a result, fracton orders show a complicated interplay between topology and geometry. Therefore, even though mobility restriction is resulted from the fracton order of the whole system, the correlation function should also contain information of mobility restriction, and such restriction would have intricate ramifications into the three types of subdimensional excitations.
First, we consider lineon excitations. The lineon-lineon correlation at long distances is expected to be anisotropic due to the one-dimensional mobility of a lineon. In real space, for lineons mobile along the direction denoted as x, the lineonlineon correlation O x,i (τ)O x, j (0) should decay more slowly along the mobile x-direction. In momentum space, we expect such anisotropy of correlation function G O x (q, τ) to be present around the Γ point, as a result of the mobility restric-tion. Then, we consider planon excitations. As demonstrated in Sec. II A, a planon is composed of a dipole of adjacent fracton excitations. There is no clear evidence that guarantees an attractive interaction between adjacent fractons. Thus, whether or not the composite is a well-defined point-like excitation at long wavenlengths is hard to be confirmed. Due to the existence of intrinsic dipole-like structure of a planon, planon-planon correlation function should be anisotropic not only in the fracton ordered phase, but also in the paramagnetic phase with strong external field h z . Hence, anisotropy of planon-planon correlation is not necessarily related to fracton orders. At last, we consider fracton excitations. As a fracton is totally immobile, the fracton-fracton correlation function at long wavelengths is expected to be isotropic. Therefore, fracton-fracton correlation is relatively simple, and we do not expect to read much information of fracton orders from such correlation functions. In addition to the correlation function, we can also expect the spectral functions of subdimensional particles to show dispersive behavior along their movable directions.
As the X-cube model itself is isotropic, the anisotropy of lineon-lineon correlation at long wavelengths suggests a hidden rotational symmetry breaking may be related to the type-I fracton order. As a concrete example, the X-cube model is obviously invariant under a 90 • spatial rotation with respect to xaxis, while the correlation functions of lineons mobile along y (and z) directions do not have this symmetry. In contrast, for pure topological orders, such as the isotropic 3D TCM, we expect the charge-charge correlation to be also isotropic. Therefore, such anisotropy of lineon-lineon correlation (or an isotropy breaking) may be recognized as a signature of Type-I fracton orders. Nevertheless, in type-II fracton orders, such as the Haah's code [27,29,118], we do not expect to find such spatial symmetry breaking due to the lack of subdimensional mobile excitations.

III. NUMERICAL RESUTLS AND ANALYSIS
In this section, we present our QMC+SAC results of lineon, fracton, and planon excitations in the X-cube model with external Zeeman fields.

A. Numerical setting and critical fields
Numerical simulation of fracton systems is notoriously difficult and arduous. We first begin with a few technique details which are related with the fracton physics. We implement our QMC simulation in the framework of stochastic series expansion (SSE) [103], where the sampling configuration is constructed by Taylor expansion of the partition function in a chosen basis. In our case, it is {σ z } for lineon and {σ x } for fracton in the simulation. Unfortunately, the four (B v,i ) and twelve-spin (A c,i ) interactions in the X-cube model cause a glassy and fragmental configuration space, which leads to a low efficiency in the Monte Carlo sampling. For example, on one hand, the local update is not efficient for QMC simulation of the X-cube model, due to the extended interactions in the Hamiltonian. On the other hand, naively using the cluster update algorithm [104,119], the cluster would rapidly extend due to the B v,i and A c,i interactions and an extremely large area of spin is always suggested to be flipped, but flipping these extremely large areas has the same effect as flipping few spins. To solve this problem, we firstly modify the cluster update algorithm to slow down the cluster growth, and apply it along with the local update in our simulation. Secondly, for system size L = 10, we warm up our Markov chains from an initial configuration from an equilibrium configuration in L = 2 QMC simulation (spatially repeated from L = 2 to L = 10), to help the larger size simulation to thermalize quickly. The details of our QMC update scheme are given in Appendix A 1, A 2 and A 3. Our numerical simulations are carried out under the periodic boundary condition with linear system size L = 6 and L = 10, and total size N = L 3 and inverse temperature β = 2L.
Note that the ground state of the X-cube model is a vacuum state of all three subdimensional excitations, which means O x , n f and n x all commute with the X-cube Hamiltonian Eq. (1) with h x = 0 and h z = 0. Therefore, to observed these subdimensional excitations, the external perturbation fields are required to introduce a fluctuation to the corresponding operator (like O x for lineons). Here, we apply the h x transverse field to see lineon excitation and h z for fractons and planons. As described in Fig. 3(a), such perturbations also lead to a first-order phase transition from the fracton phase to the paramagnetic phase, which can be viewed as a consequence of the immobility of fracton excitation [52,55,56]. To compute correlations of these excitations with better data quality, beside the improvement in the QMC update scheme discussed above, we further utilize the quantum annealing (QA) algorithm [120][121][122], in which quantum parameter, h x(z) , is slowly changed with the annealing step ∆h x(z) = 0.01 for a faster convergence to the optimal state. In the fracton phase, starting from the exactly solvable point (h x = 0, h z = 0), our simulations scan the parameter towards (h x = 0.9, h z = 0) for the lineon and (h z = 0.3, h x = 0) for the fracton cases, where we apply 2 × 10 5 Monte Carlo step for each annealing step. And for measurements in the paramagnetic phase, our simulations scan from the paramagnetic limit in x(z)-direction [PLx(z)] point towards h x = 0.9 for the h x case and h z = 0.3 for the h z . The detailed QA implementation is given in Appendix. A 3. In Fig. 3(b,c), the energy per spin e are measured in both the fracton phase and the paramagnetic phase with h x and h z field on a L = 10 system, respectively. These first-order phase transitions occur at h x ≈ 0.9 in the h x perturbing case and h z ≈ 0.3 for the h z case.

B. Correlation functions
To show mobility constraints on subdimensional excitations, we first measure their real-space normalized correlation functions defined as C n x (r) = n x,i n x,i+r − n x,i 2 n 2 x,i − n x,i 2 .
As discussed in Sec. II B, the anisotropic properties of C O x and C n x are tightly related to their own mobility in each direction. In Fig. 4, we present the numerical results of these correlation functions on a L = 10 system. The lineon-lineon correlation function, C O x (r), is anisotropic in the fracton phase. Taking O x with h x = 0.8 in the fracton phase as an example, C O x (r) decays more slowly along the x-direction with the neighboring correlation C O x (x) = 0.421(3), but fast along the yand z-direction with C O x (y) = 0.0003(6) and C O x (z) = 0.00002 (6), which is shown in Fig. 4(a-c). Such a property comes from the mobility constraint that lineons can only move along x-direction. And in the paramagnetic phase, the fracton order is completely suppressed and it will be inappropriate to interprete O x as lineon density. Therefore, the restriction of the lineon mobility also disappears, and so does the anisotropic property. As a result, C O x (r) rapidly decays to zero along all directions, which is given in Fig. 4(d-f).
Although the fracton order can also be suppressed by the h z field, the real-space correlation function of planons, C n x (r), and fracton, C n f (r), behave differently from their 1dimensional cousin due to their different mobility nature. For the planon, which is mobile in the y − z plane, its correlation function C n x decays slowly in both yand z-direction with neighboring C n x (y) = 0.251(6) and C n x (z) = 0.238 (8), but fast in the x-direction C n x (x) = 0.007(2) [see Fig. 4(m-o)]. Such a restriction of mobility is broken when the model enters the paramagnetic phase. But, different from the lineon case, the planon operator n x can be viewed as a composite of two neighboring fracton operators n f in the x-direction. Consequently, C n x (r) remains anisotropic even in the paramagnetic phase. However, it slowly decays along its intrinsic structure direction (x-direction for n x ), which is shown in Fig. 4(p-r) with neighboring correlation C n x (x) = 0.296 (8). When it comes to fracton, it is immobile in all directions. This restriction leads to the isotropic correlation function in both fracton and paramagnetic phases, which looks featureless compared to the correlation data in the lineon channel. However, since the fractons are created at four corners by a membrane operator, C n f (x) = 0.324(6) at the neighboring point in Fig. 4(gi), which are C n f (x) = 0.020(1) in the paramagnetic phase [ Fig. 4(j-l)].
In Figs. 5, 6 and 7, the first rows are in the fracton phase, with external fields h x = 0.7 − 0.9 in Fig. 5(b-d) and h z = 0.1 − 0.3 in both Fig. 6(a-c)  . From the results, we see that, for any given momentum point, the spectral functions in all three channels do not have signals when frequency is below a finite value. Therefore, the spectral functions have finite gaps in the fracton phase, meaning that the density fluctuations of all three kinds of subdimensional excitations must cost finite energy. This is expected, since all subdimensional excitations have finite single-particle gaps. Furthermore, in Figs. 5(b-d), 6(a-c), and 7(a-c), we find that, upon increasing external fields, quantum fluctuations get stronger and stronger, which causes more and more pronounced spectral signals (quantified by brightness in the ω − k plane) in the fracton phase. And in the paramagnetic phase shown in the second rows of Figs. 5, 6, and 7, the gaps for density fluctuations become larger and larger upon increasing external fields. We must stress that, since subdimensional particles are no longer well-defined excitations in the trivial paramagnetic phase, rigorously speaking, the spectral functions should no longer be interpreted as measurement of density fluctuations of subdimensional particles.
Let us first focus on the spectral function of lineons. The peak energy (i.e., the frequency of the brightest signals) of the spectral function of lineons is approximately given by ∆ ≈ 8 in both the fracton phase [ Fig. 5(b-d)] and the paramagnetic phase [ Fig. 5(e-g)]. However, the profile of the spectral function in the fracton phase clearly shows a dispersion behavior along the Γ → X, while it is flat in the paramagnetic phase. Theoretically, this dispersive behavior of density fluctuations can be traced back to the free mobility of lineons along xdirection. Moreover, in the paramagnetic phase, the profile of the spectral function is flat over the whole momentum space, but its peak energy becomes larger and larger upon increasing h x from Fig. 5(e) to (g).
Next, we perform a mean-field analysis on the spectral functions of lineons. Since the kinetic energy of lineons along dispersionless directions is strictly suppressed such that interaction effect is no longer negligible, which motivates us to focus on dispersive directions. Without loss of generality, we consider lineons mobile along x direction. At first, we regard O x,i as a density operator of lineons on site i. Then, we need to notice that the ground state of X-cube model perturbed by Zeeman fields is no longer a "vacuum" for lineons but contains finite density of lineons. Furthermore, by noticing that lineons are hard-core and mobile within one-dimensional subspace, and we expect the difference between the exchange phases of hardcore bosons and fermions to have little influence when there is few particles in 1D (thus particles can hardly exchange positions), it is phenomenologically reasonable to model lineons as 1D fermions with both positive and negative energy spectra where all states with negative energy are filled. Therefore, we propose a density-density correlation function [123,124] of lineons of the following form to fit the QMC results of spectral functions: Here, ξ ν k is a gapped dispersion relation. The indices ν, ν = ± and we require that ξ + k = −ξ − k and ξ + k > 0. n F (ξ) = (e βξ + 1) −1 is the Fermi-Dirac distribution, and β is the inverse temperature. V is the volume of the system and η is an infinitesimal real number. Besides, compared with the result in Ref. [123,124], a k-dependent matrix element has been omitted, since here only the qualitative tendency is of concern. Then, at zero temperature, χ 0 (q, ω) can be reduced to , where the second term can be safely neglected when ω > 0. Since the lineons we consider here are only mobile along x direction, the volume V reduces to the linear size L upon summing over k y and k z . To proceed further, we assume the lineon dispersion to be: x is taken from the perturbative result in Ref. [56]. By further phenomenologically considering a RPA approximation, i.e., to effectively include the interaction effect (J is a coupling constant), the dynamic structure factor (spectral function) can be obtained as A(q, ω) = −Imχ RPA (q, ω). Here, we set parameters η = 0.01, L = 2000, J = 14. We compare the RPA results and QMC results in Fig. 8, which shows a qualitative consistency between the two approaches. The QMC results in Fig. 8 are carried out on a L = 6 system also annealing from the exactly solvable point but with more bins in the QMC measurement, which help us to reduce the error in G O x and solve out a high-quality spectral function A O x . While a more standard slave-particle mean-field treatment [125,126] is desirable, we leave it to future investigation. As lineons, fractons, and planons are localized along some certain directions, it will be also interesting to probe such localizations. One of possible ways is to compute complementary entanglement properties from both real space cut and momentum space cut via techniques and observables in e.g., Refs. [127][128][129][130], which we also leave to future investigations.
After analyzing the spectral function for lineons, let us move to the spectral functions of fractons and planons. Figs. 6 FIG. 8. Comparison between the dynamic structure factors in the lineon channel obtained by mean-field+RPA method and the QMC process. The momentum range is (X → Γ) along which lineons are dispersive. The RPA results are carried out in with coupling constant J = 14. Fig. 8(a-c) are in the fracton phase with the external field h x raising from 0.7 to 0.9. Here, the QMC results are simulated at system size L = 6. and 7 respectively show the spectral functions of fracton and planon excitations with h z changing from 0.1 to 0.6 on a L = 10 system, respectively. Since h z term does not commute with the A c,i , it leads to a fluctuation of A c,i and causes the appearances of fractons and planons. For small perturbation (h z ≤ 0.3), the X-cube model stays in the fracton phase. Figs. 6(a-c) shows the fracton spectral function in the fracton phase, and Figs. 7(a-c) are that of planon. In panel (a) of Figs. 6-7, the perturbation from h z is too weak and the dynamic signal of the fracton and planon is hard to be observed. By increasing h z , the stronger fluctuations in n f and n x make their corresponding spectral signal more significant from panel (a) to (c). The profiles of spectral functions of both fractons and planons are flat with peak energy ∆ ≈ 8, which is approximately four times of single fracton energy [56]. When h z ≥ 0.3, the X-cube model enters the trivial paramagnetic phase. In this case, the profiles of spectral functions of both fracton and planon channels remain flat but the peak energy is becoming larger with stronger perturbation, which are presented in Figs. 6(d-f) and 7(d-f).
In summary, the most significant feature in the fracton phases in Fig. 5-7 is that the dispersion behaviors of these spectral functions show strong anisotropy corresponding to their mobility restriction. For example, in the Fig. 5(b-d), the profile of the spectra of lineon in the fracton phase clearly shows a dispersive behavior along the Γ → X, while it is flat in the paramagnetic phase (Fig. 5(e-g)). Meanwhile, with the help of the mean-field+RPA method, we can reproduce the profile of the lineon spectra along the dispersive direction as well. Besides, for fracton spectra as shwon in Fig. 6, we see no obvious anisotropy, that is also consistent with the totally immobile feature of fractons. As a result, we believe the anisotropic behavior is a characteristic feature of Type-I fracton orders, where partially mobile excitations exist.

IV. SUMMARY AND OUTLOOK
Before closing this paper, let us briefly summarize our numerical findings. In the lineon channels, by means of the quantum fluctuations from the h x field, lineons are excited, whose particle number becomes finite and dense in the ground state. In the fracton phase, the real-space lineonlineon correlation shows a significant value on the neighboring points along its one-dimensional mobile direction, e.g., the x-direction for O x . Meanwhile, the peak profile of the spectral function of O x also presents a dispersive behavior along the Γ → X direction. With the help of the mean-field theory and RPA method by assuming an anisotropic dispersion relation of lineons, we qualitatively recover the QMC results of the spectral function of lineons, which demonstrates an exotic connection between spectral functions and mobility constraints on dynamics of lineons. By comparing with the results in the trivial paramagnetic phase, these anisotropic properties can serve as an experimental signature of the fracton order.
In the fracton and planon channels, under the perturbation of h z field, the fluctuation of A c,i terms leads to the appearance of both fracton and planon excitations. The fracton-fracton correlation function is isotropic in both fracton and paramagnetic phase, which is due to its completely immobile nature. However, in the plannon case, its correlation function behaves anisotropically in both phases but for different reasons. It is due to its mobility constraint in the fracton phase but the intrinsic dipole-like structure in the paramagnetic case. Such a difference leads to a significant geometrical rotation of these anisotropic features, which could be a very experimental signal for the purpose of identifying the type-I, i.e., X-cube fracton order. And, for spectral functions, the spectral peak energies in both these channels in the fracton phase are nonsensitive against the increasing h z but become sensitive in the paramagnetic case.
In conclusion, with the external Zeeman-like magnetic fields pushing the X-cube model away from its exactly solvable point, we have improved QMC methods to overcome the glassy Hilbert space, and computed both real-space correlation functions and dynamic structure factors in three subdimensional excitation channels in both the fracton and the paramagnetic phases through the method of QMC+SAC in a wide parameter space. We have found how mobility constraints on local dynamical properties of subdimensional excitations is correlated to the existence of fracton topological order. This fact is nothing but the exotic interplay of topology (global topological order) and geometry (local dynamical properties) in fracton orders. As our numerical results show novel dy-namical features that are yet to be fully understood, we expect our numerical results will stimulate more efforts in both theoretical and numerical studies in fracton physics. While the X-cube Hamiltonian with magnetic fields looks intricate, it will be interesting and challenging to analytically study correlation and spectral functions in the theoretical framework of slave-particle (projective construction) techniques. Our numerical results also set a further interesting question about the connection between the single-particle dispersions of the lineon, fracton and planon channels and their density spectral functions [56]. Finally, as spectroscopy measurements have been widely performed in strongly-correlated materials, our numerical results on dynamical signature will be helpful in experimental identifications of fracton orders in spectroscopy measurements such as neutron scattering and nuclear magnetic resonance.

Diagonal update
The diagonal update is about inserting and removing the operators. In the inserting process, the whole number of possible insert positions are Here N c is the total number of the cubes and N v for the crosses. And N l is the number of the spins. We prefer a position by randomly picking a number N x in N, which determine the operator type would be inserted. If 1 ≤ N x ≤ N c , we propose to insert a H Ac,i operator at the N x cube. When and h x 0, a H h x,i operator is applied to corresponding spin. Moreover, when N c + N v + 1 ≤ N x ≤ N c + N v + N l and h x = 0, or N c + N v + N l + 1 ≤ N x ≤ N c + N v + 2N l , we propose to add a H h z,i operator on the N x spin.
After an identity operator is encountered, such insertion of the operator H D would be accepted with probability with m for the length of the configuration and n for the operator number. Noted that the value of H D not only depends on the type of operator, but also the spin states at the corresponding position.
In the removing case, we sweep the configuration and remove the diagonal operator H D with probability When it comes to the off-diagonal update process, both the local update and the modified cluster update are applied in our simulation. Firstly, there are two kinds of operators in the configuration space, which are the pure diagonal operator (H Ac,i and H h x,i ) and the quantum operator (H B v,i and H h z,i ). Due to constant term in Eq.A1, the diagonal element in H B v,i and H h z,i are no-zero. Thus, in the configuration, the quantum operator can be both diagonal and off-diagonal.
In the local update process, a leg of a quantum operator (H B v,i and H h z,i ) is selected randomly in the configuration. Then, starting from this leg, we create all the update-lines of this vertex and go along the imaginary time until finding another operator acting on the same position (see Fig. 9(a)). The spins between these two operators are proposed to be flipped (the red region in the Fig. 9(a)).
For the cluster update, also starting from a randomly picked-up vertex leg of a quantum operator, the cluster is built following the rules as below. (1) When the cluster building line meets a pure diagonal operator (H Ac,i and H h x,i ), it would go through the operator directly as Fig. 9(d) shows. (2) When the cluster building line meets a quantum operator (H B v,i and H h z,i ), it would evolve in two different ways. This line can go through the operator directly as Fig. 9(d). Or it will be reflected from the other vertex leg as Fig. 9(c). In each cluster update process, we pick 10% of these quantum operators randomly and treat them in the Fig. 9(c) way in the cluster constructing process while the others in the Fig. 9(d) ways. Within this treatment, the cluster extending would be slowed down, which make the cluster update more effective. And it is worth to note that this cluster building process would turn back into the typical cluster update with treat each quantum operators in the Fig. 9(c) way. Finally, the spins including in the cluster (the red region in the Fig. 9(b)), are suggested to be flipped.
In both of these updates, the spins included in the red region would be flipped with the acceptances probability given by Here, n +,new(old) are the number of diagonal operators with H Ac(h x ) = 4K(h x ) in the new (old) configuration. It means that the weight ratio depends on the number of the overlapvalue-changing diagonal operator.

Deducing initial configuration
In the Monte Carlo simulation, the final results should be independent of the initial configuration. However, a better choice of the initial configuration is helpful to reduce the computational cost. To find such a better choice, with the concept of super-resolution, we deduce an initial configuration of ( f L) 3 β system from an equilibrium L 3 β configuration. And here, in our case, the initial configuration of L = 10 is deduced from the L = 2 result with f = 5.
Under the framework of the SSE method, the partition function is expanded as This expansion constructs a configuration space for the Monte Carlo sampling. Such a configuration keeps two important features, the operator string n i=1 H i and the initial state |α 0 . Finding a better initial configuration by super-resolving is to find a mapping from the configuration of L 3 β to that of ( f L) 3 β, which would make fewer mutations in physic values observing from the QMC process, such as energy density and magnetization. To do so, we discuss the operator string and the initial state respectively.
For the operator string, with the periodic boundary condition, the super-resolution process can be viewed as making f 3 replicas but shifting the location of the operators. Firstly, the operators in the original configuration are marked as {τ 0 i }, and the length of the configuration is extended m → f 3 m. Then, we pick up f 3 n positions in f 3 m randomly, which are marked as {τ j }. And we place the τ 0 i operator on the positions τ f 3 (i−1)+1 to τ f 3 i for each replicas. For instance, in Fig. 10, the τ 0 1 operator is replaced on the positions τ 1 and τ 2 . In this way, we keep the fact that the operators should appear uniformly in the configuration spaces. When it comes to the initial state, noted that the ground state |GS of the X-cube model is where | ↑↑ ... ↑ denote the fully polarized state [56,112,113]. Therefore, we prefer the fully polarized state in the chosen basis as the new initial state of the larger system size configuration. Starting from this new configuration, with further Monte Carlo steps (≥ 10 4 ), a new equilibrium configuration will be reached. In our work, we applied such a process in the simulation at the exactly solvable point.
Then, due to the first-order phase transition caused by h x and h z , our simulation is also required the QA process to achieve a faster convergence to the ground state under the external Zeeman-like magnetic fields, in which the quantum parameter, h x(z) , would be exactly slowly changed and the configuration from the last parameter result would be set as an initial configuration for QMC simulation at the next parameter. Our simulation in the fracton phase of the X-cube model scans from the exactly solvable point with an annealing step ∆h x(z) = 0.01 and over 2 × 10 5 Monte Carlo steps at each annealing step [120][121][122]. And the measurements in the paramagnetic phase are from PLz(x) points in Fig. 3 with the same annealing step.
Moreover, we measure the density operators changing with the external Zeeman field, which are plotted in Fig.11. In Fig.11, the red circles are in the fracton phase while the green triangle are in the trivial phase. As it is shown, in the fracton phase, the density operators of all three kinds of subdimesional excitations are weak meaning that these subdimensional exicitations are hard to be created in the fracton phase. In contrast, these density operators have significant values in the trivial phase, which caused by the strong quantum fluctuation due to the external Zeeman fields. Meanwhile, the value of these density operators change discontinuously near the quantum critical point, which is due to the first order phase transition nature between the fracton and the trivial phase.