Anomalous diffusion, prethermalization, and particle binding in an interacting flat band system

We study the broadening of initially localized wave packets in a quasi one-dimensional diamond ladder with interacting, spinless fermions. The lattice possesses a flat band causing localization. We place special focus on the transition away from the flat band many-body localized case by adding very weak dispersion. By doing so, we allow propagation of the wave packet on significantly different timescales which causes anomalous diffusion. Due to the temporal separation of dynamic processes, an interaction-induced, prethermal equilibrium becomes apparent. A physical picture of light and heavy modes for this prethermal behavior can be obtained within Born–Oppenheimer approximation via basis transformation of the original Hamiltonian. This reveals a detachment between light, symmetric and heavy, anti-symmetric particle species. We show that the prethermal state is characterized by heavy particles binding together mediated by the light particles.


Introduction
Many-body localization (MBL) has been a major topic in condensed matter and nonequilibrium physics in the last decade [1,2,3,4,5].Presence of strong, quenched disorder characterizes the typical MBL phase and leads to conservation of local information and suppressed transport.The question under which circumstances other isolated quantum systems can form long-lasting states with preserved quantities is under debate in a wide field of physical scenarios.Apart from Anderson localization [6] as driving mechanism for MBL phenomena, a plethora of systems without quenched disorder have been investigated to provide some form of localization.After the possibility of disorder-free localization has been discussed [7,8], observations have been made in the context of glass transitions [9,10,11,12,13], mixtures of light and heavy particles [14,15,16,17,18], in systems with an extensive number of conserved quantities [19,20], in gauge theories [21,22,23], Stark potentials [24,25,26,27], in systems with strong Hilbert space fragmentation [28,29,30,31] and in an extensive variety of other implementations [32,33,34,35,36,37,38,39,40,41].More recently MBL is discussed in systems which exhibit flat bands (FBs) [42,43,44,45,46,47,48,49,50,51,52].FBs can appear in different kind of systems and can even appear as topologically protected surface states [53,54,55,56,57,58].FBs may remain stable even for strong interactions [59].Like any typical Anderson-like system they show ergodicity breaking characteristics in the strongly disordered limit, but also for very weak disorder [42,43,44] (inverse Anderson transition [60]) or completely without disorder [46,47,48,49].A recent study provides clarification on FB MBL [48].Decoupling of FB Hamiltonians reveals that conserved charges due to compact localized states in degenerate Hilbert subspaces act as effective disorder.Relations to mentioned systems with extensive numbers of conserved quantities and fragmented Hilbert spaces could be drawn.In this work we investigate transport in a typical system with FB.We show that this transport can be diffusive as well as subdiffusive depending on the flatness of the FB.In the limit of very weak dispersions a metastable, prethermal equilibrium emerges [61].Such kind of prethermal plateaus has been observed before in systems with disorderfree MBL [15,28,16,20,21,38,37,17]. We report a noteworthy correlation between FB MBL and quasi MBL in systems with light and heavy particles [16] based on the prethermalization characteristics.We further demonstrate that the prethermal plateau is caused by the circumstance that a light particles species (occupations in dispersive band (DB) states) act as confining potential for heavy particles (FB states).We show that this can be understood in terms of the Born-Oppenheimer (BO) approximation [62].First, we present our system with a single FB coexisting with two DBs and discuss occurring dynamics in sec.2.1 and 2.2.Anomalous diffusion and the arising prethermalization are described in 2.3.General characteristics of the prethermal state, its dependence on system parameters and a treatment in perturbation theory are covered in sec.3.1 and 3.2.BO approximation is carried out in sec.3.3 next to obtain a physical understanding of prethermalization in this model.The arising energy potential for heavy FB occupations within BO approximation plays an important role for the interpretation of previous results.This is studied in sec.3.4.

Model
The Hamiltonian of our system (fig.1a) represents spinless fermions on a diamond lattice with periodic boundary conditions.This and similar structures, often also denoted as rhombic lattice, have been discussed as promising in terms of localization [63,64,65,66,67,59,51].It reads: with fermionic creation (annihilation) operators ĉ † iσ (ĉ iσ ) and occupation number operators niσ = ĉ † iσ ĉiσ .σ ∈ {a, b, c} is an orbital index for the blue, red and yellow lattice sites shown in fig.1a.J terms on solid lines in fig.1a describe nearest neighbor hopping independent of the orbital index.J is also taken as unit of energy and set to 1. J ′ on dashed lines only allows hopping between neighbored a and c orbitals of adjacent cells.The latter terms lead to perturbation of the FB (fig.1b) by giving it a (a) Diamond structure width and thus dispersion.V describes a repulsive particle-particle interaction between nearest neighbors.⟨iσ, jσ ′ ⟩ stands for a summation over all nearest neighbored lattice sites {i, σ} and {j, σ ′ } along solid lines in fig.1a.⟨iσ, jσ ′ ⟩ * sums all connections along dashed lines between next-nearest neighbored a and c orbitals from adjacent cells i and j.
In the course of this article we consider three different sparsely filled systems with L cells, N = 3 • L lattice sites in total, and M ≪ N particles.The large-sized system I has 22 cells and M = 5 particles, the smaller system II contains 8 cells and is filled with M = 4 particles (except for fig.4b) and system III for interpretation purposes consists of 10 cells and M = 3 particles.The total number of particles is conserved giving the Hilbert space dimensions d = N M = (8936928, 10626, 4060) for the given systems (I, II, III).This model contains FB eigenstates inhabiting the space of anti-symmetric a & c orbital superpositions due to mirror symmetry along the central axis.Opposite signs from anti-symmetry cause an effective decoupling of neighbored cells leading to an extensive amount of locally conserved quantities.This phenomenon of compact localized states is shared by several FB lattices [68,69,70].In addition, these eigenstates are quite robust with respect to particle-particle interaction [59] and spread across the energy spectrum [49].Local conservation caused by FB states can be bypassed by adding the perturbation J ′ , as depicted in fig.2b.

