Topological Nonsymmorphic Crystalline Superconductors

# Topological Nonsymmorphic Crystalline Superconductors

Qing-Ze Wang Department of Physics, The Pennsylvania State University, University Park, Pennsylvania 16802-6300, USA    Chao-Xing Liu Department of Physics, The Pennsylvania State University, University Park, Pennsylvania 16802-6300, USA
July 24, 2019
###### Abstract

Topological superconductors possess a nodeless superconducting gap in the bulk and gapless zero energy modes, known as “Majorana zero modes”, at the boundary of a finite system. In this work, we introduce a new class of topological superconductors, which are protected by nonsymmorphic crystalline symmetry and thus dubbed “topological nonsymmorphic crystalline superconductors”. We construct an explicit Bogoliubov-de Gennes type of model for this superconducting phase in the D class and show how Majorana zero modes in this model are protected by glide plane symmetry. Furthermore, we generalize the classification of topological nonsymmorphic crystalline superconductors to the classes with time reversal symmetry, including the DIII and BDI classes, in two dimensions. Our theory provides a guidance to search for new topological superconducting materials with nonsymmorphic crystal structures.

###### pacs:
74.78.-w, 73.43.-f, 73.20.At, 74.20.Rp

## I Introduction

The research on topological superconductors (TSCs) has attracted intensive interests due to its gapless boundary excitations, known as the “Majorana zero modes”ryu2002 (); fu2008 (); roy2008 (); qi2009 (); wilczek2009 (); law2009 (); lutchyn2010 (); hasan2010 (); qi2011 (); alicea2012 (); leijnse2012 (); mourik2012 (); rokhinson2012 (); das2012 (); beenakker2013 (), with intrinsically non-local nature and exotic exchange statistics, and aims in the potential applications in low-decoherence quantum information processing and topological quantum computationsivanov2001 (); moore1991 (); nayak2008 (); alicea2011 (). The search for new topological superconducting phases and materials is a substantial step for this goal.

The first classification of TSCs (and also other topological insulating phases) was achieved by Schnyder et al.schnyder2008 () based on Altland-Zirnbauer symmetry class zirnbauer1996 (); altland1997 () for the systems with or without particle-hole symmetry (PHS), time reversal symmetry(TRS) and their combination, the so-called chiral symmetry. Later, it was realized that when additional symmetry exists in a system, new topological phases can be obtained, and the gapless edge/surface modes require the protection from additional symmetry. In particular, it has been shown that new topological insulating and superconducting phases emerge when the system has mirror symmetry and are dubbed “topological mirror insulators”teo2008 (); hsieh2012 (); shiozaki2014 (); wang2015 (); sun2015 () and “topological mirror superconductors” zhang2013 (); ueno2013 (); chiu2013 (); shiozaki2014 (), respectively. Recent work has also revealed that nonsymmorphic symmetry, including glide plane symmetry and screw axis symmetry, can lead to new topological insulating phases, as well as topological semi-metal phasesparameswaran2013 (); liu2014 (); fang2015 (); shiozaki2015 (); fang2015a (). In this work, we are interested in the role of the nonsymmorphic crystalline symmetry, mainly glide plane symmetry, in the classification of TSCs. We focus on the following three questions: (1) are there any topologically non-trivial phases that are protected by glide plane symmetry? (2) What’s the difference between glide plane symmetry and mirror symmetry in the classification of TSCs? (3) What’s the relationship between this superconducting phase and other TSCs? Below, we will first discuss the role of glide plane symmetry in the classification of superconducting gap functions, which indicates the possibility of topological superconductors protected by glide plane symmetry, thus dubbed “topological nonsymmorphic crystalline superconductors (TNCSc)”. We also construct an explicit tight-binding model in the D class with boundary Majorana zero modes and demonstrate that the existence of Majorana zero modes comes from nonsymmorphic symmetry of this model. Finally, we discuss the relationship between TNSCs and weak TSCs and generalize TNCSc to the classes DIII and BDI with time reversal symmetry for both spinless and spin- fermions.

## Ii nonsymmorphic symmetry and superconducting gap function

In this section, we will first consider the role of glide plane symmetry in the classification of superconducting gap functions. We start from a generic Bogoliubov-de Gennes (BdG) type of Hamiltonian of superconductors with nonsymmorphic symmetry in the normal states, which can be written in the momentum space as

 H=12∑k(c†k,cT−k)HBdG(ckc†T−k)

