Higher-Order Floquet Topological Insulators with Anomalous Corner States

# Higher-Order Floquet Topological Insulators with Anomalous Corner States

Biao Huang Department of Physics and Astronomy, University of Pittsburgh, Pittsburgh PA 15260, USA    W. Vincent Liu Department of Physics and Astronomy, University of Pittsburgh, Pittsburgh PA 15260, USA Wilczek Quantum Center, School of Physics and Astronomy and T. D. Lee Institute, Shanghai Jiao Tong University, Shanghai 200240, China Department of Physics and Astronomy, University of Pittsburgh, Pittsburgh PA 15260, USA Department of Physics and Astronomy, University of Pittsburgh, Pittsburgh PA 15260, USA Wilczek Quantum Center, School of Physics and Astronomy and T. D. Lee Institute, Shanghai Jiao Tong University, Shanghai 200240, China
July 26, 2019
###### Abstract

Higher order topological insulators have emerged as a new class of phases, whose robust in-gap “corner” modes arise from the bulk higher-order multipoles beyond the dipoles in conventional topological insulators. Despite rapid theoretical and experimental breakthroughs, all discussions have been constrained to the static scenario due to the lack of specific schemes to compute higher-order dynamical topological invariant. Here we provide a concrete model and explicit constructions of topological invariants for a Floquet-driven system exhibiting anomalous corner states. The bulk quadrupolar moment for the eigenstates of static Floquet operators vanishes identically, while the anomalous topological invariant associated with full-time evolution correctly describes the quantized corner charges. The signature of such a phase in cold atom experiments is discussed through corner particle dynamics and a Floquet-Bloch band tomography.

Introduction — The advent of topological insulators (TI) Hasan and Kane (2010); Qi and Zhang (2011); Ching-Kai Chiu and Ryu (2016) have revolutionized our understanding in the phases of matter by the principle of bulk-boundary correspondence. The non-trivial topology of the bulk wave-function is predicted by the principle to render robust edge states occupying a region one dimension lower than that of the bulk. Such a picture can be understood as a quantized dipole, or “first-order”, polarization of the bulk wave function which results in excessive charge at the system’s boundary. A wide class of topological systems have since been identified, including the TI Kane and Mele (2005a, b); Bernevig et al. (2006); Hsieh et al. (2009), the Weyl semimetals Xiangang Wan and Savrasov (2011); Lv et al. (2015); Xu et al. (2015); Yan and Felser (2017), quantum anomalous Hall effects Liu et al. (2016); Chang et al. (2013), and the topological superconductors Sato and Ando (2017); Zhang et al. (2018).

Analogous to electrodynamics, one naturally wonders about the generalization of dipoles to higher-order multipoles for bulk wave functions, which would modify the bulk-boundary correspondence. Such a novel class of “higher-order” topological insulators (HOTI) were constructed successfully in recent theories Benalcazar et al. (2017a, b); Kunst et al. (2018); Schindler et al. (2018a); Song et al. (2017); Langbehn et al. (2017); Trifunovic and Brouwer (); Calugaru et al. (); Queiroz and Stern (); Ezawa (), and have quickly led to experimental realizations in photonic Serra-Garcia et al. (2018); Peterson et al. (2018), electric circuits Imhof et al. () and solid state systems Wang et al. (); Schindler et al. (2018b). With vanishing dipole but quantized -th order multipoles in HOTI, both the bulk and edge exhibit gapped spectrum, while in-gap “corner” states emerge in a region being -dimensional lower than that of the bulk. Inspired by such success, rapid progress has been made towards higher order semimetals Lin and Hughes (); Ezawa (2018a, b), superconductivity Shapourian et al. (2018); Wang et al. (2018a); Khalaf (2018); Zhu (2018); Yan et al. (2018); Wang et al. (2018b), spin liquids Dwivedi et al. () and symmetry protected topological phases You et al. (); Rasmussen and Lu () in the past year.

So far, all discussions on HOTI have been focusing on static scenarios. It is known, however, that periodically driven (Floquet) systems far from equilibrium are fertile grounds for intriguing phenomena without static counterparts Nandkishore and Huse (2015); Eckardt (2017); Sacha and Zakrzewski (2017); Abanin et al. (2018). In particular, there exist “anomalous” Floquet insulators (AFI) Rudner et al. (2013); Roy and Harper (2017); Yao et al. (2017); Mukherjee et al. (2017); Peng et al. (2016); Maczewsky et al. (2016) whose static topological invariants vanish for all Floquet-Bloch bands, but edge states still emerge due to winding numbers of evolution operators with genuine time dependence. It is therefore tantalizing to explore the higher order extensions of AFI, which would open the door to a whole new set of non-equilibrium topological matters with multipole features. Yet, the current scheme for studying HOTI lacks a natural way to generalize the topological numbers built from Hamiltonian eigenstates to that from , prohibiting practical investigations of Bloch wave multipoles in dynamical systems. Clearly, an urgent need is posted to bridge the theoretical gap between HOTI and AFI, and to extend the experimental realization of AFI into the higher-order scenarios.

