Universality classes of topological phase transitions with higher-order band crossing

In topological insulators and topological superconductors, the discrete jump of the topological invariant upon tuning a certain system parameter defines a topological phase transition. A unified framework is employed to address the quantum criticality of the topological phase transitions in one to three spatial dimensions, which simultaneously incorporates the symmetry classification, order of band crossing, $m$-fold rotational symmetry, correlation functions, critical exponents, scaling laws, and renormalization group approach. We first classify higher-order Dirac models according to the time-reversal, particle-hole, and chiral symmetries, and determine the even-oddness of the order of band crossing in each symmetry class. The even-oddness further constrains the rotational symmetry $m$ permitted in a symmetry class. Expressing the topological invariant in terms of a momentum space integration over a curvature function, the order of band crossing determines the critical exponent of the curvature function, as well as that of the Wannier state correlation function introduced through the Fourier transform of the curvature function. The conservation of topological invariant further yields a scaling law between critical exponents. In addition, a renormalization group approach based on deforming the curvature function is demonstrated for all dimensions and symmetry classes. Through clarification of how the critical quantities, including the jump of the topological invariant and critical exponents, depend on the nonspatial and the rotational symmetry, our work introduces the notion of universality class into the description of topological phase transitions.


Introduction
A recently emerged issue in the research of topological insulators (TIs) and topological superconductors (TSCs) is the understanding of quantum criticality near topological phase transitions [1,2,3,4,5,6]. The topology of TIs and TSCs are characterized by the topological invariant defined from the Bloch wave function, which is calculated by different means according to the dimension and symmetry class of the system [7,8,9,10]. The discrete jump of the topological invariant upon tuning a certain parameter M at a critical point M c defines the topological phase transition. In practice, the lowenergy effective theory of these materials are well described by a Dirac Hamiltonian that respects the symmetry of the system, and the tuning parameter M represents the mass term, or equivalently the bulk gap of the energy spectrum, of the Dirac Hamiltonian. At the critical point M c , the bulk gap closes, i.e., a topological phase transition necessarily involves a gap-closing. However, one should keep in mind that the converse is not true, namely gap-closing in the energy spectrum does not necessarily involve topology.
The first step towards understanding quantum criticality of such transitions is, of course, whether any asymptotic critical behavior manifests as M → M c . Given that the system does not necessarily possess a Landau order parameter, and the topological invariant jumps discretely at the critical point, the notion of asymptotic critical behavior that unambiguously defines the transition does not seem to be obvious at a glance. Recently, a series of work on scaling and critical exponents shed a light on this issue [11,12,6]. The observation in these works is that the topological invariant, being a global property of the Bloch state, is always calculated through momentum space integration over a certain function referred to as the curvature function. If M is varied but the system remains the same topological phase, then the topological invariant remains unchanged, but the profile of the curvature function changes. Moreover, as M → M c , the curvature function gradually develops a divergence at the high symmetry point (HSP) where the gap-closing takes place. On the two sides of the critical point M + c and M − c , the divergence is of opposite sign. This divergence and sign change of the curvature function serve as the asymptotic critical behavior that unambiguously identifies a topological phase transition, in constrast to any other quantum phase transition that does not involve topology.
Another major step towards understanding the topological phase transition is the identification of correlation functions, despite the system may not possess a local order parameter. Based on the theory of charge polarization [13,14,15] and the theory of orbital magnetization [16,17,18,19,20], it is recognized that in certain symmetry classes, the Fourier transform of the curvature function yields a correlation function. In noninteracting systems, the correlation function is a measure of the overlap of Wannier functions that are a certain distance apart [6]. Because the curvature function diverges as M → M c , so does the correlation length ξ, consistent with the notion of scale invariance in the usual second-order phase transition. In addition, the power law divergence of ξ with respect to M defines the critical exponent ν. It is further unveiled that the exponent ν is determined by the order of band crossing n of the Dirac Hamiltonian at the critical point [6]. Depending on whether the bands at the critical point M = M c show linear, quadratic, or higher-order band crossing, the critical exponent and the jump of topological invariant take different values. These statements are drawn from investigating systems in a specific symmetry class, namely 2D class A in the periodic table of the topological materials, and is also found to manifest in 2D strongly interacting TIs solved by means of twisted boundary conditions [21], and weakly interacting 1D and 2D TIs solved by means of single-particle Green's functions [22].
On the other hand, previous investigation shows that an important mechanism to stabilize the aforementioned higher-order band crossing is the discrete m-fold rotational symmetry of the underlying lattice structure [23,24]. Through investigating 2D and 3D multiple Weyl points [25,26,27,28,29,30,31,32,33,34], the rotational eigenvalues are shown to be in one to one correspondence with the order of band crossing in the low-energy sector of these materials, which are described by higher-order Dirac or Weyl models [23,24]. The inclusion of additional nonspatial symmetry further constrains the rotational eigenvalues and hence the band crossing [24]. However, shall this stabilization mechanism be generally applicable to TIs and TSCs in any symmetry class, it must also take into account the fact that the order of band crossing is already constrained by the nonspatial symmetries in the symmetry classification [7,8,10].
In this article, we investigate the quantum criticality of topological phase transitions within a framework that simultaneously incorporates all the aforementioned ingredients, namely (1) symmetry classification, (2) order of band crossing, (3) m-fold rotational symmetry (for 2D and 3D), (4) correlation function, (4) critical exponents, and (5) the discrete jump of topological invariant. Our classification scheme is based on the following procedure. We first classify the higher-order Dirac models from 1D to 3D according to the time-reversal (TR), particle hole (PH), and chiral symmetry, and clarify the evenoddness of the order of the band crossing n in each symmetry class. We reveal that the order of band crossing n in certain classes also determines the jump of topological invariant ∆C per gap-closing. When rotational symmetry is further imposed in 2D and 3D, the even-oddness of n also constrains the rotational eigenvalues. As a result, the allowed rotational symmetry m depends on the dimension, symmetry class, and the system being fermionic or bosonic. The notion of Wannier state correlation function is clarified for all the 15 topologically nonotrivial dimension-symmetry classes. Moreover, we verify that the critical exponent ν of the correlation length ξ, as well as that of the edge state decay length in the topologically nontrivial phase, is always determined by the order of band crossing n, but not necessarily m or ∆C. On top of these, we will demonstrate that a previously proposed scaling law [6], as well as a renormalization group approach based on the curvature function [11,12], remain valid for all the 15 cases. All these features together give rise to a coherent picture for the quantum criticality of topological phase transitions, especially how the critical quantities {∆C, n, ν} depend on the dimension and symmetry, which we refer to as the universality classes of TIs and TSCs.
The article is organized in the following manner. In Sec. 2, we give an overview of the generic critical behavior that is universal to all the 15 dimension-symmetry classes from 1D to 3D, including the divergence of the curvature function, critical exponents, scaling laws, and the renormalization group approach. In Secs. 3 to 5, we detail the interplay between TR, PH, chiral, as well as the m-fold rotational symmetry in 2D and 3D, especially how they influence the critical quantities {∆C, n, ν}. The results are concisely summarized in Table 1. Because some classes in lower dimensions can be obtained from higher ones through a dimensional reduction [8], we start from 3D systems, and then reduce the dimension to 2D and 1D. Section 6 summarizes the features obtained within this framework. Table 1. Summary of the results in the 15 topologically nontrivial dimensionsymmetry classes from 1D to 3D, which lists the topological invariant C, jump of topological invariant across a topological phase transition per gap-closing ∆C, the order of band crossing n, the allowed m-fold rotational symmetry for fermionic F = 1 and bosonic F = 0 systems in 2D and 3D, and the corresponding Wannier state correlation functionF . 2. Generic critical behavior of higher-order Dirac models