Time evolution
The system is initially prepared in a highly localized state represented as a peaked wave packet.The propagation of this wave packet over time across the lattice is modeled with either a completely flat or slightly dispersive band.J ′ is used as parameter to continuously tune the system from a dispersive to a localized regime in order to study the transition into the FB MBL state in terms of propagation characteristics.The initial states fulfill the conditions of dynamical typicality [71,72].It states that pure initial states with similar expectation values for some observable in the beginning evolve alike in that observable over time.Simulating a single or only few of such states drawn at random can be considered typical and therefore representative for a much larger set of pure states.They are constructed in a way described in [73,74,18] and write: where |{n 1a . . .n pτ . . .n Lc } i ⟩ is the full d-dimensional basis in occupation number representation.|Ψ⟩'s are normalized such that ∥Ψ∥ = 1.Projecting the particle number operator of some lattice site {p, τ } on a random vector |Ψ r ⟩ (Gaussian distributed coefficients c i with mean 0, variance 1) generates a state where the localization maximum of the particle density forms a δ peak with ⟨n iσ ⟩ = 1 on site {p, τ }.The other lattice sites are randomly occupied with a many-body background density n bg ≈ M −1 N −1 .A depiction of such an initial state is shown as sketch on the diamond ladder in fig.2a.Averaging of observables from several random initial states might be necessary to obtain smooth results for smaller dimensions (affecting system II).Single particle analysis shows that FB states are compositions of anti-symmetric a and c superpositions with empty b orbitals [59].The peak of the initial wave packet is thus chosen on an a orbital (interchangeable with c) which maximizes overlap with the FB.A δ occupation of an a site in the single particle picture would be generated by equally adding FB and DB excitations, i.e. the FBs contribution is one half.This naive calculus is not accurate for many interacting particles, the idea however proves true and FB contributions are significant in |Ψ⟩.First we consider time evolutions of the described states |Ψ⟩ in position space.Calculations for system I are performed using Krylov subspace methods [75,76,77,78] for matrix exponentials.Smaller systems II and III are fully diagonalized and the time evolution operator is created utilizing the full eigenbasis of H. Fig. 2c-e show the initial part of the equilibration process in system I for different values of J ′ and constant interaction strength V = 1.Times t are given in dimensionless multiples of J −1 .Expectation values ⟨n iσ (t)⟩ are calculated for a single |Ψ(t)⟩.Three different dynamics can be identified: at early times t ∼ J −1 = 1 a ballistic, linear in t, spreading due to DB contributions in |Ψ⟩ proceeds; during early and intermediate times oscillations between a and c orbitals occur until those occupations relax to some FB equilibrium value n loc (noticeable double peak on {p, a} and {p, c}, see fig.2b  and c); for J ′ ̸ = 0, at times t ∼ J ′−1 the decay of the double peak in cell p sets in until some thermal value n th ≈ M N is reached on all lattice sites, see fig.2d and e.The last process is in the focus of this work.Its timescale and scaling with time is controlled by the J ′ value, and thus by the width of the perturbed FB.Fig. 2d shows an example for diffusive spreading caused by J ′ = 0.05 while J ′ = 0.005 in 2e leads to much slower subdiffusive spreading.A categorization of this process into ballistic, diffusive or subdiffusive can be made based on the exponent of its power law scaling with t [79,80,81].This classification is generalized in the next step.

