Exploring many-body localization in quantum systems coupled to an environment via Wegner-Wilson flows

Inspired by recent experiments on many-body localized systems coupled to an environment, we apply a Flow Equation method to study the problem of a disorder chain of spinless fermions, coupled via density-density interactions to a second clean chain of spinless fermions. In particular, we focus on the conditions for the onset of a many-body localized phase in the clean sector of our model by proximity to the dirty one. We ﬁnd that a many-body localization proximity effect in the clean component is established when the density of dirty fermions exceeds a threshold value. From the ﬂow equation method we ﬁnd that, similar to many-body localization in a single chain, the many-body localization proximity effect is also described by an extensive set of local integrals of motion. Furthermore, by tuning the geometry of the inter-chain couplings, we show that the dynamics of the model is ruled, on intermediate time scales, by an emergent set of quasi-conserved charges.


Introduction
The advent of cold gas experiments [1] has revitalized interest in fundamental questions of quantum thermodynamics in isolated many-body systems. One of the most intriguing avenues of research is the quest for non-ergodic phases of quantum matter. Examples range from integrable models [2,3] to quantum scars [4] and include the prominent example of ergodicity breaking by strong disorder: many-body localization (MBL) [5,6]

Flow Equation Approach For a Single Chain
The key idea of the FE approach is to introduce a family of unitary transformations, U (l), parameterised by a 'renormalization group' scale, l, and generated by the anti-hermitian operator, η(l), via the relation, U (l) = T l exp η(l)dl . The fixed point of the FE procedure in the l → ∞ limit, is a diagonal Hamiltonian with dressed couplings. Operators, O(l), flow according to the equation dO dl = [η(l), O(l)]. A customary procedure for constructing η(l) is to first separate the Hamiltonian into its diagonal, H 0 (l), and off-diagonal, V (l) parts. Then, the generator is constructed as η(l) ≡ [H 0 (l), V (l)] which guarantees vanishing off-diagonal terms at the fixed point, l → ∞ [58]. Typically, the solution of an interacting quantum many-body system via the FE approach would require a broad set of variational parameters keeping track of the nested hierarchies of multi-particles correlations.
However, in the case of MBL systems, a guiding insight in fixing the variational ansatz for the flow equations comes from the l-bit picture [59,60], which provides a method to numerically solve the flow in an efficient way: only the first leading terms describing pairwise interactions between the l-bits are retained, while higher order effects are truncated and discarded. This represents an excellent description as long as the system is strongly localized. Given this ansatz for H(l), the flow of the couplings is readily given by the solution of dH dl = [η(l), H(l)]. In other words, the flow brings the Hamiltonian of a single disordered fermionic wire (for instance, H d in Eq. (2)) into an effectively diagonal one at the fixed point This, in turn, shows that the FE method effectively brings the Hamiltonian into an l-bit basis, with couplings between the integrals-of-motion that decay in space as ∆ ij (∞) ij ∝ exp(−|i − j|/ξ). The values of h i (∞) and ∆ ij (∞) depend on the specific disorder realization. Therefore, to 3 consider disorder averaged quantities, the flow equations must be computed independently for each disorder realization. In addition to extracting the conserved charges and l-bit Hamiltonian in Eq. 1, the FE method can be used to approximate a crossover region from the MBL phase to a delocalized phase [60]. This region is identified with the parameter space where truncation error proliferates. These errors indicate the departure from an MBL phase because they indicate that the true unitary transformation must contain correlations between local degrees of freedom that are not captured by the ansatz. Since the growth of correlation between local degrees of freedom is suggestive of delocalization, the proliferation of truncation error is also indicative of a breakdown of the MBL phase. In order to measure the truncation error, one calculates the so-called 'second invariant' [60,61], a quantity conserved by the exact unitary transformation. Since the truncation breaks the unitarity of the flow, the truncation error is controlled by changes in second invariant.
By setting a small threshold for the change in the second invariant, a tight bound on the MBL phase region can be identified with the parameter space where the truncation yields error within the threshold. Such analysis performed on the single chain led gives a phase boundary consistent with exact diagonalization [60]. We discuss the second invariant in detail as it pertains to the MBL proximity effect in section 3.3.

