Robustness and eventual slow decay of bound states of interacting microwave photons in the Google Quantum AI experiment

Integrable models are characterized by the existence of stable excitations that can propagate indefinitely without decaying. This includes multi-magnon bound states in the celebrated XXZ spin chain model and its integrable Floquet counterpart. A recent Google Quantum AI experiment [A. Morvan et al., Nature 612, 240 (2022)] realizing the Floquet model demonstrated the persistence of such collective excitations even when the integrability is broken: this observation is at odds with the expectation of ergodic dynamics in generic non-integrable systems. We here study the spectrum of the model realized in the experiment using exact diagonalization and physical arguments. We find that isolated bands corresponding to the descendants of the exact bound states of the integrable model are clearly observable in the spectrum for a large range of system sizes. However, our numerical analysis of the localization properties of the eigenstates suggests that the bound states become unstable in the thermodynamic limit. A perturbative estimate of the decay rate agrees with the prediction of an eventual instability for large system sizes.


I. INTRODUCTION
In recent years, quantum simulators have become a powerful tool to investigate the non-equilibrium dynamics of quantum many-body systems.Experiments based on many different platforms have now the capability to prepare a quantum state with good fidelity and evolve it almost unitarily, preserving its coherence for sufficiently long times to observe its interacting dynamics [1][2][3][4][5][6].These advances sparked a considerable interest in the theoretical investigation of the approach to thermal equilibrium in isolated quantum many-body systems [7][8][9][10][11][12][13]. At the same time, the available quantum simulators inspired the search for "unusual" behaviors, where a system does not reach thermal equilibrium.A paradigmatic example is the phenomenon known as quantum many-body scarring, which was first discovered in an experiment with Rydberg atom arrays [14]: it was found that certain quantum states exhibit long-lived non-ergodic dynamics, while other states instead show fast thermalization.Following the initial discovery, numerous theoretical studies have been conducted to elucidate the conditions under which quantum many-body scars occur and to develop a theoretical framework that can explain and predict such phenomena [15][16][17][18][19][20][21].Other examples of exotic dynamical phenomena that have been observed in quantum simulators include, for instance, Hilbert space fragmentation [22][23][24], dynamical phase transitions [25][26][27], time crystals [28][29][30][31][32][33][34], and noise-resilient edge modes [35][36][37][38][39][40][41].
Another puzzle for our understanding of quantum many-body dynamics was recently observed in a Google Quantum AI (GQAI) experiment [42,43] based on superconducting circuits hosting microwave photons.The dynamics of the photons, which can hop between neighboring sites of a chain shown in Fig. 1(a) and interact with each other, is described by an integrable quantum circuit.The experiment showed the presence of bound states of up to 5 photons: these are exact eigenstates of the model 1d chain with extra sites forming a comb in the GQAI experiment, with the comb teeth connected to every other site of the original chain; a unit cell and labelling of sites in each cell is shown.(c) Schematic of distinct two-particle bound states, which become dimers in the very strong binding limit, in a unit cell on the comb; these give rise to three bands in the momentum space.(d) Distinct three-particle bound states, schematized as trimers, in a unit cell on the comb, producing four bands in the momentum space.
predicted by Bethe ansatz [44], where the photons are nearly adjacent and form a single collective excitation.The bound states observed in the study are enclosed in the continuous spectrum of multi-particle states, i.e., eigenstates of (mostly) distant photons that can scatter off each other.Even in the absence of a gap, the conservation laws of the integrable model protect the bound states from mixing with the underlying continuum.The work [42] demonstrated a remarkable agreement of the experimental results with the analytical solution.
The integrable circuit was then perturbed by coupling every other site of the chain with additional sites as shown in Fig. 1(b), thus breaking the integrability of the model.In the absence of any conservation laws that protect them, the exact bound states of the unperturbed model are expected to quickly decay in the continuum of multi-particle states.Nevertheless, the experiment showed the persistence of the bound states after many cycles of the circuits, even for fairly large values of the perturbation.
The origin of this robustness cannot be easily attributed to any known mechanism for slow relaxation.Thus, we do not expect usual quantum many-body scarring mechanisms to be operative in our systems with O(1) integrability breaking, and there is no evidence of even more dramatic Hilbert space fragmentation phenomena.Among other possibilities, the rigorous theory of prethermalization [37] applies only for perturbing much more special initial Hamiltonians and does not apply for perturbing generic integrable systems like in our case (except in the regime with dominant on-site terms).One could contemplate a recently discussed mechanism of weak integrability breaking of integrable models, where special perturbations effectively break the integrability only at higher order [45][46][47].However, our experience with this mechanism suggests that it would require more complicated special perturbations than in the experiment (and it would approximately preserve features in the entire spectrum and not just some specific states).While we cannot rule out such interesting and unusual explanations, we are led to consider perhaps simpler mechanisms having to do with particular quantitative properties of the specific system.Thus, an early example in Ref. [48] of unusual thermalization in a non-integrable model for some initial states was explained in Ref. [49] by proximity of these states to the boundaries (ground or ceiling states) of the many-body spectrum; that is, the dynamics was governed by being in a "special corner" of the manybody Hilbert space away from ETH-like physics expected in the middle of the many-body spectrum.
While for Floquet systems there is no notion of ground or ceiling states and hence no such proximity can be invoked, we note that the circuit and setup that models the experiment in [42] is not a truly many-body system.The puzzling robust bound states were observed for a system with a fixed number of photons (N = 3), and the Hilbert space dimension scales only polynomially with the system size, not exponentially.It is reasonable to expect that in such few-body models the decay is not as fast as for typical eigenstates in the middle of a truly many-body spectrum.Therefore, in order to understand and characterize the robustness of the bound states, it is important to perform a quantitative study for increasing system sizes, and to contrast the persistence of the bound states against the scaling expected for similar few-body models.
In this work, we examine the spectral properties of the perturbed model in this few-body regime using exact diagonalization and physical arguments.We compare the spectrum of the circuit with the one of the related Hamil-tonian system under corresponding continuous time dynamics.First, we find that in the closely related Hamiltonian case the N = 3 bound states are in fact protected by a gap and persist in the thermodynamic limit.This is true also for O(1) integrability breaking perturbations for (a range of) serendipitous choices of parameters, including the Hamiltonian parameters reasonably related to the ones in the experiment.On the other hand, in the case of the Floquet circuit the bound states are "folded" and overlap with the multi-particle continuum in the corresponding Floquet spectrum, but for the values of parameters studied in the experiment and N = 3 this folding affects only the boundaries of the Hamiltonian spectrum.Because of this and for other quantitative reasons, despite the folding, we clearly detect isolated bands in the Floquet spectrum that correspond to the unperturbed bound states, even for large values of the perturbation.Our numerical analysis suggests, however, that these bound states eventually become unstable in the thermodynamic limit.To argue this, we study the inverse participation ratio in the center of mass frame, whose scaling with system size can distinguish a bound state from a "scattering" (i.e., unbound) state.We show that the inverse participation ratio of the bound states decreases with the system size, indicating a decay in the thermodynamic limit.However, the decay is very slow: we attribute the slowness to the impossibility in our parameter regime for a 3-photon bound state to resonantly decay into a scattering state of a 2-photon bound state and a single photon; the 3-photon bound state can only decay into a scattering state of 3 isolated photons, a process with a very small (but non-zero) matrix element.Our conclusions are further corroborated by a perturbative analysis of the matrix elements and density of states, which shows a small but finite decay rate in the thermodynamic limit.
The paper is structured as follows.In Sec.II we introduce the integrable models (Hamiltonian and Floquet circuit) and the non-integrable perturbation.In Sec.III we analyse the spectra of the Hamiltonian and the Floquet circuit in the sectors with N = 1, 2, and 3 photons, for different values of the perturbation strength.In Sec.IV we study the properties of the 3-particle bound states in the Floquet circuit for increasing system sizes, and show their eventual instability in the thermodynamic limit.In Sec.V we study the matrix elements and density of states in the unperturbed model, in order to explain the instability and derive a perturbative estimate of the decay rate.We conclude in Sec.VI with suggestions for more experimental tests and broader prospects.

II. HAMILTONIAN AND FLOQUET MODELS
We will focus on a qubit model that evolves under (i) a static Hamiltonian or (ii) a Floquet circuit.In both cases, the dynamics we consider conserves the total number of qubits in the state "1" (we will refer to them as "particles") and can be expressed in terms of the two-qubit operator h j,k (w, u) acting on neighbouring sites: where w is the hopping amplitude between sites j and k and u is an interaction (attraction energy for u > 0) when particles occupy these sites.
The Hamiltonian and the Floquet circuit that we will consider are perturbations of an integrable Hamiltonian/Floquet model, respectively.We will first introduce the integrable models (Sec.II A and II B), and then consider the non-integrable perturbation (Sec.II C).

