Generating a second-order topological insulator with multiple corner states by periodic driving

TL;DR

This paper explores how periodic driving can induce second-order topological insulator phases with multiple corner states in a modified BHZ model, revealing new topological phenomena and invariant characterizations.

Contribution

It demonstrates the emergence of multiple corner states through periodic driving in a modified BHZ model, expanding understanding of Floquet topological phases.

Findings

Periodic driving induces multiple corner states.

Corner states can be labeled by eigenvalues of a specific operator.

Topological invariants like Chern and winding numbers explain phase transitions.

Abstract

We study the effects of periodic driving on a variant of the Bernevig-Hughes-Zhang (BHZ) model defined on a square lattice. In the absence of driving, the model has both topological and nontopological phases depending on the different parameter values. We also study the anisotropic BHZ model and show that, unlike the isotropic model, it has a nontopological phase which has states localized on only two of the four edges of a finite-sized square. When an appropriate term is added, the edge states get gapped and gapless states appear at the four corners of a square; we have shown that these corner states can be labeled by the eigenvalues of a certain operator. When the system is driven periodically by a sequence of two pulses, multiple corner states may appear depending on the driving frequency and other parameters. We discuss to what extent the system can be characterized by topological invariants such as the Chern number and a diagonal winding number. We have shown that the locations of the jumps in these invariants can be understood in terms of the Floquet operator at both the time-reversal invariant momenta and other momenta which have no special symmetries.

Equations

$$H(\mathbf{k})$$

$$E(\mathbf{k})$$

$$H(k_x=k_y)$$

$$a(k_x)$$

$$b(k_x)$$

$$\gamma$$

$$\mathrm{where}{\varphi }_{k_x}$$

$$H(k_x)$$

$$h_{k_x}^0$$

$$H$$

$$M(t)=\{\begin{array}{ll}{M}_1& \text{if}0<t<T/4\\ M_2& \text{if}T/4<t<3T/4\\ M_1& \text{if}3T/4<t<T\end{array}$$

$$U_F$$

$$U_F|{\psi }_j⟩=e^{-i{ϵ}_{j}T}|{\psi }_j⟩,$$

$${\mathcal{S}}_2{U}_F{\mathcal{S}}_2^{-1}=U_F^{-1}.$$

$$U_F(k_x=k_y)=\exp [i(a(k_x){\tau }^z\otimes {\sigma }^0+b(k_x){\tau }^x\otimes {\sigma }^x)],$$

$$H(\mathbf{k})$$

$$E(\mathbf{k})$$

$$U_F(\mathbf{k})=\pm {\mathbb{I}}_4.$$

$$U_F({\mathbf{k}}_{\pm })=e^{-i(M_1+M_2\pm 4t_0)(T/2){\tau }^z}.$$

$$(M_1+M_2\pm 4t_0)\frac{T}{2}=n\pi ,$$

$$U_F({\mathbf{k}}_{1,2})=e^{-i(M_1+M_2)(T/2){\tau }^z},$$

$$(M_1+M_2)\frac{T}{2}=n\pi ,$$

$$[M_1+M_2\pm 2(t_x+t_y)]\frac{T}{2}=n\pi ,$$

$$[M_1+M_2\pm 2(t_x-t_y)]\frac{T}{2}=n\pi .$$

$$U_F({\mathbf{k}}_{1,2})=U_1{U}_2.$$

$$[(M_1+t_x\cos k_x+t_y\cos k_y){)}^2$$

$$+{\Delta }_1^2({\sin }^2{k}_x+{\sin }^2{k}_y){]}^{1/2}\frac{T}{2}=n_1\pi ,$$

$$[(M_2+t_x\cos k_x+t_y\cos k_y){)}^2$$

$$+{\Delta }_1^2({\sin }^2{k}_x+{\sin }^2{k}_y){]}^{1/2}\frac{T}{2}=n_2\pi ,$$

$$t_x\cos k_x+t_y\cos k_y=-\frac{1}{2}(M_1+M_2).$$

$$\frac{\cos k_x}{t_x}=\frac{\cos k_y}{t_y}.$$

$$\cos k_x$$

Full text

I Introduction

Topological insulators (TIs) have been studied extensively for the last several years hasan ; qi . A key feature of these materials is that the bulk states are gapped but there are gapless states at the boundaries which contribute to transport and other properties at low temperatures. Further, there is a bulk-boundary correspondence, namely, the bulk bands are characterized by a topological invariant (such as the Chern number for two-dimensional TIs hasan ; qi ) which is an integer, and the number of states with a particular momentum at each of the boundaries is equal to the topological invariant. Recently, a generalization of these materials called higher-order TIs has been introduced benalcazar . For instance, a second-order TI in two dimensions is a system in which the bulk and edge states are both gapped but there are gapless states at the corners of the system peng1 ; langbehn ; song ; ezawa ; dwivedi ; miert ; schindler ; franca ; xie ; calugaru ; kang ; queiroz ; hsu ; matsugatani . An electrical circuit realization of such a system has been reported recently imhof .

Closed quantum systems driven periodically in time constitute another area that has been studied by several groups in recent years dziarmaga ; polkovnikov ; dutta ; dalessio . In particular, there has been much interest in understanding the conditions under which periodic driving can generate topological phases and boundary modes top1 ; top2 ; top3 ; top4 ; top5 ; top6 ; top12 ; top15 ; top16 ; top17 ; top18 ; top19 ; top21 ; top22 ; top24 ; top25 ; top26 . It is sometimes found that even if the time-dependent Hamiltonian lies in a nontopological phase at each instant of time, the unitary time-evolution operator for one time period (called the Floquet operator) has eigenstates which are localized near the boundaries of the system. It is therefore interesting to investigate whether higher-order TIs can also be generated by periodic driving, for example, whether such a driving can generate corner states in a two-dimensional system which has no such states in the absence of driving.

There have been some earlier studies of the generation of higher-order TIs by periodic driving bomantara ; huang ; peng2 ; rodriguez2 . Using a mirror symmetry which is present in these models, a mirror-graded winding number was derived which can predict the number of corner or hinge states which appear as a result of the driving.

In contrast to the above studies, we will study a variant of the well-known Bernevig-Hughes-Zhang (BHZ) model which, we will show, hosts corner states for certain values of the system parameters. We will then study what happens when one of the parameters is varied periodically in time and show that this can generate corner states. We will then consider an anisotropic version of the BHZ model which has not been studied before to the best of our knowledge. Unlike the isotropic model, this has a nontopological phase in which only two out of the four edges of a finite-sized square has edge states which are not topologically protected. We find that this model can also have corner states, with or without driving. Interestingly, we find that driving in certain parameter regimes can generate more than one state at each corner, unlike the time-independent model which never has more than one localized state at each corner. Following earlier papers schindler ; imhof ; top22 ; asboth ; yao ; rodriguez ; bomantara ; huang ; peng2 ; rodriguez2 , we will also study two topological invariants, a Chern number and a winding number, to see whether these can be used to understand the edge and corner states which appear in the different phases of our model. Interestingly, we find that the winding number can be defined for both the isotropic and anisotropic models, even though only the isotropic model has a mirror symmetry.

