The anomalous Floquet-Anderson insulator as a non-adiabatic quantized charge pump

Periodically driven quantum systems provide a novel and versatile platform for realizing topological phenomena. Among these are analogs of topological insulators and superconductors, attainable in static systems; however, some of these phenomena are unique to the periodically driven case. Here, we show that disordered, periodically driven systems admit an"anomalous"two dimensional phase, whose quasi-energy spectrum consists of chiral edge modes that coexist with a fully localized bulk - an impossibility for static Hamiltonians. This unique situation serves as the basis for a new topologically-protected non-equilibrium transport phenomenon: quantized non-adiabatic charge pumping. We identify the bulk topological invariant that characterizes the new phase (which we call the"anomalous Floquet Anderson Insulator", or AFAI). We provide explicit models which constitute a proof of principle for the existence of the new phase. Finally, we present evidence that the disorder-driven transition from the AFAI to a trivial, fully localized phase is in the same universality class as the quantum Hall plateau transition.


I. INTRODUCTION
Time-dependent driving opens many new routes for realizing and studying topological phenomena in manybody quantum systems. For periodic driving, Floquet theory provides a powerful framework for analyzing quantum dynamics. In particular, periodically driven crystalline systems can be characterized in terms of Floquet-Bloch band structures, analogous to the band structures of nondriven systems. Recently, an intense wave of activity has developed around exploring the possibilities of using periodic driving to realize "Floquet topological insulators," i.e., driven system analogues of topological insulators  in a variety of solid-state [23], atomic, and optical contexts [24,25].
In addition to the phenomena above, driven systems may also host unique types of robust topological phenomena, without analogues in static systems [3,[26][27][28][29][30][31][32]. A prominent example is Thouless's quantized adiabatic pump [33]: a gapped one-dimensional system transmits a precisely quantized charge when subjected to a periodic modulation, in the adiabatic limit of slow driving. More recently, a new example was discovered: a two-dimensional driven system can support chiral edge states even if all of its bulk Floquet bands have zero Chern numbers [3,27]. This situation stands in sharp contrast to that of static two-dimensional systems, where the existence of chiral edge states is intimately tied to the topological structure of the system's bulk bands, as captured by their Chern numbers [34]. A system exhibiting this anomalous behavior was recently realized using microwave photonic networks [35,36].
In this work, we show that the unique topological characteristics of periodically driven systems can lead to qualitatively new phenomena when spatial disorder is introduced. First, it is possible for robust chiral edge states to exist in a two-dimensional driven system where all bulk states are Anderson localized; we refer to such a system as an anomalous Floquet-Anderson insulator (AFAI). This situation cannot occur in the absence of driving, where the existence of chiral edge states necessarily implies that there must be delocalized bulk states at some energies [37]. Crucially, in an AFAI this relation is circumvented FIG. 1. The anomalous Floquet-Anderson insulator (AFAI), in a disordered two-dimensional periodically driven system with time-dependent Hamiltonian HðtÞ. In the AFAI phase all bulk states are localized, yet the system hosts chiral propagating edge states at all quasienergies. The nontrivial topology of the phase is characterized by a nonzero value of the winding number defined in Eq. (8).
Published by the American Physical Society under the terms of the Creative Commons Attribution 3.0 License. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI. by the periodicity of quasienergy: the edge states persist through all quasienergies, completely wrapping around the quasienergy Brillouin zone. Moreover, the combination of these novel chiral edge states and a fully localized bulk gives rise to an intriguing nonequilbrium topological transport phenomenon: quantized nonadiabatic charge pumping.
This paper is organized as follows. In Sec. II, we present a heuristic physical picture of the AFAI phase and give an overview its properties. Section III analyzes a simple solvable model that serves as a proof of principle for the AFAI phase. In Sec. IV, we discuss the topological invariant characterizing the AFAI and show that the AFAI exhibits edge modes at every quasienergy. In Sec. V, we show how the edge mode structure leads to quantized charge pumping. In Sec. VI, we conduct a numerical study of a wider class of models exhibiting the AFAI phase. We numerically demonstrate the properties discussed in Secs. IV and V. At strong disorder, we find a topological transition between the AFAI and a "trivial" Floquet insulator where all states are localized (including at the edges).

II. PHYSICAL PICTURE AND SUMMARY OF THE MAIN RESULTS
To begin, in this section we summarize the main results of this work. In particular, we describe the unique spectral characteristics of the AFAI and the novel nonadiabatic quantized pumping phenomenon that it hosts. Our aim in this section is to provide a heuristic-level picture of our findings; detailed derivations and further discussion follow in Secs. IV-VI.
The AFAI is a unique phase occurring in disordered, periodically driven two-dimensional systems. The system evolves according to a time-periodic Hamiltonian, HðtÞ ¼ Hðt þ TÞ, where T is the driving period. For describing the long-time behavior of such systems it is useful to study the stroboscopic "Floquet" evolution operator UðTÞ ¼ T e −i R T 0 dtHðtÞ . The eigenvalues of UðTÞ are phases, e −iεT , where the quasienergy ε is defined in a periodic Floquet-Brillouin zone with period Ω ¼ 2π=T.
In two spatial dimensions, disordered, periodically driven systems may exhibit a variety of phases. Some of these phases have direct analogies in nondriven systems. In cases where such analogies exist, all features of the driven system can be derived from an associated timeindependent "effective Hamiltonian" H eff , defined such that UðTÞ ¼ e −iTH eff . Topological characteristics, such as the presence or absence of chiral edge states at sample boundaries, are, in particular, captured in those cases by the effective Hamiltonian and its associated Chern numbers, just as for nondriven systems (recall that Chern numbers provide a full topological characterization of noninteracting static systems without symmetries). As a result, many phenomena exhibited by static systems can be mimicked by periodically driven systems; examples include the direct correspondence between chiral edge states and bulk Chern numbers described above, as well as disorder-induced topological transitions, as exhibited by the "topological Anderson insulator" [38] and its Floquet counterpart [11].
The AFAI phase we introduce in this work is a phase of a disordered periodically driven systems whose characteristics are qualitatively distinct from those achievable in the absence of driving. Its defining property is that all its bulk Floquet states are localized by the disorder; nevertheless, its edges support chiral edge states. This unusual situation has a number of intriguing physical consequences, as we describe below.
In the absence of driving, chiral edge states must be accompanied by delocalized states in the bulk of the system. This can be seen by considering a system in an annular geometry and tracking how its spectrum evolves as magnetic flux is threaded through the hole of the annulus. Once a full flux quantum is inserted, the Hamiltonian is equivalent to the original one and therefore its spectrum must be unchanged compared to the original one. Tracing the evolution of a given state as the flux is inserted, there are two options once a full flux quantum is reached: (1) the state returns to its original energy or (2) the state "flows" to a new energy [39]. The edge states evolve according to option (2). The only way for this spectral flow to terminate is if a delocalized bulk state is reached, connecting the upward-and downward-flowing families of states on opposite edges. Thus, we see that chiral edge states cannot exist without delocalized bulk states, as otherwise the spectral flow would have to continue up and down to infinite energies.
The above argument fails when considering a periodically driven system, where the quasienergy spectrum is periodic with a period Ω ¼ 2π=T. In this case, the flow of the edge states need not terminate in a delocalized bulk state. Instead, the flow of the edge states can "wrap" around the quasienergy zone. In this light, it appears that it may be possible to find a system that exhibits chiral edge states and at the same time has all of the bulk states localized. If this unique situation can indeed be realized, it would furthermore imply that the chiral edge states must be present at every quasienergy. Can these intriguing properties be realized in a twodimensional, disordered, periodically driven system? The current work is the first to address this question.
To prove the existence of this anomalous Floquet-Anderson insulator phase, in Sec. III we start with a simple exactly solvable periodically driven system, whose Floquet spectrum and eigenstates can be obtained explicitly. We then study its behavior in the presence of weak perturbations that break the solvability, and show that the AFAI phase is robust over a finite extent in parameter space. The nontrivial character of this phase is captured by an integervalued bulk topological invariant W that characterizes the Floquet operator of a system on a torus with threaded fluxes (Sec. IVA). The value of the bulk invariant gives the number of chiral edge states that would appear in a geometry with boundaries. In Sec. IV B, we show that the edge state spectrum itself can be directly characterized in terms of a different winding number invariant, defined for a system in a cylindrical geometry.
Our analysis is complemented by numerical simulations, which explore both the weak and strong disorder regimes. As the disorder strength is increased, we find that the system undergoes a phase transition to a trivial Anderson insulator; see Sec. VI B.
The above properties describe the single-particle characteristics of the AFAI and lead to its defining characteristic as a many-body system: robust quantized charge pumping that persists in a nonadiabatic driving regime. This behavior is in contrast to that of the one-dimensional "Thouless pump," where the pumped charge is quantized only in the limit of infinitely slow driving [33]. The setup that realizes nonadiabatic quantized pumping is illustrated in Fig. 1. We consider a strip of AFAI in which all of the states close to one edge are populated with fermions. In this situation, the total current flowing through the strip is quantized as an integer times the inverse driving period, hIi ¼ eW=T. Here, W is the bulk topological invariant discussed above, and hIi is the current averaged over many driving periods. The quantized charge pumping is a direct result of the edge structure defined above: when the fermion occupation is of the form shown in Fig. 1, the edge states on one side are completely filled, while the localization of the bulk states prevents current from flowing in the direction perpendicular to the edge. In Sec. V, we derive the relation between the quantization of the charge current and the bulk topological invariant and discuss temporal fluctuations about the quantized value.
Putting together all of the characteristics described above, we may thus define the AFAI as a disordered, periodically driven system in which all the bulk Floquet eigenstates are localized, and (i) the quasienergy independent winding number W is nonzero, (ii) chiral edge modes propagate along sample boundaries at all quasienergies, and (iii) a quantized current is pumped whenever all states along one edge are filled with fermions. As we show below, the properties (i)-(iii) all follow from each other.

III. EXPLICIT DEMONSTRATION OF THE ANOMALOUS FLOQUET-ANDERSON PHASE
In this section, we study a simple model that allows us to explicitly demonstrate the existence and robustness of the AFAI phase. We start from an exactly solvable model for a periodically driven disordered system, which exhibits localized bulk bands with zero Chern numbers and chiral edge states at all quasienergies. We then argue that these properties are robust to generic small perturbations (which break the solvability): upon adding perturbations, the bulk states remain localized and the chiral edge states persist. The model thus proves the existence of a phase of periodically driven systems in which chiral (Floquet) edge states may coexist with a fully localized bulk. This is the AFAI phase. In the next section, we discuss the robustness of the phase and the edge and bulk topological invariants that characterize its universality.
As a starting point, we consider a clean (nondisordered) system on a square lattice, introduced in Ref. [27]. The timeperiodic, piecewise-constant Hamiltonian is of the form H clean ðtÞ ¼ H n , for ½ðn − 1ÞT=5 ≤ t < ðnT=5Þ, n ¼ 1; …; 5. The square lattice is divided into two sublattices, A and B (shown as filled and empty circles in Fig. 2). During each of the first four segments of the driving, n ¼ 1; …; 4, hopping matrix elements of strength J between the A and B sublattices are turned on and off in a cyclic, clockwise fashion, as shown in Fig. 2(a): during segment n ¼ 1, 2, 3, or 4, each site in the A sublattice is connected by hopping to the site above, to the right, below, or to the left of it, respectively. In the fifth segment of the period, all of the hoppings are set to zero, and an on-site potential δ A;B is applied on the A and B sublattice sites, respectively.
We choose the hopping strength J such that ½ðJTÞ=5 ¼ ðπ=2Þ. For this value of J, during each hopping During steps 1-4, nearest-neighbor hopping is applied along the colored bonds as shown. The hopping strength J is chosen such that a particle hops between adjacent sites with probability one during each step. In the fifth step, all hopping is turned off and a random disorder potential is applied (the same potential is used for all subsequent driving cycles). (b) Over one driving period, a particle initialized on any site in the bulk returns precisely to its original position (blue arrow). Along the upper edge, a particle initialized on site ðx; 0Þ of the A sublattice shifts two sites to the right (red straight arrow) and acquires a phase e −i½π=2þV ðxþ2;0Þ T=5 . (c) With disorder, the quasienergy spectrum consists of two nonoverlapping bands of width 2δV, centered at ε ¼ AEπ=2T.
segment of the driving period a particle that starts on one of the sites hops to the neighboring site with unit probability. The on-site potential, applied only while all hopping matrix elements are turned off, is chosen to be δ A;B ¼ AEðπ=2TÞ.
With the parameter values chosen above, it is easy to find the Floquet eigenstates and quasienergies of the timedependent Hamiltonian H clean ðtÞ. The bulk spectrum consists of two flat Floquet bands with quasienergies AEðπ=2TÞ, with the corresponding Floquet eigenstates localized on either the A or B sublattice. Over each driving period, a particle initially localized on either an A or a B site in the bulk encircles a single plaquette and returns to its original position; see blue arrow in Fig. 2(b). In a cylindrical geometry, motion along the edge is also easily visualized: particles on the first row of sites in the A (B) sublattice along the upper (lower) edge shift one unit cell to the right (left), as shown by the red arrow on the upper edge in Fig. 2(b). The corresponding eigenstates are, therefore, plane waves, localized on the first row of sites in the A (B) sublattice along the upper (lower) edge. The two edges therefore host linearly dispersing chiral modes in the quasienergy gaps between the two bulk bands.
We now introduce a specific form of a time-dependent disorder potential VðtÞ, which still allows for an exact solution. The full time-dependent Hamiltonian is given by H 0 ðtÞ ¼ H clean ðtÞ þ VðtÞ. During the fifth segment of the driving period, we let VðtÞ ¼ P r V r c † r c r , where V r is drawn from a smooth distribution with support in the range ½−δV; δV, and c r is the annihilation operator at the integervalued lattice position r ¼ ðx; yÞ. During segments 1-4, VðtÞ ¼ 0. We choose δV < ðπ=2TÞ, ensuring that the gap between bands remains open.
By following the evolution of a state that is localized on a single bulk site r at time t ¼ 0, one can easily verify that this state is a Floquet eigenstate, whose quasienergy is AEðπ=2TÞ þ V r [here, þ1 (−1) refers to an initial site in the A (B) sublattice]. With periodic boundary conditions, the Floquet spectrum consists of two bands, with quasienergies in the ranges ½AEðπ=2TÞ − δV; AEðπ=2TÞ þ δV; see In a cylindrical geometry, one can similarly follow the evolution of a state that is initially localized on a site at the edge, as discussed above and in Fig. 2(b). A particle initialized at site r ¼ ðx; 0Þ on the A sublattice on the upper edge is translated by two sites to the right, and then picks up a phase e −i½π=2þV ðxþ2;0Þ T=5 . The Floquet eigenstates localized at the edge are no longer perfect plane waves, but they do remain extended with support entirely on the first row of sites along the edge.
To explicitly construct the edge eigenstates, note that the shift operation [denoted S in Fig. 2(b)] followed by the application of a local phase gives where ψ x ðtÞ is the wave function of an edge eigenstate at ðx; 0Þ. The eigenvalue relation further yields ψ xþ2 ðTÞ ¼ ψ xþ2 ð0Þe −iεT . Taking the ansatz ψ x ¼ e iðkxþδϕ x Þ yields ε ¼ k=T and δϕ xþ2 ¼ δϕ x − ðπ=2Þ − V ðxþ2;0Þ T=5. This recurrence relation for the phases is trivially solved by iteration for a system in the geometry of an infinite strip, for any value of k. For a finite cylinder with 2N sites around the perimeter, the edge state closes on itself and yields an additional nontrivial condition: Satisfying this constraint imposes a quantization condition on the allowed values of k. In the thermodynamic limit, all values 0 ≤ k < 2π become allowed and we explicitly see that chiral edge states persist at all quasienergies, in the presence of a fully localized bulk. By the definition presented in Sec. II, the Hamiltonian H 0 ðtÞ ¼ H clean ðtÞ þ VðtÞ thus realizes the AFAI phase. Clearly, the above model utilizes a very specific form of the periodic driving and of the added disorder. Nevertheless, we argue that the AFAI is a robust phase that does not require fine-tuning. To demonstrate the robustness of the phase, we now consider a generic local perturbation of H 0 ðtÞ that preserves the periodicity in time, H λ ðtÞ ¼ H 0 ðtÞ þ λDðtÞ, and show that the AFAI phase survives up to a finite value of λ.
The perturbation DðtÞ is assumed to be periodic in time and short ranged in real space, such that the matrix elements of DðtÞ vanish beyond the rth neighbor on the square lattice. For δV ¼ 0 (no disorder) and λ ≠ 0, the bulk eigenstates of UðTÞ are generically dispersive and delocalized. However, we argue that for δV > 0 and for a sufficiently small λ, all of the bulk Floquet states remain localized. To see this, we derive a time-independent effective Hamiltonian H eff λ for the Floquet problem (on the torus), with δV ≠ 0, λ ≠ 0: where T denotes time ordering. We further write the effective Hamiltonian as H eff corresponds to the unperturbed (λ ¼ 0) effective Hamiltonian, defined such that its eigenstates lie in the range ½−ðπ=TÞ; ðπ=TÞÞ. Here, we are considering a system with periodic boundary conditions; we comment on the edge states later.
The key point, which we show below, is that for sufficiently small δV and λ, the hopping matrix elements of the effective static Hamiltonian H eff ð0Þ þ D eff decay exponentially with distance. If, in addition, λ ≪ δV=Ω, then all of the Floquet eigenstates remain localized [40].
To find the effective Hamiltonian for λ ≠ 0 we need to solve for D eff . The unperturbed effective Hamiltonian is of the form where η r ¼ 0 (1) for j on the A (B) sublattice. We express D eff as a power series in λ: To find D ðmÞ eff , we expand both sides of Eq. (2) in powers of λ and compare them order by order. The details of the calculation are given in Appendix A. The results can be summarized by considering the explicit representation of the operators D ðnÞ eff as a tight-binding Hamiltonian: Using the explicit form for H eff ð0Þ given in Eq.
(3), we find for the lowest-order term with DðtÞ ¼ U 0 ðt; 0Þ † DðtÞU 0 ðt; 0Þ, where U 0 ðt; 0Þ is the λ ¼ 0 evolution operator, and E r; is the zeroth-order quasienergy difference between the states localized at sites r and r 0 . As long as E r;r 0 T is smaller than 2π for every pair of sites, which is the case for δV < ½π=ð2TÞ, the factor ½ðiE r;r 0 Þ=ðe iE r;r 0 T − 1Þ in Eq. (6) is bounded. Similarly, the matrix elements Δ ðnÞ r;r 0 are all nonsingular (see Appendix A). Under these conditions, we expect the expansion in powers of λ to converge. In Appendix A, we argue that for sufficiently small λ, the matrix elements of H eff λ decay exponentially with distance. Therefore, H eff λ has the form of a tight-binding model with random on-site potentials and weak, short-range hopping. In this context, we expect all states to remain localized up to a critical strength of λ.
Since all of the bulk states remain localized as λ is turned on, the chiral edge states that exist for λ ¼ 0 cannot disappear; the only way to remove them is by closing the mobility gap in the bulk, allowing the two counterpropagating states at the two opposite edges to backscatter into each other. Hence, we expect the chiral edge states to persist up to a critical value of λ where the bulk mobility gap closes.
In this section, we used a simple explicit model to demonstrate the existence of the AFAI, proving by example that chiral edge states may exist in a periodically driven system in which all bulk states are localized. Next, we show that this behavior, in fact, has a topological origin, and is thus much more general than the example used for this existence proof.

IV. TOPOLOGICAL INVARIANTS OF THE AFAI
In this section, we derive "bulk" and "edge" topological invariants that characterize the robust features of the AFAI spectrum for systems in periodic and open geometries, respectively. As we show, the bulk invariant, computed for an AFAI on a torus, gives the net number of chiral edge modes when the system is opened to a cylindrical geometry.

A. Bulk invariant
We now consider a generic disordered twodimensional system with a time-periodic Hamiltonian HðtÞ ¼ Hðt þ TÞ. No special form of disorder is assumed. As an additional ingredient, we also consider constant (time-independent) fluxes, Θ ¼ ðθ x ; θ y Þ, which are threaded through the torus [41]. This yields a family of time-dependent Hamiltonians HðΘ; tÞ, and their associated evolution operators, UðΘ; tÞ HðΘ;t 0 Þdt 0 . As a first step in constructing the bulk topological invariant, we define an associated, "deformed" timeperiodic evolution operator for the system on a torus: with H eff ε ðΘÞ ¼ ði=TÞ log UðΘ; TÞ. Note that, by construction, U ε ðΘ; TÞ ¼ 1. The explicit dependence on ε in the above definitions comes from the necessary choice of a branch cut for log; we use a definition such that With these definitions at hand, we can define the "winding number" [42]: In Eq. (8), we use the shorthand U ε ≡ U ε ðΘ; tÞ, and W ε is an integer, which can, in principle, depend on the quasienergy ε. Note that in order for W ε to be well defined, the quasienergy ε has to remain in a spectral gap of UðΘ; TÞ for every value of the threaded fluxes Θ (otherwise, the operator U ε is discontinuous as a function of Θ). For a large enough disordered system, almost all values of ε satisfy this requirement, since upon changing fluxes θ x and θ y , the quasienergies of the localized bulk states change only by an amount proportional to e −L=ξ , where ξ is the localization length and L is the linear system size. In contrast, the average level spacing is proportional to 1=L 2 . When all of the eigenstates of UðΘ; TÞ are localized, the invariant W ε is, in fact, independent of ε. This follows from the relation between the winding number W ε and the Chern numbers characterizing the eigenstates of UðΘ; TÞ [27], In the above, C ε 1 ;ε 2 is the total Chern number of the eigenstates with quasienergies between ε 1 and ε 2 : is a projector onto the eigenstates of UðΘ; TÞ with quasienergies between ε 1 and ε 2 . If all of the bulk eigenstates are localized, C ε 1 ;ε 2 ¼ 0 [43]. Therefore, in this case, by Eq. (9), W ε 1 ¼ W ε 2 for every pair of quasienergies. Then, we can drop the subscript ε, and refer to the winding number simply as W.

B. Edge invariant and chiral edge states
We now concentrate on the edge structure of the AFAI, considering a system in a cylindrical geometry. To reach the AFAI phase, we envision a generalization of the model in Sec. III in which disorder is added to a (translationally invariant) Floquet-Bloch system, for which all the Chern numbers of UðTÞ vanish but with one chiral edge state on each edge running through each of the bulk gaps of the quasienergy spectrum. The setup and spectrum are shown schematically in Figs. 3(a) and 3(b).
In the clean limit, there are chiral edge states in any bulk quasienergy gap with a nonzero winding number [27]. Clearly, these edge states cannot localize when weak disorder is added. Moreover, intuitively, if all of the bulk states are localized, the chiral edge states must persist even within the spectral region of the bulk states. To see this, consider inserting a flux quantum through the hole of the cylinder. As a function of the applied flux, the chiral edge states exhibit a nontrivial "spectral flow": even though the spectrum as a whole is periodic as a function of flux, every state evolves into the next state in the spectrum [ Fig. 3(d)]. The spectral flow cannot terminate in the bulk bands since all the bulk states are localized and are, hence, insensitive to the flux. Thus, there must exist a delocalized, chiral edge state at every quasienergy to "carry" the spectral flow.
To make this argument more precise, we define a topological invariant that directly characterizes the spectral flow of the edge states. To construct the edge topological invariant, we consider a cylinder that extends from y ¼ 0 to y ¼ L y , with a flux θ x inserted through the hole of the cylinder. We use tildes to denote operators for the system on the cylinder, in particular, including the Hamiltoniañ Hðθ x ; tÞ and evolution operatorŨðθ x ; tÞ.
We now isolate the topological features of the edge states by deforming the evolution operator in the regions away from the edges, such that the evolution in the bulk takes a simple universal form, while the evolution near the edges is unaffected. In particular, we change the Hamiltoniañ Hðθ x ; tÞ away from the edges such that all of the Floquet eigenstates localized at least a distance l 0 from the edges have quasienergy ε ¼ 0. (Here, l 0 is taken to be larger than the original bulk localization length.) The resulting evolution operator, which we denote byŨ ε ðθ x ; tÞ, interpolates smoothly betweenŨðθ x ; tÞ in the vicinity of the edge and U ε ðΘ; tÞ of Eq. (7) in the bulk. Note that in the latter the value of θ y can be chosen arbitrarily [44]. We give an explicit formulation of such a deformation procedure in Appendix B. Importantly, this deformation can be performed such that the bulk Floquet states remain localized throughout the entire deformation process.
The deformed Floquet operatorŨðθ x ; TÞ takes the blockdiagonal form The sub-blocksŨ ε;1 ðθ y ; TÞ andŨ ε;2 ðθ x ; TÞ correspond to sites with 0 ≤ y ≤ l 0 and L y − l 0 ≤ y ≤ L y , respectively; the unity block acts on sites with l 0 < y < L y − l 0 . The precise value of l 0 is not important, as long as it is much larger than the bulk localization length of the original evolution operator,Ũðθ x ; TÞ. The integer-valued "edge winding number" that characterizes the spectral flow of the AFAI is defined as where in the first line of Eq. (12), and below, we use the shorthandŨ j ðTÞ ≡Ũ j;ε ðθ x ; TÞ. In the second line of Eq. (12), the sum runs over all the eigenstates of U 1 ðTÞ, and ε j are their corresponding quasienergy values.
A nonzero n edge necessarily implies that there are delocalized states along the edge; if all states were localized, their quasienergies would be almost insensitive to θ x , and, hence, n edge would be zero. Note also that, since in the AFAI all the bulk states are localized, changing l 0 would not change n edge ; this amounts to adding a few localized states to the spectrum ofŨ 1 ðTÞ, and cannot change its winding number.

V. QUANTIZED CHARGE PUMPING
We now describe the quantized nonadiabatic pumping phenomenon that characterizes the AFAI phase. Consider an AFAI placed in a cylindrical geometry, as in Fig. 3(c). Fermions are loaded into the system such that in the initial state all the lattice sites are filled up to a distance of l ≫ ξ from one edge of the cylinder, and all the other sites are empty. Here, ξ is the localization length characterizing states far from the edges of the system. Below, we show that in the thermodynamic limit the current across a vertical cut through the cylinder, averaged over many driving periods, is equal to n edge , Eq. (12), divided by the driving period T. In Appendix C 2, we furthermore show that the long-time average of the pumped charge per driving period is also given by the bulk topological invariant W, which we define in Sec. IVA. The exact form in which we terminate the filled region will not matter, as long as all the sites near one edge are filled and all the sites near the other edge are empty. The system thus serves as a quantized charge pump, but unlike the quantized pump introduced by Thouless [33], there is no requirement for adiabaticity.

A. Setup
To set up the calculation of the charge pumping in the AFAI, we choose coordinates such that x is the direction along the edges of the cylinder and y is the transverse direction. We denote the initial many-body (Slater determinant) state, in which all sites up to a distance of l from the edge are filled, by jΨð0Þi. Then, the charge pumped across the line x ¼ x 0 between t ¼ 0 and t ¼ τ is given by Here, θ x is the flux through the cylinder andHðθ x ; tÞ is the corresponding Hamiltonian. For concreteness, we specify a gauge in Eq. (13), such that on the lattice every hopping matrix element that crosses the line x ¼ x 0 has a phase of e iθ x . With this choice, the current operator across the line We are interested in the average pumped charge over N periods. Below, we show that in the limit of large N, . Note that if the system is initialized in a Slater determinant of Floquet eigenstates, then the Oð1=NÞ term in Eq. (14) is absent, since then the pumped charge per period is the same in all periods: In order to compute the charge pumped per period, it is useful to express jΨðtÞi as a superposition of Slater determinants of Floquet eigenstates. As we show in Appendix C 1, when averaging the pumped charge over N periods, the contribution of the off-diagonal terms between different Floquet eigenstates decays at least as fast as 1=N. The diagonal terms yield a contribution that depends on the evolution over a single period, giving In the above, jψ j ðtÞi are the single-particle Floquet states, which evolve in time as jψ j ðtÞi ¼ e −iεt jϕ j ðtÞi [where jϕ j ðtÞi is periodic in time], and n j are the Floquet state occupation numbers in the initial state, At this point, the average current per period depends on θ x . In the thermodynamic limit, we expect this dependence to disappear. As in the case of the quantization of the Hall conductance [46], we average over θ x [47]. We therefore get Equation (16) relates the average current in a period to the spectral flow of the Floquet spectrum as the flux θ x is threaded. It is reminiscent of the expression for the edge topological invariant n edge , Eq. (12), defined in terms of the "deformed" evolution operatorŨ ε ðTÞ. For the exactly solvable model presented in Sec. III, it is straightforward to check that the current [Eq. (16)] is indeed quantized. In that model, all bulk states have ∂ε j =∂θ x ¼ 0, whereas for all the extended states along the upper edge, ∂ε j =∂θ x ¼ 2π=NT, where N is the number of unit cells along the perimeter of the cylinder. Hence, upon summing over all states localized near one edge of the cylinder, we get that Q ∞ ¼ 1, independently of how precisely the filled region is terminated (as long as all the extended states along the edge are occupied).
Below, we give a heuristic argument that, more generally, in the AFAI phase Q ∞ ¼ n edge , up to corrections that are exponentially small in l. A more rigorous (but technically cumbersome) derivation of the relation between the pumped charge and the bulk invariant is presented in Appendix C. Numerical evidence for the quantization of the pumped charge is shown in Sec. VI.

B. Argument for Q ∞ ¼ n edge
Our strategy in analyzing the pumped charge is to deform the evolution operator into the "ideal" form, U ε ðTÞ of Eq. (11), for which the pumped charge is exactly quantized, and to put bounds on the correction to the pumped charge due to the deformation. We define the deformation process according to Appendix B, with l 0 , the width of the strip beyond which the quasienergy spectrum becomes flat, chosen such that l ∼ l 0 . Clearly, for the deformed evolution operator, n j ¼ 1 for every eigenstate ofŨ 1 . Therefore, the deformed evolution operator has an exactly quantized pumped charge, equal to n edge . Now, consider the pumped charge of the original (undeformed) evolution. We can roughly divide the Floquet states that contribute to Eq. (16) into three categories.
(1) States that are localized far from occupied region, y ≫ l. For these states, n j is exponentially small, and hence their contribution to Q ∞ is negligible. (2) States that are localized near the edge, y ≪ l. These states have n j ≈ 1. Their wave functions and quasienergies, and hence their contribution to Q ∞ , are essentially unaffected by the deformation process. (3) States that are localized near the boundary between occupied and unoccupied sites, y ∼ l. For such states, n j is neither close to 0 nor to 1; however, these states are localized in the x direction (as are all the bulk states in the AFAI). Therefore, ∂ε j =∂θ x of these states is exponentially small, and they contribute negligibly to Q ∞ . As θ x varies, there are avoided crossings in the spectrum, in which the character of the eigenstates changes. For example, an eigenstate localized around y 1 ≪ l may undergo an avoided crossing with an eigenstate localized around y 2 ∼ l. When θ x is tuned to such degeneracy points, the two eigenstates hybridize strongly and do not fall into either of the categories discussed above. Such resonances affect both ∂ε j =∂θ x and the occupations n j of the resonant states. However, since the eigenstates that cross are localized in distant spatial areas, the matrix element that couples them is exponentially small. Therefore, significant hybridization requires their energies to be tuned into resonance with exponential accuracy, limiting the regions of deviation to exponentially small ranges of θ x , of order e −l=ξ . The number of such resonances increases only polynomially with the size of the system, and, therefore, for L y ≫ l ≫ ξ and L x ∝ L y , their effect on Q ∞ is exponentially small.
We conclude that, in the thermodynamic limit, all of the contributions to Q ∞ in Eq. (16) that are not exponentially suppressed are also exponentially insensitive to the deformation process. Therefore, Q ∞ ¼ n edge .

VI. NUMERICAL RESULTS
Numerical simulations substantiate the conclusions of Secs. III-V. First, we briefly summarize our main findings, and then we describe the simulations and results in more detail in the sections below. For the simulations, we use a variant of the model discussed in Sec. III, defined on a square lattice: where H clean ðtÞ is the time-dependent, piecewise-constant Hamiltonian described in Sec. III [see Fig. 2(a)]. We define D ¼ 1 2T P r ð−1Þ η r c † r c r , and take V r to be uniformly distributed in the interval ½−δV; δV. Using numerics, we are now able to study the more generic case in which the sublattice potential (denoted here by D), as well as the disorder potential are time-independent (in contrast to the model studied in Sec. III). The parameters of the model are chosen to be λ ¼ π and δ AB ¼ 0.
In the clean case (δV ¼ 0), the system exhibits an anomalous Floquet-Bloch band structure: the Chern numbers of all the bulk bands are zero, but the winding number W ε ¼ 1 for any value of ε within each of the band gaps [27]. Such a band structure is depicted in Fig. 3(b). When the disorder potential is turned on, however, the system enters the AFAI phase. Below, we show numerically that the bulk states become localized and coexist with edge states, which occur in all quasienergies. Furthermore, when the system is initialized with fermions filling all of the sites in the vicinity of one edge, while the rest remain empty, as in Sec. V, the disordered system exhibits quantized amount of charge pumped per period, when averaged over long times. Finally, we examine the behavior of the system as the strength of the disorder potential is increased. We find that when the disorder strength reaches a certain critical value, the system undergoes a topological phase transition where the winding number changes from 1 to 0. For stronger disorder, a "trivial" phase (where all bulk states are localized and there are no chiral edge states) is stabilized.

A. Localization, edge modes, and quantized charge pumping in the AFAI
The localization properties of the bulk Floquet eigenstates of Eq. (17) can be extracted from the statistics of the spacings between the quasienergy levels. For localized states, the distribution of the level spacing is expected to have a Poissonian form. In contrast, extended states exhibit level repulsion and obey Wigner-Dyson statistics [48]. To distinguish between these distributions, it is convenient to use the ratio between the spacings of adjacent quasienergy levels [49][50][51]. Choosing the quasienergy zone to be between −π=T and π=T (i.e., choosing −i log e iεT ¼ εT for −π=T ≤ ε < π=T), we label quasienergies in ascending order. We then define the level-spacing ratio (LSR) as r ¼ minfδ n ; δ nþ1 g=maxfδ n ;δ nþ1 g, where δ n ¼ ε n − ε n−1 . This ratio, r ≤ 1, converges to different values for extended and localized states, depending on the symmetries of the system. For localized states, r loc ≈ 0.39 [49], while for extended states, r ext ≈ 0.6 [51]. The latter value is obtained when one assumes that the quasienergies are distributed according to the circular unitary ensemble [51], and, in the thermodynamic limit, coincides with the value obtained by the more familiar Gaussian unitary ensemble (GUE).
Since the Floquet problem does not possess any generic symmetries such as time-reversal, particle-hole, or chiral symmetry, we expect its localization properties to be similar to those of the unitary class [52][53][54]. In analogy with the situation in static Hamiltonians in the unitary class [55], we expect that arbitrarily weak disorder is sufficient to localize all Floquet states (on the torus). However, for weak disorder, the characteristic localization length ξ can be extremely long, and easily exceeds the system sizes accessible in our numerical simulations. Therefore, the level-spacing ratio is expected to show a gradual crossover from having the characteristic of delocalized states, r ext ≈ 0.6, when ξ ≫ L, to the value that indicates localized behavior, r ext ≈ 0.39, when ξ ≪ L.
This behavior is demonstrated in Figs. 4(a)(i)-4(a)(iii), where we plot the disorder-averaged level-spacing ratio r and the density of Floquet states as a function of the quasienergy for different disorder strengths. For weak disorder, δVT ¼ 0.5, Fig. 4(a)(i) shows that the level spacing ratio is r ≈ 0.6 in any spectral region where Floquet states exist. On the other hand, Fig. 4(a)(iii) shows that, already for δVT ¼ 4, the level-spacing ratio approaches r ≈ 0.39 at all quasienergies, as expected from localized states.
Note that, as the disorder strength increases, the levelspacing ratio decreases uniformly throughout the spectrum [Figs. 4(a)(i)-4(a)(iii); the same behavior is seen at weaker values of the disorder (not shown)]. There is no quasienergy in which the LSR remains close to 0.6, corresponding delocalized Floquet eigenstates. This is consistent with the expectation that the bulk Floquet states become localized even for weak disorder, and the localization length becomes shorter as the disorder strength increases. The behavior of the LSR as a function of system size, Fig. 4(b), also shows behavior consistent with the above expectation. In contrast, if the bulk bands of the clean systems carried nonzero Chern numbers, delocalized states would persist in the bands up to a critical strength of the disorder, at which point they would merge and annihilate.
In the AFAI phase all the bulk states are localized, but the edge hosts chiral modes at any quasienergy (cf. Sec. IV). To test this, we simulate the time evolution of wave packets initialized either in the bulk or near the edge of the system. We consider the system in a rectangular geometry. The initial state, jψ 0 i, is localized to a single site, r 0 ¼ ðx 0 ; y 0 Þ. To obtain information on quasienergy resolved propagation, we investigate the disorder-averaged transmission probability, jG N ðr; r 0 ; εÞj 2 , which is a function of both quasienergy ε and the total time of evolution T f ¼ NT.
Here, the bar denotes disorder averaging. The transmission amplitude in each disorder realization, G N , is obtained by a partial Fourier transform of the real-time amplitude, Gðr; r 0 ; tÞ ¼ hrjUðtÞjψ 0 i, and is given by The real-time transmission amplitudeGðtÞ is computed numerically by a split operator decomposition. Figures 5(a) and 5(b) show jG N j 2 at different quasienergies, for initial states on the edge and in the bulk, respectively. The simulations are done for a disorder strength VT ¼ 4. At this disorder strength, the analysis of the level-spacing statistics shown in Fig. 4(a) indicates that all the bulk Floquet bands are localized with a localization length smaller than the system size. Figure 5(a) shows the value of jG N j 2 when the wave packet is initialized at the edge of the system, r 0 ¼ ð96; 18Þ. The wave packet propagates chirally along the edge. The figure exemplifies that the edge modes are robust in the presence of disorder, and are present at all quasienergies. Importantly, edge states are also observed at quasienergies where the bulk density of states is appreciable, indicating that the chiral edge states coexist with localized bulk states [the density of states in the bulk is shown in Fig. 4(a)].
In contrast, Fig. 5(b) shows jG N j 2 for a wave packet initialized in the middle of the system. The wave packet remains localized at all quasienergies, as expected if all bulk Floquet eigenstates are localized. This confirms that the model we study numerically indeed exhibits the basic properties of the AFAI phase: fully localized Floquet bulk states, coexisting with chiral edge states that exist at every quasienergy.
Next, we numerically demonstrate the quantized charge pumping property of the AFAI. Using the model described above, we numerically compute the value ofQ ∞ given by Eq. (15) for a single value of the flux, θ x ¼ 0. When computingQ ∞ , we average the charge pumped across all the lines running parallel to the y direction of the cylinder (see Fig. 3), as well as over 100 disorder realizations. In Fig. 6(a), we show the cumulative average of the pumped charge per cycle in the limit of long times,Q ∞ [cf. Eq. (15)], as a function of disorder strength. At weak disorder, when the localization length is smaller than the system size,Q ∞ is clearly not quantized. However, for disorder strength δVT ≳ 5, the value ofQ ∞ quickly tends towards unity. This agrees with the results presented in Fig. 4(a)(iii), which indicate that at this disorder strength the localization length is substantially smaller than L ¼ 70. Finite-size scaling demonstrating thatQ ∞ indeed asymptotes to unity in the thermodynamic limit is presented in the inset of Fig. 6(a).
The value of the cumulative average of the pumped charge over N periods, hQi NT =N [cf. Eq. (13)] is plotted versus N in Fig. 6(b), demonstrating its approach toQ ∞ for large values of N (i.e., at long times). As in Fig. 6(a), we average over all the lines running parallel to the y direction, and over 100 disorder realizations. We examine the asymptotic behavior of hQi NT and find a power-law behavior of the form hQi NT ¼Q ∞ þ cN −υ with υ ¼ 1.72, shown in the inset of Fig. 6(b). Note that for a single disorder realization and a single vertical cut, hQi NT is expected to exhibit an oscillatory behavior with an envelope that decays as 1=N; see Appendix C 1. This expectation is indeed confirmed by our numerical simulations, as we show in Appendix D. In contrast, Fig. 6(b) shows a power-law behavior with a power larger than 1 and no oscillations; this is clearly the result of averaging over the frequencies appearing in hQi NT =N for each disorder realization and vertical cut. The above results numerically confirm the discussion in Sec. V and conclude our numerical analysis of the AFAI phase.

B. Strong disorder transition
For sufficiently strong disorder, we expect the AFAI to give way to a topologically trivial localized phase in which the winding number vanishes. We now analyze the transition between the AFAI and this "trivial" phase. As explained above, the winding number W ε can change only if a delocalized state crosses through the quasienergy ε as disorder is added. In the AFAI phase, all of the bulk states are already localized. How does the transition between the two phases occur?
Clearly, at the transition, delocalized states must appear in the quasienergy spectrum. As disorder is increased, the delocalized states must sweep the full quasienergy zone, changing the topological invariant W ε as they do so. The transition from the AFAI phase to the trivial phase can therefore occur through a range of disorder strength δV − c < δV < δV þ c , where δV − c is the disorder strength at which the first delocalized state appears, and δV þ c is the disorder strength at which all Floquet states are again localized, and W ε ¼ 0 for all ε. Below, we support this scenario using numerical simulations, and furthermore provide evidence suggesting that the transition is of the quantum Hall universality class.
We study the same model used in Sec. VI A and examine the level-spacing ratio r as a function of disorder strength and quasienergy. For this model, our simulations indicate δV − c ≈ δV þ c , within our resolution (limited by the system size). In Fig. 7(a), we plot r, averaged over disorder realizations and all quasienergies. We see that at disorder strength δV c T ≈ 40 the level spacing ratio reaches r ≈ 0.6, indicating delocalization. On either side of this point, r approaches 0.39 as the system size increases, which indicates localization. The peak in the value of r as a function of disorder gets sharper for larger system size, which is a signature of a critical point of this transition. In Fig. 7(b), we show that at disorder strength δV c the LSR is independent of the quasienergy with r ≈ 0.6 (for disorder strengths close to V c , we also find that the LSR is independent of the quasienergy, but with r < 0.6). This indicates that all of the Floquet states have a delocalized character at this disorder strength, which leads us to conclude that δV c ¼ δV − c ≈ δV þ c . At the critical point, δV ¼ δV c , we expect the wave functions to have a fractal character [56]. This behavior is manifested in the distribution of the inverse participation ratio (IPR), P 2 ¼ P r jψðrÞj 4 . We study the distribution of the IPR, Pðlog P 2 Þ, among all the Floquet eigenstates and averaged over disorder realizations. Figure 7(c) shows the  distribution for different system sizes. We note that the shapes of the distributions for different sizes are similar, a signature of criticality. In two dimensions, the average value of the IPR at a critical point is expected to scale like hP 2 i ∼ L −D 2 , with D 2 < 2 [56]. Figure 7(d) shows the scaling collapse of all the distributions. From the collapse, we find the fractal dimension, D 2 ¼ 1.3. The inset in this figure also shows a linear scaling loghP 2 i ∼ −D 2 log L.
The critical exponent D 2 we find in our numerical simulations is close to the value found for the universality class of quantum Hall plateau transitions [56,57], D 2 ≈ 1.4, indicating that the transition from the AFAI to the trivial phase may belong to this universality class. This is natural to expect, since, like the quantum Hall transition, in the transition out of the AFAI phase a delocalized state with a nonzero Chern number must "sweep" through every quasienergy to erase the chiral edge states. We expect that the AFAI transition can be described in terms of quantum percolation in a disordered network model, similar to the Chalker-Coddington model for the plateau transitions [58]. We leave such investigations for future work.

VII. DISCUSSION
In this paper, we demonstrate the existence of a new nonequilibrium phase of matter: the anomalous Floquet-Anderson insulator. The phase emerges in the presence of time-periodic driving and disorder in a two-dimensional system, and features a unique combination of chiral edge states and a fully localized bulk. Such a situation cannot occur in nondriven systems, where the presence of chiral edge states necessarily implies the existence of delocalized bulk states where the chiral branches of the spectrum can terminate. In a driven system, the periodicity of the quasienergy spectrum alleviates this constraint, allowing chiral states to "wrap around" the quasienergy zone and close on themselves.
One of the key physical manifestations of the AFAI is a new type of nonadiabatic quantized pumping, which occurs when all states near one edge of the system are filled. It is interesting to compare this phenomenon with Thouless's quantized adiabatic pumping, described in Ref. [33].
The complementary relationship between pumping in the AFAI and the Thouless case is best revealed by first viewing the Thouless pump from the point of view of its Floquet spectrum. In Thouless's one-dimensional pump, a periodic potential is deformed adiabatically such that in each time cycle a quantized amount of charge is pumped through the system. In the adiabatic limit, the quasienergy spectrum of the pump exhibits one pair of counterpropagating one-dimensional chiral Floquet-Bloch bands, which wrap around the quasienergy Brillouin [see Fig. 8(a)]. The (nonzero) quasienergy winding number of each band gives the associated quantized pumped charge [3]. Importantly, for any finite cycle time the two counterpropagating states hybridize and destroy the perfect quantization of the charge pumped per cycle [ Fig. 8(b)].
In a strip geometry, the AFAI can be viewed as a quasione-dimensional system. As discussed in Sec. IV, the system hosts chiral edge states that run in opposite directions on opposite edges. Furthermore, as shown by the spectral flow (see Fig. 3), these counterpropagating chiral modes cover the entire quasienergy zone, analogous to the counterpropagating modes of the Thouless pump [ Fig. 8(a)]. Crucially, however, the counterpropagating modes of the AFAI are spatially separated and, therefore, their coupling is exponentially suppressed: no adiabaticity restriction is needed to protect quantization. Thus, quantized pumping at finite frequency can be achieved in the AFAI phase.
How is the AFAI manifested in experiments? First, the localized bulk and chiral propagating edge states could be directly imaged, for example, in cold atomic or optical setups. More naturally for a solid-state electronic system, the pumping current could be monitored in a two-terminal setup. Unlike the case of an adiabatic pump, where a quantized charge is pumped at zero source-drain bias, to observe quantized charge pumping in the AFAI the chiral propagating states of one edge of the system would need to be completely filled at one end of the sample and emptied at the opposite end. We speculate that this can be achieved using a large source-drain bias. A detailed analysis of such nonequilibrium transport in a two-or multiterminal setup, as well as an investigation of promising candidate systems, are important directions for future study. While the model we present in Sec. III can be implemented directly in a variety of systems, including all-optical, microwave resonator systems and even cold-atomic systems, we expect the AFAI to be realized using other classes of models, such as those based on standard band structures with a weak uniform drive.
The implications of our results go beyond those specific to the class of systems we study in this paper. As a direct generalization of our results, one can consider constructing anomalous Floquet insulators in different dimensions and symmetry classes. Floquet-Bloch band structures that generalize those of Ref. [27] can serve as a starting point for constructing such anomalous periodically driven systems. Going beyond the single-particle level, an important challenge is to understand how the properties of the AFAI change in the presence of interactions. An exciting possibility is to obtain a topologically nontrivial steady state for an interacting, periodically driven system [15][16][17]19,20]. The common wisdom dictates that a periodically driven system with dispersive modes is doomed to evolve into a highly random state that is essentially an infinite temperature state as far as any finite-order correlation functions are concerned [51,[59][60][61][62][63]. Our results on the single-particle level demonstrate that it is possible to obtain a topological Floquet spectrum with no delocalized states away from the edges of the system. It is therefore possible that such periodically driven systems can serve as a good starting point for constructing topologically nontrivial steady states for interacting, disorder (many-body localized) periodically driven systems [64]. What types of topological steady states can be obtained by this method, and what are their observable signatures, will be interesting subjects for future work.

APPENDIX A: PERTURBATIVE DERIVATION OF THE EFFECTIVE HAMILTONIAN
Here, we outline the details of the derivation we use to demonstrate the perturbative stability of the AFAI phase. We first examine the effective static Hamiltonian defined in Eq. (2), expressed as a power series in λ. We insert Eq. (4) into Eq. (2) and expand both sides in powers of λ. For the left-hand side, we obtain The right-hand side of Eq. (2), expanded in powers of λ, reads Equating Eqs. (A1) and (A2), and using e −iTH 0 eff ¼ U 0 ðT; 0Þ, we find that Z As E r;r 0 T → 2π, Δ ðnÞ r;r 0 may diverge for all n, and the expansion in λ fails. This reflects the fact that, in general, a Floquet operator whose eigenvalues are spread throughout the quasienergy zone cannot be generated by a static, local Hamiltonian.
From the form of the right-hand side of Eq. (A3), we can analyze the maximum range of the hopping matrix elements in D ð1Þ eff . We denote the maximum range of the hopping matrix elements in DðtÞ by r, where r ¼ 1 corresponds to a nearest-neighbor hop, r ¼ 2 to second neighbors, and so on. Since the matrix elements of the unperturbed evolution operator U 0 ðt; 0Þ vanish beyond second-neighbor sites on the square lattice, we find that D ð1Þ eff contains matrix elements whose range is at most r þ 4. Similarly, from Eq. (A4), Δ ð2Þ r;r 0 vanishes beyond the 2r þ 6th neighbor, and, more generally, Δ ðnÞ r;r 0 vanishes beyond range nðr þ 2Þ þ 2. Hence, the matrix elements of D eff at range nðr þ 2Þ þ 2 contain the exponentially small factor λ n .

APPENDIX B: CONSTRUCTION OF THE DEFORMED EVOLUTION OPERATOR ON THE CYLINDER
In this appendix, we construct the deformed evolution operator on the cylinder,Ũ ε ðtÞ of Eq. (11). For brevity, we suppress the appearance of the variable θ x throughout this appendix. The deformation is designed such that at t ¼ T it interpolates smoothly betweenŨðTÞ (corresponding to the cylinder) in the vicinity of the edges of the cylinder and 1 in the bulk of the cylinder. We first define the family of operators F ðsÞ ¼ P r αðy; sÞjrihrj, where αðy; sÞ ¼ 8 > > > > > > > > < > > > > > > > > : Here, we choose ξ ≪ l 1 ≪ l 2 ≪ l 0 , where ξ is the bulk localization length. Analogously to Eq. (7), the family of deformed evolution operators is defined as U ε ðt; sÞ ¼ŨðtÞ exp½itF ðsÞH eff ε F ðsÞ: Here, H eff ε is defined as in Eq. (7), i.e., H eff ε ¼ ði=TÞ log UðTÞ, where UðTÞ is the evolution operator for a full period on the torus, and ε specifies the location of the branch cut of the log in the definition of H eff ε [see discussion below Eq. (7)]. The deformed evolution operator corresponds toŨ ε ðtÞ ≡Ũ ε ðt; s ¼ 1Þ.
Importantly, all the bulk Floquet states ofŨ ε ðT; sÞ remain localized throughout the deformation process. This is because, for a sufficiently large l 2 , the bulk Floquet states barely change as s varies; only the quasienergies change. The resulting evolution operator,Ũ ε ðtÞ, depends on the choice of ε; however, the evolution operator near the edge is essentially independent of ε. Therefore, the edge invariant, Eq. (12), does not depend on ε.
Note that, strictly speaking,Ũ ε ðtÞ of Eq. (B2) is not precisely of the block-diagonal form of Eq. (11). It still has exponentially small but nonzero matrix elements connecting the different blocks. However, a second deformation can takeŨ ε ðTÞ to the form appearing in Eq. (11).

APPENDIX C: QUANTIZED CHARGE PUMPING AND THE WINDING NUMBER
In this appendix, we show that for the AFAI a nonzero value for the winding number W ε implies quantized charge pumping. As in Sec. V, we take an initial state with all sites filled in a strip of width l near one edge of the AFAI [see Fig. 3(c)], and the rest to be empty.
To calculate the time dependence of the pumped charge, we begin by deriving an expression for the instantaneous current flowing across a longitudinal cut through the cylinder, i.e., across a line parallel to the y axis. The corresponding current operator is found by first allowing a flux θ x to be threaded through the cylinder. Next, we pick a gauge where the gauge (vector) potential is nonzero only on the links connecting sites with x ¼ L x to sites with x ¼ 0. The net current flowing across the cut between x ¼ L x and x ¼ 0 is then described by the operator I x ðtÞ ¼ ∂Hðθ x ; tÞ=∂θ x , whereHðθ x ; tÞ is the Hamiltonian of the system in the presence of the flux θ x . Here, the tilde denotes the cylindrical geometry.
Below, we first (Appendix C 1) show that, when averaged over many periods, the charge pumped approaches a quantized value Q ∞ equal to the edge topological invariant n edge , expressed in terms of the spectral flow on one edge of the system [Eq. (16)]. We then show (Appendix C 2) that Q ∞ is in fact equal to the bulk topological invariant W ε , given by Eq. (8).

Quantized charge pumping from spectral flow
We start from Eq. (13), which gives the charge pumped during the time integral 0 < t < T. The initial state, which is a single Slater determinant in terms of position eigenstates, is given by a superposition of Slater determinants in terms of Floquet states: Note that A S are fixed coefficients that depend only on the initial state. Inserting this into the expression for the pumped charge in the interval 0 < t < T, Eq. (13), we get The double sum in Eq. (C2) contains both diagonal and off-diagonal terms for the contributions of single-particle Floquet states. Denoting each of the contributions by Q jk , and using i∂ t jψðtÞi ¼HðtÞjψðtÞi, the different terms in Eq. (C2) can be written as According to Floquet's theorem, the Floquet states can be written as jψ j ðtÞi ¼ e −iε j t jϕðtÞi, where jϕ j ðtÞi ¼ jϕ j ðt þ TÞi is a periodic function. Substituting this into Eq. (C3), we obtain for the diagonal terms and for the of-diagonal terms Clearly, when computing the average charge pumped over N periods, hQi NT =N, the off-diagonal terms will give a contribution that decays as 1=N, while the diagonal terms will give the contribution that does not decay with N, where n j is the probability for the jth Floquet state to be occupied, n j ¼ hΨjψ † j ψ j jΨi. Averaging over the flux values, we obtain Eq. (16), which is the result we set out to obtain in this section.

Quantized charge pumping and the winding number
We now show that Q ∞ , the charge pumped over a period averaged over long times, is equal to the winding number W, appearing in Eq. (8). We show that in the limit of a large number of cycles N, the average charge pumped per cycle contains a quantized piece equal to the winding number plus a small correction that decays with the averaging time at least as fast as 1=N.
Recall that the many-body initial state that we consider is a single Slater determinant with electrons populating all sites in the strip of width l near one edge of the cylinder. At time t, the expectation value of the current hI x ðtÞi is given by the sum of contributions from each of these single-particle states, propagated forward in time with the evolution operatorŨ ≡Ũðθ x ; tÞ for the system with the threaded flux θ x . Defining a projector P l that projects onto all sites within the strip of initially occupied sites, the current is given by Using Eq. (C7), the total charge pumped in the first cycle, hQi T ¼ R T 0 dthI x ðtÞi, is given by Rearranging and using the chain rule, Eq. (C8) becomes In the thermodynamic limit, the current hI x ðθ x Þi is expected to be insensitive to the value of the threaded flux. Thus, replacing the pumped charge by its value averaged over all θ x , and using i∂ t U ¼ HU, we obtain or equivalently, through integration by parts, InsertingŨŨ † ¼ 1 and using ð∂ λŨ † ÞŨ ¼ −Ũ † ∂ λŨ in each of the terms in the above equation gives where we use TrfP l ½A; Bg ¼ TrfA½B; P l g. We now examine the cases for which the commutator in Eq. (C12) is nonzero. DenotingÃ ≡ ðŨ † ∂ tŨ Þ, the matrix element hrj½Ã; P l jr 0 i is nonzero in the two cases: hrj½Ã; P l jr 0 i ¼ −hrjÃjr 0 i; P l jri ¼ jri; P l jr 0 i ¼ 0; hrj½Ã; P l jr 0 i ¼ hrjÃjr 0 i; P l jri ¼ 0; P l jr 0 i ¼ jr 0 i: To set up a convenient means for enforcing the conditions above, we introduce an auxiliary gauge transformation under which the single-particle states on the sites jri ≡ jx; yi transform as jri → jri; y<l; jri → e iθ y jri; l ≤ y ≤ L y : We denote the unitary operator that applies this gauge transformation as G θ y . Because Eq. (C14) defines a pure gauge transformation, the pumped charge cannot depend on the value of θ y . Therefore, we are free to average hQi T over all such gauges. Using ½G θ y ; P l ¼ 0 and defining Importantly, Eqs. (C13) and (C14) can be expressed as whereby the pumped charge becomes Thus far, we have expressed the average current using the evolutionŨ on the cylinder. Here, we aim to obtain a bulkboundary correspondence, relating the pumped charge to the evolution operator on a torus, i.e., a geometry without edges. We consider a completion of the cylinder to a torus with fluxes θ x and θ y threaded through the two holes of the torus. The torus Hamiltonian, H ≡ Hðθ x ; θ y ; tÞ, is identical to G † θ yH ðθ x ; tÞG θ y in the interior of the cylinder. The corresponding evolution operator is denoted by U ≡ Uðθ x ; θ y ; tÞ. Importantly, U and A ¼ U † ∂ t U, as well asŨ and A ¼Ũ † ∂ tŨ , are local operators for 0 < t < T. Therefore, up to corrections that are exponentially suppressed in the size of the system, i∂ θ y A ¼ i∂ θ yÃ ðθ y Þ; i.e., the derivative with respect to θ y gives an identical result in the case of a torus and a cylinder. Moreover, since i∂ θ yÃ ðθ y Þ is also a local operator, the only matrix elements of hr 0 jBðθ y Þjri contributing in Eq. (C17) are those for which r and r 0 are in the interior of the cylinder and close to the edge of the initially filled strip, i.e., y ≈ l. For these matrix elements,B (defined on the cylinder) and B ≡ U † ∂ θ x U (defined on the torus) are identical (up to corrections that are exponentially small in the size of the system). Therefore, in Eq. (C17) we can replaceÃ andB with A and B, giving where, for brevity, we denote dθ x dθ y ¼ dΘ, and unite the integrals under a single integral sign.
Inserting UU † ¼ 1 between the two terms in the trace in the equation above, and again using ð∂ λ U † ÞU ¼ −U † ∂ λ U, we get Using we get For the moment, we focus on the second term in the above equation. Integrating by parts gives I dΘ Using the cyclic property of the trace, U † U ¼ 1, and integration by parts with respect to θ x and θ y , it is possible to show that the third and fifth terms (containing the double derivatives) cancel. The second and fourth terms, using U∂ λ U † ¼ −∂ λ UU † , can be shown to give an identical contribution to the first term in Eq. (C21), but with a factor of − 1 2 . Defining the functional W½UðtÞ for a bulk evolution UðtÞ as the net charge pumped during one driving cycle, assuming initial filling of a strip of sites covering one edge, is given by It is important to note that W½U is quantized (and equal to a winding number, as discussed in Ref. [27]) only for the case where the evolution is periodic, satisfying UðTÞ ¼ Uð0Þ. For such "ideal evolutions," the second term in Eq. (C24) clearly vanishes, and, therefore, the pumped charge is quantized and given by the winding number.
For a "nonideal evolution," where UðTÞ ≠ Uð0Þ, W½U need not be an integer. However, if the initially filled strip near the edge is wide enough such that all edge states are occupied with probabilities exponentially close to 1, then the spectral flow arguments we present in Appendix C 1 indicate that the average charge pumped per cycle will yield a quantized value, with a correction that vanishes at least as fast as 1=N. As we now show, this behavior can be seen directly through further manipulations of Eq. (C24).
Consider a "continued" evolutionÛðtÞ, defined on a larger time period of 2T. We defineÛðtÞ such that it is equal to the original evolution operator UðtÞ for 0 ≤ t ≤ T, and to e iH eff ðt−TÞ UðTÞ ¼ e iH eff ðt−2TÞ for T < t ≤ 2T. As in Ref. [27], H eff ¼ ði=TÞ log UðTÞ; in the discussion below, the choice of the branch cut of the log is unimportant, as long as the system is in a localized phase. As constructed, UðtÞ is an ideal evolution in the larger period 2T; i.e.,Ûð2TÞ ¼Ûð0Þ ¼ 1.
Starting with Eq. (C24), we add and subtract the quantity W½e iH eff ðt−TÞ , i.e., Eq. (C23) with UðtÞ replaced by e iH eff ðt−TÞ . After a shift of the time variable by T, the added piece combines with W½U to give W 2 ½Û, where W 2 is defined as in Eq. (C23) with the time integration taken from 0 to 2T. The subtracted piece remains as a correction. The charge pumped over one cycle is then where the last term arises from the full derivative (second term) in Eq. (C24). Crucially, becauseÛ is 2T periodic, W 2 ½Û is a true winding number and is quantized. The issue remains to characterize the contributions of the second and third terms; below, we show that they can be neglected in the limit of a large number of pumping cycles. Consider the average charge pumped over N driving cycles, hQi NT =N ¼ ð1=NÞ To analyze this quantity we repeat the manipulations leading up to Eq. (C24). In moving from the evolution operator on the cylinder to that on a torus, see discussion above Eq. (C18), we furthermore use the fact that in the localized phase U and A ¼ U † ∂ t U remain local even at long times. Likewise, U andÃ ¼Ũ † ∂ tŨ (for the cylinder) are local in the y coordinate. In this way we find where W N ½· is defined in the same way as W½· in Eq. (C23), but with the time integration taken up to NT rather than T. The factor fðNÞ on the right-hand side of this equation arises from the term corresponding to the full derivative term in Eq. (C24); its magnitude is bounded, and therefore the ratio fðNÞ=N decays to zero as N goes to infinity. To see how the last term in Eq. (C26) vanishes for large N, consider U in terms of its spectral decomposition, UðNTÞ ¼ P n e −iε n NT P n , where P n is the projector onto Floquet state n. Each derivative contributes two terms: U † ∂ θ j U ¼ −ið∂ θ j ε n ÞNTP n þ P n ∂ θ j P n . We consider each of the four resulting terms from the product Tr½ðU † ∂ θ x UÞðU † ∂ θ y UÞ separately.
First, when both derivatives act on the quasienergies, we get ðNTÞ 2 P n ∂ θ x ε n ∂ θ y ε n . A nonzero value for these terms would imply the existence of a current that grows linearly in time, which is unphysical. Moreover, as shown in Appendix C 1, as a general rule the time-averaged current (or pumped charge) must limit to a constant plus a correction that decreases at least inversely with time. In the fully localized phase, the quasienergies fε n g are exponentially insensitive to changes in the fluxes θ x and θ y (see arguments below), and therefore these terms clearly give a vanishing contribution in the thermodynamic limit. Therefore, these terms can (and must) be dropped within the level of all other approximations of exponential accuracy employed above.
Next, when one of the derivatives acts on the quasienergy and the other acts on a projector, we get terms like NT P n ∂ θ x ε n Tr½P n ∂ θ y P n . These terms strictly vanish due to the general identity Tr½PðdPÞ=ðdλÞ ¼ 0, for any parameter λ upon which the projector P depends.
Finally, when both derivatives act on the projectors we get a nonvanishing contribution of the form P n;m Tr½ðP n ∂ θ x P n ÞðP m ∂ θ y P m Þ. Crucially, these terms do not depend on the length of the averaging interval NT. Therefore, the quantity fðNÞ in Eq. (C26) is, in fact, constant in N, and the ratio fðNÞ=N decays to zero in the long-time (large N) limit. Furthermore, in the localized phase, the contributions from projectors onto states localized far from the bonds where the gauge fields θ x;y act are exponentially suppressed. Thus, it is clear that for any fixed N the quantity fðNÞ=N remains finite in the thermodynamic limit L x ; L y → ∞.
To evaluate W N ½U in Eq. (C26) we break up the integral over the range 0 ≤ t ≤ NT into N segments of length T. Shifting the time variable within each segment to run between 0 and T, we obtain W N ½U ¼ X N−1 n¼0 W½U n ; U n ðtÞ ¼ UðtÞUðnTÞ: As discussed for the nonideal evolution UðtÞ above, the operators fU n g are not periodic in time and therefore W½U n is not quantized.
To isolate the quantized contribution to Eq. (C26), we add and subtract a "return map" contribution W½e iH eff ðt−TÞ UðnTÞ for each term W½U n . We further define the "continued" evolutionÛ n ðtÞ ¼ÛðtÞUðnTÞ, withÛðtÞ as given above. Note thatÛ n is periodic in time with period 2T [thoughÛ n ð0Þ ¼Û n ð2TÞ ≠ 1], and therefore W 2 ½Û n is separately quantized for each n. Moreover, by virtue of the fact that the winding number W 2 ½Û n is a topological invariant for periodic evolutions, its value cannot change under smooth deformations ofÛ n . In particular, we may deformÛ n →Û via the continuous transformation U n ðt; sÞ ¼ÛðtÞUðð1 − sÞnTÞ, by taking s from 0 to 1. Hence, we find that W 2 ½Û n ¼ W 2 ½Û, and, therefore, P N−1 n¼0 W 2 ½Û n ¼ NW 2 ½Û. Inserting the result above into Eq. (C26) and subtracting the appropriate return map contribution, we obtain where we combine the contributions of the return maps for all n into one term W N ½e iH eff ðt−NTÞ . Note that by shifting time arguments we can make the replacement W N ½e iH eff ðt−NTÞ ¼ −W N ½e −iH eff t . The quantity W N ½e −iH eff t is not necessarily quantized, since the unitary e −iH eff t is not a periodic function of t over the range 0 ≤ t ≤ NT. However, if the eigenstates of H eff are all localized, then we can show (see below) that W N ½e −iH eff t decays with N as 1=N (or faster). The underlying reason is that, in the localized case, the eigenstates of H eff do not flow under insertion of the fluxes θ x and θ y into the torus. For now we simply assert this claim, and prove it at the end of this section. Accepting the claim to be true, we obtain wherefðNÞ is bounded by a constant as a function of N. Equation (C29) is the result we set out to prove in this appendix. Finally, to close the loose ends, we show that W N ½e −iH eff t is bounded as a function of N. Using the spectral decomposition e −iH eff t ¼ P n e −iε n t P n , we have dΘdt e −iΔεt ε n i8π 2 TrfP n ½P m 1 ∂ θ x P m 2 ;P k 1 ∂ θ y P k 2 g; ðC30Þ with Δε ¼ ε m 1 þ ε k 1 − ε m 2 − ε k 2 , and the sum taken over the integers n; m 1;2 ; k 1;2 . To get to Eq. (C30), we use TrfP n ½P m ; P k g ¼ 0 and TrfP n ½P m ; P k 1 ∂ θ y P k 2 g ¼ 0.
When H eff is fully localized, its eigenstates do not "wrap" around the cycles of the torus, and therefore are insusceptible to the flux insertion. Therefore, up to corrections that are suppressed as expð−L=ξÞ, where ξ is the localization length, we have that (i) the eigenvalues ε n ðΘÞ are independent of the values of the fluxes and (ii) under changing the values for the fluxes, the projectors P n transform as if transforming under a local gauge transformation, P n ðΘÞ ¼ e iϒ P n e −iϒ , with e iϒ ¼ e iQ x θ x þQ y θ y . Here, Q x and Q y are projectors on sites that define the gauge transformation felt by the localized eigenstates.
Returning to Eq. (C30), we note that in order for a term in Eq. (C30) to grow with N, it must have Δε ¼ 0.
Excluding the possibility of a fine-tuned degeneracy that occurs on a finite area in flux space, the condition Δε ¼ 0 requires that either m 1 ¼ m 2 , k 1 ¼ k 2 or m 1 ¼ k 2 , m 1 ¼ k 1 . However, the contribution of the latter two cases to W N ½e −iH eff t can be shown to vanish by substituting ∂ θ α P m ¼ ie iϒ ½Q α ; P m e −iϒ in Eq. (C30) and applying straightforward algebraic manipulations. Therefore, we find that all the terms in W N ½e −iH eff t are bounded by a constant as a function of N, and thereby we obtain Eq. (C28).

APPENDIX D: CHARGE PUMPING STATISTICS
As mentioned in the main text, quantization of the charge pumped per cycle is realized for every individual disorder realization in a large system. In this appendix, we show numerical results for a single disorder realization in a finite system of size 50 × 50 lattice sites. In the main text we define the pumped charge by integrating the current across a single vertical "cut" across the system in a cylinder geometry (a cut along the y direction, cf. Fig. 3). Here, in Fig. 9 we show that the detailed time dependence of hQi NT , the cumulative average charge pumped per cycle across each single cut, displays a unique pattern of decaying oscillations and limits to one at large times [ Fig. 9(a)]. Current conservation implies that the average currents across all cuts must be equal in the long-time limit; the spread of values at large times is due to numerical discretization error. In Fig. 9(b), we show the current averaged over all vertical cuts in the same 50 × 50 system. Here, the rapid convergence to the quantized value is clearly displayed.