with

 HBdG=(h(k)−μΔ(k)Δ†(k)−h∗(−k)+μ), (1)

where is for single-particle Hamiltonian of normal states, is the chemical potential and denotes the superconducting gap function. is an annihilation operator with components and we also use () to denote each component with for spins and orbitals(lattice sites) . The superconducting gap function is related to annihilation operators by , where is the strength of attractive interactions. The BdG Hamiltonian satisfies the PHS with the PHS operator , where is the first Pauli matrix acting on the Nambu space, is an unit matrix and K is complex conjugation. The PHS (or Fermi statistics) requires the constraint for the gap function.

Next, we consider how the nonsymmorphic symmetry yields constraint on the forms of single-particle Hamiltonian and superconducting gap functions in a nonsymmorphic crystal. Here we consider the glide plane symmetry, represented by where is a mirror operator and is a non-primitive translation operator along a direction within the mirror plane. For single-particle Hamiltonian, the glide plane symmetry requires , where is the representation matrix for glide plane symmetry at the momentum and defined as appendix (). Here we emphasize that the representation matrix for glide plane symmetry takes the form , where is a phase factor due to a non-primitive translation and is the projective representation of mirror operator . For the case with only glide plane symmetry, all the projective representations are one dimensional (1D) and equivalent to the conventional representations.

The symmetry of the gap function is determined by the Cooper pair wave functions, which transform as the direct product of the representation appendix (). For the case with only glide plane symmetry, all the 1D representations can be labeled by where for spin- systems and for spinless systems. Thus, the gap function should transform as , where applies for both the spin- and spinless systems and depends on the nature of superconducting gap functionsappendix (). We will show how to classify different superconducting gap functions based on glide plane symmetry explicitly for a model Hamiltonian in the next section.

We emphasize that the superconducting gap function may preserve () or spontaneously break () glide plane symmetry. Nevertheless, similar to the case of inversion symmetryfu2010 (); yang2014 () or mirror symmetryueno2013 (), one can always re-define a glide plane symmetry operation as for the BdG type of Hamiltonian, which satisfies the condition . In this way, we can regard the BdG Hamiltonian as a semiconductor Hamiltonian with additional PHS.

Due to the existence of the glide plane symmetry , the eigenstates of the BdG Hamiltonian, , can also be chosen to be the eigenstate of , , on the glide invariant plane (GIP) in the momentum space, (mod ), where is a reciprocal lattice vector. Here is given by () for the spinless (spin-) systems and we call the eigenvalue as glide parity. Next, we look at the relationship of glide parities between one eigenstate and its partner under PHS. Direct calculation gives by using that appendix (). Therefore, and its particle-hole partner possess glide parity and , respectively. This leads to the conclusion as depicted in Fig. 1. When the gap function satisfies symmetry, for the spinless (spin-) systems, and its particle-hole partner share the same glide parity along the momentum line () on the GIP, while they have opposite glide parities along the momentum line () on the GIP. When the gap function satisfies symmetry, we find an opposite behavior for the momentum lines and , compared to the case of symmetry. We notice that the momentum line corresponds to the BZ boundary since is a primitive lattice vector of the system.

Here we emphasize different roles of glide plane symmetry and mirror symmetry for the BdG Hamiltonian of superconductivity. For the glide plane symmetry and the corresponding mirror symmetry , the GIP and the mirror invariant plane are the same. As shown in Ref. ueno2013, , the PHS either preserves the subspace with a fixed mirror parity or transforms the subspace with one mirror parity to the other. In contrast, due to the additional phase factor from the non-primitive translation of glide plane symmetry, the behaviors of PHS acting on the glide parity subspaces are always opposite for the momentum lines and . This prevents us to define a topological invariant on the whole 2D GIP since two glide parity subspaces are always “connected” to each other. However, if we limit the glide parity subspace only on the momentum line or , the PHS will either preserve the glide parity subspace or transform the subspace with one glide parity to the other, similar to the case of mirror symmetry. This immediately suggests the possibility of defining topological invariants on the 1D momentum lines or for superconductors with glide plane symmetry. Below, we will present explicitly a BdG type of model Hamiltonian with glide plane symmetry and show the existence of Majorana zero modes at the boundary. Then we will discuss bulk topological invariants and the corresponding topological classification.