The plan of this paper is as follows. In Sec. II we present the Hamiltonian of an isotropic model constructed by adding an appropriate term to the BHZ Hamiltonian; this additional term is essential to generate corner states schindler . After analyzing the different symmetries of the Hamiltonian, we will discuss how two topological invariants called the Chern number hasan ; qi and a ‘diagonal’ winding number schindler ; imhof ; rodriguez2 can be used to predict when gapless edge corner states will appear. We then discuss how to numerically study edge states by looking at a ribbon which is infinitely long but has a finite width and corner states by looking at a finite-sized square. We have shown that the corner states can be labeled by the eigenvalues ($\pm 1$) of a certain operator. In Sec. III, we study the effects of periodic driving on the isotropic model. To this end, we numerically calculate the Floquet operator and find its eigenvalues and eigenstates. We show how this modifies the regions of nonzero Chern and diagonal winding numbers; depending on the parameters, the magnitudes of the Chern number and diagonal winding number can be larger than 1. We also find that the driving can generate corner states; these states always have Floquet eigenvalues equal to $\pm 1$. In Sec. IV, we study an anisotropic model in which the hopping parameter has different values in the $x$ and $y$ directions. This model has a nontopological phase where the Chern number is zero but the winding number is nonzero; in this phase, there are edge states on only two of the edges of a finite-sized square which are not topologically protected. We again study the effects of periodic driving and show that this can generate corner states. We find that there may be more than one state at each corner if the driving frequency is low. In Sec. V we summarize our results and point out some directions for future research. Finally, we discuss in Appendix A how the locations as well as the magnitudes of the jumps in the Chern number and diagonal winding number for the periodically driven system can be understood in terms of the contributions from the time-reversal invariant momenta, and in Appendix B how the jumps in the Chern number can occur due to contributions from other momenta where there are no special symmetries.

II Time-independent isotropic model

In this section, we will study the properties of a variant of the BHZ model which is motivated by a three-dimensional model discussed in Ref. schindler, . The purpose of this is to contrast the properties of this system with that of the periodically driven system that we will study in Sec. III.

II.1 Bulk Hamiltonian

We will consider the following Hamiltonian for a system with periodic boundary conditions in which the momentum ${\bf k}=(k_{x},k_{y})$ is a good quantum number schindler

[TABLE]

where $\vec{\tau}$ and $\vec{\sigma}$ are Pauli matrices acting on the orbital and spin degrees of freedom respectively, and $\tau^{0}$ and $\sigma^{0}$ denote $2\times 2$ identity matrices in these two spaces respectively. In Eq. (1), $t_{0}$ denotes a spin-independent but orbital-dependent hopping amplitude between nearest-neighbor sites, and $\Delta_{1}$ denotes a spin-orbit coupling. We will find later that the last term with coefficient $\Delta_{2}$ is necessary in order to have corner states, with or without periodic driving. This term corresponds to orbital currents which break time-reversal symmetry oppositely in the $x$ and $y$ directions; physical systems where such a term can appear are described in Ref. schindler, .

We first study the symmetries of the Hamiltonian in Eq. (1) for the case with $\Delta_{2}=0$.

Time-reversal $\mathcal{T}$: For $\Delta_{2}=0$, we have the BHZ model bhz , a well-known example of a two-dimensional TI. The Hamiltonian is invariant under the time-reversal transformation, i.e., $\mathcal{T}H({\bf k})\mathcal{T}^{-1}=H(-{\bf k})$, where $\mathcal{T}=\tau^{0}\sigma^{y}\mathcal{K}$, and $\mathcal{K}$ denotes the complex conjugation operator. 2. 2.

Charge conjugation $\mathcal{C}$: $H({\bf k})$ has a particle-hole or charge conjugation symmetry given by $\mathcal{C}=\tau^{y}\sigma^{y}\mathcal{K}$ such that $\mathcal{C}H({\bf k})\mathcal{C}^{-1}=-H^{*}(-{\bf k})$, and a chiral symmetry $\mathcal{S}_{1}=\mathcal{CT}=\tau^{y}\sigma^{0}$ with $\mathcal{S}_{1}H({\bf k})\mathcal{S}_{1}^{-1}=-H({\bf k})$. Therefore, the system belongs to the BDI class of the Altland-Zirnbauer classification altland . There is another symmetry operator $\mathcal{S}_{2}=\tau^{x}\sigma^{z}$ which gives $\mathcal{S}_{2}H({\bf k})\mathcal{S}_{2}^{-1}=-H({\bf k})$. 3. 3.

Fourfold rotation $\mathcal{C}_{4}$: $H({\bf k})$ has a fourfold rotation symmetry about the $z$ axis, i.e., $\mathcal{C}_{4}H({\bf k})(\mathcal{C}_{4})^{-1}=H(\mathcal{C}_{4}{\bf k})$, where $\mathcal{C}_{4}=\tau^{0}e^{-i(\pi/4)\sigma^{z}}$ and $\mathcal{C}_{4}(k_{x},k_{y})=(k_{y},-k_{x})$. 4. 4.

The operator $\mathcal{P}_{1}=\mathcal{S}_{1}\mathcal{S}_{2}\tau^{z}\sigma^{z}$ commutes with the Hamiltonian. We note that $\mathcal{P}_{1}^{2}=\mathbb{I}_{4}$, the $4\times 4$ identity matrix. We will show later (see the discussion after Eq. (2)) that there is a twofold degeneracy of the energy spectrum for each value of ${\bf k}$; the two degenerate states have eigenvalues of $\mathcal{P}_{1}=\pm 1$.

We now consider the spectrum $E({\bf k})$ obtained from Eq. (1) for $\Delta_{2}=0$. Since there is a twofold degeneracy, for each value of ${\bf k}$, there are only two distinct energy bands corresponding to positive and negative energies; the exact expressions for $E({\bf k})$ will be presented below. (We will present numerical results only for the negative-energy bands, but with $\mathcal{P}_{1}$ equal to both $+1$ and $-1$. The states in these bands can be found by projecting with the operators $(\mathbb{I}_{4}\pm\mathcal{P}_{1})/2$ respectively). The gap between the two bands is shown in Fig. 1(a) as a surface plot. For a nonzero value of $\Delta_{1}$, the spectrum is found to be gapless for $M=0,\pm 2$. The region $|M|>2$ is topologically trivial and the Chern number is zero as shown in Fig. 1(b). (The Chern number is computed using the method described in Ref. fukui, ; the Berry curvature is calculated as a function of $(k_{x},k_{y})$ and the Chern number is then obtained by integrating over the Brillouin zone. We note that the bands with $\mathcal{P}_{1}=\pm 1$ have opposite values of the Chern numbers). The regions $-2<M<0$ and $0<M<2$, labeled as $A$ and $B$ in the figure, are both topological. The Chern numbers in these two regions are nonzero and have opposite signs. By studying the Hamiltonian for an infinitely long ribbon (in Sec. II.2), we find there there are robust one-dimensional gapless edge states in the topological regime as shown in Fig. 2(a).

Now switching on the $\Delta_{2}$ term, we find that the symmetries of the Hamiltonian get reduced as follows.

For $\Delta_{2}\neq 0$, the time-reversal symmetry $\mathcal{T}$ as well as the rotational symmetry $\mathcal{C}_{4}$ are broken. However the symmetry given by the product of the two, i.e., $\mathcal{C}_{4}T$ is preserved. This means that $(\mathcal{C}_{4}\mathcal{T})H({\bf k})(\mathcal{C}_{4}\mathcal{T})^{-1}=H(\mathcal{C}_{4}T{\bf k})$, where $\mathcal{C}_{4}T(k_{x},k_{y})=(-k_{y},k_{x})$. 2. 2.

