Fragile phonon topology on the honeycomb lattice with time-reversal symmetry

TL;DR

This paper investigates the topological properties of phonons in honeycomb lattice crystals with time-reversal symmetry, revealing four fragile topological phases and their physical manifestations such as edge and corner modes.

Contribution

It applies topological quantum chemistry to phonons on the honeycomb lattice, identifying four fragile topological phases and connecting them to realistic models like graphene.

Findings

Four fragile topological phonon phases identified

Edge and corner phonon modes observed in models

Graphene is close to a topologically nontrivial phonon phase

Abstract

We use the methods of topological quantum chemistry to explore the topology of phonons on time-reversal symmetric crystals with the structure of the planar honeycomb (layer group $p6/mmm$). This approach is not tied to a particular model of atomic vibrations, but is applied to the most general dynamical matrix constrained only by the symmetries of the system. We show that four distinct fragile topological phonon phases are generically possible. Truncating the dynamical matrix to third nearest neighbors yields a model that realizes the different phonon topologies, characterized by the existence of phononic edge and corner modes and by Wilson loops with winding numbers one and two. Fitting the dynamical matrix to the DFT phonon bands shows that graphene is not very far from a topologically nontrivial phonon phase.

Tables

Table 1. Table 1: Irreps at the three high symmetry points in the BZ. All irreps are 1-dimensional, except for Γ 5 subscript Γ 5 \Gamma_{5} , Γ 6 subscript Γ 6 \Gamma_{6} and K 3 subscript 𝐾 3 K_{3} , which are 2-dimensional. The little point groups at the Γ Γ \Gamma , K 𝐾 K and M 𝑀 M points are given by 6 ​ m ​ m ​ ( C 6 ​ v ) 6 𝑚 𝑚 subscript 𝐶 6 𝑣 6mm(C_{6v}) , 3 ​ m ​ ( C 3 ​ v ) 3 𝑚 subscript 𝐶 3 𝑣 3m(C_{3v}) and C 2 ​ v ​ ( 2 ​ m ) subscript 𝐶 2 𝑣 2 𝑚 C_{2v}(2m) respectively

Phase	Acoustic branch	Optical branch
Ia	${\Gamma }_6;K_1+K_2;M_1+M_2$	${\Gamma }_5;K_3;M_3+M_4$
Ib	${\Gamma }_6;K_3;M_1+M_2$	${\Gamma }_5;K_1+K_2;M_3+M_4$
IIa	${\Gamma }_6;K_1+K_2;M_3+M_4$	${\Gamma }_5;K_3;M_1+M_2$
IIb	${\Gamma }_6;K_3;M_3+M_4$	${\Gamma }_5;K_1+K_2;M_1+M_2$

Table 2. Table 2: Combinations of EBRs that reproduce the irreps at the high symmetry points for the eight branches in Table  1 . Note that these combinations are in general non-unique.

Phase	Acoustic branch	Optical branch
Ia	$A_1@1a-A_1@2b+B_2@3c$	$A_1@3c-A_1@1a$
Ib	$B_1@3c-B_2@1a$	$B_1@1a+E@2b-B_2@3c$
IIa	$E@2b-E_2@1a$	$E_2@1a$
IIb	$E_1@1a$	$E@2b-E_1@1a$

Table 3. Table 3: δ ​ ω 2 ​ ( 𝐤 , D a ) 𝛿 superscript 𝜔 2 𝐤 subscript 𝐷 𝑎 \delta\omega^{2}(\mathbf{k},D_{a}) for the little group irreps. d 𝑑 d is the dimension of the irrep. The last column indicates the C 2 ​ z subscript 𝐶 2 𝑧 C_{2z} -eigenvalues of the normal modes.

$G_{𝐤}$	Irrep	d	$\delta {\omega }^2$	$C_{2z}$
$C_{6v}$	${\Gamma }_5$	$2$	$-6\left(a_1+a_3\right)$	$+1$
	${\Gamma }_6$	$2$	$0$	$-1$
$C_{3v}$	$K_1$	1	$3\sqrt{3}{d}_2+3\left(b_1+b_3\right)$	-
	$K_2$	1	$3\sqrt{3}{d}_2-3\left(b_1+b_3\right)$	-
	$K_3$	2	$-3\sqrt{3}{d}_2$	-
$C_{2v}$	$M_1$	1	$a_1-3a_3+2b_1-4b_2$	$+1$
	$M_2$	1	$a_1-3a_3-2b_1+4b_2$	$+1$
	$M_3$	1	$-a_1+3a_3-2b_1-4b_2$	$-1$
	$M_4$	1	$-a_1+3a_3+2b_1+4b_2$	$-1$

Table 4. Table 4: Topologically nontrivial phases.

Phase	Conditions
Ia
Ib	$d_2>0$
IIa		$a_1>3a_3$
IIb	$d_2>0$	$a_1>3a_3$

Table 5. Table 5: Phonon energies in m ​ e ​ V 𝑚 𝑒 𝑉 meV at the high symmetry points.

	DFT	3NN
$\omega \left({\Gamma }_5\right)$	$192$	$192$
$\omega \left(K_1\right)$	$122$	$124$
$\omega \left(K_2\right)$	$161$	$163$
$\omega \left(K_3\right)$	$149$	$151$
$\omega \left(M_1\right)$	$163$	$161$
$\omega \left(M_2\right)$	$171$	$169$
$\omega \left(M_3\right)$	$77$	$72$
$\omega \left(M_4\right)$	$165$	$163$

Table 6. Table 6: Non-vanishing real coefficients parametrizing the harmonic potential between n ​ t ​ h 𝑛 𝑡 ℎ nth -nearest neighbors for n = 1 , … , 7 𝑛 1 … 7 n=1,\ldots,7 .

1	2	3	4	5	6	7
$a_1$	$a_2$	$a_3$	$a_4$	$a_5$	$a_6$	$a_7$
$b_1$	$b_2$	$b_3$	$b_4$	$b_5$	$b_6$	$b_7$
$-$	$-$	$-$	$c_4$	$-$	$-$	$c_7$
$-$	$d_2$	$-$	$d_4$	$-$	$d_6$	$d_7$

Table 7. Table 7: Matrices for the relevant irreps at the Γ Γ \Gamma , M 𝑀 M and K 𝐾 K points used to compute the topological index  Q c ​ o ​ r ​ n ​ e ​ r ( 6 ) subscript superscript 𝑄 6 𝑐 𝑜 𝑟 𝑛 𝑒 𝑟 Q^{(6)}_{corner} .

$C_{2z}$	${\Gamma }_5$	${\Gamma }_6$	$M_1$	$M_2$	$M_3$	$M_4$	$C_{3z}$	${\Gamma }_5$	${\Gamma }_6$	$K_1$	$K_2$	$K_3$
			1	1	-1	-1				1	1

Table 8. Table 8: Values of Q c ​ o ​ r ​ n ​ e ​ r ( 6 ) subscript superscript 𝑄 6 𝑐 𝑜 𝑟 𝑛 𝑒 𝑟 Q^{(6)}_{corner} for the different isolated branches.

Phase	Acoustic branch	Optical branch
Ia	$5/6$	$1/2$
Ib	$1/2$	$5/6$
IIa	$1/3$	$0$
IIb	$0$	$1/3$

Equations

$$U=\frac{1}{2}\sum _{i,j}{\mathbf{r}}_i^t{\mathrm{U}}_{ij}{\mathbf{r}}_j,$$

$${\mathbf{e}}_1=\frac{1}{2}\hat{\mathbf{x}}+\frac{\sqrt{3}}{2}\hat{\mathbf{y}},{\mathbf{e}}_2=-\frac{1}{2}\hat{\mathbf{x}}+\frac{\sqrt{3}}{2}\hat{\mathbf{y}},$$

