Renormalization group analysis on emergence of higher rank symmetry and higher moment conservation

Higher rank symmetry and higher moment conservation have been drawn considerable attention from, e.g., subdiffusive transport to fracton topological order. In this paper, we perform a one-loop renormalization group (RG) analysis and show how these phenomena emerge at low energies. We consider a $d$-dimensional model of interacting bosons of d components. At higher-rank-symmetric points with conserved angular moments, the $a$-th bosons have kinetic energy only along the $x^a$ direction. Therefore, the symmetric points look highly anisotropic and fine-tuned. By studying RG in a wide vicinity of the symmetric points, we find that symmetry-disallowed kinetic terms tend to be irrelevant within the perturbative regime, which potentially leads to emergent higher-rank symmetry and higher-moment conservation at the deep infrared limit. While non-perturbative analysis is called for in the future, by regarding higher-rank symmetry as an emergent phenomenon, the RG analysis presented in this paper holds alternative promise for realizing higher-rank symmetry and higher-moment conservation in experimentally achievable systems.


I. INTRODUCTION
The celebrated Noether's theorem [1] relates a conservation law to an underlying continuous symmetry. For example, in a U (1)-symmetric Hamiltonian of bosons, bosonic operatorΦ(x) is changed to e iθΦ (x) under symmetry transformations where the real parameter θ doesn't depend on coordinate x = (x 1 , x 2 , · · · , x d ) in ddimensional space. By means of Noether's theorem, one can show that the total boson number, i.e.,´d d xρ(x), is a conserved quantity, where particle number density ρ(x) =Φ † (x)Φ(x). Apparently, ρ(x) is just the zeroth order of conventional multipole expansions: in a standard electromagnetism textbook [2]. In a particle-number-conserving system, higher moment conservation, e.g., conservation of dipoles and quadrupoles, is in principle allowed. Furthermore, if the density is vector-like with multi-components, denoted as ρ = (ρ 1 , ρ 2 , · · · , ρ d ), then we can define another set of multipole expansions: where the third one is angular moment. ab··· is Levi-Civita symbol. If d = 2, 3, it can be rewritten in a compact form: ρ × x. In d = 2, ρ × x = ρ 1 x 2 − ρ 2 x 1 . Indeed, recently we have been witnesses to ongoing research progress on higher-moment conservation and the associated higher-rank version of global symmetry [3][4][5][6][7][8][9], especially in the field of fracton physics . Some typical examples of research include subdiffusive transport at late times, non-ergodicity, Hilbert space fragmentation, and spontaneous symmetry breaking [31,33,41,[68][69][70][71][72][73][74]. In a simple scalar theory, the associated higher-rank symmetry transformations are parametrized by θ(x) that is a polynomial function of x [8]. Inspired by the conventional correspondence between global symmetry and gauge symmetry, upon "gauging" higher-rank symmetry, higher-rank gauge fields can be obtained [5]. Here, the gauge fields are usually higher-rank symmetric tensor fields, which leads to generalized Maxwell equations [6] and exotic theory of spin systems in Yb-based breathing pyrochlores [75].
As a nontrivial consequence of higher moment conservation, the mobility of particles is inevitably restricted, either partially or completely. For example, it is quite intuitive that dipole conservation strictly forbids a single particle motion along all spatial directions. Such particles are called "fractons" or 0-dimensional particles [3,4]. Similarly, one can define lineons (1-dimensional particle) that are movable within a stack of parallel straight lines and planons (2-dimensional particle) that are movable within a stack of parallel planes. Regarding these strange particles as bosons, we can consider their Bose-Einstein condensation, such that the spontaneous breaking of higher-rank symmetry occurs. As a result, a class of exotic quantum phases of matter dubbed fractonic superfluids [68,69] is formed. In Ref. [69], a convenient notation dSF i was introduced to denote d-dimensional fractonic superfluids with i-dimensional particle condensation, e.g., dSF 0 with condensed fractons and dSF 1 with condensed lineons. The conventional superfluid phase corresponds to dSF d where bosons can freely move.
Higher-rank symmetric microscopic models often look quite unrealistic, highly anisotropic and fine-tuned [5]. For example, Hamiltonian does not has the usual kinetic energy term [68]. And interaction is delicately designed [69]. However, as we've known, in many condensed matter systems, symmetry has been found to be significantly enhanced at low energies. For example, Lorentz symmetry emerge in graphene which is microscopically built by non-relativistic electrons. Thus, one may wonder whether it is possible that long-wavelength low-energy limit will conserve higher moments and respect higherrank symmetry as an emergent phenomenon.
For this purpose, we may apply the traditional theoretical approach: renormalization group (RG). If there exists a phase region such that all models in the region can flow to symmetric points, we can regard higher-rank symmetry as an emergent symmetry. Theoretically, one advantage of such an emergent higher-rank symmetry is its robustness against symmetry-breaking perturbation. Practically, we expect that such a scenario holds promise for more flexible realization of exotic higher-rank symmetry and higher-moment conservation in both theoretical and experimental studies.
In this paper, we identify such a wide phase region that supports emergent higher-rank symmetry and conservation of angular moments, i.e.,´d 2 xρ × x =´d 2 x(ρ 1 x 2 − ρ 2 x 1 ) for a two-component boson fields in two dimensions. We start with a two-dimensional many-boson system in the normal state (i.e. without lineon condensation) of fractonic superfluids (denoted as 2SF 1 ). The Hamiltonian is a symmetric point in the parameter space where the a-th (a = 1, 2) component bosons only have kinetic terms along a-th axis (dubbed "diagonal " kinetic terms). There also exists a weak inter-component scattering term allowed by higher-rank symmetry. We shall perform a RG analysis in the vicinity of the symmetric point by adding symmetry-disallowed kinetic terms (dubbed "off-diagonal " kinetic terms) as a perturbation. The one-loop calculation of β-function shows that there exists a finite phase region (Fig. 4) where off-diagonal kinetic terms tend to be irrelevant under RG iteration. In other words, the high-energy model, which is not symmetric but more realistic and less fine-tuned, has a tendency to flow to the symmetric point. As a result, higherrank symmetry as well as conservation of angular moments emerges.
The remainder of this paper is organized as follows. In Sec. II, we introduce the d-component bosonic systems and its higher-rank symmetry. In Sec. III, we discuss the scaling and Feynamn rules of the d-component bosonic systems. Further, we figure out the β-functions of parameters in the systems with renormalization group (RG) analysis and depict the global phase diagram. In Sec. IV, we summarize and provide our prediction on conditions of possible realization of the systems with higher-rank symmetry.