## Iii Model Hamiltonian in the D class

Our spinless fermion model with glide plane symmetry is based on a two dimensional (2D) rectangle lattice with two sets of equivalent sites, as shown by A and B sites in Fig. 2(a) and (b). The glide plane symmetry operator is given by with a reflection along the z direction followed by a translation of along the x direction ( is a lattice constant), and relates the A sites to the B sites. The normal state Hamiltonian reads

 h(k)= ϵ(k)σ0+t3cos((kx−ϕ)a2)cos(kxa2)σ1 (2) +t3cos((kx−ϕ)a2)sin(kxa2)σ2

on the basis and , where , is a 22 unit matrix, with , , are Pauli matrices that describe the A and B sites and depends on the choice of orbitals appendix (). Furthermore, the glide plane symmetry operator on such a basis is . One can easily check that and .

As discussed above, the gap functions can be classified according to glide plane symmetry and when the glide plane symmetry for the BdG Hamiltonian is , the gap function satisfies three conditions: (PHS); (glide plane symmetry) and appendix (). The complete classificaiton of gap functions for this model Hamiltonian is discussed in the Supplemental Material appendix (). Here we only consider two typical gap functions and with the symmetries and , respectively. We take the BdG Hamiltonian (Eq. 1) with the single-particle Hamiltonian (Eq. 2) and the gap function and calculate energy dispersion of this Hamiltonian on a slab configuration. The slab is chosen to be infinite along the x direction and finite along the y direction, so that the glide plane symmetry is still preserved. The energy dispersion is shown in Fig. 2 (c) for and (d) for . In both cases, one can find two edge bands appearing in the bulk superconducting gap at one edge. However, these two edge bands cross at zero energy and give rise to Majorana zero modes at for (Fig. 2(c)), but at for (Fig. 2(d)).

The underlying physical reason of different positions of Majorana zero modes for these two cases comes from the relation between glide plane symmetry and PHS discussed in the last section. Let’s take the case of the gap function with symmetry as an example. The state and its particle-hole partner share the same glide parity at (), and thus it is possible for them to be the same state. Since PHS changes the energy of to of , the eigen energy must be zero once they are the same state. This analysis also suggests that two Majorana zero modes at must belong to different glide parity subspace, and thus no coupling is allowed between them to open a gap. In contrast, the glide parities for and are opposite at (). Thus, these two states must be different at and PHS can not require their energies to be zero. This analysis can also be applied to with symmetry and leads to the opposite conclusion. Another intuitive picture to prove non-trivial properties of 1D edge modes in Fig. 2 (c) and (d) is to consider a general one dimensional superconductor with glide plane symmetry. As shown in Ref. liu2014, ; fang2015, ; shiozaki2015, , due to the glide plane symmetry, all the bands must appear in pairs, as shown schematically by two black lines (two bands with opposite glide parities) in Fig. 2 (f). Furthermore, the PHS of superconductivity requires two additional hole bands at the negative energy, as shown by two red lines in Fig. 2 (f). Therefore, there must be even number of pairs of bands for a 1D nonsymmorphic superconductor. A single pair of bands shown in Fig. 2 (c) and (d) can only exist at the 1D boundary of a 2D system. This gives the “no-go” theorem for nonsymmorphic superconductorswu2006 ().

## Iv Bulk topological invariants and the extended Brillouin zone

The above analysis has shown that two momentum lines and play the essential role in the classification of TSCs in nonsymmorphic crystals. We can view the bulk Hamiltonian on or as a 1D Hamiltonian. For the case of with the symmetry, the Hamiltonian (Eq. 1) has PHS along the line for each glide parity subspace, thus belonging to the D class, while it has no PHS along the line for each glide parity subspace, as shown in Fig. 1 (a). Since two glide parity subspaces are decoupled, one topological invariant of the D class can be defined on the line in the glide parity subspace for a 1D Hamiltonian. In contrast, for the case of with the symmetry, one topological invariant can be defined on the line . In our example, we can re-write the BdG Hamiltonian with the eigenstates of as a basis and one can see immediately for the case with the () symmetry, the Hamiltonian is exactly equivalent to the 1D Kitaev model of p-wave superconductorskitaev2001 () in each glide parity subspace when ()appendix ().