Curvature function and correlation function
We discuss the TIs and TSCs whose low-energy effective theory is a Dirac model satisfying the following criteria: (1) The low-energy effective Dirac Hamiltonian H(k) = d(k) · Γ is constructed out of Γ a matrices that satisfies the Clifford algebra Γ a , Γ b = 2δ ab . In addition, the low-energy dispersion takes the form where the integer n signifies the order of band crossing at the topological phase transition M c = 0. The mass term M is assumed to have no momentum dependence.
(2) The model is classfied according to the TR, PH and chiral symmetries In each spatial dimension from 1D to 3D, there are 5 symmetry classes that are topologically nontrivial. We will discuss all the 15 nontrivial dimension-symmetry classes case by case in the following sections.
(3) In 2D and 3D, we consider m-fold rotation symmetric Dirac Hamiltonians and choose an eigenvalue basis in which both the Dirac Hamiltonian and the rotation operator are diagonal. In this basis, the operator of the m-fold rotation symmetry can be written as where p = 0, 1, 2...m − 1 are the rotational eigenvalues. The spin statistical factor F indicates whether the basis is fermionic F = 1 or bosonic F = 0. We assume that in each symmetry class, the basis can be either fermionic or bosonic. The rotational symmetry requires that where R m rotates the momentum k about the k z axis by 2π/m. In 1D this rotational symmetry is absent.
As will be demonstrated case by case in the following sections, the topological invariant C (or a related one that judges the topological phase transition) in all the 15 dimension-symmetry classes can be cast into the form of a D-dimensional integration over the Brillouin zone (BZ) where F (k, M ) is referred to as the curvature function, the precise form of which depends on the dimension and symmetry class. The points in the BZ satisfying k 0 = −k 0 (up to a reciprocal lattice vector) are referred to as the high symmetry points (HSPs). The curvature function is generally an even function F (k 0 + δk, M ) = F (k 0 − δk, M ) around the HSP due to certain symmetry, such as inversion or TR symmetry. While the topological invariant C remains constant within a topological phase in the parameter space of M , the profile of F (k, M ) varies with changing M , which is key to our analysis of critical behavior.
Another key ingredient in our analysis is the Wannier states |Rn defined from the Bloch state by where N denotes the number of lattice sites, andr is the position operator. We propose a correlation function derived through the Fourier transform of the curvature function, In each of the 15 dimension-symmetry classes, we will demonstrate that the correlation function always takes the form of measuring the overlap of Wannier functions.