Anomalous Diffusion
To get a better measure for the description of (anomalous) diffusion we introduce the time-dependent variance of the position operator (mean squared displacement) Σ 2 (t) [82].It stands for the broadness of a wave packet spreading over the system.The smaller Σ 2 (t) the stronger it is localized.For our system it reads: where N −1 is the many particle background density as before.Averaging over orbitals σ in each cell is necessary since a, b and c for itself aren't good quantum numbers.The cell averaged occupations n i (t) are normalized such that i n i (t) = 1 at any point in time.Expectation values ⟨n iσ (t)⟩ are calculated for |Ψ(t)⟩ from eq. ( 2).The course of Σ 2 (t) in fig.3a reproduces the dynamics shown in the time evolution of the particle density in fig.2c-e.First, we see a rapid increase from ballistic spreading due to DB contributions of the initial wave packet at short timescales followed by finite size related oscillations.These stem from circulation of the density wave through the system.After that a dynamical regime dominated by J ′ sets in.For J ′ = 0, Σ 2 (t) settles down at some constant, phenomenological FB equilibrium value Σ 2 loc after attenuation of the finite size oscillations.This state, except for small fluctuations, won't change for t → ∞.The strength of J ′ determines how the propagation process between Σ 2 loc and a thermal value Σ 2 th ≈ 1 12 (L 2 − 1) > Σ 2 loc proceeds (Σ 2 th is identical to [18]).This process follows some power law t α with the dynamical exponent α which can be extracted from slopes of the Σ 2 (t) curves in log-log graphs.In the present system α assumes values between 0 ≤ α ≤ 2, which are canonically categorized as localization for α = 0, subdiffusion between 0 < α < 1, normal diffusion for α = 1, superdiffusion between 1 < α < 2 and ballistic propagation for α = 2 [83,84,73,74,18].Note that α is double the value of given exponents in fig. 2 due to squaring of Σ 2 (t).Regular diffusion α ≈ 1 is seen for moderately small J ′ < J.For much smaller J ′ ≪ J we observe subdiffusion until a complete lack of further spreading for J ′ = 0.For small J ′ it is assumed that the simulation time frame can be cut and that observed slopes can just be extrapolated, as demonstrated in fig.3a.The resulting exponent α is then shown for different J ′ in fig.3b.Under these assumptions a smooth logarithmic scaling of α with J ′ is found.Continuous slowdowns of propagation from normal diffusion through a subdiffusive regime to localization have been observed in studies dealing with MBL transitions in disordered systems [85,86,87].Similar behavior has also been seen in quasiperiodic interacting systems [80].Logarithmic dependencies were not reported, though.The root of the fit function f (x) = 1.347 + 0.356 log 10 (x) lies at x 0 = 10 −3.781 being the point at which a phase transition between subdiffusion and localization is to be expected.Yet, this result is strongly influenced by the hard cut-off  (a) Time-dependent mean squared displacement Σ 2 (t) calculated in system I for a single initial state.Decreasing J ′ values slow down the propagation process.Data from fig. 2c-e is contained in (a).Some characteristic power law slopes for ballistics, diffusion, subdiffusion and complete lack of spreading are shown (α = 2, 1, 0.5, 0.25, 0).(b) Dynamical exponents α vs. J ′ for system I. Linear regression fitting of semi-log data shows a logarithmic scaling α ∝ log(J ′ ).(c) Σ 2 (t) in the J ′ → 0 limit for long times in system II averaged for ten initial states.Shaded area indicate a one sigma confidence interval.The separation of a single α slope into two, α pth and α th , with a prethermal plateau Σ 2 pth in between becomes apparent for very small J ′ .(d) Dynamical exponents α vs. J ′ for system II averaged for ten initial states with error bars.α splits into α pth and α th for perturbations J ′ ≲ 10 −2 .Both settle at finite values lim α pth ≈ 0.5 and lim at t = 10 3 ignoring the possibility that the propagation characteristics may change at later times.This renders determination of the diffusion exponent α in the described way inaccurate.Fig. 3a only shows a tiny segment of the whole thermalization process for small J ′ values.The propagation along J ′ doesn't set in before timescales of ∼ J ′−1 are reached and the point in time t th (see fig. 4a and c) at which Σ 2 (t) approaches Σ 2 th scales with J ′−2 .This parameter scaling with time coincides with light and heavy particle species in [16].In order to properly study the J ′ → 0 limit for bands with extremely weak dispersion a correspondingly long time evolution is necessary.This becomes impractical for systems which are not accessible via full diagonalization like system I. Hence we consider the smaller system II in fig.3c, but for much longer times.At perturbation strengths J ′ ≲ 10 −3 visible bends in the slope from Σ 2 loc to Σ 2 th appear which continue to form prethermal plateaus Σ 2 pth with decreasing J ′ .This is similar to several disorder-free systems [15,28,16,20,21,38,37,17] where plateaus of various observables are found in comparable cases of high parameter inequality and long time evolution.The duration of the system remaining in the prethermal state scales with J ′−2 .Since the momentary state of the system at Σ 2 loc and Σ 2 pth violates the eigenstate thermalization hypothesis [49] and preserves information of the initial state, their values are not universal in contrast to Σ 2 th .α is only well-defined for a single slope between Σ 2 loc and Σ  3b and d with aid of fit function f (x) = 1.347 + 0.356 log 10 (x), we find good agreement for J ′ ≳ 10 −3 between α from system I and the combination of α and α th in system II.Yet, in contrast to fig.3b where extrapolation gave a continuous decrease of α → 0 for J ′ → 0, we observe convergence of α pth and α th to finite values lim α pth ≈ 0.5 and lim α th ≈ 0.25 which both lie in the subdiffusive regime.This model does not show a continuous FB MBL transition in terms of band flatness.Localization with α = 0 for J ′ = 0 remains an isolated point and any perturbation J ′ , however small, causes subdiffusion and thermalization in the end.

