Suppressed topological phase transitions due to nonsymmorphism in SnTe stacking

We combine first principles calculations with a group theory analysis to investigate topological phase transitions in the stacking of SnTe monolayers. We show that distinct finite stacking yields different symmetry-imposed degeneracy, which dictates the hybridization properties of opposite surface states. For SnTe aligned along the [001] direction, an (even) odd number of monolayers yields a (non)symmorphic space group. For the symmorphic case, the hybridization of surface states lead to band inversions and topological phase transitions as the sample height is reduced. In contrast, for a nonsymmorphic stacking, an extra degeneracy is guaranteed by symmetry, thus avoiding the hybridization and topological phase transitions, even in the limit of a few monolayers. Our group theory analysis provide a clear picture for this phenomenology and matches well the first principles calculations.

The recently discovered topological insulators (TIs) are classified accordingly to the symmetries of its crystal lattice and/or time-reversal symmetry, which yields its topologically protected edge or surface states 1,2 . Particularly, the class of topological crystalline insulators (TCIs) have Dirac-like bands protected by space group symmetries 3,4 . The first material predicted to be a TCI was SnTe 5 , which was promptly verified experimentally 6 , and followed by other IV-VI compounds [7][8][9][10] . The monolayers of IV-VI materials have been predicted to be two-dimensional (2D) TCIs 11,12 . Both the bulk and monolayers of SnTe are classified by their mirror Chern number |n M | = 2, yielding an even number of Dirac cones in the Brillouin zone. In between these two limits, the stacking of monolayers show an intriguing nonmonotonic evolution of the band gap for stackings oriented along the [001] 12-15 , [111] 16-20 and [110] 21,22 directions. Particularly, odd and even stackings of [001] SnTe belong to distinct space groups: the first is symmorphic, while the latter is nonsymmorphic; i.e. some symmetry operations require translations by a fraction of the unit cell. It is known that nonsymmorphic groups yield an extra degeneracy with respect to its symmorphic counterpart 23,24 . Moreover, from the nonsymmorphic space groups arise newly proposed exotic topological properties [25][26][27][28][29][30][31][32] .
In this paper the distinct topological phases of odd and even stackings of SnTe monolayers are explained by the contrast between their symmorphic and nonsymmorphic crystal symmetries. First, we analyze first principles [density functional theory (DFT) [33][34][35][36] , see Methods] band structures as a function of the number N of stacked layers. Interestingly, we find that the symmorphic stacking presents two topological phase transitions, while no phase transition is seen for the nonsymmorphic stackings. To understand this disparity, we develop low energy models for the surface states of both symmorphic (odd N) and nonsymmorphic (even N) stackings. We show that, due to finite-size effects, topological surface states located at opposite surfaces hybridize differently in the symmorphic and nonsymmorphic cases. Particularly, the nonsymmorphism of the even stackings enforces an extra degeneracy, thus avoiding a band inversion. Ultimately, this protects the system from a topological phase transition. Our results provide a clear picture of finite size effects, and its respective space group, on the hybridization of surface states. Although we focus here on the [001] stacking of SnTe monolayers, equivalent considerations shall hold for other directions, and as well for any material that may present distinct space groups for different stacking heights.