II. MODEL AND SYMMETRY
The symmetric point Hamiltonian for the d-component bosonic systems in real space [69] is given by H = H 0 +H 1 with where H 0 is free Hamiltonian and H 1 is interaction Hamiltonian. Here t a = m −1 a stands for the inverse of mass along the a-th direction and µ stands for the chemical potential.Φ † a (x) andΦ a (x) are respectively creation and annihilation operators of d-component bosons, and satisfy the bosonic commutation relations. The interaction strength K ab is a symmetric matrix with vanishing diagonal elements, i.e., K aa = 0, K ab = K ba . Each term in H is invariant under both the conventional global symmetry transformations (Φ a → e iθaΦ a , θ a ∈ R) and higher-rank symmetry transformations: for each pair (φ a , φ b ) with λ ab = −λ ba ∈ R. According to the Noether's theorem, the conventional global U (1) symmetry and the higher-rank symmetry correspond to conserved total charge (particle number)Q a and conserved total angular momentsQ ab (Q ab = −Q ba ) [1,5,69]: Hereρ a =Φ † aΦa . Intuitively, the conserved quantities Q ab enforce that a single a-th component boson can only move along the a-th direction. More explanation on the classical mechanical consequence of the conservation is available in Appendix A and Ref. [69].
In the coherent-state path-integral formulation with imaginary time, the Lagrangian density L can be written as L = φ * a ∂ τ φ a + H with action S =´dτ d d xL. Here the bosonic fields φ a = φ a (x, τ ) ∈ C are the eigenvales of coherent-state operatorsΦ a (x, τ ). With the Fourier transformation of the coherent state: φ a (τ, r) = . The right one represents the Feynman loop diagram of the correction for the parameter T1.
In the first term on the r.h.s., the simplified notation φ a stands for φ a (iω n , k) which is the frequencymomentum image of the field φ a (τ, r). ω n is a bosonic Matsubara frequency and k is a momentum vector: k = (k 1 , k 2 , · · · , k d ). We use k a to denote the a-th spatial component of momentum vector k. k stands for ωn´d d k (2π) d . In the second term, since there are four pairs of frequency and momentum, we introduce a new notation φ ia to compactly represent φ a (iω ni , k i ) where the label i = 1, 2, · · · , 4. k a i stands for the a-th spatial component of momentum vector k i . The sum 1,2,3,4 denotes Other notations like 1,2 , 2 , . . . in the forthcoming texts are defined in the similar way. Besides, the kinetic energy with momentum k is defined as: ξ a = 1 2 t a (k a ) 2 − µ(In the following text, we began to introduce anisotropic kinetic energy). For momentum k i , the associated kinetic energy Before moving forward, we perturb the Lagrangian density by adding small "off-diagonal" kinetic terms that break higher-rank symmetry. They are the terms deviating the model from the symmetric point. As such, kinetic terms of both directions are present, which can be written as The kinetic parameter t a is rewritten as t aa for the notational convenience. Those off-diagonal kinetic terms with nonzero t ab (a = b) manifestly break higher-rank symmetry. Similarly, we can also understand these parameters as inverse of mass of field configuration φ a along other directions other than the a-th one: t ab = 1/m ab .In this way, the kinetic energy with momentum k can be redefined as:

A. Scaling and Feynman rules
We consider d = 2. We set the restriction of field configuration φ a as: Here κ is an arbitrary constant with the momentum dimension and Λ is the cut-off of momentum here. The high-energy part corresponds to: where the scaling parameter s > 1 and sends k to k/s. We also define s = e l with l > 0. We consider the free part of the action In the following text, we also call it fast modes and denote it with >. It is noticed that we do not take the part related to chemical potential into account since it must be relevant if we choose the kinetic part to be marginal. Suppose the scaling dimension of φ is ∆ φ , and we can assume the change of temperature and frequencies are described as: T → s −z T , and ω n → s −z ω n , n ∈ Z respectively. So the free part is scaled to To fix the momentum dependence, we need ∆ φ = d 2 + 1 = 2 and z = 2. In d = 2 case considered here, the interaction matrix is simply determined by a single parameter, i.e., K 12 = K 21 := K. For further simplifying calculation for d = 2, we assume t 1 := t 11 = t 2 := t 22 := T 0 , t 12 = t 21 := T 1 .
There are two small parameters compared to T 0 . The first one is symmetry-disallowed off-diagonal kinetic parameter T 1 which is marginal. It can also been seen as a perturbation relative to the interaction parameter KΛ 2 . The second one is the irrelevant interaction parameter KΛ 2 considered as an infinitesimal quantity compared with diagonal kinetic parameter T 0 . Hence, the following calculation will be proceeded in the perturbative regime: T 0 KΛ 2 T 1 , where Λ is the momentum cut-off.
Next, we write Feynman rules. The bare Feynman propagators are given by − µ stands for the kinetic energy of φ ia [see definition below Eq. (7)]. In addition, the interaction vertex is drawn in Fig. 1. This diagram represents the interaction term, i.e. the second term in Eq. (7), where solid lines represent the bosonic fields and vertex stands for the interaction coefficient. In the RG procedure to be performed, we apply the standard cumulant expansion and relate the mean of the exponential to the exponential of the means [76,77]: Here φ > and φ < respectively correspond to the fast modes and slow modes of bosonic fields and S int stands for interaction part of action. In this case, the notation > is to take the average over the fast modes. Therefore, after calculating the average on fast modes, we have the form of the effective action.
The form of S int > contributes to the β-functions of kinetic parameters, also named first-order correction. Further, The interaction part of action reads: Below we determine the scaling of the K matrix. To be specific, we define the new momentum and frequency as: Then the scaling of parameters K ab can be determined. With the form of interacting action (8), we obtain: Here the dimensionful δ-function of momenta is also scaled [δ(k ) = s −d δ(k)] but the Kronecker symbol of Matsubara frequencies is dimensionless and invariant upon scaling. Therefore, the scaling of K ab is given by K ab = K ab s −d . It illustrates that the K matrix at tree level is irrelevant in perturbative RG. Then, let us consider the first-order correction to T 1 : S int > contributing to the correction of kinetic parameters.
To proceed further, we should split each momentum integration into slow part (<) and fast part (>). In this way, the integration of four momenta (k 1 , k 2 , k 3 , k 4 ) is split to 2 4 = 16 combinations. According to Wick's theorem, for formulas of bare propagators, we conclude that either φ * 1a and φ 4a must be paired or φ * 2b and φ 3b must be paired. Other contractions vanish due to K ab = 0 when a = b. Therefore, in 16 combinations, only the following two are non-vanishing: • Case-I: The momenta carried by φ * 1a and φ 4a , i.e., k 1 and k 4 , are fast momenta 1 . k 1 = k 4 and ω n1 = ω n4 required by the formula of bare propagator; • Case-II: The momenta carried by φ * 2b and φ 3b , i.e., k 2 and k 3 , are fast momenta. k 2 = k 3 and ω n2 = ω n3 required by the formula of bare propagator.
These two contractions correspond to the Feynman diagram in Fig. 1. We firstly focus on case-I. When 1,4 are fast modes, the momentum conservation and Feynman rules will tell us: k 2 = k 3 , k 1 = k 4 ; ω n2 = ω n3 , ω n1 = ω n4 . Since we focus on the two-dimensional case (d = 2), we set all off-diagonal K-matrix elements as: K 12 = K 21 := K. For further simplifying calculation, we assume t 1 := t 11 = t 2 := t 22 := T 0 , t 12 = t 21 := T 1 . In this way, S int I > is given by: where k 1 and k 2 correspond to fast and slow momenta respectively. The bosonic Matsubara summation is applied: 1 β ωn In the same way, we can give the form of case-II: By observing these two actions, we find they are completely equivalent to each other by exchanging the indices (a ↔ b , 1 ↔ 2). Then we focus on the case-I and multiply it by 2, namely: where we arrange all terms into three parts: S int S int The fast momentum k 1 = (k 1 1 , k 2 1 ) takes values in the elliptic shell (Λ is the momentum cut-off and the RG flow parameter s = e l with s → 1 and l → 0): which defines the domain of integral´>. Here we should focus on the integrated domain of fast momenta. It is only determined by the fields carrying those fast momenta. To be specific, we consider the fast momentum k i carried by the fast mode φ ia or φ * ia with the label i = 1, 2, ..., 4 defined above. Its integrated domain is given by This conclusion will be used in the following text. To compute S int Using the new momentum variables, we can work out the integral´> in S int (1) > in Eq. (13) (s−1 is small enough): where we have considered infinitesimal s − 1 and the definition of s = e l . And, ξ 0Λ := T0 2 Λ 2 − µ. As a result, where S 2 := 2π T0 T1 and f B (ξ oΛ ) := 1 e βξ 0Λ −1 . Below, we will use dl to replace l since it is infinitesimal.
The second one S int > in Eq. (14) vanishes since it is an odd function of the integrated momentum k 1 . The third term in Eq. (15) is directly connected to correction of the chemical potential. We just briefly give our results here for it is not our focus: Therefore, the only term that can renormalize T 1 is S int (1) > in Eq. (18). The bare off-diagonal kinetic term is given by a =b k T1 2 (k a 2 ) 2 φ * 2b φ 2b in the frequency-momentum space. Hence, we have the βfunction of the parameter T 1 by referencing the effective action: This corresponds to Fig. 1. Besides, by referencing the term containing the chemical potential, we have the β-function of the chemical potential µ: dµ dl = 2µ − 1 2 Λ 4 S2 (2π) 2 f B (ξ aΛ )K with the term 2µ originating from the contribution from the tree-level diagrams.