Dependence on System Parameters
So far we considered a single value for interaction strength equal to the unit of energy V = J = 1.Smaller values of V lead to disappearance of the prethermal plateau.This is shown in fig.4a.A single diffusive process at t ∼ J ′−1 instead of two subdiffusive processes leads to thermalization.As discussed above, t ∼ J ′−1 poses a minimal time while stronger V prolongs the process of thermalization up to t ∼ J ′−2 , see also fig.4c.Despite having repulsive particle-particle interaction the system tends to stay localized for longer the higher V gets which has also been observed in systems with weak disorder [88].Variation of the total number of particles M gives a very similar picture in fig.4b.Prethermalization and subdiffusion can only be seen for at least M = 3 particles.M = 1 and M = 2 curves are akin to low interaction cases in fig.4a.This is surprising in the case of two particles in which an interacting system resembles interactionless dynamics.There are additional necessary conditions which occupations of flat and dispersive bands have to fulfill for prethermalization, as we will see in the following.To further study the influence of interaction strength V on thermalization, we consider the above mentioned point in time t th at which the system thermalizes, as indicated in fig.4a.Linear scaling of t th with V is found between 10 • J ′ ≲ V ≲ 10 • J.The thermalization time is prolonged between J ′−1 ≲ t th ≲ J ′−2 .Together with J ′ t th scales ∝ V J ′2 in this regime.An identical relation has been found in [16] and is a noteworthy similarity between a system with FB and a mixture of light and heavy particles (spins).th in dependence of interaction strength V for two J ′ values.Data stem from averaging over ten initial states.Presence of interaction prolongs the thermalization process between J ′−1 ≲ t th ≲ J ′−2 linearly with diminishing effects in the V → 0 and V → ∞ limits.
However, in the limit of weak and strong interaction its influence on the timescale ends.

Perturbation Theory
Due to very different energy scales and the particularly small parameter J ′ a treatment by means of perturbation theory is helpful.A scheme for similar problems in degenerate translational-invariant spin systems is e.g.described in [90].The fact that we know exact solutions of the model for the system configuration II allows a treatment based on Brillouin-Wigner expansion series [91] which is briefly summarized in the following.Perturbative J ′ contributions are separated from the system (1) which gives two relevant Schrödinger equations: and Perturbation theory sys.II (V = 1, J = 10 −5 ) Σ 2 (t) averaged for ten initial states with shaded areas as one sigma confidence interval.Without perturbation, the system remains in the FB localized state with Σ 2 (t) settling at Σ 2 loc .The prethermal plateau Σ 2 pth is reached in linear J ′ perturbation and remains stable.In second order J ′2 the full thermalization process can be reproduced.
H 0 contains the diamond structure including both J hopping and interactions V .H 1 solely describes the dashed connections between neighbored cells, see fig.1a.J ′ plays the role of the small parameter for perturbation series expansion.The time evolution operator up to the n'th order is constructed as follows: with U containing exact eigenstates |φ i ⟩ from eq. ( 5).Thus, the order of perturbation expansion is determined by the respective order of the eigenvalues in J ′ while eigenstates remain exact.The perturbed eigenvalues E (n) i up to second order read: i ⟩ from eq. ( 4) with the only difference that degenerate subspaces S of H 0 caused by the FB are diagonalized with respect to H 1 .Such a procedure is common in degenerate perturbation theory (see e.g.chapter 5.2 in [92]) and is described in the following.The subsequent matrix is calculated for all eigenstates |φ (0) i,j ⟩ ∈ S: H S 1 is then diagonalized by: from which coefficients (|ψ i ⟩) j are used to rearrange unperturbed eigenstates |φ (0) i ⟩.This yields the states |φ (0)  i ⟩ defined as: which are normalized due to normalization of |ψ i ⟩.This procedure is applied on all subspaces.|φ i ⟩ are defined as: of which |φ (1)  j ⟩ are again diagonalized and normalized with the formalism described above.Σ 2 (t) is shown in different orders perturbation theory in fig. 5.The zeroth order results are identical with the FB localized case J ′ = 0.In first order the system evolves into the prethermal state which is stable within this approximation.Single hopping processes via J ′ are incapable of leading to a further spreading of particle density.In second order the system evolves similar to the exact time evolution.It thus appears that the prethermal state is reached by first order processes in J ′ , while thermalization requires at least second order processes.The confining potential which traps the system in the prethermal state in first order can be found and explained in Born-Oppenheimer approximation next.