Results and Discussions
DFT analysis. Bulk SnTe has a rock-salt lattice and shows a small band gap at the L point. There, the conduction and valence bands are inverted due to the spin-orbit (SO) interaction, as shown in Fig. 1(a,b), yielding its 3D TCI regime characterized by the mirror Chern number |n M | = 2. Around the L point, the effective Hamiltonian 37 obeys the P31m (or D ) 3d 1 symmorphic space group. In the opposite limit, a single monolayer of SnTe is a 2D TCI 4,5,24,37 . In this case, without SO the bands are already inverted, but it is metallic around the X point of the 2D Brillouin zone, Fig. 1(c). The SO coupling opens the gap shown in Fig. 1(d), yielding the topological regime. The band structure around X obeys the symmorphic space group Pmmm (or D ) 2h 1 24 . The spinful double group of both bulk and monolayer admits doubly degenerate irreducible representations (IRREPs) with odd and even inversion characters 38 , which can be characterized by the orbital projections Sn , respectively 39 . These projections label the color code in Fig. 1(a-d) and illustrates the band inversions.  Fig. 1(f), finite-size effects arise distinctively for odd or even N. In both cases, the stacking projects a pair of distinct L points into X, Fig. 1(e). For large N, the distinction between odd or even number of layers is small. Each surface shows a pair of Dirac bands, which are highlighted by the colors in Fig. 2(a,b). However, as N is reduced, topological states localized at opposite surfaces hybridize differently for the odd/even N (symmorphic/nonsymmorphic). The hybridization becomes relevant for stackings  λ 4 , where λ ∼ v E / F g  is the penetration length of the surface states. From Fig. 1(b) we estimate the bulk gap E g ≈ 0.11 eV and Fermi velocity ℏv F ≈ 2 eVÅ, which gives λ ≈ 18 Å or N ≈ 6 layers. Therefore, the hybridization becomes relevant for  N 24, which matches the gap evolution in Fig. 3(a,b).
Next we first present these gaps and the inversion eigenvalues obtained from the DFT simulations for the symmorphic and nonsymmorphic cases for various N. Later on we introduce a group theory model for the surface states, which allow us to understand how the states from opposite surfaces hybridize as N is reduced. We find that the nonsymmorphic symmetry of even N stackings imposes an extra degeneracy at X, which protects the band crossings for all even N; e.g. see arrows in Fig. 2(b) for N = 28. In contrast, for odd N these bands hybridize as N is reduced; e.g. see arrows in Fig. 2(a) for N = 27. Within the model we will see that these distinct behaviors explain the phase transitions seen for the symmorphic case, and its absence for the nonsymmorphic stacking.
Symmorphic stacking. The ab-initio band structure for odd stacking ranging from N = 5 to 39 layers is shown in Fig. 2(a). For large N, the Dirac crossings that characterizes the TCI regime are noticeable (highlighted in blue). As N is reduced, the hybridization between the top and bottom surface states take place and the bands become massive. To analyze this process, we plot the total gap ΔE and the local gap ΔE X at X in Fig. 3(a,b). For the symmorphic case (solid green circles) two phase transitions are striking. That is, ΔE X changes sign from N = 1 to 3, and from N = 11 to 13, which is characterized by the inversion eigenvalue in Fig. 3(c). For N = 1 the monolayer is a 2D TCI, in agreement with previous results 4,5,24,37 .
To understand the transition from N = 11 to 13, the density of states for the surface bands around the Dirac crossing off the X point are projected along the [001] layers, as shown in Fig. 4. Interestingly, in all cases (symmorphic and nonsymmorphic for all N) it shows that the densities are localized at the top and bottom surfaces, within the length λ (≈6 layers) as discussed previously. Each surface state belongs to an even or odd inversion IRREP of the Pmmm (or D ) 2h 1 space group 38 . These surface states do not change their character when passing from N <11 to> 13. In fact, the sole effect of reducing N is to increase the overlap and hybridization between top and bottom states. Therefore, for any  N 6 the staking maintains its 3D TCI character defined by its bulk invariants. Consequently, the band inversions observed for odd stacking shown in Fig. 3(c) must be of a lower order, i.e. a 2D topological phase transition. An illustrative model 40   Information. This 2D topological phase transition, can be better characterized by the mirror Chern number computed using the surface state bands. Using the effective model described below, we calculate this mirror Chern number n M as a function of the gap sign (ΔE X ). We find |n M | = 2 for ΔE X < 0 and n M = 0 for ΔE X > 0. This is shown by the color code in Fig. 3(c). For the monolayer case (N = 1) it is known that |n M | = 2 12,24 .   show smooth monotonic and asymptotic behaviors in Fig. 3(a,b) (open pink circles). There is no phase transition. The band structures in Fig. 2(b) show clear hybridized Dirac bands for N > 6. Notice that the crossings at X remain closed for all even N, in contrast with the hybridization seen in Fig. 2 Fig. 2(a), and N = 40 in Fig. 2(b)], yielding gapless Dirac crossings. As we will see later on, the lack of phase transitions in this case is due to the degeneracies imposed by the nonsymmorphic symmetry. Its space group is the Pmmn (or D ) 2h 13 , and due to the Bloch phase at X, the only allowed double group IRREPs are degenerate due to time-reversal symmetry 24,38 . This yields a fourfold degeneracy with pairs of states with even and odd inversion eigenvalues. This is illustrated in Fig. 3(d). Moreover, from the PDOS shown in Fig. 4 we observe that the 3D topological character is always present for nonsymmorphic stacking, independent on the number of layers. Since there is no phase transition, its topological nature is governed by the precursor bulk 3D TCI with mirror Chern number |n M | = 2.
Effective model. Next, we develop low energy models for symmorphic and nonsymmorphic stackings, which will help us to understand the phase transitions presented above. Indeed we find that only two transitions should be expected for the symmorphic case, in agreement with ab-initio results.
Let us start by defining a basis as the surface states from a semi-infinite stacking limit. Then, we allow it to hybridize with the solutions from the opposite surface due to finite symmorphic or nonsymmorphic stackings. In all cases we consider the wave-functions around X to be composed by s, p x and p z atomic orbitals of Sn and Te, as obtained from the DFT first principles calculations.
Semi-infinite stacking. Around X, the semi-infinite stacking obeys the symmorphic Pmm2 (or C 2v ) space group. Starting with the spinless case, we find by inspection that the surface states belong to the A 1 and B 1 IRREPs of C 2v . Therefore we label the spinful eigenkets for our basis as {A 1 ; ↑, A 1 ; ↓, B 1 ; ↑, B 1 ; ↓}, where {↑, ↓} refers to the spin along z. The C 2v group is composed by a rotation C 2 (z) = iτ 3 ⊗σ z around the z axis, and mirrors M x = −iτ 3 ⊗ σ x and M y = −iτ 0 ⊗ σ y that reflect x → − x, and y → − y, respectively. Here the τ i (i = 0, 1, 2, 3) are su(2) matrices acting on the A_1/B_1 orbitals and σ j (j = 0, x, y, z) are the su(2) spin matrices. The time-reversal symmetry operator is where  is the complex conjugation operator. Up to linear order in k x and k y the effective Hamiltonian that commutes with these symmetries is where k ± = (k x , ±ik y ). The coefficients α = (α x , α y ) and β = (β x , β y ) are in-plane (σ x and σ y ) SOCs that define the anisotropic Fermi velocity near k = 0 for the Dirac cones at energies = ± Δ + Δ E e 0 0 2 1 2 , while γ = (γ x , γ y ) couples orbitals A 1 to B 1 with the same spin at finite k.
We fit this Hamiltonian to the DFT data for a large stacking (N = 59 or 60), which gives the parameters: e 0 = 0.0065 eV, Δ 0 = 0.043 eV, Δ 1 ≈ 0, α x = 1.6 eV, α y = 3.3 eV, β x = 1.3 eV, β y = 2.0 eV, γ x = 0.4, γ y = 1.5 eV. The resulting band structure is shown in Fig. 5. For small k near X the agreement is patent, while for large k parabolic terms would be necessary to improve it. Particularly, along the − X M direction one branch enters the conduction band quickly and can be disregarded, as we do not account for its coupling to the bulk in the effective model. More importantly, the spin projections (color coded) match well the DFT results. Along Γ − X, both the mirror M y and the spin σ y are preserved, allowing us to label the states with its eigenvalues M y = ±i, which match the σ y projections since [M y ,σ y ] = 0. Along − X M the mirror M x is preserved, while [H ∞ , σ x ] ≠ 0, thus we label the states with the eigenvalues of M x = ±i, while the spin σ x projections are mixed (color code) at the anti-crossing induced by the γ y term.
Finite stackings and phase transitions. The space groups of the symmorphic (D 2h 1 ) and nonsymmorphic (D 2h 13 ) stackings share the C 2v as a common subgroup. Therefore, to model the finite stackings, we consider that the top (t) and bottom (b) surface states are approximately given by the semi-infinite solutions above. That is, our basis is now labeled by . Introducing s j (j = 0, 1, 2, 3) as su(2) matrices acting on the top/bottom subspace, the C 2v operations are now For the odd symmorphic stacking, the D 2h 1 also contains the operations that connect the top and bottom surfaces: inversion I, mirror M z , and rotations C 2 (x) and C 2 (y). On the other hand, for an even number of monolayers, the D 2h 13 space group is given by the C 2v operations, plus the nonsymmorphic counterparts of I, M z , C 2 (x) and C 2 (y). Using Seitz notation, these are where →  is a nonprimitive translation by half a unit cell [see Fig. 1(f)]. Since M z = I ⋅ C 2 (z), C 2 (x) = I ⋅ M x , and C 2 (y) = I ⋅ M y (and equivalent expressions for the nonsymmorphic operations), it is sufficient to define how the inversion operator acts on our basis in both cases. For the odd stacking, the inversion transforms DFT results, it is sufficient to analyze the Γ − X direction (i.e. k y = 0). Along this direction M y is a good quantum number, therefore we can block diagonalize the Hamiltonian into M y = ±i subspaces. For M y = + i, the dominant terms are ( ) ] cases. The effective Hamiltonians for M y = −i are obtained simply replacing α x → − α x and β x → − β x . In the semi-infinite limit, α x , β x and Δ 0 are the same as in H ∞ , but deviations should be expected for a small N. The new terms δ j (j = 0, 1, 2, 3) represent the hybridization between top and bottom surface states, which couple the linear dispersions due to the finite size of the system. As the number of layers N is reduced, these hybridizations induce anticrossings δ j , which increase as indicated by the arrows in Fig. 6(a,b). Other terms allowed by symmetry only contribute to the fine tuning of the spectrum (see Supplementary Information). Three-dimensional visualizations of spectrums E(k x , k y ) are shown in Fig. 6(c) using the parameters from the semi-infinite model above, and δ 0 = δ 1 = δ 2 = 5 meV for the symmorphic stacking, and δ 3 = 5 meV for the nonsymmorphic case. A direct comparison show that the only qualitative difference between the dispersions are the hybridization at X. For δ 1 ≈ δ 0 , the gaps at X [see inset of Fig. 3(a)] are for the symmorphic and nonsymmorphic cases, respectively. For the nonsymmorphic stacking, Eq. (5) tell us that the gap at X is always Δ < E 0 X NS ( ) for any intensity of the finite-size hybridization δ 3 . Consequently, there cannot be an inversion of its surface state bands. Hence, there is no phase transition for the even stacking. This is a direct consequence of the fourfold degeneracy imposed by the nonsymmorphic symmetry, which does not allow the band crossing at X to hybridize in Fig. 6(b). On the other hand, for the symmorphic stacking, the gap at X, given by Eq. (4), changes sign with increasing δ 0 ≈ δ 1 , with the 2 , which corresponds to the N = 11 transition in Fig. 3(b). In both cases, we diagonalize the inversion operator within a degenerate subspace. The resulting eigenvalues match the DFT data in Fig. 3(c,d). Moreover, this hybridization picture agrees well with the DFT band structures from Fig. 2, and the gap evolutions shown in Fig. 3(a,b) for N ≥ 3.
The transition from N = 3 to 1 is more subtle. Indeed for a single monolayer (N = 1) the definition of top and bottom surfaces become meaningless. In order to understand this transition, first let us recall that the penetration length of the surface states throughout the bulk is λ ∼ 18 Å or ∼ N 6 layers. Therefore, for  N 6 we can safely assume that states from opposite surfaces hybridize, but does not reach all the way to the opposite surface, as seen in Fig. 4. This validates the model described above and the hybridization picture from Fig. 6. However, for ∼ N 6 the surface states lengths are larger than the stacking itself, which now goes beyond our hybridization model above. Instead, we find that as N approaches the monolayer limit, a structural phase transition occurs. This can be characterized by the lattice parameters a N . For the monolayer we obtain a 1 = 6.19 Å, while for N = 2 it jumps to a 2 = 6.31 Å, and for N = 3, a 3 = 6.33 Å. For larger N it smoothly evolves to the bulk value a ∞ = 6.40 Å. This sharp transition from N = 3 to N = 1 enhances the atomic orbital hybridizations and atomic spin-orbit couplings, which dictates this phase transition.

