Chiral Floquet Phases of Many-body Localized Bosons

We construct and classify chiral topological phases in driven (Floquet) systems of strongly interacting bosons, with finite-dimensional site Hilbert spaces, in two spatial dimensions. The construction proceeds by introducing exactly soluble models with chiral edges, which in the presence of many-body localization (MBL) in the bulk are argued to lead to stable chiral phases. These chiral phases do not require any symmetry, and in fact owe their existence to the absence of energy conservation in driven systems. Surprisingly, we show that they are classified by a quantized many-body index, which is well defined for any MBL Floquet system. The value of this index, which is always the logarithm of a positive rational number, can be interpreted as the entropy per Floquet cycle pumped along the edge, formalizing the notion of quantum-information flow. We explicitly compute this index for specific models, and show that the nontrivial topology leads to edge thermalization, which provides an interesting link between bulk topology and chaos at the edge. We also discuss chiral Floquet phases in interacting fermionic systems and their relation to chiral bosonic phases.


I. INTRODUCTION
Topological phases are typically discussed in terms of the ground state properties of gapped Hamiltonians. One of the earliest examples is the integer quantum Hall (IQH) effect , where a quantized Hall conductance is established on cooling to low temperatures that are well below the gap scale. The IQH insulator is an example of a short-range entangled (SRE) topological phase, defined as having a unique ground state on closed manifolds. [1] A well known class of SRE phases are the symmetry-protected topological (SPT) phases, which require the presence of a protecting symmetry, examples of which include the electronic topological insulators. In contrast, the IQH insulators exemplify a more basic class of 2d 'chiral' SRE phases whose nontrivial nature persists even in the absence of any symmetry. These nontrivial chiral phases are fully characterized by the existence of chiral edge modes, and can be diagnosed by their quantized thermal Hall conductance, [2] which is proportional to the chiral central charge 'c' of the effective edge conformal field theory. In particular, the chiral central charge is quantized in units of c * = 1/2 (c * = 8) for fermions (bosons) with no additional symmetry. [3][4][5] More recently it has been realized that, outside of ground state physics, topological phases can also appear in the highly excited states of a many-body system. [6][7][8] An essential ingredient in this scenario is many-body localization, which allows for a description of states in terms of conserved local integrals of motion. [9][10][11] This prevents thermalization and endows the excited states with properties similar to those of the ground states of gapped systems. [12] These distinctive properties can be observed in the dynamics of simple initial states, without the need for cooling, so that in addition to their intrinsic interest, many-body localized (MBL) topological phases also possess practical advantages in terms of realizability. Although many SPT phases have been realized in this MBL context, [6,7,13,14] the more fundamental 2d chiral topological phases have so far evaded such a realization. Indeed, it was argued in Refs. [3], [14] and [15] that such chiral phases are not permitted in excited states, because the resulting system would then be unstable to thermalization.
A key ingredient in the argument of Refs. [3], [14] and [15] is the inability of a system characterized by local integrals of motion to support a quantized thermal Hall conductance. On the other hand, it has also been understood recently that MBL systems can be stable to the introduction of a time-dependent, periodic Floquet drive, [16][17][18] and this opens up the possibility for the realization of new topological phases in MBL Floquet systems. [19][20][21][22][23] Such a periodic drive renders the thermal Hall conductance ambiguous, since the energy is not conserved, and allows for the possibility of realizing chiral MBL phases in Floquet systems.
Here we will argue that chiral phases, which are intrinsically dynamical in nature and do not have static MBL counterparts, can indeed be realized in periodically driven MBL systems. This should be differentiated from various Floquet engineering proposals, which aim at effectively realizing equilibrium phases through periodic driving. [24,25] The key observation we need is that, despite the absence of a chiral central charge or any other conserved quantity, the notion of a chiral edge is still tenable in such systems. The quantity being transported at the edge is quantum information, which is omnipresent in quantum systems and can be pumped much like conserved charge or energy. [26] Specifically, we construct a class of exactly solved spin models (termed 'bosonic chiral Floquet' models) which, despite being localized in the bulk, possess such chiral edges.
These models are bosonic analogues of the 'Anomalous Floquet-Anderson Insulator' (AFAI), introduced in the insightful works of Refs. [27][28][29][30], where unidirectional edge states appear in conjunction with an Anderson lo-calized bulk. However, the stability of these free fermion systems in the presence of interactions remains uncertain, since previous works have not identified an intrinsically many-body index that is quantized in these models; in particular, the topology of these models was linked to certain winding numbers, which are explicitly singleparticle properties. By contrast we will characterize our bosonic chiral Floquet models by an intrinsically interacting many-body topological index. The topological invariant we use was first discussed in Ref. [26] (henceforth referred to as 'GNVW'), in which the index was developed to exhaustively classify 1d quantum cellular automata. In this work, we will bridge the quantum information perspective in GNVW to the physical problems of bosonic Floquet systems in 2d. Interestingly, the topological invariant is not just an integer, as one might expect from the familiar IQH classification, but rather takes the form ν = log(p/q), where p/q is a positive rational number and depends on the dimensions of the local Hilbert spaces.
A feature of the MBL Floquet problems is the notion of a generic edge phase. In contrast to MBL Floquet SPT phases, where the edges can be localized following spontaneous symmetry breaking, the absence of symmetry requirements implies that the edge of an MBL chiral Floquet phase can never be localized. This indicates that bulk topology can enforce chaotic dynamics at the edge, which is supported by our numerical computation on a model of the chiral Floquet edge.
While our discussion will mostly focus on bosonic systems, where rigorous results are available, we will also shed some light on fermionic problems. We present strong evidence that minimal fermionic chiral Floquet phase is topologically equivalent to the bosonic counterpart. This is in stark contrast to the equilibrium quantum Hall phases, where the minimal bosonic integer chiral phase is equivalent to 8 copies of the minimal fermionic one.
The plan of attack in this paper will be as follows. We will introduce a class of exactly soluble bosonic Floquet models (the 'SWAP models'), and construct a topological model of interacting bosons that features the 1d translation operator at the edge (Sec. II). Leveraging GNVW's results on the classification of 1d locality-preserving unitary operators, we argue that this edge dynamics cannot be realized in any purely 1d system undergoing finitetime Hamiltonian evolution, and hence signifies a topological 2d bulk. Having analyzed the SWAP models, we will abstract their essential features and show that they serve as a complete set of representatives for 2d bosonic MBL chiral Floquet phases. This is achieved by first developing a sharp notion of bulk-boundary correspondence in MBL Floquet systems (Sec. III), including subtleties associated with the robustness of many-body localization in 2d, and then discussing in depth the interpretation, properties, and classification structure of the topological index in GNVW, which will be referred to as the 'chiral unitary index' and denoted by ν (Sec. IV). After addressing the formal aspects of the MBL chiral Floquet phases, we will demonstrate the explicit computability of ν via a matrix-product representation (Sec. V), discuss an experimental proposal for realization using hardcore bosons in a shaken optical lattice (Sec. VI), and elaborate on the physical consequences associated with the anomalous chiral edges of these models (Sec. VII). We will conclude by making connection between the bosonic and fermionic problems (Sec. VIII), and then discussing future directions of research motivated by the results in this work (Sec. IX).

II. BOSONIC CHIRAL FLOQUET MODELS
The goal of this section is to construct an infinite set of Floquet models in 2d bosonic spin systems, which we will refer to as 'bosonic chiral Floquet' models. For concreteness, we will first analyze in detail a model of spin-1/2's in subsection II A, in which the bulk and edge degrees of freedom (DOF) are manifestly decoupled, and the edge Floquet operator is simply given by the 1d translation operator (also frequently called a 'shift', as in GNVW). We will then generalize the model to an infinite set of bosonic chiral Floquet models in subsection II B, which show similar edge dynamics as the spin-1/2 model. We will then argue heuristically in II C that some of these models are topologically nontrivial, as their edges, the 1d translation operators, are anomalous.

A. Example of a bosonic chiral Floquet model
The model we construct can be viewed as a bosonic version of the fermionic AFAI model presented in Refs. [28] and [30], where the fermions are replaced by hardcore bosons, and each fermionic 'hop' is replaced by a bosonic SWAP gate. In anticipation of the generalizations in subsequent subsections, we we will view the hardcore bosons as spin-1/2's.
More concretely, let Λ be a checkerboard lattice with two (square) sub-lattices A and B. We will view Λ as rotated 45 degrees from the horizontal, and choose primitive lattice vectorsx,ŷ oriented in the x and y directions respectively. We can then view Λ as a crystal lattice with a two-site unit cell, with basis vectors 0 for A and (x +ŷ)/2 for B. Each site will host a spin-1/2. For concreteness, let us take a large rectangular system with size N x ×N y , so that the total number of sites is |Λ| ≡ 2N x N y .
Our Hamiltonian will be a piecewise constant function of time, with four distinct time steps, each of duration T /4. The four corresponding Hamiltonians will be denotedĤ (s) , s = 1, . . . , 4, so thatĤ(t) =Ĥ (s) for (s − 1)T /4 ≤ t ≤ sT /4. For each s, we definê By construction, all terms inĤ (s) commute and so the time-evolution operator exp(−iĤ (s) T /4) for a single time step can be readily computed: whereχ x,y is simply the 'SWAP' gate between the spin-1/2 DOF at sites x and y. Hence, the total Floquet operator can be viewed as a quantum circuit built entirely of SWAP gates, and is merely a permutation of the lattice sites. We denote it byÛ F ≡P F ≡P (4) . . .P (1) to emphasize that each step is a site permutation. Due to its permutation nature, a SWAP circuit is exactly soluble, as we discuss in Appendix A. In particular we will see below that with periodic boundary conditions (PBC), U F is just the identity.
To compute the Floquet operator, which is a site permutation, it is enough to work in the 'single-spin-flip sector', defined as the subspace spanned by the orthonormal set of single-spin-flip states {|r µ ≡Ŝ + rµ |0 }, where |0 ≡ | ↓↓ . . . ↓ , and S + rµ = (X rµ + iY rµ ) is the spinraising operator at site r µ . Restricted to this subspace, P F ≡ r µ |P F |r µ is a |Λ|-dimensional matrix, and can be efficiently analyzed. Intuitively, one can compute P F by simply 'hopping' the spin-flip following the SWAP gates in the Floquet cycle. With PBC, one sees that any spinflip circles an adjacent plaquette and returns to its starting position after the fourth step (Fig. 1a). This implies P F = 1 |Λ| , and henceÛ F =P F =1.
Although the Floquet operator is apparently trivial with PBC, it could still be topologically nontrivial and display protected, anomalous edge dynamics when the system is under open boundary conditions (OBC). Explicitly, we consider a cylindrical geometry periodic in x but open in y. This opens up two circular edges respectively consisting of the A sites at y = 1 and B sites at y = N y . The computation of the Floquet operator with OBC,Û F ≡P F , proceeds as before. In the bulk, i.e. for y = 1 and N y , one simply findsP F | Bulk =1, as none of the SWAP gates acting on the bulk sites have been affected. On the other hand, for a spin flip starting at, say, (x, y = 1) A , the SWAP gates for the third and fourth time step have been 'deleted', so in the single-spin-flip sector and henceP F permutes the sites (x, 1) A in the same way as the unit right-translation operator along the boundary, t y=1;A . Similarly, for a spin flip starting at (x, N y ) B , we have which corresponds to the action of the unit lefttranslation operator (t y=Ny;B ) −1 . Altogether, one findŝ where the bulk and edge DOF are explicitly decoupled. Hence, despite the apparently trivial bulk Floquet operator, the dynamics of the system is nontrivial at the edge, and is governed by the 1d translation operator.

