Integrability breaking and bound states in Google’s decorated XXZ circuits

,


I. INTRODUCTION
Recent advances in quantum simulators based on ultracold atoms, trapped ions and superconducting circuits [1][2][3][4][5][6][7][8] have opened a window to studying far-fromequilibrium dynamics and thermalization in isolated many-body systems [9][10][11][12].The behavior of generic thermalizing systems is described by the Eigenstate Thermalization Hypothesis (ETH) [13][14][15] which seeks to explain the process of thermalization at the level of the system's energy eigenstates.In certain systems, the ETH can break down, allowing for new types of dynamical behavior and phases of matter to emerge [16].One of the most striking manifestations of ergodicity breakdown occurs in finely tuned one-dimensional systems [17], which fail to thermalize due to their rich symmetry structure known as quantum integrability [18,19].
A paradigmatic quantum-integrable system is the spin-1/2 XXZ model, which describes the low-energy physics of certain ferromagnetic materials [20].In one spatial dimension, the model's tour de force analytic solution in the isotropic limit was presented by Bethe in the 1930s [21].One remarkable consequence of that solution was a special class of eigenstates that can be viewed as bound states of magnons -the elementary quasiparticle excitations, whose signatures were observed in spectroscopic experiments [22][23][24].However, due to the challenges of probing bound states via conventional techniques such as inelastic neutron scattering, it has been proposed [25] that local quenches [26][27][28][29][30][31][32] may provide deeper insight into the physics of bound states [33][34][35][36][37][38][39].Dynamical signatures of bound states were indeed observed in systems of 87 Rb atoms in an optical lattice, realizing an effective Heisenberg model [40].
While previous studies mostly focused on systems with continuous dynamics governed by a static Hamiltonian, it is also possible to construct equivalent Floquet models defined as a product of unitary matrices.Such models, whose quantum dynamics is intrinsically discrete, are better suited for quantum simulators which operate as a sequence of unitary gates.Quantum circuit models that correspond to the spin-1/2 Heisenberg model in the high frequency limit were studied in Refs.[41,42].Remarkably, the Floquet circuit realization was shown to be integrable for arbitrary parameters and not only in the small time step limit where it reduces via Troterrization to the Hamiltonian model [41][42][43].
The Floquet XXZ model was recently experimentally realized using a ring of superconducting qubits connected by high-fidelity fSim quantum logic gates [44].These qubits interact with each other by superconducting currents and can host excitations in the form of trapped photons.This setup has allowed for the preparation and observation of bound states of a few interacting photons, which were predicted and analytically studied in Ref. [45].One of the advantages was the possibility of controllably breaking the integrability by attaching extra qubits to the main chain and thus changing the geometry of the system.In contrast with the expectation that the bound states are protected by integrability, it was experimentally observed that these states survive even in the non-integrable regime, as previously suggested for the Hamiltonian version of the model [25].However, the robustness of the bound states was not studied in detail and the question of which mechanism protects it in the non-integrable case remains open.
In this work, we use classical simulations, based on exact diagonalization (ED) and matrix product states (MPS), to gain understanding of the experiment in Ref. [44].Specifically, we study the statistical properties of the Floquet spectrum in order to detect the transition from integrable to the non-integrable regime.We also employ time-evolving block decimation (TEBD) simulations to investigate the evolution of bound states and their robustness.In this way, we are able to reach far larger system sizes, photon numbers, and timescales compared to the quantum hardware [44].In contrast to the experiment, which has limitations due to the unwanted leakage of photons, the photon number is conserved in our study.We find that sectors with small but fixed photon number have non-thermalizing spectral properties, which affect both their level statistics and quantum dynamics.Additionally, we confirm the experimental finding that the bound states in these sectors persist beyond the integrable regime.While this effect is pronounced in dilute systems containing small photon numbers, it appears to be restricted to zero density of excitations in the thermodynamic limit.By contrast, sectors with small but finite excitation density are found to thermalize rapidly as the photon number is increased, in parallel with the fast decay of bound states.
The remainder of this paper is organized as follows.In Sec.II we introduce the Floquet XXZ model that will be the main object of our study.In Sec.III we compute the corresponding Floquet modes and investigate the statistical properties of their energy levels, including the average ratio of consecutive energy gaps, the density of states and the spectral form factor.In Sec.IV we study the evolution of bound states and their robustness to integrability breaking.We perform extrapolations to infinite system size and compare the data against the diagonal ensemble predictions, which provides information about the infinite-time limit.In Sec.V we discuss several cases beyond those studied in experiment, in particular systems with a constant filling factor and different decoration patterns, including non-symmetric ones.We summarize our results and discuss their implications in Sec.VI.Appendices provide more details about the corresponding continuous XXZ model, effects of different parameters, and the number of special eigenstates which affect the level statistics.