The model
In this section we extend the flow equation method to the system depicted in Fig. 1. We consider a system composed of two wires of interacting spinless fermions coupled via an inter-chain density-density interaction of strength ∆ I . The Hamiltonian of the system reads (cf. Fig. 1) where the sums run over N s dirty sites in the Hamiltonian H d and over δ × N s clean sites (with δ 1) in the HamiltonianH c . The fields, h i , are drawn from a uniform box distribution of variance W , i.e. h i ∈ [−W, W ]; for sufficiently large W , the chain of fermions, d i , will be in the MBL phase, and will act on the clean fermionic component, c i , as a source of disorder. Even though we study a microscopic model that contains inhomogeneities only in the on-site fields, h i , we write couplings in Eq. (2) with a generic dependence on spatial indices to emphasize that, already at the first steps of integration of the flow equations, couplings can inherit an explicit spatial dependence from the disordered fields.

Two Chain Ansatz
Similar to the FE method for single chain MBL phase, the FE method for the two-chain problem aims to construct a unitary transformation, U (l) = T l exp η(l)dl , that diagonalizes the Hamiltonian Eq. 2. In both cases, an exact calculation would require keeping track of O(2 Ns2 ) matrix elements and is therefore numerically unfeasible. As for the MBL phase in the single chain problem, the local nature of the MBL proximity effect allows one to circumvent this issue via an ansatz for the Hamiltonian, H(l) = U † (l)HU (l) at scale l of the unitary transform. The ansatz we use for the two-chain problem is H(l) = H 0 (l) + V (l) where : A : denotes Wick Ordering [51], the fieldsh c(d) are given below in Eqs. (4), and we will use the convention that the first index in ∆ I ij refers to the clean chain. In the limit l → ∞, V (l) → 0, the fixed-point Hamiltonian, H(l → ∞), is diagonal in an extensive set of l-bits localized on both the clean and dirty sites.
As customary for flow equation methods [60,51], we use Wick-ordered operators, : A :, with respect to a reference state ρ. Wick ordering reduces errors in the truncated Hamiltonian H(l) for the Hilbert space spanned by few particle excitations on top of the reference state ρ [51]. As done 5 by Thomson et al. in [60], we choose a reference state with zero entanglement between local degrees of freedom. This extreme locality condition serves a starting point for the FE unitary transformation to capture the entanglement of the MBL proximity effect. The state ρ employed is a Boltzmann distribution, ρ = 1 Z e −ΘHw , with inverse temperature Θ, chemical potentials fixing particle densities n d and n c , and Hamiltonian The choice of this state allows to easily control energy density, Θ, and particle density distribution, n d .
By Wick-ordering the Hamiltonian at flow time l = 0, the clean and dirty chains pick up effective fields, given byh where their distribution depends on the dirty chain density, n d , the inter-chain coupling ∆ I , and the disorder, W , in the dirty chain. From the expressions of the fields in Eq. (4), it is natural to observe that, if the dirty chain is sufficiently disordered and the inter-chain couplings are sizable, the clean chain will localize as result of the effective disordered field,h c i . Note, the ansatz in (3) has the notational symmetry By exploiting this symmetry, it is easy to derive flow equations for operators of the dirty chain from those of the clean one, and vice-versa. We will refer to terms (or equations) produced by such symmetry transformations using the notion C ↔ D in the following.