Generic critical behavior and universality classes
Through investigating the low-energy continuous models of all the 15 nontrivial dimension-symmetry classes from 1D to 3D, we found that the critical behavior of the curvature functions fall into two different scenarios: Peak-divergence scenario.-For linear band crossing n = 1, the curvature function F (k, M ) peaks at the gap-closing HSP. The peak gradually narrows and the height increases as the system approaches the critical point M → M c , and eventually the peak diverges and flips sign across the transition. We call this critical behavior the peakdivergence scenario. In this scenario, the peak can be well fitted by an Ornstein-Zernike form The critical behavior described above is summarized as Denoting the critical exponents of F (k 0 , M ) and ξ as the conservation of the topological invariant C div under the diverging peak yields a scaling law which is to be satisfied by any linear, isotropic Dirac model in any dimension. Because of this Ornstein-Zernike form, the Wannier state correlation function as a Fourier transform of the curvature function, as in Eq. (7), decays with correlation length ξ. Shell-divergence scenario.-For any high order band crossing n > 1, the extremum of the curvature function F max (k max , M ) forms a D − 1 dimensional shell surrounding the HSP, with a radius denoted by k max and a thickness denoted by k wid . The critical behavior is that the shell gradually reduces it radius as M → M c , and F max (k max , M ) gradually diverges and flips sigh across the transition. In other words, For the Dirac models in Sec. 2.1, the k max and k wid have the same critical exponent The conservation of the portion of the topological invariant C div in the diverging shell where Ω takes into account the angular integration, yields the same scaling law as Eq. (12). The two scenarios are depicted pictorially for 1D and 2D in Fig. 1. To be more specific, find that the low-energy continuous model for most of the symmetry classes has the topological invariant in Eq. (5) in the following scaling form Expanding the curvature function in the integrand for small k yields For the peak-divergence scenario at n = 1, the Ornstein-Zernike form is evident, which renders γ = D and ν = 1, satisfying Eq. (12). For the shell-divergence scenario at n > 1, solving for the k max from ∂ k F (k, M ) = 0 yields k max ∝ |M | 1/n , and subsequently we obtain F max (k max , M ) ∝ |M | −D/n . Consequently, the critical exponents are γ = D/n and ν = 1/n, in agreement with the scaling law in Eq. (12). This analysis indicates that the critical exponent ν for the length scale ξ are basically determined by the order of band crossing n, as can also be inferred by counting the dimension in the dispersion in Eq. (1). The edge state decay length in the topologically nontrivial phase has the same critical exponent as the correlation length. This can be seen by considering the model to be defined in the half-space x > 0, and considering the edge state with zero transverse momentum, i.e., k y = k z = 0 in 3D and k y = 0 in 2D. The problem can be calculated by projecting the higher-order Dirac Hamiltonian into real space and solve for the zeroenergy edge state [6] where Γ x and Γ M represent the Γ-matrices that multiply k n x and M , respectively, whose precise form depends on the symmetry class and dimension. Multiplying the equation by Γ x and using the ansatz The eigenvalue η is chosen such that one of the roots ξ loc is real and identifiable with the edge state decay length, which obviously also has critical exponent ν = 1/n. Thus, the correlation length and the edge state decay length are essentially the same length scale, although we should emphasize that the edge state only exists in the topologically nontrivial phase, whereas the Wannier state correlation function is well-defined in either the trivial or nontrivial phase.

Curvature renormalization group approach
The curvature renormalization group (CRG) approach has been proposed to judge topological phase transitions in a variety of systems [11,12,6,21,22]. Here we demonstrate that the method applies to all the 15 dimension-symmetry classes at any order of band crossing n and rotational symmetry m. The approach is essentially an iterative method to search for the trajectory (RG flow) in the parameter space of M along which the maximum of the curvature function is reduced. In this way, the topological phase transitions, at which the curvature function diverges, can be identified. The scaling equation demands that, at a given M , we find the new M that satisfies The mapping M → M yields the RG flow that identifies the topological phase transitions. Here δk is a small momentum displacement away from the HSP, and one chooses 0 < b < 1 for the peak-divergence scenario n = 1, whereas b > 1 for the shell-divergence scenario n > 1. The divergence of the curvature function is gradually reduced under the scaling procedure, as shown schematically in Fig. 2. We remark that a similar RG approach has also been proposed to search for the transition based on the reduced density matrix obtained from bipartition, since it shows a similar critical behavior [35].
Since the curvature function in many symmetry classes takes the generic scaling form of Eq. (17), we examine here explicitly the CRG approach applied to this scaling form. As the approach concerns only a small k region near the HSP which is set to be at the origin in our continuous models k 0 = 0, we expand the scaling form in Eq. (17) by For the peak-divergence scenario n = 1, putting Eq. (21) into the leading order expansion of the scaling equation, Eq. (20), yields the RG equation The corresponding RG flow flows away from the critical point M c = 0, with a flow rate that diverges as M → M c . For the shell-divergence scenario at n > 1, we approximate Eq. (21) by Treating δM M and ∆k = (b − 1)δk δk as small quantities and expanding to the leading order, and then putting Eq. (23) into the Eq. (20) yields the leading order RG equation The equation corresponds to an RG flow that flows away from M c = 0, with a flow rate that vanishes at M = M c . In other words, the topological phase transition manifests as an unstable fixed point at n > 1.

Topological phase transitions in three dimensions
In this section, we study topological phase transitions in three dimensions. We consider all the five symmetry classes that have nontrivial topology. The Wannier state correlation functions of class AIII, DIII, CII, and CI are constructed in a similar manner from the winding numbers, whereas in class AII it is constructed differently from the