II. MODEL
The experiment from Ref. [44] has realized a decorated ring of superconducting qubits, schematically illustrated in Fig. 1(a).If the occupancy is limited to zero or one photon per qubit, the photons can be modeled as hardcore bosons.Since we are considering a ring of qubits, we will impose periodic boundary conditions (PBCs) in our ED calculations, unless stated otherwise.The fundamental building block of the circuit is a 2-qubit fSim gate acting on pairs of adjacent qubits, where θ and β determine the nearest-neighbor hopping amplitude and phase, while φ represents the strength of interactions between neighboring qubits.The parameter β mimics the external magnetic flux threading the ring.In the following, we will primarily consider the case fSim(θ, φ, β=0) = fSim(θ, φ).Fig. 1(a) is a sketch of the model with decorations attached to every other site as in Ref. [44].The number of photons will be denoted by N and the total number of sites by L=L sites +L decor , which includes both the sites on the ring L sites and the extra sites L decor .The sketch also depicts a state with N =3 adjacent photons, which will typically be used as the initial state in our simulations.Note that there is another, similar configuration of three adjacent photons, that is simply shifted by one lattice site.This configuration is inequivalent to the one in Fig. 1(a) because it is connected to two decorations instead of one.As specified below, we will occasionally find it useful to average the results over these two initial states.In addition to the layout shown here, in Sec.V we will also consider other decoration patterns.In general, we find the dynamical properties are highly sensitive to the number of photons and the decoration pattern.
Fig. 1(b) shows the corresponding quantum circuit which consists of fSim and SWAP gates.The states of the even, odd and decoration qubits are denoted by |e i , |o i and |d i , respectively.Our classical TEBD simulations follow the layout in Fig. 1(b) and, for convenience, assume open boundary conditions (OBCs).We emphasize that the results below are insensitive to the choice of boundary conditions, as we will demonstrate good agreement between TEBD with OBCs and ED with PBCs.The circuit is defined by first applying fSim gates across all odd bonds, then across all even bonds.Since the even and odd bonds are thus not equivalent, the system is invariant to translation by two lattice sites of the main chain.Additional gates which couple to the integrabilitybreaking extra sites |d i are subsequently applied, which can further reduce the symmetry of the full system depending on the pattern of arrangement of the extra sites.The coupling parameter θ is used to tune between the integrable and non-integrable regimes, while the interaction strength φ = φ is the same, both along the main chain and between the chain and the decorations.( As shown in Appendix A, the Hamiltonian of the XXZ model corresponds to the Trotter-Suzuki expansion of this unitary in the φ, θ, θ → 0 limit.The isotropic XXX version of the model in Eq. ( 2) was first proposed in Ref. [41], while the Floquet XXZ model was formulated in Ref. [42] and analytically studied in detail in Ref. [45].The latter used Bethe ansatz to derive the dispersion of bound states containing an arbitrary number of photons.These bound states are formed by stable magnon quasiparticles, and there are two different phases depending on the ratio of θ and φ: (1) gapped phase φ > 2θ, where bound states of any photon number exist for any momentum, and (2) gapless phase φ < 2θ where the bound states are only present for a finite range of momenta.The maximal group velocity was found to decrease with the number of photons in the bound state.Quantum simulations [44] have later confirmed the analytical relations between the velocity of quasiparticles and their momentum.However, analytical solutions are not available for the non-integrable case, where the integrability is broken either by adding certain perturbations or by changing the geometry of the system.We will use classical simulations to numerically study this regime.

A. Circular orthogonal ensemble
Before we analyze in detail the Floquet spectrum of Eq. ( 2), we must understand the relevant symmetries of the model as they affect the random matrix theory ensemble describing the spectrum after breaking the integrability [47].For example, the undecorated model with PBCs is invariant under translation by two sites, due to the even and odd layers of fSim gates being applied separately.Such a circuit is also invariant under spatial inversion.However, attaching extra qubits to some of the sites reduces the symmetry of the full system.Regular patterns with decorations on every nth site will preserve some form of translation invariance, although with a larger unit cell.Furthermore, the system can be inversion-symmetric only if the decoration pattern itself is also inversion-symmetric.However, in some cases, such as that with decorations on every other site, the inversion of the decorations can be incompatible with the inversion of the main ring due to different reflection axes, so the full system only has translation symmetry, even though the decoration pattern is still inversion-symmetric.This will be discussed in more detail in Section V B.
For a general unitary matrix ÛF , the level statistics is expected to conform to the Circular Unitary Ensemble (CUE).However, as will be apparent in Secs.III A and V B, in most cases studied here we obtain the Circular Orthogonal Ensemble (COE) statistics instead.COE would trivially ensue if ÛF = Û T F , however this is not the case here for any arrangement of decorations.Our calculations show that the necessary conditions for COE level statistics are an inversion-symmetric decoration pattern and equal parameters for the fSim gates on the even and odd bonds along the ring, as defined in Eq. ( 2).Additionally, the mirror axis for inversion needs to be centred on a site, not on a bond between two sites.If R is the inversion symmetry operator which reflects the qubits along this axis, we well then have R Ûodd R = Ûeven , R Ûeven R = Ûodd and R Ûdec R = Ûdec .For simplicity, we define a modified one-cycle unitary operator The operators ÛF and Û F have the same spectrum, since they differ only by a time shift.It is now easy to see that Û F = R Û T F R. This can be understood as an additional symmetry which relates the evolution operator and its transpose, resulting in COE level statistics.Our situation is reminiscent of Ref. [48], where the Floquet spectrum was shown to have COE instead of CUE statistics if there is a transformation which connects the two steps of the Floquet unitary.
Another possibility is when the mirror axis is between two adjacent sites, leading to R Ûodd R = Ûodd and R Ûeven R = Ûeven .We will then have Û F = R Û F R, meaning that R is simply another symmetry of Û F which needs to be resolved.The level statistics in the sector where R has eigenvalue +1 is then CUE.We only find deviations from this expectation for small numbers of decorations such as two or four adjacent decorations, where the level statistics after resolving the R symmetry is somewhere between COE and CUE.However, it seems to increase towards CUE as the density of decorations or the number of photons is increased.There are also special cases which are inversion symmetric in respect to both types of mirror axes, such as the pattern with decorations on every third site.In those cases the level statistics stays COE even after resolving the R symmetry.In contrast, all non-symmetric decoration arrangements were found to exhibit CUE level statistics.

III. SPECTRAL PROPERTIES
In this section we analyze the spectrum of our unitary circuit model in Eq. ( 2).This model does not have a Hamiltonian representation in the general case, since the mapping to the XXZ model (Appendix A) is only valid in the dt→0 limit.As a consequence, the system does not have eigenstates in the usual sense.However, we can instead compute the eigenstates of the one-cycle evolution operator ÛF (2), which are known as the Floquet modes.The corresponding Floquet quasienergy spectrum is periodic with periodicity 2π/T , where T is the time length of one cycle.We set the units such that T = 1.We will investigate the properties of the Floquet modes and quasienergies from two complementary perspectives.On the one hand, we will study the level statistics and density of states, which directly derive from the quasienergies and thus tell us about the behavior of the system at very late time scales corresponding to the Heisenberg time.On the other hand, we will contrast these results against the spectral form factor, which provides information about intermediate time scales relevant for transport, such as the Thouless time.

A. Level statistics
In order to determine whether our model Eq. ( 2) is integrable or chaotic, we study the statistics of its quasienergy levels.In particular, we examine the level statistics ratio, r = min(s n , s n+1 )/ max(s n , s n+1 ), characterizing the spacing of adjacent quasienergy gaps s n = n+1 − n [49].An integrable system is expected to follow the Poisson distribution with the average value r P ≈ 0.386, while in the chaotic regime the expected distribution for our case, as explained in Sec.II A above, is the Circular Orthogonal Ensemble (COE) with r COE ≈ 0.527 [47,50].We vary the hopping amplitude θ between the main chain and the extra sites from 0 to π and plot the corresponding r (θ ).Fig. 2(a) shows the results for N = 3 photons for various chain lengths, while the extrapolation to an infinitely large system L→∞ is plotted in the inset.This result should be contrasted against the results for N = 4 and N = 5 photons in Fig. 2(b).
Turning on the coupling to the decorations is expected to break integrability, which can indeed be observed in Fig. 2 where the value of r (θ ) rapidly jumps towards r COE as soon as θ = 0.For N ≥ 4 photons, as soon as θ 0.05π, the level statistics becomes pinned to the COE value, in agreement with the usual expectation for integrability breaking in Hamiltonian systems [51].However, the case with N = 3 photons shows a visible departure from these expectations, exhibiting pronounced dips towards the Poisson value at special values of θsee Fig. 2(a).Furthermore, we find that the positions of the dips in r depend on the main chain hopping amplitude θ, but not on the interaction strength φ or the flux through the ring β, see Appendix B. No emergent In all the plots, we resolve the translation symmetry and consider only the k=0 momentum sector.
symmetry which would explain the dips at certain values of θ could be identified.Instead, we will relate the presence of dips with special structures in the density of states in Sec.III B below.We note that in all cases plotted in Fig. 2, the value of r (0) lies below the Poisson line, even though the model is known to be integrable at θ = 0.This is simply due to a large number of degeneracies present in the Floquet spectrum, which originate from the decorations and produce a peak in the probability distribution for zero level spacing.Even though there is no hopping to the extra sites when θ = 0, there are still states where one or more photons are frozen in these additional sites.A state with all photons outside the main chain has zero energy, as do some states with two separate photons on the main chain and all other photons outside.We found that completely removing the extra sites brings r (0) closer to r P .We also note that, while the hopping amplitudes inside the main chain and between the chain and the extra sites are different, θ = θ , the nearest-neighbor interaction strength is equal in both cases φ = φ .This means that the photons frozen in decorations can still interact with the other photons.
Finally, we emphasize that in Fig. 2 we assumed a fixed decoration pattern chosen in Ref. [44], where an extra site is attached to every other main-chain site, such that the system is invariant to translations by two sites.In Sec.V B we consider other decoration patterns, including the case of a single decoration, three decorations on the second, fourth and sixth site, and a random arrangement of decorations on half of all sites.All of these patterns break the translation symmetry and do not exhibit oscillations in r (θ ) that are visible in Fig. 2(a).Instead, r first reaches a plateau and then starts to decay at larger values of θ .The plateau is at r COE for inversion-symmetric patterns and at r CUE for nonsymmetric ones.We have also considered patterns which preserve some form of translation symmetry, for example those with decorations on every site, or every third, fourth or fifth site.The level statistics for these patterns shows similar properties to the previous case of decorations on every other site, with deviations from r COE for N = 3, albeit with minima and maxima in r at different locations.In contrast, no such deviations were observed for N ≥ 4.

B. Density of states
The intriguing features in the level statistics observed in Fig. 2(a) can be understood from the density of states (DOS).For example, sharp peaks in DOS signal a large number of degeneracies in the spectrum, which can decrease the value of r .In Fig. 3, we plot the normalized DOS curves for N = 3 and N = 5 photons at several values of θ that were marked by A-E in Fig. 2(a).Both photon numbers exhibit a peak at =0 when θ =0, which is explained by the previously discussed large number of zero modes due to the extra sites.This zero-energy peak is much more prominent for N = 3 and its relative height decreases with N .The results for N = 4 (not shown) are in between those for N = 3 and N = 5, with more peaks than N = 5, but still overall flatter than N = 3.As θ is increased, the DOS curves become more flat.However, several other notable peaks are present for N = 3.Although these peaks are visible at all θ , they are particularly sharp at those values where r deviates from r COE (e.g.θ ∈ [0.25π, 0.45π] and θ ∈ [0.75π, 0.85π]), see Figs. 2(a) and 3.The peaks in DOS are not present for non-symmetric patterns of extra sites, such as just one or three decorations, which will be discussed in Sec.V B.
The sharp peaks in DOS can be attributed to the existence of special eigenstates with a relatively simple structure, which can be built by combining single-photon and two-photons states.Analytical expressions for the dispersions of a single photon or N -bound photons in the integrable (non-decorated) circuit are known [45].Adding the decorations with θ =0 results in an additional zeroenergy band in the single-photon dispersion, since all the photons in the extra sites are frozen.This means that, for example, single-photon and two-photon eigenstates are still present in the three-photon spectrum at θ =0, since we can just move the remaining photons to the decorations, where they will have zero energy.In Fig. 4(a), we compare the actual two-photon Floquet spectrum (dots) with the states constructed from two single-photon states (crosses).The color scale represents the deviation from the nearest analytically constructed state.The agreement is remarkably good, which is not surprising given that the system is very dilute and only nearest-neighbour interactions are present.The two bands at the bottom of the plot are two-photon bound states, which are also in agreement with analytical expressions from Ref. [45]. Figure 4(a) is for the integrable case with no extra sites.Adding the decorations leads to the appearance of two additional bands, see Fig. 4(b) and compare with (a).The first one is a bound-state band, which corresponds to one photon in the main chain and another in adjacent decoration and is completely flat.The second one is a wider band of single-photon states corresponding to one photon inside and the other in a non-adjacent decoration.This wider band is centered around zero and has high DOS on its edges, which coincides with the peaks around ±π/3 in three-photon DOS from Fig. 3(a).Another smaller peak in DOS around −0.75π comes from the flat band of two bound photons.
We can conclude from the previous discussion that the three-photon DOS is strongly influenced by special single-and two-photon eigenstates.This effect is not so prominent in DOS for 4 or more photons, likely because the number of special states is much smaller compared to the total Hilbert space size.In Appendix C we quantify this and show that the proportion of special states for a fixed photon number N becomes asymptotically independent of the system size L.However, the saturation value still strongly depends on N , e.g., the special states comprise as many as 70% of all states for N =3 but only 1% for N =8.This analysis can now be extended to finite values of θ .The analytical expressions for the single-photon dispersion are not available in this case, but can be easily numerically computed for different coupling strengths θ .There are still three different bands, since each unit cell contains three sites.This numerical data can be used to construct three-photon bands, which correspond to three separate particles.This is a good approximation in a dilute system, even with non-zero interaction strength.In this way we obtain ten different bands, e.g., all three photons in the first band (denoted by 111), two photons in first and one in second (112), one photon per each band (123), etc.The dependence of these bands of special states on θ is shown in Fig. 5(a).As θ is increased, the bands move and cross each other.The DOS is typically higher near the edges of the bands, so we expect the DOS to be amplified when two bands cross.The edges of 222 and 123 bands overlap around θ =0.35π, which is where the level spacing ratio deviates the most from r COE .Several other bands also overlap around this point.The DOS plot for the special bands shown in Fig. 5(b) roughly corresponds to the peaks in Fig. 3(a).Therefore, as in the θ =0 case, the peaks in DOS at θ =0 are also explained by the special states which comprise a large proportion of the Hilbert space in systems with smaller numbers of photons, such as N =3.However, it is not obvious from Fig. 5(b) in which θ regions the level statistics deviates the most from values expected in chaotic systems.In particular, there are three very prominent peaks around θ =0.65π,where the level spacing distribution is actually very close to COE.We conjecture that these peaks are not narrow enough to lead to a sufficient number of degeneracies that could affect the level statistics.One might expect that the specially constructed states are a better approxima- tion for a non-interacting system and that the r (θ ) dependence would look different at smaller values of the interaction strength φ.This, however, is not the case, as shown in Appendix B, where it can be observed that the level statistics barely changes with φ.

C. Spectral form factor
The level statistics quantities considered above derive from the properties of eigenvalues of the Floquet unitary, hence they describe the behavior of the system at late times.In order to gain information about intermediate times, we study the spectral form factor (SFF) [52]: which is defined in terms of two-point correlations between Floquet quasienergies, n .As we set the time period of one unitary cycle to T =1, the time in the above equation is equal to the number of cycles, t=n c .The SFF is known to behave differently in integrable and chaotic systems, see Refs.[53][54][55][56] for some recent examples.In both cases, the SFF behavior at short times is governed by microscopic details of the system and therefore it is non-universal.After this initial transient, in integrable systems (Poisson ensemble) the SFF stays approximately constant around the value equal to the Hilbert space dimension H, K P (t) ≈ H.In nonintegrable systems, it first reaches a global minimum and, around the Thouless time t Th , it starts to grow ap-proximately linearly, according to the predictions of random matrix theory, until it saturates at H around the Heisenberg time t H ∼ H.The level statistics and density of states studied previously naturally pertain to the times of order t H , where the discreteness of the Floquet quasienergy spectrum is resolved.
The SFF is typically noisy and suffers from a lack of self-averaging [57,58].In order to smoothen its time dependence, we choose to average it over the flux through the ring β.This parameter does not qualitatively affect the level statistics, as confirmed in Appendix B. Additionally, after averaging over 100 values of β ∈ [0, π], we also compute the moving average at each time point by taking into account the nearest 60 points, which finally results in relatively smooth curves.The averaged SFFs for N = 3 photons, L = 60 + 30 sites and different values of θ are shown in the inset of Fig. 6(a).After an initial period of non-universal behavior, the SFF for θ = 0 assumes an approximately constant value, confirming that the system is integrable.In contrast, a clear linear ramp followed by saturation emerges for all studied values of θ > 0, consistent with broken integrability.We note that the SFF for θ =0 saturates at a higher value than θ >0, where the plateau is exactly as expected at H. The reason for this is a large number n 0 of zero modes in the integrable case, which increases the late-time value of the SFF to H + n 2 0 (at θ =0).
Furthermore, the Thouless time t Th can be extracted from the SFF data.This time gives us the onset of the universal behavior described by random matrix theory (i.e., the linear ramp).The COE prediction for SFF in the time window 0 < t < H is [47] K COE (t) = 2t − t ln(1 + 2t/H), (5) as shown by the dashed black curve in Fig. 6.In theory, the Thouless time could be defined as the smallest time for which K(t) = K COE (t).However, since K(t) is typically not smooth enough even after averaging, in practice we use the following criterion to determine the Thouless time [59] ln(K(t Th )/K COE (t Th )) = 0.4.
The precise value of the filtering parameter 0.4 is unimportant, as long as it is finite but not too small.In Fig. 6, we plot the extracted Thouless time together with the average level spacing ratio r for N =3 and varying θ .Interestingly, the two curves exhibit very similar features, which means that the previously observed deviations in the level statistics for N =3 leave an imprint in the thermalization properties of the system, i.e., thermalization occurs later in systems which are farther away from the non-integrable case.
The agreement of the SFF with the random matrix theory prediction K COE is not particularly good for N =3 photons.This supports our previous observation that the energy spectra in small photon number sectors have special properties, e.g., as seen in the oscillations in level statistics and non-monotonic DOS.The agreement with COE becomes better as the number of photons N increases, as can be seen for N = 5 in Fig. 6(b).The number of zero modes at θ =0 is now much smaller than the Hilbert space dimension, so the dashed horizontal lines at H and H + n 2 0 are visually indistinguishable.There is also less variance in K(t) curves for different values of θ >0, which is reflected in the almost constant value of the extracted Thouless time, as shown in the same figure.This is in line with the level spacing ratio r , which shows no oscillation with θ for this photon number but instead remains approximately constant around r COE .

IV. BOUND STATES
Thus far, we have focused on generic aspects of thermalization at the level of the entire Floquet spectrum of the decorated XXZ circuit in Eq. ( 2).However, one of the motivations behind the experiment [44] to study this particular model is the fact that its integrable version hosts a special class of ballistically propagating bound states.While such bound states are here protected by integrability, they represent only a fraction of all eigenstates and therefore it is not hard to imagine that they may persist, due to some other protective mechanism, after integrability is broken.We now examine in detail the stability of such states after decorating the circuit to break its integrability.
In the Ising limit of the Hamiltonian version of the XXZ model, J z J, an N -particle bound state corre-sponds to N adjacent spins being flipped [25], Even far from the Ising limit when J z ≥ J, the behavior of an N -particle bound state can be understood by starting from such an initial state, which is no longer an eigenstate.The same is true for the Floquet XXZ circuit, with the Hamiltonian Ising limit corresponding to φ θ.Starting from the initial state (7), the "bound state probability" (BSP) after n c cycles is given by where ) is the probability of finding photons in N adjacent sites, where the indices i and j label the sites on the main chain.Conversely, the probability of any other N -photon configuration was denoted n S .The BSP defined in this way was experimentally measured in Ref. [44] and was found to gradually decay over time even at θ =0.This decay was due to experimental imperfections rather than an intrinsic property of the model.For an ideal implementation of the XXZ circuit, the BSP drops rapidly before fluctuating around a steady finite value, as will be shown below.However, once integrability breaking terms are introduced into the Floquet circuit, there is no requirement for the N -photon bound states to continue to be stable at late times.Below, we will focus on understanding the effect of integrability breaking on BSP dynamics using TEBD simulations implemented in iTensor [46].Subsequently, we will show that other observables can reveal a signature of bound states by probing the memory of the initial configuration in Eq. (7) as the system evolves in time.Finally, we will interpret these results by examining the structure of the Floquet modes, in particular their overlap on the initial state in Eq. (7).