C. Vertex correction to K
We then turn to the vertex correction with all the contractions contained in S 2 int > : Since we only consider the loop-diagram contribution, only two kinds of contractions will not vanish under calculation. The first term S 2 int (1) > can be concluded as four cases: • Case 1: The momenta k 1 , k 2 , k 3 and k 4 are fast momenta. k 1 = k 3 , ω n1 = ω n 3 ; k 2 = k 4 and ω n2 = ω n 4 . • Case 2: The momenta k 1 , k 2 , k 3 and k 4 are fast momenta. k 1 = k 4 , ω n1 = ω n 4 ; k 2 = k 3 and ω n2 = ω n 3 .
These four cases correspond to the Feynman diagram in Fig. 2. Similarly, it can also be proved that these four cases are equal to one another by exchanging the indices. We directly take the Case 1 as an example and multiply it by 4 for simplicity. By further considering the delta functions in vertex correction, the relationships between momenta and frequencies are given by: k 1 = k 3 , k 2 = k 4 , k 1 + k 2 = k 3 + k 4 = k 1 + k 2 = k 3 + k 4 ; ω n1 = ω n 3 , ω n2 = ω n 4 , ω n1 + ω n2 = ω n3 + ω n4 = ω n 1 + ω n 2 = ω n 3 + ω n 4 . Hence, with the formula (20), we obtain that where we introduce the notation: This expression contains a sum of Matsubara frequencies n 2 : for Matsubara frequencies for bosons: ω n = 2πn β , we can obtain the result with substitutions: ξ 1 = iω n 1 + iω n 2 − ξ 1 +2 −2,a and ξ 2 = ξ 2b . The sum of n 2 can be rewritten as: where we apply f B (iω n1 + iω n2 − ξ) = > can be rewritten as: where the momenta k 1 and k 2 can be ignored in the expression since φ 1 a and φ 2 b are both slow modes. The integrated region of the momentum k 2 in S 2