3D class AIII
For 3D class AIII that has S 2 = 1, we choose the Γ-matrices [8] Γ a = α x , α y , α z , β, −iβγ 5 , The chiral operator is S = β, which together with the chiral symmetry in Eq. (2) requires that d 4 (k) = 0, whereas there is no constraint on d i (k) being even or odd in k for i = 1, 2, 3, 5. We follow the convention to choose d 5 = M as the mass term. The Hamiltonian then takes the form where we have denoted with n = n + + n − , such that the dispersion satisfies with n ∈ Z.
Using the rotational symmetry operator the chiral symmetry [C m , S] = 0 does not give any constraint on C m . The rotational symmetry in Eq. (4) demands Because g = g − iM contains the mass term that must be invariant under the rotational transformation, it is required that α p α * r = 1 and α s α * q = 1, and therefore r = p and q = s. The constraint on g then gives g(k ± , k z ) = g(k ± e ±i 2π m , k z ), which is consistent with our parametrization in Eq. (27). The parametrization of f in Eq. (27), together which is the relation between the parametrization and the rotational symmetry eigenvalues, and we have used n + − n − ∈ Z due to n = n + + n − ∈ Z. The x in Eq. (31) and throughout the article is an arbitrary integer. Evidently, Eq. (31) can be satisfied by any of the m = 2, 3, 4, 6 for either the bosonic or fermionic basis. We define the q-matrix by (not to be confused with one of the rotational eigenvalues) Writting the parametrization in Eq. (27) into cylindrical coordinate (k ⊥ , φ, k z ) by f = k n ⊥ e i(n + −n − )φ and g = k n z , and using the derivatives yields the following combination in the integration of topological invariant [8] µνρ Converting the integrand into spherical coordinates (k, θ, φ) yields The topological invariant is given by the 3D winding number obtained from the integration of the above quantity where We see that there is a topological phase transition ∆C = 0 if and only if the order of band crossing n is an odd number. Here the extra factor 2nB( n 2 , n) is an artifact of this continuous model, and the jump of topological invariant should be treated as ∆C = |C(M > 0) − C(M < 0)| = n + − n − ∈ 2Z + 1, as demonstrated by the following lattice model The resulting topological invariant is which gives ∆C = n − ∈ 2Z + 1.
The Wannier state correlation function in this class is introduced from the Chern-Simons gauge field A ab µ = u a− |∂ µ |u b− defined from the eigenstates for the filled bands We find that the integrand in Eq. (36) can be written as As a result, the topological invariant in Eq. (38) is equivalently In addition, the Chern-Simons gauge field A ab µ can be expressed as the Fourier transform of the filled-band Wannier states where |Rb is the Wannier state of filled band b localized at home cell R. We use following compact notation for the eigenstates and the corresponding Wannier states It then follows that the topological invariant in Eq. (41) can be written in terms of the Wannier states entirely in real space where we have converted the discrete sum of Bravais points into an integral. Similarly, we introduce the Wannier state correlation function from the Fourier transform of the integrand in the topological invariant These real space representations are presented diagrammatically in Fig. 3. In addition, Eq. (40) and Eq. (35) indicates that the curvature function satisfies the scaling form in Eq. (17) with D = 3, and hence the correlation length ξ has critical exponent ν = 1/n. Finally, the Chern-Simons invariant [8], equivalently the magnetoelectric polarizability, can as well be expressed in terms of Wannier states according to the same argument.

3D class DIII
For 3D class DIII that has T 2 = −1 and C 2 = 1, we follow Ref. [8] and use the same Γmatrices as in Eq. (25). The TR and PH operators are T = σ y ⊗ τ x K and C = σ y ⊗ τ y K. It follows that the TR symmetry in Eq. (2) requires Thus only the d 5 term can be the mass term d 5 = M . On the other hand, the PH symmetry requires Thus the d 4 term cannot exist, and the Hamiltonian takes the same form as Eq. (26), with the same parametrization on f and g as in Eq. (27). Given the rotational symmetry operator in Eq. (29), demanding [T, C m ] = 0 and [C, C m ] = 0 requires α p = α * s , α q = α * r . The rotational symmetry in Eq. (4) requires that from which we see that α p = α * q , and hence the rotational operator is constrained to be the of form The constraint on f in Eq. (49) renders which together with the requirement that f = d 1 − id 2 must be odd in momentum yields We see that the fermionic case F = 1 can be satisfied by any of the m = 2, 3, 4, 6, but bosonic case F = 0 can only be realized at m = 3. The topological invariant in this class is the 3D winding number discussed in Sec. 3.1, and hence the formalism of Wannier state correlation function, correlation length, and critical exponents follows that from Eq. (39) to Eq. (45).