$${\mathbf{\delta }}_A=\frac{2}{3}{\mathbf{e}}_1+\frac{2}{3}{\mathbf{e}}_2,{\mathbf{\delta }}_B=\frac{1}{3}{\mathbf{e}}_1+\frac{1}{3}{\mathbf{e}}_2.$$

$$\xi =\frac{1}{\sqrt{2}}(x+iy),\bar{\xi }=\frac{1}{\sqrt{2}}(x-iy).$$

$$U_{ij}=a_n{\text{openone}}_s+b_n{s}_x+c_n{s}_y+i{d}_n{s}_z,$$

$$D_{ab}(\mathbf{k})=\sum _{\mathbf{R}}{U}_{a,0;b,\mathbf{R}}{e}^{i\mathbf{k}\cdot \mathbf{R}}.$$

$$v(\mathbf{k}{)}^t=({\xi }_A(\mathbf{k}),{\bar{\xi }}_A(\mathbf{k}),{\xi }_B(\mathbf{k}),{\bar{\xi }}_B(\mathbf{k}))$$

$$D(\mathbf{k})=$$

$$+b_1[\left(1-\cos \frac{k_x}{2}\cos \frac{\sqrt{3}}{2}{k}_y\right){\sigma }_x\otimes s_x+\left(\cos \frac{k_x}{2}\sin \frac{\sqrt{3}}{2}{k}_y\right){\sigma }_y\otimes s_x$$

$$-\sqrt{3}\left(\sin \frac{k_x}{2}\cos \frac{\sqrt{3}}{2}{k}_y\right){\sigma }_y\otimes s_y-\sqrt{3}\left(\sin \frac{k_x}{2}\sin \frac{\sqrt{3}}{2}{k}_y\right){\sigma }_x\otimes s_y]$$

$$+a_2\left(2\cos k_x+4\cos \frac{k_x}{2}\cos \frac{\sqrt{3}}{2}{k}_y\right){\text{openone}}_4+d_2\left(4\sin \frac{k_x}{2}\cos \frac{\sqrt{3}}{2}{k}_y-2\sin k_x\right){\sigma }_z\otimes s_z$$

$$+b_2\left[2\left(\cos k_x-\cos \frac{k_x}{2}\cos \frac{\sqrt{3}}{2}{k}_y\right){\text{openone}}_{\sigma }\otimes s_x+2\sqrt{3}\left(\sin \frac{k_x}{2}\sin \frac{\sqrt{3}}{2}{k}_y\right){\text{openone}}_{\sigma }\otimes s_y\right]$$

$$+a_3\left[\left(2\cos k_x+\cos \sqrt{3}{k}_y\right){\sigma }_x\otimes {\text{openone}}_s-(\sin \sqrt{3}{k}_y){\sigma }_y\otimes {\text{openone}}_s\right]$$

$$+b_3\left[\left(-\cos k_x+\cos \sqrt{3}{k}_y\right){\sigma }_x\otimes s_x-(\sin \sqrt{3}{k}_y){\sigma }_y\otimes s_x+\sqrt{3}(\sin k_x){\sigma }_y\otimes s_y\right].$$

$${\omega }^2(\mathbf{k},D_a)={\omega }^2(\mathbf{k})+\delta {\omega }^2(\mathbf{k},D_a),$$

$${\omega }^2(\Gamma )$$

$${\omega }^2(K)$$

$${\omega }^2(M)$$

$$2\sqrt{3}|d_2|>|b_1+b_3|.$$

$$|2b_1+4b_2|+|2b_1-4b_2|<2|a_1-3a_3|.$$

$$|a_1-3a_3|>\mathrm{Max}\{2|b_1|,4|b_2|\}.$$

$$\begin{array}[]{clccl}a_{1}&=-6.03\cdot 10^{3}&&b_{2}&=1.16\cdot 10^{3}\\ a_{2}&=-0.387\cdot 10^{3}&&b_{3}&=1.12\cdot 10^{3}\\ a_{3}&=-0.114\cdot 10^{3}&&d_{2}&=-0.171\cdot 10^{3}\\ b_{1}&=-2.99\cdot 10^{3}.&&&\\ \end{array}$$

$$2\sqrt{3}|d_2|\simeq 0.31|b_1+b_3|,$$

$$|a_1-3a_3|$$

$$W(k_2)\equiv P{e}^{i{\int }_0^{2\pi }d{k}_1{A}_1(k_1,k_2)},$$

$$U=\frac{1}{2}\sum _{i,j}{\mathbf{r}}_i^t{\mathrm{U}}_{ij}{\mathbf{r}}_j,$$

$$U_{ij}=\left(\begin{array}[]{cc}a_{n}+b_{n}&-c_{n}-d_{n}\\ -c_{n}+d_{n}&a_{n}-b_{n}\\ \end{array}\right),$$

$$U_{i^{\prime }{j}^{\prime }}=V(g)U_{ij}V(g{)}^{-1}.$$

$$\theta ={\text{openone}}_{\sigma }\otimes s_x,$$

$$\theta D(\mathbf{k}{)}^{\ast }{\theta }^{†}=D(-\mathbf{k}),$$

$${\tilde{D}}_{ab}(\mathbf{k})=\sum _{\mathbf{R}}{U}_{a,0;b,\mathbf{R}}{e}^{i\mathbf{k}\cdot (\mathbf{R}-{\mathbf{\delta }}_b-{\mathbf{\delta }}_a)},$$

$$\tilde{D}(\mathbf{k})=V^{†}(\mathbf{k})D(\mathbf{k})V(\mathbf{k}),$$

$$V(\mathbf{k})=\left(\begin{array}[]{cc}e^{i\mathbf{k}\cdot\boldsymbol{\delta}_{A}}\openone_{s}&0\\ 0&e^{i\mathbf{k}\cdot\boldsymbol{\delta}_{B}}\openone_{s}\\ \end{array}\right)\,.$$

Full text

I Introduction

The honeycomb lattice has played an important role in our understanding of topological states of matter. The Haldane Haldane (1988) and Kane-Mele Kane and Mele (2005a, b) models have clarified the effects of time-reversal symmetry, or its absence, on the topology of electron hamiltonians. For phonons, rotating honeycomb lattices Wang et al. (2015a); Kariyado and Hatsugai (2015); Liu et al. (2017a) where the Coriolis force plays the role of a magnetic field and gyroscopic phononic crystals Wang et al. (2015b) have been used to study the existence of nontrivial phonon topology in the absence of time-reversal symmetry. More recently, topological phonons have been reported on the valley-mixing Kekulé deformation of the honeycomb lattice Liu et al. (2017b).

On the other hand, new insights about band theory based on symmetry indicators Po et al. (2017); Khalaf et al. (2018), band combinatorics Kruthoff et al. (2017) and Topological Quantum Chemistry (TQC) Bradlyn et al. (2017) have been recently used to classify all the nontrivial electron band topologies compatible with a given crystal structure. In particular, the methods of TQC have been used to uncover thousands of new materials with topological electron bands Vergniory et al. (2018). But these new insights have not, to the best of our knowledge, been applied to the search for nontrivial topology in phonon spectra Süsstrunk and Huber (2015, 2016); Huber (2016); Kane and Lubensky (2013); Stenull et al. (2016); Mei et al. (2012); Zhang et al. (2018); Li et al. (2018); Esmann et al. (2018); Singh et al. (2018). In this paper we adapt the methods of TCQ to the analysis of phonon bands and, as a proof of principle, use them to unveil the existence of four previously unknown topological phases on the time-reversal symmetric planar honeycomb. These phases are characterized by the existence of edge and corner phonon modes and nontrivial windings in the Wilson line spectra.