A. Hamiltonian: integrable case
The integrable chain Hamiltonian with periodic boundary conditions (assumed throughout) is simply ( The Hamiltonian H 0 can be equivalently written in terms of spin 1/2 operators, defined as which is the XXZ model in the uniform magnetic field.

B. Floquet circuit: integrable case
An integrable circuit, also known as the Floquet XXZ model, can be defined using the following two qubit gates: where fSim is the two-site unitary gate using the same notation as in the GQAI experiment [42].The gates are applied to even-odd and odd-even pairs of sites in a brickwall pattern.The unitary operator that describes the evolution over a single cycle is defined as The integrability of the model was first demonstrated in [50], and a Bethe ansatz solution was obtained in [44].The bound states of the model, that were analytically studied in [44], were then detected in the GQAI experiment [42].

C. Non-integrable perturbation
In Ref. [42] the integrable model was perturbed by adding sites as in Fig. 1(b).The additional sites are connected in the shape of "comb teeth" to every other site of the original chain.To label the sites in the new geometry, we use a composite index (R, α) where R labels the unit cell and α ∈ {1, 2, 1 ′ } labels the three sites of each unit cell.In this new geometry, we can write perturbed models for both the Hamiltonian and the Floquet case.
The Hamiltonian is given by The Floquet circuit is similarly obtained by applying the two-qubit gates on the three sets of pairs: When w ′ = 0 (θ ′ = 0 in the Floquet case), the particle number on the original chain is conserved (denoted N 1∪2 ), and the particle number on each (R, 1 ′ ) site is conserved.In what follows, we often use more crude grouping of states labelled by sectors In this case, in the sector with no particles on the 1 ′ sites, the u ′ term does not operate at all and the Hamiltonian (Floquet circuit) is equivalent to the integrable chain in Eq. (2) [Eq.( 5)] for any u ′ (ϕ ′ ).

III. SPECTRUM COMPARATIVE STUDY OF THE HAMILTONIAN AND FLOQUET SYSTEMS
In this section we present the spectra of the Hamiltonian (8) and of the Floquet circuit (9) for comparable values of parameters.In particular, for the Floquet circuit, we use the same parameters used in the GQAI experiment, namely θ = π/6 and ϕ = ϕ ′ = 2π/3.We expect this model to resemble the continuous Hamiltonian evolution for a ratio of parameters u/w = ϕ/θ = 4 and u ′ /w = ϕ ′ /θ = 4.We will first discuss the general features of the spectrum in the system with N = 1, 2, 3 particles.While the spectrum for N = 3 is rather complicated, we will nevertheless show some qualitative, distinctive features through a comparison between the Floquet and the Hamiltonian cases.We will then focus on the N = 3 bound states, by examining a narrower window of the spectrum using some observables that can specifically signal the presence of bound states.The spectra as a function of momentum k are plotted in Fig. 2 in a sector with fixed number of particles N = 1 and N = 2: for the Hamiltonian H in Eq. ( 8) we plot the spectrum E(k) (rescaled with a factor −π/6 for comparison); for the Floquet circuit we plot the quasienergies ϵ(k), defined as the complex phases of the eigenvalues of the Floquet operator in Eq. (9).The spectra are computed for increasing values of the ratio w ′ /w = θ ′ /θ.For w ′ /w = θ ′ /θ = 0, the particles on the 1 ′ sites cannot hop, so the particle number on the original chain is conserved, and the particle number on each (R, 1 ′ ) site is conserved.In this case, in the sector with no particles on the 1 ′ sites, the u ′ (ϕ ′ ) term vanishes (acts trivially) and the Hamiltonian (Floquet circuit) is equivalent to the original integrable chain for any u ′ (ϕ ′ ).The eigenstates in this sector are represented in red, and are analyzed in more detail in App.A for the Floquet case.
For w ′ /w = θ ′ /θ ̸ = 0 the states of the integrable chain hybridize with the states having non-zero occupation of the 1 ′ sites.As a result, gaps open at k = ±π in the single particle spectrum (N = 1).A rearrangement of the spectrum is observed in the two-particle sector (N = 2): as the bands in the single-particle spectrum become flatter, some gaps open in the two-particle continuum and some isolated states appear in the gaps.The bound states (three isolated bands at the top, corresponding to the three dimer configurations in Fig. 1(c)) are still observable in both the Hamiltonian and Floquet case as w ′ /w = θ ′ /θ is increased from 0 to 1: while they overlap with the two-particle continuum for some values of momentum, for other values they are protected by a gap.
The spectra for N = 1 and N = 2 have very similar features in the Hamiltonian and Floquet cases.The most notable difference is the breaking of time-reversal symmetry in the Floquet case which makes the quasienergy spectrum asymmetric for k → −k.When we instead compare the Hamiltonian and the Floquet spectra in the threeparticle sector (N = 3, Fig. 3) we observe a substantial difference: since the quasienergies are defined modulo 2π, the Floquet spectrum corresponds to a "folded" Hamiltonian spectrum.As a consequence, the 3-particle bound states, which are gapped and thus stable in the Hamiltonian case, are folded and overlap with the continuum of the Floquet spectrum and therefore they are not protected by a gap.The mixing of the bound states with the continuum in the Floquet case can lead to the decay of the bound states.Nevertheless, quantitatively this mixing can be still fairly weak, and the bound states may be visible in the spectrum even for fairly large system sizes.
These bound states are difficult to identify in Fig. 3  [Note that when w ′ /w = θ ′ /θ = 0 all states in the sector (0, 3) have zero energy and their color is RGB (0,0,0), which is black.]The system size we choose is rather small (Luc = 12 unit cells) to avoid overwhelming the plots.

in
the Floquet case, where they overlap with the continuum.To characterize these states and to discern them from the continuum, it is useful to consider quantities that are sensitive to the relative configuration of the three particles.Examples of such observables are shown in Fig. 4: for each eigenstate |ψ i,k ⟩ with momentum k, we compute the probabilities of the following configurations for the three particles in neighboring sites [trimers in Fig. 1(d)]: where |. ..⟩ k represents the (normalized) projection of the ". . ." state in the sector with momentum k, i.e., . . . . . .
An equivalent definition for P T I (and analogously for the other trimers) is (18) In Fig. 4 we plot the probabilities P T I + P T II and P T III + P T IV for each eigenstate.For w ′ /w = θ ′ /θ = 0, we observe two (dispersive) bands with large P T I + P T II , which correspond to the exact three-particle bound states of the integrable chain (see App. A).Two (flat) bands have large P T III + P T IV : they can be interpreted as bound states of the two particles on the chain localized in the potential of the particle on the extra sites (which acts as an immobile impurity), see App.B for detailed study in the Hamiltonian case.As we turn on the hopping along the teeth (w ′ /w = θ ′ /θ ̸ = 0), the bands with large P T III + P T IV acquire a rather weak dispersion, signaling a very low mobility of the trimers with one particle on the extra sites.
In the cases w ′ /w = θ ′ /θ = 0.1, 0.5, 1.0 the four bound states are still characterized by large values of  15), ( 16)].The system size is Luc = 36 unit cells.The system parameters are the same as in Fig. 3.A narrow window of quasienergies is shown focusing on the three-particle bound states, revealing their character of being primarily chain trimers or chain-extra site trimers and also where significant mixing is present.In the Hamiltonian system the three-particle bound states are isolated from the continuum and persist in the thermodynamic limit even for w ′ /w = 1, while in the Floquet case they are inside a continuum and will decay in the thermodynamic limit but apparently survive to fairly large sizes.The similarity between the bound states in the Hamiltonian and Floquet systems is notable, allowing to infer properties of the latter ones as well as of the surrounding continuum from the more simple Hamiltonian understanding.
both P T I + P T II and P T III + P T IV : this shows that even in the Floquet circuit, where the bound states are not protected by a gap, the hybridization is strong among the four bound states but quite weak between the bound states and the continuum.Even for the largest value we consider (θ ′ /θ = 1.0), the hybridization with states in the continuum is clearly visible only for one of the four states (the one with quasienergy ϵ ∼ 4).
Another useful quantity to characterize the localization properties of the particles in the bound states is the appropriate inverse participation ratio.This quantity has been extensively used as a signature of localization induced by disorder both for single particle models [51,52] and for many-body systems [53][54][55].Here, instead, we will use it to probe the localization of particles with respect to their center of mass: this localization is induced by interactions and can lead to the presence of bound states.
The idea is to find a quantity that is sensitive also to other "bound" configurations of the three particles, that are not captured by the tightest binding configurations defining P T I , P T II , P T III , P T IV .For example, a bound state can have large overlaps with some other configurations such as |. . . . ..⟩ k .However, we still expect it to have support on a small (compared to the dimension of the sector) number of such configurations.
To this end, it is useful to define the computational basis in the sector with momentum k.Similarly to the definition of Eq. ( 17), each state in the computational basis is defined from a classical configuration c (called representative) of the three particles as where T is the translation operator by one unit cell and M c ≥ 0 is the smallest positive integer such that . The state |c⟩ k can be also viewed as a normalized projection of the state |c⟩ into the sector with momentum k.The basis is obtained by taking a single representative c for each class of configurations that are related by translations.We define the momentum-resolved inverse participation ratio (IPR) for a normalized eigenstate |ψ i,k ⟩ with momentum k as where the sum runs over all the representatives (i.e., distinct classes) that define the basis.This quantity is an indicator of the localization of the three particles in the frame of their center of mass [57].For a bound state, this quantity converges to a finite value in the large L uc limit, while it scales as 1/L uc for scattering states of a two-particle bound state with a single particle, and as 1/L 2 uc for scattering states of three unbounded particles.As shown in Fig. 5, the inverse participation ratio I k takes large values for four distinct bands of bound states: these are the same states that were characterized by large values of P T I + P T II and P T III + P T IV in Fig. 4. Figure 5 shows that, while the value of I k remains large for the Hamiltonian case even up to w ′ /w = 1, the bound states in the Floquet spectrum exhibit a clear decrease of I k as θ ′ /θ goes from 0.0 to 1.0.This suggests that perturbing the integrable Floquet circuit through the inclusion of additional sites tends to unbind the original bound states of the model.

