Hierarchy of many-body invariants and quantized magnetization in anomalous Floquet insulators

We uncover a new family of few-body topological phases in periodically driven fermionic systems in two dimensions. These phases, which we term correlation-induced anomalous Floquet insulators (CIAFIs), are characterized by quantized contributions to the bulk magnetization from multi-particle correlations, and are classified by a family of integer-valued topological invariants. The CIAFI phases do not require many-body localization, but arise in the generic situation of k-particle localization, where the system is localized (due to disorder) for any finite number of particles up to a maximum number, k. We moreover show that, when fully many-body localized, periodically driven systems of interacting fermions in two dimensions are characterized by a quantized magnetization in the bulk, thus confirming the quantization of magnetization of the anomalous Floquet insulator. We demonstrate our results with numerical simulations.

Disorder plays a crucial role for stabilizing Floquet phases in closed systems. In particular, in the presence of interactions, disorder-induced many-body localization (MBL) provides a mechanism for the system to avoid uncontrollably absorbing energy from the driving field, and thereby to retain nontrivial properties at long times [38][39][40]. Importantly, the requirement of many-body localization does not preclude the system from exhibiting a variety of types of symmetry-breaking and topological order [25,26,37].
In this paper we characterize the topological properties of time-evolution in two-dimensional periodically driven systems of fermions which exhibit either full many-body localization, or a weaker form of "k-particle localization" that we define below [37][38][39][40] (see Fig. 1). Recent results suggest that this class of systems can support a nontrivial topological phase, known as the Anomalous Floquet Insulator [37] (AFI), which can be seen as the generalization of the AFAI to interacting systems (see Refs. 30 and 31). Despite being localized and insulating, the AFI features nontrivial circulating currents in the bulk, which in the noninteracting case (the AFAI) give rise to quantized orbital magnetization [30]. In a geometry with boundaries, the AFI supports thermalizing chiral edge states coexisting with a localized bulk [31,37]. The existence AFI as a stable many-body state of matter rests on the existence of MBL; even if MBL does hold out to infinite times, the phenomenology of the AFI is expected to persist for at least exponentially long times.
The motivation of our work is to determine the topological invariant(s) that characterize the AFI. Focusing on the topological characterization of the micromotion of particles in the bulk (i.e., the dynamics which take place within each driving period), we uncover two main results.
As our first result, we confirm that, like the AFAI, the AFI is characterized by a quantized magnetization density in regions of the bulk where all states are occupied, as schematically depicted in Fig. 1a. Specifically, the magnetization density is quantized as µ 1 /T where T denotes the driving period, and µ 1 is an integer characterizing the topological phase. This quantization is protected by many-body localization, and µ 1 cannot change under any deformation of the system that preserves MBL.
As the second major finding of our work, we uncover a rich new structure of topological invariants that emerges in the interacting case: while periodically driven systems of noninteracting fermions in two dimensions (such as the AFAI) may be characterized by a single invariant µ 1 , their interacting counterparts are characterized by a family of integer-valued topological invariants µ 1 , µ 2 , . . .. The invariant µ encodes information about the contribution to the time-averaged magnetization from -particle correlations. Hence, interactions allow for a richer topological structure in the system.
The topological protection of the invariant µ relies on a less restrictive notion of localization than the conventional notion of MBL. Specifically, µ is well-defined and topologically protected when all Floquet eigenstates with up to k particles are localized for some k ≥ l. We term this notion of localization "k-particle localization." (a) The anomalous Floquet insulator (AFI) is characterized by drive-induced circulating motion of particles in the bulk. Nontrivial topology is revealed in a quantized, nonzero magnetization density within regions where all states are filled, given by m = µ 1 T , where µ1 is a nonzero integer. (b) With sufficiently strong interactions, a new class of interaction-induced topological phases can emerge, which we term correlation-induced anomalous Floquet insulators (CIAFI's). CIAFI phases are characterized by a quantized, nonzero contribution to the magnetization from -particle correlations. Such correlations can for example arise due to immobilization of many-particle bound states, as depicted in the figure. (c,d) Topological phase transition between the AFI and a CIAFI phase with µ2 = 2 obtained from numerical simulations of a driven Hubbard-like model (see Sec. IV for details). (c) Contribution to the time-averaged magnetization in the system due to two-particle correlations, S2 (see Sec. I for definition and relationship with µ2), as a function of the interaction strength V . (d) The correlation length ξ in the system diverges for interaction strength V comparable to the hopping J, indicating a topological transition between AFI and CIAFI phases.
Many-body localization corresponds to k-particle localization in the limit where k and the system size goes to infinity, while allowing the particle density to be finite in the thermodynamic limit. While the existence of MBL in more than one dimension is still a subject of debate [41], k-particle localization for finite k is well established in any dimension [42]. It is likely that systems exhibiting k-particle localization, even if not fully MBL, may still display long-lived transient phenomena: delocalization in such systems must be induced by k+1-particle correlated processes, whose rates are expected to be exponentially suppressed in k for sufficiently weak interactions.
Our results above show that k-particle localized Floquet systems of interacting fermions in 2D are characterized by k independent topological invariants, µ 1 , . . . µ k . When one or more of the higher-order invariants are nonzero, the system is in a new, strongly-correlated, intrinsically nonequilibrium phase that is topologically distinct from any noninteracting system, including the (noninteracting) AFAI. We term this class of phases Correlation-Induced AFIs (CIAFIs). Here we consider a broader notion of the term "phase" than for equilibrium systems; in the sense we consider here, a phase charac-terizes the structure of the Hamiltonian of the isolated system, independently of the particular state of the system (and in particularly, independently of particle density and temperature).
We present a family of models which interpolate from the AFI phase to a CIAFI phase with a nonzero value of µ 2 , and demonstrate the existence of a nontrivial CIAFI phase in the model through numerical simulations [see Fig. 1 The arguments leading to the identification of the higher-order invariants µ can in principle also be applied to bosonic systems where the total number of bosons is conserved (e.g., as in systems of bosonic atoms in optical lattices). Hence AFI and CIAFI phases also exist for k-particle localized bosonic systems. However, for simplicity, in this paper, we consider fermionic systems only.
The rest of the paper is organized as follows. In Sec. I, we summarize the main results of this paper. In Sec. II we briefly review the structure of the Floquet operator in many-body and k-particle localized systems, and of the orbital magnetization operator. In Sec. III we use the time-averaged magnetization density operator to identify a set of topological invariants {µ } that characterize the AFI phase, and show that nonzero values of the invariants give rise to a quantized magnetization density in regions where all sites are occupied (Sec. III D). In Sec. IV we present a family of models that realize both the AFI and CIAFI phases, and support our conclusions with numerical simulations of these models. We conclude with a discussion in Sec. V.

I. SUMMARY OF MAIN RESULTS
We begin by summarizing the main results of this paper. We consider a two-dimensional periodically driven systems of interacting fermions, which is k-particle (or many-body) localized due to disorder [43]. To characterize the topology of the system, we quantify the circulating motion of particles in the bulk. This circulating motion can be captured through the time-averaged magnetization density operator of each plaquette p in the Heisenberg picture,m p . The magnetization densitym p measures the total time-averaged current that circulates around the plaquette; see Sec. II for a definition of this operator and a review of its properties. From its intrinsic properties, we show that the trace ofm p defines a family of topological invariants for the system. Specifically, the trace ofm p in the -particle subspace, Tr m p , for each = 1, . . . k, must take the same value for each plaquette in the system; this value cannot change under any smooth deformation of the parameters of the system that preserves k-particle localization. Hence Tr m p for each = 1, . . . k constitutes a topological invariant of the system. The intrinsic invariants µ 1 . . . µ k described in the introduction are constructed by forming system-size independent, integer-valued combinations of the (system size dependent) invariants Tr 1mp , . . . Tr kmp ; see Sec. III C for further details.
To illustrate the physical meaning of the invariants {µ }, consider first the case where the system holds a single fermion, initially located on site i in the lattice (we assume, without loss of generality, that each site holds a single orbital). When all single-particle Floquet eigenstates are localized, the particle will remain confined near site i at all times. However, the driving field may cause the particle to undergo circulating motion, as schematically depicted in the bottom left of Fig. 1(b). This circulating motion gives rise to a nonzero long-timeaveraged (orbital) moment,M i . For both single-and many-particle systems (which we consider below), the total time-averaged magnetic moment can be computed as the integral of magnetization density over the entire lattice, pm p a 2 . Ref. [31] showed that the sum ofM i over all single-particle states, S 1 ≡ iM i , is quantized as an integer times A/T , where A denotes the area of the system; this integer defines µ 1 . As an implication, magnetization density is quantized in the bulk of the system in regions where all states are occupied.
We now consider the dynamics resulting from initializing the system in a two-particle state where sites i and j are occupied. We letM ij denote the total longtime-averaged magnetization of the system resulting from this initialization. In the absence of interactions, one can verify thatM ij =M i +M j . However, with interactions present,M ij generically differs fromM i +M j when sites i and j are close to each other. The deviation can be measured by the "magnetization cumulant" C ij ≡M ij − (M i +M j ). In Sec. III below, we show that, when all 1-and 2-particle states are localized, the sum of C ij over all distinct two-particle configurations, S 2 ≡ i<j C ij , must be quantized, as an integer µ 2 times A/T . The number µ 2 cannot change under any perturbation that preserves localization of states with 1 and 2 particles. Thus, µ 2 is a topological invariant protected by 2-particle localization, and characterizes the contribution to the magnetization associated with 2-particle correlations. The higher-order invariants, µ for > 2, are defined analogously to µ 2 from higher-order "cumulants" of the magnetization (see Sec. III C for details), and µ is protected under any perturbation that preserves -particle localization. We term the class of phases characterized by nonzero values of the higher-order invariants (i.e., µ for > 1) as correlation-induced anomalous Floquet insulators (CIAFIs). The AFI phase is the MBL extension of the noninteracting AFAI, where all higher-order invariants must be zero, and can thus only be characterized by a nonzero value of µ 1 . Hence the CIAFI phases are distinct from the AFI.
In Sec. IV we present a model that realizes a CIAFI phase with µ 2 = −2. The model consists of spin-1/2 fermions on a bipartite square lattice with Hubbard-like on-site interactions and disorder, subject to the 5-step driving protocol of the canonical AFAI model [16,30,31] [see Fig. 3(a)].
As discussed in Sec. IV, and shown numerically in Fig. 1(c), the strength of the Hubbard-type interactions, V , controls the topological phase of the model [see Fig. 1(b)]: when interactions are absent (V = 0), the system is in the AFAI phase with µ 1 = 2, while all higherorder invariants take value zero [31]. When interactions are weak, but finite, our numerical results indicate that many-body localization persists, and hence the system remains in the AFI phase with µ 1 = 2 (here the factor of 2 accounts for the two spin species). In particular, the values of all higher-order invariants must remain zero [S 2 = 0, see Fig. 1(c)]. However, when interactions are much stronger than the tunneling rate between the sites, J, they act to block tunneling to or from doubly-occupied sites, resulting in nonzero values of C ij for such configurations. We demonstrate that this effect drives the model into a CIAFI phase with µ 2 = −2 (S 2 = −2A/T ). In Fig. 1(d), we confirm that the transition between the AFI and CIAFI phases in this model is accompanied by a divergence of the localization length of the two-particle states of the system.

II. MANY-BODY AND k-PARTICLE LOCALIZATION IN PERIODICALLY DRIVEN SYSTEMS
The main result of this work is to characterize the topological properties of time-evolution in two-dimensional periodically-driven k-particle (or many-body) localized fermionic systems. As a preliminary step, in this section we review the structure of the Floquet operator in such systems.
The system we study is a two-dimensional lattice systems of interacting fermions, of physical dimensions L × L, subject to periodic driving. While our results apply to any type of lattice, below we assume for simplicity that the system is defined on a square lattice with lattice constant a and (time-dependent) nearest-neighbor tunneling. The time evolution of the system is described by the time-periodic Hamiltonian H(t) = H(t + T ), where T is the driving period. To avoid complications from the coexistence of thermalizing chiral edge states and a localized bulk [37], we focus on the case where the system is defined on a torus, such that no edges are present [44].
A. Structure of Floquet operator in many-body localized systems We first review the structure of the Floquet operator when the system is many-body localized, i.e., when any state of the system exhibits localized behavior in the thermodynamic limit. The concepts we introduce here also form a basis for our discussion of the more general case of k-particle localization (Sec. II B).
When the system is MBL, it has a complete set of emergent local integrals of motion [39,40,45, 46] (LIOMs), {n a }. The LIOMs form a mutually commuting set of quasilocal operators that are individually preserved by the stroboscopic evolution of the system [47]. The number of independent LIOMs in the localized system is given by the dimension D 1 of the system's single-particle Hilbert space. For spinless fermions with one orbital per site, we have D 1 = L 2 /a 2 . The LIOMs {n α } may thus be labelled by a single index α which runs from 1 to D 1 .
To make the discussion more concrete, the LIOMs can be identified from the system's Floquet operator [39], U (T ). The Floquet operator is defined as the evolution operator of the system, U (t) ≡ T e −i t 0 dt H(t) , evaluated for a time interval corresponding to one complete driving period T . Here T denotes the time-ordering operation, and we work in units where = 1 throughout. Analogously to nondriven systems, the stroboscopic timeevolution (i.e., the time-evolution at integer multiples of the driving period T ) is conveniently expressed in terms of the eigenstates of the Floquet operator, {|ψ n }, known as Floquet eigenstates. These satisfy U (T )|ψ n = e −iεnT , where ε n has units of energy and is known as quasienergy. Note that each quasienergy ε n is only defined modulo the driving frequency Ω ≡ 2π/T . The stroboscopic timeevolution is hence equivalent to that generated by the static effective Hamiltonian, H eff ≡ n ε n |ψ n ψ n |, since In the many-body localized regime, the effective Hamiltonian takes the form (1) Each coefficient ε α1...a (referred to as a quasienergy coefficient in the following) is associated with a particular combinationn α1 . . .n α k formed from the D distinct LIOMs, and has units of energy. Each sum α1...α in Eq. (1) runs over all D combinations of distinct LIOMs, where a b denotes the binomial coefficient. The above form of the Floquet operator implies that each LIOMn α is preserved by the stroboscopic evolution of the system, and thus the operators {n α } are integrals of motion.
We now review some important properties of the LIOMs which we use in the following. Firstly, each LIOMn α can be written in the form of a fermionic counting operator:n α =f † αfα . Heref α is a (dressed) quasilocal fermionic annihilation operator, constructed from the original lattice annihilation and creation operators {ĉ i } and {ĉ † i }, respectively, as: jĉ kĉlĉm + · · · , whereĉ i annihilates a fermion on site i in the lattice. Through the identification of the LIOMs with fermionic counting operators, we note that αn α gives the total number of fermions in the system.
Another crucial property of the LIOMs is that each LIOMn α has its support localized around a particular location r α in the lattice. Specifically, the magnitude of the coefficient ψ α i1...i decreases exponentially with the distance s from any of the sites i 1 , . . . i to r α : ψ α i1...i ∼ e −s/ξ f , where the length scale ξ f sets the spatial extent of the LIOMs. Similarly to the LIOMs, the quasienergy coefficients {ε α1...α } also exhibit localized behavior. Specifically, ε α1...α decays as e −d/ξε , where d is the distance between any two of the LIOM centers r α1 . . . r α k ; here ξ ε is another localization length scale (not necessarily identical to ξ f , see Ref. 48).
As is evident above, MBL systems may be characterized by several distinct localization lengths [48]. In particular, the LIOM expansion above establishes two length scales, ξ f and ξ ε . In the following, we will make use of an additional relevant length scale, ξ l , which characterizes the spread of time-evolved operators. B. k-particle localization As we explained in the introduction, the topological classification we develop in this work applies to a more general class of systems than those exhibiting full MBL; specifically, the invariants we identify can be defined for any system that is k-particle localized for some nonzero k. As defined in the introduction, k-particle localization is understood as the situation where all Floquet eigenstates holding particles for = 1, . . . k are localized. In the remainder of this paper we will make use of similar notation, such that always refers to a specific particlenumber sector, while k refers to the "degree of localization" of the system: i.e., k is defined as the integer such that Floquet eigenstates in the system with k or fewer particles are localized, while at least one Floquet eigenstate with k + 1 particles is delocalized.
For k-particle localized systems, we expect a LIOM decomposition and effective Hamiltonian H eff as defined in Eq. (1) can be written to describe the evolution in Fock space of up to k particles, with the expansion truncated to kth order. Full MBL can be seen as a special case of k-particle localization; specifically, MBL can be understood as the k → ∞ limit of k-particle localization where the localization length of the truncated LIOM expansion described above remains bounded for all k.

III. TOPOLOGICAL INVARIANTS OF THE TIME EVOLUTION
In this section, as the main result of our work, we characterize the micromotion of k-particle localized systems (which includes the case of MBL as described above). We show that such systems may exhibit non-trivial micromotion, featuring steady-state circulating currents at long times. We characterize these circulating currents by analyzing the time-averaged magnetization density operator of the system. From this analysis we identify a set of topological invariants µ 1 . . . µ k that characterize the steady-state circulating currents that the system may support.  4)]. In many-body localized systems, the time-averaged current passing through a cut C is determined by the difference between the currents circulating around the cut's two end-points, p and q. The currents circulating around plaquette p are measured by the magnetization density operatormp. b) Ampere's law on the lattice. The difference in magnetization densities between two adjacent plaquettes p and q gives the currentĪpq on the bond between them.
In a stepwise fashion, below we consider the dynamics of a k-particle localized system in the -particle subspace for each = 1, . . . k (allowing k to be infinite for fully MBL systems). This approach ensures that our our results do not rely on full MBL to be valid, while still applying to this class of systems if such exist.