B. General classes of bosonic chiral Floquet models
The only place we used the spin-1/2 nature of the model in the previous analysis was Eq. (3), which can be easily generalized to the case of spin (p − 1)/2 sites for any p ≥ 1. In other words, we can replace the two-dimensional site DOF with p-dimensional ones and obtain an infinite class of bosonic chiral Floquet models. At the lower edge (y = 1), each of these models gives the edge unitary operatorŶ =t (p) , the unit righttranslation operator for a chain with p dimensional site Hilbert spaces. Note, however, that spin-0 DOF (p = 1) are special, as in that case the full many-body Hilbert space is one-dimensional and hence any model is necessarily trivial. (For simplicity, we will drop the superscript (p) when the Hilbert space dimension is not emphasized, witht understood to bet (p) for some p > 1.) A further generalization is obtained by taking any of the above models and performing the four time steps in reverse. This results in a mirror image model, with the edge unitary operatorŶ now given by the reverse translations. Even more generally, one can consider a 'stacked' model with two independent layers: a model with sites of dimension p stacked on top of the mirror image of a model with sites of dimension q. This givesŶ =t (p) ⊗t (q) , where we lett (q) ≡ (t (q) ) −1 for clarity.
After the construction of this infinite set of models, a natural question is to ask if these models are all distinct, i.e. given positive integers p, q, p , q , is the (p, q) model equivalent to the (p , q ) model? At this stage, however, it is not even clear what the word 'equivalent' means, since the usual notion of bulk-boundary correspondence in equilibrium systems cannot be applied to our Floquet setup. To even pose the question, one must first develop a corresponding notion in Floquet systems.
We will undertake this task in Sec. III, where we show that the bulk and boundary dynamics can be systematically decoupled in any MBL Floquet system. Granted such decoupling, one can classify MBL Floquet systems by studying smooth deformations among their boundaries. Before proceeding with this analysis, however, it is instructive to first examine the classification in a heuristic fashion, using the fact that the (p, q) models constructed here showcase explicit bulk-edge decoupling. Specifically, in Sec. IV, we will see that in general the (p, q) model is equivalent to the (p , q ) model whenever p/q = p /q . Before that, we will present below a heuristic argument suggesting that the (p, 1) models are nontrivial for p > 1, while the (p, p) models are always trivial.

C. Anomaly of the translation operator
Our intuition with the bulk-boundary correspondence dictates that we should identify obstructions to 'smoothly' deforming the boundary operatorŶ to the trivial oneŶ =1, i.e. obstructions in interpretingŶ as the finite-time evolution of a local Hamiltonian defined only near the boundary. We sayŶ is 'locally genreated' when no obstruction is present, and is 'anomalous' whenever such interpretation is impossible. Such an anomalous boundary signifies, at a heuristic level, the presence of a topologically nontrivial bulk. As the (p, 1) models giveŶ =t (p) , we will argue they are nontrivial by showingt (p) is anomalous when p > 1. The actual proof of the claim is deferred to section IV, where we define the quantized 'chiral unitary index' proposed in GNVW and discuss its implications.
We will look for an obstruction to interpretingt (p) as a locally generated unitary. For our purpose, the class of 'locally generated' unitaries is basically the same as the class of finite-depth quantum circuits of local unitaries (FDLUs), [31] so one should look for an obstruction to designing an FDLU that implementst in a 1d system. [32] Sincet (p) is a site permutation operation, it can be built using only local SWAP gates, as we show in Fig. 1b. However, in this construction L − 1 layers of SWAPs are used on a ring of size L, and therefore the circuit depth diverges in the thermodynamic limit. This suggests the impossibility of simulatingt (p) using any FDLU, which, when proven, will lead to our claim that the (p, 1) models are topologically nontrivial.
For the (p, p) models, however, the edge operator Y =t (p) ⊗t (p) can be realized in a purely 1d system using a nearest-neighbor SWAP circuit of depth two (Fig. 1c). This construction immediately shows that the (p, p) models are topologically trivial for all values of p.

III. BULK-BOUNDARY CORRESPONDENCE IN MBL FLOQUET SYSTEMS
We have seen that the bosonic chiral Floquet models, featuring the 1d translation operators at the edge, are suggestively nontrivial. We will now provide a strong theoretical grounding for this claim. The first step is to establish the notion of bulk-boundary correspondence in Floquet systems.
The bulk-boundary correspondence is frequently applied in the classification of gapped equilibrium phases of matter, where the ability to cleanly distinguish between bulk and boundary DOF in the presence of a gap allows one to classify distinct bulk phases in terms of the anomalous properties of their edges, which are symmetry-protected to be gapless in d = 2. In a Floquet system, however, the notion of a ground state, and hence an excitation gap, is meaningless. Nonetheless, if a Floquet system is MBL, so that one can identify mutually commuting quasi-local DOF, [9][10][11] a version of the bulk-boundary correspondence becomes possible. While similar approaches have been adopted for 1d systems in Refs. [19][20][21][22][23], in the following we generalize the previous discussions and present a self-contained construction applicable to Floquet systems in any spatial dimension d, and then specialize to d = 2 and apply it to the bosonic chiral Floquet models we constructed.

A. MBL Floquet systems in any dimension
We begin by first clarifying what we mean by 'MBL' and 'distinct' in our Floquet context. We will call a Floquet system MBL if its Floquet operatorÛ F can be written in the bulk as a product of mutually commuting quasi-local unitary operators. Here, we say an operator is 'quasi-local' if its nontrivial action is exponentially localized. (See precise definition of 'MBL' in Eq. (10) below). Two MBL Floquet systems are said to be distinct if it is impossible to interpolate between their Hamilto-niansĤ 0 (t) andĤ 1 (t) within the class of MBL Floquet systems. More precisely, we sayĤ 0 (t) andĤ 1 (t) are distinct if for any continuous interpolation the corresponding family of Floquet operatorsÛ F (s) will fail to be MBL for some s.
Next we proceed to explain the asserted bulk-boundary correspondence in detail. For concreteness, we work on a large boundary-less 2d geometry, say a torus. Write the Hamiltonian asĤ where each local termĤ r (t) =Ĥ r (t+T ) is nontrivial only within a finite neighborhood b(r) surrounding lattice site r. Consider the Floquet operator which is uniquely-defined up to an arbitrary choice in the reference time t = 0. We sayÛ F is MBL when it can be written as a product of mutually commuting quasi-local unitariesÛ r :Û Here eachÛ r is a quasi-local operator approximately localized within a finite neighborhood b (r) of lattice site r, and [Û r ,Û r ] = 0. In the 'l-bit' notation commonly used to describe MBL systems, eachÛ r can be represented asÛ r = exp(−iF r ({τ r })), where {τ r } is a complete set of Hermitian, mutually-commuting, quasi-local and conserved operators (l-bits), and F r is a local real function.
Identifying the exact form of the l-bits for a given system is a nontrivial task, but in the following discussion we will not need to use this l-bit representation explicitly. Before we move on to develop the notion of bulkboundary correspondence in an MBL system, we pause to discuss the distinction between our use of the phrases 'MBL Floquet systems' and 'MBL Floquet phases'. To define a phase, one should demonstrate that the defining properties are robust, i.e. they are stable under all small physical perturbations. Assuming many-body localization is robust, one can view MBL systems as the outof-equilibrium analogues of gapped equilibrium systems, and define a phase as a collection of systems that can be smoothly deformed into each other while maintaining their MBL nature. Although we implicitly assume such robustness of 2d many-body localization throughout this work when discussing MBL Floquet phases, we expect our results to hold in some long pre-thermal regime independent of the robustness of many-body localization. [33][34][35][36] B. Bulk-boundary correspondence To discuss bulk-boundary correspondence in an MBL Floquet system, it is useful to first define two length scales that are important in our MBL Floquet problem: the 'localization length' ξ and the 'Lieb-Robinson length' LR . Physically, ξ is the maximum radius of b (r), the neighborhood of sites on which anyÛ r can act nontrivially. Equivalently,Û r can be viewed as a local operator acting in b (r) tensored with the identity on the sites outside b (r). Due to the mentioned exponentially decaying tails of the l-bits, ξ so-defined is in principle unbounded. Instead, we will relax our definition and allowÛ r to differ from the identity outside b (r) by a small error. This allows for a finite ξ, with the corresponding error bounded by O(e −r /ξ ) for sites at a distance r away from b (r).
Another important length scale is LR ≡ T v LR , where v LR is the maximum Lieb-Robinson velocity [37] ofĤ(t) for 0 ≤ t ≤ T . The intuition behind this 'Lieb-Robinson length' LR is that since the Floquet period T is finite andĤ(t) is local, there should be a finite causal 'lightcone' beyond which no communication, i.e. information exchange, is possible. Using the Lieb-Robinson theorem, this intuition is formalized as follows: for any localized operatorÔ,Û † FÔÛ F can act nontrivially only on sites within a distance LR from the support ofÔ, again up to exponentially small corrections. Note that, for an MBL system, ξ and LR are not independent: From Eq. (10), one sees that LR ≤ 2ξ, the maximum possible distance between any two points lying in the same neighborhood b (r). Now consider a large region D in our lattice (a disk in 2d, or a ball in d > 2), with diameter much bigger than 2ξ. There are two unitary operators that can be naturally associated to the system on the ball D. One is the Floquet unitary constructed fromĤ D , the Hamiltonian truncated to D: The other is just the truncation of the original Floquet unitaryÛ F of the full system to the ball D: Note that while U F can be constructed for any Floquet system, to construct U F one needs to invoke the MBL assumption. By the Lieb-Robinson bound, these two uni-tariesÛ F andÛ F must agree in the interior of D, with an error that decreases exponentially for sites away from the boundary. In other words, their mismatch, is a unitary operator approximately localized on sites within a distance 2ξ from the boundary of D (Fig. 2). By construction, 2ξ is much smaller than the linear size of D, and thereforeŶ is effectively a unitary operator acting on the (d − 1)-dimensional boundary region of D. Note thatÛ F is defined using the time-dependent Hamil-tonianĤ D (t), and henceŶ can capture information that is missing in the bulk Floquet HamiltonianĤ F . In particular, this includes the physical consequences of a possibly nontrivial 'micro-motion' within the Floquet period.
The key observation now is that even thoughÛ F is MBL everywhere by construction,Û F may not be MBL at the boundary, whereÛ F andÛ F do not necessarily agree. Hence,Ŷ is not MBL in general. Rather, we can only argue thatŶ is locality-preserving, in the sense that conjugating byŶ takes local operators to nearby local operators (as follows from the Lieb-Robinson bound). Crucially, there is no guarantee that a locality-preserving unitary can be interpreted as the finite-time evolution operator of a local Hamiltonian. Specifically, in the case of Y , there is no guarantee thatŶ can be interpreted as a finite-time local Hamiltonian evolution in d − 1 dimensions. Whenever there is an obstruction in such reinterpretation, we sayŶ is anomalous and this signifies the presence of a topologically nontrivial bulk. We emphasize here that the mentioned anomaly makes no reference to symmetries, and is thus a 'stronger' anomaly than those associated with the boundary of MBL-SPT Floquet phases recently studied in Refs. [19][20][21][22][23].