IV. EVENTUAL INSTABILITY OF THE 3-PARTICLE BOUND STATES IN THE FLOQUET MODEL
The decrease of the inverse participation ratio as a function of w ′ shown in the previous section, Fig. 5 at fixed L uc = 36, suggests that the localization length in the center of mass frame tends to grow with the perturbation strength.If this localization length diverges in the thermodynamic limit L uc → ∞, the 3-particle bound state is unstable.
In order to probe the eventual decay of the 3-particle bound state, we examine the scaling of the inverse participation ratio with the system size, for θ ′ = θ, where the perturbation is of the same order as the unperturbed parameters and the bound states are expected to be least robust of the parameters θ ′ ≤ θ considered in the previous sections.In Fig. 6 (upper panel), we plot the inverse participation ratio of the energy eigenstates in the k = 0 sector (computed as described in Sec.III) for different numbers of unit cells L uc .Despite the fairly large perturbation strength, the inverse participation ratio I k=0 still shows four clearly visible peaks that correspond to the four bound states of the unperturbed model.For large L uc , however, it is possible to notice the effect of the perturbation, which mixes the bound states with the underlying continuum and smears out the peaks.To quantify this effect, in Fig. 6 (lower panel) we plot the height of the peak (defined as the maximum of I k=0 in an appropriate energy window) as a function of L uc for the four bound states.The fluctuations with the system size are still very large, indicating that the results are still very sensitive to the finiteness of the level spacings in the spectrum and to their statistical fluctuations.Nevertheless, all the data shows a decreasing trend with system size, for all the four bound states.These results suggest that the bound states will ultimately decay in the thermodynamic limit, and that the decay is very slow, leading to persistent 3-particle bound states for numerically and experimentally accessible system sizes.
Of the four bound states, we observe that one of them (the one with quasienergy ϵ ∼ 4.0032) exhibits significantly faster decay of I k=0 with the system size (Fig. 6, lower panels).We will now argue using numerical experiments that this faster decay is caused by the proximity in the spectrum with scattering states of a two-particle bound state with a single particle (which we will refer to as "2+1 continuum" below).Note that strictly speaking this label refers to states present in the integrable model at ϕ ′ = 0 in the sector (N 1∪2 , N 1 ′ ) = (3, 0), while the states highlighted in Fig. 7 at ϕ ′ = ϕ are their descendants.As shown in Fig. 7(a), the 2+1 continuum is characterized by a larger probability of trimer configurations and larger I k compared to the scattering states of individual particles (which we will call "1+1+1 continuum").The edge of the 2+1 continuum is very close to the lowest bound state, suggesting that these states are  responsible for the larger hybridization.In Fig. 7(b) we also consider a different choice of parameters, for which the spectra of the two lowest bound states are partially enclosed in the 2+1 continuum.In this case, these two states show a fast decay of I k=0 with the system size L uc , while the other bound states are more resilient (Fig. 8).
These results suggest that the decay is fast when the 3-particle bound state can decay in the 2+1 continuum, but is very slow when it can only decay in the 1+1+1 continuum.This explains the apparent robustness observed in the experiment [42].

V. PERTURBATIVE ANALYSIS OF EVENTUAL INSTABILITY
To address the question of an eventual instability, it is useful to study the matrix elements of the perturbation (i.e., the hopping along the teeth) in the basis of the unperturbed eigenstates of the Floquet model with ϕ ′ = ϕ, θ ′ = 0.The perturbation of the Floquet circuit is where θ ′ ≪ 1 is the small perturbative parameter, and V is defined as It is adequate for our purposes to use intuition from the Hermitean perturbation theory treating the Floquet quasienergies as the unperturbed energies and θ ′ • V as the perturbation.
We label the four bound states of the θ ′ = 0 model in the sector with total momentum k = 0 as |ψ n ⟩ with n = 0, 1, 2, 3: |ψ 0 ⟩, |ψ 2 ⟩ have quasienergies ϵ 0 ≈ 4.0512 and ϵ 2 ≈ 4.3614, respectively, and belong to the sector with (N 1∪2 , N 1 ′ ) = (2, 1); |ψ 1 ⟩, |ψ 3 ⟩ have quasienergies ϵ 1 ≈ 4.2657 and ϵ 3 ≈ 4.4755, and belong to the sector with (N 1∪2 , N 1 ′ ) = (3, 0).We generally expect a bound state |ψ n ⟩ to be unstable to a perturbation V if the "Fermi's golden rule rate" , where V nj ≡ ⟨ψ n |V |ϵ j ⟩ is the matrix element connecting the n-th bound state with the state |ϵ j ⟩ in the continuum and E n is the energy of the bound state.From the numerical study we know that the 3-particle bound states of interest to us do not overlap in energy with the 2+1 states (i.e., scattering states of a 2-particle bound state and a particle).Then conservation of energy and momentum implies that the bound states can only decay in the 1+1+1 continuum for k = 0.For these continuum states, from a simple counting argument we expect a density of states ∝ L 2 uc .The matrix element V nj = ⟨ψ n |V |ϵ j ⟩ between a state in the 1+1+1 continuum and a bound state (i.e., a localized state in the center of mass frame) can similarly be estimated from a simple argument: the state |ϵ j ⟩ of the three particles can be approximated as a (properly symmetrized) product of three plane waves, while |ψ n ⟩ is a single plane wave (with k = 0); the matrix element is non-zero only when the three particles are next to each other; taking into account the normalizations of the plane waves, we get that the matrix element scales as We then expect the product between the density of states and |V nj | 2 to yield an O(1) rate in the thermodynamic limit.
Note that a similar argument for a decay into the 2+1 continuum would give a density of states ∝ L uc and a matrix element , resulting, again, in a finite rate in the thermodynamic limit.In Figs. 9 and 10, we numerically check our prediction for the scaling of the matrix elements and of the density of states for the decay into the 1+1+1 continuum.We plot the matrix elements |V nj | (including only the ones that are non-zero) multiplied by L uc .As a measure of the density of states, we consider pairs ϵ j , ϵ j+1 of nearby levels in the quasienergy spectrum (including only the states with non zero |V nj |) and we plot the inverse spacing (∆ϵ) −1 = (ϵ j+1 − ϵ j ) −1 multiplied by L −2 uc as a function of the average quasienergy ϵ = (ϵ j+1 +ϵ j )/2.We find that, for all bound states n = 0, 1, 2, 3, both |V nj | • L uc and (L 2 uc ∆ϵ) −1 show a good data collapse, with no systematic dependence on system size [58], confirming our predictions that |V nj | ∝ L −1 uc and (∆ϵ) −1 ∝ L 2 uc .The scaling is further analysed in Fig. 11(a) by plotting the average of |V nj | 2 for each n = 0, 1, 2, 3 over the energy windows plotted in Figs. 10 and 9, as a function of system size.The results of the fits show a dependence |V nj | 2 ∝ Figure 10.Similar analysis as in Fig. 9, but for the bound states n = 0 (left) and n = 2 (right), which belong to the sector with (N1∪2, N 1 ′ ) = (2, 1).In the case of degeneracies, we consider a single |ϵj⟩, proportional to the projection of V |ψn⟩ on the degenerate subspace.The data collapse shows that both the matrix elements (top) and the density of states (bottom) obey the predicted scaling with the system size for 1+1+1 continuum.