A. Characterization of micromotion
To characterize the micromotion of k-particle localized systems, in this subsection we consider the dynamics within the subspace of states holding particles, where ≤ k. Naively, one might expect that the timeaveraged current density in this subspace always vanishes due to localization. Indeed, there can be no net flow of charge across any closed curve. However, for an open curve (or "cut"), as schematically depicted in Fig. 2a, a nonzero time-averaged current may run across the cut due to uncompensated local circulating currents around the curve's endpoints. The total current circulating around a point in a given plaquette is precisely the magnetization density in this plaquette.
To establish this relationship in more rigorous terms, we consider the total time-averaged current that passes through a cut C between plaquettes p and q in the lattice, as depicted in Fig. 2a. The operator I C (t) measuring the current through the cut C is given by where I b denotes the bond current operator on bond b (restricted to the -particle subspace) [49], and the sum runs over the set B C of all bonds that cross the cut C [see Appendix A for an explicit definition of I b (t)]. Note that I b (t), and thereby I C (t), depends on time in the Schrödinger picture due to the explicit time-dependence of the Hamiltonian H(t).
To characterize the circulating currents in the system, we seek the long-time-averaged expectation value of the current I C for an arbitrary initial -particle state, |ψ . Here we introduce the notation O ≡ lim τ →∞ 1 τ τ 0 dt ψ(t)|O(t)|ψ(t) to indicate the timeaveraged expectation value in the state |ψ(t) . The timeaveraged current I C may equivalently be computed in the Heisenberg picture as I C = ψ|Ī C |ψ , where |ψ denotes the initial many-body state of the system, and I C denotes the long-time-average of the current operator I C in the Heisenberg picture: where U (t) denotes the system's time-evolution operator as defined above. For later, we define O ≡ lim τ →∞ 1 τ τ 0 dt U † (t)O(t)U (t) for any operator O. As argued above, the time-averaged currentĪ C across cut C can only have nonzero expectation value due to localized circulating currents at the cut's two endpoints, p and q. This implies thatĪ C only depends on the details of the system near plaquettes p and q. In Appendix A we verify this intuition, by proving that the opera-torĪ C only has support near the two endpoints of the cut C. Specifically, assuming only k-particle localization and conservation of charge, we show that, within theparticle subspace, where ≤ k,Ī C must take the form where the operatorm p has its full support (up to an exponentially small correction) within a distance ξ l from plaquette p, and similarly form q . Here ξ l is a finite, systemsize independent length scale measuring the spread of operators in the system (within the -particle subspace): specifically, for any time-periodic operator A(t) with a finite region of support R, the long-time averageĀ (when restricted to the -particle subspace) is a local integral of motion with support within a finite distance ξ l from R (up to an exponentially small correction) [50]. Crucially, the operatorm p in Eq. (4) is the same for any cut with an endpoint in plaquette p. Thus, Eq. (4) uniquely defines the operatorm p for each plaquette p in the system, up to a correction exponentially small in system size. Specifically, let plaquette q be separated from plaquette p by a distance d, of order the system size, L. In this case,m p can be identified uniquely from the terms ofĪ C which have support nearest to plaquette p, up to a correction of order O(e −d/ξ l ) ∼ O(e −L/ξ l ).
For each plaquette p,m p may be defined from Eq. (4) as described above by considering a cut of length ∼ L (up to an exponentially small correction). The set of operators {m p } obtained in this way then obey Eq. (4) for any two plaquettes in the lattice. In particular, when the plaquettes p and q are adjacent, Eq. (4) implies that m p −m q =Ī pq , whereĪ pq measures the time-averaged current on the bond separating plaquettes p and q, as schematically depicted in Fig. 2b. This relationship is the time-averaged lattice version of Ampere's law, which relates the current density, j, to the magnetization density, m: j = ∇ × m (see Ref. 30). We thus identify the operatorm p as the time-averaged magnetization density in the system at plaquette p [51]. As the above discussion shows, the time-averaged magnetizationm p measures the total current circulating around plaquette p.