C. Application to the bosonic chiral Floquet models
We will now apply the formal construction above to the bosonic chiral Floquet models we described in Sec. II, corresponding to d = 2. For simplicity, we will specialize the discussion to the p = 2 model presented in Sec. II A, although it generalizes immediately to the general classes of (p, q) models discussed in Sec. II B.
Recall that, with PBC, the Floquet operator of the model is simplyÛ F =P F =1. Now imagine appending to the original Floquet evolution a fifth time step, corresponding to the time evolution operatorÛ (5) . The bulk Floquet operator is then modified intoÛ F =Û (5)PF = U (5) , and in particular the system is MBL ifÛ (5) is MBL. In fact, even without this fifth time step, i.e. even if we had setÛ (5) =1, we would have also obtained a localized model, since, as noted above, the time evolution operator for the first four time steps is simplyP F =1. The purpose of the disorder introduced in the fifth time step is to render the MBL property robust against small arbitrary perturbations of the Hamiltonian, and hence allow our system to represent a stable Floquet MBL phase. The particular form of the disorder is unimportant, as long as the resulting unitary can be deformed into the model we constructed while remaining MBL, i.e. we stay within the same MBL Floquet phase. In the following, we will choose a particular form ofÛ (5) for concreteness and convenience.
First we replace time intervals T /4 with T /5 in the clean system and simultaneously increase the strength of the Hamiltonian such that the single-step evolution operator,Û (s) for s = 1, . . . , 4, is unchanged (Eq. 3), and then consider the disordering Hamiltonian where the coupling constants J x ∈ [0, 5π/T ] are independent random variables. Again, H (5) is exactly soluble, witĥ Note that with this choice ofÛ (5) , the model also has a U (1) symmetry generated by spin rotation about the z-axis. This symmetry is inessential to our construction. Once again, one can consider a more general disordering Hamiltonian breaking this symmetry, and as we will argue in Sec. IV the nontrivial nature of the phase survives.
To compute the edge operatorŶ as defined in Eq. (13), we have to consider two Floquet operators defined on an open geometry: the operatorÛ F corresponding to the evolution of the truncated Hamiltonian, andÛ F , the direct truncation of the bulk Floquet operator. Although we previously chose the edge to be the boundary of a disk D, in the present construction it is more convenient to take the cylindrical geometry adopted in Sec. II, which corresponds to taking OBC in one of the two torus directions.
To evaluateÛ F , recall from Eq. (6) that, with this OBC, the first four time steps in the Floquet cycle give rise to the evolution operatorP F =t y=1;A ⊗ (t y=Ny;B ) −1 . Since the disordering Hamiltonian Eq. (14) is on-site, it is unaffected by the change in boundary condition, and therefore we find By the same token,Û F =Û (5)P F , whereP F is the truncation of the PBC permutation operatorP F to the cylinder defined by the OBC, i.e. we should discard all terms in P F that cross the boundary between y = 1 and y = N y . However, asP F =1 it is again unaffected by the truncation. This givesÛ F =Û (5) , and thereforê giving the same edge operator as in the clean model in Sec. II. Note that the simplicity ofŶ hinges on our choice of a simple on-site disordering HamiltonianĤ (5) . For a more general disordering Hamiltonian,Û (5) andÛ (5) will only cancel exactly in the bulk, and modifyŶ by a unitary operator bi-local on the two circular edges. Nevertheless, the tensor-product factorization in Eq. (17) will still hold with an accuracy O(e −Ny/ξ ). We will see in Sec. IV that, as long as we stay in the same MBL Floquet phase, such modification does not affect the nontrivial nature of the model, which is captured by the chiral unitary index of a single edge.
To further support the claim that the MBL chiral Floquet model, characterized by the edge operatorŶ =t at a single edge, is in a nontrivial phase, it is instructive to study its behavior as we tune the system to a manifestly trivial fixed point. More explicitly, consider two 2d MBL Floquet systems defined on the same cylinder geometry, described by the Floquet operatorsÛ F (0) and U F (1) with respectively anomalous and trivial edges, saŷ Y (0) =t andŶ (1) =1. Imagine a smooth family of Floquet systems {Û F (s) : 0 ≤ s ≤ 1} interpolating between the two. By our claim,Û F (s) will necessarily become delocalized at some s c if the MBL chiral Floquet model is indeed nontrivial. To illustrate this, we define a family of models by randomly deleting a fraction s of the SWAP gates applied in the Floquet period. WhileÛ F (s = 0) is simply the original chiral Floquet model, very few sites are permuted as s → 1, and henceÛ (s = 1) =Û (5) and Y (1) =1. As shown in Fig. 3, the system indeed fails to be MBL around s = 0.5, consistent with the argument above. Plotted against the right axis is the localization length ξ. ξ was computed for N × N systems with PBC, where N = 60, 80, 100 and 120. As shown in the plot, the system is localized for small and large values of s, but is delocalized around s = 0.5. Accompanying the delocalization is a change in the edge behavior, as can be seen in the data on W, plotted against the left axis. W is defined as follows: The edge operatorŶ for each disorder realization is computed following the discussion near Eq. (17), using a 100 × 100 system in the cylindrical geometry periodic in x but open in y. Regardless of s,Ŷ can be interpreted as an element of the site permutation group, and hence can be factorized into mutually commuting terms, with each term being an oriented 'loop' in the lattice (Appendix A). The edge operator at a single edge is identified by restricting to the 'loops' that are strictly contained within the lower half of the system, and W is defined as the net winding number of these 'loops' about the cylinder axis. W = 1 and 0 respectively indicate a chiral and non-chiral edge operator. Note that W is not welldefined near s = 0.5, where ξ and N become comparable, and therefore W is not quantized in that region.

IV. THE CHIRAL UNITARY INDEX
In the previous sections we reduced the problem of classifying 2d bosonic MBL Floquet systems to that of classifying 1d locality-preserving unitary operatorsŶ , and we constructed an infinite set of models, labeled by two positive integers (p, q), which realizesŶ =t (p) ⊗t (q) . Having provided circumstantial evidence that some of these models, e.g. the (p, 1) models with p > 1, are nontrivial, we now argue that the quantized index defined in GNVW classifies bosonic MBL chiral Floquet phases in 2d via bulk-edge correspondence. We will call this index the 'chiral unitary index ' and denote it by ν. Roughly speaking, ν(Ŷ ) is a measure of the imbalance in quantum information pumping at the edge under the evolution governed by the locality-preserving unitary operator Y . As we will see below, ν takes the form log(p/q), where p, q are relatively prime positive integers.
In particular, in this section we will justify the use of the adjective 'chiral', show that ν is insensitive to the shape of the edge, and argue that it is in fact robust to any smooth deformation preserving the MBL nature of the bulk.

A. One-dimensional quantum information transport
The goal of this section is to develop a classification of 1d locality-preserving unitary operators, which we denote byŶ . Recall that by 'locality-preserving', we mean that conjugation byŶ takes local operators to nearby local operators. Above we have argued heuristically that certain locality-preserving unitaries, liket, are anomalous, i.e. cannot arise as the finite-time evolution of a 1d local Hamiltonian. Precisely such a class of operators was studied in GNVW, which showed that the anomaly associated withŶ is exhaustively captured by a quantized index, with the translation operators playing the role of generators for nontrivial indices.
To get some intuition on what makes a localitypreserving operatorŶ anomalous, it is useful to first recall the equilibrium classification of 2d gapped quantum phases. In the equilibrium setting, the effective low energy dynamics of the 1d edge is typically described by a conformal field theory, which is partially characterized by its 'chiral central charge'. The chiral central charge describes the chiral flow of energy along the edge at finite temperature -roughly, it counts the number of chiral edge modes. When non-zero, it signals an anomaly, i.e. an obstruction to realizing the corresponding field theory as the low energy theory of a truly 1d microscopic system.
In the present Floquet setting, instead of a 1d effective low energy field theory, we are given a locality-preserving unitaryŶ governing the edge dynamics. The analogue of the anomaly in this setting, namely the inability to real-izeŶ as the Floquet operator of a truly 1d system, turns out to also have an interpretation in terms of chiral transport. While the appearance of chirality may have been expected given the vital role played by the translation operator in the previous sections, the transport interpretation is now more subtle, as there is no conserved charge or even energy.
To understand what is being transported, consider again the unit translation operatort on a 1d ring of length L 1, and letρ be the density operator of any quantum state defined on the ring. For any local mea-surementÔ x , the transformed stateρ ≡tρt † gives results obeying Tr(ρ Ô x+1 ) = Tr(ρÔ x ), where Tr denotes the trace over a set of orthonormal bases in the physical Hilbert space. Stated in words, the measurement results ofρ at x + 1 are identical to those ofρ at x. Such unidirectional, perfect correlation is independent of the choice and structure ofρ, and is indeed a property of the operatort itself. Correlation is information -what is pumped in a state-independent manner is therefore, suggestively, quantum information, which are present in our Floquet setting even when both charge and energy are not conserved.
To formalize this idea, we now give a precise definition of the chiral unitary index of ν(Ŷ ).

B. Definition of the chiral unitary index ν
Let us work with a finite, but arbitrarily large, 1d ring whose lattice sites, labeled by x, host bosonic 'spin' Hilbert spaces H x of dimension p x . For concreteness, pick an orthonormal basis {|i x } for each H x , and let A x be the set of local operators at a site x: Since the elements of A x , being operators, can both be added and multiplied, A x forms an algebra. An explicit orthonormal basis for Additionally, one can consider more general operator algebras A S consisting of operators acting on a collection of sites S: Now letŶ be a locality-preserving unitary operator in this 1d system. As described in the previous section, the degree to whichŶ preserves locality is quantified by a Lieb-Robinson length LR . Recall the defining property of LR is that for any local operatorÔ,Ŷ †ÔŶ can act nontrivially only on sites at most a distance LR from the support ofÔ, up to exponentially small corrections.
We now set up a 'flow gauge' that measures how the locality-preserving unitaryŶ pumps quantum information across a spatial cut in the system. Specifically, consider two contiguous intervals of sites L and R, residing immediately to the left and right of the cut respectively. These intervals host independent operator algebras A L and A R , in the sense that [â L ,â R ] = 0 for allâ L ∈ A L andâ R ∈ A R . For brevity, we will say A and B commute and write [A, B] = 0 when they are independent in this sense.

captures the extent to which
A L becomes entangled with A R under the action ofŶ , and therefore reflects the amount of quantum information pumped from left to right across the cut. [38] Similarly, a measure of how [Y (A R ), A L ] = 0 captures the flow of quantum information from right to left. To fully characterize the chiral character ofŶ , we should look at the difference between these two measurements.
To ground the discussion quantitatively, we will need to introduce a measure of the extent to which two sets of operators A and B, acting on a common Hilbert space H Λ with dim(H Λ ) = p Λ , are independent, i.e. we seek a measure that is minimized when the two sets of operators commute ([A, B] = 0), and is maximized when the two sets coincide (A = B). A natural candidate is the overlap of the respective orthonormal operator bases, which we define below. Let A and B be respectively spanned by the bases {ê a ij : i, j = 1, . . . , p a } and {ê b lm : l, m = 1, . . . , p b }, and define where the trace Tr Λ is restricted to H Λ . In Appendix C, we show that η so-defined enjoys the desired properties: (i) it is independent of the orthonormal basis chosen; With all these preparations we can now define the index. We let A L and A R be as above, but with the extra condition that they have sizes ≥ LR . Let H Λ be the Hilbert space associated with a very large finite region Λ containing both L and R (for example Λ could be the whole 1d system, assumed to be a finite ring). We then define the chiral unitary index ν(Ŷ ) as where ind(Ŷ ) is equivalent to the index defined in Eq. (45) of GNVW, and we will refer to it as the 'GNVW index'. Although there is no formal distinction between ν and ind, we nonetheless introduce ν, i.e. takes the logarithm of the GNVW index, such that the connection of the index with the more familiar equilibrium results will be more transparent, as we elaborate on below.

C. Index of the translation operator
To develop some intuition for the chiral unitary index ν(Ŷ ), or equivalently, its exponential, the GNVW index ind(Ŷ ), let us first compute it for two simple cases. First, ifŶ =1 is the identity operator, then clearly ν(1) = log(1) = 0, since [A L , A R ] = 0. Second, in the case that the spin Hilbert spaces H x all have the same dimension p x = p > 1, we can letŶ =t be the unit right-translation operator. The Lieb-Robinson length oft is just 1 and therefore intervals L and R of length 1 -i.e. single sites -are already sufficiently large to be used in the index computation in Eq. (20).
which is non-zero as long as p > 1, and so as claimed the translation operator has a nontrivial index. In addition, by a similar computation one sees that ν(t † ) = − log p, consistent with the intuition thatt andt † pump quantum information in opposite directions.
Pictorially, the GNVW index of the translation operator corresponds to the transport of the entire on-site p-dimensional Hilbert space H x across the cut. Note also that the computation of ν(t) above is independent of x, the location of the cut. In addition, one can check that expanding the sizes of L and R will leave ν(t) invariant. Hence, ν(t) is independent of the arbitrariness in defining A L and A R , as one would expect for a well-defined index.
Sincet 2 brings the p 2 -dimensional Hilbert space H x−1 ⊗ H x across the cut, we should expect ν(t 2 ) = log(p 2 ) = 2 log p. In addition, for a system stacked with an identical copy of itself,t ⊗t also brings a p 2dimensional Hilbert space across the cut, and similarly we expect ν(t ⊗t) = 2 log p. These observations suggest that, as claimed, the chiral unitary index is additive (or equivalently the GNVW index is multiplicative) under both composition (multiplication) and tensor product of the unitaries [26] (Appendix C):