The charge conjugation and one of the other symmetries act as before: $\mathcal{C}H({\bf k})\mathcal{C}^{-1}=-H^{*}(-{\bf k})$ and $\mathcal{S}_{2}H({\bf k})\mathcal{S}_{2}=-H({\bf k})$. The operator $\mathcal{S}_{1}$ does not describe a symmetry of the system.

For $\Delta_{2}\neq 0$, we cannot use $\mathcal{P}_{1}$ to distinguish between the two degenerate bands since $\mathcal{P}_{1}$ does not commute with the Hamiltonian. It is therefore not possible to calculate the Chern number when $\Delta_{2}\neq 0$.

We note that in the presence of a $\Delta_{2}$ term, the $\mathcal{C}_{4}T$ symmetry gaps out the edge states but gapless corner states appear for certain values of $M$. The system is then called a second-order TI.

The energy spectrum of the Hamiltonian in Eq. (1) can be found as follows. Since the four matrices $\tau^{z}\otimes\sigma^{0}$, $\tau^{x}\otimes\sigma^{x}$, $\tau^{x}\otimes\sigma^{y}$ and $\tau^{y}\otimes\sigma^{0}$ anticommute with each other and the square of each of them is equal to $\mathbb{I}_{4}$, the spectrum can be found by taking the square of $H({\bf k})$. This gives

[TABLE]

(The fact that the Hamiltonian is a $4\times 4$ traceless matrix and there are only two possible energy levels implies that each energy level must be twofold degenerate). Equation (2) implies that the bulk gap can only vanish at one of the four momenta ${\bf k}=(0,0),~{}(0,\pi),~{}(\pi,0)$ and $(\pi,\pi)$; further the energy at one of these momenta is zero only if (i) $M=\pm 2t_{0}$ or (ii) $M=0$ and $\Delta_{2}=0$. These parameter values give the locations of topological phase transitions as we will see below. At all other values of the parameters, the spectrum will be gapped at all momenta.

Based on the existence of mirror symmetries, certain mirror-graded topological invariants were derived in Refs. schindler, , imhof, and rodriguez2, . While our isotropic model has a similar mirror symmetry, the anisotropic model discussed in Sec. IV does not have that symmetry. Nevertheless, we will see that both models allow us to define a topological invariant called ‘diagonal’ winding number (which is motivated by the mirror-graded topological invariants mentioned above); we use the word ‘diagonal’ since this winding number is defined on one of the diagonals in the Brillouin zone as we will see. In order to derive this winding number, we observe that the Hamiltonian in Eq. (1) simplifies considerably along the diagonals $k_{x}=\pm k_{y}$. For instance, along $k_{x}=k_{y}$, the Hamiltonian takes the form

[TABLE]

which, interestingly, does not depend on $\Delta_{2}$. Next, we see that by rotating by an appropriate angle using the matrix $\sigma^{z}$, we can transform $(\sigma^{x}+\sigma^{y})/\sqrt{2}\to\sigma^{x}$. We then obtain

[TABLE]

This Hamiltonian only depends on two matrices, $\tau^{z}\otimes\sigma^{0}$ and $\tau^{x}\otimes\sigma^{x}$, which are Hermitian and square to $\mathbb{I}_{4}$. If we take the coefficients of the two matrices,

[TABLE]

to be the coordinates of a point in a two-dimensional plane, we obtain a closed curve as $k_{x}$ goes from $-\pi$ to $\pi$. Assuming that the system is gapped for all values of $(k_{x},k_{y})$ which includes the line $k_{x}=k_{y}$ as a special case, we see that $(a(k_{x}),b(k_{x}))$ will not be equal to $(0,0)$ for any value of $k_{x}$. We can therefore define the winding number $\gamma$ of the closed curve around the origin; this will be a topological invariant since it will not change if any of the parameters are changed slightly (assuming that such changes do not make the curve pass through the origin). Mathematically, we have

[TABLE]

We note that the existence of the winding number is crucially dependent on the fact that on the diagonals $k_{x}=\pm k_{y}$, the Hamiltonian essentially reduces to a sum of only two anticommuting matrices. If the Hamiltonian contained a third anticommuting matrix such as $\tau^{x}\otimes\sigma^{z}$ or $\tau^{y}\otimes\sigma^{0}$, it would not be possible to define the diagonal winding number. We also observe that the diagonal winding number does not depend on the value of $\Delta_{2}$, unlike the Chern number which can only be calculated if $\Delta_{2}=0$.

A parametric plot of $b(k_{x})$ versus $b(k_{x})$ is shown in Fig. 1(d) for the case when $M=0.1$, $t_{0}=1$, and $\Delta_{1}=1$. We can see easily that the diagonal winding number is 1. This is consistent with the lattice calculation for this set of parameters where we find one localized state at every corner.

Figure 1(c) shows the diagonal winding number for the model described in Eq. (1). Comparing this with Fig. 1(b), we see that the regions of nonzero values of the Chern number and diagonal winding number coexist in this model when $\Delta_{2}=0$.

II.2 Edge states for an infinitely long ribbon

We will first study the effect of the $\Delta_{2}$ term on the edge states in this system. To this end, we consider a strip of the material which is infinitely long in the $x$ direction and has a finite width (with $N_{y}$ sites) in the $y$ direction. This means that $k_{x}$ is a good quantum number and for each $k_{x}$, we effectively have a chain with $N_{y}$ sites which is related to its neighboring chains by factors of $\exp(\pm ik_{x})$. At each site there are four degrees of freedom corresponding to the two orbitals ($d$ and $f$ which we denote by $\tau^{z}=\pm 1$) and two spins ($\uparrow$ and $\downarrow$ denoted by $\sigma^{z}=\pm 1$). The creation operator for the $d$ orbital at the $n_{y}$-th site of the chain for spin $s$ is denoted by $d^{\dagger}_{n_{y},s}$ where $s=\uparrow,\downarrow$. Similarly, for the $f$ orbital, the corresponding operator is $f^{\dagger}_{n_{y},s}$. Hence the spinor for each value of $k_{x}$ has $4N_{y}$ components and the Hamiltonian $H(k_{x})$ is a $4N_{y}\times 4N_{y}$ matrix. By keeping the $k_{x}$ part of the Hamiltonian as it is, but discretizing in the $y$ direction, we obtain the following Hamiltonian,

[TABLE]

where $\mathbb{I}_{N_{y}}$ is the $N_{y}\times N_{y}$ identity matrix.

In our numerical calculations, we will set $t_{0}=1$ and express all the other parameters in units of $t_{0}$. Diagonalizing Eq. (7a), we find $4N_{y}$ energy eigenvalues and eigenvectors for each momentum $k_{x}$, such that each eigenvector is a $(4N_{y})$-component spinor. The energy levels as a function of momentum $k_{x}$ are shown in Fig. 2 for both $\Delta_{2}=0$ and $\Delta_{2}\neq 0$. We have taken $M=1$ and $\Delta_{1}=0.5$.