3D class AII
For 3D class AII that has T 2 = −1, we use the same representation of the Γ-matrices as in Eq. (25), and the TR operator T = iσ y ⊗ IK following Ref. [8]. The TR symmetry requires that and hence only d 4 = M can be the mass term, and all others must be odd functions of momentum. Demanding the rotational operator in Eq. (29) to commute with the TR operator [C m , T ] = 0, renders α p = α * q and α r = α * s , and consequently C m = diag(α p , α * p , α r , α * r ). We may write the Hamiltonian into the form where f = d 1 − id 2 and g = d 3 − id 5 . The rotational symmetry in Eq. (4) demands that We parametrize f and g by With this parametrization in mind, the constraint on g in Eq. (55) demands that p = r, thus the rotational operator takes the form of Eq. (50). Consequently, we arrive at the same constraint to the rotational eigenvalues as Eq. (52).
The Z 2 invariant C for 3D lattice models in class AII can be constructed from the Pfaffian of the filled-band m-matrix m αβ (k) = u αk |T |u βk at HSPs [36,37], provided the Pfaffian is a real number [36,37] where k 0,i is the i-th HSP. The energy eigenstates for the filled bands E − (k) = −d are The Pfaffian that enters the Z 2 topological invariant in Eq. (57) is The critical behavior of the Pfaffian is that as M → M c , although the Pfaffian at k 0 remains at ±1, its second derivative diverges (or equivalently the Laplacian diverges since the model is isotropic around each k 0 ). This motivates us to consider an additional topological invariant C that uses the Laplacian of the Pfaffian as the curvature function This choice of curvature function is sound because it integrates to a topological invariant, which can be seen by considering the integration over the BZ of a D-dimensional cubic lattice Moreover, since Pf[m] is usually an even function, its Laplacian ∇ 2 k Pf[m] is also an even function. Note that since C is always zero, it is not directly related to the Z 2 invariant in Eq. (57).
We introduce the Wannier state correlation function for TR-invariant systems by considering the Fourier transform of the curvature function in Eq. (60). In terms of the Wannier state |Rn of the n-th filled-band localized at home cell R, introduced in Eq. (6), the derivatives of the Pfaffian take the form Universality classes of topological phase transitions with higher-order band crossing 16 The Fourier transform of the Laplacian of the Pfaffian is theñ where . This corresponds to a correlation function that measures the overlap of the Wannier state of the first filled-band centered at R and that of the second filled-band centered at the origin, as a matrix element of R 2 T .
The Laplacian of the Pfaffian in our continuous model in spherical coordinates is given by whose Fourier transform reads At any n, the correlation function is a decaying function of argument R/ξ, with a correlation length ξ ∝ |M | − 1 n that has critical exponent ν = 1/n. At n = 1, the correlation function monotonically decays, whereas at n > 1, the correlation function oscillates and decays, as shown in Fig. 4. These features remain true for all the correlation functions discussed in the present work.

3D class CII
For 3D class CII that has T 2 = −1 and C 2 = −1, the minimal model is a 8 × 8 Dirac model [7]. We choose the chiral representation, in which the seven Γ-matrices are given by [8] Γ a = Γ a 4×4 ⊗ η x , for a = 1 ∼ 4 , where Γ a 4×4 are those in Eq. (25). Γ 6 = I σ ⊗ I τ ⊗ η z = S is used to implement the chiral symmetry. The Hamiltonian written in the basis of the other six Γ-matrices is block-off-diagonal We will denote σ i , τ i , and η i as the Pauli matrices for the spin, particle-hole, and valley grading. The TR and PH operators may be chosen as T = iσ y K ⊗ I τ ⊗ I η and C = −iσ y K ⊗ I τ ⊗ η z . The block-off-diagonal part of the Hamiltonian in Eq. (67) expressed in terms of the d-vector yields The TR symmetry requires that Consequently, the mass term is d 4 = M . The PH symmetry also gives the same condition. We parametrize the 8 × 8 rotational operator by with α p i = e i 2π m (pi+ F 2 ) . The coexistence of TR and rotational symmetry [C m , T ] = 0 requires that α p 1 = α * p 2 ≡ α p . The coexistence of PH and rotational symmetry [C m , C] = 0 also requires the same thing. As a result, the α p in Eq. (70) takes the form α p = diag(α p , α * p ). We now examine the rotational symmetry. To simplify the calculation, we rewrite the Hamiltonian in Eq. (67) by Universality classes of topological phase transitions with higher-order band crossing 18 The rotational symmetry in eq. (4) gives the constraint The first two equations in Eq. (72) yields α p = α q = α r = α s , so the final form of the rotational operator is Using the parametrization such that the dispersion satisfies Eq. (28), the constraint on f yields because n = n + + n − ∈ 2Z + 1 due to Eq. (69). For fermions F = 1, the condition in Eq. (75) can be by any of the m = 2, 3, 4, 6, where as for bosons F = 0 it can only be satisfied at m = 3. We now introduce a winding number and the corresponding Wannier state correlation function. Defining the 4 × 4 q-matrix from the off-diagonal part of the Hamiltonian as [7] we found that µνρ Tr q † ∂ µ qq † ∂ ν qq † ∂ ρ q = 0, and hence the simple version of the winding number in Eq. (36) vanishes. However, if we add a Γ 5 4×4 (defined in Eq. (25)) into the definition then it yields a non-vanishing winding number C (see Eq. (36)) where B(a, b) is the Beta function. Although this winding number does not correspond to the true Z 2 invariant in this class, the integrand in Eq. (78) satisfies all the critical properties of the curvature function discussed in Sec. 2.2. In addition, from the eigenstates of the filled bands the Berry connection A ab

3D class CI
For 3D class CI that has T 2 = 1, C 2 = −1, we use the minimal 8 × 8 model that writes the Hamiltonian in the chiral basis [7] where the 4 × 4 matrices {α i , β, γ 5 } are those in Eq. (25). The D matrix expressed in terms of the d = (d 1 , d 2 , d 3 , d 4 ) vector is with f = d 1 − id 2 and g = d 3 + id 4 . The 8 × 8 TR and PH operators are T = I ⊗ I ⊗ σ x K and C = I ⊗ I ⊗ (−iσ y )K. The TR and PH symmetry require Consequently, only d 4 = M can be the mass term. Now we consider the 8 × 8 rotational operator in Eq. (70). The coexistence of TR and rotational symmetry, and the coexistence of PH and rotational symmetry imply α p = α * r and α q = α * s , so the rotational operator is constrained to take the form We use Eqs. (81) and (82) to write the Hamiltonian into the form Writing down each component in the 8 × 8 matrix explicitly and using the definition in Eq. (70), we see that rotational symmetry in Eq. (4) requires To satisfy the Dirac form of the dispersion, we proceed to parametrize with n = n + + n − ∈ 2Z + 1 due to Eq. (83). Since g is invariant under rotation g(k) = g(R m k), it requires that α p 2 α q 1 = α * p 1 α * q 2 = α p 1 α q 2 = 1. A detailed analysis shows that these conditions cannot be satisfied at any nonzero m for either fermions F = 1 or bosons F = 0, hence 3D class CI is not compatible with the m-fold rotational symmetry.
We proceed to discuss higher-order Dirac model n ≥ 1 without considering the rotational symmetry. Defining the q-matrix the winding number is calculated from the integration of which is twice of that in Eq. (34) for 3D class AIII. As a result, the winding number defined in Eq. (36) is always an even number C ∈ 2Z. The notion of correlation function and critical exponents then follows the discussion in Sec. 3.1.

