Dynamical preparation of Floquet Chern insulators

Realizing topological insulators is of great current interest because of their remarkable properties and possible future applications. There are recent proposals, based on Floquet analyses, that one can generate topologically nontrivial insulators by periodically driving topologically trivial ones. Here we address what happens if one follows the dynamics in such systems. Specifically, we present an exact study of the time evolution of a graphene-like system subjected to a circularly polarized electric field. We prove that, for infinite (translationally invariant) systems, the Chern number is conserved under unitary evolution. For systems with boundaries, on the other hand, we show that a properly defined topological invariant, the Bott index, can change. Hence, it should be possible to experimentally prepare topological states starting from non-topological ones. We show that the chirality of the edge current in such systems can be controlled by adjusting the filling.

where H S is time-independent, H 1 (t) is time periodic with period T , and the amplitude f (t) is given by We restrict our analysis to non-interacting fermionic Hamiltonians, for which a complete characterization of the topological phases exists [16]. We assume that the initial state |ψ S 0 is the ground state of the static Hamiltonian H S . At time t > 0, the time-dependent wavefunction is |ψ(t) . We are interested in situations in which the undriven system is described by a topologically trivial static Hamiltonian H S and the driving is such that the Floquet Hamiltonian H F is topologically non-trivial. We stress that the fact that the Floquet Hamiltonian is topological does not imply that the wavefunction |ψ(t) is also topological. Hence, it is natural to ask whether it possible, under unitary evolution, to transform the ground state of a topologically trivial static Hamiltonian |ψ S 0 into the ground state of a topologically non-trivial Floquet Hamiltonian |ψ F 0 . If not, how "close" can |ψ(t) be to |ψ F 0 ? We focus on spinless fermions in a honeycomb lattice with nearest-neighbor hopping J and a staggered sublattice potential ∆ subject to a circularly polarized electric field E ac (t) = E 0 f (t) [− cos (Ω t) , sin (Ω t)] [6,7]. In the electromagnetic gauge, in which the vector potential is zero, the time-dependent Hamiltonian is given by Eq. (1) with where the sum in H S extends over nearest neighbor sites, α ∈ {1, 2} indicates the sublattices A and B, respectively, n i,α = c † i,α c i,α are the site number operators, U E ( r, t) = −e r · E ac (t) is the electrostatic potential energy, and e is the electric charge.
We work in the high-frequency limit in which the driving frequency is larger than the single-particle bandwidth [17], i.e., Ω > 6 J. Moreover, in order to obtain a non-trivial high-frequency limit, we scale the electric field with the frequency of the driving [18], eaE 0 ∼ Ω where a is the lattice spacing. Our parameters are: and are chosen so that the (effective) Floquet Hamiltonian H F is topological. The staggered sublattice potential ∆ is introduced to make direct connection with the experiment in Ref. [17], and to ensure that the edge modes that are not topological in nature are gapped out. The period of the driving is T = 2π/Ω and we consider ramping times τ /T ∈ [0−10 3 ]. We stress that this choice of parameters is relevant for the recent experimental realization of the Haldane model in cold atoms [17]. In Ref. [17], τ = 20 ms and 1/T = Ω/(2π) = 4 kHz so that τ /T = 80. However, the loading procedure there was more complex than the linear ramp considered here. Infinite System. We first consider the infinite (translationally invariant) system which, at half-filling, can be mapped onto a collection of independent pseudo spin-1 2 by going to momentum space. The Hamiltonian H = k H k and the density matrix ρ = |ψ ψ| = k ρ k are (we take = 1 in what follows): Here, 1 2×2 is the 2 × 2 identity matrix, σ are the Pauli matrices, and S k and B k are three dimensional, timedependent, vectors fields defined in the two-dimensional Brillouin Zone (BZ) (see Fig. 1a). For a pure state, the vector S k has unit length and the Chern number (Ch) of the state is simply the number of wrappings of the pseudo spin configuration around the Bloch sphere [5]: Here the integral extends over the BZ. In the ground state, the pseudo spin configuration is parallel to the pseudo magnetic field. This does not need to be the case out of equilibrium. The exact equation of motion is: which is simply the precession of the pseudo spin S k around the pseudo magnetic field B k . With this mapping, the ground states |ψ S 0 and |ψ F 0 obtained by filling the valence bands of H S and H F are represented by the pseudo spin configurations S S k and S F k , respectively. The overlap of the two ground states is: For the parameters chosen [see Eq. (see Figs. 1b and 1c). This implies that there is at least one k-point in the BZ for which the vectors S S k and S F k are anti-parallel and, as a consequence, the overlap of the ground states is identically zero, ψ S 0 |ψ F 0 2 = 0. Despite this fact, |ψ S 0 has a large occupation on the valence band of H F . This is because the pseudo spin configurations S S k and S F k differ significantly only close to the Dirac points, while they are almost identical in the rest of the BZ (see Fig. 1d). The normalized occupation (Occ) of the Floquet valence band is 3 is the area of the BZ. We can now consider the dynamical process by which the periodic driving is turned on. In principle, the Chern number inherits a time-dependent from the timedependence of the pseudo spin configuration S(t) obtained by integrating the equation of motion (6) subject to the initial condition S k = S S k . However, a straightforward calculation shows that this is not the case. This follows from the fact that ∂ t Ch can be written as: If S k (t) and B k (t) are sufficiently smooth vector fields in the BZ then it follows that the expression above is identically zero (see supplementary material). From Eq. (6) one can see that, under a smooth pseudo magnetic field B k (t) and an initially smooth pseudo spin configuration, S k (t) remains smooth. We can therefore formulate a nogo theorem as follows: If the initial pseudo spin configuration is smooth (at least C 1 ) in the Brillouin zone and the pseudo magnetic field is smooth (at least C 2 ), then the Chern number is conserved under the unitary evolution generated by the pseudo magnetic field.
We should stress that the smoothness condition on B k (t) is not very restrictive. For example, the band structure of graphene is singular at the two Dirac points, but the pseudo magnetic field configuration: is analytic in the BZ and satisfies the condition of the theorem. The statement that the Chern number cannot change independently of the time-evolution considered is similar to the result that, under unitary evolution, a pure state remains pure. Both results do not predict the exact wavefunction at the end of a dynamical process but constrain the possible outcomes. A straightforward application of our result is that the time-evolved wavefunction has Chern number equal to zero and therefore it must be orthogonal to the Floquet ground state (which has Chern number one): System with Boundaries. Experimental systems have boundaries, so here we address what happens when translational invariance is broken. We consider a finite system with open boundary conditions (see inset in Fig. 2a). The system contains 928 lattice sites evenly divided into A and B sublattices. To characterize the topological properties of systems with broken translational symmetry one cannot rely on the Chern number. We use two complementary indicators: (i) the cumulative local density of states CLDOS(E) = E −∞ dε LDOS(ε) and (ii) the Bott index [9]. The Bott index is a topological invariant that is the equivalent of the Chern number for system with boundaries. In equilibrium, the Bott index at some energy E is computed by projecting some special matrices onto the eigenstates of the Hamiltonian with energy ε < E. This assumes that the eigenstates with energy ε < E are fully occupied while the eigenstates with energies ε > E are empty. We have extended the Bott index definition to non-equilibrium situations by taking into account the non-equilibrium occupation of the eigenstates (see Methods). We stress that, while the CLDOS depends exclusively on the Hamiltonian, the Bott index is a property of the state as it depends, in a crucial way, on the occupation of the energy levels. In particular, in out of equilibrium situations, the Bott index inherits a time-dependence from the time dependence of the state. In Fig. 2b and 2c we show these two indicators for the static Hamiltonian H S and the Floquet Hamiltonian H F , respectively. The CLDOS(E) of H S (both at the center of the sample and along the edges) has a plateau around E = 0 signifying that there are no states at E = 0, i.e., the system is gapped. Moreover, the Bott index is identically zero indicating that both H S and its ground state |ψ S 0 are topologically trivial. On the contrary, H F has edge states inside the bulk gap, as shown by the fi- nite (zero) slope of the CLDOS at the edge (center) for E ≈ 0. The existence of topologically protected edge modes is confirmed by the Bott index. In fact, in equilibrium, the Bott index at some energy, Bott(E), is equal to the number of edge states at that energy. As one can see in Fig. 2b, Bott(E 0) = 1, indicating that the ground state of H F at half filling (E = 0) is topologically nontrivial.
In order to study the adiabatic turning on of the periodic driving, we solve the time-dependent Schrödinger equation [19,20] subject to the initial condition |ψ(t = 0) = |ψ S 0 (see Methods). At stroboscopic times t n = n T during the time evolution, we monitor the Bott index and the overlaps of |ψ(t) with both the initial state and the Floquet ground state: Their behavior, for ramping time τ = 80 T , are shown in Fig. 2a. One can see that overlap with the initial state decays to zero rapidly while the overlap with the ground state of the Floquet Hamiltonian increases. For t n = n T > τ , the electric field has reached its final value and the overlap with |ψ F 0 becomes independent on n since |ψ F 0 is an eigenstate of the evolution operator over a period, i.e. ψ F 0 |ψ(t n+1 ) Interestingly, for the parameters chosen, at t ≈ τ −11T (i.e., slightly before the electric field is fully on) the Bott index jumps from zero and becomes one, indicating the wavefunction has acquired a topological character. This is in stark contrast with the behavior of the infinite system in which the Chern number is strictly conserved during the time-evolution. We also note that the overlap with the Floquet ground state |ψ F 0 increases non-monotonically with time. This suggests that the final overlap ψ F 0 |ψ (t = ∞) 2 can be increased by using more sophisticated ramping protocols. For example, by instantaneously quenching the amplitude of the electric field from its value when the overlap ψ F 0 |ψ (t) 2 has a local maximum to its final value. In the inset of Fig. 2a, we plot the value of the overlap with the Floquet ground state at the end of the ramp, i.e., ψ F 0 |ψ (t = ∞) 2 , for different ramping times τ and observe that, as expected, it generally increases with increasing τ .
To make further contact with experiments, we investigate the current that flows through the sample under driving which is, in principle, a measurale quantity [21]. We note that the physical current is different from the current one would obtain using the Floquet Hamiltonian [18,22] (see also Supplementary Information). The former connects only nearest neighbor sites, while the latter can flow between far away lattice sites if there are longer range hopping terms in H F . Contrary to the overlap, the current continues to evolve for t > τ , when the electric field has already reached its final value. After the end of the ramp, the current is not conserved locally since the site occupancies are not stationary. We therefore wait for the site occupations to reach their asymptotic values (up to small oscillations within a driving period) before measuring the (instantaneous) current, which is then averaged over 20 periods to obtain the data shown in Fig 3b. We note that the current is, as expected, concentrated along the edges but contains both clockwise and counterclockwise components.
Discussion. The exact non-equilibrium dynamics of many-body noninteracting fermions during the switching on of a periodic perturbation appears to be nontrivial in systems with nontrivial topology. The two topological invariants studied here, the Chern number in translationally invariant (thermodynamically large) systems and the Bott index in systems with boundaries, exhibit completely different behavior. The Chern number is conserved while the Bott index is found to change under unitary evolution. The exact conservation of the Chern number might appear surprising since the naive expectation is that, during any non-equilibrium process, some excitations are generated and the final state, which corresponds to a partially filled valence and conductance band need not have a quantized Chern number. However, this argument does not take into account that, under unitary evolution, each quasi-momentum k is in a coherent superposition of the valence and conduction band. It is precisely this coherence that prevents the Chern number from changing. Hence, within the conditions discussed, it is not possible to adiabatically transform a topologically trivial into a topologically nontrivial wavefunction. This is something that might also occur in interacting systems as topological phase transitions are first order (discontinuous) transitions [23,24]. Our results open many new interesting questions: (i) What is the nature of the topological transition in systems with boundaries? (ii) What is the dynamics of the edge states [25] across those transitions? (iii) What happens in the presence of interactions [26,27] and/or a coupling to a bath?

METHODS
Bott index for out-of-equilibrium systems. Consider a single-particle Hamiltonian (defined by a matrix H) on a lattice. Given the two diagonal matrices X i,j = x i δ i,j and Y i,j = y i δ i,j , where x i and y i are the coordinates of the i th lattice site, we defined two unitary matrices: where L x,y are the linear dimensions of the system. Let P be the projectors onto the eigenstates up to energy E, i.e., P ≡ <E | |. In equilibrium, the Bott index at energy E is defined as [9]: , whereŨ x = P U x P andŨ y = P U y P are matrices U x and U y projected onto the states up to energy E. We generalize the definition of the Bott index by taking into account the non-equilibrium occupation of the states. The Bott index of the occupied states is obtained by replacing the matricesŨ x andŨ y with the matrices U x (t) and U y (t): where P (t) is the projector on states occupied at time t. For example, if we had P = |α α| + |β β| at t = 0, at time t the projector becomes (for a system of non-interacting particles as the one we are considering) P (t) = |α(t) α(t)| + |β(t) β(t)|, where |α(t) = U (t)|α and |β(t) = U (t)|β are the time-evolved states. Numerical simulations for system with boundaries. The time-dependent Hamiltonian is given by Eq. (1). Because of its noninteracting character, this problem can be efficiently solved in the single-particle basis [19,20]. We denote as H S and H 1 (t) the N s × N s matrices (N s being the number of lattice sites) that represent the static and time-dependent parts of the Hamiltonian in real space.
The evolution operator over a cycle is given by: where t j = j T N and U (t j+1 , t j ) is computed using a second-order Trotter-Suzuki decomposition [28][29][30]: Since H S is time-independent, e −iδt H S needs to be computed only once (this is done by diagonalizing H S ). This leaves the computation of e − iδt 2 H1(t+δt/2) , from the already diagonal H 1 (t), to be computed at each time step. By exact diagonalization of U (T, 0), we obtain the Floquet eigenstates and eigenvalues: from which the single-particle Floquet Hamiltonian can be explicitly built as H F = l |ψ l ε l ψ l | where ε l = T θ l . The lowest energy single-particle eigenstates of H F (H S ) are then collected into a rectangular matrix W F (W S ) of size N s × N p , where N p is the number of particles in the ground state (at half filling N p = N s /2).
The time evolution of the many-particle system is obtained by multiplying the matrix W S from the left with the square matrix U (T, 0) of size N s × N s . The overlaps between many-particle wavefunctions can also be easily computed as determinants of products of matrices such as W F and W S , and their adjoints [19,20]. Moreover, the c † i c j elements of the equal-time single-particle density matrix are given by the i, j element of the square matrix W 0 W † 0 of size N s × N s (here the overline indicates complex conjugation). From the knowledge of the equal-time single-particle density matrix many observables (such as the current) can be easily computed (see Supplementary  Information).

I. ACKNOWLEDGMENT
We thank M. Bukov, T. Iadecola, M. Kolodrubetz, T. C. Lang, A. Polkovnikov, and K. Sun for illuminating discussions. We thank Kolodrubetz for bringing to our attention Ref. [25] and T. C. Lang for sharing his software used to plot currents. This work was supported by the Office of Naval Research.
where we have used the distribution property of the cross product ∂ x (S × B) = S x × B + S × B x . Now we can apply the Binet-Cauchy identity to obtain: In a similar way we get: Putting all together, and carrying out the cancellations, we get: So far this expression is completely general. Now we use that the initial state is a pure state: from which it follows that S · S = 1, i.e., S is a unit vector for any point in the BZ. We then observe that: Therefore we arrive at Eq. (8) in the main text. We observe that: If the vector field B is smooth the mixed derivative commute, i.e., B x,y = B y,x , and we arrive at: If ∂ x (B y · S) and ∂ y (B x · S) are continuous, we can use the periodicity of B and S in the Brillouin zone to obtain: We note that, up to this point, we have simply shown thatĊh(t) = 0 if S(t) and B(t) are sufficiently smooth and S(t) represents a pure state, i.e., S · S = 1. However, to prove that the Chern number is conserved at all times, we still need to prove that under time evolution (i) a pure state remains pure and (ii) a smooth pseudo-spin configuration remains smooth. To verify that this is indeed the case, we look into the equation of motion (6). We note that, under this equation, the length of the vector S is conserved, i.e., ∂ t (S · S) = 0, and therefore the condition (i) is verified. We also note that the equation of motion is a linear differential equation. If S(t = 0) is smooth in k x , k y and B(t) is smooth in k x , k y then S(t) remains smooth at all times. Therefore, the condition (ii) is verified. This concludes the proof of the theorem.
Edge States in the system with boundaries In Fig. 4 we show some single-particle Floquet eigenstates at energies ε ≤ 0 within and outside the bulk gap (see Fig. 2b). As expected, the wavefunctions with energies within the bulk gap are localized on the edges. It is interesting to note that these states have much larger occupation on the zigzag edges (left and right) than on the armchairs (top and bottom) edges. The chiral nature of the edge states can be probed by preparing a single-particle δ-like wavefunction in the center of the left and right edge. The subsequent time-evolution shows that some particle density diffuse towards the center of the sample while the rest remains localized on the edges and moves clockwise. The fraction of charge that persists on the edge is determined by the overlap between the initial δ-like wavefunction and the edge state and it is expected to scales as 1/L where L is the length of the edge.