The bulk states form the two continuous energy bands shown in blue. These two bands are separated by an energy gap whose value depends on both $M$ and $\Delta_{1}$ consistent with Fig. 1(a). Let us first consider Fig. 2(a) where we have taken $\Delta_{2}=0$. Here the bulk gap hosts gapless edge states (red solid lines) which go from one band to the other and are related to the Chern number of the infinite system by the bulk-boundary correspondence; we thus have a two-dimensional TI. On switching on a $\Delta_{2}$ term, these edge states become gapped as shown in Fig. 2(b). Since each of these edge states lie within the same band, i.e., they do not connect the two bands, they are topologically trivial. This is consistent with the fact the Chern numbers are zero for the infinite system with $\Delta_{2}\neq 0$. A useful topological invariant for $\Delta_{2}\neq 0$ is the diagonal winding number $\gamma$ which is explained in Sec. II.1. However, the bulk-boundary correspondence for this case has to be analyzed using a system which is finite in both the directions, such as a finite-sized square lattice, as discussed in the next section.

II.3 Corner states for a finite-sized square lattice

In order to find the corner states of the system, we consider a finite-sized square lattice lying in the $x-y$ plane with $N_{x}$ and $N_{y}$ sites along the $x$ and $y$ directions respectively. The total number of lattice points is $N=N_{x}\times N_{y}$. Since there are two orbital ($d$ and $f$) and two spin ($\uparrow$ and $\downarrow$) degrees of freedom at each site, we arrive at a $4N\times 4N$ Hamiltonian,

[TABLE]

where ${\langle n,n^{\prime}\rangle}$ denotes nearest neighbors, $\epsilon_{n,n^{\prime}}=i$ or 1 and $\xi_{n,n^{\prime}}=\pm 1$ for nearest neighbors along the $x$ and $y$ directions respectively.

The different terms in Eq. (8) can be understood as follows. $M$ acts as a staggered chemical potential for the two orbitals. $t_{0}$ is the amplitude for nearest-neighbor hopping that keeps both the spin and orbital the same but differs in sign for the two orbitals. $\Delta_{1}$ flips both the spin and orbital degrees of freedom and also depends on the direction of hopping via $\epsilon_{n,n^{\prime}}$. $\Delta_{2}$ describes a hopping that flips the orbital but keeps the spin the same, and this term depends on the direction of hopping through $\xi_{n,n^{\prime}}$. Diagonalizing this Hamiltonian gives the $4N$ energy eigenvalues and eigenvectors, all of which turn out to be doubly degenerate.

For our numerical calculations, we have considered a square lattice with $25$ sites in each direction, i.e., $N_{x}=N_{y}=25$. The parameter $t_{0}$ is set to unity and all other parameters and the energy are expressed in units of $t_{0}$. The results for $M=1$, $\Delta_{1}=1$ and $\Delta_{2}=0.1$ are shown in Fig. 3. The plot of energy eigenvalues versus the eigenvalue index is shown in Fig. 3(a). The edge states (red) and corner states (blue) are clearly separated in energy from the bulk states (black). As is clear from the inset in this figure, there are four states which are very close to zero energy. They become degenerate in the thermodynamic limit; these fourfold degenerate states are found to exist only at the corners of the square. One such state which lives at the corner labeled $3$ is shown in Fig. 3(b). Similar states at zero energy exist at each of the four corners. We find that the decay length of these states is much larger along the edges than along the diagonal direction into the bulk. This can be understood as follows. Typically, the decay length of a boundary state is inversely proportional to the gap of the corresponding bulk states. In our system, the bulk gap is much larger than the gap of the edge states (the latter is proportional to $\Delta_{2}$); we can clearly see this difference in Fig. 3(a). Hence the decay length into the bulk is much smaller than the decay length along the edges.

We have found that the four-component spinors corresponding to these corner states have an interesting structure. We recall that the Hamiltonian in Eq. (1) changes sign under a transformation by $\mathcal{S}_{2}=\tau^{x}\sigma^{z}$. We therefore expect the space of zero-energy states to remain invariant under the action of $\mathcal{S}_{2}$; in particular, the corner states should be eigenstates of $\mathcal{S}_{2}$. We find that the states localized at corners labeled $1$ and $3$ have eigenvalues of $\mathcal{S}_{2}$ equal to $+1$ while those at corners $2$ and $4$ have eigenvalues $-1$. We can understand the relation between the eigenvalues of $\mathcal{S}_{2}$ for the states at the different corners as follows. First, Eq. (1) implies that a rotation by $\pi$ about the $z$ axis (which transforms $x\to-x$ and $y\to-y$) can be performed by the operator $\sigma^{z}$. This transformation interchanges the corners (1,3) and corners (2,4). Since $\sigma^{z}$ commutes with $\mathcal{S}_{2}$, this implies that the eigenvalues of $\mathcal{S}_{2}$ for the states at corners 1 and 3 must be identical, and similarly for the eigenvalues for the states at corners 2 and 4. Second, we see from Eq. (1) that the reflection $x\to-x$ and $y\to y$ can be performed by the operator $\sigma^{y}$. This interchanges the corners (1,2) and corners (3,4). Since $\sigma^{y}$ anticommutes with $\mathcal{S}_{2}$, this implies that the eigenvalues of $\mathcal{S}_{2}$ for the states at corners 1 and 2 must have opposite signs, and similarly for the eigenvalues for the states at corners 3 and 4.

The edge states, whose energies are shown in red in Fig. 3(a), are found to be localized along all the edges of the system. However, these are not protected topologically as explained in the analysis of the ribbon geometry in Sec. II.2 (they do not go from one bulk band to the other).

To recapitulate, in this section we have analyzed a topologically nontrivial two-dimensional system which, in certain parameter regimes is found to have robust zero-energy “corner states”. These can be connected to the bulk system using a topological invariant called the diagonal winding number.

It is useful to understand why the $\Delta_{2}$ term in Eq. (1) gaps out the edge states and produces zero-energy corner states. Basically, the $\Delta_{2}$ term provides a mass term for the edge states and therefore gaps them out. However, the mass term has opposite signs for edges running along the $x$ and $y$ directions, due to the $\cos k_{x}-\cos k_{y}$ structure. Since the mass terms changes sign on going around a corner from an edge along $x$ to an edge along $y$, a zero-energy state appears at each corner schindler . This mechanism for the appearance of zero-energy corner states is different from the one discussed in Ref. matsugatani, . That work considers a model of a three-dimensional TI; the addition of a weak magnetic field gaps out the surface states but produces gapless states which are localized along the vertical hinges. Then the coupling in the $z$ direction is smoothly turned off; this produces a set of decoupled two-dimensional TIs each of which has zero-energy corner states which are the remnants of the hinge states. We also note that our model lies in the symmetry class BDI while the one discussed in Ref. matsugatani, is in the class AIII.

III Periodically driven isotropic model

Having understood the properties of the time-independent model of a second-order TI in various parameter regimes, we now proceed to study what happens when the system is driven periodically in time. In particular, we will study the effect of varying the parameter $M$ in the Hamiltonian. We will study the topological properties of this driven system both in momentum space (i.e., for a system with periodic boundary conditions) as well as for a finite-sized square lattice using a time-evolution operator, and we will again use the Chern number and diagonal winding number to characterize the system bomantara ; huang ; rodriguez2 ; peng2 .

First, let us consider the bulk system, i.e., we work with the momentum space Hamiltonian $H$ given in Eq. (1) with the parameter $M$ varying in time as,

[TABLE]