L α
uc with α in the range −2.03, −1.49, in rough agreement with the expect scaling |V nj | 2 ∼ L −2 uc for matrix elements to the 1+1+1 continuum states.In Fig. 11(b) we plot the density of states ρ n (computed as the number of states in the same energy windows, divided by the width of the windows) as a function of L uc : the results agree with the expected scaling ρ n ∝ L 2 uc .The rate of decay is then computed as and Γn is plotted in Fig. 11(c).The rate is approximately constant, with no evident dependence on system size (except for a mild increase, which is visibly present only for small L uc ).However, the values of the rates Γn for different bound states span many orders of magnitudes.For n = 1, 2, 3 they are smaller than 10 −5 , suggesting that the bound states may persist at finite size for large values of the perturbation (or equivalently, for large system size with a fixed perturbation strength).[The corresponding physical decay rates Γ 1,2,3 < 10 −5 × (θ ′ ) 2 , for smallish θ ′ ≤ π/6, give lifetimes that can exceed 10 5 Floquet cycles.]The n = 0 has the largest decay rate (Γ 0 ≈ 4 • 10 −3 × (θ ′ ) 2 , significantly larger than the other three), and hence it should exhibit a more enhanced decay with system size for θ ′ ̸ = 0, in agreement with the results of Fig. 6.We remark, however, that the four bound states studied in Sec.IV for θ ′ = θ cannot be simply attributed in a one-to-one correspondence to the unperturbed states |ψ n ⟩, with n = 0, 1, 2, 3 at θ ′ = 0, because of the strong hybridization of the bands visible at the studied θ ′ ̸ = 0, in particular for n = 1, 2, 3; hence one should not use such perturbative estimates literally for all θ ′ of interest.
Nevertheless, it is suggestive, to use the estimate of the decay rates to interpret the decrease of the IPR peaks of Fig. 6 with increasing L uc .The Wigner-Weisskopf theory for the decay of a state in quantum-mechanics relates the decay rate with the width of the resonance in the frequency domain.At finite system size, this treatment can break down because the spectrum is discrete: in order to observe the finite width of a resonance, we need the level spacing to be much smaller than the width.Extrapolating our estimates of Γ n to the value θ ′ = θ = π/6 of Fig. 6 (well beyond the perturbative regime), we find Γ 0 ≈ 10 −3 and Γ 1,2,3 < 3 • 10 −6 .From our fits for the density of states close to the four peaks for θ ′ = θ = π/6 we are finding that ρ n ≈ f n L 2 uc with f n ≈ (0.04, 0.04, 0.03, 0.02), so the average level spacing ρ −1 n becomes of the order of Γ n for L uc ≈ (Γ n f n ) −1/2 , which results in L uc ≈ 150 for n = 0 and L uc > 2700 for n = 1, 2, 3.This estimate confirms that larger system sizes are needed in order to clearly observe the decay for the three rightmost peaks in the inverse participation ratio in Fig. 6.If the system size is not sufficiently large, this measure is sensitive to the fluctuations of the individual eigenergies, resulting in the noisy dependence observed in Fig. 6.

VI. CONCLUSIONS
In this work, we analysed the robustness of the bound states observed in the Google Quantum AI experiment [42].We compared the Hamiltonian and the Floquet spectrum, showing that the bound states in the N = 3 sector, which are protected by a gap in the Hamiltonian case, overlap (fold) with the other edge of the spectrum when trying to connect to the Floquet case.This is consistent with the direct study in the Floquet case where the bound states are surrounded by continuum states.We characterized the bound states by studying their overlaps with the trimer configurations and their inverse participation ratio resolved in sectors of total momentum.Our results suggest that similar many-body spectroscopic techniques as the ones applied to observe the bands of exact bound states in the integrable circuit can be used to detect the bands of approximate bound states in the perturbed circuit and to measure properties such as the maximal band velocity and their microscopic structure.For example, depending on the band, the particles in the bound state either reside primarily on the chain as in the trimers T I and T II in Fig. 1, or have one particle on the extra sites as in the trimers T III and T IV in the same figure.Such more detailed dynamical and structural properties of the bound states in fact change significantly as one varies θ ′ from 0 to θ (without much effect on their apparent robustness, in part because the bound states primarily mix among themselves) and could be probed directly in experiments.
For θ ′ = 0, particles located on the extra sites act as impurities in an integrable model on the chain.Recently, exact spatially bound states inside a continuous spectrum have been proposed in some integrable Hamiltonian models in the presence of an impurity [59,60].Interestingly, one of the 3-particle bound states in our modelling of the experimental system at θ ′ = 0, namely n = 2 state with ϵ ≈ 4.3614 in Sec.V from the sector (N 1∪2 , N 1 ′ ) = (2, 1), appears to correspond to a similar instance in the Floquet setting, see discussion in App.C (where this state is referred to as ϵ ≈ −1.92 from 2π shift).It would be interesting to investigate possible existence and stability of such states more broadly in the Floquet XXZ model with an impurity, both theoretically and experimentally.
For θ ′ ̸ = 0, while the bands of bound states are clearly visible even for fairly large system sizes, our finite size scaling analysis shows that the bound states tend to decay for increasing L uc .The decay is more rapid for one of the four bound states, due to the proximity with the 2+1 continuum.For other values of the parameters, other bound states can become similarly more unstable: we anticipate that such a difference in the robustness of the bound states can be probed in the same experimental apparatus, by preparing different initial states [such as the trimer configurations T III and T IV in Fig. 1 that would decay much faster for modified parameters as in Figs.7(b) and 8].A numerical analysis of the matrix elements and of the density of states in the unperturbed model confirms the presence of small but finite decay rates for all the bound states.
Our explanation for the current experiments on the non-integrable model is thus a quantitative few-body one and does not require true many-body unusual thermalization.An interesting question for future work is the possibility of so-called weak integrability-breaking perturbations for the Floquet XXZ model, and for Floquet integrable models in general.These perturbations, which can be systematically constructed to preserve integrability up to a given order in the perturbation strength, have been studied in the context of Hamiltonian integrable models.Understanding their possible structure in Floquet circuits would allow for the experimental verification of slow dynamics in digital quantum devices, such as the one used in the Google Quantum AI experiment.
Note added -A study by Hudomal et.al. [61], which appeared shortly after our work, also studies the bound states and the spectral properties of the GQAI experiment.Reference [61] focuses on different properties of the spectrum and studies the time evolution using TEBD simulations, and reaches different conclusions about the robustness of the bound states in the infinite time and infinite system size limit, which however may be hindered by their available time scales/sizes compared to potentially very long decay times.
(Note that the momentum k has exactly the same meaning as in the main text, since T 2 corresponds to translation by one unit cell.)Then, the states |γ, k⟩ are eigenstates of the Floquet operator U 0 with quasienergies ϵ = 2γ − k (mod 2π) . (A4) The model was solved in Ref. [44] using Bethe ansatz, and the following dispersion relation of the ℓ-particle bound state was found for generic ℓ cos (ϵ ℓ-string (k) − χ) = cos 2 (α) − sin 2 (α) cos(k), (A5) where , (A7) Our numerical study below of the spectra of U 0 and the bound states is consistent with these predictions.In Figs. 12 and 13 we plot the spectrum of the operators V 0 and U 0 for the sectors with N = 2 and N = 3 particles.For N = 2 (Fig. 12), we see that the band of bound states is gapped in the spectrum of V 0 , while one branch of the band in the spectrum of U 0 is not gapped for k ≈ 0: the "folding" procedure that gives the spectrum of U 0 from the one of V 0 brings part of the bound state band in the continuum of 1+1 states.However, the presence of a gap in V 0 implies that the bound states are robust to sufficiently small (but finite) perturbations of the Floquet operator that preserve the "brickwall" structure and the total number of particles (i.e., to the perturbations of U 0 that correspond to perturbations of the V 0 operator).This protection of the bound states is not manifest when one only considers the momentum-resolved spectrum of U 0 in the presence of the overlap with the continuum because such a view does not take into account the nontrivial conserved quantity V 0 , namely [U 0 , V 0 ] = 0; this additional "symmetry" is fully taken into account when one considers V 0 together with T 2 .
Similar arguments hold for N = 3: as we see in Fig. 13, part of the bound state band is gapped around k = 0 in the spectrum of V 0 .This part of the band is robust against any (small) perturbation of V 0 , despite the absence of a gap in the spectrum of U 0 .Note, however, that the comb teeth perturbation considered in the main text following the GQAI experiments does not preserve the brickwall structure that is crucial for the reduction of the full problem to V 0 , and the above protection does not operate in this case.
understanding of their characters.While the 3-particle bound states cease to be sharp in the Floquet system of interest, for reasons discussed in the main text the Hamiltonian system with modest folding in the quasienergy space still provides reasonable approximations, allowing much of the intuition about the spectral properties and bound states of the Hamiltonian case to be transferred to the Floquet case.
Hamiltonian spectra shown in this Appendix in Fig. 14 and top panels of Figs. 15 and 16 are the same as in the corresponding Hamiltonian panels in Figs. 2 and 3, except that they are not multiplied by (−π/6) and are instead presented here in the Hamiltonian energy units.At the expense of some repetition, this significantly simplifies referencing features in the spectrum and also puts the bound states at the bottom of the spectrum, allowing us to readily use intuition from familiar perturbation theory/effective Hamiltonian tools near ground states.In quantitative demonstrations below, we use the same parameters w = 1, u = u ′ = 4w, varying w ′ , as in the main text.
For easy reference, we first list the exact dispersions of the 2-particle (N = 2 sector) and 3-particle (N = 3 sector) bound states in the unperturbed Hamiltonian H 0 , Eq. ( 2), translated from the known results for the XXZ chain [62][63][64][65][66]: where k chain refers to the natural momentum on the chain with translation symmetry by one site [67].For u = 4w, the 2-string and 3-string bound states are separated from continuum states by finite gaps at each k chain ; the top of the 2-string bound state band at k chain = π happens to coincide with the bottom of the 1+1 continuum at k = 0, while the top of the 3-string bound state band lies significantly below the bottom of the corresponding 2+1 continuum (which lies below the 1+1+1 continuum).
Turning to the perturbed problem with the additional sites (the comb lattice in Fig. 1), we start with some general remarks.The conserved total number of particles in the system is denoted N .When w ′ = 0, the particle number on the original chain is separately conserved and is denoted N 1∪2 ; furthermore, the particle number on each (R, 1 ′ ) site is conserved (we often use a more crude conserved number N 1 ′ = N − N 1∪2 to group states).In this case, in the sector with N 1 ′ = 0, the u ′ term does not operate at all and the Hamiltonian is equivalent to the original integrable chain for any u ′ .Thus, in this sector, the u ′ part by itself does not break the integrability, and we sometimes refer to the model with w ′ = 0 as integrable.On the other hand, still keeping w ′ = 0, in sectors with N 1 ′ > 0, the occupations of the 1 ′ sites do not fluctuate but create static potentials (−u) < 0 for the chain particles on the 1 sites connected to the occupied 1 ′ sites.
The model for the particles on the chain is non-integrable because of these static "impurity potentials", although for small N 1 ′ it is possible to use some intuition/results from the integrable model away from the impurities.In the case of total of one particle in the system, the u and u ′ terms do not operate at all and the problem reduces to a single-particle problem with hopping amplitude w along the chain and w ′ on the 1-1 ′ links connecting to the extra sites.This problem is easily solved and has three bands

