Formation of robust bound states of interacting microwave photons

Systems of correlated particles appear in many fields of modern science and represent some of the most intractable computational problems in nature. The computational challenge in these systems arises when interactions become comparable to other energy scales, which makes the state of each particle depend on all other particles1. The lack of general solutions for the three-body problem and acceptable theory for strongly correlated electrons shows that our understanding of correlated systems fades when the particle number or the interaction strength increases. One of the hallmarks of interacting systems is the formation of multiparticle bound states2–9. Here we develop a high-fidelity parameterizable fSim gate and implement the periodic quantum circuit of the spin-½ XXZ model in a ring of 24 superconducting qubits. We study the propagation of these excitations and observe their bound nature for up to five photons. We devise a phase-sensitive method for constructing the few-body spectrum of the bound states and extract their pseudo-charge by introducing a synthetic flux. By introducing interactions between the ring and additional qubits, we observe an unexpected resilience of the bound states to integrability breaking. This finding goes against the idea that bound states in non-integrable systems are unstable when their energies overlap with the continuum spectrum. Our work provides experimental evidence for bound states of interacting photons and discovers their stability beyond the integrability limit.

Photons that propagate in vacuum do not interact with each other; however, many technological applications and the study of fundamental physics require interacting photons.Consequently, realizing quantum platforms with strong interactions between photons constitutes a major scientific goal 10,11 .In this regard, superconducting circuits are promising candidates since they provide a configurable lattice where microwave photons can be confined to a qubit site, hop between the sites, and interact with each other.Each site can host localized electromagnetic oscillations and hence be occupied with a discrete number of microwave photon excitations.The tunability of coupling elements allows photons to hop between the sites, and the non-linearity of qubits leads to interaction between the photons.The zero-and single-photon occupancies of qubits are used as the |0 and |1 states in quantum information processing.Here we also confine the dynamics to zero or single occupancy for a given qubit, the so-called hard core boson limit, and show that microwave photons can remain adjacent and form coherent bound states.
The advent of quantum processors is giving rise to a paradigm shift in the studies of correlated systems [12][13][14][15][16] .While theoretical studies of condensed matter models were focused on Hamiltonian systems for many decades, high-fidelity quantum processors commonly operate based on unitary gates rather than continuous Hamiltonian dynamics.This experimental access to periodic (Floquet) unitary dynamics opens the door to a plethora of nonequilibrium phenomena 17 .Since such periodic dynamics often cannot be described in terms of a local Hamiltonian, established results are fewer and far in between [18][19][20] .For instance, until recently, there was no theoretically known example of bound state formation for interacting Floquet dynamics.
Integrable models form the cornerstone of our understanding of dynamical systems and can serve to bench-mark quantum processors.A relevant example of an interacting integrable model is the 1D quantum spin-1/2 XXZ model which is known to support bound states [2][3][4][5]21 . Recntly, the shared symmetries of the spin-1/2 XXZ Hamiltonian model with its Floquet counterpart led to a proof for the integrability of the XXZ Floquet quantum circuits 22,23 .Later, Aleiner obtained the full spectrum for these Floquet systems and provided analytical results for bound states 24 .The advantage of using quantum processors in studying these models becomes apparent when going beyond the integrability limit, where the classical counterpart of the circuit shows chaos, and analytical and numerical techniques fail to scale favorably.
To define systems with bound states, consider a chain of coupled qubits and the unitary evolution Û of interacting photons on this array.We divide the computational space of all bitstrings with n ph photons into two sets: one set T composed of all bitstrings in which all photons are in adjacent sites, e.g.|00...011100...00 ; the other set S contains all other n ph bitstrings, e.g.|00...101001...00 .A bound state is formed when the eigenstates of the system can be expanded as the superposition of bitstrings mainly in T and with smaller weight in S. Therefore, for any initial state |ψ 0 ∈ T the photons remain adjacent at all future times |ψ = Û |ψ 0 , which implies that almost every projective measurement returns a bitstring in T (Fig. 1a).
The emergence of a thermodynamic phase or the formation of a bound state in Floquet dynamics seems rather implausible at first sight.In a closed Floquet system there is no notion of lowest energy, a key concept in equilibrium physics.Therefore, the energy minimization that commonly stabilizes bound states in e.g.atoms does not hold.In the absence of interactions and in 1D, photons hop independently and the evolution can be mapped to that of free fermions.In this limit, obviously, no bound state can  be formed.The key question of bound state formation is whether the effect of kinetic energy (hopping) that moves photons away from each other could be balanced by interactions.In Fig. 1b, we provide a plausibility argument to illustrate this point.Consider two photons that are initially occupying adjacent sites, in the low kinetic energy regime where maximum one hopping event occurs in the span of a few cycles.In the spirit of Feynman path formulation, the probability of a given configuration at a later time can be obtained from summing over all possible paths that lead to that configuration with proper weights.When photons are in adjacent sites, they accumulate phase due to the interaction.In the three depicted paths, the accumulated phases are different, thus leading to destructive interference.Hence, the interactions suppress the probability of unbound configurations and facilitate the formation of bound states.The control sequence used to generate unitary evolution in our experiment consists of a periodic application of entangling gates in a 1D ring of N Q = 24 qubits (Fig. 1c).Within each cycle, 2-qubit fSim gates are applied between all pairs in the ring.In the 2-qubit subspace, {|00 , |01 , |10 , |11 }, this gate can be written as where θ and β set the amplitude and phase, respectively, of hopping between adjacent qubit lattice sites, and the conditional-phase angle φ imparts a phase on the |11 state upon interaction of two adjacent photons.In the supplementary information, we show that we can achieve this gate with high fidelity (∼ 1%) for several angles.
In the following, we will denote fSim(θ, φ, β = 0) as fSim(θ, φ).The qubit chain is periodically driven by a quantum circuit, with the cycle unitary: In the limit of β = 0 and θ, φ → 0, this model becomes the Trotter-Suzuki expansion of the XXZ Hamiltonian model.
To quantify to what extent photons remain bound together, we prepare an initial state with n ph photons at adjacent sites and measure the photon occupancy of all sites after each cycle with approximately 5,000 repetitions.In Fig. 2a we plot the average photon occupancy (1− Ẑj )/2 on each site j as a function of circuit depth for the fSim angles θ = π/6 and φ = 2π/3.Since the fSim gates are excitation number conserving, all data are post-selected for the bitstrings with the proper number of excitations, which allows us to mitigate errors induced by population decay.While n ph = 1 is not a bound state, it provides a benchmark, where we can clearly see the quantum random walk of a single particle and its familiar interference pattern.For n ph = 2, we observe the appearance of two wavefronts: the fastest one corresponds to unbound photons, whereas the other one corresponds to the 2-photon bound state.For n ph > 2, the concentration of the population near the center indicates that the photons do not disperse far, but instead stay close to each other.In the supplementary information, we also present the situation where the initial photons are not adjacent, in which case the system tends toward a uniform distribution.
To extract the wavefront velocity, we select the measured bitstrings in which the photons remain adjacent, i.e. in T , and discard the ones in S. In panel c, we present the spatially and temporally resolved probabilities of the "center of photon mass" (CM, Fig. 2b) of these T bitstrings.With this selection, the first panel in c shows a very similar pattern to the single-particle propagation in a, highlighting the composite nature of the bound state.The propagation velocities of the bound states can now be easily seen, and as expected, the larger bound states propagate more slowly.The wavefronts propagate with constant velocity, indicating that the bound photons move ballistically and without effects of impurity scattering.The extracted maximum group velocities of the bound states, v max g (Fig. 2d), match very well with that corresponding to the analytical dispersion relations derived in ref 24 , which take the same functional form for all n ph : where α and χ are functions of n ph , θ, and φ (see supplementary information for exact forms).
In order to characterize the stability of the bound state, it is useful to consider the evolution of the fraction of bitstrings in which the photons remain adjacent, n T /(n T + n S ) (where n T (S) is the number of bitstrings in T (S)), which reflects contributions from both internal unitary dynamics as well as external decoherence (Fig. 2e).In the absence of dephasing, n T should reach a steadystate value after the observed initial drop.However, we observe a slow decay which we attribute to the dephasing of the qubits, since the data is post-selected to remove T 1 photon loss effects.A remarkable feature of the data is that the decay rate for various n ph values is the same, indicating that this decay is dominated by bond breaking at the edges of the bound state.
To show that the bound photons are quasiparticles with well-defined momentum, energy, and charge, we study the spectrum of the bound states using a many-body spectroscopy technique 25 .We measure the energy of the bound states by comparing their accumulated phase over time relative to the vacuum state |0 ⊗N Q .This is achieved by preparing n ph adjacent qubits in the |+X -state and measuring the following n ph -body correlator that couples the bound states with the vacuum state: for all sets of n ph adjacent qubits (Fig. 3a).This protocol is based on measuring the Green function of the system.While the correlator above is not Hermitian, it can be reconstructed by measuring its constituent terms (e.g. and summing these with the proper complex pre-factors.We note that since C j,n ph only couples the n ph -photon terms to the vacuum, the initial product state used here serves the same purpose as an entangled superposition state |000..00 + |00..0110..00 .By expanding these states in the momentum basis (k-space), it becomes evident that C j,n ph contains the phase information needed to evalu-ate the dispersion relation of the n ph bound states: where |k and α k are bound n ph -photon momentum states and their coefficients, respectively.
Fig. 3b shows the real and imaginary parts of the correlator for the case of two photons.While the real space data displays a rather intricate pattern (Fig. 3b), conversion to the energy and momentum domain through a 2D Fourier transform reveals a clear band structure for both the single-particle and the many-body states (Fig. 3c).The observed bands, which are defined modulo 2π/cycle due to the discrete time translation symmetry of the Floquet circuit, are in agreement with the predictions of Eq. 3, as illustrated in colored dashed curves.The bands shift when the photon number increases, as expected from the higher total interaction energy.Moreover, they become flatter, a characteristic feature of increased interaction effects.
In order for a bound state to form, the interaction energy must be sufficiently high compared to the kinetic energy of the particles.In particular, bound states are only expected to exist for all momenta when φ > 2θ 24 .To explore this dependence on φ/θ, we also measure the band structure for n ph = 2 in the weakly interacting regime (θ = π/3, φ = π/6; Fig. 3d), which exhibits very different behavior from the more strongly interacting case studied in Fig. 3c: while no band is observed for most External magnetic fields can shift the energy bands and reveal the electric pseudo-charge of the quasi-particles constituting the band.We produce a synthetic magnetic flux Φ that threads the ring of qubits by performing Zrotations with angles ±Φ/N Q on the qubits before and after the two-qubit fSim gates, resulting in a complex hopping phase β = Φ/N Q when a photon moves from site j to j +1 26 .As a consequence, the eigenstates are expected to attain a phase (n ph β)•j, effectively shifting their quasimomentum by n ph β.Fig. 3e displays the flux dependence of the two-photon band structure, exhibiting a clear shift in momentum as Φ increases.In Fig. 3f, we extract the shift for n ph = 1−5 and observe excellent agreement with the theoretical predictions 24 .Crucially, the momentum shift is found to scale linearly with n ph , indicating that the observed states have the correct pseudo-charge.
Generally, bound states in the continuum are rare and very fragile, and their stability rely on integrability or symmetries 27,28 .Familiar stable dimers, such as excitons in semiconductors, have energy resonances in the spectral gap.In the system considered here, the bound states are predicted to almost always be inside the continuum due to the periodicity of the quasi-energy.Our results shown in Figs. 3 demonstrate an experimental verification of this remarkable theoretical prediction in the integrable limit and constitutes our first major result.
Next we probe the stability of the bound states against integrability breaking.Fermi's golden rule suggests that any weak perturbation that breaks the underlying symmetry will lead to an instability and a rapid decay of the bound states into the continuum.We examine the robustness of the n ph = 3 bound state by constructing a quasi-1D lattice where every other site of the 14 qubit ring is coupled to an extra qubit site (Fig. 4a).The extra sites increase the Hilbert space dimension and ensure that the system is not integrable.We implement the circuit depicted in Fig. 4b with three layers of fSim gates in each cycle.The first two layers are the XXZ ring dynamics with the same parameters used in Fig. 2: θ = π/6 and φ = 2π/3.In the third layer we also use φ = 2π/3 but vary the swap angle θ to tune the strength of the integrability breaking perturbation.
Fig. 4c shows the probability of measuring three-photon T -bitstrings as a function of time for various θ angles.g, Momentum averaged quasi-energy spectra for varying θ fitted with Lorentzian.The bound state peak slowly disappears with the increase of θ .
In the limit of small θ , where the integrability breaking is weak, the system shows a slowly decaying probability, similar to the unperturbed (integrable, θ = 0) results presented in Fig 2 .In Fig. 4d, we show the dependence of this probability on perturbation strength after two fixed circuit depths.For strong perturbations, the integrability breaking washes out the bound state and the probability rapidly decays to the equiprobable distribution in the full Hilbert space of 3 photons (T +S).However, the surprising finding is that even up to θ = π/6, which corresponds to perturbation gates identical to the gates on the main ring, i.e. a strong perturbation, there is very little decay in n T .This observation demonstrates the resilience of the bound state to perturbations far beyond weak integrability breaking for n ph = 3.We further confirm this finding by performing spectroscopy of these states, which shows the presence of the n ph = 3 bound states up to large perturbations (Fig. 4e).By fitting the momentum averaged spectra (Fig. 4g), we extract the θ -dependence of the half-width of the band (Fig. 4f).Indeed, we find that the bandwidth is insensitive to θ up to very large perturbation.
These observations are at odds with the expectation that non-integrable perturbation leads to the fast decay of bound states into the continuum.One known exception is many-body scars, where certain initial states exhibit periodic revivals and do not thermalize 29,30 .Moreover, in the case of weak integrability breaking, robustness to perturbations can result from quasi-conserved or hidden conserved quantities 31,32 .However, the resilience observed here extends well beyond the weak integrability breaking regime typically considered in such scenarios.Alternatively, the presence of highly incommensurate en-ergy scales in the problem can lead to a very slow decay in a chaotic system due to parametrically small transition matrix elements, a phenomenon called prethermalization 33,34 .Our experiment finds the survival of an integrable system's feature -bound states -for large perturbation and in the absence of obvious scale separation, which may point to a new regime arising due to interplay of integrability and prethermalization.
The key enabler of our experiment is the capability of tuning high fidelity fSim gates to change the ratio of kinetic to interaction energy, as well as directly measuring multi-body correlators C j,n ph , both of which are hard to access in conventional solid state and atomic physics experiments.Aided by these capabilities, we observed the formation of multi-photon bound states and discovered a striking resilience to non-integrable perturbations.This experimental finding, although still observed for computationally tractable scales, in the absence of any theoretical prediction, constitutes our second major result (Fig. 4).A proper understanding of this unexpected discovery is currently lacking.

I. Quantum processor details and coherence times
The experiment is performed on a quantum processor with similar design as that in Ref. [? ].The qubits are transmons with tunable frequencies and interqubit couplings.Figure S1a shows the single-qubit relaxation times of the 24 qubits used in the experiment, where a median value of T 1 = 16.1 µs is found.The dephasing times T * 2 , measured via Ramsey interferometry, are shown in Fig. S1b and have a median value of 5.3 µs.Lastly, The T 2 values after CPMG dynamical decoupling sequences are also shown in Fig. S1b and have a median of 17.8 µs.

II. 2-Qubit fSim gates A. fSim calibration
The floquet unitaries used in the experiment are composed of alternating layers of fSim(θ, φ, β) gates which are defined as: where θ is the SWAP angle and φ is the conditional phase, and β is phase accumulated resulting from hopping.For open chains, β is not gauge invariant and can be ignored, but for closed chains, non-zero β values lead to a total flux threading the closed chain.fSim(θ, φ, β) describes the unitary form output by a DC pulse bringing the fundamental frequencies ω 1 and ω 2 of two transmons into resonance and turning on their interqubit coupling g for a given time duration t p , as illustrated in Fig. S2a.During t p , resonant interaction between the |10 and |01 states of the two transmons leads to population transfer and a finite θ.Additionally, dispersive interaction between the |11 and |02 (as well as |20 ) states of the two qubits gives rise to a finite conditional phase φ.
Due to the frequency detunings of the qubits during the DC pulse, the fSim unitary also includes additional singlequbit Z rotations and is more generally described by: where γ, α and β are complex phases incurred by the single-qubit Z rotations.These single-qubit phases are calibrated and reduced to nearly zero using the technique of Floquet calibration described in our previous works [? ??].Here we focus on the tuning and calibration of the two-qubit angles θ and φ. Figure S2c and Figure S2d show experimentally obtained values of θ and φ as functions of pulse parameters t p and g max .In these measurements, the technique of unitary tomography [? ] is used to estimate the angles.We have also enforced a Gaussian filter with time constant 5 ns on the rising/falling edges of the pulse on g to ensure adiabatic 3 ) (square), ( 1 6 , 1 2 ) (triangle) and ( 13 , 1 6 ) (diamond).c, Experimentally measured θ as a function of pulse length tp and maximum interqubit coupling gmax.The approximate pulse parameters for the fSim gates in panel b are indicated by their corresponding symbols.d, Similar plot as panel c but with φ shown instead.e, θ and φ as functions of tp for a fixed gmax/2π of 36 MHz.Solid lines show linear fits.f, θ and φ as functions of gmax for a fixed tp of 5 ns.Solid lines show a linear fit ∝ gmax to θ(gmax) and a quadratic fit ∝ g 2 max to φ(gmax).
evolution with respect to the |11 → |20 and |11 → |02 transitions.This is important to minimize leakage.We observe that θ shows a series of maxima/minima corresponding to values of t p and g max where |01 is fully transformed to |10 or returned back to |01 .On the other hand, φ increases monotonically over t p and g max until it reaches a maximum value of π where it is wrapped by 2π and becomes −π.Given the dependence of both unitary angles on both pulse parameters t p and g max , independent control of θ or φ is not possible with a single pulse parameter.As such, past works have added a resonant pulse between the |11 and |02 states to enact a pure CPHASE gate, thereby enabling full tunablity over θ and φ [? ].The additional pulse, however, significantly increases the complexity of quantum control and is also prone to leakage.In this work, we have chosen to perform fSim gates directly using the single pulse in Fig. S2a.Our approach relies on the different scaling of θ and φ with the pulse parameters, as illustrated by Fig. S2e and Fig. S2f.Here we observe that while θ and φ both scale linearly with t p , the scaling with g max is different for the two angles: whereas θ scales linearly with g max , φ ∝ g 2 max due to the fact that dispersive shift of the |11 state by the |02 and |20 states is proportional to g 2 /∆, where ∆ is the frequency difference between |11 and the |02 (|20 ).The difference in scaling implies that it is possible to achieve a desired combination of θ and φ by choosing a particular "contour" in the 2D space of (t p , g max ) where θ has the target value, then increasing (decreasing) g max while decreasing (increasing) t p until φ attains the target value as well.
Practically, φ and θ are calibrated via a simple gradient descent method: we start with an initial guess (t 0 , g 0 ) for the pulse parameters (t p , g max ) based on the 2D scan shown in Fig. S2c and Fig. S2d.The corresponding values of φ and θ are then accurately determined via Floquet calibration [? ??] which we denote as φ 0 and θ 0 .We then calibrate the fSim angles at (t p , g max ) = (t 0 + δt, g 0 ) and (t p , g max ) = (t 0 , g 0 + δg).The results allow us to approximate the following gradient matrix: A new set of pulse parameters (t 1 , g 1 ) are then computed from the gradient matrix and the deviations from target fSim angles, (∆t, ∆g) = (t c − t 0 , g c − g 0 ), via: The fSim angles (t p , g max ) are then measured at the new pulse parameters (t 1 , g 1 ) and the process is repeated.Generally only two gradient descents are sufficient to reach control errors on the level of ∼20 mrad for both θ and φ.
When repeating the fSim gate n times in Floquet calibration, the error accumulates as:

B. fSim gate control error
The angular errors, which are measured using periodic Floquet calibration as outlined in section II A, are displayed in Fig. II A along with the measured 2-qubit Pauli error for the four pairs of θ,φ studied in our work.The angular errors of the fSim gate can be combined into an overall control error by calculating the coherent gate infidelity: where e P is the Pauli error, D = 4 is the size of the computational subspace of the gate, and U target(actual) is the target (actual) unitary.Inserting the unitary matrix in eq.S2, we find the control error in terms of the angular errors: For the four angle pairs in our study, we find median control errors of c (θ = π/6, φ Qubit pair  C. Choice of parameters θ and φ The choice of the angles θ and φ of the fSim gate is dictated by several considerations, with regards to both the physics of the bound state and the experimental parameters of the fSim gate.In the main text, we present several angles as they have differents properties: • First and foremost, the behavior is significantly different for the gapped φ > 2θ regime and the gapless regime.
We have predominantly focused on the gapped regime where bound states exist for all momenta, simplifying the analysis.
• Secondly, as shown in Fig. S4b, the amount of dispersion in the spectrum increases with increasing (decreasing) values of θ (φ).Hence, the ratio φ/θ should not be too large in order to maximize the visibility of the dispersion.
• Higher values of φ cause the bound state to be more localized, as intuitively expected from the stronger interactions.
• Finally the quality of the fSim gate depends on the angles, since strong interaction requires longer gates time or larger interaction strength (see section II).
FIG. S5.Pulse sequences for determining the spectra of a, a single excitation, and b, two excitations for 4 qubits.
Multi-excitation example.Next, we consider the photon conserving dynamics of two excitations in a 4 qubit system, with the computational basis |1100 , |1010 , |1001 , |0011 , |0101 , |00110 .We start by placing Q3 and Q4 in a |0 + |1 superposition, which gives rise to a superposition of the vacuum and a two-excitation state, as well as undesired single excitation states: The appearance of single excitation terms is not desired and could have been avoided by using proper entangling gates to arrive at |0000 + |0011 as a more relevant initial state.However, as our experimental results show, populating (wrong) manifolds with fewer number of excitations, is not harmful since the two-qubit lowering operator does not couple these manifolds to the vacuum.The existence of these undesired states does, however, reduce the signal contrast, due to distributed probabilities.The evolution of the initial states results in: The two-qubit lowering operator acting on Q3 and Q4, can also be written as Hence, the presence of terms in other manifolds does not leading to wrong answers and ψ t | σ − 3 σ − 4 |ψ t provides the desired 2-photon spectrum.The number of independent Pauli strings that need to be measured scales exponentially with the number of excitations in the manifold, and hence this method is not scalable to large photon number systems.4), we post select the outcome bitstrings that preserve the number of photons.This post-selection is justified by the fact that the fSim gate is an excitation preserving gate for any angle.The observed decay is due to the T1 decay of the qubits during the circuit.

