Piezoelectricity and topological quantum phase transitions in two-dimensional spin-orbit coupled crystals with time-reversal symmetry

Finding new physical responses that signal topological quantum phase transitions is of both theoretical and experimental importance. Here, we demonstrate that the piezoelectric response can change discontinuously across a topological quantum phase transition in two-dimensional time-reversal invariant systems with spin-orbit coupling, thus serving as a direct probe of the transition. We study all gap closing cases for all 7 plane groups that allow non-vanishing piezoelectricity, and find that any gap closing with 1 fine-tuning parameter between two gapped states changes either the Z2 invariant or the locally stable valley Chern number. The jump of the piezoelectric response is found to exist for all these transitions, and we propose the HgTe/CdTe quantum well and BaMnSb2 as two potential experimental platforms. Our work provides a general theoretical framework to classify topological quantum phase transitions, and reveals their ubiquitous relation to the piezoelectric response.


Numbers of FTPs and Effective
In this section, we derive Eq. (3) in the main text via linear response theory from Eq. (2) in the main text, which is equivalent to the derivation in Supplementary Ref. [1]. The derivation is done with the natural unit c = = 1 and the metric (−, +, +).
To apply the linear response theory, we start from an action S that includes the electronic effective model and the leading order effect of the infinitesimal strain. Since the current is present in Eq.
(2) , we should include the U (1) gauge field that accounts for the electromagnetic field. With the U (1) gauge field, the action reads where k µ = (ω, k) µ , A µ and u ij and ψ follow the same Fourier transformation rule, G 0 (k) = [ω − h 0 (k)(1 − i )] −1 is the time-ordered Green function without the electron-strain coupling, the chemical potential is chosen to be the zero energy, and M ij is the matrix coupled to the strain tensor u ij . To the leading order, the linear response is given by the following effective action where f ij,µν = − 1 2 and the absence of the Chern-Simons term AdA is due to the T symmetry. With Supplementary Eq. (2) and Eq. (2) , we can use the condition that u ij is uniform to derive the expression of the PET, resulting in To further derive Eq. (3) , we define h(k, u ij ) = h 0 (k) + u ij M ij and G(k, u ij ) = [ω − h(k, u ij )(1 − i )] −1 as the Hamiltonian and Green function with the electron-strain coupling, respectively. Using ∂ kµ G −1 = ∂ kµ G −1 0 and ∂ uij G −1 = −M ij , we can revise Supplementary Eq. (4) to Define X µ = (ω, k i , u jk ) and then the above equation can be further transformed to where µνρ is the Levi-Civita symbol. Integrating out ω in the above equation with the Wick rotation gives Eq. (3) . Although the derivation here is done for = c = 1, all the expressions of γ ijk and the resultant Eq. (3) stay the same after converting to the SI unit as they carry the right unit for the PET in 2+1D.
Finally, we would like to discuss the effect of the identity term of h 0 in Supplementary Eq. (2) when h 0 is a two band model. In general, the Hamiltonian can always be split into the identity part and the traceless part as h 0 (k) = m 0 (k)1 + h traceless 0 (k). The eigenvalues of h 0 (k) then read m 0 (k) ± ε(k), where ±ε(k) are two eigenvalues of h traceless 0 (k) with ε(k) > 0 chosen without loss of generality. As the model is gapped and the Fermi energy (E = 0) is chosen to lie inside the gap, we have ε(k) > |m 0 (k)| ≥ 0. Since the poles of G 0 are at ω = [m 0 (k) ± ε(k)](1 − i ), integrating ω along (−∞, ∞) in f ij,µν of Supplementary Eq. (2) gives the same result as integrating ω along (−∞ + m 0 (k)(1 − i ), ∞ + m 0 (k)(1 − i )) owing to the absence of poles in between the two paths. As a result, we can directly neglect the identity term of a two-band insulating h 0 in f ij,µν of Supplementary Eq. (2).

Supplementary Note 2. Details on PET for Each PG
The discussion on the electronic effective model and FTP of the gap closing between two non-degenerate states has some overlap with Supplementary Ref. [2]. However, the topological property and PET jump of the gap closing between two insulating states have not been discussed in Supplementary Ref. [2].

A. PG p1
In the main text, the effective Hamiltonian for scenario (ii) of p1 is derived in a non-Cartesian coordinate system, which is not convenient for the generalization to other PGs with more crystalline symmetries. Thus, we re-derive the effective Hamiltonian in the Cartesian coordinate system, as given by (see Supplementary Note 3 A) Here we only perform the unitary transformation on the bases of the Hamiltonian and do not rotate the momentum or the coordinate system. Correspondingly, the PET jump across the direct TQPT at m = 0 can be derived as Supplementary Eq. (7)-(8) resemble the conclusion for p1 in the Results and are useful for the discussion of the other 6 PGs with non-vanishing PET. We would like to discuss more about the GC and PET for p1. In the first scenario, all TRIM have no essential differences and the gap closing always happens between two Kramers pairs unless more parameters are finely tuned. Therefore, there is no need to further classify this scenario into finer cases, and the codimension for the gap closing is 5, indicating that this scenario cannot be direct TQPT 3 . According to the main text, no finer classification is needed for the second scenario either, the codimension of the gap closing scenario is 1, and it is indeed a direct TQPT that changes the Z 2 index and leads to the PET jump.