within a single time period $[0,T]$; we then continue this periodically by taking $M(t+T)=M(t)$ for all $t$. We are interested in studying the system stroboscopically, i.e., at times $t=\mathcal{N}T$ where $\mathcal{N}$ runs over all integers. We will calculate the quasienergy spectrum by numerically diagonalizing the Floquet operator $U_{F}$ (defined in Eq. (10) below) for each momentum ${\bf k}$. We find that each of the two quasienergy bands is twofold degenerate.

The time-evolution operator is given by

[TABLE]

where $T$ is the time period of the driving and $H_{1}$ and $H_{2}$ are the Hamiltonians with $M=M_{1}$ and $M_{2}$ respectively. (The symbol $\mathscr{T}_{t}$ denotes the time-ordered product). Since we will study both momentum and real space systems, the Hamiltonian takes two different forms for these two cases, as will be explained below. Further, since $U_{F}$ is a unitary operator, its eigenvalues are complex numbers with unit magnitude, i.e.,

[TABLE]

where $\epsilon_{j}$’s are the quasienergies and are defined modulo $2\pi/T$. We define the first Floquet Zone such that $\epsilon_{j}\in[-\pi/T,\pi/T]$ and $|\psi_{j}\rangle$’s are the corresponding Floquet eigenstates.

The driving protocol described in Eq. (9) has been chosen to satisfy a particular symmetry of the Floquet operator $U_{F}$. We saw earlier that $\mathcal{S}_{2}=\tau^{x}\sigma^{z}$ satisfies $\mathcal{S}_{2}H\mathcal{S}_{2}^{-1}=-H$; this is true in both momentum space and real space. Equation (10) then implies that

[TABLE]

This implies that if $|\psi_{j}\rangle$ is an eigenstate of $U_{F}$ with eigenvalue $e^{-i\epsilon_{j}T}$, then $\mathcal{S}_{2}|\psi_{j}\rangle$ is an eigenstate of $U_{F}$ with eigenvalue $e^{+i\epsilon_{j}T}$. Thus the eigenvalues of $U_{F}$ must appear in complex conjugate pairs. Next, if there is an eigenstate $|\psi_{j}\rangle$ with eigenvalue $e^{-i\epsilon_{j}T}=\pm 1$, $\mathcal{S}_{2}|\psi_{j}\rangle$ must be the same as $|\psi_{j}\rangle$. Hence $|\psi_{j}\rangle$ must be an eigenstate of $\mathcal{S}_{2}$, and the eigenvalue must be $\pm 1$ since $(\mathcal{S}_{2})^{2}$ is the identity operator. We will see below that corner states always appear with $e^{-i\epsilon_{j}T}=\pm 1$; hence they must also be eigenstates of $\mathcal{S}_{2}$ with eigenvalue $\pm 1$. Further, the arguments given in Sec. II.3 regarding the relative eigenvalues of $\mathcal{S}_{2}$ at the four corners are valid for the periodically driven system as well. Namely, a rotation by $\pi$ about the $z$ axis is performed by $\sigma^{z}$ while a reflection of the $x$ coordinate is performed by $\sigma^{y}$, and these commute and anticommute respectively with $\mathcal{S}_{2}$. Combining these together, we see that if corner 1 has $(n_{+},n_{-})$ states with eigenvalues of $\mathcal{S}_{2}$ equal to $(+1,-1)$ respectively, corner 3 will have the same number of states with those eigenvalues, while corners 2 and 4 will have the number of eigenvalues of the two kinds interchanged to $(n_{-},n_{+})$.

Having computed the Floquet operator $U_{F}({\bf k})$, we investigate if there is any correspondence between the bulk states and the corner states that appear due to the driving. To choose an appropriate time period for studying the Floquet problem, we have to compare the time period $T$ with an intrinsic timescale of the system, which, in this case can be taken to be $1/t_{0}$. For fast driving, i.e., when the time period $T\ll 1/t_{0}$, the system does not have time to respond to the changing Hamiltonian and the properties are not very different from the static system. On the other hand, the system has interesting behavior for intermediate driving frequency, i.e., when $T$ is comparable to $1/t_{0}$.

To study the intermediate frequency regime, we fix $M_{1}=2.5$ and $M_{2}=3.5$ and numerically find the eigenvalues and eigenvectors of $U_{F}({\bf k})$ for different values of $T\approx 1/t_{0}$. We then use these eigenvectors to compute the topological invariants. The results for this are shown in Fig. 4. The two topological invariants i.e., the Chern number and diagonal winding number are shown for the cases $\Delta_{2}=0$ and $\Delta_{2}=0.1$ in Figs. 4(a) and 4(b) respectively.

In Fig. 4(a) we have shown the Chern number for the case $\Delta_{2}=0$. The locations and magnitudes of the jumps in the Chern number can be understood by studying the time-evolution operator at some special points in the Brillouin zone. The Chern number changes abruptly when $U_{F}({\bf k})$ becomes equal to $\mathbb{I}_{4}$ at one of the four time-reversal invariant momenta $(0,0)$, $(0,\pi)$, $(\pi,0)$ and $(\pi,\pi)$. The magnitude and sign of the jump is determined by the momentum point where the time-evolution operator becomes $\mathbb{I}_{4}$. A detailed explanation of this is given in Appendix A.

We can also define the diagonal winding number for the periodically driven system. Along one of the diagonals, say, $k_{x}=k_{y}$, we saw in Sec. II.1 that, after an appropriate transformation, the Hamiltonian becomes a linear combination of two matrices, $\tau^{z}\otimes\sigma^{0}$ and $\tau^{x}\otimes\sigma^{x}$, at any time $t$. The commutator of these two matrices gives a third matrix, $\tau^{y}\otimes\sigma^{x}$, such that the three matrices form a closed Lie algebra. The Floquet operator $U_{F}(k_{x}=k_{y})$ defined in Eq. (10) must therefore be an exponential of a linear combination of these three matrices. Next, we observe that $\tau^{z}\otimes\sigma^{0}$ and $\tau^{x}\otimes\sigma^{x}$ are symmetric matrices while $\tau^{y}\otimes\sigma^{x}$ is antisymmetric. Since $H_{1}$ and $H_{2}$ are symmetric, $U_{F}(k_{x}=k_{y})$ given by Eq. (10) must also be a symmetric matrix. Hence it must be the exponential of a linear combination of only $\tau^{z}\otimes\sigma^{0}$ and $\tau^{x}\otimes\sigma^{x}$, and not $\tau^{y}\otimes\sigma^{x}$. Thus we can write

[TABLE]

where $(a(k_{x}),b(k_{x}))$ are fixed uniquely by demanding that $0<\sqrt{(a(k_{x}))^{2}+(b(k_{x}))^{2}}<\pi$; we can demand this if $U_{F}(k_{x}=k_{y})$ is not equal to $\pm\mathbb{I}_{4}$ for any value of $k_{x}$ top15 . We now take $(a(k_{x}),b(k_{x}))$ to be the coordinates of a point in a two-dimensional place, thus defining a closed curve as $k_{x}$ goes from $-\pi$ to $\pi$. We then define the winding number of the closed curve around the origin as in Eq. (6).