Born-Oppenheimer Interpretation
The Born-Oppenheimer (BO) approximation [62] takes a different approach to tackle a scenario like (4) and (5) in which energy scales of system parts are magnitudes apart.It becomes a good approximation for high mass differences similar to perturbation theory when J ′ /J becomes small.Yet, instead of solving only a single Schrödinger equation for the unperturbed part of the system and applying an expansion series to approximate the influence of the perturbation on its eigensystem, it separately solves two Schrödinger equations for both parts in parametrical dependency on each other.We discuss the simplest form of BO approximation for our model neglecting non-adiabatic effects.One prerequisite to utilize this approximation is the existence of differently fast modes which is clearly fulfilled in our model.The large energy ratio between J and J ′ results in vastly different timescales of occurring dynamics.The dynamics of symmetric DB modes is fast at times ∼ J −1 while anti-symmetric FB modes occur at times ∼ J ′−1 .DB states can be interpreted as a light particle species compared to heavy particles from FB states similar to electrons and nuclei in molecular and solid-state physics.Yet, while electrons and nuclei obey lepton and baryon number conservation and occupy distinct subspaces of the Hilbert space, this cannot be said for DB and FB states.In the present (b) BO approximation in system II √ 2J as main junction.J ′ is neglected since J ′ ≪ J.The heavy subsystem contains − orbitals connected by J ′ .Both subsystems interact only via V without particle interchange.(b) Comparison of Σ 2 (t) in the original system (solid lines) and in BO approximation after transformation (dashed lines) for a few different V and J ′ parameters.Σ 2 (t) was averaged over 50 initial states with shaded area as confidence interval.BO data has been smoothed using [89].The courses of the original system, including prethermalization and thermalization at proper timescale, are reproduced within BO approximation with good accuracy.Thermalization times t th in BO approximation fit in with original data, see rectangles.
representation of (1) in fig.1a orbitals a and c can both be occupied by DB and FB states.In order to formulate parametrical dependencies within BO approximation it is crucial to have two separate sets of basis states for both, light and heavy, particle species.A decoupling of light DB modes and heavy FB modes can be achieved using following canonical transformations based on their symmetry [93,48]: They transform the Hamiltonian H from eq. ( 1) to H with: where H0 describes the light part of the system including interaction terms and J ′ • H1 is the heavy part.niσ = d † iσ diσ are again occupation number operators.The transformation changes orbital degrees of freedom in basis states from σ ∈ {a, b, c} to symmetry related σ ∈ {+, b, −} and results in an altered graph shown in fig.6a (see also fig.3b in [48]).+ and b orbitals form a light triangular subsystem ‡ in which √ 2J is the dominant hopping term between sites.The J ′ term in this light system is negligible due to presence of large J. Thus, H1+ in ( 13) can be omitted.The heavy subsystem consists of − orbitals with J ′ as single connection between sites.Both subsystem are effectively decoupled and only interact via V with each other, yet there are no transitions possible which is similar to [16].Occupations in the light and heavy part of the system remain conserved.Transforming states from the {a, b, c} to the {+, b, −} basis can be performed by means of the following operator: ς is defined as the number of particles in a and c orbitals.The conservation of orbitals allows expansion of an eigenstate | φi ⟩ of the full transformed Hamiltonian H into the following product ansatz: with |χ − j ⟩ being a state in occupation number representation from the heavy subsystem and |ϕ + i [χ − j ]⟩ is a corresponding d + -dimensional light eigenstate.The latter one parametrically depends on the configuration of heavy particles which is encoded in |χ − j ⟩.Total dimensions of each subsystem for up to M − heavy and up to M + light particles are given as: ‡ Note that a structure like our light system itself hosts another FB for the specific values of J ′ = ±1 which has also been discussion in terms of FB localization, see e.g.[65,59,94,95] The full Hilbert space dimension d can be restored via The initial step to address the full H (13) in BO approximation consists of first solving the light Schrödinger equation d − times: for all heavy configurations |χ − j ⟩.Interaction terms in H0 cause a repulsive force between light and heavy particles which is reflected in ε The potential seen by the light system resembles finite δ peaks on b orbitals next to heavy particle locations.The light system will quickly go into the groundstate assuming that the vast difference in timescales enables it to instantaneously adjust and relax with regard to any change of heavy particles configurations §.Thus, for the second step of BO approximation only groundstate energies ε + 0 [χ − j ] are taken into account.The second Schrödinger equation for J ′ • H1 reads: The light system remains untouched.Groundstate energies of the light system {ε + 0 [χ − j ]} for all heavy configurations form a d − -dimensional potential landscape through which heavy particles travel with hopping strength J ′ .A few scenarios for system II with different V and J ′ values are calculated in BO approximation in fig.6b by deploying the heavy eigensystem from eq. ( 18) for time evolution: where U contains heavy eigenstates |ϕ − i ⟩.The initial state |Ψ(t 0 )⟩ for BO approximation in {a, b, c} basis is chosen to be in the FB localized phase after relaxation of DB states at some time t 0 < J ′−1 .The operator (14) transforms where time evolution is continued using eq.( 19).Both prethermalization and the process to a final state at thermalization time t th can be observed in this simplified model with reduced complexity and degrees of freedom.The dynamics in the approximated model is solely caused by the movement of heavy particles in the BO potential ε + 0 [χ − j ] which is studied in detail next.