Second Invariant and Phase Boundary Analysis
From the ansatz in Eq. 3, we can derive the FE generator η(l) = [H 0 (l), V (l)]. Then, by matching the truncated terms in the Heisenberg equation of motion, dH(l)/dl = [η(l), H(l)], the truncated flow equations, a set of first order differential equations for the couplings, can be derived where the β functions are order three polynomials in the couplings Γ and their forms discussed in detail in section 5. An l-bit Hamiltonian, H(l → ∞), is retrieved by numerically evolving Eq. 7 with initial conditions given by the bare physical couplings, evolved for large l. Deep in the MBL proximity effect phase, these differential equations describe a unitary transform to a diagonal Hamiltonian H(l → ∞), and the unitary transformation described by η(l), along with the Hamiltonian H(∞), can be used to predict dynamics of relevant observables [60,51,62]. When either chain delocalizes, the l-bit Hamiltonian ansatz will be an insufficient representation of the effective Hamiltonian, and the couplings η(l) and H(∞) cannot be used to make 6 predictions. To detect the breakdown of the MBL proximity effect ansatz, we monitor the extent that the truncated flow equations, Eq. 7, break unitarity. For this goal, we employ a quantity known as the second invariant, see, for instance, previous work in Refs. [60,61]. The second invariant is the p = 2 case of a class of many invariants of the FE unitary transformation given by Tr[H(l) p ]. It is particularly easy to calculate for the spin systems and is given by: We can then quantify the error made by a given ansatz by computing the change in the second invariant: If δI is small, then the MBL proximity effect ansatz in Eq. 3, and the approximations discussed above, represent a reliable description and can be used to compute dynamics and the local conserved l-bits. On the other hand, when δI is large, we have an indication that the ansatz fails and that we cannot use the generator η(l) nor the l-bit Hamiltonian H(∞) to make predictions. By identifying a threshold for δI, we can find a tight bound on the phase boundary for the MBL proximity effect. While the choice of threshold is arbitrary, by making it stringently small, one can ensure that below that threshold the MBL proximity effect is properly captured. On the other hand, if it is above that threshold, we must conclude that 1) the system is delocalized or, 2) it is localized in an operator basis not captured by the ansatz. If 2) is the case, then, the operator basis must contain either non-local operators or operators capturing stronger correlations. In either case, a reasonably chosen threshold should yield an approximate boundary for the MBL proximity effect.