In this work, we explicitly construct the models and topological invariants for such a higher order Floquet topological insulator (HOFTI) exhibiting anomalous quadrupoles. Within this phase, the static nested polarization, constructed by replacing the Hamiltonian as in previous work with Floquet operator , vanish identically. But anomalous “corner” states still arise and are described by the dynamical quadrupoles contained in the full evolution of . The key step is to involve a Hermitian mapping of , which does not change the topological property of evolution operators while allowing for projections of onto the system’s boundary. Nested Wilson loops constructed from the evolution operators projected to the boundary correctly capture the dynamical edge topology and predict the quantized charge accumulating at the system’s corners, coinciding with numerical results. Further, a cold atom realization of a HOFTI is discussed, together with its signatures of anomalous corner states in detections. Our work paves the way for a systematic study on non-equilibrium topological matters of higher multipole nature.

Models, symmetries, and phase diagrams — As a minimal model, we consider a binary drive with Hamiltonians in two driving sectors illustrated in Fig. 1. To set the time reversal invariant point at time , we write them into a 3-step driving with period , where the Hamiltonian reads

 H(k,t)=⎧⎪⎨⎪⎩γh1,t∈[0,T/4];λh2k,t∈[T/4,3T/4];γh1,t∈[3T/4,T]. (1)

Here are hopping constants, and the (dimensionless) instantaneous Hamiltonians written in momentum space read

 h1 =τ1σ0−τ2σ2, h2k =coskxτ1σ0−sinkxτ2σ3−coskyτ2σ2−sinkyτ2σ1, (2)

where are Pauli matrices spanning the basis for 4 sublattices  111The direct product of Pauli matrices can be understood as, i.e. . as shown in Fig. 1. The evolution operator is then

 U(k,t)=Pτe−i∫t0H(k,τ)dτ, (3)

where denotes the path-ordering of time . Such a model enjoys high solvability sup () and rich phase diagrams including both the normal and anomalous Floquet topological phases.

The dynamical model defined in Eqs. (1)–(3) satisfies all of the time reversal , particle-hole , and chiral symmetries  Roy and Harper (2017); Yao et al. (2017)

 T−1U(k,t)T =U∗(−k,−t), T=τ0σ0; C−1U(k,t)C =U∗(−k,t), C=τ3σ0; S−1U(k,t)S =U(k,−t), S=τ3σ0. (4)

And . Here is complex conjugation , and are unitary matrices. Thus, the system belongs to the BDI class which holds no topological indices for the conventional first-order Floquet TIs Roy and Harper (2017); Yao et al. (2017) in two dimensions. That means the bulk dipoles always vanish, and the in-gap modes in an open-boundary system would be attributed to higher order multipoles.

We are interested in the system’s characters at spectroscopic time , with being integers. The bulk spectrum of the Floquet operator : can be obtained as sup ()

 exp(iEk±) =exp(±iEk)=fk±i√1−f2k, fk=cosEk =cos(√2γ)cos(√2λ) −coskx+cosky2sin(√2γ)sin(√2λ), (5)

with each quasi-energy band being two-fold degenerate. The gap closes when , giving the topological phase boundaries

 √2λ=±√2γ+mπ,m∈Z. (6)

Therefore, one can divide the irreducible phase diagram into four distinct regions as shown in Fig. 2(a) 222Note that the axis lines are not phase boundaries except for the special points satisfying Eq. (6).. To characterize the topological properties, two numerical results are presented for each phase 333The lattice sizes in Fig. 2(c) are for the spectrum, and for transverse polarization . To break the full degeneracy of corner states, we put in a static, infinitesimal onsite chemical potential biases for sublattices 1, 2, 3 and 4. The choice is to fix the signs of polarizations at and ..

First, exact diagonalization of the real-space Floquet operator with open boundary conditions are shown in Fig. 2(c). Enforced by particle-hole symmetries, the boundary modes only exist at quasienergy and/or , represented by the red dots. Each set of boundary modes within one bulk gap involves four eigenstates localizing at the four corners of the sample respectively. Fig. 2(b) illustrates the amplitudes of one such mode.

Second, we compute the transverse polarization in a semi-infinite system (described by ), which represents the shift of wave-function centers towards -direction away from lattice sites Benalcazar et al. (2017a, b); sup (),

 Wy(x,ky) =Pk′ye−i∮kydk′yAy(x,k′y)=∑j|νj(x,ky)⟩e2πiνj(x)⟨νj(x,ky)|, py(x) =1Ly∑m,j,ky|[νj(x,ky)]mum(x,ky)|2νj(x). (7)

Here is the Wilson loop with the Bloch-wave Berry connection , whose non-Abelian components . in the first equation is the base-point for the loop integration . ’s are eigenstates of in lower or upper bands, which render ’s of opposite signs. means path-ordering. If contains bulk quadrupoles, there would appear non-zero near the edges as it represents “nested” polarization resulting in corner modes. For parameters near phase transitions, corner modes could be more extended into the bulk, but the half-system total is quantized to be half-integer (modulo integer) in the topological quadrupolar phase.