D. Interpretation of the index
We have motivated the definition of the chiral unitary index, Eq. (20), by quantifying how a localitypreserving unitary operator redistributes quantum information. This is exemplified by the index computation for the translation operator, which gives ν(t (p) ) = log p, the maximum von Neumann entropy supported by a plevel system. These facts make it natural to identify the chiral unitary index with the flow of information entropy per Floquet cycle along an 1d edge of the 2d system.
Such unidirectional transport of physical quantities is a defining feature of chiral phases of matter. For instance, equilibrium chiral phases are accompanied by a quantized nontrivial thermal Hall conductance. Can we then view the MBL chiral Floquet phases as being in one-to-one correspondence to their equilibrium counterpart, with the only difference being a transport of quantum information in lieu of energy quanta? We will attack this problem by first posing a seemingly different one: Conventional wisdom holds that chiral phases with counter-propagating edge modes can 'cancel' each other and reproduce a trivial phase. How does this intuition generalize to the bosonic MBL chiral Floquet phases?
To answer this, it is instructive to first consider the (p, p) model, which should be trivial by our intuition on 'edge cancellation'. To see this, we simply note that it has a trivial index: ν t (p) ⊗t (p) = ν(t (p) )+ν(t (p) ) = log p− log p = 0. Physically, this implies the operatort (p) ⊗t (p) is not anomalous, i.e. it can be well approximated by an FDLU, as we have shown in Sec. II C.
Importantly, not all translation operators are equal. For instance, consider the (2, 3) model, which has index ν(t (2) ⊗t (3) ) = log(2) − log 3 = 0. Physically, this amounts to the observation that the information capacity of a qutrit (p = 3) is larger than that of a qubit (p = 2). Hence, two chiral unitary operators annihilate each other if and only if they are 'equal and opposite'. The phenomena of 'opposite' operators annihilating each other is familiar from the equilibrium classification of SRE chiral phases, but the presence of logarithms in the current Floquet situation is novel.
Since the (p, q) model has index log(p/q), under stacking of these models their chiral unitary indices form (log Q + , +), the group of log-positive-rational number under addition. [26] This suggests that the classification in our MBL Floquet setting is much richer than that of their equilibrium counterparts. Nonetheless, if one is restricted to quantum spin systems with isomorphic p-dimensional site Hilbert spaces, the realizable indices are reduced to a subgroup of the general case. For instance, if p is prime then the classification is reduced to (log p Z , +) (Z, +), similar to the equilibrium result -this follows immediately from the alternative definition of the index given in section 7.2 of GNVW. Slightly more generally, if p has n p prime factors, then each prime factor generates one 'dimension' of a 'lattice' and the classification becomes (Z np , +).

E. Robustness of the index
While we have focused on the (p, q) model and their edges,Ŷ =t (p) ⊗t (q) , to develop a physical picture for the chiral unitary index ν, the narration thus far is missing a crucial ingredient: the quantization of ν. This quantization is vital for its role in the classification, since it is needed to ensure stability of the classification to any continuous change of the system, e.g. coming from a change in the precise shape of the 1d boundary chosen in the definition of Sec. III, or from generic perturbation to the time-periodic Hamiltonian. To see this explicitly, if this quantization was absent, one could only claim that the group (log Q + , +) classified the special class of (p, q) models, but could not view it as the general classification of 2d bosonic MBL chiral Floquet phases.
This quantization property is nontrivial, since a priori a local dressing of the translation operator could already interfere with the transport of quantum information. Specifically, although it is not at all obvious from its definition in Eq. (20), the GNVW index ind(Ŷ ) turns out to be a positive rational number for any locality-preservingŶ . In GNVW, a proof of this statement is based on an alternate but equivalent definition of the index as the ratio of the dimensions of two certain finitedimensional operator algebras, which is manifestly a positive rational number. Mathematically-minded readers are encouraged to consult section 7 of GNVW for further details of this argument. In the remainder of this subsection, however, we will instead provide a more physical argument for the quantization.
We first present the main idea behind the derivation.
Intuitively, small physical deformations can be viewed as a dressing of the original operator by an FDLU, which is built from local unitaries. However, a local unitary cannot be chiral: Suppose on the contrary that a local unitary operatorÛ I , defined on the finite 1d interval I with boundaries ∂ L and ∂ R , is chiral. Then under evolution governed byÛ I , quantum information is gradually depleted from (say) ∂ L and accumulates at ∂ R . This is at odds with unitarity, since the dimensions of the local Hilbert space at the two edges are fixed. Hence, we see that ν(Û I ) = 0, and from Eq. (22) we conclude that any FDLU has ν = 0 and therefore the index is robust. Note that the contradiction above is evaded by locality-preserving unitaries defined on a boundary-less geometry, say on a ring or on an infinite line. These are precisely the geometries for which the translation operatort is well-defined, and indeed the geometries that arise as the boundary of a 2d system. With this picture in mind, we proceed to show that Eq. (20) is invariant under smooth physical deformations. Specifically, consider two locality preserving uni-tariesŶ (0) andŶ (1), together with a smooth interpo-lationŶ (s) between them, parametrized by s ∈ [0, 1]. Assume that the Lieb-Robinson lengths LR (s) ofŶ (s) are all uniformly bounded by some LR : LR (s) ≤ LR . Thenγ(s) ≡Ŷ (s)Ŷ (0) † is a family of unitary operators contractible to the identity. For small δs, we can expand whereĥ(s) is a Hermitian operator. Sinceγ(s) is also locality-preserving,ĥ(s) has to be local. Thereforeγ (1) can be viewed as generated by a local Hamiltonian evolution, and is hence well approximated by an FDLU. Aŝ Y (1) =γ(1)Ŷ (0) and the chiral unitary index ν is additive under composition (Eq. (22)), it remains to show that ν = 0 for any FDLU. Using the composition property once again, one simply needs to show the triviality of a single layer of (necessarily commuting) local unitaries. In this situation, the only local unitary that can potentially lead to a nontrivial index is the one entangling A L and A R , i.e. the unitary that is sliced by the cut. To see that ν = 0 for such a unitary operator, we present below an argument originally given in GNVW. LetÛ LR be the local unitary sitting at the cut. We can take L and R as large as we please such thatÛ LR . . , p α } as a set of orthonormal basis for H α , α = L, R, we can writeÛ where (U LR ) Computing explicitly (assuming index summation con-vention), one finds in the last line of Eq. (25) amounts to an index relabeling. As all the indices are summed over, the expression is unchanged and hence ν(Û LR ) = 0.
We have thus shown that the chiral unitary index is invariant under small continuous deformations. As a corollary of the above arguments, we also see that the index is independent of the location of the cut. However, since the composition property, Eq. (22), is not explicitly derived in this work, our arguments do not constitute as a proof of the index quantization -for a rigorous proof of this fact, see section 7 of GNVW. Additionally, GNVW proves that if ν(Ŷ ) = 0,Ŷ must necessarily be an FDLU. Combined, these statements imply that any locality-preservingŶ can be continuously connected to a (p, q)-edge for some relatively prime positive integers p, q, and so the stacked (p, q) models constructed in the previous section form a complete set of representatives for bosonic MBL chiral Floquet phases in 2d.

V. NUMERICAL INDEX COMPUTATION VIA MPUS
We have argued that the quantized chiral unitary index implies the existence of 2d bosonic MBL chiral Floquet phases, and that the SWAP models we presented serve as prototypes. However, since the class of problems under consideration is strongly interacting, one may expect that a general method to numerically compute the index is unavailable -in fact even specifying the edge unitary operatorŶ in a concrete manner is technically challenging. Contrary to this expectation, we will demonstrate below that the index is in fact numerically computable through the use of matrix-product representations.
One of the most successful frameworks for representing a 1d quantum operator is through the use of a 'matrix product' representation, [39,40] in which an operatorÔ is written aŝ ix+2,jx+2 · · · × | · · · i x i x+1 i x+2 · · · · · · j x j x+1 j x+2 · · · |, where for each x and i x , j x = 1, . . . , p x , M ix,jx is a χ x × χ x+1 matrix, with the maximum value of {χ x } known as the bond dimension. The ellipses may extend to infinity for an infinite chain, terminate at two ends for an open chain, or be subjected to a trace in the bond space for a system with periodic boundary condition. In this work we will focus exclusively on matrixproduct unitaries (MPUs). In addition, in discussing MPUs it is often convenient to introduce a diagrammatic representation, where each tensor at a site is represented by a box with two 'physical' and two 'virtual' legs (Fig. 4a). The physical legs are attached to 'bras' and 'kets', corresponding respectively to in-coming and out-going physical states; the virtual legs are contracted with those of the neighboring sites, and provide 'room' for quantum information transport, which is necessary when the unitary is not on-site.
FDLUs in 1d are well-described by MPUs with finite bond dimensions, since each local gate in the circuit can be readily represented as an MPU. [39,40] However, to simulate the edge of a 2d MBL Floquet phase, we must relax ourselves from locally generated unitaries to locality-preserving unitaries. The prototypes of locality-preserving, but not locally generated, 1d unitaries are the translation operators. These turn out also to admit a simple MPU representation. To see this, assume dim(H x ) = p for all sites x, and let the bond dimension be also p. Now consider the tensor T [x] ix,jx α,β = δ ix,α δ jx,β , which gives an MPU with bond dimension χ = p. Pictorially, the delta functions act as 'connectors' of 'quantum information pipes', with the pairing of indices leading to a 'pipe' that connects the in-coming states at x to the out-going ones at x + 1 (Fig. 4b). The tensor T [x] therefore defines the unit right-translation operator, as one can explicitly verify.
Having discussed the MPU representations of both the FDLUs and the translation operators, we now make connection with the study of 2d MBL Floquet phases. Recall that the GNVW classification is exhaustive, i.e. any such unitary can be written as the product of an FDLU and a (stacked) translation operator with a suitable index. Since both the FDLUs and translation operators admit MPU representations, the GNVW result implies all locality-preserving 1d unitaries can be efficiently simulated by MPUs. Further using the bulk-boundary correspondence we established in Sec. III, one sees that the MPU representation provides a universal framework for describing the edge of any 2d MBL Floquet system.
In addition, by representing the edge operatorŶ of a 2d MBL Floquet system as an MPU, one can also compute the chiral unitary index ν of the system within the MPU formalism. To leverage the power of this framework, we recast the index formula Eq. (20) into the MPU language. Although a concrete derivation, which we present in Appendix D, will necessary involve a fair bit of technicality, its diagrammatic representation is quite intuitive, as we discuss below.
As discussed, the in-coming and out-going legs of an MPU can be loosely viewed as the inlets and outlets of a pipeline, with each of the openings labeled by their location on the chain. Since ν is a measure of the net quantum information flow, in the pipeline analogy the index computation amounts to quantifying the rate at which an incompressible fluid in the pipeline flows across a spatial cut. This rate can be extracted by suitably comparing the inflow and outflow across the cut. For instance, to measure the left-to-right flow, one can imagine closing all the openings of the pipe, except for the inlet at x and the outlet at x + 1 . When some test fluid is pumped into the pipe at x, it can either get trapped inside the pipe, or emerge out at x + 1. The portion that flows out from the outlet must have passed through the cut between sites x and x + 1, and therefore its volume reflects the capacity of the pipe across that cut. In the MPU language, the analogue of 'closing pipe openings' is to contract the free-hanging legs, i.e. to take traces. This is shown in Fig. 4c , which enters the index computation via Eqs. (19) and (20), takes a form similar to the pipeline picture presented above.t In the analogy above, the inlet at x − 1 and the outlet at x + 1 could also be connected, and our measurement was blind to the flow through this connection, which also passes through the specified cut. Physically, the range of such 'connection' is the range of quantum information redistribution, which is what we defined as the Lieb-Robinson length LR . To accurately evaluate the flow across the cut, one should therefore increase the number of inlets and outlets open on the two sides until all possible connections have been exposed, i.e. one needs to choose the intervals L and R in the index computation to be at least as big as the Lieb-Robinson length. This can be observed in Fig. 5, where we show the numerical results on the computation of the chiral unitary index using the MPU formula in Appendix D.