B. Topological invariance of Tr kmp
We now show that, for each value of = 1, . . . k, the trace ofm p in the -particle subspace, Tr m p , takes the same value for all plaquettes in the system. Subsequently (in Sec. III B 1) we show that this universal value is quantized as an integer multiple of 1/T , z . Periodically driven k-particle localized systems of fermions in two dimensions are thus characterized by the k integer-valued topological invariants z 1 . . . z k .
We prove the topological invariance of Tr m p through a simple line of arguments. First, Eq. (4) implies: Using the cyclic property of the trace and Recall from Eq. (2) that the current operator I C (t) is given by a sum of bond current operators. Noting that any bond current operator I b (t) is by construction traceless (see Appendix A), we conclude that Tr Ī C = 0. Hence we find: This relation holds for any pair of plaquettes in the lattice. Therefore, for a given disorder realization, Tr m p must take the same universal value for all plaquettes in the system. We now show that the universal value of Tr m p is a topological invariant of the system in the thermodynamic limit (L → ∞) [52]. Consider perturbing H(t) within some subregion R of the system (by a small but finite amount), in such a way that -particle localization is preserved. Before and after the perturbation, Tr m p only depends on the details of the system around the plaquette p, up to an exponentially small correction (due to the exponentially decaying tails of the LIOMs). Hence, for a plaquette p located a distance of order L/2 from the region R, Tr m p0 may only change by an amount of order e −L/2ξ l due to the perturbation. Since Tr m p is given by the same value for all plaquettes in the system, Tr m p must remain unaffected by the perturbation even for plaquettes within the region where the system is perturbed, R. Thus, Tr m p is unaffected by any local perturbation that preserves -particle localization, up to a correction exponentially suppressed in system size. We conclude that Tr m p is a topological invariant of the system, protected by -particle localization.
In the following, it is convenient to parameterize the topologically-invariant value of Tr m p by a dimensionless number; we hence let z denote the value of Tr m p in units of the inverse driving period, such that Tr m p = z /T .

Quantization of z
Here we show that the dimensionless invariant z must take an integer value for each . To do this, we use an approach that generalizes the one employed for the noninteracting case in Ref. 30. This subsection provides a summary of the proof, while full details are given in Appendix B.
To begin, we consider the total time-averaged magnetization operator,M ≡ pm p a 2 . Since Tr m p takes the value z /T for all plaquettes in the system, we have To establish the quantization of z , we proceed in two steps. First, we obtain Tr M from the response of the system to the insertion of the weak uniform magnetic field B 0 = 2π/L 2 that corresponds to one flux quantum piercing the torus (note that the flux quantum is given by 2π in the units we employ): we show that, in the thermodynamic limit, whereŨ (T ) denotes the Floquet operator of the system in the presence of the magnetic field B 0 , and | · | denotes the determinant within the -particle subspace. Subsequently, we show that the determinants |Ũ | and |U | must be identical (see also Ref. 30); this implies that Tr (M )B 0 T equals an integer multiple of 2π. Using B 0 = 2π/L 2 along with Eq. (7), we conclude that z must be an integer. To obtain Eq. (8) (which forms the first step in our derivation), we show that the magnetic moment of each -particle Floquet eigenstate, |ψ n , gives the response of its quasienergy, ε n , to the addition of the weak magnetic field B 0 . Lettingε n denote the perturbed quasienergy level in the one-flux system associated with |ψ n (see the following for details, and, in particular, for a discussion of the perturbation-induced resonances), we show in Appendix B thatε Specifically, the sum ofε n − ε n over all -particle Floquet states satisfies where O(e −L/ξ ) denotes some (dimensionfull) correction which goes to zero as e −L/ξ in the thermodynamic limit. We obtain Eq. (8) from Eq. (10) by multiplying with −iT , taking the exponentials on both sides and recalling that |Ũ (T )| = exp(−i nε n T ) and likewise for U (T ).
Eq. (10) can be obtained through first-order perturbation theory in B 0 . In Appendix B, we provide a rigorous derivation of this result, along with an exact definition of the one-to-one relationship between the quasienergy levels of the one-and zero-flux systems which Eq. (10) implicitly requires. (In particular, we give the prescription for uniquely identifyingε n for each "unperturbed" quasienergy level ε n .). Here we summarize the arguments: near the region of support of |ψ n [53], the Hamiltonian of the one-flux system,H(t), is given by where θ b denotes the Peierls phase on bond b induced by the magnetic field B 0 , and I b (t) denotes the bond current operator (see Sec. III A and Appendix A). Note that there is a gauge freedom in choosing the Peierls phases; we choose them to be of order 1/L 2 near the region of support of |ψ n (such that the subleading correction in the above expansion ofH(t) can be neglected in the thermodynamic limit).
In the thermodynamic limit L → ∞, one may naively expect that the quasienergy spectrum of the one-flux system can be obtained through a first-order perturba- However, note that the convergence of such an expansion to first order is only ensured if the ratio between the matrix elements of δH in the Floquet eigenstate basis and the corresponding quasienergy level spacings, r mn ≡ ψ m |δH(t)|ψ n /(ε m − ε n ), is much smaller than 1 for all choices of -particle Floquet eigenstates m and n. While the perturbation δH(t) is of order L −2 , the many-body level spacing in the -particle subspace is of order Ω/(L 2 ), where Ω ≡ 2π/T denotes the angular driving frequency. Hence, in the thermodynamic limit r mn can potentially be much larger than 1 for certain choices of m and n. However, in Appendix B we provide a careful analysis that confirms our initial expectation: with a probability that goes to 1 in the thermodynamic limit (for each between 1 and k), r nm goes to zero for all choices of m and n. This result arises because states where ψ n |δH|ψ m is nonvanishing must be spatially close, and hence experience local level repulsion.
The above discussion shows that the quasienergy level corresponding to the state |ψ n in the one-flux system, ε n , is captured by first-order perturbation theory with respect to δH(t). Expanding the quasienergyε n to first order in δH(T ), we obtaiñ along with the fact that in a Floquet eigenstate the timeaveraged expectation value over one period is identical to the long-time average, we find whereĪ b denotes the long-time average of the bond current I b (t) in the Heisenberg picture (see Sec. III A).
Recall from Eq. (4) (see also Fig. 2b) where p b and q b denotes the two adjacent plaquettes separated by the bond b, such that b is oriented counterclockwise with respect to p b [49]. Inserting this result into Eq. (12), we note that each plaquette in the lattice appears four times exactly (namely once for each of the four bonds bounding the plaquette). Rearranging the terms from a sum over bonds to a sum over plaquettes, we thus find ε n − ε n ≈ − p ψ n |m p |ψ n (θ bp,1 + θ bp,2 + θ bp,3 + θ bp,4 ). (13) where b p,i denotes the lattice bond that constitutes the ith edge of plaquette p (counted in clockwise order starting from the positive x-direction), and θ bp i gives the Peierls phase acquired by traversing the bond counterclockwise with respect to p. The sum of Peierls phases θ bp,1 + θ bp,2 + θ bp,3 + θ bp,4 hence gives the flux through plaquette p, and hence yields exactly B 0 a 2 for each plaquette. Eq. (9) follows by usingM ≡ p a 2m p . The rigorous derivation in Appendix B shows that the correction to the approximate equality in Eq. (9) scales with system size as L −4 , and hence is subleading in thermodynamic limit (recall that B 0 ∼ L 2 ). We subsequently use the LIOM structure of the Floquet operator in Eq. (1) to show that, remarkably, these individual corrections approximately cancel out when summed over all -particle states, yielding an exponentially suppressed net correction, which scales with system size as e −L/ξ . This establishes Eq. (10), and thereby also Eq. (8).
What remains to be shown is that U (T ) andŨ (T ) have identical determinants in the -particle subspace. We show this using the approach from Ref. 30: the determinant of any time-evolution operator can be found from the time-integrated trace of the Hamiltonian [17]: which can be straightforwardly verified using the spectral decomposition of U (t). Identifying the integrand in the right-hand side above as Combining this with Eq. (7) and using that B 0 = 2π/L 2 , we conclude that z must be an integer.