MBL Proximity Effect
In this section, we present numerical results, for system sizes unattainable with exact diagonalization, that establish the validity of using an l-bit Hamiltonian to describe the MBL proximity effect. . We study the model introduced in section 3.1 for two equal length chains of length N s = 24 (48 total sites), and numerically solve the differential flow equations, Eq. 7. For this model, the initial couplings are given as: , and clean and dirty number densities are computed with respect to the Wick ordering reference state, n We focus on the limit in which the disordered system would be strongly localized and vary the inter-chain coupling, clean chain hopping strength and reference state parameters. Therefore, we set W = 60, ∆ d = J d = 0.1, and vary the parameters ∆ I , J c , n d (µ d ) and Θ. By setting ∆ c = 0.1, we also focus our attention to the limit in which the clean intra-chain coupling is weak.
The exact form of the truncated flow-equations are given in AppendixC and discussed in section 5. For a fixed configuration of h d i , the truncated flow equations are numerically evolved for a sufficiently long flow-time such that 1) the hoppings, J c(d) ij (l), have become sufficiently small, and 2) there is no appreciable change in the flow of any other coupling. The evolution is repeated for different random instances of h d i , and we present the disorder average of the asymptotic (l → ∞) couplings.
In analogy to a single disordered chain, we define an effective disorder parameter as W c = ∆ I /2J c and work in a limit in which the clean chain is expected to be strongly localized: ∆ I = 45, J c = 0.1, Θ = 0.3 and n d = 0.5 (i.e. W c = 225). We choose such a strong effective disorder to benchmark the method and isolate the effects of varying different parameters. Solving the numerical flow equations (see AppendixD for details on numerical implementation), we find that the density-density couplings, ∆ c(d) ij , are exponentially suppressed in |i − j|, as it occurs in the applications of the Wegner flow to single disordered chains [59,60]. In Fig. 2a, we show the decay in space of the disorder-averaged, asymptotic, density-density couplings, ∆ c ij (l → ∞), on a logarithmic scale, and they illustrate the onset of an MBL phase in the clean chain. As discussed below, the change in the second invariant for these parameters is small for the majority of disorder realizations and thus confirms the validity of the MBL proximity effect ansatz employed in this ansatz.
The top panel of Fig. 2 shows that by decreasing the inter-chain coupling, the final densitydensity couplings between the l-bits present a slower decay in space suggesting a departure from the MBL proximity phase. The effective disorder parameter, W c = ∆ I /2J c , can be used to compare with the disordered Heisenberg chain (a prototype of MBL), which shows a transition at W/J = 4. By considering the second-invariant, we find that the truncation produces minimal error for W c 10 and the MBL proximity is well-established. Note that while we benchmark the method with W c = 225, we found the MBL proximity effect to be consistent with a l-bit ansatz for a reasonable effective disorder strength of W c > 10. While for W c 10, the error grows with decreasing W c and suggests that somewhere in the range W c 10 the system undergoes a transition to a delocalized phase. In this limit, we have found that the final density-density couplings for the dirty-chain,∆ d ij , are still strongly localized while those for the clean-chain are not. This suggests that the source of truncation error is due to the clean-chain becoming delocalized.
The bottom panel of Fig. 2b is one of the most interesting results of our analysis. Here, different curves correspond to different fermionic densities of the dirty component in Hamiltonian (2), with fixed total fermionic density, n tot ≡ n d + n c = 0.5. This variation of n d follows a similar logic to the experiment in Ref. [63], where a complementary situation has been considered (the melting of an MBL phase by coupling to a clean bath). There, the delocalizing effect of the clean component on the dirty component has been experimentally observed in a mixture of collisionally coupled ultra-cold bosons in a two-dimensional optical lattice. Above a certain critical density of bosons, the clean component acts as an ergodic bath and destroys the features 8  of the MBL phase in the dirty sector. Complementary, we find that a critical density of dirty fermions is required in order for the MBL systems to be sufficiently large to entail localization in the clean component. The analysis of the second invariant identifies that the MBL proximity effect is well-established for n d > 0.25, and suggests that for some value of n d less than 0.25, the clean chain goes through a delocalization transition. It is important to note that we are unable to identify with accuracy the point of transition since our ansatz fails close to it (see also Ref. [60]).
We have also studied the effect of increasing the clean-chain hopping, J c and the energy density parameterised by the inverse temperature, Θ of the reference state ρ. We found that the lbit ansatz, Eq. 3, becomes inefficient for large clean chain hopping, J c > 0.5, and at large energy densities, Θ < 0.05. In these limits, the clean chain couplings, ∆ c ij , begin to delocalize while the dirty chain couplings, ∆ d ij , are unaffected. This dependence of localization on the hopping strength is similar to a standard MBL system (the system delocalizes at strong hopping), while  the dependence on the energy density of the dirty chain is novel. At low energy density, the dirty chain charge distribution in the reference state, n d k , and, correspondingly, the effective clean disorder fields,h c k , are strongly disordered, and the clean chain localizes. While for high energy density, the reference state has no disorder in the dirty chain densities, and the clean chain delocalizes. Extrapolating these results, we expect that the localization of the clean chain depends on the disorder of the dirty chain charge distribution.
We summarize our results in the portrait of Fig. 3, which shows the region of the Θ-n d plane where the change in the second invariant is expected to be smaller than our chosen threshold δI c = 0.1. In addition to depicting the trends just discussed, it shows that the dirty chain densities of the reference state must be strongly disorder to compensate for a weaker inter-chain coupling ∆ I (W c ), in order to induce MBL in the clean sector.