B. PGs p1m1, c1m1 and p1g1
In this part, we study three PGs, p1m1, c1m1, and p1g1, all of which are generated by a mirror-related symmetry U and the lattice translation. U is a mirror operation for p1m1/c1m1 and a glide operation for p1g1. The difference between p1m1 and c1m1 lies on the directions of the primitive lattice vectors relative to the mirror line, which is not important for our discussion here. Without loss of generality, we choose the mirror or glide line to be perpendicular to x, labelled as m x or g x , respectively. The glide operation is thus denoted as g x = {m x |0 1 2 }, where "0 1 2 " represents the translation by half the primitive lattice vector along y. The U symmetry in these three PGs requires γ xxx = γ xyy = γ yxy = γ yyx = 0, whereas the PET components γ xxy , γ xyx , γ yxx , γ yyy are allowed to be nonzero. For the symmetry analysis here, the PET behaves the same under the glide and mirror operations since u ij is considered in the continuum limit.
In order to classify the gap closing scenarios, we define the group G 0 for a gap closing momentum k 0 such that G 0 contains all symmetry operations that leave k 0 invariant. Since G 0 can include the TR-related operation, it can be larger than the little group of k 0 . Based on G 0 , we obtain in total 4 gap closing scenarios for these three PGs: (i) the gap closing at TRIM (G 0 contains T ), (ii) G 0 contains U but not T , (iii) G 0 contains UT but not T , (iv) G 0 contains no symmetries other than the lattice translation, which we refer to as the trivial G 0 . As summarized in Tab. 1 in the main text, the TQPT exists in scenario (iii) and (iv), which can lead to the jump of symmetry-allowed PET components.

I. Scenario (i): TRIM
In scenario (i), the gap closing requires 3 (5) fine-tuning parameters (FTPs) for p1m1 and c1m1 if m x is (is not) in the G 0 . (See Supplementary Note 3 B.) For p1g1, the TRIM (Γ, X, Y and M ) are split into two classes according to the value of g 2 x : Γ, X with g 2 x = −1 and Y, M with g 2 x = 1. The gap closing at Γ, X needs 3 FTPs since g x behaves the same as m x , while the gap closing at Y, M needs only 1 FTP if it happens between two Kramers pairs with opposite g x eigenvalues. However, such gap closing at Y, M is in between two g x -protected gapless phases with codimension 0, where the bands with opposite g x eigenvalues cross with each other at momenta other than Y, M as shown in Supplementary Fig. 1. Therefore, there is no direct TQPT between two gapped phases in scenario (i).
The same situation occurs for scenario (ii). In scenario (ii), the gap closes at two different momenta ±k 0 that are invariant under the U operation, meaning that the bases at ±k 0 can have definite U eigenvalues. The gap closing between the two bases with the same U eigenvalues requires 2 FTPs, as discussed in Supplementary Note 3 B. When the gap closes between two bands with opposite U eigenvalues, the system always enters a stable U-protected gapless phase with 0 codimension. (This case is not the same as the scenario (i) since only one side is guaranteed to be gapless.) Thus, the gap closing cases cannot be direct TQPTs.
In scenario (iii), the gap closing occurs at two different momenta ±k 0 that are invariant under UT , as shown by the orange dashed lines in Fig. 1 (b) and (c) in the main text.
For p1m1 and c1m1 (U = m x ), (m x T ) 2 = 1 suggests that we can have m x T=K at k 0 by choosing the appropriate bases and the band touching point at k 0 should typically occur between two non-degenerate bands. We further take T=iσ y K by choosing the appropriate bases at −k 0 , and thus the two-band effective models h ± (q, u) at ±k 0 can be given by Supplementary Eq. (7) with extra constraints v 0 = ξ a1,xy = ξ a1,yx = ξ y,xx = ξ y,yy = 0 (9) for a 1 = 0, x, z. As a result, only 1 FTP m is needed for the gap closing (m = 0), and only one single Dirac cone exists in half 1BZ at the transition, leading to the change of the Z 2 index. Based on Supplementary Eq. (8), the jump of symmetry-allowed PET components across this TQPT can be derived as For p1g1 with U = g x , since (g x T ) 2 = 1 at (k x , 0) and (g x T ) 2 = −1 at (k x , ±π), we have two different gap closing cases. When the gap closes at (±k 0,x , 0), the algebra relation involving g x T is the same as m x T , e.g. (g x T ) 2 = (m x T ) 2 = 1, and thus the effective Hamiltonian can be chosen to be the same as that for p1m1 and c1m1, leading to 1 FTP, Z 2 index change, and the same form of PET jump. On the contrast, due to (g x T ) 2 = −1 at (±k 0,x , ±π), the gap closing needs 4 FTPs and thus no TQPT can occur in this case. (See Supplementary Note 3 B.)
There is no change of Z 2 index for this scenario, since two Dirac cones exist in half 1BZ when the gap closes and the CN of contracted half 1BZ can only change by an even number. Nevertheless, scenario (iv) can still be "topological" in the context of valley Chern number (VCN) as elaborated in the following. Due to the Dirac Hamiltonian form shown in Supplementary Eq. (7), the Berry curvature is peaked at each valley k 0,1,2,3 for a small m and can be captured by the electronic part of the corresponding effective Hamiltonian. Then, we can integrate the Berry curvature given by the effective model and get the VCN 4,5 for each valley as N ki = −η i sgn(v x v y )sgn(m)/2 with i = 0, 1, 2, 3. The values of η i at different valleys are related by the TR and U symmetries, both of which flip the sign of the Berry curvature. Thus, we have η 0 = η 3 = 1 and η 1 = η 2 = −1. It should be pointed out that the Berry curvature integral is not over the entire 1BZ and the VCN at each valley thus does not need to be an integer. Nevertheless, the change of VCN across the gap closing is defined on a closed manifold and must be an integer number, given by ∆N ki = −η i sgn(v x v y ) as varying m from 0 − to 0 + . For the convenience of further discussion, we can define the VCN of the whole system 5 as N val = i η i N ki = −2sgn(v x v y )sgn(m), and the change of the VCN becomes ∆N val = −4sgn(v x v y ) = 4∆N + with the factor 4 for the four valleys. Therefore, if we restrict all the valleys to be far apart in the momentum space, the change of the VCN is a well-defined topological invariant and this gap closing scenario is a TQPT.
In principle, tuning parameters may merge different valleys at some high symmetry momentum, e.g. the valleys at k 0 and Uk 0 merged at the mirror or glide line. Therefore, without the constraint of well-defined valleys, two phases with different VCNs can share the same band topology and thus can be adiabatically connected. It means the topology characterized by VCN is "locally stable" 6 , though globally unstable. Nevertheless, we restrict all valleys to be well-defined in our discussion and refer to the gap closing scenario as a TQPT.
Next we study the change of the PET components at this TQPT, which can be split into two parts: ∆γ (0) originating from ±k 0 and ∆γ (1) given by ±Uk 0 . ∆γ (0) equals to Supplementary Eq. (8) since the effective models at ±k 0 are the same as Supplementary Eq. (7). Owing to the mirror or glide symmetry, ∆γ (1) is related to ∆γ (0) as As a result, we obtain the non-zero jump of symmetry-allowed PET components ∆γ ijk = ∆γ C. PG p3 PG p3 is generated by 3-fold rotation C 3 and the lattice translation. Owing to C 3 , the PET only has two independent components γ xxx and γ yyy as Again, we classify the gap closing for p3 according to G 0 , resulting in three different scenarios: (i) G 0 contains T , (ii) G 0 contains C 3 but not T , and (iii) G 0 is trivial. Here we do not have a scenario for G 0 containing C 3 T but no T , since (C 3 T ) 3 is equivalent to T . As summarized in Tab. 1 in the main text and elaborated in the following, in any of the above scenarios, there are gap closing cases between gapped states that need only 1 FTP, change the Z 2 index, and lead to the discontinuous change of symmetry-allowed PET components.