C. Cumulant basis of invariants
The above discussion shows that k-particle localized systems are characterized by the k independent, integervalued topological invariants z 1 . . . z k . Here z gives the trace of the magnetization density operator in theparticle subspace (in units of the inverse driving period). However, each z depends on the size of the system, and thus is not an intrinsic property of the system. For instance, in noninteracting systems, z scales as L 2( −1) , where L is the physical dimension of the system [54]. In this subsection we construct linear combinations of the invariants z 1 . . . z k that give an equivalent set of system size independent invariants µ 1 . . . µ k that characterize the intrinsic topological properties of the system.
The intrinsic invariants µ 1 . . . µ k can be expressed as the cumulants of the magnetization operator, as discussed in Sec. I. To illustrate, consider the time-averaged magnetic moment,M ≡ p a 2m p , of a state where two particles are initialized on sites i and j, which we denotē M ij . The average of the total magnetic moment, taken over all 2-particle states, is given by 1 D2 (z 2 L 2 /T ), where D denotes the dimension of the -particle subspace. For each i and j, we writeM ij =M i +M j + C ij , where, as in Sec. I,M i denotes the time-averaged magnetization of the system holding a single particle initially located at site i. From this definition of C ij , we find where we used that Tr M = z L 2 /T for = 1, 2. The right hand side is evidently an integer multiple of 1/T . We take this integer to be our definition of the intrinsic invariant µ 2 . Note that µ 2 gives the mean value of S i ≡ j =i C ij over all sites i (recall that C ij = C ji ). Importantly, due to the fact that the two particles only influence each other's motion when they are within a localization length of one another, the cumulant C ij is only significant for O(ξ 2 l /a 2 ) choices of j for each i. The mean value of S i is therefore an intrinsic quantity, which does not depend on the system size; in particular, it remains finite in the thermodynamic limit. In the noninteracting case, C ij = 0, and µ 2 = 0. Thus, µ 2 gives the contribution to the magnetization from 2-particle correlations.
We extend this definition to higher numbers of particles, by expandingM in terms of the fermionic annihilation and creation operators, SinceM preserves the number of particles, we havē (16) Without loss of generality, we take M i1...i k ;j1...j k to be nonzero only if i 1 < i 2 . . . < i k and j 1 > j 2 . . . > j k , such that each independent combination of creation and annihilation operators appears only once in the above sum. We see that the expectation value ofM in a singleparticle state |i ≡ĉ † i |0 (where |0 denotes the vacuum state) is given by M i;i . We thus identify M ii =M i , whereM i was defined above. Likewise, in the twoparticle-state |ij ≡ĉ † iĉ † j |0 (where i < j), the expectation value ofM is given by M i;i + M j;j + M ij;ji . We thus identify M ij;ji = C ij . The higher-order cumulants can be defined in a similar fashion, such that C i1,...i = M i1...i ;i ...i1 . Note that the long-time average of an operator in the Heisenberg picture, such asM , must be diagonal in the Floquet eigenstate basis; for example, M i;j is diagonal in the basis of single-particle Floquet eigenstates.
Due to localization and the locality of interactions (see above), the coefficient C i1...i can only be nonzero if all sites i 1 . . . i are spatially close (on the scale of ξ l ). Thus, through arguments analogous to those below Eq. (15), for intrinsic quantity of the system. This motivates us to define the -th intrinsic invariant as: To relate µ to the invariants z 1 . . . z k , we take theparticle trace in Eq. (16).
(this can be verified from combinatorial arguments), where D 1 = L 2 denotes the dimension of the system's single-particle subspace, we find where we used Tr M = z L 2 /T . By induction, one can verify that each µ is an integer. First, by the definition above, µ 1 equals z 1 , and hence is an integer. For To further elucidate the physical meaning of the intrinsic invariant µ , we express it in terms of the LIOMs that were introduced in Sec. II. Since the long-time average of any Heisenberg picture operator is diagonal in the basis of Floquet eigenstates [55], the operatorm p must be an integral of motion [56]. This requiresm p to take the following form in terms of the of the LIOMs {n α } that we introduced in Eq. (1): Here, for each term involving a products of LIOMs, the sum α1...α runs over the D1 distinct combinations of LIOM indices α 1 . . . α . Due to the finite support of the operatorm p , we note that the coefficient m p α1...α vanishes as e −d/ξ l , where d is the distance from the plaquette p to the center of the most remote of the LIOMs α 1 . . . α .
Taking the -particle trace in Eq. (19) and using Tr [n α1 . . .n αν ] = D−ν −ν , we find Comparing with Eq. (18) for each = 1 . . . k, we find Note that µ is independent of the choice of plaquette p. From the expression above, it is evident that µ characterizes the intrinsic topological properties of the system. Since the magnetization coefficients {m p α1...α } vanish when the distance from any of the LIOM centers r α1 . . . r α to plaquette p becomes large, the right-hand side of Eq. (21) is independent of system size in the thermodynamic limit. In essence, µ captures the contribution of -body correlations to the magnetization density.

D. Quantized magnetization density in fully occupied regions
As a final part of this section, we show that the values of the invariants µ 1 . . . µ k can be measured directly from the magnetization density within a region of the system where all sites are occupied. In particular, for the AFI (which is fully MBL and for which only µ 1 takes nonzero value), the magnetization density is given by µ 1 /T . Consider preparing the system in an -particle state |Ψ R (where ≤ k) by filling all sites in some finite region of the lattice, R, of linear dimension d, with all sites outside R remaining empty (here we assume this requires fewer than k particles). For a plaquette p located deep within the fully occupied region, we find the time-averaged magnetization density as m p = m p R , where we introduced the shorthand O R ≡ Ψ R |O|Ψ R . Using the expansion ofm p in Eq. (19), we thus find: (22) To analyze the sum, we note that, for a LIOMn a whose center r a is located deep within the filled region R, all sites wheren a has its support are occupied. Thuŝ Here the correction arises from the exponentially decaying tail ofn α outside the filled region. For terms in the above equation where the centers of all the LIOMs α 1 . . . α ν are located near the plaquette p, the above result implies that n α1 . . .n αν R = 1 + O(e −d/ξ l ), since all of the LIOMŝ n α1 . . .n αν are located deep within the initially occupied region. For all remaining terms in Eq. (22), one or more LIOMs α 1 . . . α ν are located outside the filled region, and thus reside at least a distance ∼ d from the plaquette p. In this case, the coefficient m p α1...αν is exponentially small in d/ξ l [see the discussion below Eq. (19)]. For both categories of terms we can thus set Ψ R |m p α1...ανnα1 . . .n αν |Ψ R = m p α1...αν , at the cost of a correction of order e −d/ξ l . Doing so, we obtain Using Eq. (21), we identify the -th sum above as the invariant µ /T . Recalling that Ψ R |m p |Ψ R = m p , we thus find: The above discussion thus shows that the magnetization density deep within the filled region is given by the (convergent [58]) sum of the invariants {µ }. In particular, for the AFI, where only µ 1 is nonzero, m p = µ 1 /T . We note that the individual invariants µ 1 . . . µ k may be extracted from the dependence of the magnetization density on the particle density in the system. Specifically, for a random initial state with a uniform, finite particle density ρ, the expectation value n α1 . . .n αν , averaged over all choices of LIOMs, is given by ρ ν . Hence, at finite particle density ρ, the average magnetization density in the system is given by m p ≈ 1 T ν=1 µ ν ρ ν . The values of the individual invariants µ ν can thus be extracted from a fit of m p as a function of ρ.

IV. SPECIFIC MODEL AND NUMERICAL SIMULATIONS
In this section we present a simple model for a periodically driven system of interacting fermions in two dimensions, which realizes either the AFI or a CIAFI phase. The model was briefly discussed in Sec. I. We first consider the limit of weak interaction. In this regime we argue that the system realizes the AFI phase with µ 1 = 2. Subsequently, we show that, in the limit of strong interactions, the model is characterized by a quantized, nonzero value of the "two-particle cumulant" of the magnetization density, consistently with a CIAFI phase characterized by µ 2 = −2. To support our conclusions, we provide numerical simulations of the model in the to above regimes.
The model we consider consists of fermions with spin-1/2 living in a two-dimensional bipartite square lattice with periodic boundary conditions. The Hamiltonian is given by where H dr (t) describes piecewise-constant, timedependent hopping, H dis denotes a disorder potential, while H int describes an on-site interaction between the fermions. The driving protocol, which is contained in H dr (t), is divided into five segments, as depicted in ηT /4, while the fifth segment has duration (1 − η)T ; the parameter η is a number between 0 and 1 which controls the localization properties of the model (see below). In the first four segments, H dr (t) turns hopping on for the four different bond types in a counterclockwise fashion, as indicated in Fig. 3a, while H dr (t) = 0 in the fifth segment. More specifically, in the j-th segment (where j ≤ 4), Hereĉ r,s annihilates a fermion on site r with spin s, and the vectors {b j } are given by b 1 = −b 3 = (a, 0) and b 2 = −b 4 = (0, a). The r-sum above runs over all sites in sublattice A of the bipartite square lattice. We set the tunneling strength to J = 2π ηT , such that, in the absence of disorder and interactions, H dr would generate a perfect transfer of particles across the active bonds in each of the first four segments. The parameter η controls how rapidly the "hopping π-pulses" are applied (and thereby how strong they are relative to the disorder and interaction potentials), and thus controls the localization properties of the model; smaller η yields stronger localization (see Ref. 37).
The disorder and interaction terms H dis and H int are constant throughout the driving period and are given by For each site, w r takes a random value in the interval [−W, W ], andρ r,s ≡ĉ † r,sĉr,s denotes the occupancy on site r. The parameter V has units of energy and denotes the strength of the interactions. Note that when V J, tunneling is effectively blocked between doubly-occupied and vacant sites. As we show below, this blocking leads to a nonzero value of the higher-order invariant µ 2 .
To characterize the topological properties of the model, we consider the dynamics of particles in the two limits of weak and strong interactions. Below we demonstrate how these two regimes drive the model into the AFI phase with µ 2 = 2 and a CIAFI phase with µ 2 = −2, respectively. We substantiate these conclusions with numerical simulations in Sec. IV A.
In the absence of interactions, V = 0, the model in Eq. (24) reduces to two decoupled copies of the AFAI model from Ref. 31. When interactions are weak, but nonzero, Ref. 37 suggests that the phase remains MBL (i.e., non-thermalizing). Since the model should be connected to the non-interacting AFAI, we hence expect the system to be in the AFI phase [37] with winding number µ 1 = 2 (see also discussion in Sec. I). The factor of 2 arises from the extra species of fermions introduced due to the spin-1/2 degree of freedom.
We now show that the model above is in a CIAFI phase with µ 2 = −2 in the limit of strong interactions, V → ∞. To see this, we consider the time-averaged magnetic mo-mentM ij (see Sec. III C) that results when initially occupying two single-particle states i and j, where each choice of i or j corresponds to a particular site and spin. Recall that tunneling is blocked when the first particle is located on, or tunnels to, a site occupied by the second particle. Hence, doublons (i.e., states where two particles occupy the same site) remain frozen in place, implying thatM ij = 0 if i and j correspond to the same site being occupied. For all other initial configurations, interactions effectively do not affect the dynamics, and one can verify thatM ij =M i +M j , whereM i denotes the timeaveraged magnetic moment in the single-particle state i. As a result, the "cumulant" C ij ≡M ij −M i −M j takes value −2a 2 /T when the initialization ij corresponds to a doublon configuration, and value zero for all other 2-particle initializations (see Sec. III C for definition of C ij ). We recall from Sec. III C that µ 2 = S 2 T /L 2 , where S 2 ≡ i<j C ij T /L 2 . Since there are L 2 /a 2 distinct doublon configurations, where L denotes the physical dimension of the lattice, we find that S 2 = −2L 2 /T . Thus, µ 2 = −2 in the limit W = 0, V → ∞. From the discussion in Sec. III, we expect the quantization of µ 2 to persist for finite disorder, W , and finite (but large) values of the interaction strength, V .
The discussion above shows that the model in Eq. (24) is characterized by two distinct values of the invariant µ 2 in the limits where V = 0 and V → ∞, respectively. Due to the robust quantization of µ 2 , which is protected by 2particle localization, we hence conclude that the system supports two distinct topological phases that arise when V J and V J, respectively. The transition between the phases is separated by a critical point, V c [42]: when V is increased past V c in the thermodynamic limit, the localization length in the two-particle sector should diverge at V = V c , while µ 2 changes abruptly from 0 to −2.