Born-Oppenheimer Potential
The potential ε + 0 [χ − j ] originates from the interaction between the light and heavy subsystem and plays a crucial role regarding prethermalization.While interpretations can be extended to larger systems we restrict the following examination of ε + 0 [χ − j ] to § Non-adiabatic corrections would need to be considered, if heavy particle dynamics causes excitations of the light subsystem, see e.g.[96,97,98].Such excitations are expected to become important for increasing values of J ′ /J.   the simplified system III with L = 10 cells and a total of M = 3 particles for the sake of phenomenological discussion.The number of particles in each subsystem is further fixed to a single light particle, M + = 1, and two heavy particles, M − = 2.A single light particle is required to mediate interaction between heavy particles, while at least two of those are necessary to create nontrivial ε + 0 [χ − j ].The potential for a single heavy particle would be just flat due to translational invariance of the heavy subsystem because of periodic boundary conditions.ε + 0 [χ − j ] depends on interaction strength V and on the relative distance r between heavy particles in |χ − j ⟩: which is why ε + 0 (r, V ) is used as notation from now on.Fig. 7a presents a few characteristic cases for interaction and distance dependence of ε + 0 (r, V ).Considering the inset, a generic r-dependent course of ε + 0 (r, V ) for a high constant value of V almost follows Hooke's law.It is minimal for a clustered configuration of heavy particles with r = 1 and becomes maximized for fragmented r = L/2.In terms of interaction-dependency, BO potentials grow continuously with V , where ) holds for all V > 0, see main part of fig.7a.For strong interactions limit values ϵ 1 and ϵ 2 are asymptotically approached.They ).V = 0 data has been smoothed via [89].The theoretical limit 1 4 (L + 1) for even distributions is marked, but not fully reached due to finite size corrections.(b) Relative distance r pb (t) in system II with same parameters.(b) shares the plot legend with (a).Data were averaged over 10 initial states with shaded area as confidence interval.V = 0 data has been smoothed using [89].Inset: mean squared displacement Σ 2 (t) data corresponding to the main figure .can be derived from calculating the groundstate energy of the light particle being in a box since heavy particles act as infinite barriers in the V → ∞ limit.ϵ 1,2 read (e.g.eq.(B.12) in [92]): The mass of the particle is 1/2 √ 2 in accordance with H0 , eq. ( 13).The smaller ϵ 1 describes a single box with interior size of 2(L − 2), while ϵ 2 is the energy of a particle distributed over two boxes with sizes 2(L/2 − 1).Factor 2 takes two orbitals per unit cell into account (+ and b).A single wide box means lower momentum uncertainty and therefore costs less energy than multiple smaller boxes.The light system generally favors clustered arrangements of the heavy subsystem in terms of energy minimization.This is true for all V > 0, and leads to an effective attraction between heavy particles caused by ε + 0 (r, V ) which we identify as mechanism behind prethermalization in this model.Another way to picture this attraction is presented in fig.7b where an energy landscape ε + 0 (r, V ) in dependence of both heavy particle positions x 1,2 is shown for constant V .
As before there is a minimum at distance r = |x 1 − x 2 | = 1 and a maximum at r = L/2 which are translational invariant due to periodic boundary conditions.The stability of the prethermal state in first order J ′ (fig.5) can be explained by examining possible paths heavy particles are able to pass within this energy landscape.In first order J ′ two heavy particles are trapped at a single r = 1 spot in a valley of stability if J ′ is smaller than the potential barrier: Prethermalization can thus be understood as heavy particles forming an immobile bound pair within the BO potential.Here we find an explanation for the necessity of at least M = 3 particles to observe a prethermal plateau, as seen in fig.4b.At least two heavy particles are required for the bound pair and a single light particle is needed to mediate interactions and to create the BO potential in the first place.Second order J ′2 hopping, however, enables bound pairs to travel along the valley of stability what led to the observed full delocalization in terms of mean squared displacement Σ 2 (t) at times t ∼ J ′−2 .This circumstance of moving bound pairs is reflected in the time evolution of relative distance r pb (t) between heavy particles.It is calculated via: r pb (t) is related to r from eq. ( 20) but takes periodic boundaries into account and is defined for general states | Ψ(t)⟩ in the transformed {+, b, −} basis.Only heavy subspaces with M − ≥ 2 are considered.In presence of attraction due to ε + 0 (r, V ), r pb (t) remains constant since movements within the valley of stability in fig.7b do not change relative distance between heavy particles.This is demonstrated in fig.8a in an idealized case of system III in BO approximation.The initial state of the heavy subsystem consists of two heavy particles next to each other.Time evolution is performed by means of eq. ( 19).r pb (t) grows in absence of ε + 0 (r, V ) for V = 0 in first order J ′ meaning that individual heavy particles are able to move freely through the system, which is in contrast to the V = 1 case.Fig. 8b for system II with identical parameters confirms former results without utilizing BO approximation.Here, we have chosen initial states that are peaked wave packets in the original {a, b, c} basis similar to eq. (2).To account for an inhomogeneous initial distribution in terms of relative distance, they possess an additional density peak at an a orbital in cell p + 1 next to the original peak in cell p.These initial states are transformed into {+, b, −} basis via (14) and evolved using regular time evolution of H. Unlike before, early values of r pb (t) are larger than 1 due to a broader distribution of the initial states caused by the random background in (2).Yet again, an increase of relative distance only occurs for V = 0.It remains constant for V = 1 with only minuscule fluctuations after t > J ′−1 .Regarding the inset of fig.8b, this is in contrast to Σ 2 (t) which shows the known behavior from fig. 3a and 4a.Σ 2 (t) increases in both cases up to thermal Σ 2 th .Full delocalization in terms of Σ 2 (t) occurs despite presence of the bound pairs in the heavy subsystem.Hence, thermalization in the parameter regime of small perturbation J ′ and sufficiently strong V does not imply breaking up of these bonds.This is similar to Hydrogen gas where molecules are stable at room temperature, albeit considered being thermal.This picture of heavy particle binding can also qualitatively explain the behavior of the thermalization time t th as a function of V in fig.4c In general, the barrier height ∆ε + 0 (V ) from eq. ( 22) increases linearly in V for V ≪ 1 and saturates to a constant value for V ≫ 1 like ε + 0 (r, V ) in fig.7a.The crossover occurs around V ≈ 1, because ∆ε + 0 (V ) is determined by the ground state of the light particles inside the static potential of a heavy particle configuration, which is proportional to V .Concerning the thermalization time t th we can now distinguish three regimes: 1. weak interaction V , when ∆ε + 0 (V ) < J ′ , 2. intermediate interaction, when ∆ε + 0 (V ) > J ′ , but V < 1, and 3. strong interaction, when ∆ε + 0 (V ) > J ′ and V > 1.In the weak interaction regime heavy particles do not bind, because their kinetic energy J ′ is larger than the binding energy ∆ε + 0 (V ).Thus, thermalization is dominated by single heavy particle hopping, which means t th ∝ 1 J ′ .In the intermediate interaction regime heavy particles bind to pairs.In this case thermalization occurs via second order hopping processes through the barrier ∆ε + 0 (V ), which leads to t th ∝ ∆ε + 0 (V ) J ′2 ∝ V J ′2 .Finally, in the strong interaction regime ∆ε + 0 (V ) becomes constant.Therefore the second order hopping processes lead to t th ∝ 1 J ′2 .These three regimes and their scaling with V and J ′ are in nice agreement with the behavior seen in fig.4c.Within the picture of heavy particle binding we can also qualitatively understand, why the dynamical exponents α pth and α th differ by a factor of two, as seen in fig.3d.For small values of J ′ the thermalization of this system occurs two well separated steps: on the timescale of t ∼ J ′−1 a single heavy particle moves through the system until it meets another heavy particle to form a bound pair.This process is subdiffusive with an anomalous exponent α pth .Within the prethermal plateau the pair remains immobile until a timescale of t ∼ J ′−2 is reached.Thermalization finally proceeds as a correlated hopping process second order in J ′ .The timescale of pair diffusion is thus squared with respect to single particle diffusion, which leads to half an exponent α th ≈ α pth /2.