A. Dynamics of bound state probability
The BSP dynamics for N =3, 4 and 5 photon bound states is presented in Figs.7(a)-(c) for various strengths of the integrability-breaking θ .By increasing θ from 0 to π/2, the decorations become more strongly coupled to the main chain and the bound states are eventually destroyed.However, at intermediate θ the BSP does not appear to decay to zero, even after many cycles.This is true even when θ is comparable in size to the natural energy scale along the chain, θ ≈ θ.For larger bound states, θ introduces large, slow oscillations into the BSP that are independent of system size.The origin of these oscillations will be explained in Sec.IV C.
Typically, an infinitesimally small perturbation is sufficient to destroy integrability in the thermodynamic limit L→∞ and infinite time limit t→∞.We access these limits by extrapolating the numerical data for the BSP via We also averaged over two inequivalent initial states (7).Inset: Average BSP for different photon numbers N at fixed L=300+150.Data was averaged over 150 cycles and obtained using TEBD with bond dimension χ = 320 for 3 ≤ N ≤ 12, θ = π/6, φ = 2π/3 and varying θ .(e) Diagonal ensemble prediction for the probability to remain in a bound state, averaged over two possible initial configurations.Data in this panel is obtained using ED with PBCs, on system sizes N =3, L=30+15 (with Hilbert space dimension dim=14189); N =4, L=20+10 (dim=27404); N =5, L=14+7 (dim=20348); N =6, L=12+6 (dim=18563).Inset: Extrapolation to infinite system size for a translation invariant initial configuration with k=0 momentum.
two methods: time-averaging the TEBD results and evaluating the diagonal ensemble predictions from ED data.The latter directly takes the t→∞ limit by assuming that the off-diagonal elements of the density matrix average out to zero [15,50].These two methods have different advantages and limitations.While the TEBD method allows us to study dynamics in very large systems, these simulations become computationally more expensive as the evolution time increases, which limits the total number of cycles.On the other hand, the diagonal ensemble prediction provides information about the BSP at infinite time, however, it requires a computation of the complete eigenspectrum using ED, which limits the maximal system size.The total Hilbert space size is constrained by the amount of RAM available for diagonalization, while our implementation relies on 128-bit integers to represent basis configurations, which limits the maximal number of sites L ≤ 128, irrespective of the photon number N .In principle, the latter restriction can be lifted using a more flexible encoding of the basis states, at the cost of sacrificing some of the computation efficiency.
For the TEBD time average of the BSP, we consider 100 cycles between cycle n c = 20 and cycle n c = 120 for a variety of system sizes ranging from L = 20 + 10 to L = 300 + 150.In this way, we exclude the data at very short times which may be impacted by non-universal effects.By fitting the average BSP at each system size according to BSP(L, θ ) = α(1/L) + BSP ∞ we extrapolate to L → ∞ and obtain the result plotted in Fig. 7(d).This procedure was repeated for several photon numbers N .In each case, the initial state was chosen according to Eq. ( 7), which is not translation-invariant.As discussed in Sec.II, there are two such inequivalent configurations and our results are averaged over both.We find that the bound states are robust for a finite range of θ which decreases as the size of the bound states increases.
We also address how the robustness changes as the bound states increase in size, but continue to be dilute relative to the total system size, N/L 1.We calculate the BSP for bound states between sizes N = 3 and N = 12, averaged over n c = 150 cycles, to find BSP(N, θ ).The results for a fixed number of sites L = 300 + 150 are plotted in the inset of Fig. 7(d), where it can be seen that the BSP(N, θ ) curves are starting to converge for larger values of N .These results suggest that large but dilute bound states continue to be robust.The diagonal ensemble results, which directly access the infinite time limit t→∞ of the BSP, can be seen in Fig. 7(e).These results are consistent with the extrapolated TEBD results, suggesting small bound states are robust for a finite range of θ .In particular, the N =3 bound states appear to be robust up to values of the integrability breaking that are comparable to the onchain hopping terms.As N becomes larger the bound states appear to become less robust but both Fig. 7(e) and Fig. 7(d) suggest the N =4, N =5 and N =6 bound states are robust for a finite range of θ .
For our diagonal ensemble calculations in Fig. 7(e) we also averaged over two inequivalent initial configurations (7).Since these states break translation invariance, we have to work in the full Hilbert space and therefore cannot obtain enough data points for reliable system-size scaling.However, if we form a translation-invariant initial state, we can restrict to the k=0 momentum sector and reach much larger system sizes.This allows us to extrapolate the diagonal ensemble value for BSP to L→∞, as shown in the inset of Fig. 7(e).These results suggest that the bound states of N =3, 4 and 5 photons are robust at moderate values of θ , which are clearly in the non-integrable regime according to the level statistics in Fig. 2, even in the infinite time and infinite size limit.
Although the results for translation-invariant and noninvariant initial states in Fig. 7(e) and inset are qualitatively similar, there are also some minor differences.Most notably, the BSP decays quadratically with θ for the symmetric state and linearly for the non-symmetric initial state, suggesting that symmetric states are more robust to integrability breaking.This is not surprising, given that Floquet modes have well defined momenta due to the overall symmetry of the system, and therefore a translation symmetric state can have a higher overlap with a single mode than a non-symmetric one.Nevertheless, these differences appear to rapidly diminish with an increase in N .