A. Numerical simulations
Here we substantiate the discussion above through numerical simulations of the model: we first consider the limit of weak interactions, and show that the (quantized) average magnetic moment per particle remains unaffected by the nonzero interaction strength, as our analytical discussion predicts for an AFI phase with µ 1 = 2. Subsequently, we show that that the model is characterized by a quantized, nonzero value of the invariant µ 2 , when V is large, demonstrating that the system is in a CIAFI phase, distinct from the µ 1 = 2, µ 2 = 0 AFI phase.

Weak interactions: AFI phase with µ1 = 2
We first present data from simulations of the model described above, in the limit of weak interactions. We consider a single disorder realization of the model with parameters W = 2π/T , V = 0.1 W , and η = 1/16. From Ref. 37, we expect the model is many-body localized with these parameters. Since the model is obtained by adding weak interactions to a model of the AFAI with winding number 2 (see Refs. 30 and 31; here the factor of 2 arises because of the spin degeneracy), we moreover expect the system to be in the µ 1 = 2 AFI phase (i.e., with µ = 0 for > 1).
To probe the topology of the system, we compute the mean magnetic moments of random time-evolved 4particle states in a lattice of 6 × 6 sites. The long-time averaged magnetic moment, introduced in Sec. III, is defined asM = p a 2m p . The mean expectation value ofM , averaged over randomly chosen -particle states (i.e., states chosen randomly from a given orthonormal basis) is given by M 0 [ ] ≡ D −1 Tr M , where the binomial coefficient D counts the number of possibleparticle states in the system of D = 2L 2 single-particle states (here the factor of 2 arises due to the spin degeneracy, and L = 6 for the case we consider). Using that Tr M = z L 2 /T , along with Eq. (18), we can express M 0 [ ] in terms of the topological invariants µ 1 . . . µ : particles, our expectation that µ 1 = 2 while µ = 0 for > 1 hence would lead to corresponding to an average magnetic moment per particle of a 2 /T . This result was previously established for the noninteracting limit of the model (where the system is in the AFAI phase) [30]. The discussion above hence shows that the quantized average magnetic moment per particle in the AFAI is unaffected by interactions, as long as the system remains in the AFI phase.
To compute M 0 in the simulation, we pick as initial states 1972 random configurations of four particles located on individual sites. We evolve each initialization for 5,000 driving periods with a fixed disorder realization (the same for all initial states). Fig. 3b shows the particle density in the resulting final state for one of the realizations, after evolution for 5,000 periods. White dots and arrows indicate the corresponding initial configuration of occupied sites and spins. Note that the particle density remains non-uniform and confined near the initial location of the particles, consistent with many-body localization. We compute the time-averaged magnetic moment M for each of the 1972 states, using the time-averaged bond-currents. The 1972 values of M we obtained in this way are plotted in the histogram in Fig. 3c. Fig. 3d shows the time-averaged bond currents and magnetization density in the system for the same state used in Fig. 3b, used to calculate the magnetization. The distribution of M obtained from these initializations was found to have mean 3.999997 a 2 /T and standard deviation δM = 0.001a 2 /T , resulting in a standard deviation of the mean at δM/ √ 1972 ≈ 0.00003a 2 /T . This result is consistent with a µ 1 = 2 AFI phase [see Eq. (27)].

Strong interactions: CIAFI phase with (µ1, µ2) = (2, −2)
We now demonstrate that strong interactions drive the model into a CIAFI phase with µ 2 = −2. These data were briefly discussed in Sec. I. Here we present them in further detail.
To show that large interaction strength drives the model into the CIAFI phase, we keep W and η fixed, but vary V . We moreover consider a single disorder realization with 18×18 sites. For each value of V we consider, we obtained the time-evolution over 1000 driving periods for between 179 and 324 randomly chosen initializations where the two particles were located on particular sites and had distinct spins [59].
To establish the existence of a phase transition between the AFI and CIAFI phase, we considered the localization length in the system. We measured this using the inverse participation ratio of the density in the final state that resulted from each of the initializations we considered, P ≡ ( r |ρ r | 2 ) −1 , where ρ r = s=↑,↓ ĉ † r,sĉr,s denotes the particle density on site r in the final state. When each particle is localized on a particular site, P takes the value 1/4 (in the case of a doublon configuration) or 1/2. In contrast, P = L 2 /4 indicates full delocalization (corresponding to ρ r = 2/L 2 for all r). More generally, P can effectively be seen as 1/4 times the number of sites where the final state has support. This motivates us to define the effective localization length of the system, ξ IPR , as the average value of √ 4Pa 2 obtained from the initializations we probed.
In Fig. 1d, we plot the above localization length of the system, ξ IPR , as a function of V . As is evident in the figure, the localization length remains small for small values of V . This indicates that the µ 1 = 2 AFI phase at V = 0 remains stable for finite values of the interaction strength, as was also suggested by the results in Sec. IV A 1. In the range between V = J and V = 10J, the localization length diverges, consistent with a phase transition. For V 10J, the localization length becomes small again, indicating the system has transitioned back into a stable phase. The localization length appears to remain small as V goes to ∞; we hence expect this new phase to be the µ 2 = −2 CIAFI phase.
To verify the existence of two distinct phases (namely the µ 1 = 2 AFI and the µ 1 , µ 2 = 2, −2 CIAFI phases), we computed the sum III C or I for definition of these quantities). In Fig. 1c, we plot the value of this sum. The data shows a clear transition between µ 2 = 0 to µ 2 = −2 in the range V = J to V = 10J, where the localization length diverges. This further supports the existence of a µ 1 , µ 2 = 2, −2 CIAFI phase for strong interactions, which is distinct from the AFI phase.