Conclusion
To summarize, we first studied general transport dynamics in a system with a FB which gradually receives dispersion.Near the FB MBL limit, i.e.FB dispersion of exactly zero, diffusion becomes slower changing into a subdiffusive regime.Dynamical exponents seem to continuously decrease logarithmically for weaker dispersion.However, when we look at fully diagonalizable systems prethermalization is revealed.The MBL phase remains an isolated point in this model.Any small dispersion causes subdiffusion and thermalization at sufficiently long times.The phenomenon of a prethermal state can be understood in a transformed Hamiltonian in terms of BO approximation.Occupations of FB states behave like heavy particles which are bound by the attractive potential mediated by light DB particles.This binding of heavy particles characterizes the prethermal phase.The BO approximation also sheds light on the thermal state which is made of delocalized bound particles in parameter regimes of strong interaction and weak dispersion.The behavior of the thermalization time and the ratio of two dynamical exponents can be understood as well.The success of BO approximation in explaining many of the observed results regarding subdiffusion and prethermalization gives opportunity for application in other disorderfree systems.A decoupling of degrees of freedom on an eigenstate level as done in this work, may certainly ease approaching thermodynamic limits in many models which are otherwise inaccessible.

Figure 1 .
Figure 1.(a) Sketch and (b) single particle dispersion relation of the examined diamond ladder.