For the finite-sized system on a $25\times 25$ lattice, we find that corner states only appear when the diagonal winding number of the bulk system is nonzero, implying that there is a nontrivial second-order bulk boundary correspondence. The time-evolution operator $U$ for the lattice model is given by the time-ordered product given in Eq. (10) with the Hamiltonian $H$ of the form given in Eq. (8) with the parameter $M$ given by Eq. (9) for the two halves of the cycle, i.e., $H_{1}=H|_{M=M_{1}}$ and $H_{2}=H|_{M=M_{2}}$. For a finite-sized square lattice, $U_{F}$ is a $4N\times 4N$ matrix where $N=N_{x}\times N_{y}$. The factor of $4$ is from the two spin and two orbital degrees of freedom at each lattice site. Diagonalizing this $U_{F}$ gives $4N$ quasienergy eigenvalues denoted by $\epsilon_{j}$ in Eq. (11) and each eigenvector $|\psi_{j}\rangle$ is a $4N$-component spinor.

Since each eigenvalue $\exp({-i\epsilon_{j}T})$ is a complex number of unit magnitude, it can be represented as a point on the unit circle. Figure 5(a) shows a plot of $\cos(\epsilon_{j}T)$ versus $\sin(\epsilon_{j}T)$. The points marked in red (shown more clearly in the inset) correspond to the eigenstates which are localized at one of the corners. One such state is shown in Fig. 5(b). These corner states are fourfold degenerate and lie close to $\sin(\epsilon_{j}T)=0$. For this particular choice of parameters, these quasienergies lie at $\cos(\epsilon_{j}T)=-1$. They may also lie at $\cos(\epsilon_{j}T)=+1$ for a different set of parameters. For any choice of parameters, as long as we are in the higher-order topological sector, i.e., the diagonal winding number is nonzero, these corner states exist and are always separated from the bulk and edge states by a finite gap which is proportional to $\Delta_{2}$.

It is interesting to note that corner states can appear even when both $M_{1}$ and $M_{2}$ are larger than 2 as we can see from Fig. 4(b). Thus the periodic driving can generate corner states even when the instantaneous Hamiltonian does not have corner states at any time; as Fig. 1(c) shows, the time-independent Hamiltonian has no corner states if $M>2$.

IV Anisotropic model

We now consider a variation of the model discussed so far by making the hopping amplitude $t_{0}$ anisotropic, i.e., taking the hoppings along $x$ and $y$ directions to be $t_{x}$ and $t_{y}$ respectively. In the case of the static system, the Hamiltonian in Eq. (1) becomes

[TABLE]

(In the anisotropic model, the coefficients of $\cos k_{x}$ and $\cos k_{y}$ in the last term in Eq. (14) could, in principle, be different, namely, $\Delta_{2x}\neq\Delta_{2y}$. However, we will assume that they are equal, and we will call the coefficient $\Delta_{2}$ as before. This will allow us to define a diagonal winding number as before). This Hamiltonian has the same symmetries as the one in Eq. (1) except for the ones which involve the $\mathcal{C}_{4}$ transformation. Further, unlike the isotropic model, the energy spectrum does not remain the same if we interchange $k_{x}$ and $k_{y}$.

We can find the energy spectrum of Eq. (14) in the same way as given in Eq. (2). We find that

[TABLE]

This implies that the bulk gap can vanish only if (i) $M=\pm(t_{x}+t_{y})$ or (ii) $M=\pm(t_{x}-t_{y})$ and $\Delta_{2}=0$. These give the locations of topological phase transitions.

We can now follow the same procedure as outlined in Sec. II.1 to find the topological invariants, i.e., the Chern number and the diagonal winding number, corresponding to the cases with $\Delta_{2}=0$ and $\Delta_{2}\neq 0$ respectively. Note that the transformations given in Eqs. (3-5) which lead to a point in a two-dimensional plane and hence to a winding number continue to work in this anisotropic system, as long as the system is gapped at all momenta. The phase diagram and the values of the topological invariants are shown in Fig. 6.

In the isotropic case, as discussed in Sec. II.1, at $M=0$ there is a transition from one topologically nontrivial phase to another (regions labeled $A$ and $B$ in Fig. 1(b)). However, on introducing an anisotropy by taking $t_{x}\neq t_{y}$, the two topologically nontrivial regions are separated by an intermediate region where the Chern number is zero. (The boundaries of all the regions can be found by finding the values of $M$ where the energy eigenvalues of Eq. (14), with $\Delta_{2}=0$, vanish at one the four time-reversal invariant momenta). For the case $t_{x}=1$ and $t_{y}=2$, we see in Fig. 6(a) that the intermediate phase, labeled as $II$, lies in the region $|M|\leq 1$, whereas the phases labeled $I$ and $III$ are topologically nontrivial with Chern numbers $C=\pm 1$ respectively. The width $W_{M}$ of phase $II$ depends on the hoppings as $W_{M}=2||t_{x}|-|t_{y}||$.

On studying the Hamiltonian on a square lattice with $t_{x}\neq t_{y}$, we find that the edge states in the topological and nontopological phases behave very differently from each other. In regions $I$ and $III$, the edge states are present along all the edges of the system and are topologically protected. However, in region $II$, the edge states exist only on the edges parallel to the $x$ direction (when $|t_{x}|<|t_{y}|$) and are not topologically protected, even though they are gapless (if $\Delta_{2}=0$) and they lie inside the bulk gap. This is consistent with the value of the Chern number which is zero in region $II$. If $|t_{y}|<|t_{x}|$, the edge states are found only on the edges parallel to the $y$ direction. The system can therefore be described as a weak topological insulator in region $II$ for $\Delta_{2}=0$ hasan ; qi ; fu .

For $\Delta_{2}=0$, the different phases of the system can be distinguished from each other by the Chern number as shown in Fig. 6(a). If $\Delta_{2}\neq 0$, the Chern number is ill-defined in all the phases, but they can be distinguished from each other by the diagonal winding number (whose value does not depend on $\Delta_{2}$). Figure 6(b) shows the winding number $\gamma$ for the anisotropic system with $\Delta_{2}=0.3$. We see that $\gamma=1$ for $|M|\leq 3$. The winding number does not differentiate among the phases $I$, $II$ and $III$ which are topologically different from one another as shown by the Chern number. We thus see that the regions of nonzero Chern and diagonal winding numbers are not identical in the anisotropic model.

The corner states are obtained by diagonalizing the Hamiltonian in real space on a square lattice. The expression for this Hamiltonian is similar to Eq. (8) with the hoppings along $x$ and $y$ directions being $t_{x}$ and $t_{y}$ respectively. Figure 7 shows the edge states in regions $I$ and $II$ for $\Delta_{2}=0$. We see that in region $I$, the edge states lie on all the edges while in region $II$, they only lie on the edges parallel to the $x$ direction (for $|t_{x}|<|t_{y}|$). This is consistent with the Chern numbers in these two regions. Figure 8 shows the energy eigenvalues (the corner states lie at zero energy) and a corner state in region $II$ for $\Delta_{2}=0.3$.

We now study the effect of periodically driving the system by varying the parameter $M$ between two values $M_{1}$ and $M_{2}$ both of which lie in the nontopological regime, i.e., in region $II$ of Fig. 6. The values of $T$ where topological transitions occur depend on the values of $M_{1}$ and $M_{2}$ chosen in Eq. (9).