Topological phase transitions in two dimensions
Among the five topologically nontrivial symmetry classes in two dimensions, as we address in this section, the topological invariant in class A, C, and D are calculated from the same Berry curvature formula, and so follows the same correlation function. On the other hand, the class AII and DIII are described by the Z 2 invariant and the correlation function derived accordingly.

2D class A
The minimal model for a 2D class A model is described by the Hamiltonian The rotational operation reads C m = diag(α p , α q ), and hence the rotational symmetry in Eq. (4) requires The parametrization yields since there is no nonspatial symmetry in this class to constrain the order of band crossing, one has n = n + + n − ∈ Z. This condition can be satisfied by any m for either bosons or fermions.
The topological invariant C in 2D class A is given by the integration of the Berry curvature Ω(k x , k y ). The valance band Berry connection as A = u k |i∇ k |u k , which yields the Berry curvature The derivatives may be converted into the polar coordinates by which yields the Berry curvature The topological invariant is calculated by converting the integration over the BZ into an integration over the entire polar plane The above invariant is normalized such that the change of topological invariant across the topological phase transition ∆C = C(M > 0) − C(M < 0) = n + − n − is consistent with ∆C = 1 for the usual linear Dirac model (n + , n − ) = (1, 0). The Fourier transform of the Berry curvature can be expressed in terms of the Wannier states by [38,39,40,6] which measures the overlap of the Wannier functions centering at R and at the origin 0, sandwiched by the factor R x r y − R y r x . Within our continuous model formalism where J 0 (Rk) is the 0 th -order Bessel function. Using the Berry curvature in Eq. (96), we find that the Fourier transform at any {n + , n − } gives a decaying function λ(R) with a correlation length proportional to ξ ∝ |M | −1/n ∝ |M | −ν . Moreover, because {n + , n − } are positive numbers, and hence n ≥ |n + − n − |, we have the inequality i.e., the order of band crossing can only be greater than or equal to the jump of topological invariant. Finally, we remark that out of all the correlation functions introduced in the present work, onlyF 2D (R) is gauge invariant (because Berry curvature is gauge invariant) and hence measurable. It can be measured by, for instance, performing a Fourier transform of the Berry curvature measured in the cold atom experiments [41,42,43], from which the critical behavior discussed above may be verified.

2D class C
We now discuss 2D class C that has only particle-hole symmetry C 2 = −1. The minimal model is a 2 × 2 Dirac model with the PH operator C = σ y K. Using the rotational operator C m = diag(α p , α q ), the requirement [C, C m ] = 0 yields α q = α * p , and hence the rotational operator has to take the form C m = diag(α p , α * p ). On the other hand, the PH symmetry requires d i (k) = d i (−k), i.e., all the three components are even in momentum. Without loss of generality we choose d 3 to be the mass term, and parametrize the Hamiltonian by The requirement that d i (k) is even in momentum dictates n ∈ 2Z, indicating that the order of band crossing can only be even. The rotational symmetry in Eq. (4) renders and consequently must be satisfied. As a result, when the system is fermionic F = 1, it can only have m = 3 fold rotational symmetry, whereas the bosonic case F = 0 can be satisfied with any m. The topological invariant is calculated from the same Chern number as in 2D class A in Sec. 4.1. Consequently, the Wannier state correlation function and the critical exponent of the correlation length are introduced exactly in the same way, using the formalism from Eq. (94) to (100). However, it should be emphasized that in 2D class C here the jump of topological invariant ∆C, order of band crossing, and inverse of the critical exponent are all even numbers, in contrast to 2D class A in which they can be either even or odd.

2D class D
The calculation of 2D class D, where the PH symmetry satisfies C 2 = 1, is similar to that in 2D class C presented in Sec. 4.2. We follow the formalism in Ref. [8] and use the PH operator C = σ x K. The PH symmetry requires We thereby choose d 3 = M as the mass term. Using the Hamiltonian and the lowest order expansion of f in Eq. (101), we see that the order of band crossing is odd n ∈ 2Z+1. The commutation of rotational operator C m = diag(α p , α q ) with the PH operator [C, C m ] = 0 yields α q = α * p . As a result, the rotational eigenvalues take the form C m = diag(α p , α * p ). Equation (4) renders the following condition on f Therefore, must be satisfied. If the system is fermionic F = 1, this condition can be satisfied in any m-fold rotational symmetry, whereas if the system is bosonic F = 0, then only m = 3 is allowed. The topological invariant is calculated from the Chern number in Sec. 4.1, and the Wannier state correlation function and the critical exponent of the correlation length are introduced by the formalism from Eq. (94) to (100). However, it should be emphasized that in 2D class D here the jump of topological invariant ∆C, order of band crossing, and inverse of the critical exponent are all odd numbers.