One particle
which are shown in Fig. 14.The symmetry E ↔ −E around the zero energy as well as the origin of the entire E = 0 band are due to the bipartite hopping nature of this single-particle problem: There are no on-site potentials and the hoppings connect only sites on different sublattices A ≡ {(R, 2), (R, 1 ′ )} and B ≡ {(R, 1)}.Since there are twice as many A sites as there are B sites, L uc zero-energy states residing entirely on the A sublattice.

Two particles
Figure 15(a-d) shows the momentum resolved exact diagonalization results in the Hamiltonian units, with fixed u = u ′ = 4w, where we vary w ′ from w ′ = 0 to w ′ = w; we set w = 1.We do not have an exact solution for general w ′ in this case, but we can understand rough features and in particular the 2-particle bound states by developing from certain controlled limits, which we now describe.
a. w ′ = 0, general w The eigenstates can be divided into three main groups (sectors) defined by specifying the (conserved) number of (N 1∪2 , N 1 ′ ) = (2, 0): There are L uc (2L uc − 1) such states and they are marked red in panel (a) in Fig. 15.This sector is equivalent to the unperturbed integrable model on the original chain.The corresponding states in Fig. 15(a) represent simply folding of the original chain spectrum to the new Brillouin zone.We clearly see the 2-particle continuum spanning energy window (−4w, 4w) in the thermodynamic limit.We are mainly interested in the lowest-energy states forming two bands of the 2-string bound states.In the original chain, the Hamiltonian is invariant under translation by a single site, and the bound states form a single band in the corresponding Brillouin zone, with dispersion given by Eq. (B1).This band is here folded to the new Brillouin zone, resulting in two distinct bands.
(N 1∪2 , N 1 ′ ) = (1, 1): There are 2L 2 uc such states and they are marked green in panel (a) in Fig. 15.These can be further subdivided into subgroups labelled by a location (R 0 , 1 ′ ) of the one particle on the 1 ′ sites, which remains completely localized since w ′ = 0.The spectra are identical for different R 0 ; each energy level is hence repeated L uc times in the full spectrum and gives a flat band in the energy vs momentum plot in Fig. 15(a).
For a fixed R 0 , we have a problem of one particle hop-ping on the original chain of 2L uc sites with the hopping amplitude w and a single attractive "impurity potential" (−u ′ ) < 0 felt by the particle when it is on the site (R 0 , 1).In this hopping problem, we expect one localized state near the impurity potential and 2L uc − 1 delocalized states.The delocalized states span energy window of (−2w, 2w) in the thermodynamic limit.
One can easily solve for the localized state ψ localized (j) = Ce −κ|j−j0| in the thermodynamic limit, where j labels sites as in the original chain and j 0 is the corresponding label of the site (R 0 , 1).The localized state energy and the rate of the wavefunction decay per lattice site are In the full system, this can be viewed as a 2-particle bound state with one of the particles immobile on the (R 0 , 1 ′ ) site and the other residing on the chain sites but dynamically bound to the immobile particle.For our numerical parameters u ′ = 4w we obtain 472 and e −κ = √ 5 − 2 ≈ 0.236; thus, the localization length is 0.693 of the original chain lattice spacing, which means we have a fairly compact bound state, and the above expression for ϵ localized is very accurate even for relatively small sizes.
As mentioned earlier, this analysis gives identical spectra for each of the L uc possible locations R 0 , which results in L uc -fold degeneracy for each found eigenvalue.Each such eigenvalue gives rise to a completely flat band when the full spectrum is resolved in momentum.Looking at the green states in Fig. 15(a) marking the present sector, we see the corresponding dense set of flat bands in the energy window (−2w, 2w) for the delocalized states and the flat band near ϵ localized for the 2-particle bound states.It is a numerical accident for the specific parameters that this energy is very close to where the two red 2-particle bound state bands meet, which happens at ϵ 2-string (k chain = π/2) = −u − 2w 2 /u = −4.5 [using Eq. (B1)].Note that in our case where the dominant interaction binding energies are taken to be the same, u ′ = u (motivated by the GQAI experiments), we expect all 2-particle bound states (i.e., chain-chain and chainextra site) to be roughly in the same ballpark, while the precise band locations depend on further dynamical details from the hopping energy.
(N 1∪2 , N 1 ′ ) = (0, 2): There are L uc (L uc − 1)/2 such states where both particles are on the 1 ′ sites.These states all have energy 0 and are marked blue in panel (a) in Fig. 15.They are not important in our considerations below focusing on the bottom of the spectrum near where the 2-string bound states reside.
b. Small w ′ , general w: qualitative considerations Now we add non-zero but small w ′ .The w ′ term does not act within the above sectors but connects the (N 1∪2 , N 1 ′ ) = (2, 0) and (1, 1) sectors [and also the (1, 1) and (0, 2) sectors].This roughly explains why in Fig. 15 (a-d) outside of the (−2w, 2w) energy window the (2, 0) sector states are not strongly affected for w ′ = 0.1 and w ′ = 0.5 (even all the way to w ′ = 1), except close to the boundary of this window and to the bound states from the (N 1∪2 , N 1 ′ ) = (1, 1) sector.The latter bound states (which we refer to as chain-extra site bound states) couple with the bound states from the (N 1∪2 , N 1 ′ ) = (2, 0) sector (chain-chain bound states) and with the bottom of the two-particle continuum from the (N 1∪2 , N 1 ′ ) = (2, 0) sector.There are some additional features in the middle of the spectrum that arise with increasing w ′ , but these are not of direct interest to us and are not studied further.
From now on we focus on the 2-string bound states.One of the effects of adding w ′ on the chain-chain bound states is that they can lower their energy by virtual processes where one of the particles hops off the chain onto a 1 ′ site and back.This lowering of the energy is roughly w ′2 /u and is visible in Fig. 15 (b-d) for the red 2-string bound state bands.On the other hand, the chain-extra site bound states, in the limit of very compact bound states, do not have this mechanism and to this order their energy would remain unchanged; such tendency is visible in Fig. 15 (b-d) for the green 2-string bound states.However, the chain-extra site bound states can hybridize with the chain-chain bound states with amplitude O(ww ′ /u) particularly when their energies are close, which happens near the wavevector π at w ′ = 0 and moves to smaller wavevectors as w ′ increases and the chain-chain bound states move to lower energies.This hybridization with the moving central location is visible in the progression in Fig. 15 (a-d) as we increase w ′ , while for w ′ = 1 the descendants of the chain-extra site bound states no longer overlap with the descendants of the chain-chain bound states.
From the figure, we see that the predominantly green 2string bound states survive for momenta sufficiently away from zero even for w ′ = 1, while they are in the continuum of states for momenta close to zero and will not survive in the thermodynamic limit.For the predominantly red 2-string bound states, only states in the upper band with momenta close to zero enter the continuum spectrum and will not survive in the thermodynamic limit (although their decay rate is likely very small), while the rest of this upper band and all of the lower band 2-string bound states clearly survive in the thermodynamic limit protected from the continuum by gaps at the corresponding momenta.
c. Perturbative treatment for w, w ′ ≪ u, u ′ Some of the above qualitative arguments can be made more precise by taking the limit w, w ′ ≪ u, u ′ , which we discuss here for completeness.In the absence of the hoppings, the lowest-energy states are the following dimer states specified by the particle locations on the full system (chain + extra sites), and depicted in Fig. 1(c): In the second-order perturbation theory, first we have diagonal corrections, which for later convenience we write in the ket-bra notation for the dimer states associated with location R as defined above: Here we explicitly see the claimed lowering of the energies of the chain-chain dimers D I , D II via virtual fluctuations involving w ′ hops, while no such lowering is present for the chain-extra site dimers D III .Next, in the same order, the chain-chain dimers can hop along the chain with amplitude w 2 /u.Explicitly, in the above notation,