In Fig. 9 (a), the Chern number is calculated for the anisotropic Hamiltonian in Eq. (14) for $\Delta_{2}=0$. The jumps in the Chern number at some specific values of the parameter $T$ signify topological transitions and can be understood as discussed earlier. We note that $U_{F}({\bf k})=\pm~{}\mathbb{I}_{4}$ at these points and the magnitude of the Chern number can exceed 1 unlike the static system. It turns out that the jumps in the Chern number occurring for $T<4.4$ are due to contributions from the time-reversal invariant momenta, while the jumps occurring for $T>4.4$ are due to contributions from other momenta which have no special symmetries. This is explained in Appendices A and B respectively. We also see some sharp fluctuations in the Chern number at $T=4.81$ and $4.91$; the reasons for these are explained in Appendix B.

We can compute the diagonal winding number $\gamma$ from the wave functions of the quasienergy states of the Floquet operator; the value of $\gamma$ does not depend on $\Delta_{2}$. The winding number also jumps at some specific values of $T$ signifying second-order topological transitions.

Finally we consider a finite-sized square with $30\times 30$ sites and periodically drive the anisotropic lattice Hamiltonian. From the momentum-space calculations, we find regions with $\gamma=1$ and 3, as shown in Figs. 10(a) and 10(b). For the lattice calculations, we will work with the driving period, i.e., $T=3$, and demonstrate that the diagonal winding number 3 indeed counts the number of corner states. The corner states appear in this system only for those parameter values for which the winding number of the Hamiltonian in Eq. (14) is nonzero. For low-frequency driving, i.e., when the time period of driving is large compared to the timescale $1/t_{0}$, we sometimes find more than 4 corner states in the system which become degenerate at zero quasienergy in the thermodynamic limit. Figure 11 shows that for some values of the parameters, there are a total 12 corner states (3 at each corner), all lying at zero quasienergy in the thermodynamic limit. All the three states shown at corner 3 have the eigenvalue -1 for the operator $\mathcal{S}_{2}=\tau^{x}\sigma^{z}$.

V Discussion

We first summarize the results obtained in this work. We considered a variant of the BHZ model in which there is an additional term (with a coefficient $\Delta_{2}$) which breaks the $\mathcal{C}_{4}$ and $\mathcal{T}$ symmetries separately but is invariant under their product. We first considered an isotropic model and studied how its properties change as a parameter $M$ is varied. The system can be characterized by two topological invariants called the Chern number and diagonal winding number. The Chern number can be calculated only when the parameter $\Delta_{2}=0$, while the winding number can be calculated for any value of $\Delta_{2}$ and it does not depend on the value of $\Delta_{2}$. When $\Delta_{2}=0$, there are no corner states. If the Chern number is nonzero, there are gapless edge states; this system is then in a usual TI phase (called first-order). When $\Delta_{2}\neq 0$ and the winding number is nonzero, there are gapped edge states and gapless states which are localized at the corners of a finite-sized square lattice; the system is then in a second-order TI phase. We have shown that corner states have an interesting structure, namely, they are eigenstates of the operator $\tau^{x}\sigma^{z}$. We then studied the effect of periodic driving of the parameter $M$. We find that the driving can generate corner states in certain ranges of parameters where the corresponding time-independent Hamiltonian has no such states. We have shown that the number of corner states is given by the winding number.

Next, we have studied an anisotropic version of the model where the hopping amplitudes in the $x$ and $y$ directions are different. Once again, we can define a Chern number if $\Delta_{2}=0$ and a diagonal winding number for any value of $\Delta_{2}$. An interesting feature of this model is that when $\Delta_{2}=0$, there is a nontopological phase in which the Chern number is zero and there are states localized along only two out of the four edges of a finite-sized square. When $\Delta_{2}\neq 0$ and the winding number is nonzero, gapless corner states appear. We again find that periodic driving of $M$ in this anisotropic model can generate corner states and their number is given by the winding number.

For the periodically driven models (both isotropic and anisotropic), we sometimes find that there can be more than one state localized at each corner of the system. We have shown that the locations of the jumps in the Chern number and diagonal winding number can be understood as the points where the Floquet operator becomes proportional to the identity matrix at one of the time-reversal invariant momenta and, in the case of the Chern number, at other momenta where there are no special symmetries; these jumps signify topological transitions. Although there has been no explicit discussion in the earlier literature of momenta with no special symmetries which can lead to jumps in the Chern number, the role of such momenta has been alluded to in Refs. top25, , top26, and bomantara, .

We have shown that the diagonal winding number can also predict the number of states at each corner for the periodically driven system. A possible direction for future study would be to look for a topological invariant which can predict the number of corner states with eigenvalues of $\tau^{x}\sigma^{z}=\pm 1$ separately. For instance, an interesting topological invariant to study would be the winding number of the one-dimensional states which are localized at one particular edge of an infinitely long ribbon with a finite width; these states would be labeled by a momentum $k_{x}$ (as shown in Fig. 2), and an expression similar to Eq. (6) could then be used to calculate the winding number. At first sight, such a numerical calculation appears challenging because a finite width allows hybridization between the states at the opposite edges and we would therefore generally find states which are superpositions of states localized at the two edges. However, if these can be projected to one of the edges in a way which varies smoothly with $k_{x}$, it may be possible to calculate the corresponding winding number. This may provide a way of understanding the number of corner states.

Acknowledgments

R.S. and A.D. thank Adhip Agarwala for useful discussions. A.D. acknowledges funding from SERB, DST, India for NPDF Research Grant No. PDF/2016/001482. D.S. thanks DST, India for Project No. SR/S2/JCB-44/2010 for financial support.

Appendix A Contributions to Chern and diagonal winding numbers from

time-reversal invariant momenta

The jumps in the Chern number shown in Figs. 4(a) and 9 (a) can be understood by studying the time-evolution operators at some specific points in the Brillouin zone. These are the points where

[TABLE]

First, let us study the anisotropic model and $\Delta_{2}=0$. Consider the Hamiltonian at the momenta ${\bf k}_{+}=(0,0)$ and ${\bf k}_{-}=(\pi,\pi)$. At these two momenta the $\Delta_{1}$ and $\Delta_{2}$ terms vanish. Therefore the Hamiltonians in the two halves of the cycle, $H_{1}$ and $H_{2}$ commute with each other and the time-evolution operator given in Eq. (10) can be simplified to

[TABLE]

This gives the condition in Eq. (16) if

[TABLE]

where $n$ is an integer. For all our calculations in the isotropic case we have set $t_{0}=1$, $M_{1}=2.5$ and $M_{2}=3.5$. For ${\bf k}_{+}$, the values of $T$ obtained from Eq. (18) are given by $n\pi/5$ Similarly, for ${\bf k}_{-}$, $T$ is found to be $n\pi$. These are precisely the points in Fig. 4(a) where the Chern number jumps; all the jumps are by $\pm 1$ since there is a contribution from either ${\bf k}_{+}$ or ${\bf k}_{-}$ but not both.

Similarly, the points ${\bf k}_{1}=(0,\pi)$ and ${\bf k}_{2}=(\pi,0)$ also contribute. In this case Eq. (10) reduces to

[TABLE]

which satisfies the condition in Eq. (16) if

[TABLE]

where $n$ is an integer. This implies that there should be jumps in the Chern number at $T=n\pi/3$. Note that the contributions at these values of $T$ come from both ${\bf k}_{1}$ and ${\bf k}_{2}$. This explains the jumps of $\pm 2$ in the Chern number at these values of $T$ in Fig. 4(a).

A similar analysis can be carried out for the anisotropic model where $t_{x}\neq t_{y}$. We then find that Eq. (18) gets modified to