VI. EXPERIMENTAL PROPOSAL
Having constructed explicit models and identified a topological index for bosonic MBL chiral Floquet phases in 2d, we now propose an experimental setup that realizes the ν = log 2 phase using hardcore bosonic ultracold atoms loaded onto a shaken optical lattice.
In the hardcore limit, each site Hilbert space is twodimensional and can be viewed as an effective spin-1/2. Inspired by previous experiments in creating more complicated lattice geometries, [41,42] we will consider a pair of short-and long-wavelength optical lattices. Two square lattices (red and blue), formed by two pairs of retro-reflected laser beams, are used, and a deep vertical lattice is employed to render the system quasi-2d (Fig. 6a). To create the checkerboard geometry, we take the ratio of the wavelengths to be √ 2, and rotate the two lattices by 45 degrees. We suppose the frequencies of the laser beams are chosen to be respectively blue-and reddetuned from the atomic resonance, and therefore the potential minima of the blue and the maxima of the red are energetically favored.
Suppose the blue lattice is deep such that the system is strongly Mott-insulating when the red lattice is ignored. The role of the shaken red lattice is to periodically 'turn on' the targeted bonds, which enhances tunneling of particles. Intuitively, the tunneling will be perfect when the shaking profile is appropriately designed, akin to the action of a SWAP gate.
More specifically, we consider a five-step protocol similar to our construction of bosonic chiral Floquet models in Sec. II A and Sec. III C. In particular, the last step was again introduced only for incorporating disorder into the system. This can be achieved by turning on a speckle potential, which creates a disordered on-site chemical potential. [43,44] Similarly, the first four steps will resemble the SWAP models. In these steps, each of duration T /5, the red lattice is ramped on and then off, with the phase of the lasers set to align the potential maxima (energy minima for the atoms) to the locations indicated in Fig. 6b. The potential maxima of the red lattice mark the 'on' bonds, while all the other bonds are considered as 'off'. In a crude approximation, we neglect the tunneling across the 'off' bonds, and in each step s = 1, . . . , 4 we consider the Bose-Hubbard Hamiltonian H s (t) ≈ r∈RL sĤr (t), wherê Proposed optical lattice setup for realizing the bosonic MBL chiral Floquet phase using ultracold hardcore bosons. (a) We consider two square lattices, indicated respectively by blue circles and red diamonds, that are rotated by 45 degrees with respect to each other. We assume the blue lattice is deep and the system is strongly Mott insulating without the red lattice. The role of the red lattice is to offset the energy barrier between pairs of sites, and thereby turning 'on' the bonds between them, indicated by purple lines. (b) The red lattice is both amplitude-and phase-modulated, and sequentially turns on different sets of bonds in the first four steps of the Floquet cycle. modulated as the lattice is ramped on and off, though that of U (t) is relatively unimportant as we will take the hardcore limit anyway.
By construction, all the terms inĤ s (t) commute, and we simply solve the dynamics governed by the two-site systemĤr(t). In the hardcore limit U (t) → ∞, the effective Hamiltonian in the low-energy subspace, with the particle-number basis (|00 , |10 , |01 , |11 ), is a 4 × 4 matrix for which the time evolution operator can be easily com-puted: where φ(t) = t 0 dt J(t ) is the integrated rotation angle in the one-particle subspace. While a more careful estimation will require the computation of J(t) for realistic lattice parameters, it should be possible to design an appropriate amplitude-modulation profile realizing φ(T /5) = π/2, which gives a 'perfect' tunneling. As such, U ḡ r (T /5) is equivalent to the SWAP gate up to the phase rotation R ϕ ≡ diag(1, −i, −i, 1).
Thanks to the topological nature of the problem, this phase rotation does not affect the realization of the model. To see this explicitly, note that the phase rotation above can be written aŝ with the number operatorn µ = |1 µ 1 µ | in the hardcore limit.R ϕ can be interpreted as the evolution of a Hamiltonian diagonal inn µ . Since all the number operators commute,Û 2d F remains exactly soluble. In particular, the Floquet operator (together with the fifth, disordering step) takes the MBL form, with {n rµ } being the l-bits and the functions F j only coupling l-bits in the same or adjacent unit cells.
The argument above can be readily extended to a smooth family of phase rotations connectingR ϕ to the identity, and hence the proposed system is indeed in the same ν = log 2 MBL chiral Floquet phase as the SWAP model in Sec. II A. By the same token, restoring the tunneling between the neglected bonds should have a mild effect provided that the system remains MBL, which is anticipated when the blue lattice is sufficiently deep.
While lattice shaking techniques have already been demonstrated in 2d, [45] simulating MBL models in ultracold atom experiments is still an active research frontier. Nonetheless, several recent works have already demonstrated their feasibility, [46][47][48] and realization of our proposal could be soon within reach.

VII. PHYSICAL CONSEQUENCES
In this section, we discuss the physical consequences and experimental signatures associated with the bosonic MBL chiral Floquet phases. We will focus on the thermal behavior of the edge in Sec. VII A, and the transport of quantum information in Sec. VII B.