V. Supplementary data
The measurement becomes more susceptible to T1 errors as the size of the bound state increases, thus increasing the necessary number of repetitions to construct the statistics.For each number of photons, we initialize the system with adjacent excitations (blue) or with excitations separated by a few vacancies (orange).In the case where the excitations are initially adjacent, we find that the distribution stays concentrated in the part of Hilbert space with few vacancies.When the initial state contains vacancies, on the other hand, the evolution tends toward exploring the entire Hilbert space and approaches the equiprobable distribution represented by a dashed line.Dashed lines: theoretical prediction.We observe very good agreement with the theoretically predicted band structure across all the three angle pairs that satisfy φ > 2θ, required to have bound states at all momenta.As predicted by theory, the bands are found to shift with increasing interaction strength (φ), and the width of the dispersion increases with increasing θ.In the regime where φ < 2θ, we observe a two-photon bound state for momenta near k = ±π, while no bound states are observed for higher n ph in this regime.The latter is likely due to the fact that the overlap between the initial product state and the n ph -bound state scales as 2 −n ph , thus causing a reduction in the signal-to-noise ratio at high n ph .Dashed black vertical lines: theoretically predicted momentum threshold for the existence of bound states.for n ph = 1 − 5. c,d, Same as a,b, but Fourier transformed to momentum space.e, Absolute value of the multi-particle correlator Fourier transformed to momentum space.f-h, Linecuts at k = 5π/6 for c, d, and e, respectively.Black curves show fits on the form A k=5π/6,t = αe (iω−1/τ )t (single fit for all three plots at each photon number).Inset in center panel: extracted decay rate τ −1 (in units of cycles −1 ) as a function of photon number.