One key concept in the TQC framework is that of (elementary) band representation Zak (1980, 1981); Cano et al. (2018a). Roughly speaking, a band representation describes how the set of atomic orbitals in a crystal transform under the symmetry operations of the space group of the crystal. As the number of atoms in an ideal crystal is infinite, a band representation is an infinite dimensional representation of the space group. The Bloch-wave combinations $\sum_{\mathbf{R}}\exp(i\mathbf{k}\cdot\mathbf{R})\phi_{i}(\mathbf{r}-\mathbf{R})$ of atomic orbitals foot (1) transform according to a representation of the little group $G_{\mathbf{k}}$ Bradley and Cracknell (1972). Thus a band representation induces little group representations at all the points in the Brillouin zone (BZ).

A band representation that can not be written as the sum of two band representations is an elementary band representation (EBR). With some exceptions that are well understood and tabulated Michel and Zak (2001); Cano et al. (2018a), an EBR is generated from a set of orbitals that transform under an irreducible representation (irrep) $D$ of the local site symmetry group $G_{w}$ of a maximal symmetry Wyckoff position (WP) $w$. The corresponding EBR will be denoted by $D@w$.

Any band representation is either elementary or can be written as a sum of EBRs. If a subset of bands is separated by a gap from the other bands and does not transform as a band representation, then the subset does not have an atomic limit and is topologically nontrivial Bradlyn et al. (2017). This gives an efficient method to identify topologically nontrivial subsets of bands. In practical terms, one considers the irreps describing the transformation of the subset of bands at the high symmetry points in the Brillouin zone. If the collection of little group irreps can not be obtained from a sum of EBRs, then the subset of bands necessarily has nontrivial topology.

This method is applied in Section II to the phonon bands on the honeycomb, and a dynamical matrix that incorporates all the symmetry-allowed couplings including third nearest neighbors is presented in Section III. The phase diagram for the dynamical matrix is explored in Section IV, where it is shown that, depending on the values of the coupling constant, the model can describe four distinct fragile topological phonon phases. The fact that the dynamical matrix includes all the couplings compatible with the system symmetries allows for very precise fits to real materials. As an example, we place graphene on the phase diagram and show that it is not very far from one of the four fragile topological phases. The nontrivial topological phases are characterized by their Wilson loop winding in Section V, where it is shown that zero winding number is forbidden by symmetry for any isolated two-band branch. Edge and corner modes are analyzed in Section VI, and the conclusions and discussion can be found in Section VII.

II Band representations for phonons

The concept of band representation is easily extended to the atomic vibrations of a crystal Michel et al. (1995); Walker and Zak (1995). Instead of atomic orbitals, one considers the set of vectors giving the displacements from equilibrium of all the atoms in the crystal and their transformation properties. This defines the “mechanical” band representation. Intuitively, the mechanical band representation can be understood as an electron band representation with (spinless) $p$-orbitals, since $p$-orbitals, like atomic displacements, transform according to the “vector representation” Bradley and Cracknell (1972) $V$. Thus the mechanical band representation $M$ is induced from the vector representation at each occupied Wyckoff position $w$ in the crystal i.e., $M=\sum_{w}V@w$.

An important difference between electron and phonon band representations concerns the number of degrees of freedom and their location. Electronic orbitals are functions defined throughout space, and we can consider many orbitals in an atom, while an atomic displacement is described just by a single vector at the atom location. As a consequence, the number of phonon bands is always equal to three times the number of atoms in the primitive unit cell. In this regard, the electromagnetic fields considered in the recent application of TQC to photonic crystals Blanco de Paz et al. (2019) are closer to the electron wave-functions than to the crystal vibrations studied here. Another property of phonon bands without parallel in electrons is the existence of three acoustic bands, that must satisfy the constraint $\lim_{\mathbf{k}\rightarrow 0}\omega(\mathbf{k})=0$. The acoustic modes at $\mathbf{k}=0$ represent uniform translations $\mathbf{t}$ of all the atoms in the crystal and transform according to the vector representation $V$.

II.1 Mechanical band representation for the planar honeycomb

The layer group for the planar honeycomb is LG 80 ($p6/mmm$), which corresponds to the space group SG 191 ($P6/mmm$). In-plane $(x,y)$ and off-plane $(z)$ vibrations of the lattice are respectively even and odd under reflections on the horizontal mirror plane, and decouple in the harmonic approximation. Thus, for practical purposes, it is sufficient and simpler to consider the subgroup SG 183 ($P6mm$), while treating in-plane and off-plane vibrations separately. Note, however, that the horizontal mirror plane symmetry is not broken. This means, in particular, that we are considering a suspended lattice as opposed to a lattice on a substratum.

The atoms in the honeycomb lattice are at the Wyckoff position $2b$, with site symmetry group isomorphic to the point group $3m$ ($C_{3v})$. The vector representation for $3m$ is reducible, $V=A_{1}(z)+E(x,y)$ Bradley and Cracknell (1972). This implies that the mechanical band representation can be written as the sum of two EBRs induced respectively from $A_{1}$ and $E$ at $2b$, i.e., $M=A_{1}@{2b}+E@{2b}$, that describe the transformation properties of the two off-plane and four in-plane bands.

The application BandRep at the Bilbao Christallographic Server (BCS)  Aroyo et al. (2006a, b, 2011) gives the irreps induced by any EBR at the high symmetry points in the BZ, and also tells whether the EBR is decomposable Bradlyn et al. (2018); Vergniory et al. (2017). An EBR is decomposable if the corresponding band can be split into two sets separated by a gap. It is known that a split EBR always gives rise to a topologically nontrivial phase Bradlyn et al. (2017); Bouhon et al. (2018). According to the BCS, $M=A_{1}@{2b}$ is indecomposable. In that sense the two off-plane phonon bands are nontopological and will not be considered further in this letter. Note, however, that they host a massless Dirac point at the $K$-point. On the other hand, $E@{2b}$ is decomposable and, as a consequence, the four in-plane bands can split into two disconnected sets. Table 1 gives, for all the subsets consistent with the compatibility relations, the irreps at the three high symmetry points in the BZ. The acoustic branches are identified by noting that the vector representation at the $\Gamma$ point is given by $V=\Gamma_{1}(z)+\Gamma_{6}(x,y)$ and therefore $\Gamma_{6}$ is associated with uniform in-plane translations.

The next step is to determine which of the branches in Table 1 can transform as band representations. To this end, we try to obtain sums of EBRs that induce the same little group irreps at the high symmetry points in the BZ. It turns out that this is possible only for the optical branch of phase IIa and the acoustic branch of IIb. As shown in Table 2, in order to obtain the irreps of the remaining branches we must subtract some EBRs. This means that only the optical branch of phase IIa and the acoustic branch of IIb may transform as band representations, whereas the remaining six branches must be topologically nontrivial. This is consistent with the general result that says that, when a set of bands transforming as an EBR splits, at least one of the resulting subsets of bands must have nontrivial topology Po et al. (2018); Cano et al. (2018b); Bradlyn et al. (2019). The negative integral coefficients in Table 2 are indicators of fragile topology, where the subset of bands may be trivialized by the addition of trivial bands Po et al. (2018).