[TABLE]

while Eq. (20) changes to

[TABLE]

Equations (21-22) explain the locations of the jumps in the Chern number in Fig. 9 (a), where we have set $M_{1}=-0.9$, $M_{2}=-0.45$, $t_{x}=1$ and $t_{y}=2$.

We now consider the case when $\Delta_{2}\neq 0$. The behavior at the momenta ${\bf k}_{\pm}$ is the same as in the case $\Delta_{2}=0$, because at these points the $\Delta_{2}$ term vanishes, and the Hamiltonians for the two halves of the cycle continue to commute. However, this is not true at the momenta ${\bf k}_{1,2}$ since the $\Delta_{2}$ term survives at these momenta. The time-evolution operator must therefore be written as

[TABLE]

We would again like this to be equal to $\pm\mathbb{I}_{4}$ as in Eq. (16). However, this would imply that $U_{1}$ and $U_{2}$ must be inverses of each other (possibly up to a sign) and hence must commute. This contradicts the fact that $U_{1}$ and $U_{2}$ do not commute at ${\bf k}_{1,2}$. This contradiction implies that these two momenta cannot contribute to a change in any topological invariant. Indeed we see in Fig. 4(b) that the diagonal winding number jumps when $U_{F}({\bf k})=\pm~{}\mathbb{I}_{4}$ at ${\bf k}={\bf k}_{\pm}$.

We can summarize the results presented here as follow. For $\Delta_{2}=0$, there are jumps in the Chern number whenever $U_{F}({\bf k})=\pm\mathbb{I}_{4}$ at any of the four momenta $(0,0)$, $(\pi,\pi)$, $(0,\pi)$ and $(\pi,0)$. For any value of $\Delta_{2}$, there are jumps in the diagonal winding number whenever $U_{F}({\bf k})=\pm\mathbb{I}_{4}$ at any of the two momenta $(0,0)$ and $(\pi,\pi)$.

Appendix B Contributions to the Chern number from other momenta

We will now study if Eq. (16) can be satisfied at momenta which differ from the four time-reversal invariant points and have no special symmetries. We will only consider the case $\Delta_{2}=0$ here. We would like to know whether the Floquet operator given in Eq. (10) can be equal to to $\pm\mathbb{I}_{4}$, when $H_{1}({\bf k})$ and $H_{2}({\bf k})$ are given by Eq. (1) or (14) and $M=M_{1}$ and $M_{2}$ respectively. It is clear that if ${\bf k}$ is not at one of the four time-reversal invariant points, $e^{-iH_{1}({\bf k})T/2}$ and $e^{-iH_{2}({\bf k})T/2}$ will generally not commute with each other. Hence Eq. (10) can be equal to $\pm\mathbb{I}_{4}$ only if $e^{-iH_{1}({\bf k})T/2}$ and $e^{-iH_{2}({\bf k})T/2}$ are separately equal to $\pm\mathbb{I}_{4}$. Equation (15) implies that this will happen if

[TABLE]

where $n_{1}$ and $n_{2}$ are positive integers. We now see that since there are two conditions to be satisfied, it may be possible to vary the two quantities $k_{x}$ and $k_{y}$ to satisfy Eqs. (24), provided that $T$ is large enough. We also see that such solutions for $(k_{x},k_{y})$ will appear in groups of four since these equations remain unchanged if $k_{x}\to-k_{x}$ or $k_{y}\to-k_{y}$. Hence the Chern number may be expected to change by $\pm 4$. This is precisely what we observe in Fig. 9 (a) at $T\simeq 4.49$.

We can derive the precise value of $T$ where the Chern number jumps by 4 in Fig. 9 (a) as follows. Given that $M_{1}=-0.9$, $M_{2}=-0.45$, $\Delta_{1}=1$, $t_{x}=1$ and $t_{y}=2$, we find numerically that there are peaks in the Berry curvature near $k_{x}=\pm 1.44$ and $k_{y}=\pm 1.30$ when $T$ is close to $4.49$. If we substitute these values of $k_{x}$, $k_{y}$ and $T$ in Eq. (24), we find that $n_{1}$ and $n_{2}$ are close to 1. Given this information, we can analytically derive the values of $k_{x}$, $k_{y}$ and $T$ which satisfy those equations with $n_{1}=n_{2}=1$ exactly. We first observe that $(M_{1}+t_{x}\cos k_{x}+t_{y}\cos k_{y})^{2}$ must be equal to $(M_{2}+t_{x}\cos k_{x}+t_{y}\cos k_{y})^{2}$. Since $M_{1}\neq M_{2}$, this implies that

[TABLE]

Next, the minimum value of $T$ where Eqs. (24) will be satisfied will correspond to the maximum value of $(M_{1}+t_{x}\cos k_{x}+t_{y}\cos k_{y})^{2}+\Delta_{1}^{2}(\sin^{2}k_{x}+\sin^{2}k_{y})$ as a function of $(k_{x},k_{y})$, subject to the condition in Eq. (25). Solving this maximization problem, we obtain a second condition

[TABLE]

Using Eqs. (25) and (26), we find that

[TABLE]

which gives four possible values of the momentum in the first Brillouin zone,

[TABLE]

Note that the values of $(k_{x},k_{y})$ depend on the parameters $t_{x},~{}t_{y},~{}M_{1}$ and $M_{2}$, unlike the time-reversal invariant momenta. Substituting Eq. (28) in Eqs. (24) gives

[TABLE]

which agrees well with the value where the Chern number jumps by 4 in Fig. 9 (a).

The general solution to Eqs. (24) can be found as follows. Since $M_{1}\neq M_{2}$, these equations imply that

[TABLE]

Substituting this in one of the equations in Eqs. (24), and using the identity $\sin^{2}k_{x}+\sin^{2}k_{y}=2-\cos^{2}k_{x}-\cos^{2}k_{y}$ and Eq. (30), we obtain a quadratic equation for $\cos k_{x}$. This generally gives two solutions for $\cos k_{x}$ and therefore for $\cos k_{y}$. We therefore get four possible values of $k_{x}$ and two values of $k_{y}$ in each case, giving a total of eight different solutions for $(k_{x},k_{y})$. (In special cases, these eight solutions become degenerate and reduce to four solutions).

We find that as $T$ increases beyond the value given in Eq. (29), there are more and more values of $(k_{x},k_{y})$ for which Eqs. (24) are satisfied for various integer values of $n_{1}$ and $n_{2}$. This gives rise to an increasing number of jumps in the Chern number.

In addition to the above momenta where the Floquet eigenvalues are exactly equal to $\pm 1$, we find that there are sometimes large numbers of momenta where the Floquet eigenvalues come very close to (but not exactly equal to) $\pm 1$, and the quasienergy gaps become very small. As a result, the Berry curvature becomes very large near those points; this gives rise to sharp fluctuations in the Chern number, since this is obtained by integrating the Berry curvature over the Brillouin zone. We see examples of this in Fig. 9 (a) at $T=4.81$ and $4.91$. These fluctuations in the Chern number are numerical artifacts. Namely, we have numerically studied the vicinity of the small quasienergy gaps using a very fine resolution in the momentum $(k_{x},k_{y})$, and we have found that although the gaps are very small, they are not exactly zero. The very small gaps lead to large fluctuations in the numerically calculated value of the Berry curvature and therefore of the Chern number.