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

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

Wei Chen and Andreas P. Schnyder Department of Physics, PUC-Rio, Rio de Janeiro 22451-900, Brazil Max Planck Institute for Solid State Research, Stuttgart 70569, Germany
###### Abstract

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, -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 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.

## 1 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 at a critical point defines the topological phase transition. In practice, the low-energy effective theory of these materials are well described by a Dirac Hamiltonian that respects the symmetry of the system, and the tuning parameter represents the mass term, or equivalently the bulk gap of the energy spectrum, of the Dirac Hamiltonian. At the critical point , 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 . 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 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 , 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 and , 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 , 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 defines the critical exponent . It is further unveiled that the exponent is determined by the order of band crossing of the Dirac Hamiltonian at the critical point[6]. Depending on whether the bands at the critical point 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 -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) -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 even-oddness of the order of the band crossing in each symmetry class. We reveal that the order of band crossing in certain classes also determines the jump of topological invariant per gap-closing. When rotational symmetry is further imposed in 2D and 3D, the even-oddness of also constrains the rotational eigenvalues. As a result, the allowed rotational symmetry 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 , but not necessarily or . 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 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 -fold rotational symmetry in 2D and 3D, especially how they influence the critical quantities . 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.

## 2 Generic critical behavior of higher-order Dirac models

### 2.1 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 is constructed out of matrices that satisfies the Clifford algebra . In addition, the low-energy dispersion takes the form

 E±(k)=±d=±(k2n+M2)1/2, (1)

where the integer signifies the order of band crossing at the topological phase transition . The mass term is assumed to have no momentum dependence.

(2) The model is classfied according to the TR, PH and chiral symmetries

 TH(k)T−1=H(−k), CH(k)C−1=−H(−k), SH(k)S−1=−H(k). (2)

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 -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 -fold rotation symmetry can be written as

 Cm=diag(αp,αq...),αp=ei2πm(p+F2), (3)

where are the rotational eigenvalues. The spin statistical factor indicates whether the basis is fermionic or bosonic . We assume that in each symmetry class, the basis can be either fermionic or bosonic. The rotational symmetry requires that

 CmH(k)C−1m=H(Rmk), (4)

where rotates the momentum about the axis by . In 1D this rotational symmetry is absent.