From the numerical results, one can readily identify that 1⃝ and 2⃝ in Fig. 2 are the (normal) topological and trivial insulating phases respectively. The in-gap corner modes in 1⃝ arise from non-trivial polarization appearing at the system’s edge. In fact, 1⃝ and 2⃝ are smoothly connected to the static phases in Refs. Benalcazar et al. (2017a, b), as when in 1⃝ (or in 2⃝ ), the Floquet model Eq. (1) is described by a static Hamiltonian (or ) in Eq. (2), which holds non-trivial (or trivial) Bloch-wave quadrupoles Benalcazar et al. (2017a). Further, the phase 3⃝ is also a normal topological one, but with opposite polarization compared with 1⃝ and the corner modes appear at gap. The representative at is equivalent to up to a global gauge transformation by .

The most interesting phenomenon occurs in phase 4⃝ . The Floquet operator reduces to with the representative parameters , which seems like a topological trivial one. Indeed, vanishes identically in this case, indicating that the static Floquet operator involves no quantized quadrupoles. However, the corner modes do show up in both the gaps. The contradiction signals that the anomalous corner modes in phase 4⃝ result from the quantized quadrupoles associated with the full dynamics of throughout a period, as we will discuss next.

Topological invariant — To define a topological number with time being an independent parameter, we need to periodize the evolution operator in time such that functions as in the parameter space. This can be constructed by using the return map Rudner et al. (2013); Roy and Harper (2017); Yao et al. (2017):

 Uε(t)=U(t)e−iH(ε)efft/T, (8)

where , and denotes the branch cut when taking the logarithm. Here ’s are eigenstates of the Floquet operator , and therefore we have by construction because . Choices of branch cuts determine the non-analytic point in quasienergy spectrums, which is the gap that we check the possible existence of corner modes.

Let us focus on class BDI related to our model above. In the presence of chiral symmetry, at half evolution period takes block diagonal/off-diagonal forms depending on branch cuts  Ching-Kai Chiu and Ryu (2016); Yao et al. (2017)

 Uε=0(T2)=(0U+U−0), Uε=π(T2)=(U+00U−), (9)

where are unitary matrices. We emphasize that these operators still carry the information of the full evolution through and is not simply a static one at half period.

Up to now, the procedures are standard for an AFI Rudner et al. (2013); Roy and Harper (2017); Yao et al. (2017). Our major result is the scheme to analyze the winding number of the unitary within the subspaces of certain Wannier bands. To do so, we introduce the tool Hermitian operator 444The similar “Hermitian map” is introduced in Roy and Harper (2017) to discuss the first order Floquet topological insulators, and the one-to-one correspondence between Hermitian maps and , and the equivalence of topology for and for the Hermitian map are proved there.

 Htool=(0U†−U−0). (10)

Note , the eigenvalues of the Hermitian matrix are . Thus, it introduces a map of , where denotes torus and is the spatial dimension. The solution of can be written as

 |m±⟩=1√2(±|αm⟩U−|αm⟩), (11)

where is an arbitrary normalized U(N) spinor , and denotes eigenvalues . Since the is a trivial reproduction of the properties of the U(N), we can focus on one of the two branches, say, the branch, and simplify the notation . The eigenvector serves to introduce the Berry connection of within the subspace spanned by selected :

 i[Aμ]mn=⟨m|∂kμ|n⟩=12⟨αm|(U†−∂kμU−)|αn⟩+⟨αm|∂kμ|αn⟩. (12)

Without any subspace projection, one natural choice is , where the subscript denotes the -th element of the U(N) spinor, and ranges over the whole U(N) space. Then we recover the usual Berry connection for the unitary matrix

 i[Aμ]mn=12(U†−∂kμU−)mn, (13)

where the on the right-hand-side denotes matrix indices.

With the aid of Eq. (12), we can compute the higher order winding number using the following procedure in the discrete Brillouin zone.

(1) Calculate the (first-order) Wilson loop with “bare” Berry connections defined in Eq. (13),

 Fx,k =12(I+U†−,kU−,k+Δkex),Δkx=2πLx, Wx,k =Fx,kFx,k+Δkxex⋯Fx,k+(Lx−1)Δkxex. (14)

For our model, is an SU(2) matrix. It represents the polarization of Bloch waves towards -direction.

(2) Diagonalize the first order Wilson loop ,

 Wx,k|νj(k)⟩=ei2πνj(ky)|νj(k)⟩. (15)

The phases are the Wannier band spectrum, and ’s carry the information of edge topology.

(3) Obtain the nested Wilson loop as

 ~F(j)yk =12⟨νj(k)|U†−,kU−,k+Δkyey|νj(k)⟩ +⟨νj(k)|νj(k+Δkyey)⟩−12,Δky=2πLy, ~W(j)yk =~F(j)yk~F(j)yk+Δkyey⋯~F(j)yk+(Ly−1)Δkyey. (16)