The different phases in Table 1 are associated with the ordering of the phonon frequencies at points $K$ and $M$. For example, in order to obtain phase Ia, $\omega(K_{1})$ and $\omega(K_{2})$ must be lower than $\omega(K_{3})$, and $\omega(M_{1})$ and $\omega(M_{2})$ must also be lower than $\omega(M_{3})$ and $\omega(M_{4})$. This must be so in order to avoid band crossings between the two subsets, as such band crossings would cause the subsets to reconnect Bradlyn et al. (2018). For that reason, whenever $\omega(K_{3})$ is between $\omega(K_{1})$ and $\omega(K_{2})$, the four bands are interconnected and we are in a nontopological phase. The frequencies at different points in the BZ are obtained by diagonalizing the dynamical matrix Maradudin and Vosko (1968).

III Dynamical matrix for in-plane modes

The harmonic potential energy of the ions can be written as

[TABLE]

where $i,j$ run over all the atoms in the lattice, $\mathbf{r}_{i}=(x_{i},y_{i})$ is the displacement of atom $i$ from equilibrium and $U_{ij}=U_{ji}^{t}$ is the matrix of force constants between atoms $i$ and $j$. Assuming that the atoms $i,j$ are nth-nearest neighbors, the $2\times 2$ matrix $\mathrm{U}_{ij}$ is parametrized by four real coefficients $(a_{n},b_{n},c_{n},d_{n})$, often with a single set of coefficients for all nth-neighbor pairs. As explained in Appendix A, some of these coefficients may be forced to vanish by the symmetries of the lattice.

The basis vectors for the Bravais lattice are given by

[TABLE]

and are shown in Fig. 1 together with their reciprocal lattice vectors $\mathbf{g}_{i}$, which satisfy $\mathbf{g}_{i}\cdot\mathbf{e}_{j}=2\pi\delta_{ij}$. Sites on the $A(B)$ sublattice are at positions $\mathbf{R}+\boldsymbol{\delta}_{A(B)}$, where $\mathbf{R}$ denotes a lattice translation and

[TABLE]

Then each atom is specified by a composite label $i\equiv(a,\mathbf{R})$ that gives the unit cell $\mathbf{R}$ and the sublattice $a=A,B$. In what follows it will be convenient to make a change of coordinates from linear to circular polarizations

[TABLE]

As shown in Appendix A, in these coordinates, the matrix $U_{ij}$ can be written as

[TABLE]

where $s_{i}$ are Pauli matrices for phonon ‘spin’, with the upper (lower) components giving the amplitude of the right (left) circular polarization. The dynamical matrix is obtained as the Fourier transform of the matrix of force constants

[TABLE]

Note that the vector indices $(x,y)$ or $(\xi,\bar{\xi})$ are implicit.

The dynamical matrix acts on 4-component vectors

[TABLE]

and can be expressed in terms of Kronecker products of Pauli matrices $\sigma_{i}$ for the two sublattices, and $s_{i}$ for phonon ‘spin’. The result, up to third nearest neighbor interactions is given by

[TABLE]

[TABLE]

Note that all the coefficients in this expression are real. A detailed analysis of the symmetries of the dynamical matrix can be found in Appendix B.

It is interesting to note that, if we neglect third nearest neighbor interactions, the dynamical matrix (III) is closely related to the Kane-Mele hamiltonian Kane and Mele (2005b) for spinful $p_{z}$ orbitals in graphene, with the two phonon circular polarizations playing, to some extent, the role of electron spin. In fact, comparing (III) with the Kane-Mele hamiltonian shows that the terms proportional to $a_{1}$ have the structure of the Dirac hamiltonian, while the term with $d_{2}$ coincides with the Haldane spin-orbit coupling in the Kane-Mele model. On the other hand, the functions of ($k_{x},k_{y}$) in the terms proportional to $b_{1}$ are identical to the ones appearing in the Rashba spin-orbit coupling in the Kane-Mele model, although the matrix structures are different. Actually, writing the couplings in real space shows that the Rashba term in the Kane-Mele model differs from the $b_{1}$-term just by a factor of the imaginary unit $i$. That both terms manage to be time-reversal invariant reflects the fact that electron spin-$1/2$ and phonon spin-$1$ transform differently under time reversal.

Another consequence of the different transformation properties of electron and phonon spin, this time under space symmetries, is that the Rashba coupling breaks reflection symmetry by the horizontal mirror plane, while the $b_{1}$-term does not. Finally, one last manifestation of the differences between spin-$1/2$ and spin-$1$ systems is the absence of Kramers degeneracy in the phonon spectrum, where we find four different frequencies at the time-reversal symmetric $M$-point, instead of the two doubly degenerate energies of the electron spectrum.

As mentioned above, in-plane displacements transform according to the vector representation, just like $(p_{x},p_{y})$ orbitals do. As a consequence, the dynamical matrix for in-plane modes has to be closely related to the hamiltonian for spinless $(p_{x},p_{y})$-orbitals in graphene. This is actually the case, as the model in Ref. Cano et al., 2018b coincides with the dynamical matrix (III) if we take $(a_{1},b_{1})\sim(t_{\sigma}\pm t_{\pi})$, $d_{2}\sim x$, where $(t_{\sigma},t_{\pi},x)$ are the couplings in Ref. Cano et al., 2018b, and set all the other couplings and the on-site term (proportional to $\openone_{4}$) in (III) to zero. Note, however, that in that limit phases IIa and IIb, to be discussed in the next section, are out of reach and therefore were not mentioned in Ref. Cano et al., 2018b.

IV Phase diagram for the planar honeycomb

The dynamical matrix can be diagonalized analytically at the three high symmetry points in the BZ, and irreps can be assigned to the corresponding eigenmodes. The result may be written as

[TABLE]

where $D_{a}$ denotes an irrep of the little group $G_{\mathbf{k}}$, $\omega^{2}(\mathbf{k})$ is common to all the bands at the $\mathbf{k}$-point and $\delta\omega^{2}(\mathbf{k},D_{a})$ is specific to each irrep. We obtain

[TABLE]

while $\delta\omega^{2}(\mathbf{k},D_{a})$ is given in Table 3, where the last column shows the $C_{2z}$-eigenvalues of the normal modes at the $C_{2z}$-invariant points $\Gamma$ and $M$.

As mentioned above, in order to be in a nontrivial phase $\omega(K_{3})$ should be higher or lower than $\omega(K_{1})$ and $\omega(K_{2})$. A look at the frequencies in Table 3 shows that this implies

[TABLE]

A second necessary condition in order to have two disconnected sets of bands is that $\omega(M_{1})$ and $\omega(M_{2})$ are both higher or lower than $\omega(M_{3})$ and $\omega(M_{4})$. This is equivalent to

[TABLE]

For $|b_{1}|>2|b_{2}|$ this reduces to $2|b_{1}|<|a_{1}-3a_{3}|$ while, for $|b_{1}|<2|b_{2}|$, Eq. (12) is equivalent to $4|b_{2}|<|a_{1}-3a_{3}|$. The two cases can be combined into

[TABLE]

Assuming that Eqs. (11,13) are satisfied, Table 4 identifies the nontrivial topological phases. If either condition in Eqs. (11,13) fails to be satisfied, we will be in a nontopological phase where all four bands are interconnected. It is important to note that, unlike the hamiltonian $H(\mathbf{k})$ for electron bands, the dynamical matrix has to satisfy several stability conditions to prevent the existence of imaginary phonon frequencies. For instance, imposing $\omega^{2}(\Gamma)>0$ requires $a_{1}+a_{3}<0$.