More insights about this system can be obtained from the view of the extended Brillouin zone (BZ)lee2008 (), which has been widely used in the field of iron pnictide superconductors. For nonsymmorphic crystals, all the eigenstates of the Hamiltonian can be labeled by the eigenvalues of glide operators, defined as , in which is called “pseudocrystal momentum”lee2008 () and defines the extended BZ. For our model, the glide plane symmetry operation only involves translation by along the x direction, and thus the extended BZ for is doubled along the x direction (), compared to the conventional BZ for , as shown in Fig. 2 (e). Since we have , this suggests that the pseudocrystal momentum is related to momentum by when and with when . Here the sign of is determined by keeping in the region and in the region . As a result, the BdG Hamiltonian can also be rewritten as for and for and in the extended BZ. Here is the BdG Hamiltonian in the subspace with glide parity . For our model Hamiltonian, corresponds to the two by two Hamiltonian defined in the Supplemental Materialappendix (). For the case of , the form of is given by , where . We notice that if we take the hopping parameters and along the x direction to be zero, this Hamiltonian exactly corresponds to the 1D Kitaev chain with one Majorana zero mode at the open boundarykitaev2001 (). With the hopping along the x direction, all the 1D Kitaev chains are coupled along the x direction, so Majorana zero modes at the end of the chains couple to each other and expand into a band. This corresponds to the weak TSCsteo2010 (); hughes2014 (), which is in analogy to weak topological insulatorsfu2007 (). The PHS requires for the band of Majorana zero modes. Therefore, zero energy states can only appear for and ( is periodic in ), which both correspond to in the conventional BZ. In contrast, for the case of , the gap function comes from the so-called pairing for two electrons with the momenta and to form a Cooper pairyang1989 (); hu2012 (); hu2013 (); hao2014 (); wang2015g (); appendix () ( for our model). In this case, the PHS requires for the Majorana band, leading to the zero energy states at . This analysis based on the extended BZ is consistent with our previous results and show explicitly the relationship between TNSCs and weak TSCs.

## V Discussion and conclusion

The above results for TNCSc can be directly generalized to the systems with spin- and with additional time reversal (TR) symmetry. For spin- systems, since in the glide parity is given by , there is an additional minus sign when considering how the glide parity of an eigenstate of the BdG Hamiltonian transforms under PHS. This leads to the consequence that the topological invariant can be defined at () for the systems with the () symmetry. According to the standard topological classification, TR symmetry can change the symmetry class from the class to for spinless systems and for spin- systems. To see how it affects the classification of TNCSc, we consider an example of a spin- system in the class with the symmetry. If we take a state with glide parity where , the glide parity of its PHS partner has been shown to be and the glide parity of its TR partner is also , where is used and TR operator is with and the second Pauli matrix acting on spin space. One can see that chiral symmetry exists in each glide parity subspace for any momentum. In addition, at the momentum line , PHS and TRS also exist in each glide parity subspace. Therefore, the symmetry class is DIII for the momentum line and AIII for other momentum lines () in each glide parity subspace. This leads to classification at , classification at and classification at other momentum lines for the whole BdG Hamiltonianschnyder2008 (); appendix (), in sharp contrast to the classification of topological mirror superconductors in the DIII classzhang2013 (). This classification leads to the existence of edge flat bands in the DIII class for TNCSc (See Supplemental materials appendix ()). The spinless and spin- TNCSc in classes D, DIII and BDI are also studied in the Supplemental Materialappendix (). Nonsymmorphic symmetry is known to exist in several classes of superconducting materials, including iron pnictide superconductorskamihara2008 (); stewart2011 (); dagotto2013 (); cvetkovic2013 (); qing2012 (); he2013 (); tan2013 (), BiS-based layered superconductorsmizuguchi2012 (); mizuguchi2012a (); usui2012 (); yazici2013 (); jha2013 (); xing2012 (); yang2013 (); lin2013 (); dai2015 (), and heavy fermion superconductorsstewart1984 (); stewart2011 (), e.g. UPtjoynt2002 (), UBeott1983 (). Our topological classification of TNCSc can be directly applied to these systems to search for realistic topological superconducting materials.