3 FIG. 2 .
FIG. 2. Trajectory of bound photons.a, Time-and site-resolved photon occupancy on a 24-qubits ring for photon numbers n ph = 1 − 5. To measure a n ph -photon bound state, n ph adjacent qubits are prepared in the |1 state.b, Schematic and example of bitstrings in T and S. Center of mass is defined as the center of n ph adjacent occupied sites.c, Evolution of the center of mass of n ph -bound states.Each trajectory is similar to the single photon case, highlighting the composite nature of the bound states.d, Extracted maximum group velocity from the trajectory of the center of mass.Black line: theoretical prediction.e, Decay of the bound state into the single excitations continuum due to dephasing.For all panels, θ = π/6 and φ = 2π/3, and the trajectories are averaged over all possible initial states.Data are post-selected for number of excitations equal to n ph .
FIG. 3. Band structure of multi-photon bound states.a, Schematic of circuit used for many-body spectroscopy.n ph adjacent qubits are prepared in the |+ -state, before evolving the state with a variable number of fSim gates.The phase of the bound state is probed by measuring the correlator σ + i ..σ + i+n ph −1 for all sets of n ph adjacent qubits.b, Real (top) and imaginary (bottom) parts of the n ph = 2 correlator.c, Band structure for n ph = 1 − 5 (top to bottom), obtained via a 2D Fourier transform in space and time of the n ph -correlators.Color scale: absolute square of the Fourier transform, |A k,ω | 2 .Dashed curves: theoretical prediction in Eqn. 3. d, Band structure for n ph = 2 in the weakly interacting (φ < 2θ) regime, displaying the emergence of a bound state only at momenta near k = ±π.Dashed black lines: theoretically predicted momentum threshold for the existence of the bound state (see supplementary information).e, Flux dependence of the n ph = 2 band structure, displaying a gradual momentum shift as the flux increases (Φ0 = 2πNQ).Orange circles and dashed line indicate the peak position of the band.f, Extracted momentum shifts as a function of flux for n ph = 1 − 5 (top to bottom), indicating that the rate of shifting scales linearly with the photon number of the bound states, i.e. the pseudo-charge q of each bound state is proportional to its number of photons.Colored lines: theoretical prediction.