Putting everything together, we have an effective Hamiltonian
By going to momentum space, we obtain a 3 × 3 matrix that is easy to diagonalize numerically and explore the evolution of the 3 bands, e.g., as one varies w ′ relative to w.The results are shown in Fig. 15(e-h) and capture roughly the behavior seen for the full problem and discussed qualitatively earlier: As w ′ /w increases from 0 to 1, the chain-chain dimer bands move to lower energies while at the same time particularly the higher-energy one of them hybridizes significantly with the chain-extra site dimer band, as seen in the mixing of the colors of these bands in Fig. 15.This model does not capture all detailed features seen in the figure, e.g., the apparent very small gaps between the two lowest bands at k = ±π, presumably because of an inaccurate treatment of the sizable w used in the figure; however, interestingly, at w ′ = w it gives a completely flat topmost band capturing the nearly flat uppermost 2-string bound state band in the figure.
We do not consider further details here since in any case this model does not include the eventual (small) overlaps with the two-particle continuum, where we have to use the full problem ED results.
3. Three particles a. w ′ = 0, general w The eigenstates can be divided into four main groups defined by specifying the number of particles on the chain N 1∪2 and on the extra sites N 1 ′ : (N 1∪2 , N 1 ′ ) = (3, 0): There are L uc (2L uc − 1)(2L uc − 2)/3 such states and they are marked red in panel (a) in Fig. 16.This sector is equivalent to the unperturbed integrable model on the original chain.The corresponding states in Fig. 16(a) represent simply folding of the earlier chain spectrum to the new Brilloin zone.We are primarily interested in the lowest-energy states forming two bands of the 3-string bound states representing the single band of these in the original chain, with dispersion given by Eq. (B2), folded to the new Brillouin zone.
(N 1∪2 , N 1 ′ ) = (2, 1): There are L 2 uc (2L uc − 1) such states and they are marked green in panel (a) in Fig. 16.These can be further subdivided into subgroups labelled by a location (R 0 , 1 ′ ) of the one particle on the 1 ′ sites, which remains completely localized under such Hamiltonian.The spectra are identical for different R 0 ; each energy is hence repeated L uc times in the full spectrum and gives a flat band in the energy vs momentum plot in Fig. 16(a).
For a fixed R 0 , we obtain a problem with two particles on the original chain but with an attractive "impurity potential" (−u ′ ) < 0 on one site j 0 corresponding to (R 0 , 1).Away from the impurity we have We are mainly interested in the effect of the impurity on the 2-string bound states.We can roughly model these as dimers [covering sites (j, j + 1)] hopping on the original chain [hops (j, j + 1) ↔ (j + 1, j + 2)] with amplitude w dimer and background energy εdimer , which we can estimate by fitting the exact dispersion of the 2-string bound states in the integrable model, Eq. (B1), as The dimer feels the attractive impurity at j 0 for two of its positions, (j 0 − 1, j 0 ) and (j 0 , j 0 + 1), and the prob-lem is mathematically equivalent to a point particle hopping on a lattice with potential (−u ′ ) on two neighboring sites.On an infinite lattice and in this model of a rigid dimer, we can solve analytically for exponentially localized states and obtain energies The first localized state is always present and is symmetric around j 0 , while the second localized state is present if u ′ > 2w dimer (which is satisfied in our problem) and is anti-symmetric around j 0 .These localized dimer states can be viewed as 3-particle bound states where one of the particles resides on the extra sites; we will often refer to these also as 3-string bound states.The two lowest green flat bands in Fig. 16(a) [close to the red weakly dispersive bands of the 3-string bound states from the (N 1∪2 , N 1 ′ ) = (3, 0) sector] correspond to these states viewed as bands once we include all different R 0 .Estimating w dimer = 0.25 from the 2-string dispersion, we can estimate the splitting of about 2w dimer ≈ 0.5, which is somewhat larger than the actual splitting of ≈ 0.25 between the corresponding green bands in the figure; the inaccuracy is likely due to crude modelling of the 2-string bound state by the above dimer picture.
We will not consider any other localized states in this group, which will be below the delocalized continuum of states but significantly above the 3-particle bound states of interest.Of interest for connecting with the Floquet case are the highest-energy states in this group, which are near the energy ≈ 4: If we take the Hamiltonian spectrum and multiply it by "time" π/6 as a rough estimate to connect with the GQAI Floquet experiment as done in the main text, the 3-string bound states from the sector (N 1∪2 , N 1 ′ ) = (3, 0), upon "folding" modulo 2π in the Floquet quasienergy space, would land among the states that are in the window ≈ [3,4] in the Hamiltonian spectrum in this figure.The nature of these states from the (N 1∪2 , N 1 ′ ) = (2, 1) sector is as follows: there is one particle on 1 ′ and two particles on the chain in scattering states staying away from the 1 ′ particle and from each other.We expect that matrix elements of the w ′ perturbation between these and the 3-string bound states from the (N 1∪2 , N 1 ′ ) = (3, 0) sector are very small, since the w ′ term would break such a 3-string state into a particle on a 1 ′ site and a nearby dimer, and both the proximity to the "impurity" and the proximity of the two particles on the chain has low probability in the described scattering states.In Sec.V we use such characters of the states involved to obtain scaling of these matrix elements with L uc .These matrix elements are not relevant at all in the Hamiltonian problem where the 3-string bound states and these states are separated by a large energy difference, but they are important for understanding robustness of the 3-string bound states in the Floquet problem.As discussed further in Sec.V, the density of states is also an important factor in estimating the decay rate of the bound states, requiring a quantitative study as presented in the main text.
(N 1∪2 , N 1 ′ ) = (1, 2): There are L 2 uc (L uc −1) such states and they are marked blue in panel (a) in Fig. 16.Here two particles are immobile on some 1 ′ sites, say (R 0 , 1 ′ ) and ( R0 , 1 ′ ), while the remaining particle is moving on the chain and sees the immobile particles as impurity potentials (−u ′ ) at j 0 and j0 corresponding to (R 0 , 1) and ( R0 , 1).Away from the impurities we have free propagation with the band covering energy window [−2w, 2w], while the impurities will localize some states out of this band.Since u ′ is sufficiently large, we expect that the two impurities will lead to two localized states out of the band, whose energies will depend somewhat on specific relative position of the impurities but will be independent of the overall shift of R 0 and R0 , leading to flat bands.These localized-state energies are visible in the spectrum around energy ≈ −4.5.Neither of the blue states are important for understanding the 3-string bound states in the Hamiltonian and the Floquet cases.
(N 1∪2 , N 1 ′ ) = (0, 3): There are L uc (L uc −1)(L uc −2)/6 such states and they are marked black in panel (a) in Fig. 16.These have zero energy and are not of much interest for the study of the 3-string bound states.
b. Small w ′ , general w: qualitative considerations We now focus solely on the 3-string bound states.One of the important effects of adding w ′ seen in Fig. 16 is that the energies of the 3-string states whose particles reside on the chain go down with increasing w ′ , while the energies of the 3-string states that have one particle on the extra sites remain essentially unchanged.We can understand this simply as follows.For the former 3-string states, which we will refer to as chain 3-string states [see also pictures T I and T II in Fig. 1(d) in the tight trimer limit], the non-zero w ′ allows virtual fluctuations involving hopping of one of the particles off the chain onto a nearby extra site site and back, leading to lowering of the energy.On the other hand, for the latter 3-string states, which we will refer to as dimer-1 ′ 3-string bound states [see also pictures T III and T IV in Fig. 1(d)], such virtual fluctuations are not available when the dimer is tightly bound to the particle on the 1 ′ site.The antisymmetric dimer-1 ′ 3-string bound state [schematically, antisymmetric combination of T III and T IV , see next subsection for more details] is separated from the chain 3string bound states already at w ′ = 0, and the separation only increases with adding w ′ .From the evolution in panels (a)-(d) in Fig. 16 and the essentially unchanged bright green color of the corresponding band even at w ′ = w, we conclude that the character of this state remains essentially unchanged.On the other hand, the symmetric dimer-1 ′ 3-string bound state energy at w ′ = 0 is close to the lowest energy of the upper chain 3-string band near k = ±π, and as the latter moves down upon increasing w ′ the two bands overlap and mix particularly near momenta where their energies are close.By the time w ′ reaches value of 0.5 the two bands are already separated and stay separated afterwards, also pushing a bit away from each other by level repulsion.We can then conclude that the darker-green band is roughly the symmetric dimer-1 ′ 3-string bound state with small admixture of fluctuations to the chain 3-string bound state.Both green bands remain essentially flat since moving such a dimer-1 ′ 3-string bound state requires at least four whops and two w ′ -hops, i.e., high order perturbation theory in the hoppings relative to the interactions.
Finally, when the chain 3-string bound states are well separated from the dimer-1 ′ 3-string bound states, we can understand the effect of w ′ on the former in more detail as follows: We start with the picture of a trimer hopping on the chain.In the presence of the extra sites, the trimer has two inequivalent positions: one where both ends of the trimer are over extra sites and the other where the middle of the trimer is over an extra site.In the former case, virtual fluctuations lower the energy of the trimer by 2w ′2 /u while in the latter case they lower the energy by only w ′2 /u.Thus, we can model the effect of small w ′ on the trimer as a potential −3w ′2 /(2u)+(−1) j w ′2 /(2u).This has both a uniform part shifting everything down in energy and a staggered part that will open a gap of roughly w ′2 /u at ±π/2 in the original chain Brilloin zone.Folded to the new Brilloin zone, we have a picture roughly similar to the two red bands with the gap near the corresponding Brilloin zone boundaries.The above picture gives an estimate of the gap between the bands as w ′2 /u = 0.25 at w ′ = w = 1, u = 4, which is somewhat larger but is still fairly close to the observed gap in Fig. 16(d).The inaccuracy is likely due to approximations when modelling the 3-string states by rigid trimers and also due to a larger admixture of the dimer-1 ′ in the upper band, as seen in the difference of the colors between the upper and lower red bands, which is not treated in the above trimer model.
To summarize, in this Hamiltonian system, besides the exact numerical results showing stable 3-string bound states well-separated from the continuum, we also have a fairly complete qualitative picture of the 3-string bound states.Note that even though we expect no decay of the 3-string bound states at w ′ = 1 similar to the integrable w ′ = 0 case, there are still significant quantitative differences in the 3-string propagation properties.For example, the maximal band velocity that determines the propagation wavefront is significantly smaller in the perturbed case.We expect similar quantitative effects also for (approximate) 3-particle bound states in the Floquet system, which could be tested in experiments.
c. Perturbative treatment for w, w ′ ≪ u, u ′ Some of the intuition including w ′ perturbatively relied on the pictures of tightly bound dimers or trimers.This is formally valid in the regime of small w ≪ u, and for w, w ′ ≪ u, u ′ we can perform the corresponding calculations systematically, which we do here for completeness.In the absence of the hoppings, the lowest energy states are the following trimer states specified by the particle locations on the chain + extra sites system [see Fig. 1(d)]: The states T I , T II form a degenerate manifold, and so do T III , T IV .The corresponding two manifolds are separated in energy if u ′ ̸ = u and would be treated separately in this case.On the other hand, they are degenerate if u ′ = u and should be treated together in this case.
Adding small w and w ′ , the first perturbative corrections appear in quadratic order.First, there are diagonal corrections appearing from virtual process processes where one of the particles hops away from the other two and then comes back, obtained by simply examining available such virtual moves: Next, at this order the above T III and T IV at the same R 0 get connected with matrix elements h eff III,IV = h eff IV,III = − The T III -T IV block is diagonalized by considering symmetric and anti-symmetric combinations, obtaining Finally, when u ′ = u, which we assume from now on, we need to also consider connections betweem T II and T III , T IV , which appear at this order: We see that at this order the T I is not coupled with the rest of the states and has the energy, Next, the anti-symmetric combination of T III and T IV also decouples and has the energy, sizes shows much larger fluctuations for this state than for the other bound states (see inset in Fig. 17).These fluctuations may be attributed to a weak hybridization with states in the continuum.The same fluctuations can be observed in the matrix elements V nj studied in Sec.V.
For completeness, we note that some states in the (N 1∪2 , N 1 ′ ) = (1, 2) sector (in blue), with quasienergy around ϵ ≈ 2.2, also exhibit size-independent I k=0 , even though they are not 3-particle bound states in the same sense as above.In this sector, two of the particles are completely localized on the extra sites and then serve as impurity potentials for the third particle that resides on the chain.This particle can be either in an extended state on the chain giving I k=0 ∼ L −1 uc (these states show collapse in the I k=0 • L uc panel in Fig. 17), or it can be localized on one of the impurities giving size-independent I k=0 ; it is the latter states, whose details also depend on the relative location of the two impurities, that show up near ϵ ≈ 2.2.