###### Acknowledgements.
We would like to thank X. Dai, X.Y. Dong, Ken Shiozaki, Fan Zhang and Jiangping Hu for the helpful discussions.

Note added. - After finishing this paper, we notice a paper on arxivhao2015 (), which concerns possible topological superconducting phases in monolayer FeSe and potential relation to nonsymmorphic symmetry. We also notice another recent paper on arxivshiozaki20152 () about topological classification of TNCSc based on the twisted equivariant K-theory.

## Appendix A Bogoliubov-de Gennes Hamiltonian and the glide plane symmetry

We start from a generic mean-field Hamiltonian of superconductors with nonsymmorphic symmetry in the normal states, which can be written in the momentum space as

 H =∑k,α,βc†k,α(hα,β(k)−μ)ck,β+∑k,α,β12(Δ†α,βc−k,αck,β+Δα,βc†k,αc†−k,β) (3) =12∑k(c†k,cT−k)HBdG(ckc†T−k)

with

 HBdG=(h(k)−μΔ(k)Δ†(k)−h∗(−k)+μ), (4)

where is for single-particle Hamiltonian of normal states, is the chemical potential and denotes the superconducting gap function. is an annihilation operator with components and we also use () to denote each component with for spins and orbitals(lattice sites) . The superconducting gap function is related to annihilation operators by , where is the strength of attractive interactions. The Bogoliubov-de gennes (BdG) Hamiltonian satisfies the particle-hole symmetry (PHS) with the PHS operator , where is Pauli matrix acting on the Nambu space, is an unit matrix and K is complex conjugation. The PHS (or Fermi statistics) requires the constraint for the gap function.

glide plane symmetry can be expressed as with a mirror operator and a non-primitive translation operator along the direction within the mirror plane. We have and where is representation of space group in the little group of wave-vector . Further, , where is the projective representation of mirror operator .

Since the normal state has the glide plane symmetry, . Explicitly,

 ∑k,αβgc†k,αg−1hαβ(k)gck,αg−1 =∑k,α,βD∗k,αγ(g)c†gk,γhαβ(k)Dk,βλ(g)cgk,λ =∑k,γλc†k,γ(∑αβD∗g−1k,αγ(g)hαβ(g−1k)Dg−1k,βλ(g))ck,λ

Therefore, we have

 ∑αβD∗g−1k,αγ(g)hαβ(g−1k)Dg−1k,βλ(g)=hγ,λ(k) (5)

More elegantly, the normal state Hamiltonian under glide plane symmetry satisfies

 D†k(g)h(k)Dk(g)=h(gk). (6)

We also have .

The pairing terms in Eq. 3 fulfill the requirement of glide plane symmetry. However, the gap functions are not necessary to respect the glide plane symmetry.

 g(∑k,αβΔαβ(k)c†k,αc†−k,β)g−1=∑k,αβγλD∗k,αγ(g)D∗−k,βλ(g)Δαβ(k)c†gk,αc†−gk,β=∑k,γλ~Δγλ(k)c†k,γc†−k,λ

, where is the transformed gap function. Thus, we arrive at . This indicates that the gap function transforms according to the decomposition of the direct product of the representation .

If we choose the basis as the eigenstates of glide plane symmetry operator, all representations are reduced to one dimension(1D) and where for spin- systems and for spinless systems. Thus, . We can write down the requirement of gap function under nonsymmorphic symmetry in a compact way

 ~Δ(gk)=D†k(g)Δ(k)D∗−k(g)=ηΔ(gk) (7)

where applies for both spin- and spinless systems.

The matrix form of the BdG Hamiltonian is shown in Eq. 4 in the basis of Nambu space . According to Eq. 7, we can treat the BdG Hamiltonian as a semiconductor Hamiltonian with additional PHS and re-define the glide plane symmetry operator as

 Gη(k)=(Dk(g)00ηD∗−k(g)) (8)

Next we check how the BdG Hamiltonian transforms under glide plane symmetry by considering possible gap functions required in Eq. 7. For the case with , we have

 G−1η(k)HBdG(k)Gη(k) =(D−1k(g)00ηD∗,−1−k(g))(h(k)−μΔ(k)Δ(k)−h∗(−k)+μ)(Dk(g)00ηD∗−k(g)) (9) =(h(gk)−μD−1k(g)Δ(k)ηD∗−k(g)h.c.−h∗(−gk)+μ) =(h(gk)−μη2Δ(gk)h.c.−h∗(−gk)+μ) =(h(gk)−μΔ(gk)h.c.−h∗(−gk)+μ) =HBdG(gk)