int
(1) > is given by T0 k a 2 , we limit our integrated region on (Λ/s) 2 ≤k 2 ≤ Λ 2 . Therefore, the contribution from the first type of contraction is given by: where the parameter A takes the form of: where we use analytic continuation here and replace the imaginary frequencies with ω 1 and ω 2 . It is necessarily assume that the frequencies and momenta on the external lines can be ignored compared with ξ 0Λ and Λ. It is obvious that the expression (24) is positive. In other words, A > 0. More calculations on A are shown in Appendix B. By comparing the form of (23) with the bare interaction vertex, we find that it will correct the interaction parameter K effectively. According to the effective action:S ef f = S 0 + S int − 1 2 ( S 2 int − S int 2 ), we have the form of the β-function of parameter K: where the term −2K originates from the rescaling of slow momenta carried by the slow modes in (9) above. We have the general form of K(l) and T 1 (l): with the separatrix: The second type of contraction S 2
The Feynman diagrams related are presented in Fig. 3. After integrating the fast momenta, all of these terms will be turned into terms with constant coefficients. Therefore, integration from the two connected diagrams will be explained as corrections of other possible terms such as g|φ a | 4 . We do not care about this parameter for it does no contribution to the parameters K and T 1 . For simplicity, we do not present more details here.

D. Global phase diagram
In summary, we can see the first-order corrections only correct the kinetic terms on directions other than the a-th one of field configurations φ a [see definition below Eq. (7)]. In this way, the kinetic energy along the a-th axis will not be influenced by the contraction. We can directly prove that the higher-order correction will still contribute nothing to the parameter T 0 . Hence, we have the β-function for parameter T 1 corresponding to Fig. 1, given by Eq. (19). We can safely come to the conclusion that if K > 0, the parameter t 1 will reduce in the RG flow. This requires the parameter K to be positive. Only with positive K can parameter T 1 flow to zero in the RG analysis. Nevertheless, at this step, it is insufficient to tell whether the parameter T 1 is irrelevant here since the elements in the matrix K will also flow to zero. What we need is further calculation on vertex correction, which has been given in Eq. (25).
A numerical result of beta functions (19) and (25) is shown in Fig. 4, where the red line is the separatrix (28). In the region below the line, all T 1 (0) will flow to zero in RG flow, indicating the emergence of higherrank symmetry in low-energy physics since those terms violating the higher-rank symmetry will vanish in the low-energy physics.
In addition, after comparing the parameters T 1 (0) and K(0)Λ 2 near the separatrix, we have T 1 (0) K(0)Λ 2 T 0 with the approximation ξ0Λ 2k B T 1, which can be realized in cold-atom systems. This tells us our calculation is consistent with our initial assumption above.
From the phase diagram, we can safely come to the conclusion that the higher-rank symmetry will emerge after integrating the high-energy modes since off-diagonal kinetic parameters are just perturbations to the system. However, it should be noted that our RG analysis is only valid near the initial point l = 0. We here also do not take higher-order terms into account. Therefore, our calculation and results may be invalidated when l → ∞. The reason we do not just finish our RG analysis on the axis: T 1 (0) = 0 is we only consider the lower-order corrections. If we start with an initial point on the horizontal axis: T 1 (0) = 0, we have S 2 → ∞ in dT1 dl , which requires more higher-order terms since they contain expressions with higher order of S 2 . We have to take all of them into account for accuracy, which is beyond the current perturbation calculations. Although our calculation cannot take all the higher-order corrections into account, it also exhibits the tendency of RG flow giving the accurate prediction near the region with l = 0. Alternatively, the exact flat "band" along one direction for a boson, in analogy to bosonic Landau level, leads to strong correlation effect even though interaction is not large. We also comment that the perturbation calculations fail completely if we instead study models [68] with conserved dipole moments, where bosons are fractons that are immobile towards all directions. In such models, the Hamiltonian is intrinsically non-Gaussian with no kinetic terms at all. But analytic difficulty becomes much smaller when fracton condensation is considered as what Ref. [68] did.
Above discussion is all about the β-functions of parameters of 2-component bosons in two-dimensional space. Similarly, we can extend the case to the d-dimensional space with d > 2. We here give the forms of βfunctions of d-component bosonic systems for reference. The β-functions of t ab (a = b) and K ab (a = b) is given by ab − dK ab , where S d represents the volume of d-dimensional ellipsoid: They give the same type of RG flow phase diagram as Fig. 4.