FIG. 4 .
FIG. 4. Resilience to integrability breaking.a, Schematic of the 14 qubits chain with 7 extra sites in red to break the integrability.b, Integrability is broken via an extra layer of fSim-gates (red) between the chain and the extra qubits, with φ = φ and a gradually varied θ .c, Decaying probability of remaining bound for different swap angles θ .Similar to Fig. 2e, the bound state decays into the continuum due to the dephasing.d, Probability of remaining bound after 20 and 40 cycles as θ is swept.e, Spectroscopy of the n ph = 3 bound state for different θ .Note that the bound state survives even for θ = θ.f, Half width of the momentum averaged spectra (g) as a function of θ .The gray line indicates the result for the chain without the extra qubits.g, Momentum averaged quasi-energy spectra for varying θ fitted with Lorentzian.The bound state peak slowly disappears with the increase of θ .

6 FIG
FIG. S2. a, Schematic illustration of the DC pulse used to realize a fSim gate.b, Schematic plot showing the four sets of fSim angles used in this work:(θ/π, φ/π) = ( 1 3 , 5 6 ) (filled circle), ( 1 6 , 2 3 ) (square), ( 1 6 ,12 ) (triangle) and (13 , 1 6 ) (diamond).c, Experimentally measured θ as a function of pulse length tp and maximum interqubit coupling gmax.The approximate pulse parameters for the fSim gates in panel b are indicated by their corresponding symbols.d, Similar plot as panel c but with φ shown instead.e, θ and φ as functions of tp for a fixed gmax/2π of 36 MHz.Solid lines show linear fits.f, θ and φ as functions of gmax for a fixed tp of 5 ns.Solid lines show a linear fit ∝ gmax to θ(gmax) and a quadratic fit ∝ g 2 max to φ(gmax).
FIG.S4.a, Band structure of the different bound states for θ = π/6 and φ = 2π/3.b, c, Width of the n ph -photon band structure as a function of the parameters of the fSim gate.All these results are calculated using the exact solution given in eq.S12.
FIG. S6.Post-selection for Figure 2: In the trajectory experiments shown in the main text (Figures2 and 4), we post select the outcome bitstrings that preserve the number of photons.This post-selection is justified by the fact that the fSim gate is an excitation preserving gate for any angle.The observed decay is due to the T1 decay of the qubits during the circuit.The measurement becomes more susceptible to T1 errors as the size of the bound state increases, thus increasing the necessary number of repetitions to construct the statistics.
FIG. S7.Trajectory histogram:Complementary data for Figure2in the main text with angles θ = π/6 and φ = 2π/3.For each number of photons, we initialize the system with adjacent excitations (blue) or with excitations separated by a few vacancies (orange).In the case where the excitations are initially adjacent, we find that the distribution stays concentrated in the part of Hilbert space with few vacancies.When the initial state contains vacancies, on the other hand, the evolution tends toward exploring the entire Hilbert space and approaches the equiprobable distribution represented by a dashed line.
FIG. S8.Many-body spectroscopy for various fSim angles.a, Band structures analogous to those shown in Fig. 3c in the main text for n ph = 1 − 5 (top to bottom) and fSim angles θ = π/3, φ = 5π/6 (a), θ = π/6, φ = 2π/3 (b), θ = π/6, φ = π/2 (c), θ = π/3, φ = π/6 (d).Dashed lines: theoretical prediction.We observe very good agreement with the theoretically predicted band structure across all the three angle pairs that satisfy φ > 2θ, required to have bound states at all momenta.As predicted by theory, the bands are found to shift with increasing interaction strength (φ), and the width of the dispersion increases with increasing θ.In the regime where φ < 2θ, we observe a two-photon bound state for momenta near k = ±π, while no bound states are observed for higher n ph in this regime.The latter is likely due to the fact that the overlap between the initial product state and the n ph -bound state scales as 2 −n ph , thus causing a reduction in the signal-to-noise ratio at high n ph .Dashed black vertical lines: theoretically predicted momentum threshold for the existence of bound states.
FIG. S9.Decoherence of multi-particle correlator.a,b, Real (a) and imaginary (b) parts of the correlator Cj,n ph = Π j+n ph −1 i=j σ + i

2 FIG 5 FIG. S11 .
FIG. S10.Flux dependence of band structures for n ph = 1 − 5 (a-e, respectively).Momentum shifts were extracted by convolving the band structures with that at Φ = 0 (Φ = 0.2Φ0 for n ph = 5 due to more clear structure), summing over the energy axis, and finding the maximum.Colored dots indicate the corresponding extracted peak positions of the bands.
FIG.1.Bound states of photons.a, In a 1D chain of qubits hosting bound states, an initial state with adjacent photons evolves into a superposition of states in which the photons remain bound together.
b, Interactions between photons can lead to destructive interference for paths in which photons do not stay together, thus suppressing separation.c, Schematic of the gate sequence used in this work.Each cycle of evolution contains two layers of fSim gates that connect the even and odd pairs respectively.The fSim gate has three controllable parameters that set the kinetic energy (θ), the interaction strength (φ) and a synthetic magnetic flux (β).The median gate fidelity, measured with cross-entropy benchmarking, is 1.1% (see supplementary information).