Note also that phases Ia and Ib can be modeled by a dynamical matrix including only up to second nearest neighbors, whereas third nearest neighbor interactions are required to obtain phases IIa and IIb. Specifically, a minimal model that describes qualitatively the topologically nontrivial phases with WL winding number one (Ia and Ib, see next section) can be obtained by keeping only the set of parameters $(a_{1},b_{1},d_{2})$, with $|a_{1}|>2|b_{1}|$ and $2\sqrt{3}|d_{2}|>|b_{1}|$, with $a_{1}$ negative for stability. Phases IIa and IIb can be obtained by keeping also $a_{3}$, with $|a_{1}-3a_{3}|>2|b_{1}|$ and $3a_{3}<a_{1}$. Note that this last condition requires $3|a_{3}|>|a_{1}|$ which, as shown at the end of this section, would require increasing $a_{3}$ by a factor of twenty in the case of graphene. Thus having $3|a_{3}|>|a_{1}|$ is probably unrealistic for a material crystal, but may be attainable for mechanical systems or metamaterials Maldovan (2013).

IV.1 Locating graphene on the phase diagram

Here we will try to place graphene on the phase diagram for the in-plane phonon system of the planar honeycomb. To this end, we fit the DFT phonon spectrum reported in Ref. Natterer et al., 2015 to the one obtained from the third-nearest neighbor dynamical matrix (III).

Table 5 compares the phonon frequencies at the high symmetry points computed by DFT in Ref. Natterer et al., 2015 with the ones obtained by diagonalizing the third-nearest neighbor dynamical matrix (III) for the following set of parameters, in units of $(meV)^{2}$

[TABLE]

The parameters have been obtained by using the explicit formulae in Table 3 to fit the in-plane bands in Fig. 2 of Ref. Natterer et al., 2015 at the high symmetry points. As these are not labeled by irreps, some trial and error is necessary before the fit can be completed. The goal is to get the best possible agreement at the high symmetry points, while maintaining the overall qualitative features of the spectrum. This leads to the unique assignment of irreps at the high symmetry points given in Table 5. As shown in Table 4, the energies at the high symmetry points agree within a few percent, and the spectrum in Fig. 2 (left) closely resembles the one in Ref. Natterer et al., 2015.

All the in-plane bands in graphene are interconnected and the spectrum is nontopological. But now that we have parametrized the spectrum as in (14), we can try to see how far graphene is from becoming topological. Plugging (14) into (11) yields

[TABLE]

while for (13) the result is

[TABLE]

Thus we see that, whereas inequality (13) is on the verge of being satisfied, we are further from satisfying (11), although the couplings are at least the correct order of magnitude.

How much do we need to vary the parameters to get into a topologically non-trivial phase? Assume that we increase the strength of $d_{2}$ in the same proportion that we decrease $b_{1}$ and $b_{3}$, i.e., we take $d_{2}\rightarrow(1+x)d_{2}$, together with $b_{1}\rightarrow(1-x)b_{1}$ and $b_{3}\rightarrow(1-x)b_{3}$. Then we see that the required inequalities are satisfied for $x\gtrsim.52$. In other words, graphene can enter a topologically nontrivial phonon phase if some couplings are changed by about $50\%$. As graphene can be stretched by up to $20\%$ and coupling constants are often exponentially sensitive to atomic distances, this is not a wholly unrealistic possibility. Fig. 2(right) shows the in-plane phonon spectrum for $x=0.7$. As the frequency of $\omega(K_{3})$ is greater than $\omega(K_{1})$ and $\omega(K_{2})$, we are in Phase Ia. Note that phases IIa and IIb are clearly out of reach, as they would require increasing the strength of $a_{3}$ by a factor of twenty, from $|a_{3}|\simeq 0.11\cdot 10^{3}$ to more than $|a_{1}|/3=2\cdot 10^{3}$.

V Wilson loop winding and symmetries

Winding in the Wilson loop (WL) spectrum for a subset of isolated bands is a topological invariant. A winding that can not be eliminated by any perturbation that respects the symmetries of the system and does not close a gap guarantees that the subset of bands has nontrivial topology. We will consider a $\mathbf{g_{1}}$-directed Wilson loop Wilczek and Zee (1984); Soluyanov and Vanderbilt (2011); Zak (1982, 1989); Alexandradinata et al. (2016); foots (3) defined by

[TABLE]

where $P$ indicates that the integral is path-ordered and $A_{1}(\mathbf{k})_{ij}=i\langle u_{i}(\mathbf{k})|\partial_{k_{1}}u_{j}(\mathbf{k})\rangle$ is the non-abelian Berry potential built from the normal modes $u_{i}(\mathbf{k})$ of a subset of isolated bands. We take $k_{1}$ along $\mathbf{g_{1}}$ (see Fig. 1) and $k_{2}$ along $\mathbf{g_{2}}$. The eigenvalues of $W(k_{2})$ are of the form $e^{2\pi ix_{1}(k_{2})}$, where $x_{1}(k_{2})$ is the position of the center of a hybrid Wannier function Kivelson (1982); Marzari and Vanderbilt (1997) along $\mathbf{e_{1}}$. As the base-point moves along the $k_{2}$-axis from $\Gamma$ to $M$ and back to $\Gamma$, the Wannier centers move along the 1-dimensional unit cell, as shown in Fig. 3 foot (7).

According to Table 2, all the branches but the optical one of phase IIa and the acoustic one of IIb are necessarily topological, as they do not transform as band representations. On the other hand, the irreps of the optical branch of phase IIa and the acoustic branch of IIb at the high symmetry points of the BZ are such that they might transform according to band representations and could be trivial. However, as shown in Fig. 3, the WL of the acoustic branch of IIb has winding number two, which implies that the branch has nontrivial topology. In fact, we find that the WLs of all eight branches in Table 2 have non-zero winding numbers, equal to one for phases Ia and Ib, and two for IIa and IIb. In other words, there are no trivial bands in the topological phases of this phonon system.

The windings of the WLs can be understood from the $C_{2z}$-eigenvalues of the normal modes at the two $C_{2z}$-invariant points, $\Gamma$ and $M$. As shown in Table 3, the eigenvalues are $+1$ for $\Gamma_{5}$, $M_{1}$ and $M_{2}$, and $-1$ for $\Gamma_{6}$, $M_{3}$ and $M_{4}$. Comparing with Table 1, we see that the two bands in each branch of phases Ia and Ib have the same $C_{2z}$-eigenvalues at the $\Gamma$-point, and opposite to those at the $M$ point. The converse is true for phases IIa and IIb, where the $C_{2z}$-eigenvalues at the $\Gamma$ and $M$ points are equal. Having opposite $C_{2z}$-eigenvalues at $\Gamma$ and $M$ forces the Wilson bands to wind, while having the same eigenvalues is compatible with zero or, more generally, even winding numbers Alexandradinata et al. (2014). The crossings at generic $k_{2}$ and Wannier center $x_{1}=\pm 1/2$ in Fig. 3 are protected by $C_{2z}\mathcal{T}$ invariance Bouhon et al. (2018); Bradlyn et al. (2019); Song et al. (2018), where $\mathcal{T}$ is the time-reversal operation. A more involved analysis including the role of the $C_{3z}$-rotation symmetry Bouhon et al. (2018); Bradlyn et al. (2019); Song et al. (2018) shows that the allowed WL winding numbers are of the form $3n\pm 1$ ($n\in\mathbb{Z}$), see Appendix C for details. Thus zero winding is excluded, which explains the absence of trivial bands in the topological phases of this system.

We close this section by pointing out that our analysis of Wilson loop winding rests on the possibility of defining Wannier functions for phonons. In fact, the possibility of defining localized basis for lattice vibrations was established long ago by W. Kohn Kohn (1973) in a simple one-dimensional model, while a modern general treatment based on the extension of the concept of position operator to phonons can be found in Ref. Giustino and Pasquarello, 2006.

VI Edge and corner modes