V. DISCUSSION
In this work, we characterized the topological properties of periodically driven systems of interacting fermions in two dimensions. We established that the quantized magnetization of the AFAI persists in its interacting generalization, the anomalous Floquet insulator (AFIs). As a second result, we identified a new class of intrinsically-correlated nonequilibrium phases, namely the correlation-induced anomalous Floquet insulators (CIAFIs). The topological invariants characterizing the CIAFIs are encoded in the multi-particle correlations of the time-averaged magnetization density. While this work focused on driven fermionic models and their bulk topological invariants, our discussion can be readily extended to bosonic systems with particle number conservation.
Importantly, the topological protection of the CIAFIs does not require full many-body localization, but rather relies on k-particle localization, where the system is localized for any finite number of particles up to a maximum number, k. The existence of k-particle localization is well-established [42]. Since the existence of the CIAFI does not rely on full many-body localization, we may expect the behavior described above to be manifested via experimental signatures in the prethermal dynamics of systems which eventually thermalize at long times. Searching for other models that give rise to nontrivial values of these invariants and characterizing the physical properties that they imply will be interesting directions for future studies.
We demonstrated that CIAFIs may be realized in a tight-binding model with Hubbard type-interactions subject to a stepwise driving protocol. Recently, a noninter-acting version of such a model was experimentally realized with ultracold atoms in optical lattices [60]. The CIAFI phases may be achieved in a similar experiment by adding Hubbard-type interactions to the system. We expect this type of interactions is natural to implement with ultracold atoms in optical lattices. Thus, we speculate that experimental realization of CIAFI phases is feasible with current experimental platforms.
At this point it is not clear whether the CIAFI phases are compatible with MBL, i.e., if they can exist in the thermodynamic limit of L → ∞ and k → ∞. (For finite k, localization is possible, and the physics described above is rigorously applicable.) In particular, we expect that CIAFI phases will exhibit dynamics strongly dependent on the initial state. In the model of Sec. IV, initial states where some large region R is doubly occupied would support chiral edge states moving around such regions. If the initial state contains such "internal edges," they may thermalize and serve as a weak heat bath for the remainder of the system. Next, if the density of filled regions R in the system is increased, we expect that at some point thermalizing internal edges will form a connected network, destroying localization. In contrast, initial states without filled, connected regions are expected to be much more stable, since there are no direct thermalization processes which involve few nearby particles; thermalization, if it occurs at all, will proceed either due to rare thermal inclusions, or due to multi-particle tunneling into, e.g., a state with "internal edges." After the initial posting of this work, another preprint independently classified the bulk topological properties of two-dimensional MBL systems, when particle number conservation was present [61]. Interestingly the classification in Ref. 61 did not contain the CIAFI phases, suggesting that CIAFI phases and MBL may be incompatible. A definite answer for this question, however, remains lacking, and will be an interesting direction for future studies. In any case, the features above suggests that CIAFI phases (rigorously established for finite particle number) may provide a versatile playground for studying the interplay of weak thermalizing baths and MBL regions, which is expected to give new insights into the stability of MBL in 2d.
The topological classification we developed in the present work relied on particle number conservation. Chiral phases of spins and bosons without particle number conservation, which are close relatives of the AFAI (with higher-order invariants being zero, µ = 0, > 2), were considered in Ref. 29. It was shown that, when many-body localized, such phases are characterized by a quantized topological index which describes the pumping of quantum information along the edge over one driving period. Such an index arises from the rigorous classification of anomalous local unitary operators in onedimensional systems, developed by Gross et al [62]. It will be an interesting direction of future studies to investigate whether the bulk classification of the present work can be generalized to systems where particle conservation is not present.
In the future, it will moreover be interesting to investigate how thermalization is manifested in experimentally realistic situations for the CIAFI phases, and what the corresponding time scales are. With k-particle localization present (for some large k), thermalization must be driven by correlated processes involving more than k particles. It is natural to expect that such thermalizing process will be parametrically slow, and therefore signatures of the CIAFI phases (and the AFI), such as quantization of magnetization, would be observable even if MBL is eventually destroyed. A systematic study of such thermalization timescales will be an interesting question for future studies, with significance beyond the context of topological phases we considered here.
alizations where resonances occur between two or more sites separated by a distance comparable to the system size, L. Disorder realizations supporting such accidental resonances do not meet the conditions for k-particle or many-body localization, as defined in Sec. II. However, for a randomly chosen disorder realization within the k-particle localized region of parameter space, the probability that the -particle quasienergy spectrum (for each ≤ k) features any such an accidental resonance goes to zero in the thermodynamic limit L → ∞ [42]. In the following, we assume that the disorder realization under consideration does not feature such accidental resonances; within the k-particle localized regime of parameter space, this assumption holds with probability 1 in the thermodynamic limit.
[44] As for the k = 1 (i.e., single-particle) special case [30,31], we expect that k-particle localization in the bulk can coexist with delocalized edge states [42]. A detailed study of the interplay between bulk localization and delocalized edge states in the case of full MBL is left for future work; some aspects have been discussed in Refs. 37  [52] For a finite system, the fact that m p α 1 ...α is exponentially insensitive to the details of the system far away from the plaquette p means that it may only change by an amount of order e −L/ξ when the system size is increased. This implies that the sum Σα 1 ...α m p α 1 ...α is given by its value in the thermodynamic limit, up to a correction of order e −L/ξ .
[53] Here the region of support is understood as the region of the lattice where the particle density is significant in the state |ψn . See Appendix B for further details. [57] To see this, note thatnα|ΨR = (1 − fαf † α )|ΨR . The operator f † α is a polynomial in {cα} and {c † α }, where each term has the net effect of creating one fermion in the region around LIOM a. Since all sites near the LIOM a are occupied for the state |ΨR , f † α |ΨR = 0, and thuŝ nα|ΨR = |ΨR .
[58] To see that the sum in Eq. (23) converges, note that the coefficient mα 1 ...α is exponentially suppressed in d/ξ, where d is the distance from any of the LIOM centers rα 1 . . . rα to the plaquette p. The number of distinct LIOMs whose centers are located within a radius ξ l from the plaquette p is of order ξ 2 l /a 2 , where a is the lattice constant in the system. Therefore, the coefficient m p α 1 ...α vanishes exponentially when ξ 2 l /a 2 . Recalling that µ ≡ Σα 1 ...α m p α 1 ...α must take integer value for each , we thus conclude that µ equals zero when ξ 2 l /a 2 .
[59] More precisely, the initializations were divided into three classes, that contained initializations where the two particles were located on the same site, adjacent sites, and all other "on-site" initializations, respectively. The average magnetization was found through the sum of the obtained mean values of M within each class, weighted according to the number of states in the class. In this appendix we establish that the time-averaged current that passes through a cut C between two plaquettes p and q is determined by two quasilocal operators, m p andm q , with support centered at p and q, respectively [see Eq. (4) and Fig. 4]. By considering two plaquettes separated by a distance much longer than the localization length, this provides a prescription for uniquely identifying the magnetization density operatorm p (up FIG. 4. a) Schematic depiction of the argument showing that time-averaged current through a cut C between to plaquettes p and q only depends on the cut's two end-points. Specifically, since there can be no accumulation of charge over time in the region between the cuts C and C , the same current must pass through the two cuts, and thusĪC =Ī C for any two cuts C and C between the plaquettes p and q. b) The vanishing divergence of current implies thatĪC pq +ĪC qr =ĪC pr .
to exponentially small corrections in the distance, which can be of order the system size).
We recall from the main text that the operator corresponding to current through the cut C is given by where I b denotes the bond current operator on bond b, and the sum runs over all bonds that cross the cut C.
The goal of this Appendix is to find the time-averaged expectation value of the current, I C , resulting from some given initial state |ψ . As in the main text, we use O ≡ lim τ →∞ 1 τ τ 0 dt ψ(t)|O(t)|ψ(t) . The timeaveraged expectation value of the current I C may equivalently be computed in the Heisenberg picture as I C = ψ|Ī C |ψ , where |ψ denotes the initial state of the system. Here, as in the main text, for any Schrödinger picture operator O(t) [such as I C (t)],Ō denotes the timeaverage of the current I C in the Heisenberg picture, The time-averaged current operatorĪ C is thus obtained by transforming the time-dependent operator I C (t) in Eq. (A1) with evolution operator U (t), and integrating over time as in Eq. (A2). To explore the properties ofĪ C , we consider the timeaveraged current for a different cut, C , between the same two plaquettes p and q, see Fig. 4a. We note that I C (t) − I C (t) =Ṅ R (t), where N R measures the number of particles in the region R between cut C and C (shaded region in Fig. 4). Importantly, since N R is bounded by the number of sites in the region R, the long-timeaveraged value of Ṅ R must vanish. We thus conclude that I C = I C . Since this holds for any initial state |ψ , we conclude thatĪ As a next step, we note from Eqs. (A1)thatĪ C = b∈B CĪ b , whereĪ b denotes the time-averaged current on bond b [see Eq. (A2)]. We note that the operator I b (t) is local, with support only on the sites connected by the bond b. For many-body localized systems, this implies that the operatorĪ b is a localized integral of motion, with support within a distance ∼ ξ l from the bond b, up to an exponentially small correction [50]. Hence, I C is given by a sum of terms, each of which only has support within a region of radius ξ l , centered at a point along the cut C.
The requirements thatĪ C is given by a sum of local terms as described above, while at the same time taking the same value for all cuts between plaquettes p and q [Eq. (A3)], significantly constrains the form thatĪ C can take. In particular, this implies thatĪ C = I(p, q), where the operator I(p, q) only depends on the locations of the two plaquettes p and q (and not on the details of the cut C). Moreover, for any cut between plaquettes p and q, I(p, q) is given by a sum of terms which only have support in a region of width ξ l around the cut. For any site located a distance larger than ξ l from both plaquettes p or q, we can find a cut that remains separated from the site by a distance larger than ξ l . Therefore the support of operator I(p, q) can only include sites within a localization length of the endpoints p and q. Hence, we write: where A 1 (p, q) has its full support within a region of width ξ l around plaquette p, and A 2 (p, q) has support around plaquette q. The operators A 1 (p, q) and A 2 (p, q) depend only on the locations of plaquettes p and q, respectively. By letting the cut from p to q go through an arbitrary plaquette r on the torus (see Fig. 2b), we conclude from the arguments above the I(p, r) + I(r, q) = I(p, q). This implies A 1 (p, r)+A 2 (p, r)+A 1 (r, q)+A 2 (r, q) = A 1 (p, q)+A 2 (p, q).
(A5) The only terms on the left hand side with support near plaquette r are the terms A 2 (p, r), and A 2 (r, q), while none of the terms on the right-hand side have support near plaquette r. We thus conclude that A 2 (p, r) = −A 1 (r, q) for any choice of two plaquettes p and q. Hence we may write A 1 (r, q) = A(r), and A 2 (p, r) = −A(r) for some function A(r) which only depends on the location of plaquette r and has its full support near plaquette r. Using this in Eq. (A4), we find Identifying A(p) =m p , we thus conclude that Eq. (4) holds.
Appendix B: Derivation of Eq. (10) Here we derive Eq. (10), which is used to establish the integer quantization of the topological invariant z .
To recapitulate, we consider a k-particle localized system, where k may be infinite in the case of full MBL. For a given ≤ k, we consider the -particle Floquet eigenstates of the system, {|ψ n }, with corresponding quasienergies {ε n }, and letε n denote the perturbed quasienergy corresponding to ε n when the weak uniform magnetic field B 0 = 2π/L 2 is inserted that results in one flux quantum piercing the torus (see below for details). The goal of this Appendix is to establish two results. First, we show that for each -particle Floquet eigenstate, |ψ n ,ε n = ε n − B 0 ψ n |M |ψ n + O(L −5/2 ). (B1) Here, and in the remainder of this Appendix, O(L −p ) indicates a correction which goes to zero at least as fast as L −p [64]. (I.e, in the following, we only indicate how rapidly corrections decrease with system size.) Secondly, we show that, when summed over all -particle Floquet states, the corrections of order L −5/2 in Eq. (B1) approximately cancel out, yielding a net correction which is exponentially suppressed in system size: where O(e −L/ξ ) likewise indicates a correction that goes to zero as e −L/ξ in the thermodynamic limit. Eqs. (B1) and (B2) implicitly require that, for each quasienergy level ε n of the (unperturbed) zero-flux system, it should be possible to identify a unique quasienergy levelε n of the (perturbed) one-flux system which satisfies Eq. (B1). In Sec. B 4 below, we confirm that such a complete one-to-one identification is possible for all but a set of disorder realizations which has measure zero in the thermodynamic limit.
As noted in the main text, Eq. (B2) does not follow trivially from first-order perturbation theory in the weak magnetic field B 0 : under a continuous perturbation of the system, the system's quasienergy spectrum undergoes exponentially many avoided crossings due to resonances between many-body Floquet eigenstates separated by a large distance in Fock space. Hence, first-order perturbation theory breaks down for the system. Instead, we establish Eq. (10) with an alternative approach, using the localization properties of the many-particle Floquet eigenstates.
In order to follow this approach, we use a succession of auxiliary results which are not discussed in detail in the main text, but are crucial for the proof of Eqs. (B1) and (B2). The line of arguments proceeds as follows: we first show explicitly how the uniform magnetic field B 0 can be implemented in the system (Sec. B 1). Subsequently, in Sec. B 2 we show that, for a given finite region S of the lattice, it is always possible to choose a gauge where the HamiltonianH of the one-flux system resembles the Hamiltonian H of the zeroflux system locally within S, and likewise for the Floquet operatorsŨ and U (Sec. B 3). Using this result, we demonstrate in Sec. B 4 that the Floquet eigenstates and quasienergies, {|ψ n } and {ε n }, are robust to the perturbation caused by inserting of the weak uniform magnetic field B 0 , such that the one-to-one identification described above is possible. From these auxilliary results, we prove Eq. (B1) in Sec. B 5, and finally use Eq. (B1) along with the LIOM structure of the system to establish Eq. (B2) (Sec. B 6).
For the sake of brevity, throughout this Appendix we will work with a fixed degree of localization and particle number, unless otherwise noted. Thus, in the following, k and are fixed constants that refer to the system's degree of localization and to the number of particles in the system, respectively. We take ≤ k in the discussion below.