Second Invariant
Above we used the second invariant, δI, to identify when the truncated flow equations preserve the unitarity of the exact Wegner-Wilson flow and to justify the MBL-proximity effect ansatz, Eq. 3. Because the flow equation transformation depends on the disorder realization, δI varies from sample-to-sample. The left panel of Fig. 4 shows the distribution of δI for a disorder strength where the MBL-proximity effect ansatz is valid for the majority of disorder realizations, while the right panel shows the distribution for a system where the same ansatz fails for the majority of disorder realizations. In order to distinguish between these two situations, we can compute the median of δI (we don't use the mean because it is artificially biased by the few trials with large second invariant weight). As shown in Fig. 5, the median δI shows that the MBL proximity effect ansatz becomes worse for decreasing ∆ I and n D . Here we see that for W c > 10 and for n d > 0.25, the median second invariant is small and relatively unaffected by changes in W c and n d , demonstrating the validity of the MBL proximity effect ansatz. While for small W c < 10 and small n d < 0.25, the error made by truncation is large and suggestive of a transition to delocalization somewhere below these values. The large sample-to-sample variation of the second invariant suggests the presence of regions not captured by the MBL proximity effect ansatz, and future work may attempt to reduce the second invariant for these disorder realizations by improving the ansatz.

Truncation Error for the MBL Proximity Effect Ansatz
In the previous section we have shown that the l-bit ansatz, Eq. 3, accurately describes the MBL proximity effect phase and that the truncated flow equations, Eq. 7, imply a small error in approximating the exact flow equation unitary transformation, U (l). In this section, we analyze the approximations made by the truncation in Eq. 3, and we discuss, in section 5.2, the physics of the terms contributing to the truncated flow equations. The first type of operators dropped are the n > 2 body terms such as the three body scattering, : c † i c † j c † k c i c j c k :. As long as the integrals of motion do not contain n > 3 body operators with significant weight, then truncating these terms 11 will not produce significant error in the integrals of motion, FE unitary transformation, or l-bit Hamiltonian. This is confirmed by the small second invariant presented in the previous section.
It is important to note that, despite dropping these n-body scattering operators, the ansatz does not ignore all n-body correlations: while at scale l the few body terms are not n > 3 body correlated in the transformed basis, they do contain n > 3 correlations in the physical basis (i.e. U (l) : c † i (l)c j (l) : U † (l) contains n-body operators). In addition to dropping n > 3 body scattering operators from the l-bit ansatz, we drop the off-diagonal terms : n c k c † i c j : and : c † k c l c † i c j :, which we will call correlated hopping and full two-body scattering (F.S) respectively. Including these terms requires keeping track of O(N 3 s ) (O(N 4 s ) for F.S.) number of couplings and significantly increases the computational resources required. To identify the error produced by dropping these terms we highlight how they are produced as the flow evolves.
We identify 7 distinct operators by the 7 sums shown in Eq 3: (see AppendixA for explicit forms for the remaining operators). We then classify contributions to the generator by the type of off-diagonal operator appearing in the commutator: η = [H 0 , J] = η h + η ∆ + η I where: These commutators are computed using rules for Wick ordering [51] and yield: where the coefficients Γ and F are given in AppendixB. The form of the generators are either a hopping operator, : c † j c i :, a correlated hopping (C.H) operator, : n c k c † i c j : or an inter-chain correlated hopping (C.H.I) operator : n c k c † i c j :. It will be important for quantifying the error implied by our truncation to notice that each of the generators is proportional to J c ij or J d ij . In addition, the η ∆ generator is also proportional to ∆ c(d) ij . Taking the commutator [η(l), H 0 (l) + V (l)] yields contributions contained both inside and outside the ansatz, H(l), and are summarized in Table. 1. The operators outside the ansatz are dropped and produce errors proportional to their coefficients. We expect the majority of these coefficients to be small because we study the MBL proximity effect in a limit that the couplings 12  Besides these operators, there are still a few that appear linear in a small coupling and could produce larger error. For example, [η I ,ĥ c ] produces an inter-chain correlated hopping operator, n k c † i c j , which has a coefficient proportional to J c (∆ I ) 2 . While this off-diagonal operator is not small, it is initialized to zero and only affects the diagonal Hamiltonian after commuting with a generator that is also proportional to J c . Therefore, its effect on the diagonal Hamiltonian will remain small as long as J c remains small. This is confirmed by the small change in the second invariant presented above.
This completes our analysis of the error produced by the truncation in the Ansatz, Eq. 3. In summary, we have discussed how we expect that a small error will be produced in our truncation scheme, as long as ∆ c(d) and J c(d) are initialized to small values. We then referenced results in section 4, which demonstrate small truncation error via a small change in the second invariant, to confirm such expectations.

Truncated Flow Equations
In the previous section, we have sketched the derivation of the FE Heisenberg equation of motion, dH dl = [η(l), H(l)], and discussed the error produced by the truncation of the ansatz. In this section we focus on operators in [η(l), H(l)] that contribute to the ansatz and truncated flow equations (Eq. 7). We first focus on the contribution in first row, first column of table 1. For the clean chain it produces a term: and therefore contributes to the evolution of J c ij : This is the primary contribution evolving the off diagonal terms to 0, and is responsible for the intuitive physics discussed above. If we ignore the other contributions to dJ c ij /dl then the evolution of J ij is: Thus, the stronger the disorder in the effective fieldsh c i , the faster the off diagonal terms decay. In addition to producing terms in the β functions that removes the off diagonal couplings J c ij , the generator η h renormalizesĥ c(d) and generates off diagonal hoppings J  k . Eq. 15 shows that the contribution from the commutator [η h ,ĥ c ] removes off diagonal couplings, while Eq. 17 captures new terms produced by the rotation by η h . Similar physics occurs for the generators η ∆ and η I , which are constructed to remove hoppings that change energy via the density-density interaction. In a strong interacting limit, the exact unitaries produced by these generators will generate a Hamiltonian describing doublon and domain wall propagation [64]. If disorder is also strong, these quasi-particle excitation may also localize, realizing a novel MBL of correlated quasi-particles. Unfortunately, in order to capture these effects, one needs to keep track of computationally demanding correlated hopping operators [64] dropped by our ansatz.
While such considerations offer promising prospects for future work, they also have direct consequences for the contributions we include in the truncated flow equations. Since the generators η ∆ and η I transform the hopping operators,Ĵ c(d) , into a set of correlated hopping operators that commute with density-density interactions [64], the truncation above yields a transformation which simply removes the hopping operators without producing the correlated hopping operators they transform into. If these correlated hopping operators are responsible for delocalization, then removing them would produce an artificial localization. To avoid this false localization, we remove the contribution to d dl J c(d) ij coming from [η ∆ , ∆] and [η I , ∆ I ] (respectively, second row, third column; and third row, forth column; of table 1). Ignoring these contributions only produces small error for the same reason dropping the correlated hopping operators produces small error: Figure 6: The dirty chain couples to the clean chain every δ = 3 sites. The emergent integrals of motion are illustrated with different colors: n d k (blue), n c f,r=0 (red) and N f (green).
the error in J ij is proportional to J c(d) ij but its contribution to the l-bit Hamiltonian is (J The small error is numerically confirmed by a small second invariant as discussed above.
The remaining contributions from η ∆ and η I are the ones in the second column of table 1 and describe delocalization processes produced by density-density interactions. A characteristic contribution is: which produce a contribution to the evolution of ∆ c ij as: This contribution captures how the truncated flow equations break the unitary character of the FE transform in a delocalized limit. When disorder is small, J ij remains finite longer during the flow equation evolution and ∆ c ij has a longer time to grow according to the contribution in Eq. 19. This growth produces larger truncation error because, as discussed in the previous section, truncation error is only small when ∆ c ij is small. This concludes our analysis of the physical content of the contributions to the truncated flow equations. The full set of truncated flow equations used in our numerics are reported in AppendixC.

Engineering the geometry of the inter-chain couplings
We now discuss novel effects arising by tuning the coupling geometry. In the geometry of Fig. 6 each fermion of the dirty chain is coupled, every δ sites, to a fermion of the clean chain. This new geometry can still be studied using analogous flow equations to those employed above. Since the clean chain is δ times longer than the dirty chain, we can label the dirty chain with f = 0 . . . N s − 1, and conveniently reference the sites of the clean chain (k = 0 . . . N s δ − 1) with r, using k = f δ + r. f labels the dirty sites, and r = 0 . . . δ − 1 is the number of sites away from the coupled site. We can now explicitly write the initial inter-chain coupling as ∆ f,r,f = ∆ I δ f,f δ r,0 . This leads to an initial clean-chain effective field ofh c f,r = ∆ I n d f δ r,0 . 15 With this important modifications, we can straightforwardly evolve the couplings using the same truncated flow equations discussed in the previous sections. We show evolution of few of them in Fig. 7. The right panel of Fig. 7 shows the suppression of the hopping between a coupled site f, r = 0 and an uncoupled site f, r = 1, while the left panel of Fig. 7 shows the hopping between two uncoupled sites, f, r = 1 and f = f, r = 2, remaining constant. This is consistent with the expectations given by Eq. 16: for a particle to hop on to a coupled site its energy must change by (h c i −h c j ) ≈ ∆ I n d , while such a change of energy is not required for a particle hopping between two uncoupled sites.
With the hopping between uncoupled sites remaining constant, Eq. 19 predicts the divergence of the associated density-density couplings. This is depicted in the left panel of Fig. 7 and explains the failure of the MBL proximity effect ansatz. Instead of modifying the ansatz, we propose to modify the generator, η(l) of the unitary transformation. We define a modified generator η = [H 0 , V ], where we choose V to only include hoppings to coupled sites: Using such a generator, one can employ the same ansatz as above, but the transformation now results in a novel fixed point Hamiltonian describing transport between uncoupled sites and conserved charges on coupled and dirty sites (n c f,r=0 and n d f respectively). In addition, the new generator produces a next-nearest neighbor hopping across the coupled site (i.e J f,r=δ−1,f =f +1,r=1 ). In the proceeding section we derive this hopping rate as 1 τn = [J c (l = 0)] 2 /h c (l = 0), which in the limit of strong inter-chain coupling, ∆ I , is smaller than the other timescales in the system. In this limit, relaxation occurs in two steps: first, on times scales shorter then τ n , transport is blocked by the coupled sites, and second, on times scales longer then τ n , charge is allowed to diffuse across the coupled sites. In the first step, when τ τ n , the system relaxes to a state in which the charge on the bunches of uncoupled sites, N f = δ r=1 n c f,r , is conserved. While on longer times, charge on the uncoupled sites can fully relax via unconstrained transport. In the following two sections, we further investigate this novel behavior by first, in section 6.1, deriving the time, τ n , separating the two relaxation steps, and second, in section 6.2, investigating and deriving the Hamiltonian that governs short time relaxation.

Separation of Time Scales
To derive an estimate of the next-nearest neighbor hopping rate, we first assume ∆ W . This guarantees that the flow of the dirty chain reaches a steady state before there are significant changes in the clean one. We can then treat the clean chain as a single chain with an effective fieldh f,r . We write the new generator as where with J the strength of the hopping on to the coupled site. The first term in η f will suppress hopping between the coupled site and its left neighbor, while the second term will enforce the 16 same on the right neighbor. Since [η f , η f ] = 0 for δ > 2, we can focus on a single coupled site and its neighbor.
We will label the coupled site with 0 and its left and right neighbor sites with −1 and +1. The hopping and effective field couplings will then flow as where J 2 is the magnitude of the next-nearest neighbor hopping, J 2 (l) = J f,r=δ−1,f =f +1,r=1 (l). The flow of these couplings do not depend on the flow of the density-density coupling and can thus be solved independently. We use the assumption that J h and note that the flow of J is much faster than the flow of the other couplings. Thus, assumingh constant, we can approximate the flow of J(l) as Approximatingh 0 as constant, we find Thus, τ n =h 0 J 2 is the characteristic time when relaxation crosses over to full transport and eventually to thermalization. A meaningful separation of time scales therefore requiresh 0 J 2 .
In the following section we will discuss the form of the effective Hamiltonian describing the first stage of relaxation.

Effective Hamiltonian at intermediate times: τ τ n
As discussed above, relaxation in the novel geometry with large inter-chain coupling, occurs in two stages: first, during intermediate times, the model relaxes to a state in which the clean-charge distribution on the uncoupled clusters is approximately conserved, while, on longer times, the clean-charge relaxes to a homogeneous distribution. The Hamiltonian describing the first relaxation process is obtained by dropping the next nearest neighbor hoppings from the Hamiltonian, H(l → ∞). This Hamiltonian, has 3 types of conserved charges as depicted in Fig. 6: the first type, n d k , are the conserved charges on the dirty chain, the second type n c f,r=0 are the conserved charge on the coupled site and N f = δ r=1 n c f,r is the total conserved charge on an uncoupled cluster. For δ > 2 these charges do not determine the dynamics of the charge distribution within an uncoupled cluster, and we must consider the interplay between the intra-cluster tunneling and inter-cluster density-density coupling.
There are two possibilities for such interplay: the density-density coupling between two neighboring sets of uncoupled sites is smaller than J 2 , or it is larger: • In the first case, the density-density coupling can be accurately dropped from the intermediate time effective Hamiltonian. This leads to each set of uncoupled sites, labeled by f , evolving completely independently on intermediate times. The dynamics can be described as the evolution of an effective spin, The local map between the N f fermions on δ − 1 sites and the spin can be performed by identifying the basis states labeled by the eigenvalues of n c f,r =0 with the basis states labeled by the eigenvalues of L z f . Operators that are polynomial in the densities will then be mapped to operators that are polynomial in L z . The remaining terms in the Hamiltonian describe tunneling within a set of uncoupled sites with all the same f . They describe transition between the L z f basis states and are thus described by polynomials in L x f and L y f . • In the second case, when the density-density interaction between the uncoupled cluster is relevant, the local emergent spins will be coupled. Since the hopping operators at a site f commute with those at a site f , a Jordan-Wigner string is not required to correctly reproduce spin statistics, and the coupled Hamiltonian can be written as: where the function F depends on the intra-chain coupling, ∆ c , and the function R depends on h c k , J c ij , ∆ I ij , and ∆ c ij . In general, if the dirty chain or coupled sites have a disordered distribution of charges, the local operators, R f , in the Hamiltonian will be disordered too. 18 The issue of whether the system is fully localized on intermediate times, will then depend on any integrability present in this intermediate time Hamiltonian, or on the impact of disorder on R f .
In the first case, the intermediate time Hamiltonian can be diagonalized by independently diagonalizing the Hamiltonian of effective spins L f . In the second case, when the spins are coupled, further analysis is required to explore the dynamics at intermediate times and will be the subject of Sec. 6.4.

Density-Density interactions between uncoupled clusters
To determine if the effective spins L f are coupled or not, we compute the magnitude of the density-density interaction between two uncoupled clusters. We focus again on one coupled site, labeled by r = 0, and its neighboring sites, labeled by r = ±1 (for any f ). The flow equation equations for the density-density couplings then becomes These coupled differential equations describe a rotation in a three-dimensional space at an instantaneous rate 2J(l) 2 . Given that ∆ c −1,1 (l = 0) = 0, the system (28) can be solved and yields where: l 0 dl 2J 2 (l ) = J 2 (l = 0) h 2 0 (l = 0) (1 − e −2h 2 0 l ).
Therefore, the amplitude of the rotation in such three-dimensional parameter space is small in J 2 /h 2 0 .
We are now in place to discuss which of the two possibilities discussed in the previous section is realized. If J 2 (l = ∞) ∆ −1,1 (l = ∞), then an interacting Hamiltonian describes the intermediate time dynamics while, if the inequality is not satisfied, a non-interacting spin chain will describes the intermediate time dynamics. Given the assumption J h, this inequality simplifies to h ∆. Thus, for the approximation made in ansatz Hamiltonian above, we must choose h > ∆ and conclude that the intermediate time Hamiltonian describes a set of independently evolving spins.
Alternatively, we could assume the bare Hamiltonian has a next-nearest neighbor coupling of the order ∆ −1,1 (l = 0) ≈ ∆ 0,1 < h. In this case the rotation in ∆ c ij space, described by Eq. 28, would still be of a small angle, but away from an initial vector with ∆ −1,1 (l = 0) already greater than J 2 (l = ∞). Intermediate time dynamics would then be described by a set of coupled emergent spins of size L f . 6.4. Explicit form of the Hamiltonian for δ = 3 As an example, we now can consider the δ = 3 case in which there are two uncoupled sites for each dirty site f , and discuss the effective Hamiltonian governing the intermediate time dynamics. The local Hilbert space for these two sites is 4 dimensional and the basis vectors can be labeled by the different ways in which 2 sites may be occupied with particles (the label 1 indicates an occupied site) The local Hamiltonian on these sites reflects the block diagonal structure enforced by the conserved charges:  where∆ L ,∆ R , and∆ R+L are functions linear in the operators n d i and n c f =f and depend on the intra and inter-chain couplings, and fields h c , at the flow time l = ∞. For δ = 3 the conserved charge N f has eigenvalues 0, 1, and 2 that correspond to the three blocks in Eq. 32. This block structure can be represented by two trivial spin-zero subspaces and one spin-half subspace.
We consider the case that N f = 1 for each f , so that the local Hilbert space for the block of interest will be spin-half. The mapping to spin-halves can be preformed via