As will be demonstrated case by case in the following sections, the topological invariant (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 -dimensional integration over the Brillouin zone (BZ)

 C=C(M)=∫BZdDk(2π)DF(k,M), (5)

where 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 (up to a reciprocal lattice vector) are referred to as the high symmetry points (HSPs). The curvature function is generally an even function around the HSP due to certain symmetry, such as inversion or TR symmetry. While the topological invariant remains constant within a topological phase in the parameter space of , the profile of varies with changing , which is key to our analysis of critical behavior.

Another key ingredient in our analysis is the Wannier states defined from the Bloch state by

 |unk⟩=∑Re−ik⋅(^r−R)|Rn⟩,|Rn⟩=1N∑keik⋅(^r−R)|unk⟩, (6)

where denotes the number of lattice sites, and is the position operator. We propose a correlation function derived through the Fourier transform of the curvature function,

 λR=∫dDk(2π)Deik⋅RF(k,M). (7)

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.

### 2.2 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 , the curvature function peaks at the gap-closing HSP. The peak gradually narrows and the height increases as the system approaches the critical point , and eventually the peak diverges and flips sign across the transition. We call this critical behavior the peak-divergence scenario. In this scenario, the peak can be well fitted by an Ornstein-Zernike form

 F(k0+δk,M)=F(k0,M)1+ξ2δk2. (8)

The critical behavior described above is summarized as

 limM→M+cF(k0,M)=−limM→M−1cF(k0,M)=±∞, limM→Mcξ=∞. (9)

Denoting the critical exponents of and as

 |F(k0,M)|∝|M−Mc|−γ,ξ∝|M−Mc|−ν, (10)

the conservation of the topological invariant under the diverging peak

 Cdiv=F(k0,M)(D∏i=1∫ξ−1−ξ−1dki1+ξ2k2i)=F(k0,M)ξD×O(1)=const., (11)

yields a scaling law

 γ=Dν, (12)

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 , the extremum of the curvature function forms a dimensional shell surrounding the HSP, with a radius denoted by and a thickness denoted by . The critical behavior is that the shell gradually reduces it radius as , and gradually diverges and flips sigh across the transition. In other words,

 limM→M+cFmax(kmax,M)=−limM→M−1cFmax(kmax,M)=±∞, limM→Mckmax=0,limM→Mckwid=0. (13)

For the Dirac models in Sec. 2.1, the and have the same critical exponent

 |F(k0,M)|∝|M−Mc|−γ,k−1max∝k−1wid∝ξ∝|M−Mc|−ν. (14)

The conservation of the portion of the topological invariant in the diverging shell

 Cdiv=Ω∫kmax+kwid/2kmax−kwid/2dkkD−1F(k,M)≈ΩkwidkD−1maxFmax=const. (15)

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

 C=Ω∫∞0dkkD−1MkD(n−1)(M2+k2n)(D+1)/2. (16)

Expanding the curvature function in the integrand for small yields

 F(k,M)=MkD(n−1)(M2+k2n)(D+1)/2≈Sgn(M)|M|DkD(n−1)(1+D+12k2nM2). (17)

For the peak-divergence scenario at , the Ornstein-Zernike form is evident, which renders and , satisfying Eq. (12). For the shell-divergence scenario at , solving for the from yields , and subsequently we obtain . Consequently, the critical exponents are and , 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 , 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 , and considering the edge state with zero transverse momentum, i.e., in 3D and in 2D. The problem can be calculated by projecting the higher-order Dirac Hamiltonian into real space and solve for the zero-energy edge state[6]

 (d⋅Γ)ψ=[knxΓx+MΓM]ψ=[(−i)n∂nxΓx+MΓM]ψ=0, (18)

where and represent the -matrices that multiply and , respectively, whose precise form depends on the symmetry class and dimension. Multiplying the equation by and using the ansatz , with an eigenstate of , we obtain

 ξn=in−2ηM. (19)

The eigenvalue is chosen such that one of the roots is real and identifiable with the edge state decay length, which obviously also has critical exponent . 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.

### 2.3 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 and rotational symmetry . The approach is essentially an iterative method to search for the trajectory (RG flow) in the parameter space of 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 , we find the new that satisfies

 F(k0+δk,M)=F(k0+bδk,M′). (20)

The mapping yields the RG flow that identifies the topological phase transitions. Here is a small momentum displacement away from the HSP, and one chooses for the peak-divergence scenario , whereas for the shell-divergence scenario . 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 region near the HSP which is set to be at the origin in our continuous models , we expand the scaling form in Eq. (17) by

 F(δk,M)≈Sgn(M)|M|DδkD(n−1)(1−D+12δk2nM2). (21)

For the peak-divergence scenario , putting Eq. (21) into the leading order expansion of the scaling equation, Eq. (20), yields the RG equation

 dMdℓ≡M′−Mδk2=1−b22∂2kF(k,M)|k=0∂MF(k,M)|k=0=(1−b2)(D+1)2D1M. (22)

The corresponding RG flow flows away from the critical point , with a flow rate that diverges as .

For the shell-divergence scenario at , we approximate Eq. (21) by

 F(δk,M)≈Sgn(M)|M|DδkD(n−1), F(bδk,M′)=F(δk+Δk,M+δM)≈Sgn(M)|M+δM|D(δk+Δk)D(n−1). (23)

Treating and as small quantities and expanding to the leading order, and then putting Eq. (23) into the Eq. (20) yields the leading order RG equation

 dMdℓ≡δM(Δk/δk)=(n−1)M. (24)

The equation corresponds to an RG flow that flows away from , with a flow rate that vanishes at . In other words, the topological phase transition manifests as an unstable fixed point at .

## 3 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 invariant.

### 3.1 3D class AIII

For 3D class AIII that has , we choose the -matrices[8]

 Γa=(αx,αy,αz,β,−iβγ5), αi=(0σiσi0),β=(100−1),γ5=(0110). (25)

The chiral operator is , which together with the chiral symmetry in Eq. (2) requires that , whereas there is no constraint on being even or odd in for . We follow the convention to choose as the mass term. The Hamiltonian then takes the form

 H(k)=∑i=1,2,3,5diΓi=⎛⎜ ⎜ ⎜⎝ccccg′ff∗−g′∗g′∗ff∗−g′⎞⎟ ⎟ ⎟⎠, (26)

where we have denoted

 f=d1−id2=kn++kn−−, g′=d3−iM=g−iM=knz−iM, (27)

with , such that the dispersion satisfies

 E±(k)=±d=±√|f|2+|g|2+M2=±√k2n+M2, (28)

with .

Using the rotational symmetry operator

 Cm=diag(αp,αq,αr,αs),αp=ei2πm(p+F2), (29)

the chiral symmetry does not give any constraint on . The rotational symmetry in Eq. (4) demands

 αpα∗sf(k±,kz)=αrα∗qf(k±,kz)=f(k±e±i2πm,kz), αpα∗rg′(k±,kz)=αsα∗qg′(k±,kz)=g′(k±e±i2πm,kz). (30)

Because contains the mass term that must be invariant under the rotational transformation, it is required that and , and therefore and . The constraint on then gives , which is consistent with our parametrization in Eq. (27). The parametrization of in Eq. (27), together with gives

 n+−n−=xm+(p−q)∈Z, (31)

which is the relation between the parametrization and the rotational symmetry eigenvalues, and we have used due to . The in Eq. (31) and throughout the article is an arbitrary integer. Evidently, Eq. (31) can be satisfied by any of the for either the bosonic or fermionic basis.

We define the -matrix by (not to be confused with one of the rotational eigenvalues)

 q(k)=−1d(diσi−iM)=−1d(g−iMff∗−g−iM). (32)

Writting the parametrization in Eq. (27) into cylindrical coordinate by and , and using the derivatives

 ∂∂kx=cosϕ∂∂k⊥−sinϕk⊥∂∂ϕ,∂∂ky=sinϕ∂∂k⊥+cosϕk⊥∂∂ϕ, (33)

yields the following combination in the integration of topological invariant[8]

 ϵμνρTr[q†∂μqq†∂νqq†∂ρq]=−M12k2(n−1)⊥kn−1z(n+−n−)n2(k2n+M2)2. (34)

Converting the integrand into spherical coordinates yields

 ϵμνρTr[q†∂μqq†∂νqq†∂ρq] =−M12(ksinθ)2(n−1)(kcosθ)n−1(n+−n−)n2(k2n+M2)2. (35)

The topological invariant is given by the 3D winding number obtained from the integration of the above quantity

 C = π3∫d3k(2π)3ϵμνρTr[q†∂μqq†∂νqq†∂ρq] (36) = {−(n+−n−)2Sgn(M)[2nB(n2,n)]ifn∈2Z+1,0ifn∈2Z,

where is the Beta function. We see that there is a topological phase transition if and only if the order of band crossing is an odd number. Here the extra factor is an artifact of this continuous model, and the jump of topological invariant should be treated as , as demonstrated by the following lattice model

 f=d1−id2=(sinkx−isinky)n−,g=d3=sinn−kz, d5=M+coskx+cosky+coskz. (37)

The resulting topological invariant is

 C = π3∫BZd3k(2π)3ϵμνρTr[q†∂μqq†∂νqq†∂ρq] (38) = {n−forM\apprge−3andn−∈2Z+1,0forM\apprle−3orn−∈2Z,

which gives .

The Wannier state correlation function in this class is introduced from the Chern-Simons gauge field defined from the eigenstates for the filled bands

 (39)

We find that the integrand in Eq. (36) can be written as

 ϵμνρTr[q†∂μqq†∂νqq†∂ρq]=8ϵμνρTr[AμAνAρ] =12Md4ϵμνρ∂xdμ∂ydν∂zdρ. (40)

As a result, the topological invariant in Eq. (38) is equivalently

 C=8π3∫BZd3k(2π)3ϵμνρTr[AμAνAρ]. (41)

In addition, the Chern-Simons gauge field can be expressed as the Fourier transform of the filled-band Wannier states

 Aabμ=−i∑Reik⋅R⟨0a|^rμ|Rb⟩=−i∑Re−ik⋅R⟨Ra|^rμ|0b⟩, (42)

where is the Wannier state of filled band localized at home cell . We use following compact notation for the eigenstates and the corresponding Wannier states

 |u−⟩=(|u1−⟩|u2−⟩),|R⟩=(c|R1⟩|R2⟩). (43)

It then follows that the topological invariant in Eq. (41) can be written in terms of the Wannier states entirely in real space

 C=i8π3∫d3R1∫d3R2ϵμνρTr[⟨0|^rμ|R1⟩⟨R1|(^rν−R1ν)|R1+R2⟩ ×⟨R1+R2|(^rρ−R1ρ−R2ρ)|0⟩], (44)

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

 =i∫d3R1∫d3R2ϵμνρTr[⟨0|^rμ|R1⟩⟨R1|(^rν−R1ν)|R1+R2⟩ ×⟨R1+R2|(^rρ−R1ρ−R2ρ)|R⟩]. (45)

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 , and hence the correlation length has critical exponent . Finally, the Chern-Simons invariant[8], equivalently the magnetoelectric polarizability, can as well be expressed in terms of Wannier states

 CS3=−π∫d3k(2π)3ϵμνρTr[Aμ∂νAρ+23AμAνAρ] =iπ∫d3R1ϵμνρTr[⟨0|^rμ|R1⟩⟨R1|R1ν(^rρ−R1ρ)|0⟩] −i2π3∫d3R1∫d3R2ϵμνρTr[⟨0|^rμ|R1⟩⟨R1|(^rν−R1ν)|R1+R2⟩ ×⟨R1+R2|(^rρ−R1ρ−R2ρ)|R⟩]. (46)

according to the same argument.

### 3.2 3D class DIII

For 3D class DIII that has and , we follow Ref. [8] and use the same -matrices as in Eq. (25). The TR and PH operators are and . It follows that the TR symmetry in Eq. (2) requires

 di(k)=−di(−k)fori=1,2,3,4,d5(k)=d5(−k). (47)

Thus only the term can be the mass term . On the other hand, the PH symmetry requires

 di(k)=−di(−k)fori=1,2,3,di(k)=di(−k)fori=4,5. (48)

Thus the term cannot exist, and the Hamiltonian takes the same form as Eq. (26), with the same parametrization on and as in Eq. (27). Given the rotational symmetry operator in Eq. (29), demanding and requires , .

The rotational symmetry in Eq. (4) requires that

 α2pf(k±,kz)=(α∗q)2f(k±,kz)=f(k±e±i2πm,kz), αpαqg(k±,kz)=α∗pα∗qg(k±,kz)=g(k±e±i2πm,kz), (49)

from which we see that , and hence the rotational operator is constrained to be the of form

 Cm=diag(αp,α∗p,αp,α∗p). (50)

The constraint on in Eq. (49) renders

 ei2πm(2p+F)kn++kn−−=ei2πm(n+−n−)kn++kn−−, (51)

which together with the requirement that must be odd in momentum yields

 n+−n−=xm+2p+F∈2Z+1, (52)

We see that the fermionic case can be satisfied by any of the , but bosonic case can only be realized at . 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).

### 3.3 3D class AII

For 3D class AII that has , we use the same representation of the -matrices as in Eq. (25), and the TR operator following Ref. [8]. The TR symmetry requires that

 di(k)=−di(−k)fori=1,2,3,5,d4(k)=d4(−k), (53)

and hence only 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 , renders and , and consequently . We may write the Hamiltonian into the form

 H=⎛⎜ ⎜ ⎜⎝MgfMf∗−g∗g∗f−Mf∗−g−M⎞⎟ ⎟ ⎟⎠, (54)

where and . The rotational symmetry in Eq. (4) demands that

 αpαrf(k±,kz)=f(k±e±i2πm,kz), αpα∗rg(k±,kz)=g(k±e±i2πm,kz). (55)

We parametrize and by

 f=kn++kn−−,g=knz. (56)

With this parametrization in mind, the constraint on in Eq. (55) demands that , 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 invariant for 3D lattice models in class AII can be constructed from the Pfaffian of the filled-band -matrix at HSPs[36, 37], provided the Pfaffian is a real number[36, 37]

 (−1)C=∏iSgn(Pf[m(k0,i)]), (57)

where