Continuum states
In Fig. 17 we observe three types of 2+1 continuum states that show collapse in the I k=0 • L uc panel.First, we find a set of states in the (N 1∪2 , N 1 ′ ) = (3, 0) sector (red) with quasienergies from 0.8 to −2.6 + 2π.Another set of states is in the (N 1∪2 , N 1 ′ ) = (2, 1) sector (green) with quasienergies from 1.9 to −3 + 2π: these can be of two subtypes, one where the 2-particle bound state has one particle localized on an extra site (the third particle roams on the chain), and the other where the 2-particle bound state roams on the chain (the third particle is localized on an extra site); we do not try to distinguish these here.Finally, most of the states in the (N 1∪2 , N 1 ′ ) = (1, 2) sector (blue) with quasienergies from −1 to 1 show collapse in the I k=0 • L uc panel: more precisely, these states are not 2+1 states but instead have two particles completely localized on the extra sites and the third particle extended around the chain, which we have already mentioned in the previous subsection.
Turning to the states that show collapse in the I k=0 • L 2 uc panel, we see the 1+1+1 continuum, that spans the full quasienergy range in the (N 1∪2 , N 1 ′ ) = (3, 0) sector (red), and a range of quasienergies from −2.1 to 2.1 in the (N 1∪2 , N 1 ′ ) = (2, 1) sector (green) (strictily speaking, the latter are not 1+1+1 continuum but are from flat bands obtained by constructing momentum eigenstates from degenerate states where one of the particles is completely localized on an extra site while the other two are in extended states on the chain).
Note that in the quasienergy range from 1.2 to 2 we find many eigenstates in the (N 1∪2 , N 1 ′ ) = (2, 1) sector with values of I k=0 that are intermediate between the 2+1 and 1+1+1 continuum: since this sector is not integrable, the states in the 2+1 can decay in the 1+1+1 continuum (if they have the same quasienergy), so the scattering states will be a mixture of 2+1 and 1+1+1 states.We have not tried to understand these states in any detail since they are far from the 3-particle bound states of main interest to us.

Figure 1 .
Figure 1.(a) 1d chain hosting the integrable Hamiltonian and Floquet systems.(b)1d chain with extra sites forming a comb in the GQAI experiment, with the comb teeth connected to every other site of the original chain; a unit cell and labelling of sites in each cell is shown.(c) Schematic of distinct two-particle bound states, which become dimers in the very strong binding limit, in a unit cell on the comb; these give rise to three bands in the momentum space.(d) Distinct three-particle bound states, schematized as trimers, in a unit cell on the comb, producing four bands in the momentum space.

Figure 2 .
Figure 2. Spectrum of the Hamiltonian model (left side in each two-panel group) and Floquet circuit (right side) in the oneparticle (N = 1, top panels) and two-particle (N = 2, bottom panels) sectors.The spectra are shown for comparable values of parameters (w = 1, u = u ′ = 4, θ = π/6, ϕ = ϕ ′ = 2π/3), and for increasing values of w ′ /w and θ ′ /θ taken to be the same in each two-panel group.The Hamiltonian energies E are multiplied by an appropriate time evolution factor of (−π/6) for an approximate match to the Floquet circuit step.These still fit within a 2π period and hence can be compared directly with the Floquet quasienergies ϵ without requiring any folding.For the unperturbed models with w ′ /w = 0 and θ ′ /θ = 0, different colors represent states belonging to different (N1∪2, N 1 ′ ) sectors, as shown in the legend.For the perturbed models, (N1∪2, N 1 ′ ) are not good quantum numbers, and the eigenstates can have non-zero components in all the sectors: The color mixing (red/green for N = 1 and red/green/blue for N = 2) represents the squares of the norms of the projections of each eigenstate into the sectors [i.e., RGB color (255 • P (1, 0), 255 • P (0, 1), 0) for N = 1 and (255 • P (2, 0), 255 • P (1, 1), 255 • P (0, 2)) for N = 2, where P (N1∪2, N 1 ′ ) is the squared norm of the projection of the eigenstate in the (N1∪2, N 1 ′ ) sector].The system size is Luc = 120 unit cells for N = 1, and Luc = 20 unit cells for N = 2.