Implementation of magnetic flux
Here we discuss how the magnetic flux is implemented. The system we consider consists of interacting fermions on a lattice with the geometry of a torus, of dimensions L × L. The Hamiltonian of the system (in the absence of a flux) takes the form whereĉ i annihilates a fermion on site i in the lattice.
Here the first term contains both hopping and on-site potentials, including disorder, with J ij (t) = J * ji (t), while the term H int accounts for interactions. We allow both parts of the Hamiltonian to be time-dependent, with periodicity T . To simplify the discussion, we consider the case of a square lattice model with nearest-neighbour hoppings, and a density-density interaction described by H int = 1 2 i,jρ iρj V ij (t), whereρ i =ĉ † iĉ j and V ij (t) = V ji (t) is real. In the general case of a quasilocal Hamiltonian, the results below can also be derived using similar arguments.
In this subsection we are interested in finding the HamiltonianH(t) of the system when the uniform magnetic field B 0 = 2π L 2 is inserted, corresponding to one flux quantum through the surface of torus. Having assumed H int (t) to consist of density-density interactions, only the first term in Eq. (B3) is affected by the magnetic flux. The HamiltonianH(t) thus takes the form: Here, the Peierls phases {θ ij }, with θ ij = −θ ji , must ensure that the total phase acquired by traversing a closed loop on the torus is given by B 0 A S (mod 2π), where A S is the area enclosed by the loop [65].
There are (infinitely) many distinct configurations of the phases {θ ij } that satisfy this condition, corresponding to different choices of gauge for the one-flux Hamilto-nianH(t). As the starting point for the following discussion, we consider the following Landau-type gauge: let θ x i denote the Peierls phase for hopping along the bond in the positive x-direction from site i (and similarly let θ y i be the Peierls phase for hopping in the positive y-direction), and give them the values: Here x i and y i denote the coordinates of site i (defined with branch cut outside S 0 ), and δ ij denotes the Kronecker delta symbol, such that δ xi,L takes the value 1 if x i = L, while δ xi,L = 0 for all other values of x i . Recall that a is the lattice constant. The phases θ y i ensure that a trajectory encircling a plaquette acquires a phase of B 0 a 2 , if the trajectory does not cross the branch cut of the x-position operator between x = L and x = 0. The phase θ x i , which does not appear in the Landau gauge in an open geometry, is necessary to ensure that the phase is also given by B 0 a 2 (mod 2π) for trajectories encircling plaquettes across the branch cut.
The goal of the following is to show that we can choose another gauge where B 0 only weakly perturbs the Hamiltonian within a particular finite region of the lattice, S, which consists of one or more non-overlapping diskshaped regions, S 1 , . . . S N , whose combined area, A S , is much smaller than L 2 . We reach such a gauge through the following transformation to the one-flux Hamiltonian with the gauge choice as prescribed in Eq. (B4): 0 y i for sites i within subregion S n , and (x (n) 0 , y (n) 0 ) denotes the center of subregion S n . In this case, one can verify that, for sites within subregion n the Peierls phases resulting from this transformation take the following values:, The later holds since the branch cut of the x-coordinate does not intersect S. Since S n has disk geometry and is centered around (x A S for sites i within subregion S n . Hence we confirm that the Peierls phases are all of order √ A S a/L 2 for bonds within S, and thus much smaller than 1 in the limit A S L 2 specified above.

Response of the Hamiltonian
An important result we will use extensively in the following is that, for large systems, the insertion of the uniform field B 0 only weakly perturbs the system, up to a gauge transformation. To see this, we consider the action of the perturbation induced by B 0 , δH(t) ≡H(t) − H(t) (in the particular gauge we consider), on a state |ψ with an arbitrary number of particles, where all particles are located in the finite region S that was introduced in the previous subsection.
As a first step, we note that δH(t)|ψ = δH(t)P S |ψ , where P S projects into the subspace where all particles are located within S. Using thatĉ i P S = 0 if site i is located outside S, we find The Peierls phases {θ ij } are as given in Eq. (B5) above. Below, we establish an upper bound for the spectral norm [66] of δH(t)P S , δH(t)P S . To do this, we use that M ≤ Tr(M † M ), such that where K ij ≡ J ij (e iθij − 1), and we suppressed timedependence for brevity. Since θ ij = 0 for i = j, terms above are only nonzero when i 1 = i 2 and j 1 = j 2 . Thus, We now estimate the maximal scale of the right hand side above. We recall from the discussion in the end of Subsection B 1 that the Peierls phases {θ ij }, as given in Eq. (B5), are of order √ A S a/L 2 or smaller for bonds within the region S. This implies that the value of each non-vanishing term in the sum in Eq. (B8) is of order J 2 A S a 2 /L 4 or less, where J denotes the typical scale of the (off-diagonal) tunneling coefficients {J ij }.
To estimate the number of non-vanishing terms in the sum we recall, from the assumptions made in the beginning of subsection B 1, that the tunneling coefficients J ij only couple nearest-neighbor pairs of sites in the lattice. Hence, for each choice of the index j, J ij may only be non-vanishing for four choices of the index j. These considerations show that there are only of order A R /a 2 non-vanishing terms in the sum above. Using that each non-vanishing term has norm of order J 2 A R a 2 /L 4 , we find that δHP S 2 A 2 S J 2 L −4 . Here a b indicates that a is smaller than b, or of order b. Thus we conclude that In the sense of the operator norm, the difference between the Hamiltonians with and without one flux quantum uniformly piercing the entire torus decays to zero with the inverse of the total system area, when acting on states confined to the region S, and with a judicious choice of gauge.
a. Action on a localized state Using the above result, we now show that a gauge exists where δH is small when acting on states which are not strictly confined to the region S of the lattice, but rather only exponentially localized. Specifically, we consider a state |ψ , whose full support is exponentially confined to a region S which consists of one or more diskshaped subregions of radius r, with the probability of finding a particle a distance s from the center of the nearest subregion decaying as e −s/ξ l when s > r.
To conveniently quantify the extent to which particles are confined within a subregion of the lattice, for each j = 1, 2 . . ., we let |ψ j denote the component of the wavefunction |ψ where the outermost particle is located in the distance interval between (j − 1)a and ja from the nearest subregion of S. Specifically, |ψ j ≡ (P j − P j−1 )|ψ , where P j denotes the projector onto the states where all particles are located within a distance ja from the center of the nearest subregion of S. From this definition one can verify that |ψ = ∞ j=1 |ψ j . Moreover, the using that P j P k = P min(j,k) , it follows that the components are mutually orthogonal: ψ j |ψ k = 0 for j = k. From the definitions above, the probability for finding finding a particle more than a distance ja from the center of R is given by ψ|(1 − P j )|ψ = ∞ j =j+1 ψ j |ψ j . Since the left hand side must be of order e −ja/ξ for ja > r, and each term in the right hand side is positive, we must have ψ j |ψ j e −ja/ξ l for j > r/a.
We now use the above result to obtain a bound for the state δH|ψ . Inserting |ψ = ∞ j=1 |ψ j , and using P j |ψ j = |ψ j one can verify that |ψ = P R |ψ + j>r/a P j |ψ j , where P S ≡ P r/a denotes the projector into the subspace where all particles are located within the region S (for convenience we assume r to be an integer multiple of the lattice constant a). Using this result along with the triangle inequality and Eq. (B10), we hence obtain: The considerations from Sec. B 2 show that we may choose a gauge forH such that δHP S JA S /L 2 , and δHP j A 2 Sj J/L 2 for any choice of j, where A Sj ∼ (ja) 2 denotes the area of the region projected into by P j . Using that j>j0 j 2 e −j/k ∼ j 2 0 e −j0/k when j 0 k, one can then verify that where A S ∼ r 2 denotes the area of the region S. Thus, since r ξ l , we find

Response of the Floquet operator
We now show that, for any region S in the lattice that consists of one or more disk-shaped subregions, it is possible to find a gauge, the Floquet operators of the one-and zero-flux systems,Ũ (T ) and U (T ), have nearly identical actions states |ψ localized within S:Ũ (T )|ψ ≈ U (T )|ψ . Here the state is said to be localized within S if the probability of finding a particle a distance s from the center of the nearest subregion os S decays as e −s/ξ l for s > r, where r denotes the radius of S.
First, we note that (U −Ũ )|ψ = (Ũ † U − 1)|ψ . This follows from the unitarity ofŨ , since for any state |Ψ , |Ψ = Ũ † |Ψ . Using thatŨ Using that |Ψ = Ũ † |Ψ along with the triangle inequality, we thus find We now use that U (t) is local at all times 0 ≤ t ≤ T , due to the finite Lieb-Robinson velocity v of the system. The locality implies that, for the state U (t)|ψ , the probability of finding a particle a distance s from the center of S decays exponentially when s r. Using the result in Eq. (B12) from the previous subsection, we thus find Using this in the inequality in Eq. (B14), we conclude Thus, (Ũ − U )|ψ JT A S /L 2 . The result in Eq. (B16) shows that, with a judicious choice of gauge, the Floquet operators of the one-and zero flux systems give nearly identical results when acting on a localized state. In this sense, the insertion of a uniform magnetic field B 0 only weakly modifies the Floquet operator for large systems.

Response of Floquet eigenstates and quasienergy spectrum
We now show that, in the subspace with k or fewer particles, the quasienergy spectrum and Floquet eigenstates of k-particle localized systems are robust to perturbations, and only weakly affected by the insertion of the uniform magnetic field B 0 .
In this subsection, it is useful to use notation that relates the quasienergies and Floquet eigenstates to the LIOM decomposition in Eq. (1) (which is valid in the subspace of up to k particles, which we consider): in the following we thus let |Ψ α1...α ≡f † α1 . . .f † α |0 denote the Floquet eigenstate of the system for which only LIOMs α 1 . . . α take value 1 (see Sec. II A for definition off † α ), and let E α1...α denote the corresponding quasienergy.
Using this cutoff length, we show below that for each finite ≤ k, where k denotes the system's degree of localization (which is infinite for MBL systems), the -particle Floquet eigenstates {|Ψ α1...α } ofŨ can be labeled such that, for each choice of LIOMs (identified by the LIOM indices α 1 . . . α ), Eq. (B17) thus shows that, in the thermodynamic limit, each eigenstate ofŨ is identical to an eigenstate of U , up to gauge transformation and a vanishingly small correction, while Eq. (B18) shows that their associated quasienergies similarly are identical up to a vanishing correction. This establishes the one-to-one correspondence of the quasienergy levels of the zero-and one-flux systems that we summarized below Eq. (B2). Due to the possibility that the field B 0 induces a resonance between two Floquet eigenstates of U , disorder realizations do exist where one (or more) of the eigenstates ofŨ is a significantly hybridized combination of two eigenstates of U . In this case, Eq. (B17) will hold for most but not all Floquet eigenstates of the system. However, as we show here, the set of disorder realization where such a resonance-induced breakdown of Eq. (B17) occurs has measure zero in the thermodynamic limit. In this way, Eqs. (B17) and (B18) hold for almost all disorder realizations, in the thermodynamic limit.
To establish Eqs. (B17) and (B18), we first consider the case = 1 (i.e., we establish the relationships for each single-particle Floquet eigenstate). Subsequently, in a stepwise fashion, we generalize this result to states with particles, for each = 2, . . . k. a. Single-particle eigenstates Here we establish the relationships in Eqs. (B17) and (B18) for the single-particle case. We assume that k-particle localization is robust to perturbations, and thusŨ also describes a k-particle localized system (we assume k ≥ 1). Thus, in particular, each single-particle eigenstate |Ψ ofŨ has its full support within a finite disk-shaped region S of linear dimension d, with the probability of finding the particle a distance s outside S decaying as e −s/ξ l .
Due to its finite region of support, each single-particle eigenstate ofŨ , |Ψ , may only overlap significantly with Floquet eigenstates whose corresponding LIOM centers are located within a distance ∼ ξ l from S. To exploit this fact, we introduce a system-size dependent length scale d ξ l , which acts as an effective length cutoff for the region of support of a LIOM. The length d must be much smaller than L, but can otherwise be taken to be arbitrarily large, as long as d/L vanishes in the thermodynamic limit. From the considerations above it follows that |Ψ only overlaps with the finite number Floquet eigenstates, |Ψ α1 . . . |Ψ α N 1 , whose LIOM centers are located within a distance d from S, (up to a correction exponentially small in d/ξ l : For the purposes of the following, it is convenient to order the indices n according to the value of the overlap, such Note that the sequence of LIOM indices α 1 . . . α N1 depends on the choice of |Ψ ; this dependence is taken to be implicit below, for the sake of brevity. We now show that |Ψ only overlaps significantly with one of the eigenstates |Ψ α1 . . . |Ψ α N 1 , while the total weight from all other eigenstates gives a negligible contribution. To show this, note that |Ψ αn and |Ψ are eigenstates of U andŨ , respectively, and hence whereẼ is the quasienergy associated with |Ψ . Since |Ψ is exponentially well localized within S, Eq. (B16) Combining these two inequalities with Eq. (B20), we find We now consider two implications of the above inequality. Firstly, Eq. (B19) implies | Ψ α1 |Ψ | 2 1/N 1 − O(e −d/ξ l ) (c.f. the labelling of the states {|Ψ αn }). Thus, Secondly, we note that, for a random choice of |Ψ , the typical spacing between the N 1 quasienergy levels {E n } is of order ∆E ∼ W/N 1 , where W denotes the width of the single-particle quasienergy spectrum (when the quasienergy spectrum has no gaps, W = 2π/T ). In this case, only one of the quasienergies {E αn } (namely E α1 ) is close enough toẼ for Eq. (B21) to allow a significant value of Ψ n |Ψ . Thus, |Ψ ≈ |Ψ 1 for a typical choice of |Ψ . We now prove that |Ψ ≈ |Ψ 1 for any choice of |Ψ in the system (except for a measure-zero set of disorder realizations in the thermodynamic limit). To establish this result, we first note We now establish a lower bound for |E αn − E α1 |, using the fact the quasienergy levels of nearby states E α1 and E αn repel each other, and that |Ẽ − E α1 | satisfies the bound of Eq. (B22). Specifically, note that the Floquet eigenstates |Ψ 1 and |Ψ n have their support within a distance d from each other. The quasienergies E α1 and E αn are hence subject to local level repulsion when the quasienergy difference δE ≡ |E n − E 1 | is much smaller than the scale of matrix elements between them with respect to the kinetic part of the Hamiltonian (i.e. δE Je −d/ξ l ). In the limit where δE Je −d/ξ l , the probability distribution p(δE) for δE should thus resemble the Wigner-Dyson distribution for the Circular unitary ensemble (CUE) [67]: Using the above result, we now compute the expected number of pairs of nearby single-particle eigenstates |Ψ αi and |Ψ αj in the entire system, for which |E αi − E αj | is smaller than some given (small) value δE 0 . Here "nearby" refers to the eigenstates |Ψ αi and |Ψ αj having their centers located within a distance ∼ d from each other, such that they may potentially overlap with the same eigenstate ofŨ . Noting that there are O(L 2 N 1 /2a 2 ) distinct pairs of nearby eigenstates (where a denotes the lattice constant), we have where we used N 1 ∼ (d/a) 2 . Thus, in the limit where δE 0 Je −d/ξ l , We recall we may take d arbitrarily large as long as d/L → 0 in the thermodynamic limit. In the following, it is convenient to let d scale with system-size as d ∼ We conclude that, in the thermodynamic limit, there are zero pairs of Floquet eigenstates |Ψ αi and |Ψ αj with LIOM centers within a distance d ∼ 1 2 ξ log(L/a) from each other whose quasienergies differ by less than a LIOMsn α andn β are located within a distance d from the centers ofñ 1 andñ 2 . As a result, there are only of order N 2 ∼ 2d 2 /a 2 2 choices of distinct LIOMs α, β for which |Ψ αβ can significantly overlap with |Ψ .
Using the same arguments as for the single particle case (Sec. B 4 a) one can show that, for all but a measure-zero set of disorder realizations in the thermodynamic limit, there exists a unique two-particle eigenstate |Ψ αβ of U for each two-particle eigenstate |Ψ ofŨ such that (up to a gauge transformation) Separated LIOMs -Next, we consider the case where the two excited LIOMsñ 1 andñ 2 are separated by a distance ∆r larger than d. In this case, the LIOM structure of the Floquet operatorŨ [Eq. (1) in the main text] implies that, up to an exponentially small correction in the distance ∆r/ξ l , |Ψ may be written as a direct product of two single-particle eigenstates |Ψ α and |Ψ β . Here α and β refer to the labeling of the single-particle eigenstates ofŨ that was established in the previous subsection. Letting S α and S β denote the two non-overlapping regions of linear dimension d where the states |Ψ α and |Ψ β respectively have their support (up to a correction exponentially small in d/ξ l ), we have [68]: |Ψ = |Ψ α Sα ⊗ |Ψ β S β ⊗ |0 + O(e −d/ξ l ). (B35) where we used ∆r > d. Here |Ψ S denotes the restriction of the state |Ψ to the Fock space of the region S (defined from the projection of |Ψ into the subspace with no particles outside region S). The state |0 refers to the vacuum in the complementary region to S α and S β . Since the two particles in the state |Ψ are separated by a distance much larger than d, the regions S α and S β do not overlap. We recall that Eq. (B17) was already proven to hold for the single-particle case. Thus |Ψ α (the eigenstate in the presence of one flux quantum piercing the system) is approximately identical to a single-particle eigenstate |Ψ α of the zero-flux system's Floquet operator U (for all but a measure zero set of disorder realizations). Specifically, up to a gauge transformation, |Ψ α = |Ψ α + O(L −2 ). The eigenstate |Ψ α moreover has its full support in the same region S α as |Ψ α , up to a correction exponentially small in d/ξ l . Letting V α be the unitary operator that generates the transformation to the gauge in which Eq. (B17) holds for |Ψ α , we have where we used that we may take d ∼ 1 2 ξ l log(L/a), such that the correction O(e −d/ξ ) scales with system size as L −1/2 in the thermodynamic limit. Using the relation (B36) for the states |Ψ α Sα and |Ψ β S β in Eq. (B35), we hence obtain |Ψ = V α V β |Ψ α Sα ⊗ |Ψ β S β ⊗ |0 + O(L −1/2 ). (B37) Due to the LIOM structure of the Floquet operator U (Eq. (1) in the main text), |Ψ α Sα ⊗ |Ψ β S β ⊗ |0 is identical to the Floquet eigenstate |Ψ αβ of the zero-flux system, up to a correction of order e −d/ξ l . Since the product of the two gauge transformations V α and V β is itself a gauge transformation, we thus conclude that, up to a gauge transformation: The two cases we considered above show that, in the thermodynamic limit, and for all but a measure zero set of disorder realizations, each two-particle eigenstate |Ψ ofŨ is identical to a unique eigenstate of U , up to a gauge transformation, and a correction of order O(L −1/2 ). We may thus label the two-particle eigenstates ofŨ such that Eqs. (B17) and (B18) hold with = 2, and for each choice of the LIOM indices α 1 and α 2 .
c. -particle-eigenstates We finally consider the general case of an -particle eigenstate |Ψ ofŨ , where is smaller than or equal to the system's degree of localization, k. For this situation, we can apply the same structure of arguments as for the two-particle case: due to the LIOM structure of the one-flux Floquet operatorŨ , each -particle state is constructed by "exciting" LIOMsñ 1 . . .ñ . We split our line of arguments into two cases, depending on whether or not the LIOMsñ 1 . . .ñ can be divided into clusters separated from each other by distances greater than d.
In the case where the excited LIOMs can be divided into clusters in the way above, |ψ can be written as a direct product of eigenstates ofŨ with fewer than k particles, up to a correction of order e −d/ξ l . Following the same line of arguments as for the analogous twoparticle case, the relationships (B17) and (B18) can then be demonstrated to hold for this class of eigenstates using the fact that Eq. (B17) and (B18) hold for eigenstates with fewer than particles.
In the case where all LIOMs are located in the same cluster, we note that |ψ only significantly overlaps with eigenstates {|Ψ α1...α } where the centers of all the LIOMsn α1 . . .n α are located in the region S, consisting of all sites with a distance d from any of the excited LIOM'sñ 1 . . .ñ . There only exist a finite number of eigenstates N with this property. Specifically, N d 2 /a 2 counts the number of distinct configurations of k LIOMsn α1 . . .n α whose centers are located within S. Crucially, N only depends on the number of particles, , and d, and is independent of system size.
Using the same arguments as for the single-particle case, we then find that, for all but a measure zero set of disorder realizations in the thermodynamic limit, there exists a unique eigenstate |Ψ α1...α of U such that (up to a gauge transformation), where, as we described in the beginning of this Appendix, O(L −p ) denotes term scaling with system size as L −p in the thermodynamic limit (see Footnote [64]). In addition, when the LIOMs are located within a distance d from the same point,Ẽ Thus, Eqs. (B17) and (B18) hold for the -particle case in the thermodynamic limit, for any = 1, . . . k.

Relationship between magnetization density and quasienergy
Having established the auxiliary results in Secs. B 1-B 4, we are now ready to prove Eq. (B1), which is the first main goal of this appendix. To recapitulate, we seek to show that, for each -particle Floquet eigenstate, |ψ n with quasienergy ε n , the associated quasienergy for the one-flux system,ε n (see Sec. B 4 for details), satisfies ε n = ε n + B 0 ψ n |M |ψ n + O(L −5/2 ), whereM denotes the time-averaged magnetization operator (see Sec. III B 1 of the main text), and, as in Sec. B 4 above, O(L −p ) denotes a correction of order λL −p or less, where λ is some system-size independent energy scale that does not play a role for our discussion. In this step of the derivation it is useful to define a region of support, S n , for each Floquet eigenstate |ψ n . Specifically, for each Floquet eigenstate, |ψ n , and for some length scale d L, we let S n denote the smallest region of the lattice that ensures the centers of all nonzero LIOMs in the state |ψ n , α 1 . . . α , are located within a distance d from the boundary of S n . The region of support S n may consist of one or several disconnected diskshaped subregions of linear dimension d, and has area A Sn ≤ π d 2 . As in Sec. B 4, when taking the thermodynamic limit L → ∞ in the following, we let d increase logarithmically with system size as d ∼ 1 2 ξ l log(L/a). To establish Eq. (B41), for a given Floquet eigenstate |ψ n , we letŨ be the one-flux Floquet operator in a gauge where Eq. (B16) holds within S n , and let |ψ n denote the eigenstate ofŨ corresponding to |ψ n through Eq. (B17). Noting that |ψ n and |ψ n are eigenstates of U andŨ , respectively, we have