B. Memory of the initial state
The BSP is not the only local observable that reveals the unusual behavior of the bound initial states at finite θ .Persistent nonthermalizing behavior can also be seen in the site occupation, n i = ni .Since the integrability breaking decorations make up one third of the total sites on the chain, we would expect a third of the photons to be located on them after a short time when the system has sufficiently thermalized.Instead, we find that this is only true at larger θ ( π/3).The fraction of photons located on the decorations as θ is varied shows very similar behavior for bound states of different sizes.
In the integrable case, larger bound states propagate more slowly due to their small group velocity [44,45].This behavior appears to persist in the nonintegrable model.A large fraction of photons in the bound state remain in the vicinity of their initial sites even after many cycles.In Fig. 8, we show n i (n c ) for an N =12 bound state, demonstrating this robust nonthermalizing behavior.We can quantify this by calculating n init = i∈initial sites n i , the occupation of all the sites initially occupied by a photon.The average of n init for different size bound states can be seen in the inset of Fig. 8. From this perspective, the bound states appear to grow more robust as they increase in size.This at first seems to be in contradiction with the BSP results from Fig. 7.However, the BSP measures the overlap with the N -photon bound state as a whole, while n init also captures the case when the bound state loses photons from the edges while its core stays robust and does not move away significantly from its initial position.This is precisely what happens for large-N bound states.Since we are considering hardcore bosons, the photons from the middle of the bound state can only hop to the decorations and back, while the photons at the edge can move further away along the chain and become detached from the rest.