Conclusions
In conclusion, we have shown that the band structure of the symmorphic and nonsymmorphic stackings evolve differently as a function of N. This illustrates how the specific space group of a finite size sample is essential to understand its properties within the nanometer scale. Particularly, we show that the nonsymmorphic symmetry imposes an extra protected degeneracy that avoids an hybridization at the X point, thus yielding a smooth evolution of the band structure with the number of layers. Consequently, there is no phase transition in this case, and its topological nature is governed by the precursor bulk 3D TCI regime, even for samples consisted of only a few layers. In contrast, the hybridization of the symmorphic stacking yields two 2D phase transitions. One is due to a band inversion induced by the hybridization, and the other occurs as the system approaches the monolayer limit. Our space-group-dependent hybridization picture provides a clear understanding of these two cases and matches well the DFT data. Moreover, this analysis can be easily extended to other finite-sized materials that may present a variety of distinct space groups.

Methods
Density functional theory. The first principles calculations are performed using the density functional theory (DFT) within the generalized gradient approximation (GGA) for the exchange and correlation functional 35 . Fully relativistic j-dependent pseudopotential within the projector augmented wave method 36 has been used in the non-colinear spin-DFT formalism self-consistently. We use the Vienna Ab initio Simulation Package (VASP) 33,34 , with plane wave basis set with a cut-off energy of 400 eV. The Brillouin zone is sampled by using a number of k-points such that the total energy converges within the meV scale. All stackings have been fully relaxed. The optimized force criteria for convergence was less than 0.01 eV/Å. A lattice parameter of 6.19 Å was obtained for the monolayer, and 6.31 Å for the bilayer. For larger stackings, the lattice parameter varies smoothly as a function of the number of layers, and converges to the bulk value (6.40 Å) for 20 layers.
Group theory analysis. The effective models were built using the method of invariants 23,24,41 . In each case we consider the little group  that keeps the TRIM of interest invariant. To define the relevant irreducible representations (IRREPs) for our basis, we consider the orbitals that contribute to band structure around the Fermi level, as obtained in the DFT simulation. Matrix representations for each IRREP can be found by inspection. The spin degree of freedom is then included via a direct product with its su(2) operators to form the double group 38 . The character tables for each double group are shown in the Supplementary Information.