Figure 3 .
Figure 3. Spectrum of the Hamiltonian model (left side in each two-panel group) and Floquet circuit (right side) in the three-particle sector (N = 3).The spectra are shown for comparable values of parameters (w = 1, u = u ′ = 4, θ = π/6, ϕ = ϕ ′ = 2π/3), and for increasing values of w ′ /w = θ ′ /θ.The Hamiltonian energies are multiplied by the time evolution factor of (−π/6), and this gives a window wider than 2π and would require folding to compare with the Floquet quasienergies.For clarity, instead of such folding of the Hamiltonian spectrum, we show the Floquet circuit quasienergies repeated with period 2π.The Floquet Brillouin zone ϵ ∈ [−π, π) is marked with grey dashes on the quasienergy scale.For easier comparison with the Hamiltonian case, copies of the quasienergy spectrum are shown beyond the first Floquet Brillouin zone.For the unperturbed models with w ′ /w = 0 and θ ′ /θ = 0, different colors represent states belonging to different (N1∪2, N 1 ′ ) sectors, as shown in the legend in the leftmost group.For the perturbed models, the color mixing red/green/blue represents the square of the norms of the projections of each eigenstate in the sectors with (N1∪2, N 1 ′ ) = (3, 0), (2, 1), (1, 2) respectively; specifically, the color in RGB is (255 • P (3, 0), 255 • P (2, 1), 255 • P (1, 2)).[Note that when w ′ /w = θ ′ /θ = 0 all states in the sector (0, 3) have zero energy and their color is RGB (0,0,0), which is black.]The system size we choose is rather small (Luc = 12 unit cells) to avoid overwhelming the plots.

Figure 4 .
Figure 4. Top panels: Probability PT I + PT II of having the three particles in consecutive sites along the chain [Eqs.(13), (14)].Lower panels: Probability PT III + PT IV of having one particle on an extra site and the other two right next to it on the chain [Eqs.(15), (16)].The system size is Luc = 36 unit cells.The system parameters are the same as in Fig.3.A narrow window of quasienergies is shown focusing on the three-particle bound states, revealing their character of being primarily chain trimers or chain-extra site trimers and also where significant mixing is present.In the Hamiltonian system the three-particle bound states are isolated from the continuum and persist in the thermodynamic limit even for w ′ /w = 1, while in the Floquet case they are inside a continuum and will decay in the thermodynamic limit but apparently survive to fairly large sizes.The similarity between the bound states in the Hamiltonian and Floquet systems is notable, allowing to infer properties of the latter ones as well as of the surrounding continuum from the more simple Hamiltonian understanding.

Figure 5 .
Figure 5. Inverse participation ratio (for fixed momentum k).As explained in the text, this is calculated over basis states used in the momentum-resolved ED, which can be essentially viewed as describing configurations of the particles relative to their center of mass.The system size is Luc = 36 unit cells.The system parameters and the quasienergy window is the same as in Fig.4.The IPR detects both types of the three-particle bound states equally well and provides a good measure of the degree of localization of the particles in the bound state.Detailed study for varying system sizes in the Floquet system at θ ′ = θ (showing the strongest decay of the bound states in this figure) and k = 0 is performed in Fig.6.

Figure 6 .
Figure 6.Top panel: Inverse participation ratio of Floquet eigenstates with momentum k = 0 for the system with parameters θ = θ ′ = π/6, ϕ = ϕ ′ = 2π/3.The states are organized by the quasienergy on the horizontal axis, and the prominent (approximate) bound states stand out with their large IPR.Multiple system sizes are shown.Lower panels: Maximum of I k=0 in an energy window [ϵ − ∆, ϵ + ∆] (shaded areas in the top panel), with ∆ = 0.05, as a function of the system size.The result of a linear fit, performed to bring out the overall decreasing trend, is shown with a dashed line in each panel.

Figure 7 .
Figure 7. Spectrum of the perturbed Floquet circuit focusing on the continuum states near the 3-particle bound states: The 2+1 continuum is characterized by larger probability of trimer configurations (PT I + PT II and PT III + PT IV ) and larger I k than the 1+1+1 continuum.Note that the quantities are the same as the ones plotted in Figs. 4 and 5 but the color scales are saturated to a maximum value of 0.01 to bring out the 2+1 states more.(a) For the choice of parameters considered in the experiment, the bound states are separated in the spectrum from the 2+1 continuum, but one of them is significantly closer.(b) For a different choice of parameters, two out of four bound states lie in the 2+1 continuum (for some values of k).

Figure 8 .
Figure 8. Top panel: Inverse participation ratio of Floquet eigenstates with momentum k = 0 for the system with parameters θ = θ ′ = π/6, ϕ = ϕ ′ = 8π/15 from Fig. 7(b).For small system size (Luc = 10), four distinct eigenstates stand out with a large value of I k=0 .Lower panels: Maximum of I k=0 in an energy window [ϵ − ∆, ϵ + ∆] (shaded areas in the top panel), with ∆ = 0.05, as a function of the system size.The two leftmost panels (the ones which correspond to the eigenstates overlapping with the 2+1 continuum) show a clear decay of the height of the peak for increasing system size.The other two show an extremely weak decay with Luc.

Figure 9 .
Figure 9. Top: absolute value of the matrix elements Vnj (connecting the bound state |ψn⟩ with the eigenstate ϵj) multiplied by the system size Luc, as a function of the quasienergy ϵj.The data collapse shows that |Vnj| ∝ L −1 uc consistent with 1+1+1 continuum.Bottom: inverse level spacing of neighboring pairs of Floquet eigenstates (∆ϵ) −1 multiplied by L −2 uc , as a function of the average quasienergy of the pair.Only states that have a non-zero matrix elements with |ψn⟩ are considered.The data collapse shows the predicted scaling for the density of states (∆ϵ) −1 ∝ L 2 uc .The vertical grey lines indicate the energy ϵn of the bound state.The Floquet system parameters are ϕ = ϕ ′ = 2π/3, θ = π/6, θ ′ = 0, i.e., the "integrable point".The data are plotted for the bound states with n = 1 (left) and n = 3 (right), which belong to the sector with (N1∪2, N 1 ′ ) = (3, 0).

Figure 12 .
Figure 12.Spectrum of the operators V0 (left) and U0 (right) in the N = 2 sector.The colormap indicates the momentumresolved inverse participation ratio I k (the scale is saturated to a maximum value of 0.10 for more clear view).The bands of 2-particle bound states have large I k (dark blue), while the 1+1 continuum states have small I k (light green).

Figure 13 .
Figure 13.Spectrum of the operators V0 (left) and U0 (right) in the N = 3 sector.The colormap indicates the momentumresolved inverse participation ratio I k (the scale is saturated to a maximum value of 0.15).The bands of 3-particle bound states have the largest I k (dark blue), followed by the 2+1 continuum (light blue), while the 1+1+1 continuum states have small I k (light green).

Figure 14 .
Figure 14.Single particle (N = 1) spectrum of H for w = 1, u = u ′ = 4.The color mixing represents the weights on the different sectors of the unperturbed model w ′ = 0, as described in the caption of Fig. 2.

Figure 15 .
Figure 15.(a-d) Two particle (N = 2) spectrum of H for w = 1, u = u ′ = 4, Luc = 20 .The color mixing represents the weights on the different sectors of the unperturbed model w ′ = 0, as described in the caption of Fig. 2. (e-h) Zoom-out of the region containing the 2-particle bound states.Semitransparent lines: 2-particle bound states (and nearby other states) in the spectrum of H for the same parameters as in (a-d) but larger Luc = 36 used for smoothness.Solid lines: Spectrum of the effective Hamiltonian Eq. (B9) derived for w, w ′ ≪ u, u ′ .The exact spectrum and the results of the effective model have almost perfect overlap for all the three bands of bound states.

Ĥeff hop = − w 2 u R
|D II (R)⟩⟨D I (R)| + H.c.+ |D I (R)⟩⟨D II (R + 1)| + H.c. .Note that the chain-extra site dimers cannot hop by themselves at this order [for u ̸ = u ′ , such hopping can appear only at O(w ′2 w 4 /u ′5 )].Finally, for u ′ = u, assumed from now on, where we need to do degenerate perturbation theory involving all D I (R), D II (R), D III (R), the D III (R) can convert to D I (R) or D II (R) and vice versa: the 2-particle continuum (covering energy window of [−4w, 4w] = [−4, 4] for w = 1) as well as the band of 2-string bound states (energy window of [−u − 4w 2 /u, −u] = [−4, −5] for u = 4w used here).The attractive impurity will lead to appearance of some localized states out of these.

Figure 16 .
Figure 16.(a-d) Three particle (N = 3) spectrum of H for w = 1, u = u ′ = 4, Luc = 12.The color mixing represents the weights on the different sectors of the unperturbed model w ′ = 0, as described in the caption of Fig. 3. (e-h) Zoom-out of the region containing the 3-particle bound states.Semitransparent lines: bands of the 3-particle bound states for the same parameters as in (a-d) but larger Luc = 36.Solid lines: (flat) bands of the 3-particle bound states from solving the second-order effective Hamiltonian, Eqs.(B12)-(B14).In panel (e), the lowest flat band of the effective Hamiltonian (E = −2u−2w 2 /u = −8.5) is three-fold degenerate: one band belongs to the (N1∪2, N 1 ′ ) = (2, 1) sector and the others are in the (N1∪2, N 1 ′ ) = (3, 0) sector (only the first one is visible, in green).