Here the Wilson loop has projected Berry connection defined in Eq. (12), where is replaced by the Wannier wave functions . The nested Wilson loop involves the simultaneous polarization of Bloch waves towards the and direction, which leads to the corner states. In our case, and are two U(1) numbers and therefore no further diagonalization is needed 555If is a matrix, one needs to diagonalize it similar to Eq. (15) and obtain the phase argument as for in step (2).. The nested polarization can be obtained as

 ~P(j)xy=1Lx∑kx~p(j)y(kx); ~p(j)y(kx)=12πarg(~W(j)yk) (17)

Note that an arbitrary in can be taken because the phase factor does not depend on the base point after the (path-ordered) integration over in Eq. (16).

We apply the above procedure and compute the quadrupole strength of our model, as shown in Fig. 3. We see that for phase 4⃝ , both the and gap involves two gapped Wannier bands (denoted as ). Therefore, each Wannier band carries its own topological number. Indeed, identically for all up to numerical accuracy, and therefore the quadrupole strength for both gaps. Thus, correctly captures the corner charges as shown in Fig. 2(c) for phase 4⃝ .

Similar to the first order Floquet TI’s Rudner et al. (2013); Roy and Harper (2017); Yao et al. (2017), one can write the relation between dynamical and static quadrupoles as

 ~P(ε1)xy−~P(ε2)xy=~P(ε2

where , are dynamical quadrupoles computed above for gaps respectively, and is the static quadrupole computed by for the band sandwiched by and . This explains that the static quadrupoles in phase 4⃝ vanish for all Floquet-Bloch bands (see Fig. 2(c)) because , as shown in Fig. 3. One can also verify this relation for other phases, where gaps with (or without) corner modes correspond to (or ).

Experimental proposals — First, we point out that the photonic experiments on the static HOTIs Serra-Garcia et al. (2018); Peterson et al. (2018) can be generalized to our HOFTI case with little modifications. The and in Eqs. (1) have already been realized in these experiments, and one may extend them by performing repeated quenches between the two Hamiltonians resembling the situation for first order AFIs Mukherjee et al. (2017); Maczewsky et al. (2016); Cheng et al. (2018). The corner modes can be observed directly through the photon intensity at sample corners. In contrast to the static HOTI which has only one in-gap frequency peak for corner modes, our theory predicts that HOFTI would exhibit two such peaks both between and aside of the bulk modes.

Second, we discuss the signatures of HOFTI in cold atom experiments 666Note Ref Benalcazar et al. (2017a) has already elaborated the realization of with -flux Hofstadter models in the presence of optical superlattices, and the Floquet driving in our case corresponds to repeated quenching the relative lattice depth for superlattices. Thus, we only focus on the detection features missing in the previous literature. . Here, the particle conservation could lead to novel dynamics showing the localization of corner modes, and a time-of-flight (TOF) density mapping (band tomography) could reveal the trivial static polarization in the HOFTI phase. See the results of simulations in Fig. 4.

A “corner” can be engineered via the microscope methods associated with digital mirror device Tai et al. (2017); Quelle et al. (2017), which gives rise to a step-like potential barrier up to single site accuracy. With high enough , it can be represented by the open boundary conditions in our simulations. We consider the initial states with particles concentrating around the sample corner 777Here, as an illustration, we consider uniform initial particle distribution at the corner unit cells (totally 36 sites). None-uniform distributions do not change the qualitative features. Experimentally, the initial state can be prepared by an off-center harmonic trap rendering the corner sites lowest chemical potentials Aidelsburger et al. (2015), or direction moving fermion particles to the edge via microscope-based optical tweezers Endres et al. (2016), which overlap with both the corner and bulk eigenstates. Since the bulk spectrum is dispersive, systems with/without corner states would have finite/zero density remaining at the corner after long-time evolution in an in situ imaging, as shown in Fig. 4(a)–(e).

To distinguish corner states in HOFTI from the static ones in HOTI, one can apply band tomography and map out all the eigenstates for , from which the bulk nested polarization can be backed up. We briefly describe the procedures below and leave more details in SM sup (). and therefore its eigenstates can be parameterized by three angles in . Working in the regimes for phase 1⃝ or for 4⃝, where eigenstates connect smoothly to those in static systems, one can populate the “lowest two” Floquet bands even in the slow-driving situation. After equilibration, the lattice depth is ramped up and the system evolves under static chemical potentials for sublattices . Finally, momentum density is measured after TOF. As expected, three profiles of : could fully reveal  888In fact, there are infinitely many choices to back up the angles. We demonstrate one such choice in sup (). It is clear from Fig. 4(f) that the weak momentum dependence of the angles in phase 4⃝ signals the vanishing static quadrupoles.

Discussions and generalizations — We investigate in detail a model in the symmetry class BDI showing the Floquet anomalous corner state, and illustrates its signatures in cold atom experiments. Moreover, a general framework is constructed to compute the anomalous quadrupoles in such a genuine dynamical system. The three steps related to Eqs. (14)–(17) can be readily generalized to systems with more bands and with higher order multipoles. Also, since we only assume the presence of chiral symmetry (which gives Eq. (9)), the procedure can be directly applied to symmetry class AIII, DIII, CI and CII for dimensions with indices. Without chiral symmetry, one generally needs to tackle for all , but the mapping of in Eq. (11) remains valid, which allows for a projection of to subspaces by constraining . It will be interesting for future work to explore such a generalization, i.e. to class A with quadrupoles in 3D, and to provide a complete list of procedures to obtain higher order dynamical multipole in all symmetry classes and dimensions.

Acknowledgment — This work is supported by AFOSR Grant No. FA9550-16-1-0006, ARO Grant No. W911NF-11-1-0230, and NSF of China Overseas Scholar Collaborative Program Grant No. 11429402 sponsored by Peking University.

## References

• Hasan and Kane (2010) M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82, 3045 (2010).
• Qi and Zhang (2011) X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys. 83, 1057 (2011).
• Ching-Kai Chiu and Ryu (2016) A. P. S. Ching-Kai Chiu, Jeffrey C.âY. Teo and S. Ryu, Rev. Mod. Phys. 88, 035005 (2016).
• Kane and Mele (2005a) C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95, 146802 (2005a).
• Kane and Mele (2005b) C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95, 226801 (2005b).
• Bernevig et al. (2006) B. A. Bernevig, T. L. Hughes,  and S.-C. Zhang, Science 314, 1757 (2006)cond-mat/0611399v1 .
• Hsieh et al. (2009) D. Hsieh, D. Qian, L. Wray, Y. Xia, Y. S. Hor, R. J. Cava,  and M. Z. Hasan, Nature 452, 970 (2009)0902.1356v1 .
• Xiangang Wan and Savrasov (2011) A. V. Xiangang Wan, Ari M. Turner and S. Y. Savrasov, Phys. Rev. B 83, 205101 (2011).
• Lv et al. (2015) B. Q. Lv, H. M. Weng, B. B. Fu, X. P. Wang, H. Miao, J. Ma, P. Richard, X. C. Huang, L. X. Zhao, G. F. Chen, Z. Fang, X. Dai, T. Qian,  and H. Ding, Phys. Rev. X 5, 031013 (2015)1502.04684v3 .
• Xu et al. (2015) S.-Y. Xu, I. Belopolski, D. S. Sanchez, C. Guo, G. Chang, C. Zhang, G. Bian, Z. Yuan, H. Lu, Y. Feng, T.-R. Chang, P. P. Shibayev, M. L. Prokopovych, N. Alidoust, H. Zheng, C.-C. Lee, S.-M. Huang, R. Sankar, F. Chou, C.-H. Hsu, H.-T. Jeng, A. Bansil, T. Neupert, V. N. Strocov, H. Lin, S. Jia,  and M. Z. Hasan, Science 1, 10 (2015)1508.03102v2 .
• Yan and Felser (2017) B. Yan and C. Felser, Annual Review of Condensed Matter Physics 8, 337 (2017).
• Liu et al. (2016) C.-X. Liu, S.-C. Zhang,  and X.-L. Qi, Annu. Rev. Condens. Matter Phys. 7, 301 (2016)1508.07106v1 .
• Chang et al. (2013) C.-Z. Chang et al.Science 340, 167 (2013).
• Sato and Ando (2017) M. Sato and Y. Ando, Rep. Prog. Phys. 80, 076501 (2017)1608.03395v3 .
• Zhang et al. (2018) P. Zhang, K. Yaji, T. Hashimoto, Y. Ota, T. Kondo, K. Okazaki, Z. Wang, J. Wen, G. D. Gu, H. Ding,  and S. Shin, Science 360, 182 (2018)1706.05163v2 .
• Benalcazar et al. (2017a) W. A. Benalcazar, B. A. Bernevig,  and T. L. Hughes, Science 357, 61 (2017a).
• Benalcazar et al. (2017b) W. A. Benalcazar, B. A. Bernevig,  and T. L. Hughes, Phys. Rev. B 96, 245115 (2017b)1708.04230v2 .
• Kunst et al. (2018) F. K. Kunst, G. van Miert,  and E. J. Bergholtz, Phys. Rev. B 97, 241405 (2018)1712.07911v3 .
• Schindler et al. (2018a) F. Schindler, A. M. Cook, M. G. Vergniory, Z. Wang, S. S. Parkin, B. A. Bernevig,  and T. Neupert, Sci. Adv. 4, eaat0346 (2018a).
• Song et al. (2017) Z. Song, Z. Fang,  and C. Fang, Phys. Rev. Lett. 119, 246402 (2017)1708.02952v5 .
• Langbehn et al. (2017) J. Langbehn, Y. Peng, L. Trifunovic, F. von Oppen,  and P. W. Brouwer, Phys. Rev. Lett. 119, 246401 (2017)1708.03640v2 .
• (22) L. Trifunovic and P. Brouwer,  1805.02598v1 .
• (23) D. Calugaru, V. Juricic,  and B. Roy,  1808.08965v1 .
• (24) R. Queiroz and A. Stern,  1807.04141v2 .
• (25) M. Ezawa,  1806.03007v1 .
• Serra-Garcia et al. (2018) M. Serra-Garcia, V. Peri, R. Süsstrunk, O. R. Bilal, T. Larsen, L. G. Villanueva,  and S. D. Huber, Nature 555, 342 (2018)1708.05015v1 .
• Peterson et al. (2018) C. W. Peterson, W. A. Benalcazar, T. L. Hughes,  and G. Bahl, Nature 555, 346 (2018)1710.03231v1 .
• (28) S. Imhof, C. Berger, F. Bayer, J. Brehm, L. Molenkamp, T. Kiessling, F. Schindler, C. H. Lee, M. Greiter, T. Neupert,  and R. Thomale,  1708.03647v1 .
• (29) Z. Wang, B. J. Wieder, J. Li, B. Yan,  and B. A. Bernevig,  1806.11116v1 .
• Schindler et al. (2018b) F. Schindler, Z. Wang, M. G. Vergniory, A. M. Cook, A. Murani, S. Sengupta, A. Y. Kasumov, R. Deblock, S. Jeon, I. Drozdov, H. Bouchiat, S. GuÃ©ron, A. Yazdani, B. A. Bernevig,  and T. Neupert, Nat. Phys. 14, 918 (2018b)1802.02585v1 .
• (31) M. Lin and T. L. Hughes,  1708.08457v1 .
• Ezawa (2018a)
• Ezawa (2018b) M. Ezawa, Phys. Rev. B 97, 155305 (2018b)1802.03571v1 .
• Shapourian et al. (2018) H. Shapourian, Y. Wang,  and S. Ryu, Phys. Rev. B 97, 094508 (2018)1711.02122v2 .
• Wang et al. (2018a) Y. Wang, M. Lin,  and T. L. Hughes,  (2018a), 1804.01531v2 .
• Khalaf (2018) E. Khalaf, Phys. Rev. B 97, 205136 (2018)1801.10050v4 .
• Zhu (2018)
• Yan et al. (2018) Z. Yan, F. Song,  and Z. Wang, Phys. Rev. Lett. 121, 096803 (2018)1803.08545v2 .
• Wang et al. (2018b) Q. Wang, C.-C. Liu, Y.-M. Lu,  and F. Zhang,  (2018b), 1804.04711v1 .
• (40) V. Dwivedi, C. Hickey, T. Eschmann,  and S. Trebst,  1803.08922v1 .
• (41) Y. You, T. Devakul, F. J. Burnell,  and T. Neupert,  1807.09788v2 .
• (42) A. Rasmussen and Y.-M. Lu,  1809.07325v1 .
• Nandkishore and Huse (2015) R. Nandkishore and D. A. Huse, Annu. Rev. Condens. Matter Phys. 6, 15 (2015).
• Eckardt (2017) A. Eckardt, Rev. Mod. Phys. 89, 011004 (2017).
• Sacha and Zakrzewski (2017) K. Sacha and J. Zakrzewski, Reports on Progress in Physics 81, 016401 (2017).
• Abanin et al. (2018) D. A. Abanin, E. Altman, I. Bloch,  and M. Serbyn, arXiv preprint arXiv:1804.11065  (2018).
• Rudner et al. (2013) M. S. Rudner, N. H. Lindner, E. Berg,  and M. Levin, Phys. Rev. X 3, 031005 (2013).
• Roy and Harper (2017) R. Roy and F. Harper, Phys. Rev. B 96, 155118 (2017)1603.06944v2 .
• Yao et al. (2017) S. Yao, Z. Yan,  and Z. Wang, Phys. Rev. B 96, 195303 (2017)1708.05993v2 .
• Mukherjee et al. (2017) S. Mukherjee, A. Spracklen, M. Valiente, E. Andersson, P. Ãhberg, N. Goldman,  and R. R. Thomson, Nat. Commun. 8, 13918 (2017)1604.05612v2 .
• Peng et al. (2016) Y.-G. Peng, C.-Z. Qin, D.-G. Zhao, Y.-X. Shen, X.-Y. Xu, M. Bao, H. Jia,  and X.-F. Zhu, Nat. Commun. 7, 13368 (2016).
• Maczewsky et al. (2016) L. J. Maczewsky, J. M. Zeuner, S. Nolte,  and A. Szameit, Nat. Commun. 8, 13756 (2016)1605.03877v1 .
• (53) The direct product of Pauli matrices can be understood as, i.e. .
• (54) See supplemental materials for details.
• (55) Note that the axis lines are not phase boundaries except for the special points satisfying Eq. (6).
• (56) The lattice sizes in Fig. 2(c) are for the spectrum, and for transverse polarization . To break the full degeneracy of corner states, we put in a static, infinitesimal onsite chemical potential biases for sublattices 1, 2, 3 and 4. The choice is to fix the signs of polarizations at and .
• (57) The similar “Hermitian map” is introduced in Roy and Harper (2017) to discuss the first order Floquet topological insulators, and the one-to-one correspondence between Hermitian maps and , and the equivalence of topology for and for the Hermitian map are proved there.
• (58) If is a matrix, one needs to diagonalize it similar to Eq. (15) and obtain the phase argument as for in step (2).
• Cheng et al. (2018) Q. Cheng, Y. Pan, H. Wang, C. Zhang, D. Yu, A. Gover, H. Zhang, T. Li, L. Zhou,  and S. Zhu, Nat. Commun.  (2018), 1804.05134v2 .
• (60) Note Ref Benalcazar et al. (2017a) has already elaborated the realization of with -flux Hofstadter models in the presence of optical superlattices, and the Floquet driving in our case corresponds to repeated quenching the relative lattice depth for superlattices. Thus, we only focus on the detection features missing in the previous literature.
• Tai et al. (2017) M. E. Tai, A. Lukin, M. Rispoli, R. Schittko, T. Menke, D. Borgnia, P. M. Preiss, F. Grusdt, A. M. Kaufman,  and M. Greiner, Nature 546, 519 (2017).
• Quelle et al. (2017) A. Quelle, C. Weitenberg, K. Sengstock,  and C. M. Smith, New J. Phys. 19, 113010 (2017).
• (63) Here, as an illustration, we consider uniform initial particle distribution at the corner unit cells (totally 36 sites). None-uniform distributions do not change the qualitative features. Experimentally, the initial state can be prepared by an off-center harmonic trap rendering the corner sites lowest chemical potentials Aidelsburger et al. (2015), or direction moving fermion particles to the edge via microscope-based optical tweezers Endres et al. (2016).
• (64) In fact, there are infinitely many choices to back up the angles. We demonstrate one such choice in sup ().
• Aidelsburger et al. (2015) M. Aidelsburger, M. Lohse, C. Schweizer, M. Atala, J. T. Barreiro, S. Nascimbene, N. Cooper, I. Bloch,  and N. Goldman, Nature Physics 11, 162 (2015).
• Endres et al. (2016) M. Endres, H. Bernien, A. Keesling, H. Levine, E. R. Anschuetz, A. Krajenbrink, C. Senko, V. Vuletic, M. Greiner,  and M. D. Lukin, Science 354, 1024 (2016).
• Hauke et al. (2014) P. Hauke, M. Lewenstein,  and A. Eckardt, Phys. Rev. Lett. 113, 045303 (2014).
• Fläschner et al. (2016) N. Fläschner, B. Rem, M. Tarnowski, D. Vogel, D.-S. Lühmann, K. Sengstock,  and C. Weitenberg, Science 352, 1091 (2016).
• Tarnowski et al. (2017) M. Tarnowski, M. Nuske, N. Fläschner, B. Rem, D. Vogel, L. Freystatzky, K. Sengstock, L. Mathey,  and C. Weitenberg, Phys. Rev. Lett. 118, 240203 (2017)1703.02813v2 .

Supplemental Material

Supplemental Materials for “Higher-Order Floquet Topological Insulators with Anomalous Corner States” Biao Huang W. Vincent Liu July 26, 2019

## Appendix S-1 Algebraic details for evolution operators

### S-1.1 Floquet operator

First, we consider the bare evolution operator. Set which means . When :

 U(t)=cos(√2γt)−isin(√2γt)h1√2; (S1)

when ,

 U(t) =(cos(√2λ(t−1/2))−isin(√2λ(t−1/2))h2k√2)(cos(γ/√2)−isin(γ/√2)h1√2); (S2)

and finally, when ,

 U(t) =(cos(√2γ(t−3/2))−isin(√2γ(t−3/2))h1√2)(cos(√2λ)−isin(√2λ)h2k√2)(cos(γ/√2)−isin(γ/√2)h1√2). (S3)

The Floquet operator is the evolution operator at the end of a full period:

 UF=U(T)=(cos(√22γ)−isin(√22γ)h1√2)(cos(√2λ)−isin(√2λ)h2k√2)(cos(√22γ)−isin(√22γ)h1√2). (S4)

Using the relations

 h1h2k+h2kh12=(coskx+cosky)τ0σ0, h1h2kh12=coskyτ1σ0+sinkxτ2σ3−coskxτ2σ2+sinkyτ2σ1, (S5)

we have the analytical form for the Floquet operator as

 UF=fkI+i(g1kΓ1+g2kΓ2+g3kΓ3+g4kΓ4), (S6)

where

 fk =cos√2γcos√2λ−sin√2γsin√2λcoskx+cosky2 g1k =1√2sin√2λsinky g2k =1√2(sin√2γcos√2λ+cos√2γsin√2λcoskx+cosky2−sin√2λcoskx−cosky2) g3k =1√2sin√2λsinkx g4k =−1√2(sin√2γcos√2λ+cos√2γsin√2λcoskx+cosky2+sin√2λcoskx−cosky2) (S7)

and are Dirac matrices being anticommuting with each other , and ’s are real numbers satisfying . That means , and therefore its eigenstates, can be parameterized by three angles in addition to the quasienergy : . The eigenvalues of can be obtained via , which gives in our case ,

 exp(iEk±) =exp(±iEk)=fk±i√1−f2k, fk=cosEk =cos(√2γ)cos(√2λ)−coskx+cosky2sin(√2γ)sin(√2λ), (S8)

with each band being two-fold degenerate. The gap closes at when

 fk=±1,⇒√2|λ|=√2|γ|+mπ. (S9)

This gives the topological phase boundaries in the main text.

### S-1.2 Periodized evolution operator

To obtain the periodized evolution operators, one needs to compute the eigenstates of the Floquet operator. Note that all the matrices are Hermitian ones, and the eigenstates of are the same as those of the Hermitian matrix : , where

 M=1sinEk4∑j=1gjΓj=(0V†V0), V=1sinEk(g4σ0+i(g1σ1+g2σ2+g3σ3)). (S10)

We have put in the factor such that the matrix is an unitary one . The eigenstates then can be easily constructed as

 (S11)

with corresponding eigenvalues in Eq. (S8). The denotes two degenerate bands with the same eigenvalue.

For our purposes, we only need the return map at half period . So the periodized evolution operator is only computed at this point. In Eq. (S8), we choose . Thus, when the branch cut , the two eigenvalues of the return map are and . If the branch cut is , we have eigenvalues and . Thus,

 ε=0: e−iH(0)eff/2=e−iEk/2(|Ek+↑⟩⟨Ek+↑|+|Ek+↓⟩⟨Ek+↓|)−eiEk/2(|Ek−↑⟩⟨Ek−↑|+|Ek−↓⟩⟨Ek−↓|), ε=π: e−iH(π)eff/2=e−iEk/2(|Ek+↑⟩⟨Ek+↑|+|Ek+↓⟩⟨Ek+↓|)+eiEk/2(|Ek−↑⟩⟨Ek−↑|+|Ek−↓⟩⟨Ek−↓|).

Here

 ⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩|Ek+↑⟩⟨Ek+↑|+|Ek+↓⟩⟨Ek+↓|=12(I+(0Q†Q0)),|Ek−↑⟩⟨Ek−↑|+|Ek−↓⟩⟨Ek−↓|=12(I−(0Q†Q0)) Q=1sinEk(g4+ig3ig1+g2ig1−g2g4−ig3) (S12)

Thus, it can be simplified as

 ε=0:e−iH(0)eff/2=⎛⎜⎝−isinEk2σ0Q†cosEk2QcosEk2−isinEE2σ0⎞⎟⎠; (S13)

In summary, we can apply Eqs. (S7), (S8), (S12) and (S13) to compute the periodized evolution operator at half period :

 h1=τ1σ0−τ2σ2,h2k=coskxτ1σ0−sinkxτ2σ3−coskyτ2σ2−sinkyτ2σ1, U(T2)=(cos√2λ2−isin√2λ2h2k√2)(cos√2γ2−isin√2γ2h1√2) Uε=0(T2)=U(T2)e−iH(0)eff/2,Uε=π(T2)=U(T2)e−iH(π)eff/2. (S14)

## Appendix S-2 Numerical algorithms

### S-2.1 Static nested polarization

To be self-contained, we briefly review the procedure to compute static nested polarization in previous literature. Consider a static Hamiltonian or a time-independent Floquet operator , where is the period of a drive (a fixed number). The construction of topological invariants is based on their eigenstates

 H|Ek⟩=Ek|Ek⟩, or UF|Ek⟩=eiEk|Ek⟩. (S15)

Consider a 4-band model, where two gapped bands (with energy/quasi-energy) are doubly degenerate respectively, and the eigenstates can be denoted as for upper bands, and similarly for lower bands. For a system with filled lower two bands, the first order Wilson loop can be obtained as

 [Fxk]mn=⟨−Ekm|−Ek+Δkxex,n⟩, Wx,k=FxkFx,k+Δkxex⋯Fx,k+(Lx−1)Δkxex, Δkx=2πLx, (S16)

where is the lattice sites along and denotes the two filled bands. Then, one can obtain the eigenvalues and eigenstates of the Wilson loop

 Wx,k|νj(k)⟩=e2πνj(ky)|νj(k)⟩, (S17)

where denotes two Wannier bands for edge wave functions. The edge Wannier wave function is constructed by

 |wj(k)⟩=∑m[νj(k)]m|−Ekm⟩, (S18)

where denotes the -th element of the 2-component spinor . Finally, the nested polarization is obtained through