In order to compute the edge modes we have imposed periodic boundary conditions in the $x$-direction (along $\mathbf{e_{1}}$ in Fig 1) and open boundary conditions on the two zig-zag edges. The results are presented in Fig 4, where $k_{x}$ goes from [math] to $2\pi$. The edge spectrum in phase Ib (IIb) is qualitatively identical to the spectrum in phase Ia (IIa). The edge modes present several peculiarities. First of all, besides the usual gap-crossing ‘optical’ edge modes, there are also ‘acoustic’ edge modes below the bulk acoustic branch. Secondly, the edge mode spectrum is gapped, as shown in the close-ups in Fig 4. This is typical of topologically fragile phases Po et al. (2018), and is consistent with the negative coefficients in the linear combinations of EBRs in Table 2.

Note that if we set $b_{1}=b_{2}=b_{3}=0$ in Eq. (III) the dynamical matrix becomes diagonal in ‘spin’ space, i.e., right and left-handed modes decouple form each other. In that case, the edge modes become helical modes and the gaps in their spectrum disappear. As the $b_{i}$ bulk couplings are allowed by all the symmetries in the system, this is a sign of fragile topology. In the gapless limit $b_{i}=0$, acoustic and optical edge modes have opposite helicities: if the acoustic right moving edge modes have right circular polarization, then then optical right moving edge modes have left circular polarization and vice versa. Notice also that the edge mode spectrum is more complicated in the IIa and IIb phases, possibly reflecting the winding number two in the Wilson line spectrum.

Lastly, there are also corner modes, as shown in Fig 5. This is again to be expected in topologically fragile phases, and is considered a signature of higher order topology Benalcazar et al. (2017, 2017); Wang et al. (2018); Wieder and Bernevig (2018). The corner modes in Fig 5 have been computed for a parallelogram cut along the unit vectors in Fig 1, with $40$ unit cells per side and open boundary conditions along the resulting zig-zag edges. The corner modes appear at the $60$ degree angles. For different values of the parameters we have also observed additional corner modes at the $120$ degree angles, sometimes buried in the continuum. The slow damping rate of the corner modes along the edges is a reflection of the smallness of the gap in the edge mode spectrum. Corner modes are present in all four topological phases, but they are easier to observe in phases Ia and Ib, i.e., in those with winding number one.

Concerning the topological stability of the corner modes, please see the analysis in Appendix D where an attempt is made to relate the observed corner states to the recently proposed filling anomaly Benalcazar et al. (2019). Although the existence of non-vanishing secondary topological indices for this system suggests that at least some of the numerically computed corner states must be topologically protected, our analysis is inconclusive as to which ones are topologically robust.

VII Discussion

Topological phononics is a very young field going back to less than a decade ago. In spite of the activity displayed, of which Refs. Wang et al. (2015a); Kariyado and Hatsugai (2015); Liu et al. (2017a); Wang et al. (2015b); Liu et al. (2017b); Süsstrunk and Huber (2015, 2016); Huber (2016); Kane and Lubensky (2013); Stenull et al. (2016); Mei et al. (2012); Zhang et al. (2018); Li et al. (2018); Esmann et al. (2018); Singh et al. (2018) are but a small sample, only a handful of systems with topological phonon spectra are presently known, most of them mechanical contraptions, and virtually no material realizations have been described. The results reported in this paper show that the methods of TQC can be adapted to a systematic search for topologically non-trivial phonons, with the crystal structure as input data. As these methods are susceptible of automatic implementation, this may eventually enlarge enormously our knowledge of materials which are the phononic analogues of topological insulators. foots (8)

One new insight from the analysis in this paper is that any material with the structure of the honeycomb has the potential to host topological phonons while keeping all its symmetries intact. This is, for phonons, the analogue of the Kane-Mele mechanism for electrons Kane and Mele (2005a, b). Indeed, the dynamical matrix studied in this paper can be considered as the fragile topology version of the Kane-Mele model. This should be contrasted with other approaches in the literature, where one tries to open a gap in the phonon spectrum by breaking time reversal symmetry, either through Coriolis forces Wang et al. (2015a); Kariyado and Hatsugai (2015); Liu et al. (2017a) or by interactions with external magnetic fields Holz (1972); Strohm et al. (2005); Sheng et al. (2006).

There are some open questions for the future. One is the precise definition of fragile topology for phonons. When the irreps of an isolated subset of bands can be obtained as a difference $BR_{1}-BR_{2}$ of band representations, the addition of a trivial band that transforms according to $BR_{2}$ ‘trivializes’ the fragile topology Po et al. (2018). In the case of electrons, the states transforming as $BR_{2}$ can, in principle, be found in the crystal, maybe as core orbitals deep in the valence bands, maybe forming high energy conduction bands. Obviously this is not the case for phonons, where the number of bands is limited by the crystal structure. Another question is the possible connectivity constraints derived from the existence of acoustic bands. This problem does not arise in graphene, where in- and off-plane modes decouple, but may be an issue with other systems.

Acknowledgements.

It is a pleasure to thank A. Bernevig, B. Bradlyn, J. Cano and Z. Song for useful discussions and suggestions. This work has been supported in part by Spanish Science Ministry grant PGC2018-094626-B-C21 (MCIU/AEI/FEDER, EU) and by Basque Government grant IT979-16.

Appendix A Parametrization of the dynamical matrix

The symmetry constraints on the matrix of force constants for the planar honeycomb have been considered in Ref. Falkovsky, 2008. For the sake a completeness and to set the notation we give here a brief summary of the analysis for the in-plane modes that, according to the analysis in Section II, are the ones that can exhibit nontrivial topology. The harmonic potential energy for the in-plane modes is given by

[TABLE]

where $i,j$ run over all the atoms in the lattice, $\mathbf{r}_{i}=(x_{i},y_{i})$ is the displacement of atom $i$ from equilibrium and $U_{ij}=U_{ji}^{t}$ is the matrix of force constants between atoms $i$ and $j$.

Assuming $i,j$ are nth-nearest neighbors, the $2\times 2$ matrix $\mathrm{U}_{ij}$ is parametrized by four real coefficients $(a_{n},b_{n},c_{n},d_{n})$. Pick one nth-nearest neighbor $j$ to the atom $i$ and parametrize $U_{ij}$ as

[TABLE]

where the parameters are real by time reversal symmetry. Then, if $(i^{\prime}j^{\prime})$ is another pair of nth-nearest neighbors such that $\mathbf{d}_{i^{\prime}j^{\prime}}=V(g)\mathbf{d}_{ij}$, where $V(g)$ is the vector representation for the symmetry operation $g$ that transforms $ij$ into $i^{\prime}j^{\prime}$, we will have

[TABLE]

If the number of nth-nearest neighbors is sufficiently large, we may have several subsets of pairs $(ij)$ unrelated by symmetry. In that case we would need more that one set of parameters $(a_{n},b_{n},c_{n},d_{n})$. For the honeycomb lattice this will happen whenever the number of nth-nearest neighbors is greater than $12$, which is the number of elements in $6mm$. On the other hand, the number of independent parameters for given $n$ may be less that four. This happens if the two atoms $(i,j)$ are left invariant (exchanged) by a symmetry element $g$, since then Eq. (20) (combined with $U_{ij}=U_{ji}^{t}$) becomes a constraint. For the honeycomb lattice, one possibility is to have a mirror $m_{\parallel}$ along the link $\mathbf{d}_{ij}$. Assuming for simplicity that $\mathbf{d}_{12}=d_{12}\mathbf{\hat{x}}$, under the action of $m_{\parallel}$ $(x_{i},y_{i})\rightarrow(x_{i},-y_{i})$ for $i=1,2$. This implies $c_{n}=d_{n}=0$. Another possibility is to have a mirror $m_{\perp}$ perpendicular to $\mathbf{d}_{ij}$ through the midpoint of the link. This exchanges the two atoms and reverses the $x$ components of the displacements, which implies $c_{n}=0$. The effect of these constraints is summarized in Table 6 for $n$ up to seven. It is important to realize that, in order to impose a constraint, the parallel mirror $m_{\parallel}$ has to pass through the two atoms. A mirror plane that is parallel to $\mathbf{d}_{ij}$ but does not go through $ij$ does not impose any constraint on $U_{ij}$; it merely relates $U_{i^{\prime}j^{\prime}}$ to $U_{ij}$, where $(i^{\prime}j^{\prime})$ are the images of $(ij)$ by the mirror plane. Failing to appreciate this point would lead us to set $d_{2}=0$, as in Ref. Falkovsky, 2008.