IV. SUMMARY AND OUTLOOK
In this paper, we study how higher-rank symmetry (5) and angular moment conservation (6) emerge at low energies through a RG analysis. In other words, we identify them as emergent phenomena rather than strict properties of microscopic models. A phase diagram is given in Fig. 4, in which a wide parameter region is found to support emergent phenomena. Despite the limitation of oneloop perturbation techniques, we argue that emergence occurs in the deep infrared regime. On the other hand, by regarding higher-rank symmetry as emergent symmetry, our work opens a door to new way of thinking on realization of such unconventional symmetry and higher moment conservation in more realistic models, e.g., in simple frustrated spin models near symmetric points. Recently, some higher-moment conserving 1D spin systems have been found to support anomalously slow, subdiffusive late-time transport [74,78]. Thus, it is interesting to ask how emergent conservation of angular moments affect the late-time transport. Besides, the Hamiltonian can be reformulated on a square lattice. In an optical lattice of cold-atom experiments, one may choose lattice constant l = 700nm, Λ ∼ 1 l = 1.43 × 10 7 m −1 , m 0 = 1 t0 = 8.22 × 10 −34 kg, T = 10K. The condition t0Λ 2 −2µ 2k B T ∼ 10 3 where the chemical potential is negative indicates the possibility for simulating the 2-component bosons on the optical lattice. For accuracy and extension, it will be interesting to further consider the correction from higherorder term with more loops in the Feynman diagram. Finally, it will be interesting to study symmetry-protected topological phases (SPTs) with such emergent higherrank symmetry.
In the main text, we introduce the following conserved quantities (a, b = 1, 2, · · · , d ; a < b) in d-dimensional space: For each a, the conservation ofQ a requires that all ath component bosons are always in the d-dimensional space. The conserved quantitiesQ ab leads to the mobility restriction: a single a-th component boson is only allowed to move along a-th directions. To be more specific, we consider d = 2 and the following three quantities are conserved: Q 1 andQ 2 are the usual conserved charge (particle number) of bosons of the 1st and 2nd components, respectively. The conserved quantityQ 12 is the total angular moments formed by bosons. SupposeQ 1 = N 1 and Q 2 = N 2 . We can use Dirac function to expressρ 1 and ρ 2 in the eigenbasis: It means that bosons of the 1st (2nd) component are located in x i (x j ) with i = 1, 2, 3, · · · , N 1 (j = 1, 2, 3, · · · , N 2 ). Then, Q 12 are reduced to: From this expression, one can conclude that, if we move a 1st component boson, in order to keep Q 12 invariant, the boson is only allowed to be movable in the 1st direction such that its coordinate x 2 is unchanged. Of course, we can collectively move bosons of both components, such that the change in N1 i x 2 i can be canceled out by the change in N2 j x 1 j such that Q 12 is unaltered. This scenario is beyond the single particle movement and gives nontrivial effects when inter-component interaction (Kterm) is involved.