where is used. This gives us the form of symmetry transformation for the BdG Hamiltonian.

Due to the glide plane symmetry, all the eigenstates of the BdG Hamiltonian at the glide invariant plane (GIP) can be also expressed as the eigenstates of glide plane symmetry and the corresponding eigenvalues are dubbed “glide parity”, as discussed in the main text. To show the glide parities of a state and its particle-hole partner, we take an example of the case with the symmetry. If the gap function preserves the glide plane symmetry, the BdG Hamiltonian commutes with on GIPs. This indicates that one can simultaneously block diagonalize and with a set of common eigenvectors. Each block owns a glide parity with for spinless(spin-) systems.

We start from that , i.e. , where C is the PHS operator. One can pick up a common eigenstate of and with glide parity . Its particle-hole partner is denoted as . Then we have . Therefore, and its PHS partner possess glide parity and , respectively.

## Appendix B Model Hamiltonian

In this section, we will show how to construct a tight-binding model for the TNCSc. We use to denote Lwding orbital at , where is the position of unit cell and is the position of atom site or orbital inside a unit cell. The Bloch wave function is defined as by using the linear combination of atomic orbitals(LCAO). It should be emphasized that the phase factor in our construction does not include the position . This Bloch wave function can be written as where . It is easily checked that . Thus, the Bloch theorem holds for such a choice of LCAO. Under such a choice of LCAO, we have that , where is a reciprocal lattice vector. Any Hamiltonian on such a basis satisfies that , i.e.

 H(k+P)=H(k). (10)

Similarly, for the gap function, , i.e.

 Δ(k+P)=Δ(k) (11)

.

Next let us take and ( and ) to present creation and annihilation operators of () and show the real space form of the tight-binding model discussed in the section “Model Hamiltonian” of the main text. The normal state tight-binding Hamiltonian for the lattice structure reads

 He= m0∑i,s=A,Bc†i,sci,s+[∑i,s=A,B(t12c†i+dx,sci,s+t22c†i+dy,sci,s)+∑it32(c†i,Aci,B+c†i,Aci−dx,B)+H.c.] (12)

where denotes index of unit cells, are primitive lattice vectors along x and y direction, A and B denote two inequivalent atom sites and H.c. represents their conjugation parts. We further obtain a tight-binding model in the momentum space by performing an unusual Fourier transformation and , where is the position of the i unit cell. Such a Fourier transformation simplifies the Hamiltonian and leads to with is a reciprocal lattice vector.

The gap functions for configurations mentioned previously need to satisfy three conditions: (1) PHS

 ΔT(k)=−Δ(−k); (13)

(2) glide plane symmetry

 D†k(g)Δ(k)D∗−k(g)=ηΔ(gk); (14)

and (3)

 Δ(k)=Δ(k+G). (15)

These three conditions allow us to classify all the possible gap function for this model Hamiltonian. Let us define the gap functions , where , is a complex function of and are four 22 matrices, as shown in the second column in Table LABEL:pair.

The second and third column in Table LABEL:pair shows the “parity” of matrices in the sense of PHS and glide plane symmetry. Due to the PHS (13), the parity of is determined by the parity of from the second column, which is listed in the fourth column. For 2D system, if . For the last two columns, the existence of is determined by Eq. (14) and in the third column. Let us take the matrix with as an example. For , since glide operation does not act on , we obtain and there is no constraint on . But for , Eq. (14) and together requires , leading to . Thus, no term is possible for . Similar analysis can be applied to other matrices. Based on the parity of on the fourth column, we can get possible polynomials of , as listed in Table LABEL:pairf.

The parameters for the calculation of energy dispersion of the BdG Hamiltonian are shown in Table LABEL:par.

## Appendix C Hamiltonian in the extended Brillouin zone for G±

In this section, we will analyze our model Hamiltonian in the extended Brillouin zone. For the case with the symmetry, . The eigenvalues are . The eigenvectors are , , and