These are the only possible constraints on the planar honeycomb lattice for $i\neq j$. For $i=j$, the $C_{3z}$ invariance about each atom implies $U_{ii}=a_{0}\openone$. Note, however, that $a_{0}$ is not an independent parameter, due to the global translation invariance of the crystal, which implies $\sum_{j}U_{ij}=0$. This guarantees that atom $i$ does not experience any force when all the atoms in the crystal are given a uniform displacement $\mathbf{r}_{i}=\mathbf{t}$ and is behind the existence of acoustic phonon branches. This is a peculiarity of phonon dynamics without analog in electron hamiltonians, where on-site energies are independent parameters. For the planar honeycomb lattice, as a the result of global translation invariance, we have $a_{0}=-3a_{1}-6a_{2}-3a_{3}-\ldots$

Appendix B Symmetries of the dynamical matrix

The dynamical matrix is invariant under all the symmetries of the time-reversal symmetric honeycomb. The time-reversal operation may be written as $\mathcal{T}=\theta\mathcal{K}$, where $\mathcal{K}$ denotes complex conjugation and the unitary matrix $\theta$ is given by

[TABLE]

reflecting the fact that time-reversal exchanges the two circular polarizations but not the sublattices. Invariance of the dynamical matrix under time-reversal requires

[TABLE]

since $\mathcal{T}$ reverses the sign of the momenta.

By construction, the dynamical matrix $D(\mathbf{k})$ is periodic in reciprocal space, $D(\mathbf{k}+\mathbf{K})=D(\mathbf{k})$, where $\mathbf{K}$ is any vector in the reciprocal lattice. This follows trivialy from the definition (6). In order to consider the invariance of the dynamical matrix under point group operations, it is convenient to consider the alternative definition

[TABLE]

which is not periodic but is related to $D(\mathbf{k})$ by a unitary transformation

[TABLE]

where

[TABLE]

Then, invariance under a point group operation $R$ requires

[TABLE]

where $M(R)$ denotes the corresponding matrix in the mechanical representation. In-plane displacements of the atoms are invariant under the horizontal plane and for all practical purposes we can use the layer group $p6mm$. This is an infinite group that includes all the lattice translations, but we can restrict ourselves to the point group elements, as translation invariance has been implemented by taking the Fourier transform in Eq. (6). The point group $6mm$ can be generated by the sixfold counterclockwise rotation $C_{6z}^{+}$ and a vertical plane $m_{x}$ perpendicular to the $OX$ axis. The mechanical representation for these generators is given by

[TABLE]

where we have used the fact that in-plane displacements transform according to the $\Gamma_{6}$ irrep. The presence of $\sigma_{x}$ and $\openone_{\sigma}$ in the equation above reflects the fact that the sublattices are exchanged by $C_{6z}^{+}$ but not by $m_{x}$. The $\Gamma_{6}$ matrices for the generators are given by

[TABLE]

Imposing Eq. (26) for the two generators of the point group $6mm$ guarantees the invariance of the dynamical matrix under the whole layer group $p6mm$.

Appendix C Constraints on Wilson loop winding

C.1 $C_{2z}$-eigenvalues and WL winding

We mentioned in Section V that the winding of the WLs can be understood in terms of the $C_{2z}$-eigenvalues of the normal modes at the two $C_{2z}$-invariant points $\Gamma$ and $M$. Here we fill in the details of the argument following the ideas in Ref. Alexandradinata et al., 2014 and taking into account that, for the in-plane modes of the honeycomb lattice, $C_{2z}$ plays the role of inversion symmetry. The $\mathbf{g_{1}}$-directed WL is defined by

[TABLE]

For fixed $k_{2}$, this formula can be interpreted as the WL for a 1-dimensional system with reciprocal primitive cell $k_{1}\in[0,2\pi]$ along $\mathbf{g_{1}}$. Then, for values of $k_{2}$ such that the 1-dimensional system is $C_{2z}$-invariant, i.e., for $k_{2}=0,\pi$, the number $N_{(-1)}$ of $-1$ eigenvalues of the WL is given by Eq. (1) in Ref. Alexandradinata et al., 2014, namely

[TABLE]

where $n_{(\pm)}(0)$ and $n_{(\pm)}(\pi)$ are the numbers of normal modes with $C_{2z}$ eigenvalue equal to $\pm 1$ at $k_{1}=0$ and $k_{1}=\pi$ respectively. Similarly, the number of complex conjugate pairs of eigenvalues $(\lambda,\lambda^{*})$ of the WL is given by Alexandradinata et al. (2014)

[TABLE]

Now, as Fig. 6 shows, for $k_{2}=0$ the 1-dimensional WL goes through $\Gamma$ and $M^{\prime\prime}$, and we have

[TABLE]

This is equal to $2$ for phases Ia and Ib, where we have two modes with $C_{2z}$-eigenvalues $\eta$ at $\Gamma$ and $-\eta$ at $M$, and to [math] for phases IIa and IIb, where the eigenvalues are equal to $\eta$ at $\Gamma$ and $M$. Note that $M^{\prime\prime}$ and $M$ are related by a unitary $C_{3z}$-rotation that does not change the $C_{2z}$ spectrum. On the other hand, by Eq. (31) the number of complex conjugate pairs of eigenvalues $(\lambda,\lambda^{*})$ of the WL is zero in all cases. As the WL eigenvalues are equal to $e^{2\pi ix_{1}(k_{2})}$, where $x_{1}(k_{2})$ is the hybrid Wannier function center along $\mathbf{e}_{1}$, this implies $x_{1}(\Gamma)=1/2$ for phases Ia and Ib and $x_{1}(\Gamma)=0$ for phases IIa and IIb. This is exactly what is observed in Fig. 3.

For $k_{2}=\pi$ the 1-dimensional WL goes through $M$ and $M^{\prime}$, and we have

[TABLE]

irrespective of the topological phase. As above, Eq. (31) implies that the number of complex conjugate pairs of eigenvalues $(\lambda,\lambda^{*})$ of the WL is zero in all cases. Thus $x_{1}(M)=0$ for all the topological phases, in agreement with Fig. 3.

For phases IIa and IIb the hybrid Wannier functions sit at the same point for $k_{2}=0$ and $k_{2}=\pi$, and we could have non-winding WLs. This is obviously impossible for phases Ia and Ib, where the WLs necessarily must wind.

C.2 Crossings at generic points protected by $C_{2z}\mathcal{T}$