Figure 2 .
Figure 2. (a) Depiction of an initial wave packet with peak in particle density ⟨n pa (t = 0)⟩ = 1 at site {p, a} (red) and random many-body background n bg ≈ M −1 N −1 on all other sites (dark blue).(b) After inner-cell leveling between a and c orbital occupations to some finite value n loc (bright red, see top panel in (c) for comparison), particle density is constrained and only able to escape via J ′ .(c) Numerical results of time evolving a single initial state in system I with J ′ = 0. Absence of J ′ hopping causes FB localization recognizable by a static inhomogeneity in particle density in cell p. (d,e) Transport away from the FB localized state can be in a diffusive or subdiffusive regime depending on the strength of J ′ .Please note the different t axes.
dynamical exp.α 1.347 + 0.356 log 10 (J ) dynamical exp.α 1.347 + 0.356 log 10 (J ) Figure 3.(a) Time-dependent mean squared displacement Σ 2 (t) calculated in system I for a single initial state.Decreasing J ′ values slow down the propagation process.Data from fig.2c-e is contained in (a).Some characteristic power law slopes for ballistics, diffusion, subdiffusion and complete lack of spreading are shown (α = 2, 1, 0.5, 0.25, 0).(b) Dynamical exponents α vs. J ′ for system I. Linear regression fitting of semi-log data shows a logarithmic scaling α ∝ log(J ′ ).(c) Σ 2 (t) in the J ′ → 0 limit for long times in system II averaged for ten initial states.Shaded area indicate a one sigma confidence interval.The separation of a single α slope into two, α pth and α th , with a prethermal plateau Σ 2 pth in between becomes apparent for very small J ′ .(d) Dynamical exponents α vs. J ′ for system II averaged for ten initial states with error bars.α splits into α pth and α th for perturbations J ′ ≲ 10 −2 .Both settle at finite values lim

Figure 4 .
Figure 4. (a) Σ 2 (t) for a range of V values with fixed J ′ = 10 −5 averaged for ten initial states.Shaded areas show a one sigma confidence interval.Rectangles indicate determination of thermalization times t th for (c).Single particle oscillations occur at times t ≲ V −1 .Prethermalization and subdiffusion is absent for small interaction V .(b) At least M = 3 particles are necessary to observe a prethermal plateau.For M = 1 and M = 2 a single regular diffusive process leads to thermalization.Data has been computed and averaged over 1 (M = 1), 1000 (M = 2), 100 (M = 3), 10 (M = 4) initials states.Data for M = 1 and M = 2 have been smoothed via Savitzky-Golay filtering[89].(c) The point in time t th when Σ 2 (t) approaches the thermal value Σ 2 th in dependence of interaction strength V for two J ′ values.Data stem from averaging over ten initial states.Presence of interaction prolongs the thermalization process between J ′−1 ≲ t th ≲ J ′−2 linearly with diminishing effects in the V → 0 and V → ∞ limits.

Figure 6 .
Figure 6.(a) Changed graph after basis transformation.+ and b orbitals form the light subsystem with√ 2J as main junction.J ′ is neglected since J ′ ≪ J.The heavy subsystem contains − orbitals connected by J ′ .Both subsystems interact only via V without particle interchange.(b) Comparison of Σ 2 (t) in the original system (solid lines) and in BO approximation after transformation (dashed lines) for a few different V and J ′ parameters.Σ 2 (t) was averaged over 50 initial states with shaded area as confidence interval.BO data has been smoothed using[89].The courses of the original system, including prethermalization and thermalization at proper timescale, are reproduced within BO approximation with good accuracy.Thermalization times t th in BO approximation fit in with original data, see rectangles.

10 0 10 3 Figure 8 .
Figure 8.(a) Relative distance r pb (t) in BO approximation for a clustered initial state with r = 1 in presence (V = 1) and absence (V = 0) of ε + 0 (r, V ).V = 0 data has been smoothed via[89].The theoretical limit 1 4 (L + 1) for even distributions is marked, but not fully reached due to finite size corrections.(b) Relative distance r pb (t) in system II with same parameters.(b) shares the plot legend with (a).Data were averaged over 10 initial states with shaded area as confidence interval.V = 0 data has been smoothed using[89].Inset: mean squared displacement Σ 2 (t) data corresponding to the main figure.
2 th if J ′ ≳ 10 −2 holds, see fig.3d.Since emergence of Σ 2 pth leads to separation of the propagation process, α is split into α pth for the first slope between Σ 2 loc and Σ 2 pth and α th for Σ 2 pth and Σ 2 th .When comparing both results from fig.