C. Floquet modes
In order to understand the robustness of bound states to integrability breaking, we compare them with the Floquet modes of the system, i.e., the eigenstates of ÛF .In Figs.9(a)-(b) we plot the overlap of the all the eigenstates with the initial bound state for N =3 photons, contrasting the integrable case with that of θ = 0.2π.The color scale represents the number of pairs of neighboring occupied sites for each state.For example, the value of 2 corresponds to a three-photon bound state, 1 to two neighboring and one separate photon, while 0 implies three separate photons.The integrable case, θ = 0, in Fig. 9 contains the bound states (red points).The energies of these sectors overlap, but there is no mixing between the states.The overlapping energies are a consequence of the periodic Floquet spectrum.In the Hamiltonian XXZ model these sectors are separated by energy gaps.As θ is increased, the system becomes non-integrable and the sectors start to mix.At θ = 0.2π, Fig. 9(b), the bound states have almost merged with the bulk, but still remain visible.This is no longer the case after θ = 0.3π, which is consistent with the results of Figs.7(d)-(e) where it was shown that the bound states are robust only up to this point.However, the bound states survive at values of θ which are clearly in the non-integrable regime from the point of view of the level statistics, recall Fig. 2(a).
We have also repeated these calculations for a translation-invariant initial state in the k=0 momentum sector.The main difference is the reduction in number of bound Floquet modes.There are two such states in the k=0 sector and their energies and overlaps with the initial state show almost no dependence on the system size L.This is related to the fact that the BSP does not significantly decay with the system size.
The overlap plots for larger numbers of photons display similar features to N =3.At θ =0, there are N separate sectors which are defined by the number of pairs of adjacent photons.The bound Floquet modes slowly mix with the other sectors as the coupling to the decorations in increased.The case of N =5 and θ =0.1π is shown in Fig. 9(c).Here we see three prominent towers of states, which are related to the oscillations in the BSP in Fig. 7(c).These oscillations can be attributed to photons from the bound state hopping onto the extra sites and back.The towers appear as soon as θ =0 and persist until approximately θ ≈0.3π.The distance between the towers and therefore the oscillation frequency depends approximately linearly on θ .Moreover, the shape and height of the towers depend on the number of decorations attached to the initially occupied sites.
Intriguingly, we find that the towers are more prominent for the five-photons initial bound state with two decorations on the second and fourth site [Fig.9(c)] than for the state with three decorations on the first, third and fifth site.Upon closer inspection, some towers of highoverlap states can also be discerned for N =3 in Fig. 9(b).However, they are not as well differentiated as for N =5, and do not appear to be equally spaced in energy, which is the reason why the oscillations in BSP at θ =0.2π are irregular, see Fig. 7(a).Interestingly, the towers do not become better resolved with increasing the photon number and the case of N =5 actually features the sharpest towers and corresponding oscillations in BSP and various local observables, such as the number of photons in the extra sites.As a side note, similar looking towers of states are often found in systems that host "quantum many-body scars" [60][61][62], however it is not clear whether similar physics occurs in the present case.