Having the same hybrid Wannier centers at $k_{2}=0$ and $k_{2}=\pi$ is compatible with non-winding WLs, but also with even winding number. For phases IIa and IIb we find the the WLs have winding number two for all the disconnected branches. As shown in Fig. 3, this involves crossings at generic values of $k_{2}$, and we might worry at their stability against small perturbations. Here we present for phonons the analog of the arguments given in Refs. Bouhon et al., 2018; Bradlyn et al., 2019; Song et al., 2018 for the stability of generic crossings for electrons. Note that for spin-$1/2$ systems $C_{2z}^{2}\!=\!\mathcal{T}^{2}\!=\!-1$, whereas for bosons we have $C_{2z}^{2}=\mathcal{T}^{2}=1$.

The combined operation $C_{2z}\mathcal{T}$ leaves $(k_{1},k_{2})$ invariant, and the WL loop must satisfy the constraint

[TABLE]

Writing the unitary operator $W(k_{2})$ in terms of the ‘Wannier hamiltonian’ $H_{W}$

[TABLE]

the constraint becomes

[TABLE]

where the extra minus sign is due to the fact that $C_{2z}\mathcal{T}$ is an antiunitary operation

[TABLE]

where $U$ is a unitary matrix and $\mathcal{K}$ denotes complex conjugation. In order to determine the form of the matrix $U$ we note that, according to Eq. (34), for generic $k_{2}$ the WL is invariant under the magnetic group $2^{\prime}$ that has only two elements, $2^{\prime}=\{E,C_{2z}\mathcal{T}\}$, where $E$ is the identity operation. Now, according to Ref. Bradley and Cracknell, 1972, the only single-valued irreducible corepresentation for the group $2^{\prime}$ is 1-dimensional, with the unitary matrix for $C_{2z}\mathcal{T}$ given by $D(C_{2z}\mathcal{T})=\pm 1$. The two signs give unitarily equivalent corepresentations, and we may take the $+$ sign without loss of generality. As we are considering two-band WLs, these must transform as the 2-dimensional corepresentation obtained by taking two copies of the 1-dimensional irreducible corepresentation. This means that $U=\openone_{2}$ and Eq. (36) reduces to

[TABLE]

Writing $H_{W}$ as a linear combination of Pauli matrices

[TABLE]

we see that the constraint (38) implies

[TABLE]

Had we taken two copies of the 1-dimensional irreducible representation with different signs for $D(C_{2z})$, the constraint would read

[TABLE]

with solution

[TABLE]

that is related to (40) by a unitary transformation. In any case, due the periodicity of the Wannier hamiltonian eigenvalues, there will be a crossing whenever $a(k_{2})=n\pi$, and this may happen for generic values of $k_{2}$. As we can not add another Pauli matrix, small perturbations will merely shift the position of the crossings, but will not be able to remove them. Note that $a(k_{2})=n\pi$ corresponds to $x_{1}(k_{2})=n/2$, i.e., the protected crossings will take place for $x_{1}=0,1/2$. Thus the generic crossings at $x_{1}=1/2$ in Fig 3 are protected by $C_{2z}\mathcal{T}$ and the winding of the WL is the winding of the function $a(k_{2})$.

C.3 The role of $C_{3z}$-invariance

As shown in Refs. Bouhon et al., 2018; Bradlyn et al., 2019; Song et al., 2018, in order to extract the consequences of $C_{3z}$-invariance on WL winding, one has to choose the paths of integration with some care. In this subsection we apply to our system the results in Ref. Bouhon et al., 2018, to which we refer the reader for a detailed presentation. There it is shown that the WL for a single isolated band along the green dashed path in Fig. 6 is given by

[TABLE]

where $\xi_{n}^{\mathbf{k}}$ is the $C_{nz}$-eigenvalue of the normal mode at the high symmetry point $\mathbf{k}$. As the loop encloses one sixth of the BZ, one can use the symmetries of the system Bouhon et al. (2018) to show that the total WL winding number for the single isolated band is given by

[TABLE]

It is important to note that, in order to obtain this result, one has to use other symmetries of the system besides $C_{3z}$. In their absence, as emphasized in Ref. Bradlyn et al., 2019, $C_{3z}$-invariance alone may not be enough to guarantee the existence of nontrivial windings for a particular WL, even though nontrivial windings will show up in other, appropriately chosen WLs.

In order to apply Eq. (43) to our system, where we have isolated pairs of bands, we turn off all the $b_{i}$ couplings. As mentioned in Section VI , this effectively decouples the right- and left-polarized sectors, and we end up with two 2-band systems with identical frequency spectra. This is shown in Fig. 7, which should be compared with Fig. 3. Note, in particular, that setting $b_{1}=0$ does not change the WL winding numbers. That is to be expected, as we can vary the values of the $b_{i}$ couplings adiabatically without closing the global gap that separates the acoustic and optical branches. This means that we can use Eq. (43) to compute the winding number for the isolated bands in the decoupled left and right-polarized systems, knowing that the results are valid also for $b_{i}\neq 0$. See Ref. Bouhon et al., 2018 for a thorough discussion of this point and further examples with spinful electrons.

As the dynamical matrix can be diagonalized analytically at the three high symmetry points, it is straightforward to compute the required $C_{nz}$-eigenvalues. We find that the WL for phases Ia and Ib is given by

[TABLE]

which implies $\gamma=\pm\pi/3+2\pi n$, where $n\in\mathbb{Z}$. By Eq. (44), this translates into

[TABLE]

The windings in Fig. 7 (top) are obtained by setting $n=0$. Similarly, the result for phases IIa and IIb is

[TABLE]

which implies $\gamma=\pm 2\pi/3+2\pi n$, yielding

[TABLE]

where the windings in Fig. 7 (bottom) are again obtained by setting $n=0$. Combining (46) and (48) gives $w=3n\pm 1$, which is the result quoted at the end of Section V. Thus zero winding number is incompatible with $C_{3z}$-invariance and any isolated branch, acoustic or optical, must be topologically nontrivial.

Appendix D Corner modes and filling anomaly

According to Ref. Benalcazar et al., 2019, corner-localized quantized features are protected by secondary topological indices that can be computed from the symmetry eigenvalues at the high symmetry points of the BZ of systems with $C_{n}$-symmetry. Those indices are a manifestation of the filling anomaly that, for electron bands, can be understood as a mismatch between the number of electrons required for charge neutrality and to preserve the crystal symmetry. This concept has been recently applied to other types of bands, such as those in photonic Wu et al. (2020); He et al. (2019) and acoustic crystals. Lin et al. (2020); Peri et al. (2019)

The appropriate secondary topological index for a $C_{6}$-symmetric crystal is given by Benalcazar et al. (2019)

[TABLE]

with

[TABLE]

where $\#\Pi_{p}^{(n)}$ counts the number of $e^{2\pi i(p-1)/n}$ eigenvalues of $C_{n}$ at the high symmetry point $\Pi$ and we take the $\Gamma$ point as a reference. The eigenvalues should be counted for all the bands below the gap or, more generally, for an isolated subset of bands. Eqs. (49) and (50) can be easily evaluated with the help of the matrices in Table 7 and the irreps in Table 1. The results are summarized in Table 8. Note that, according to Table 2, the two subsets of bands with vanishing $Q^{(6)}_{corner}$ can be induced from a Wannier function at the center of the hexagon, wich is the Wigner-Seitz primitive cell that preserves the $C_{6}$-symmetry of the honeycomb lattice.

Table 8 shows that, for at least one subset of bands in each phase, the topological index $Q^{(6)}_{corner}$ takes a non-vanishing value. This suggests that some of the observed corner modes must be topologically protected, although pointing out which ones may require a more detailed analysis that we leave for future work. Also for the future is left the task of clarifying the meaning of concepts such as charge neutrality and fractional charges in a phononic context. See Ref. van Miert and Ortix, 2020 for a very recent illustration of the subtleties involved in elucidating the topological character of corner modes.