A. Topology-enforced edge thermalization
To this end, it is instructive to first recall the phenomenology associated with SRE topological phases of matter in equilibrium. In 2d, a bulk with SPT order is always paired with an anomalous edge, which cannot be both symmetric and gapped, i.e. the edge is either gapless, or breaks a symmetry.
In our non-equilibrium setting, the role of an excitation gap is now played by many-body localization, and therefore the edge of a nontrivial 2d MBL Floquet SPT phase cannot be both symmetric and MBL, i.e. if symmetries are not broken, the edge must be thermal. [49] When the protecting symmetry is broken, however, the boundaries can be localized again. In particular, some of the proposed phases feature boundaries that become localizable when the discrete time-translation group is reduced to a subgroup, say when sub-harmonic terms are added to the original drive. Interestingly, such symmetry breaking might happen spontaneously and lead to discrete time crystals [19,[50][51][52][53] (in Ref. [19] the analogous phase was termed π spin glass).
In fact, all the previously-proposed topologically nontrivial MBL Floquet phases are protected by some symmetries, and therefore the boundaries of these systems can always be localized by symmetry breaking. The bosonic MBL chiral Floquet phases we discovered, however, correspond to the first established examples that are nontrivial even in the absence of any symmetries, and hence the non-localizability is robust, provided that the system remains Floquet. [54] In particular, the anomalous edge is stable against the introduction of sub-harmonic terms. Such robust non-localizability of the chiral edge can be summarized as follows: If ν(Ŷ ) = 0, then for any 1d local Hamiltonian evolutionF (an FDLU), the locality-preserving operatorFŶ is never localized. This holds even when arbitrarily strong disorder is incorporated into the system viaF .
To substantiate this claim, we consider a 1d ring of spin-1/2's. LetF h ≡ e −iĤ h be the evolution operator of a disordered 1d Heisenberg Hamiltonian, given bŷ whereŜ i denotes the spin operators at site i, h is the strength of on-site random magnetic fields, and w i,α is a random number within the range [−1, 1]. It is known thatĤ h is in a thermal phase for h < 2.5, and an MBL phase for h > 2.5. [55] This is confirmed by looking at the statistics of the level spacing between adjacent energy eigenstates, δ n ≡ n+1 − n . As shown in Fig. 7a, the level statistics ofĤ h is Poissonian when h = 8 > 2.5, confirming that it is MBL. In addition, for this class of models the 'r ratio', defined by r = min(δn,δn+1) max(δn,δn+1) , is known to converge to 0.39 and 0.6 when the system is respectively MBL and thermal. [56,57] This can be seen from the data shown in Fig. 6b, in which the r ratio changes from 0.6 to 0.39 as h is increased, and is consistent with an MBL phase transition at h = 2.5. [55] In contrast, when combined with the unit righttranslation operatort, the operatorF ht is never localized Chiral edge thermalization. The presence or absence of localization is studied using exact diagonalization for a small spin-1/2 chain with size N ≤ 11 sites, and level statistics is extracted using 200 disorder realizations. Red circles and blue squares respectively indicate level statistics ofĤ h and i log e −iĤ ht . (a) The level spacing δn ≡ n+1 − n, where n is the n-th eigenvalue of the Hamiltonian for one disorder realization, is known to show distinct statistical properties when the system is MBL or thermal. The plot shows level statistics for h = 8, and the probability density functions of the normalized level spacing δ/ δ are well-fitted to the Poisson distribution (red line) and the generalized unitary ensemble (GUE; blue dashed line) respectively forĤ h and i log e −iĤ ht . This indicatesĤ h is localized but i log e −iĤ ht is thermal. (b) The localization of the system can be further quantified by studying the 'r ratio', defined as r ≡ min(δn,δ n+1 ) max(δn,δ n+1 ) . r is known to converge to 0.39 (red line) and 0.6 (blue dashed line) respectively for Poisson and GUE level statistics, [56,57] as observed at large h. The inset shows the convergence as a function of system sizes 8 ≤ L ≤ 11 at a strong disorder of h = 8. regardless of h. Such distinction between the behavior of F h andF ht can be readily seen in the level statistics, as we show in Fig. 7. In particular, althoughF h is MBL at a strong disorder strength of h = 8 > 2.5, the level statistics ofF ht agrees with that of a thermal phase (i.e., the r ratio remains at 0.6 in the entire range of h). [57] This suggests that the nontrivial bulk topology of a bosonic MBL chiral Floquet phase can lead to a robust chaotic behavior at the edge.

B. Unidirectional transport and quantum communication
Aside from connection to ergodicity, akin to its equilibrium counterpart a chiral edge operator also gives rise to signature unidirectional transport. For instance, in a system with a conserved U(1) charge, one expects an edge operatorŶ , with ν(Ŷ ) > 0, to generate a right-moving current proportional to the average charge localized at the edge. [30] However, we stress that the presence of such a conserved charge only bears witness to the chiral nature of the edge, and is not necessary for the existence of a chiral Floquet phase, which was linked to the chiral unitary index ν defined in the absence of any symmetries. For this reason, a more fundamental characterization of the chiral Floquet phases are in terms of the chiral transport of quantum information, which we investigate below. As a thought experiment, imagine Alice and Bob, sitting respectively at (0, 0) and (L, L), are separated by an experimental setup realizing the bosonic MBL chiral Floquet model discussed in Sec. II A. Suppose the system is initialized into a trivial product state, say with all spin-1/2's pointing 'up'. At time t = 0, Alice can locally create quantum entanglement by entangling a bulk and an edge qubit on her side, respectively at (1, 1) and (0, 0), and prepare the initial state where | ⇑ denotes the 'all-up' state for all sites other than those explicitly written out (Fig. 8a). After one Floquet period (t = T ), the bulk DOF only pick up an on-site random phase, whereas all the edge DOF are 'advanced' by one unit along the edge (Fig. 8b), and hence the state evolves into where e −iφ is an unimportant phase. Therefore, after 2L Floquet periods (Fig. 8c), one finds i.e. the 'pumping' of quantum information manifests as a transfer of local entanglement across the entire system. Hence, Alice and Bob now share an entangled pair of qubits, transported by the 'one-way information highway' at the edge of a 2d bosonic MBL chiral Floquet system, with which they can utilize for quantum communication.

VIII. FERMIONIC CHIRAL FLOQUET PHASES
We have so far focused on 2d bosonic systems, where we established rigorous results using the GNVW results and an MPU reformulation. However, the basic construction of the SWAP models is very closely related to the noninteracting fermionic AFAI models in Refs. [27][28][29][30]. In particular, the AFAI model features a chiral edge mode in the clean single-particle Floquet spectrum, and in fact we show in Appendix E that, viewed as a many-body operator, the AFAI Floquet evolution acts as the fermion unit translation operator along the edge. As our heuristic arguments apply equally well to fermions, this strongly suggests the AFAI models are also stable to interactions.
In light of these parallels, it is natural to ask how the bosonic and fermionic chiral Floquet phases are related: Is the fermionic phase equivalent to a bosonic one with a particular chiral topological index ν? In equilibrium, we know that chiral phases are characterized by their chiral central charge c, which measures the quantized thermal Hall conductance (aka gravitational anomaly) of the bulk. It is known that without intrinsic bulk topological order, the minimal chiral complex fermion phase is an IQH insulator with c C = 1, whereas the minimal chiral bosonic phase is the E8 state with c B = 8. [3] Hence, in equilibrium, one needs to combine 8 of the minimal complex fermion phases together to equate to a minimal bosonic phase. How is this relation modified in the Floquet setting, where continuous chiral flow of heat is replaced by discrete chiral pumping of quantum information?
Before diving into the details, let us anticipate the results. First, observe that the fermionic chiral Floquet phase (e.g. the AFAI model) pumps a two-state qubit of quantum information (0 or 1 fermion occupation) along the edge. This is the same quantum information capacity as the spin-1/2 bosonic model, which has index ν = log 2. Hence, we expect the 'conversion rate' between the minimal chiral Floquet phases of complex fermion and hardcore bosons to be 1 : 1. We will see that this is indeed the case, which is in stark contrast to the equilibrium result of c B : c C = 8 : 1. In the arguments we will assume that, even in the presence of fermions, the chiral Floquet phases are classified by an index that is additive under stacking. We will also encounter a Majorana fermion version of the model, which is loosely speaking the square-root of the minimal complex fermion chiral Floquet phase, and hence can be interpreted as having edges that pump a fractional amount of quantum information.
To establish this, we will first formally decompose complex fermions into pairs of real (Majorana) fermions, and then relate the Majorana chiral edges to complex fermions and subsequently hardcore bosons. Consider again a checkerboard lattice and a four-step driving protocol similar to the one described in Sec. II A, where instead of hardcore bosons, we let each site x host a complex fermion with creation operatorĉ † x . We then define the Majorana fermions viâ and consider a drive that acts nontrivially on χ but trivially onχ. The model construction and analysis, detailed in Appendix F, is essentially identical to before, with only minor modifications that do not affect the topological nature of the model. At a single edge, the clean model again featuresŶ =t M , the unit right-translation operator for the χ-Majorana fermions: where we let x indexes the site along the edge, andt M acts trivially onχ. We can also apply the same driving protocol to theχ fermions, which is effectively the same as stacking two copies of chiral Floquet models of Majorana fermions. Ifχ are driven in the opposite chirality compared to χ, at the edge we will have a pair of counter-propagating Majorana translation operators, which can arise from a purely 1d local Hamiltonian evolution and the bulk is therefore trivialized (Appendix F). The more interesting case is when χ andχ are driven with the same chirality. This gives the edge operatorŶ =t MtM , wheretM denotes the unit right-translation operator forχ. This gives and hence we identifyt MtM as the translation operator for the complex fermion, i.e. the conversion ratio between Majorana and complex fermions is 2 : 1, as one would expect. If we let η M , η C respectively count the copies of unit right-translation operators for Majorana and complex fermions (left-translations are counted as negatives), we can label any edge by (η M , η C ), which is additive under stacking. The two-to-one conversion between Majorana and complex fermions described above implies (2n + m, 0) ∼ (m, n) for any pair of integers n, m. Now we add hardcore bosons, or equivalently spin-1/2's, to the discussion. Their chiral unitary indices are always given by ν = η B log 2, where η B ∈ Z can again be viewed as a count of the number of bosonic chiral edge operators. The edge characterization is now extended to the tuple of integers (η M , η C , η B ). To connect fermionic and bosonic models, we will establish (3, 0, 0) ∼ (1, 0, 1), i.e. starting with an edge with three copies of chiral Majorana operators, we can convert a pair of them (corresponding to two-dimensional site Hilbert spaces) to a chiral hardcore-boson edge (Fig. 9). Combined with the conversion between Majorana and complex fermions, we have (1, 1, 0) ∼ (1, 0, 1) ⇒ (0, 1, 0) ∼ (0, 0, 1). This implies that, in our MBL Floquet setting, chiral edges of complex fermions and hardcore bosons can be converted to one another in a 1 : 1 ratio, which is in stark contrast to the equilibrium result of 8 : 1. Note, however, that one cannot directly convert a fermionic model to a bosonic one, since doing so will violate fermion-parity conservation. As we will see, in canceling (0, 1, 0) by stacking with (0, 0, −1), it is crucial that a trivial pair of counterpropagating Majorana chiral edges is present. Hence, the statement (0, 1, 0) ∼ (0, 0, 1) is to be interpreted in the sense of stable equivalence.
To complete this discussion, it remains to show (3, 0, 0) ∼ (1, 0, 1), which we will demonstrate using a relabeling argument. Let the three co-propagating Majorana translation operators be acting onχ i for i = 1, 2, 3, where the site label is suppressed. We can formally recast the site Hilbert space as that arising from a hardcore boson,τ , tensored with a new Majorana fermion,μ, by defining the following operators (on each site): where the operators {τ i } verify the algebra of Pauli matrices and conserve fermion parity. Importantly, [μ,τ i ] = 0, and henceμ andτ can be regarded as independent degrees of freedom. Restoring the site indices, the edge evolution simply translates:μ x →μ x+1 andτ x →τ x+1 . We can therefore reinterpret the edge as a chiral Majorana fermion together with a decoupled chiral hardcore boson, which implies (3, 0, 0) ∼ (1, 0, 1). Such equivalence implies the fermionic chiral Floquet phases are as stable as their bosonic counterparts. While we have not ruled out the logical possibility that the bosonic chiral phases are trivialized in the presence of idle fermions, it is highly unlikely on physical grounds. Nonetheless, in order to rigorously establish the stability of fermionic chiral phases, one must rule out this possibility by a nontrivial extension of the GNVW classification to fermionic algebras. We will leave this as a challenge for future works. Finally, we comment on the physical stability of the Majorana phase. While one can construct MBL Majorana models (Appendix F), experimentally relevant fermionic systems exhibit charge conservation, and cannot realize the Majorana phase without spontaneous superfluidity or superconductivity. Since such phases feature Goldstone modes (superfluids) or introduce long range interactions (superconductors), they appear to be at odds with many-body localization. Hence, it is more natural to consider the Majorana phases as appearing in fractionalized systems with topological order. [58] We leave these generalizations to future works.

IX. DISCUSSION
We showed that the edge dynamics of bosonic 2d systems in an MBL chiral Floquet phase can exhibit a chiral flow of quantum information. These chiral dynamics are anomalous in that they cannot be achieved by any local Hamiltonian evolution in 1d, and are characterized by the quantized chiral unitary index taking the form ν = log(p/q), where p/q ∈ Q + is a positive rational number. [26] These chiral Floquet phases rely on the robustness of MBL phases in 2d, so far an unproven conjecture. However, we reiterate that even in the absence of stable 2d many-body localization, our description will nonetheless provide an accurate universal characterization of parametrically long time pre-thermal dynamics in strongly disordered systems.
Aside from a formal classification via bulk-boundary correspondence, we have also constructed a full set of representative exactly soluble models, developed the concrete numerical framework for index computation, proposed an experimental realization using hardcore ultracold bosons in a shaken optical lattice, elaborated on the physical consequences arising from the nontriviality of the phase, and explored their fermionic counterparts.
Our results suggest that interaction, localization, chirality and periodic driving conspire to realize a novel phase of topological matter -a new tune that cannot be played if anyone in the quartet is missing. As chiral phases are known to be the 'root' of many topological phases, the present work stimulates numerous further inquires. We highlight three particular aspects below: (i) For bosonic systems, we were able to rigorously establish the stability and completeness of the chiral floquet classification via the formal machinery of GNVW.
For fermionic systems, we demonstrated the topological equivalence of chiral edges of certain fermionic and bosonic models. However, classification of fermionic (or more general anayonic) Floquet chiral phases, with the same rigor as the purely bosonic case, will require a nontrivial extension of the GNVW results to fermionic or general anyonic operator algebras, which we leave as an important open problem.
(ii) The classification of topological phases is generally richer in the presence of extra symmetry constraints, say global spin rotation symmetries. Recently, much progress has been made in the classification of Floquet SPT phases of interacting bosons, where once again MBL was invoked as an essential ingredient. [19][20][21][22][23] It was proposed that the analogue of cohomology classification in this context of Floquet SPTs in d spatial dimensions can be achieved by identifying the phases with elements of the cohomology group H d+1 (G × Z, U (1)), where G is a unitary internal symmetry group and Z denotes the discrete time translation. [21,22] Although similar cohomology classifications are known to be incapable of capturing the equilibrium chiral phases, [31] we point out in Appendix G an interesting potential link between the two in our Floquet setting. In addition, recall that the edge of an equilibrium SPT phase can be pictured as a pair of counterpropagating chiral edge modes symmetry-protected from cancellation. [31] What, if any, are the new phases of matter that arise from a symmetry-enriched version of the present work?
(iii) In Sec. VII B, we pointed out how the unidirectional transport of quantum information at the edge of a bosonic MBL chiral Floquet phase can be utilized for entanglement sharing, which is a basic ingredient for any quantum communication protocols. However, the idealized description there, which involves the exact addressing of the l-bits, cannot be directly applied to a system with generic disorder. Away from this idealized limit, there will almost certainly be a logarithmically slow spreading of entanglement into the MBL bulk (as for the dephasing in any MBL time evolution), and also to some extent back to the edge. This is in contrast to the ballistic spread in a thermalizing system, and partially quantifies the universal aspect of the 'robustness' to imperfect register addressing. Yet, how the thermal nature of the edge impacts the protocols is less clear. For instance, in the presence of disorder the chiral information channel is noisy, and would likely decohere the entangled pair. Can this decoherence be efficiently echoed away? Since the entanglement transport is a manifestation of the nontrivial chiral unitary index, can one leverage the quantized nature of the index to design a more robust quantum communication channel?
We leave these questions for future work.
Note added -After posting our manuscript, a related work by Harper and Roy appeared, [59] which gives arguments for the stability of the chiral edge of bosonic phases with particle-number conservation, and suggests a similar classification by rational numbers.
While properties [1] and [2] can be viewed as defining properties of the SWAP gate, [3] is a statement on locality and is more subtle. Indeed, since fermionic operators are never quite local, the discussion here is restricted to bosonic (spin) systems. Besides, a decomposition similar to Eq. (A2) is not unique. Yet, given any such bi-local decomposition, one can rewrite it in the stated symmetric form using properties [1] and [2].
We aim to establish that the group generated by the set of SWAPs, a subgroup of the group of unitary operators on H Λ , is equivalent to the symmetric group S |Λ| . Loosely, one simply notes that S |Λ| is the group of site permutations, and S |Λ| is known to be generated by pairwise exchanges, which are frequently called 'transpositions'. At the 'quantum level', these transpositions are given by SWAPs, and hence the claimed equivalence.
More concretely, we introduce the set of objects that are permuted in our context. For each site i ∈ Λ, we introduced the local observable algebra, A i , which is simply the set of local quantum operators, nontrivial only at site i, endowed with addition (over C) and multiplication. Now consider the set {A i : i ∈ Λ}, which for p > 1 is a finite set of |Λ| objects. We can therefore define the symmetric group S |Λ| as the (abstract) group permuting A i 's. Now, recall that S |Λ| is generated by (a subset of) the set of transpositions, and one can check from the listed properties thatτ i,j A iτi,j = A j andτ i,j A jτi,j = A i . Hence, we arrive at the stated correspondence.
This correspondence implies that any circuitĈ of SWAP gates corresponds to an element π(Ĉ) ∈ S |Λ| , where π is the isomorphism sendingτ i,j to the corresponding transposition in (the abstract group) S |Λ| . Two circuitsĈ andĈ are equivalent iff π(Ĉ) = π(Ĉ ). Hence, to solveĈ, one simply reduces the corresponding permutation π(Ĉ) to a simple form, as we will illustrate in the next subsection. In addition, note that the only condition we have imposed on the site Hilbert space is p > 1, and therefore the analysis of a SWAP circuit is largely independent of p.

Solutions of the SWAP circuits
As discussed in the main text, the restriction of a SWAP circuit to the single-spin-flip sector is only a |Λ|dimensional matrix, and hence can be solved efficiently. The mentioned correspondence between the group of SWAP circuits and the symmetric group implies the full many-body solution can be immediately inferred from this computation: One can show that the restriction to the single-spin-flip sector gives nothing other the 'natural permutation representation', which is a faithful representation of S |Λ| . Once we have determined π(Ĉ) using this representation, the entire spectrum and eigenstates ofĈ are determined, as we demonstrate below.
Recall that any element in S |Λ| can be written in a unique 'cycle decomposition', which follows from the natural permutation action of S |Λ| on the set of |Λ| objects. For instance, the notation (a, b, c, d)(e, f ) refers to a permutation a → b → c → d → a and e ↔ f . In addition, note that the cycles are 'disjoint', i.e. each cycle acts on a different subset of the set of symbols and therefore different disjoint cycles commute. A cycle involving only two objects, like (e, f ), is a 'swap' of the two and is called a 'transposition'. As discussed, S |Λ| is generated by the set of transpositions, and indeed one can verify that (a, b, c, d) = (a, b)(b, c)(c, d).
Using the established isomorphism, one can computê C explicitly via the cycle decomposition. For instance, suppose π(Ĉ) = (a, b, c, d)(e, f ), then Hence, the disjoint cycle decomposition of π(Ĉ) corresponds to a factorization ofĈ into commuting pieces, and any SWAP circuit on Λ can be viewed as a single layer of commuting unitary operators (which are not necessarily local). In addition, observe that π −1 ((a, b, c, d)) is nothing but the translation operator for the four-site 'ring' with sites a, b, c and d. As all the eigenvalues and eigenstates of the translation operator (for any number of sites and p) can be readily computed,Ĉ can be exactly solved.
Appendix B: The 'diluted' SWAP model and mapping to classical percolation In this appendix, we discuss in details the construction and analysis of the 'diluted' SWAP model, whose properties are shown in Fig. 3 of the main text.
Recall that we define the 'diluted' SWAP model as a bond-diluted version of the model introduced in Sec. II A, where any of the SWAP gates are removed at random with probability s, and left intact with probability 1 − s. Intuitively, the chiral phase should survive for a very small dilution probability s 1, whereas when most bonds are removed, s ≈ 1, one expects a trivial phase.
This defines a family of disordered model, labeled by 0 ≤ s ≤ 1, that connects the chiral bosonic model (s = 0) to a trivial model (s = 1). As detailed in Appendix A, a quantum circuit composing only of SWAP gates defined on a lattice Λ is exactly soluble, due to its equivalence to a permutation in the symmetric group S |Λ| . One can therefore also analyze the entire family of diluted SWAP models efficiently. Our main goal is to study the disorder average of ξ, the localization length of the Floquet operator (with PBC) defined in Sec. III of the main text. As the fifth time step in our model only gives rise to an on-site phase factor, ξ is determined by the permutation factorP F coming from the first four steps of the Floquet cycle. We have established in Appendix A that the permutation factor can be factorized into mutually commuting pieces via the cycle decomposition. Translating into the physical setup, each factor in the decompositionÛ F = rÛ r corresponds to a disjoint cycle, and each 'cycle' r is a collection of sites r = (r 1 , r 2 , . . . , r lr ), with the following interpretation: after each Floquet period, a spin flip localized at site r i goes to site r (i+1) (mod l r ). The radius of the neighborhood b (r) is then simply max({|r i −r| : i = 1, . . . , l r }), wherer = 1 lr lr i=1 r i is the average location of the sites in the cycle. The localization length ξ is then defined as the maximum of the radii of b (r) as r runs through all the disjoint cycles. As shown in Fig. 3, ξ diverges around s = 1/2 and signals the breakdown of many-body localization.
While ξ is a general diagnostic for many-body localization in any system, in our present construction one can consider another length scale that serves as a proxy of ξ. Pictorially, one can imagine tracking the position of a spin flip as it moves under the drive. Following the discussion above, a spin flip must come back to its starting position after l r Floquet period, where l r is the length of the cycle the spin flip belongs to. Hence, the trajectories of spin flips in the long-time limit define loops on the lattice, and the length of these loops serves as another length scale differentiating between localized and delocalized behaviors. As we will demonstrate below, utilizing this 'loop' picture one can formally map the diluted SWAP model to a classical bond percolation model, where bonds are occupied and unoccupied with probability s or (1−s) respectively. The mapping implies the loop lengths will diverge at s c = 1/2, which is consistent with the divergence of ξ .
The mapping proceeds as follows: For any diluted sequence of SWAP gates (see e.g. Fig. 10b), one can track the orbits (closed trajectories) of a spin flip initially residing on a particular site (Fig. 10c). In the chiral phase, edge trajectories circumnavigate the boundary (green dashed line), whereas bulk trajectories form small loops (blue and red lines). In the trivial phase, the chiral edge trajectories are absent, and all bulk loops are again small. To expose the connection to percolation, it is useful to utilize a dual dense-packed loop representations of classical bond percolation clusters. Consider a dense packed loop configuration in which each bond has an accompanying pair of loop segments, that are either oriented parallel to the bond (Fig. 10d, top panel) if the bond is occupied, or perpendicular to the bond (Fig. 10d, bottom panel) otherwise. Connecting up the loop segments produces a set of closed loops (Fig. 10e), which one can easily verify are in one-to-one correspondence with the trajectories of spin flips in the corresponding diluted SWAP model. As indicated by small arrows in (Fig. 10e), a spin flip initial on a given site executes an orbit corresponding to the loop either down and to the left (A sublattice, filled circles) or up and to the right (B sublattice, open circles).
This mapping immediately shows that the diluted SWAP model has two distinct phases with (s < 1/2) and without (s > 1/2) chiral edges, which survive over a finite range of dilution probability, s. The chiral edge is lost sharply at s c = 1/2, in which the model is characterized by classical percolation exponents associated with a single diverging length scale l ∼ (s−1/2) −ν with ν = 4/3.
However, we caution that the classical percolation character of the transition is likely special to the perfect SWAP gates used in the model, and that more generic quantum perturbations, such as disordering the local strengths of the Hamiltonian, will likely induce quantum fluctuations that change the nature of the transition. An analogous example is that of the quantum Hall plateau transitions, which can be viewed in a loose sense as percolation of quantum Hall droplets. In that case, tunneling between chiral edge states of the droplets causes the universal scaling properties to differ from those of classical bond percolation.

Appendix C: The GNVW index
Here we provide an explicit derivation of some of the claims on the properties of the chiral unitary index ν, and discuss why our exposition is non-rigorous. Note that all the claims have been discussed in GNVW, and so are not original.
We will first show a series of results concerning the 'algebra overlap' η. To this end, recall the definition (with summation convention on repeated indices) where Λ denotes a sufficiently large, but finite, set of sites such that any operatorsâ ∈ A andb ∈ B can only act nontrivially within Λ.
η has the following properties: (1) η is independent of the choice of Λ.
The normalization p a p b /p 2 Λ is precisely chosen for this.
so indeed η is invariant.  By the arguments in the previous claims, we can choose a basis and take Λ = a, and see that Tr a (ê a jm ) 2 = p 2 a . (C4) Next we proceed to discuss properties of the chiral unitary index ν: This is an immediate consequence of the following property of the trace: Tr(a 1 ⊗ a 2 ) = Tr 1 (a 1 )Tr 2 (a 2 ).
The proof of this claim is actually somewhat involved, and we refer the readers to GNVW. In particular, we comment that this property, which is central for the group structure of the classification, is actually quite subtle, and deserves further elaboration.
As we have emphasized throughout, in obtaining a sensible definition of the chiral unitary index ν, it is important that the Lieb-Robinson length LR is much smaller than the system size L. However, this assumption is unpleasing from a mathematical point of view, since whenever L is finite it is easy to find a finite product of locality-preserving unitary operators, sayt L/2 with L even, that violates this condition. Therefore, the classification ν ∈ (log Q + , +) could only be sensible if the thermodynamic limit L → ∞ is taken first. Yet, the many-body Hilbert space is not a mathematically welldefined concept in the thermodynamic limit, and so the formalism presented in the present work is non-rigorous. GNVW overcomes this difficulty by introducing the notion of 'quantum cellular automata', which are automorphisms of the C*-closure of the observable algebras on the infinite chain. This can be viewed as a regularization of the problem by focusing exclusively on local operators (observables), as is natural in a physics problems, and studying only the transformations (automorphisms) among local operators. When these transformations are explicitly written out in a finite system, they correspond to conjugation by locality-preserving unitary operators, and hence our emphasis on the Lieb-Robinson length. As can be seen in the definition, our goal will be to compute η(Ŷ (A L ), A R ) and η(A L ,Ŷ (A R )) when a locality-preservingŶ is given as an MPU. Similar to discussions of matrix-product states (MPS), it will be most convenient to introduce a graphical representation of the equations, as was done in (Fig. 4a). While we annotated each box by a 'M [x] ' in Fig. 4a to emphasize that the tensors are site-dependent, as one would expect in a disordered system, in the following we will drop the annotation for clarity of the diagrams. We will also represent a basis of the operator algebra at site x,ê (D1) Now consider evaluating the index by specifying a cut between sites x and x + 1, and using In evaluating η(Ŷ (A x ), A x+1 ), the nontrivial task is to com- where the repeated indices are traced over. We will discuss below how to evaluate this quantity in the MPU language.
To this end, see that (Fig. 4c) where the vertical dashed line with an arrow indicates the location of the cut, and we have chosen Λ to be the entire system, assumed to be a finite ring of size L 1. Note the † annotation, indicating that the tensors appearing are those corresponding toŶ † . Also note that the tensors at different sites are generally distinct, and thereforeŶ is not translationally invaraint.
One can combine all the tensors for sites other than x and x + 1 into one single tensor, which only has legs in the 'bond space'. Representing this combined tensor as a shaded box, we rewrite where in the second equality sign we performed a singular-value decomposition (SVD) on the shaded tensor. In practice, for sufficiently large L there is only one singular value to sum over -far away from the cut (with distances measured in the Lieb-Robinson length), the operatorŶê ix+1jx+1 is simply the identity (in principle with exponential accuracy), and therefore the evaluation of the trace can be done on a sufficiently large open interval containing the cut. Since such intervals correspond to the OBC, there should only be a single singular value in the SVD of the shaded tensor.
Finally, we evaluate (with summation convention) where we used T and * to indicate that the tensors are those forŶ T andŶ * . In addition, note that there is a further transposition in the bond space of the 'block' of tensors on the left of the figure, as indicated by the flip of the arrow on the dashed line.
One can therefore evaluate ν(Ŷ ) by contracting the tensor in Eq. (D4). However, as discussed in the main text one should choose the intervals L and R to be at least as large as the Lieb-Robinson length in order to obtain the correct index. This can be done by a direct generalization of the above discussion, but with the sites x and x+1 replaced by a collection of sites. As a technical remark, we also note that, for stability, it is important to restore the normalization factor in η in the actual numerical evaluation, for if not the magnitude of the numerical values will diverge as the size of the quantum Hilbert space involved.

Two-site MPUs
To prepare for the construction of example MPUs used in the numerical computation of the chiral unitary index, we first discuss the construction of the building blocks used: random two-site MPUs.
Consider a general two-site unitary operator where U i1,i2 j1,j2 is a p 1 p 2 dimensional unitary matrix. We use the notation that the upper indices are for the leftspace, and the lower for the right. We seek to rewrite it in a MPU form, where we define tensors M

[x]
ixjx , x = 1, 2 and i x , j x = 1, . . . , p x such that i2,j2 |i 1 i 2 j 1 j 2 |. (D6) Since we have a two-site problem, M [1] i1,j1 is a 1×χ dimensional matrix for each pair of i 1 , j 1 , and M [2] i2,j2 is χ × 1 dimensional. Similar to the canonicalization of MPS, we group the tensors into column and row (block-) vectors, and define a p 2 i,j can be found by a matrix factorization ofŨ , whose entries are fully determined by U . Note the crucial fact thatŨ is related to U by a partial transpose, and generally cannot be transformed into one another via any row-column rearrangement (indeed their dimensions do not even match when p 1 = p 2 ). This is important, for if the latter was true the MPU 'Schmidt weights' would all be 1/ √ χ, which is absurd if we are claiming full generality in our construction.
To construct a generic two-site MPU, we (i) generate a random p 1 p 2 -dimensional unitary matrix, (ii) rearrange the matrix elements as per Eq. (D8), and (iii) perform a matrix factorization (say SVD) to find the tensors M [x] i,j . In practice, however, this procedure leads to an MPU representation of an FDLU with an undesirable large bond dimension. To see why, simply note that the number of non-zero singular values in step (iii) above will generically be min(p 2 1 , p 2 2 ), and so the bond dimension of the resulting MPU will scale exponentially with the circuit depth and p. To reduce the computation difficulty, we further engineer the MPU construction process to reduce the resultant bond dimension. Note that we lose the full generality of the constructed MPU once we bring the bond dimension down.
We will do this in two steps: First we will present a particular construction that, while controlling the bond dimension, does not represent a class of two-site MPUs containing the 'generic' gates. Next we will generalize this simple construction to restore generality. Again, we emphasize that generality is only restored at the cost of allowing a bond dimension scaling as p 2 .
The first construction is loosely an analogue of designing a controlled-phase gate. Let {P i | i = 1, . . . , χ 1 } and {Q j | j = 1, . . . , χ 2 } be two sets of orthogonal projectors in H 1 and H 2 respectively, such that χ1 i=1P i =1 1 and χ2 j=1Q j =1 2 . For (slightly) more generality we should choose the basis for the projectors (i.e. the vector subspaces) arbitrarily for the different gates. This can be achieved by multiplying the MPU below, defined with a simple projector basis, by a random on-site unitarŷ U 1 ⊗Û 2 .
We consider the matrices {L i,j ,R i,j | i = 1, . . . , χ 1 ; j = 1, . . . , χ 2 } satisfyinĝ where the omitted cases (e.g. i = i , j = j forL i,j ) are unconstrained. These requirements can be fulfilled by consideringL i,j that act only in theP i subspace, and similarly forR i,j inQ j . Then we claim is an MPU with bond-dimension χ 1 χ 2 (the actual bond dimension maybe even smaller after canonicalization). Indeed, check that Note thatL i,j is j-dependent -if not we will just end up with an MPU of bond dimension 1.
In the above construction, no matter how we partition our on-site Hilbert spaces we will never recover the fully generic two-site MPU. To improve this, we observe that the key points of the above construction are the following: The second condition reduces to the usual unitary condition when χ 1 = χ 2 = 1, and the previous construction corresponds to the 'trivial' solution to this condition, akin to restricting oneself to the on-site unitary subgroup of the two-site unitaries. It is now clear how to generate more general MPU while restricting their bond dimensions -we simply use the more general solution, as discussed near Eq. (D8), for the second condition. Explicitly, we now consider where r i = rank(P i ) and r j = rank(Q j ), i.e. the dimensions of the vector subspaces defined by the projectors, and we demand This gives an MPU with bond dimension χ = i,j min(r 2 i , r 2 j ), which is restricted by how we partition the site Hilbert spaces. At the same time, we recover the fully general case when χ 1 = χ 2 = 1, at the cost of restoring the bond dimension min(p 2 1 , p 2 2 ). The discussion above can be readily generalized to MPUs that act on more than two sites, simply by a regrouping of the DOF into two 'super' sites. In particular, it will be interesting to explore if the above approach in reducing bond dimension, which manifestly preserves unitarity, offers any practical advantage in time evolution simulation of 1d quantum systems. Fig. 5 To demonstrate the computability and quantization of the chiral unitary index ν, we consider some example MPUs taking the form in Fig. 11.

Construction of MPUs used in
The lowest layer (unshaded boxes) represent a 'base' MPU with a known index. Explicitly, we take it to be either the identity operator, or (copies of) the translation operators for different site Hilbert space dimensions p. As discussed in the main text, they can be represented as MPUs.
The upper layers (shaded boxes) represent an FDLU multiplied to the 'base' MPU. While we show an FDLU of depth two in Fig. 11, in certain cases we considered FDLUs of depth four (see Table I). Each shaded box is a random two-site unitary, for which we have detailed their MPU representation in the previous subsection. Note also that, despite the appearance of the figure, the system is not translationally invariant as each shaded box represents a different random gate. ...
... Due to the simple architecture of the circuit, one can also directly read off their Lieb-Robinson lengths: One simply follows the paths upward in Fig. 11 and try to go as far as possible in one direction. For instance, if the 'base' is the identity operator, with our FDLU of depth two an operator local at site 0 can only be nontrivial in the interval [−2, 2], so LR = 2. In addition, the 'Lieb-Robinson light-cone' is strict for such circuits, i.e. the statement thatŶÔ xŶ † is nontrivial only within a distance LR from x is exactly true.
The settings for the various example MPUs used in Fig. 5 are summarized in Table I. As can be seen in the last column of the table, generic two-site MPUs were used in the FDLU for all sub-figures other than (f). For (f), we use two-site MPUs of the form in Eq. (D13), where for each gate (a shaded box in Fig. 11) we randomly pick either the left or the right site as being labeled by '1', and the other by '2'. We consider χ 1 = 3 and χ 2 = 1, where the six-dimensional site Hilbert space for site 1 is partitioned into 6 = 2 + 2 + 2. Since the 'base' MPÛ t (2) ⊗t (3) has bond dimension 6, this brings the total MPU bond dimension down from 6 3 = 216 to 6 × (3 × 4) = 72. Finally, note that in Figs. 5 (b) and (c), the evaluated index for certain cuts converge earlier than the other. This is due to an even-odd effect, imprinted by the arbitrary choice of the 'starting point' of the first layer of FDLU in Fig. 11, and is merely an unimportant artifact of our simple circuit architecture.  Fig. 11. L denotes the total length of the ring used in the computation, and '#. cut' denotes the number of cuts at which the chiral unitary index was computed. p denotes the dimension of the site Hilbert space, 'base' and 'depth' respectively denote the choice in the 'base MPU' and the depth of the FDLU, as illustrated in Fig. 11. LR indicates the Lieb-Robinson length of the MPU. 'Generic' indicates whether or not fully generic two-site MPUs are used in the FDLU.
L #. cut p Base Depth LR Generic (a) 128  In this appendix, we review the key properties of the free fermion AFAI proposed in Refs. [27, 28, and 30]. In equilibrium systems, a chiral edge mode is anomalous, i.e. it cannot exist as a purely 1d system. To see this, imagine a clean single-band 1d model, with the band structure E k satisfying dE k /dk > 0 for all values of k. This implies E π > E −π , violating the periodicity of the Brillouin zone. At the edge of a 2d bulk, however, the contradiction is resolved, as the chiral mode can terminate in the bulk bands, which will necessarily carry nonzero Chern numbers. For a Floquet system with the time-periodic singleparticle Hamiltonian H(t) = H(t + T ), however, the instantaneous Hamiltonian plays a secondary role, and the long-time behavior is governed by the single-particle Floquet operator U F ≡ T e −i T 0 dtH(t) . The eigenvalues of U F are phases, and can be expressed as e −i nT . { n } are known as the quasi-energies, and are only defined modulo 2π/T . Interestingly, since the energy direction is also periodic, the previous argument forbidding a stand alone chiral mode in 1d is now invalid. For instance, the straight chiral mode k = k/T gives π = π/T = −π + 2π/T , which seems to give a legitimate quasi-energy band.
Nonetheless, such a chiral mode is still anomalous. [28] This can be seen by noting that the winding number W = i 2π dk Tr U † ∂ k U is a topological invariant quantized to be an integer. [28] While the single-band model with k = k/T gives W = 1, the identity operator U = 1 corresponds to W = 0. For a purely 1d system, there exists a smooth family of time evolution operators U (t) interpolating between U (0) = 1 and U (T ) = U F . Since the winding number is robust against smooth deformation, the winding number of any purely 1d Floquet system is necessarily zero. [28] Hence, a chiral mode, signified by a non-zero W, remains anomalous, and can only arise as the boundary of a 2d system. In contrast to the equilibrium case, such chiral modes can terminate themselves, and so their presence no longer implies non-zero Chern numbers of the bulk bands. The bulk of the system is therefore Anderson localizable, giving rise to AFAI.
While the definition of the winding number W can be adapted for disordered systems, [30] its generalization to interacting systems is unclear. As a first step towards such generalization, it is instructive to look instead at the many-body Floquet operator corresponding to a chiral edge mode. For the straight chiral mode, one findŝ which is nothing other than the unit translation operator. This simple observation bridges the fermionic AFAI model with the chiral bosonic model we constructed, which also features the (bosonic) translation operator at the edge. Although our current formalism does not immediately apply to fermionic systems, it is suggestive that the fermionic translation operator is as anomalous as its bosonic counterpart. Therefore, to rigorously establish the stability of the fermionic AFAI against introduction of interactions, it only remains to prove that the fermionic translation operator is anomalous, which we leave as an important open question.
phases. From the equilibrium results, one expects that Floquet systems with continuous unitary symmetry G are now (partially) classified by H d+1 Borel (Z × G, U(1)) = H d+1 Borel (G, U(1)) × H d Borel (G, U(1)), with a similar interpretation as before. Of particular interest to us here is, as stated, when G = U(1), for which one finds [31] H d−1 Borel (U(1), U(1)) = H d (CP ∞ , Z) = Z for d even 0 for d odd . (G4) As such, we have i.e. all of them are Z, but the physical interpretation alternates between a 'strong' and a 'pumped' SPT. Now we focus again on d = 1, for which the Z classification above is expected to correspond to a 'pumped' MBL Floquet SPT phase in 1d, in which particles (carrier of the U(1) charge) are unidirectionally pumped in one direction in every Floquet cycle. More explicitly, we relate the above Z classification to projective representations of U(1) × Z. Let (φ, n) ∈ U(1) × Z, with group multiplication in additive notation, and let [(φ, n)] be the (projective) representation of (φ, n). (In the more common multiplicative notation, one writes e iφ ∈ U(1); in our notation here φ is defined mod 2π.) Consider the factor system ω defined via: [(φ 1 , n 1 )][(φ 2 , n 2 )] =ω((φ 1 , n 1 ), (φ 2 , n 2 ))[(φ 1 + φ 2 , n 1 + n 2 )], with the following co-cycle: which is indeed projective iff ν = 0. Now let us rephrase the above using a more physical notation. LetQ be the Hermitian charge operator generating U(1), and using the bulk-boundary correspondence in Sec. III, we write our Floquet operator, which generates the Z factor, for an open 1d spin chain as a product of three unitary operatorŝ where we assumeŷ L (ŷ R ) has a finite support localized on the left (right) edge, andû B corresponds to the bulk operator. Following the above discussion, we consider a system for which the following commutation relations are realized: Note that the first line above is nothing but a translation of Eq. (G8) into the more physical notation. In addition, we assumeû B is MBL with the local charge operators {Q x } being the l-bits, such that one can writê u B = exp(−i r F r ({Q x })) with F r being local real functions. As a consequence, there is a stronger version of the commutation relation: e iφQxû B e −iφQx =û B for all sites x.
By the localization assumption, we can consider the local charge operatorQ L , which is the truncation ofQ to the support ofŷ L . Now consider any stateρ of the system, with an average charge q L ≡ Tr ρQ L on the implying an accumulation of charges on the left edge when ν = 0, which is independent of the choice ofρ.
As argued, such symmetry charge accumulation is common for all 'pumped' bosonic MBL Floquet SPT phases in 1d, protected by a unitary abelian on-site symmetry. Unlike the previous case with G = Z 2 , however, the U(1) charges are unbounded, and so the nontrivial nature of the edge cannot be neutralized by adding subharmonic terms to the Floquet drive. In the long-time limit, therefore, infinite opposite charges will accumulate at the two ends of the chain, regardless of the initial state. Though we haven't proven that such systems are indeed anomalous, they are absurd on physical grounds -while the site Hilbert spaces could be finite-dimensional in the bulk, the edge Hilbert spaces are necessarily infinitedimensional in order for the Floquet operator to remain unitary.
It is therefore suggestive that the cohomology classification of MBL Flouqet phases might be invalid when G = U(1), or more generally when the group G is continuous, despite we have not provided any concrete evidence leading to this claim. For instance, one might argue that we should start with a system having all site Hilbert spaces being infinite-dimensional, say in a quantum rotor model. In that setup, however, it is unclear why the bulk Floquet operator cannot be restricted into a finite dimensional subspace using U(1) conservation.
Another possible resolution is that the discussion is not self-consistent. For instance, it could be that such a system can never be MBL, and hence does not exhibit a bulk-edge decoupling. Given the symmetry group U(1) is abelian, it is rather surprising if there is such an obstruction to many-body localization. Yet, we have already seen a similar obstruction -we argued that the chiral unitary operatorsŶ (ν(Ŷ ) = 0) is anomalous, and hence cannot be localized. In fact, the above discussion, concerning the unidirectional pumping of charges in a stateindependent manner, can be viewed as a U(1)-enriched version of our discussion of chiral unitary operators in 1d. This serves as a potential link between cohomological SPT classifications and chiral phases, and exploration of this link will be an interesting direction for future studies.