V. FINITE DENSITY OF EXCITATIONS AND OTHER DECORATION PATTERNS
Up to this point, we have mostly focused on the experimental setup of Ref. [44], restricting to the case with integrability-breaking decorations on every other site and small numbers of photons N ≤6.In this Section, we consider other cases that were not studied previously.In particular, we will now fix the filling factor instead of fixing the total photon number, which will allow us to investigate convergence to the thermodynamic limit by growing both the system size and the number of excitations, as conventionally done in the literature.Moreover, we will explore other decoration patterns, including those with decorations on every n-th site where the system still retains translational invariance, as well as completely random patterns which break all the symmetries.

A. Fixed filling factor
When we extrapolate systems with small but fixed photon numbers to the infinite number of sites, they become infinitely dilute.By contrast, the thermodynamic limit is conventionally taken by keeping the density constant.Thus, we introduce the filling factor ν = N/L sites and study properties of our circuit as both N and L sites are simultaneously increased such that ν remains constant.
As done previously for fixed N , we average the BSP (8) over a certain number of cycles and investigate its dependence on θ .In particular, the averaging was done between cycles n c =100 to n c =200, after the initial drop in bound state survival.The time-averaged value is then extrapolated to the thermodynamic limit, while keeping ν fixed, using a quadratic fit in 1/L, see Fig. 10.Here we set the filling factor to ν = 1/10, but the results for other ν values are similar.In contrast to dilute systems with fixed photon numbers, the average BSP now drops to very small values already at θ ≈ 0.05π.Due to a very slow decay rate of the BSP for very weak integrability breaking θ 0.05π, we expect Fig. 10 to provide only an upper bound, as the extrapolated value of the BSP would likely be smaller with an access to a longer time window.However, the TEBD calculations become significantly more time consuming with increasing number of cycles, which is a limiting factor in very large systems.In the inset of Fig. 10 we show the evolution of BSP for N =100 photons in system size L=1000+500, and several values of θ obtained using TEBD.Unlike the case of fixed but small photon numbers, the large-N bound states are less resilient to integrability breaking by coupling to the extra sites.The BSP quickly decreases with the number of cycles, as can be seen in the inset, where we only show the relatively small coupling strengths θ ∈ [0, 0.1π].The extrapolated values in Fig. 10 point to the conclusion that very small photon number sectors we studied up to this point have unconventional properties, which are not shared by thermodynamically large sectors with non-zero filling factor ν.