I. Scenario (i):TRIM
There are 4 TRIM in scenario (i), namely three M points related by C 3 and one Γ point, as labeled in Fig. 1 (f) in the main text. G 0 of each individual M point only contains T and the lattice translation, and thus the gap closing at M needs 5 FTPs, same as the gap closing at TRIM for p1.
When the gap closes at Γ point as shown in Fig. 1 (f) in the main text, G 0 also contains C 3 with C 3 3 = −1. Due to [C 3 , T ] = 0, the Kramers pairs can be classified into two types according to the C 3 eigenvalues: one with (e −iπ/3 , e iπ/3 ) and the other with (−1, −1). The gap closing between the Kramers pairs of the same type requires more than 1 FTPs, 3 for (e −iπ/3 , e iπ/3 ) type and 5 for (−1, −1) type, as discussed in Supplementary Note 3 C.
The gap closing with 1 FTP happens between the TR pairs of different types, for which the minimal four-band effective Hamiltonian in the bases (e −iπ/3 , e iπ/3 , −1, −1) reads where h p3,0 is the electron part and h p3,1 describes the electron-strain coupling τ 's and σ's are Pauli matrices that label two different Kramers pairs and two components of each Kramers pair, respectively, m is the gap closing tuning parameter, and the bases are chosen such that T= − iτ 0 σ y K. This gap closing is certainly a TQPT since it changes the numbers of IRs of the occupied bands, meaning that the two gapped states separated by this gap closing cannot be adiabatically connected. When v 1 = v 2 = 0, we can define an effective inversion symmetry P = τ z σ 0 for the electron part of Supplementary Eq. (13), P h p3,0 (−k) P † = h p3,0 (k), and thus the gap closing of h p3,0 (k) with v 1 = v 2 = 0 changes the Z 2 index according to the Fu-Kane criteria 7 since the parity of the occupied band changes. The existence of non-zero v 1 , v 2 terms that break P cannot influence the Z 2 topology change, since (i) the Z 2 topology does not rely on the effective inversion symmetry, and (ii) additional gap closing away from Γ is forbidden at m = 0 as long as the v 1,2 terms are restored adiabatically. Therefore, within a certain range of v 1,2 , the codimension-1 gap closing at m = 0 is a direct TQPT that changes the Z 2 index.
The remaining question is if the codimension-1 gap closing at m = 0 is always a Z 2 transition. To answer this question, note that we can always assume the transition at m = 0 is Z 2 for a parameter region S 1 of v i 's in h p3,0 and non-Z 2 for the other parameter region S 0 of v i 's. Since the same form of the Hamiltonian can not correspond to Z 2 and non-Z 2 transitions simultaneously, the intersection of S 0 and S 1 is empty. Now we suppose both S 0 and S 1 are codimension-0 subspaces of the v i parameter space (not the whole parameter space since only v i 's are included while m is excluded). Then, the boundary of S 0 , labeled as ∂S 0 , is a codimension-1 subspace of v i parameter space, and the special transition at (m = 0, v i ∈ ∂S 0 ) is a codimension-2 transition, as shown in Supplementary Fig. 3.
Patching the Hamiltonian with m = ± with positive and infinitesimal gives a Hamiltonian that lives on a closed manifold. This Hamiltonian has Z 2 trivial and nontrivial ground states when v i 's are in S 0 and S 1 , respectively, as shown in Supplementary Fig. 3. As v i 's change from S 0 to S 1 passing through ∂S 0 , the patched Hamiltonian must experience a gap closing at generic k points that changes Z 2 ( Supplementary Fig. 3), since there is always an energy gap at Γ for m = 0. As discussed with more detail in the following (see scenario (iii)), the gap closing at generic k points surely changes the Z 2 index, is codimension-1, and simultaneously happens at six momenta. The gap closing can only happen either for m = part or for m = − part of the patched Hamiltonian but not both, since if the gap closes twice, the Z 2 index would be changed back. It means, the codimension-1 hypersurface for the gap closing at generic k (red line in Supplementary Fig. 3) touches the codimension-1 hypersurface for gap closing at Γ (m = 0 line in Supplementary Fig. 3) just from one side of m but not passing through. As mentioned above, the touching part at (m = 0, v i ∈ ∂S 0 ) is a codimension-2 transition, owing to the assumption that both S 0 and S 1 are codimension-0 subspaces of v i parameter space.
At the touching, we must have the six generic gap closing points merging at Γ. Otherwise, we should expect the red line in Supplementary Fig. 3 to pass through the m = 0 line instead of stopping, since the gap closing process is local in the momentum space and different gap closing cases cannot influence each other if they happen at the different momenta. However, the merging process cannot be codimension-2 since moving a generic gap closing point Fig. 3. This figure shows the typical phase diagram around (m = 0, vi ∈ ∂S0) (the black dot), within the assumption that both S0 and S1 are codimension-0 subspaces of the vi parameter space. The green line is (m = 0, vi ∈ S0) which does not change Z2 index, the blue line is (m = 0, vi ∈ S1) that changes Z2 index, and the system closes the gap at six generic k points on the red line. Without loss of generality, we choose the red line to touch the m = 0 line from the positive m side. A,B, and C label three different Hamiltonians given by patching the two effective models with m = ± at different values of vi. A is Z2 trivial, C is Z2 non-trivial, and B closes the gap at the six generic k points for in m = part.
to a specific momentum while keeping the gap closed requires at least 3 FTPs (two to move the momentum and one to close gap). Therefore, S 0 and S 1 cannot be both codimension-0 subspaces of v i parameter space, and the codimension-1 gap closing at Γ can only be Z 2 or non-Z 2 but not both. Since we already show that the Z 2 transition at Γ can be codimension-1, the codimension-1 gap closing at Γ should always change the Z 2 index.
The gap closing momenta in scenario (ii) are K and K in Fig. 1 (g) in the main text. Since these two momenta are related by T , we only need to derive the effective model at one momentum, say K, and the other one can be obtained using T . At K, the C 3 symmetry has three possible eigenvalues −1, e ±iπ/3 due to C 3 3 = −1. If the gap closing is between two states with the same C 3 eigenvalues, it cannot be TQPT since the fixed gap closing momentum leads to 3 FTPs for the gap closing. (See Supplementary Note 3 C.) There are three cases for two states with different C 3 eigenvalues: (e −iπ/3 , e iπ/3 ), (e iπ/3 , −1), and (−1, e −iπ/3 ). The effective models in the three cases are equivalent since the representations of C 3 in these cases can be related to each other by multiplying a phase factor e ±i2π/3 . Therefore, we focus on the first case, of which the effective model at K (after an appropriate unitary transformation) is given by h + in Supplementary Eq. (7) with where a = 0 or z. Similarly, by choosing the appropriate bases at K such that T=iσ y K, the effective model at K is given h − in Supplementary Eq. (7) with the parameter relation listed above. As a result, the gap closing between states with different C 3 eigenvalues needs 1 FTP and changes the Z 2 index since half 1BZ contains one Dirac cone (K or K ). Based on Supplementary Eq. (8), the jump of independent PET components across this TQPT (varying m from 0 − to 0 + ) has the non-zero form where ∆N + = −sgn(v 2 ) = −1.

III. Scenario (iii): trivial G0
In scenario (iii), there are six gap closing momenta, labeled as ±k 0 , ±C 3 k 0 and ±C 2 3 k 0 , as shown by red crosses in Fig. 1 (h) in the main text. The effective Hamiltonian at ±k 0 are exactly the same as Supplementary Eq. (7) since the two momenta are related by T and no more symmetries are involved. Therefore, the gap closing scenario needs 1 FTP, and the contribution to the PET jump from the gap closing at ±k 0 is the same as Supplementary Eq. (8), noted as ∆γ (0) ijk . The effective models at ±C 3 k 0 and ±C 2 3 k 0 can be obtained from those at ±k 0 by C 3 and C 2 3 operations, respectively, whose electronic parts are also in the Dirac Hamitlonian form. The contracted half 1BZ then contains three Dirac cones at the gap closing and its CN must change by an odd number, indicating the change of Z 2 index. Furthermore, the contributions to the jump of PET components from the gap closing at ±C 3 k 0 and ±C 2 3 k 0 are ∆γ i j k and ∆γ i j k , respectively, owing to the symmetry. As a result, the jump of independent PET components is given by ∆γ ijk = ∆γ (0) ijk + ∆γ (1) ijk + ∆γ (2) ijk , which has the nonzero form D. PG p31m and PG p3m1 Both PGs p31m and p3m1 are generated by the lattice translation, the three-fold rotation C 3 , and a mirror symmetry which we choose to be m x without loss of generality. The difference between the two PGs lies on the direction of the mirror line relative to the primitive lattice vector: the mirror line is parallel or perpendicular to one primitive lattice vector for p31m or p3m1, respectively. C 3 and m x span the point group C 3v , which leads to for the PET, and thus γ yyy serves as the only independent symmetry-allowed PET component. We classify the gap closing scenarios into 4 types according to G 0 : (i) G 0 contains T , (ii) G 0 contains at least one of the three mirror symmetry operations in C 3v (again labeled as U = m x , C 3 m x , or C 2 3 m x ) but no T , (iii) G 0 contains the UT but no T , and (iv) G 0 is trivial. As summarized in Tab. 1 in the main text, all gap closing cases between gapped states with 1 FTP change either Z 2 index or the VCN, and lead to the jump of symmetry-allowed PET components. When the gap closes at Γ point ( Fig. 1 (i-j) in the main text), the generators of G 0 besides the lattice translation are C 3 , m x and T , and there are still two types of Kramers pairs characterized by the C 3 eigenvalues as those in Supplementary Note 2 C I. Owing to the extra mirror symmetry here, the number of FTPs for the gap closing between the same type of Kramers pairs becomes 2 for (e −iπ/3 , e iπ/3 ) type and 3 for (−1, −1) type as discussed in Supplementary Note 3 D. Therefore, we still only need to consider the gap closing between different types of Kramers pairs. For the convenience of the later material discussion, we choose the bases as (e −iπ/3 , −1, e iπ/3 , −1). One can always choose the TR symmetry and mirror symmetry to be represented as T= − iσ y τ 0 K and m x= − iσ x τ 0 . In this case, the effective Hamiltonian can be derived by imposing the m x on Supplementary Eq. (13), leading to The form of the Hamiltonian then reads where u 2 = u xx + u yy , u ± = u xx − u yy ± i(u xy + u yx ). The above Hamiltonian shows that the gap closing at Γ needs only 1 FTP, which is m. As discussed in Supplementary Note 3 D, this gap closing cannot drive a gapped phase into a mirror-protected gapless phase, and therefore can separate two gapped states. Similar to the discussion in Supplementary Note 2 C I, the gap closing changes the Z 2 index when tuning m from 0 − to 0 + , indicating a TQPT. When v 2 = 0, an analytical expression for the jump of independent PET component can be obtained from Supplementary Eq. (16) and Supplementary Eq. (21), which reads With parameter values v 3 = v 6 = 1eVÅ and ξ 1 = ξ 2 = 2ξ 4 = ξ 5 = 1eV , the numerical results ( Supplementary  Fig. 2(b)) for non-zero v 2 still show a PET jump across TQPT.
II. Scenario (ii): U ∈ G0 and T / ∈ G0 Scenario (ii) can be further divided into two classes depending on whether G 0 contains C 3 . When G 0 does not contain C 3 , the gap closing either requries more than 1 FTP or drives the system into a mirror-protected gapless phase with 0 codimension, similar to Supplementary Note 2 B II.
Only when the gap closes at K, K for p31m, G 0 contains C 3 . In this case, G 0 contains the group C 3v , which has one 2D irreducible representation (IR) and two different 1D IRs when acting on the states. The gap closing between the states furnishing the same IR requires 3 FTPs, similar to the case for two states with the same C 3 eigenvalue in Supplementary Note 2 C II. If the gap closes between the doubly degenerate states furnishing the 2D IR and a state furnishing a 1D IR, the system with a fixed carrier density cannot be insulating on both sides of the gap closing because the number of occupied bands is changed. If the gap closes between two states that furnish different 1D IRs, the mirror-protected gapless phase must exist on one side of the gap closing as the two states must have opposite mirror eigenvalues. Therefore, there is no direct TQPT between the insulating phases in scenario (ii).

III. Scenario (iii): UT ∈ G0 and T / ∈ G0
In scenario (iii), the gap closing cases are again divided into two different classes depending on whether G 0 has C 3 . We first discuss the class without C 3 , which happens for the gap closing at UT invariant momenta except K, K for p3m1. As shown in Fig. 1 (k-m) in the main text, the total number of inequivalent gap closing momenta is six, including ±k 0 , ±C 3 k 0 , and ±C 2 3 k 0 . Without loss of generality, we choose k 0 such that −m x k 0 is equivalent to k 0 . Then, the effective models at ±k 0 are the same as the corresponding models in Supplementary Note 2 B III, i.e. Supplementary Eq. (7) with the parameter relation Supplementary Eq. (9), indicating 1 FTP for the gap closing. Since the effective models at ±C 3 k 0 and ±C 2 3 k 0 are related to those at ±k 0 by C 3 and C 2 3 operations similar to Supplementary Note 2 C III, the jump of PET components can be derived by substituting Supplementary Eq. (9) into Supplementary Eq. (19), resulting in Moreover, since three Dirac cones exist in half 1BZ when the gap closes, the Z 2 index changes at the gap closing, making it a TQPT. The class that G 0 includes C 3 can only happen when the gap closes at K and K for PG p3m1, as shown in Fig. 1  (m) in the main text. We can choose C 3 and m x T as the generators of G 0 besides the lattice translation. Similar to Supplementary Note 2 C II, we first study K and derive the model at K by choosing the right bases such that T=iσ y K. The states at K can be labeled by C 3 eigenvalues, −1 and e ±iπ/3 given by C 3 3 = −1. Since (m x T ) 2 = 1 and C 3 m x T = m x T C −1 3 , the gap closing typically happens between two non-degenerate states, labeled by the C 3 eigenvalues as (λ 1 , λ 2 ), and we can always choose m x T=K. The λ 1 = λ 2 case cannot correspond to TQPT since 2 FTPs are needed for the gap closing as discussed in Supplementary Note 3 D, while the λ 1 = λ 2 case requires only one FTP for the gap closing similar to Supplementary Note 2 C II. Since the matrix representations of C 3 and m x T are equivalent for the three choices (λ 1 , λ 2 ) = (e −iπ/3 , e iπ/3 ), (e iπ/3 , −1), and (−1, e −iπ/3 ), they have the same effective models and we only consider the first choice. With all the above conventions and simplifications, the effective models at K and K can be given by those for Supplementary Note 2 C II with an extra constraint ξ y,yy = 0 brought by m x T . As a result, the Z 2 index does change when the gap closes, and the jump of PET components can be derived from Supplementary Eq. (18) with the above extra constraint, which reads IV. Scenario (iv): trivial G0 In scenario (iv), the gap closes simultaneously at twelve inequivalent momenta, namely ±k 0 , ±m x k 0 , ±C 3 k 0 , ±C 3 m x k 0 , ±C 2 3 k 0 and ±C 2 3 m x k 0 in Fig. 1 (n-o) in the main text. The effective model around k 0 can be chosen as h + in Supplementary Eq. (7), and the models around other gap closing momenta can be further obtained by the symmetry. Although this gap closing scenario only needs 1 FTP, it cannot induce any change of Z 2 index since there is an even number (six) of Dirac cones in contracted half 1BZ. However, the gap closing can change the VCN when the twelve valleys are well defined according to Supplementary Note 2 B IV, e.g. N k0 can change by ±1, and thus is a TQPT in the sense of the locally stable topology.
We split the change of PET components for this scenario into 3 parts: γ (0) from ±k 0 and ±m x k 0 , γ (1) from ±C 3 k 0 and ±C 3 m x k 0 , and γ (2) from ±C 2 3 k 0 and ±C 2 3 m x k 0 . Since the contribution to γ (0) contains two Kramers pairs that are related by m x , same as Supplementary Note 2 B IV, γ (0) equals to Supplementary Eq. (11). C 3 symmetry then gives ∆γ i j k and ∆γ (2) i j k , similar to Supplementary Note 2 C III. As the result, the total change of PET can be obtained from ∆γ = ∆γ (0) + ∆γ (1) + ∆γ (2) , which is propotional to the change of the VCN of the system with ∆N val = 12∆N + .

E. 10 PGs with 2D Inversion or C 2
The PET jump cannot exist in 10 PGs that contain C 2 or inversion, including p2, p2mm, p2mg, p2gg, c2mm, p4, p4mm, p4gm, p6, and p6mm. This conclusion can be drawn from the symmetry analysis of PET. Since both C 2 and inversion transform (x, y) to (−x, −y), γ ijk = −γ ijk is required for those 10 PGs, leading to the vanishing PET. Early study 6,8 also shows that a stable gapless phase can exist in between the QSH insulator and the NI when C 2 exists. In this gapless regime, 2D gapless Dirac fermions are locally stable and can only be created or annihilated in pairs.

Supplementary Note 3. Numbers of FTPs and Effective Models for the Gap Closing
The discussion on the gap closing between two non-degenerate states has some overlap with Supplementary Ref. [2].

A. PG p1
This part has been studied in Supplementary Ref. [3].

I. Not TRIM
When the gap closes at k 0 that is not a TRIM, the two-band model near the gap closing to the leading order of q = k − k 0 in general takes the form if C x q x + C y q y + M = 0. We choose C x × C y = 0 since it can be satisfied without finely tuning anything (or equivalently in a parameter subspace with 0 codimension). In this case, the gap closes when M lies in the plane spanned by two vectors C x and C y . Therefore, the codimension for the gap closing is 1 since only the angle between the vector M and the (C x , C y ) plane needs to be tuned. Next, we derive Eq.
Here C x × C y = 0 gives non-zero v x and v y . With a shift of k 0 by ( , the model is further simplified to the electronic part of h + in Supplementary Eq. (7). Finally, we define the q x + v 0 q y and q y to be q 1 and q 2 , respectively, to get Eq. (5) , which represents the most generic form of the Hamiltonian.

II. TRIM
In this part, we count the number of FTPs for the gap closing at the TRIM. Owing to the Kramers' degeneracy, every band at the TRIM is doubly degenerate, and we use the name "Kramers pair" to label the two states related by the TR symmetry. We consider the gap closing between two Kramers pairs |1, ± and |2, ± , where T |i, + = −|i, − can always be chosen by the unitary transformation. As a result, the mass term for the effective model at the TRIM reads where the bases are (|1, + , |1, − , |2, + , |2, − ) and all parameters are real. Since the momentum is fixed at TRIM (−k = k + G with G a reciprocal lattice vector), none of the terms in the above equation can be canceled by shifting the momentum. Therefore, 5 FTPs are needed for the gap closing.

I. Scenario (i): TRIM
If G 0 does not contain U, which can occur on the edge of 1BZ for c1m1, the situation is the same as the TRIM in Supplementary Note 3 A, which requires 5 FTPs. When G 0 contains U, we should discuss the U = m x case (p1m1 and c1m1) and the U = g x case (p1g1), separately.
For p1m1 and c1m1, since m 2 x = −1, two states of one Kramers pair have opposite mirror eigenvalues ±i, labeled by m x |i, ± = ±i|i, ± . On the bases (|1, + , |1, − , |2, + , |2, − ), the effective model around the gap closing between two Kramers pairs can be given by Supplementary Eq. (29) with ∆ 1 = ∆ 2 = 0, since the bases with different mirror eigenvalues cannot be coupled by the mass terms. As a result, 3 FTPs are needed for such gap closing scenario.
For p1g1, g 2 x = −1 at Γ and X and the number of FTPs for the gap closing is thus the same as the above case, which is 3. At Y and M , g 2 x = 1 and two states of one Kramers pair have the same g x eigenvalue, 1 or −1. In this case, the gap closing between two Kramers pairs with the same g x eigenvalue needs 5 FTPs, which is the same as the TRIM scenario in Supplementary Note 3 A. On the other hand, between two Kramers pairs with opposite g x eigenvalues, only 1 FTP needs to be tuned to close the gap at Y or M , since the off-diagonal terms (∆ 0,1,2,3 ) in Supplementary Eq. (29) are all forbidden.
In scenario (ii), there are two gap closing momenta ±k 0 that are related by the TR symmetry. Therefore, we only need to consider one of them, say k 0 , to derive the number of FPTs for the gap closing. At k 0 , the states can be labeled by the eigenvalues of U. If the gap closing between two states with the same U eigenvalues, the effective model can be described by Supplementary Eq. (27) with |C x | = 0. The gap closes if and only if C y q y + M = 0, realizable by making two vectors M and C y parallel. Such realization needs 2 FTPs, e.g. the two components of the projection of M on the plane perpendicular to C y .
When the gap closes between two states with different U eigenvalues, the effective model along the U-invariant line (q x = 0) reads h(q) = E 0 (q y ) + (m 0 + Cq y + B 0 q 2 y )σ z , which, by shifting the k 0,y , can be simplified to h(q) = E 0 (q y ) + (m + Bq 2 y )σ z . The gap for this Hamiltonian keeps closing when mB ≤ 0, indicating a stable gapless phase protected by U with 0 codimension.
In this scenario, we here only consider the (UT ) 2 = −1 case, where each band at the gap closing momentum is doubly degenerate. We can define the UT pair as the two degenerate states related by UT , in analog to the Kramers pair defined in Supplementary Note 3 A. Similar as Supplementary Eq. (29), there are 5 mass terms for the gap closing between two UT pairs. However, the case here is different from the TRIM scenario in Supplementary Note 3 A, since q x does not change under UT and thus the corresponding terms have the same form as the mass terms in Supplementary Eq. (29). One of the five mass terms can then be canceled by shifting k 0,x , resulting in 4 FTPs for the gap closing.

I. Scenario (i):TRIM
We first discuss the gap closing at Γ between two Kramers pairs of the same type. If the bases have C 3 eigenvalues (e −iπ/3 , e iπ/3 , e −iπ/3 , e iπ/3 ), the mass term of the effective model is given by Supplementary Eq. (29) with ∆ 1 = ∆ 2 = 0 since the bases with different C 3 eigenvalues cannot be coupled, resulting in 3 FTPs for the gap closing. If the bases have C 3 eigenvalues (−1, −1, −1, −1), the effective model equals to Supplementary Eq. (29) that has 5 FTPs for the gap closing. Now we discuss the construction of the effective model for the bases (e −iπ/3 , e iπ/3 , −1, −1). The form of the effective model, Supplementary Eq. (13), is given by the tensor product of the bases in the same IR listed in Supplementary Tab. 1. Note that the matrix representation and the bases for the E IR are not Hermitian. It means given two copies of (E, +) or (E, −) IR, say (τ + (σ x − iσ y ), τ + (σ x + iσ y )) and (k x − ik y , k x + ik y ) furnishing (E, −) IR, the coefficients used for the tensor product can be complex, e.g.
II. Scenario (ii): C3 ∈ G0 and T / ∈ G0 Here we consider the gap closing between two states with the same C 3 eigenvalues at K or K . In general, the mass terms at one gap closing momentum are m x σ x + m y σ y + m z σ z . Since the gap closing momentum is fixed, none of the three mass terms can be canceled by shifting the momentum, and hence there are 3 FTPs for the gap closing. When the two Kramers pairs carry C 3 eigenvalues as (e −iπ/3 , e iπ/3 , e −iπ/3 , e iπ/3 ), the effective model equals to Supplementary Eq. (29) with ∆ 1 = ∆ 2 = 0 before considering m x , similar to the correspond case in Supplementary Note 3 C. As m x= − iσ x for each Kramers pair, the ∆ 3 is also forbidden, resulting in 2 FTPs for the gap closing. On the other hand, if C 3 eigenvalues are all −1, the effective model equals to Supplementary Eq. (29) before considering m x , similar to the correspond case in Supplementary Note 3 C, and including m x makes ∆ 2 = ∆ 3 = 0, leading to 3 FTPs for the gap closing.
The construction of the effective model for the bases (e −iπ/3 , −1, e iπ/3 , −1) is the same as that for the (111) HgTe/CdTe quantum well, which is discussed in Supplementary Note 5. Next we show that the gap closing at Γ in this case cannot drive a gapped system to the mirror protected gapless phase. Since the three mirror lines are related by the C 3 symmetry, we only need to consider one of them, say k x = 0 that is invariant under m x . The eigenvalues along this line read with α, β take ±. E ±β bands cross at Γ and belong to the same set of connected bands. The mirror eigenvalue of the E αβ band is −αi, and then the mirror protected gapless phase happens when E ++ crosses with E −− or E +− crosses with E −+ . Both band crossings require the same condition since they are related by the TR symmetry. However, the above equation has no solution when m = 0 and v 2 3 +v 2 6 = 0. It can be seen from the inequality (a + b) 2 + c 2 + (a − b) 2 + c 2 > 2|b|, which holds unless c = 0 and |a| ≤ |b|. Therefore, without finely tuning more parameters to realize v 2 3 + v 2 6 = 0, a gapped system remains when the sign of m flips.
II. Scenario (iii): UT ∈ G0 and T / ∈ G0 Here we discuss the case when the gap closes at K and K for PG p3m1 and between two states with the same C 3 eigenvalues. Before considering m x T , the mass terms at K are m x σ x + m y σ y + m z σ z since C 3 does not provide any constraints and the fixed gap closing momentum cannot be shifted to cancel any of them. Since m x T can be chosen as σ 0 K, m y is forbideen and the remaining two mass terms serve as the 2 FTPs for the gap closing.

Supplementary Note 4. VCN in Tight-Binding Model
In this section, we discuss the quantization and physical meaning of the VCN change in a tight-binding (TB) model with p1m1. We consider a square lattice and each unit cell only contains one atom. Without loss of generality, we set the lattice constant to 1, and choose the mirror symmetry as m y . On each atom, we include a spinful s and a spinful p y orbitals, meaning that the bases can be labeled as |R, α, s with R the lattice vector, α = s, p y for orbital, and s =↑, ↓ for spin. The bases with specific Bloch momentum can be obtained by the following Fourier transformation Then, the representations of the symmetries read m y |k, α, s = |m y k, α , s (τ z ) α α (−iσ y ) s s T |k, α, s = | − k, α, s (iσ y ) s s , where τ and σ are Pauli matrices for orbital and spin. With on-site terms and nearest-neighbor hopping terms, we choose the following symmetry-allowed expression for the Hamiltonian where d 1 = m + 2t 1 cos(k x ) + 2t 2 cos(k y ) The eigenvalues of h(k) read ± d 2 For concreteness, we choose t 1 = t 2 = t 3 = λ 1 = 1/2 and m = 4/5 and assume the model is half-filled (two occupied bands). In this case, the gap closes only when we tune λ 2 to zero, as shown in Supplementary Fig. 4(a), and the gap closing points sit at k = (± arccos( −8+5 √ 3 10 ), ±5π/6) ≈ (±1.50, ±2.62), belonging to the VCN scenario (iv) for p1m1. When the gap is small but nonzero, the positions of valleys can be determined numerically by locating the peaks of Berry curvature ( Supplementary Fig. 4(b)), which are close to the gap closing points.
Finally, based on the TB model Supplementary Eq. (34), we discuss the quantization and physical consequence of the VCN change across the λ 2 = 0 gap closing. Without loss of generality, we take the valley in the k x,y > 0 quarter of 1BZ as an example to discuss the quantization. The VCN of this valley can be calculated by integrating the Berry curvature over the k x,y > 0 quarter of 1BZ. As shown in Supplementary Fig. 4(c), although VCN is not quantized on any side of the gap closing, the change of VCN across the gap closing is an integer, consistent with the effective-model analysis in the main text. According to Supplementary Ref. [9], one physical consequence of the quantized VCN change is the gapless domain-wall mode in a domain wall structure that consists of the two different gapped phases separated by the gap closing, like Supplementary Fig. 4(d). As shown in Supplementary Fig. 4(e), the VCN change for each valley matches the number of gapless domain-wall modes around that valley.

Supplementary Note 5. (111) HgTe/CdTe Quantum Well
In this section, we provide more details on the analysis of the HgTe QW. Before going into details, we first introduce some basic properties of the QW. Both HgTe and CdTe have the standard zinc-blende structure, similar to most II-VI or III-V compound semiconductors. The crystallographic space group of both compounds is F43m (space group No. 216). In the QW, HgTe serves as a well while Hg 1−x Cd x Te serves as the barrier. Similar to early experimental and theoretical studies 10-14 , we use x = 0.7.
A. d-induced PET jump for E = 0 To describe the TPQT, we project the 6-band Kane model onto the bases (|E 1 , + , |H 1 , + , |E 1 , − , |H 1 , − ) via second order perturbation 13 and get the following 4-band model where the values of the parameters are listed Supplementary Tab. 3, k 2 = k 2 1 + k 2 2 , k ± = k 1 ± ik 2 , and k 1 and k 2 are the momenta along (1, −1, 0) and (1, 1, −2), respectively. Compared to the celebrated Bernevig-Hughes-Zhang model 10 , we have an additional k-linear term A 2 due to the reduction of the full rotational symmetry to C 3 rotation symmetry. In the Supplementary Eq. (36), the TQPT shown in Fig. 2 (a) of the main text occurs at m = 0. To show the jump of the symmetry-allowed PET components at the gap closing, we need to introduce the electron-strain coupling h (1) ef f based on the symmetry: where u 2 = u 11 + u 22 and u ± = u 11 − u 22 ± i(u 12 + u 21 ). This electron-strain coupling is in the most general symmetryallowed form to the leading order of u ij , which definitely includes the IB terms, ξ 3 and ξ 4 . With Supplementary Eq. (36) and Supplementary Eq. (37), the independent PET component γ 222 can be derived analytically as resulting in the PET jump as .
Based on Supplementary Eq. (38) and ξ 1,2,3,4 = 1eV (comparable to those in Supplementary Ref. [5]), we plot the γ 222 of the function of the width in Fig. 2 (b) of the main text, which shows a jump around d = 65Å.
B. E-induced PET jump for fixed d After including the electric field, we project the modified Kane model onto the bases (|E 1 , + , |H 1 , + , |E 1 , − , |H 1 , − ) via second order perturbation and get the following 4-band model Compared with Supplementary Eq. (36), the above Hamiltonian has three extra IB terms D 1,2,3 brought by the electric field. In fact, it is now in the most general symmetry-allowed form up to the second order of the momentum for the HgTe/CdTe QW along the (111) direction. In addition, the parameter m (mass term) can also be controlled by electric field. In the contrast to (001) QW, the constant (k-independent) IB terms in Supplementary Ref. [11] 37), the parameter expression, and ξ 1,2,3,4 = 1eV comparable as those in Supplementary Ref. [5], the PET jump can be calculated.
From Supplementary Tab. 5, the most general symmetry-allowed Hamiltonian without the electron-strain coupling can be derived to the k 2 order, resulting in Supplementary Eq. (40). As shown in Supplementary Tab. 5, u ij behaves the same as the k 2 term, and thereby the electron-strain coupling has the same form as the k 2 term in Supplementary Eq. (40).