Now the representation in Eq. (107) indicates that the Hamiltonian takes the form
where we have defined f and g, and we assume their leading order expansions are Note that TR invariance requires so we use d 3 = M as the mass term. The oddness in Eq. (110) requires the expansion in Eq. (109) to satisfy n ∈ 2Z + 1 and + + − = 2Z + 1. Equation (4) yields the following constraints on f and g The last two equations of the above equation can never be satisfied at the same time, so we conclude that g = d 4 − id 5 = 0. Thus rotational symmetry requires that Γ 4 and Γ 5 terms cannot exist. The constaint on f then gives which can be satisfied by any m for either fermionic or bosonic cases. Denoting d = d 2 1 + d 2 2 + d 2 3 , the four eigenstates in a gauge choice that makes the Pfaffian real and equal to ±1 at HSPs are where the upper signs are for the occupied states {1, 2} which have eigenenergies Thus the sign of d 3 /d at HSPs determines the topology of the system according to Eq. (57).
Using the Pfaffian in Eq. (114), its Laplacian is The correlation function in our continuous model renders which displays the same behavior as other correlation functions.

2D class DIII
The discussion of 2D class DIII follows that in 3D class DIII with a dimensional reduction k z = 0, using the same Γ-matrices and TR and PH operators described in Sec. 3.2. The analysis from Eq. (47) to (50) remains true in 2D. We now translate the 3D Hamiltonian in Eq. (26) into 2D and parametrize f and g by First let us look at the parametrization of g . The constraint on g is g (k ± ) = g (k ± e ±i 2π m ).
m ( + − − ) k + + k − − , and hence + = − = /2, which implies g = k − iM only depends on the module of momentum, which contradicts the requirement that the d 3 (k) in g = d 3 − iM must be odd in momentum. Thus we conclude that the d 3 term cannot exist, and we have eventually g = −iM . The oddness of d 1 and d 2 again yields the same constraint between the n + − n − and the rotational eigenvalues as in Eq. (52), which can be satisifed by any m for fermions, and m = 3 for bosons.
Using the eigenstates in Eq. (39), with setting d 3 = 0 as discussed above, the Z 2 topological invariant can be constructed from the Pfaffian of the m-matrix using the formalism in Sec. 4.4. The TR operator in Sec. 3.2 acting on the filled-band eigenstates yields The Pfaffian of the m-matrix is then So we see that the sign of M determines the Z 2 index in our continuous model. Obviously, because of the Z 2 invariant, we always have ∆C = 1. The rest of the correlation function, critical exponents, and universality class follows that of Sec. 4.4 for 2D class AII.

Topological phase transitions in one dimension
In one dimension, the topological invariant in class BDI, AIII, and D are described by the same Zak phase. The class DIII is described by the Z 2 invariant calculated from the TR operator. On the other hand, the topology of class CII is given by a winding number of its own kind. We will detail all these calculations in this section.

1D class BDI
Following Ref. [8], we model the 1D class BDI systems, which have T 2 = 1 and C 2 = 1, by T = σ z K and C = σ x K. The TR symmetry requires On the other hand, the PH symmetry requires Combining the TR and PH symmetry, we see that the d 2 term cannot exist, and we may use d 3 = M as the mass term. The oddness of d 1 implies that we can parametrize the Hamiltonian by with an odd order of band crossing n ∈ 2Z + 1.
Denoting the valance band eigenstate of the Hamiltonian in Eq. (122) by |u − , the corresponding Berry connection actually vanishes A(k) = u − |i∂ k |u − . Thus for the sake of constructing the topological invariant and introducing the correlation function, we rotate the Hamiltonian by That is, we rotate the basis such that the mass term resides in the σ y component. The eigenstate in this basisH(k)|ũ − = E(k)|ũ − is related to the original basis by a rotation |ũ − = V k |u − , and gives a Berry connectioñ Using the parametrization in Eq. (122), we obtain the Berry connectioñ .
The topological invariant is the Zak phase calculated from the integration of Berry connection, which in our continuous model reads The topological invariant constructed from the Wannier state reads, according to the theory of charge polarization [13,14], The correlation function is given by the Fourier transform of the Berry connectioñ which measures the overlap of the Wannier state centering at the origin and that centering at R, sandwiched by a position operatorr. These Wannier state representations are shown diagrammatically in Fig. 5. For n = 1, the Berry connectionÃ(k) displays an extremumÃ max at the HSP k = 0, whereas for any n > 1, the extremum is located at a finite momentum k max , i.e., showing a double peak structure in momentum space. It is straight forward to obtain indicating a critical exponent α = 1/n. Moreover, the half-width-at-half-maximum scales as k HW HM ∝ |M | 1/n , implying the correlation length ξ ∝ |M | −1/n and a critical exponent ν = 1/n.

1D class AIII
For 1D class AIII that satisfies the chiral symmetry S 2 = 1, we follow Ref. [8]. The chiral operator is S = σ y , and hence the chiral symmetry demands that the d 2 term must vanish. We choose d 3 to be the mass term and write However, there is no restriction on the order of band crossing n ∈ Z being even or odd.
Since the Hamiltonian takes exactly the same form as that for 1D class BDI discussed in Sec. 5.1, the formalism therein also applies. We make the same rotation of basis by Eq. (123), and introduce the Berry connection and Zak phase. The resulting Wannier state correlation function takes the same form as in Eq. (128). The critical exponent of the correlation length is again ν = 1/n, except one should keep in mind that in 1D class AIII here n can be either even or odd.