B. Other decoration patterns
Finally, we consider the level statistics and BSP dynamics for different patterns of decorations.In Fig. 11 we show the average level spacing ratio r (θ ) for N =3 photons and five different types of decoration arrangements.The first one [Fig.11(a)] is only a single decoration, while the second one [Fig.11(b)] consists of three decorations attached to sites 2, 4 and 6.Both of these patterns break translation symmetry, which limits the system sizes we can reach.Unlike the level statistics in Fig. 2(a), where r was oscillatory for N =3, here we observe no such oscillations.However, r first increases to the COE value and then starts to decrease around θ ≈ 0.3π.The decrease does not happen for larger photon numbers, with the r (θ ) curve becoming flat already at N =4.Thus, we conclude that N =3 displays anomalous level statistics properties, irrespective of the decoration pattern.We have also examined the DOS for these other patterns of extra sites (not shown).The DOS distribution for the cases of one and three decorations differs from Fig. 3(a) in that it loses the sharp peaks, but the overall shape stays approximately the same and becomes flatter with increasing number of photons.The peaks likely disappear because the new decoration pattern is no longer translation-invariant, and the eigenstates which were previously degenerate are no longer related by symmetry, hence they generally have different energies.
Note that the r plateau is more pronounced and closer to the COE value for three decorations compared to a single decoration.This trend continues as we add more decorations.In Fig. 11(c) we show several random patterns where the number of extra sites is kept equal to half the number of sites inside the ring.We again observe similar behavior, with an initial plateau followed by a decrease in r .However, the plateau is now at the CUE value r CUE ≈ 0.597 instead of r COE ≈ 0.527.As explained in Sec.II A, this is due to all of the studied random patterns breaking inversion symmetry, unlike the previous cases of one and three decorations.It might seem surprising that r is not monotonic with system size in Fig. 11(c), but this is simply due to choosing completely different patterns for each system size and could be avoided by averaging over several random patterns for each L. The DOS distribution for the cases in Fig. 11(c) again has no peaks for θ =0 and is noticeably flatter than the previously considered patterns (data not shown).
We have also considered two examples of periodic patterns, one with decorations attached to every site of the main ring [Fig.11(d)] and the other with decorations on every third site [Fig.11(e)].The first case is invariant to translations by two sites and inversion which swaps the i-th site (decoration) with (L sites − i)-th [(L decor − i)-th], so these symmetries must be resolved in order to obtain the correct level statistics.We note that the full system is not inversion symmetric for the usual periodic pattern with decorations on every other site, even though the arrangement of decorations itself is.This is due to first applying the fSim gates on odd bonds and then on even bonds, Eq. ( 2).There is no reflection axis which simultaneously preserves both the decoration pattern and the order of even and odd fSim gate layers.Similar to previous results in Fig. 2(a), for N =3 we again observe deviations from r COE at certain values of θ .However, the r (θ ) curve is now symmetric around θ = π/2 with an integrable point in the middle.This is similar to the case in Fig. 14(a) in Appendix B, where r (θ ) is symmetric around θ=π/2.As can be seen from Eq. ( 1), this value of the hopping amplitude corresponds to a photon moving to the neighboring site with probability 1, so it is not surprising that this is a special case.
For the second periodic pattern [Fig.11(e)], the symmetries of the full system are translation by six sites and inversion which preserves this arrangement of decorations.In this case, we also observe oscillations in r (θ ), but the local minima and maxima are at different values of θ compared to the other patterns.As before, all oscillations disappear for N =4 or more photons.For all studied periodic patterns with decorations on every n-th site, the DOS still exhibits pronounced peaks for N =3 and to some extent for N =4.
We have also investigated the robustness of bound states for various decoration patterns.In Fig. 12 we plot the averaged BSP for N =3 photons, L sites =300 sites on the main chain averaged over n c =100 cycles between cycles n c =50 and n c =150.For a decoration on every second site we average over the two possible N =3 photon initial states around the center of the chain.We perform a similar averaging over the three possible initial states for the case of a decoration on every third site.For randomly allocated decorations we instead average over 5 different random patterns with the initial bound state at the center of the chain.These results suggest the robustness of the bound states is more dependent on the density of decorations L decor /L sites than on the ac-tual pattern.The bound states survive for larger values of θ when the number of decorations is smaller.The average BSP decays over a similar range of θ for random patterns with L decor /L sites =1/2 and for the periodic case of the same site density.For the case of decorations on every site there is a peak at θ =π/2 which corresponds to the near integrable point as seen in Fig. 11(d).As before, the bound states become increasingly less robust as the number of photons grows.

VI. CONCLUSIONS AND DISCUSSION
We have performed systematic classical simulations of the Floquet XXZ circuit on a 1D chain with integrability breaking decorations.This study was motivated by the recent Google experiment [44] which realized the same model on a ring of superconducting qubits and investigated the dynamics of its bound states.Surprisingly, the bound states were observed to be resilient to integrability breaking perturbations in the form of extra qubits attached to the ring.We have analyzed the level statistics of this model and simulated the dynamics of bound states, confirming that some of these states indeed survive in certain parts of the non-integrable regime, even for an infinite number of qubits and at infinite time.In contrast to much previous work, the focus of Ref. [44] and our own was on dilute systems containing few excitations, which have rarely been discussed in the literature (see, however, Ref. [63]).As we have demonstrated, such systems are amenable to classical simulations in large numbers of qubits, providing useful benchmarks for future studies on improved quantum hardware.
One of our most significant findings is that small but fixed photon number sectors show unusual properties in several respects.In particular, the robustness of bound states depends on the photon number, with larger states decaying more rapidly as the coupling to the integrability-breaking extra sites is increased.Moreover, the energy spectrum for N =3 photons has unusual level statistics, which deviates from the expectation for a chaotic system even for strong couplings to the integrability-breaking decorations.This was attributed to the presence of special eigenstates in the energy spectrum.These eigenstates have a relatively simple structure, which is related to one-photon and twophoton states, while the rest of the photons are located in the decorations.The proportion of such states is large enough only in sufficiently dilute systems.When the decoration pattern is periodic, some of these eigenstates are related by translation and are therefore degenerate in energy, which results in prominent peaks in the DOS and affects the level statistics.The deviations in level statistics were shown to leave an imprint in the dynamics of bound states by slowing down the thermalization.
Additionally, we have investigated systems with constant filling factors and their extrapolation to the thermodynamic limit.Such systems are no longer dilute and our findings indicate that they do not support stable bound states when integrability is broken.This is in stark contrast with the bound states in very dilute systems with small photon numbers, such as the one studied in experiment [44].Moreover, we have explored other decoration patterns, including both periodic and non-periodic ones.A brief summary of all considered systems is given in Table I.Our calculations suggest that the peaks in the density of states disappear when the pattern is not periodic, which destroys the translation symmetry of the full system.This is likely a consequence of certain eigenstates no longer being degenerate.We also find that the inversion symmetry of the decoration patterns (or lack of it) influences the level statistics.In particular, inversionsymmetric patterns are consistent with COE and nonsymmetric with CUE statistics in the chaotic regime.Deviations from these values were observed only for N =3 photons and were found to diminish as the number of photons is increased.However, our results do not indicate a link between the irregularities in level statistics and the robustness of bound states, although both properties are most prominent in dilute systems.For example, nonperiodic decoration patterns result in level spacing ratios consistent with random matrix theory, implying that the integrability is indeed fully broken, while the few-photon bound states remain robust in that regime.
One advantage of classical simulations performed in this work is direct access to the system's properties at finite energy densities.Thus, the model considered in this work would be a useful platform for benchmarking quantum algorithms that target states at finite energy density [64].Moreover, it would be interesting to explore other models which host bound states, e.g., the chiral Hubbard model [45], and investigate if such models exhibit similar behavior in relation to the density of excitations and integrability breaking by changing the geometry of the system, as described in this work.
Note added: After the completion of this work, we became aware of Ref. [65] which also studied the stability of bound eigenstates in the special case of N =3 photons and decorations on every second qubit.Based on perturbative arguments and the scaling of inverse participation ratio, Ref. [65] concluded that N =3 eigenstates slowly lose their bound state character in the L→∞ limit.Our finite-size scaling analysis above suggests that the bound state probability remains finite in this limit for N =3, however this cannot rule out the possibility of a much larger length scale, at which all dynamical signatures of bound states would ultimately disappear at infinite time.

Figure 1 .
Figure 1.(a) Sketch of the XXZ circuit model with L=14+7 sites (7 unit cells).Filled dots denote a bound state of N =3 photons.Here, the integrability-breaking decorations are attached to every other site.(b) An example of the corresponding quantum circuit with L=4+2 sites.The circuit consists of fSim (boxes) and SWAP gates (vertical lines).The alternating layers of gates acting on even/odd bonds and decorations are denoted by blue, red and green color, respectively, matching the unitaries Ueven, U odd and U dec in (a).Our classical MPS simulations in iTensor[46] follow this diagram and assume open boundary conditions.

Figure 2 .
Figure 2. Statistics of the Floquet quasienergies.Average ratio of consecutive energy gaps r for different values of θ with fixed θ = π/6, φ = 2π/3, β = 0.The horizontal dashed lines are r P ≈ 0.386 and r COE ≈ 0.527.(a) N = 3 photons for different system sizes indicated in the legend.The relevant Hilbert space dimensions range from dimL=20+10 = 406 to dimL=80+40 = 7021.The vertical lines A, B, C, D and E mark several values of θ (0.10, 0.35, 0.65, 0.75 and 0.95) that will be studied later in more detail.Inset: linear extrapolation to L→∞.(b) The case of N = 4 photons and (inset) N = 5 photons, with the largest Hilbert space dimensions dimL=40+20 = 24405 and dimL=20+10 = 14253, respectively.In all the plots, we resolve the translation symmetry and consider only the k=0 momentum sector.

5 Figure 3 .
Figure 3. Density of states (DOS) for different values of θ , normalized by the average over the entire energy spectrum [the labels A-E are defined in Fig. 2(a), while 0 denotes the integrable case θ = 0].The main panel corresponds to N = 3 photons in a system size L = 80 + 40, while the inset shows N = 5 in L = 20 + 10.In both cases, θ = π/6, φ = 2π/3.

Figure 4 .
Figure 4. (a) Comparison of the actual dispersion of twophoton states for the integrable case θ =0 with no extra sites (dots) and the theoretical prediction for two separate noninteracting photons (crosses) and two bound photons (dashed lines).The color scale corresponds to the deviation between each dot and the closest cross.(b) Same as (a) but with added extra sites, while θ =0.

Figure 5 .
Figure 5. (a) Bands of special states constructed from numerically-obtained single-photon bands for θ ∈ [0, π].Here we consider the case of N = 3 separate photons.The bands are denoted by the number of photons in the first (1), second (2) and third (3) single-photon band.(b) Corresponding DOS.The bright regions can be related to the peaks in Fig. 3.

Figure 6 .
Figure 6.(a) Comparison of the level statistics ratio r and the Thouless time extracted from the spectral form factor (SFF) for N = 3 photons, L = 60 + 30.Inset: SFF time series corresponding to the main plot for several values of θ [the labels A-E are defined in Fig. 2(a)].The data was averaged over β and smoothened by a moving average (see text).Note the logarithmic scale on both the x and y axis.The horizontal dashed lines mark the saturation values, H and H + n 2 0 , while the vertical line at nc = H is the Heisenberg time.The dashed black curve is the COE prediction for the linear ramp.The agreement with COE becomes better as N increases.(b) Same for N = 5 and L = 14 + 7, which shows a much clearer linear ramp in the SFF and better agreement with the random matrix theory at late times.

Figure 8 .
Figure 8. Site occupation of a N =12 photon bound state after some number of cycles indicated on the legend.Data is obtained by TEBD for at system size L = 300 + 150, bond dimension χ = 320 for θ = π/6, φ = 2π/3 and θ = π/4.Inset: ninit/N for bound states of sizes N = 3 to N = 12 at different values of θ .Data is obtained by TEBD for 150 cycles and the same parameters as in the main plot.

Figure 9 .
Figure 9. Overlap of the initial state with the eigenstates of the one-cycle evolution operator.(a)-(b): N = 3 photons in system size L = 30 + 15 for θ = 0 and θ = 0.2π, respectively.(c) N = 5 photons in system size L = 14 + 7 with θ = 0.1π.The color scale is the number of pairs of adjacent occupied sites.

Figure 10 .
Figure 10.Time-averaged BSP extrapolated to infinite system size.Filling factor is fixed to ν=N/Lsites=1/10.For comparison, we also replotted the case of fixed N = 3 from Fig.7.Inset: BSP dynamics for L=1000 + 500 and N =100.All data is obtained by TEBD with bond dimension χ=256.

3 Figure 12 .
Figure 12.Time averaged BSP over 100 cycles for different decoration patterns in the N = 3 excitation sector for Lsites = 300.The parameters given are θ = π/6, φ = 2π/3 and χ = 256.In the case of decorations on every second (third) site the results were averaged over two (three) non-equivalent initial bound-state configurations.The random pattern results were averaged over five different patterns with decorations on half of the sites, L decor /Lsites = 1/2.