1D class DIII
The 1D class DIII formalism can be obtained from 3D class DIII in Sec. 3.2 by a dimensional reduction k y = k z = 0. Since rotational symmetry is absent in 1D, the issue now is reduced to how to parametrize f = d 1 − id 2 and g = d 3 − iM in Eq. (26) using 1D momentum k x = k, provided d i (k) = −d i (−k) for i = 1, 2, 3. Given the oddness of d i , we should parametrize f = a 1 k n − ia 2 k n , g = a 3 k n − iM , n ∈ 2Z + 1 , such that the dispersion is The Fourier transform of which gives where |R1 denotes the Wannier state of band 1 centered at R, and |02 denotes the Wannier state of band 2 centered at the origin.

1D class D
For 1D class D with C 2 = 1, we follow Ref. [10] to use the PH operator C = σ x K, which requires that which is antisymmetric. The Z 2 topological invariant is calculated from the sign of the Pfaffian of X(k 0 ) at the two HSPs k 0 = 0 and k 0 = π in the original lattice model. In our continuous model so there is a topological phase transition at M c = 0. Alternatively, the topological invariant can be calculated from the Berry connection in the Majorana basis as in Sec. 5.1. The eigenstates in the Majorana basis calculated from Eq. (142) is The Berry connection of the filled band in this gauge is After using d 2 = k n , we obtain the same form of A(k) as in Eq. (125), and thus the same Zak phase C as in Eq. (126). We see that ∆C = 1 regardless of the order of band crossing n. The correlation function follows from Eq. (128) which measures the overlap of Wannier function sandwiched by a position operator. Because the basis of this Wannier state is the Majorana fermions, the correlation function has been referred to as the Majorana-Wannier state correlation function, as previously calculated for the 1D Kitaev chain which is an explicit example for n = 1 in this class [47]. It should be reminded that although we formulate this secion in the fermionic language, our formalism also applies to bosonic models.
The TR and PH operators are T = σ y ⊗ IK and C = I ⊗ τ y K. The TR and PH symmetries require the mass term to be d 3 = M , and the only allowed kinetic term is d 2 (k) = −d 2 (−k). Thus a general higher-order Dirac Hamiltonian in this class takes the form with n ∈ 2Z + 1. The topological invariant is given by the winding number [48] C = 1 2i 2π 0 dk 2π Tr σ y ⊗ τ y H −1 (k)∂ k H(k) , The filled-band eigenstates read To construct the correlation function, we consider a modified Berry connection such that the integrand in Eq. (149) can be written as meaning that the topological invariant can as well be expressed as an integration of Tr A ab . The Fourier transform of Tr A ab yields a correlation functioñ that measures the overlap of Wannier functions weighted by σ y ⊗ τ yr .

Conclusions
In summary, we investigate the quantum criticality of topological phase transitions for all the topologically nontrivial symmetry classes from 1D to 3D within a unified framework. Our approach generalizes the symmetry classification [7,8,9,10] to higherorder Dirac models in all physically relevant cases, and reveals the following features due to the interplay between the topological invariant C, the {T, C, S} symmetries, order of band crossing n, and the m-fold rotational symmetry in 2D and 3D: (1) The {T, C, S} symmetries constrain the order of band crossing to be even, odd, or integer in each dimension × symmetry class. (2) The even-oddness of the jump of topological invariant ∆C at the critical point may be due to the even-oddness of the band crossing, or may also be because the formula for the topological invariant C only allows certain integer values. (3) The even-oddness of the band crossing generally gives a constraint for the rotational eigenvalues. As a result, it is possible that only certain rotational symmetries m are compatible with a specific symmetry class, which may also depend on the system being fermionic or bosonic. (4) The topological invariant C (or a related one) takes the form of an integration over a curvature function. The Fourier transform of the curvature function always represents a correlation function that measures the overlap of two Wannier states that are a certain distance apart, or the overlap of multiple Wannier states in a compact form. (5) The critical exponent ν of the correlation length ξ, as well as that of the decay length of the edge state, are determined by the order of band crossing, but not necessarily ∆C or m. (6) The critical exponent ν and that of the extremum of the curvature function γ are not independent, but satisfy a scaling law owing to the conservation of the topological invariant. (7) The CRG approach based on the deformation of the curvature function can be used to judge topological phase transitions in all the dimensions and symmetry classes that have been investigated.
In the sense of clarifying how the {T, C, S, C m } symmetries influence the critical quantities {∆C, n, ν, γ}, the present work brings the notion of universality class into the research of topological phase transitions, as summarized in Table 1. The identification of Wannier state correlation function further indicates that, despite lacking a local order parameter, one may still construct a correlation function at the wave function level to characterize the quantum criticality of the system, both in the topologically trivial and nontrivial phase. In addition, the CRG approach demonstrates that the concept of scaling is a process of deforming local curvature that leaves the global topology unchanged, in contrast to the coarse graining process for Landau order parameters.A significant byproduct of our analysis is the complete classification of Dirac models with TR, PH, chiral, and m-fold rotation symmetry. Importantly, this is also of relevance for the study of band crossing points in semimetals or superconductors that are protected by m-fold rotation. That is, symmetry-protected band crossings in semimetals or superconductors with m-fold rotation symmetry described by the considered Dirac models with momentum-dependent mass terms. We anticipate that this framework we follow can be generally applicable to study how other kinds of spatial symmetries, in conjunction with the symmetry classification according to {T, C, S}, influence the quantum criticality of topological phase transitions. Applications of this kind to a wide range of spatial symmetries, such as reflection or a specific point group symmetry, await further investigations.