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, which host excitations in the form of microwave photons, are promising candidates as they provide a configurable lattice in which a discrete number of photons can be confined to a qubit site, hop between the sites and interact with each other. The tunability of coupling elements enables 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 \(\left|0\right\rangle \) and \(\left|1\right\rangle \) 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 methodological shift in the studies of correlated systems^{12,13,14,15,16}. Whereas theoretical studies of condensed matter models were focused on Hamiltonian systems for many decades, high-fidelity quantum processors commonly operate on the basis of unitary gates rather than continuous Hamiltonian dynamics. This experimental access to periodic (Floquet) unitary dynamics opens the door to a plethora of non-equilibrium phenomena^{17}. Because such periodic dynamics often cannot be described in terms of a local Hamiltonian, established results are few and far 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 benchmark quantum processors. A relevant example of an interacting integrable model is the one-dimensional (1D) quantum spin-½ XXZ model, which is known to support bound states^{2,3,4,5,21}. Recently, the shared symmetries of the spin-½ XXZ Hamiltonian model with its Floquet counterpart led to a proof for the integrability of the XXZ Floquet quantum circuits^{22,23,24}. Later, Aleiner obtained the full spectrum for these Floquet systems and provided analytical results for bound states^{25}. The advantage of using quantum processors in studying these models becomes apparent when going beyond the integrability limit, where analytical and numerical techniques fail to scale favourably.

To define systems with bound states, consider a chain of coupled qubits and the unitary evolution \(\hat{U}\) of interacting photons on this array. We divide the computational space of all bitstrings with *n*_{ph} photons into two sets: one set \({\mathcal{T}}\) is composed of all bitstrings in which all photons are in adjacent sites, for example, \(\left|00…011100…00\right\rangle \); the other set \({\mathcal{S}}\) contains all other *n*_{ph} bitstrings, for example, \(\left|00…101001…00\right\rangle \). A bound state is formed when the eigenstates of the system can be expanded as the superposition of bitstrings mainly in \({\mathcal{T}}\) and with smaller weight in \({\mathcal{S}}\). Therefore, for any initial state \(\left|{\psi }_{0}\right\rangle \in {\mathcal{T}}\) the photons remain adjacent at all future times \(\left|\psi \right\rangle =\hat{U}\left|{\psi }_{0}\right\rangle \), which implies that almost every projective measurement returns a bitstring in \({\mathcal{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, for example, 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 for 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 in which a maximum of 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, two-qubit fSim gates are applied between all pairs in the ring. In the two-qubit subspace, \(\{\left|00\right\rangle ,\left|01\right\rangle ,\left|10\right\rangle ,\left|11\right\rangle \}\), this gate can be written as

$${\rm{fSim}}(\theta ,\phi ,\beta )=\left(\begin{array}{cccc}1 & 0 & 0 & 0\\ 0 & \cos \theta & {\rm{i}}{{\rm{e}}}^{{\rm{i}}\beta }\sin \theta & 0\\ 0 & {\rm{i}}{{\rm{e}}}^{-{\rm{i}}\beta }\sin \theta & \cos \theta & 0\\ 0 & 0 & 0 & {{\rm{e}}}^{{\rm{i}}\phi }\\ \end{array}\right),$$

(1)

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 \(\left|11\right\rangle \) state on interaction of two adjacent photons. In the Supplementary Information we show that we can achieve this gate with high fidelity (approximately 99%) for several angles. In the following, we will denote \({\rm{fSim}}(\theta ,\phi ,\beta =0)\) as \({\rm{fSim}}(\theta ,\phi )\). The qubit chain is periodically driven by a quantum circuit, with the cycle unitary:

$${\hat{U}}_{{\rm{F}}}=\prod _{{\rm{even}}\,{\rm{bonds}}}{\rm{fSim}}(\theta ,\phi ,\beta )\prod _{{\rm{odd}}\,{\rm{bonds}}}{\rm{fSim}}(\theta ,\phi ,\beta ).$$

(2)

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-\langle \hat{{Z}_{j}}\rangle )/2\) on each site *j* as a function of circuit depth for the fSim angles *θ* = π/6 and *ϕ* = 2π/3. Because the fSim gates are excitation number conserving, all data are postselected for the bitstrings with the proper number of excitations, which allows us to mitigate errors induced by population decay. Although *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 two-photon bound state. For *n*_{ph} > 2, the concentration of the population near the centre indicates that the photons do not disperse far, but instead stay close to each other. In the Supplementary Information we also present the situation in which the initial photons are not adjacent, in which case the system tends towards a uniform distribution.

To extract the wavefront velocity, we select the measured bitstrings in which the photons remain adjacent, that is, in \({\mathcal{T}}\), and discard the ones in \({\mathcal{S}}\). In Fig. 2c, we present the spatially and temporally resolved probabilities of the ‘centre of photon mass’ (CM, Fig. 2b) of these \({\mathcal{T}}\) bitstrings. With this selection, the first panel in Fig. 2c shows a very similar pattern to the single-particle propagation in Fig. 2a, 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}_{\,\text{g}}^{\text{max}}\) (Fig. 2d), match very well with those corresponding to the analytical dispersion relations derived in ref. ^{25}, which take the same functional form for all *n*_{ph}:

$$\cos (E(k)-\chi )={\cos }^{2}(\alpha )-{\sin }^{2}(\alpha )\cos (k),$$

(3)

where *E* is the quasi-energy, *k* is the momentum, and *α* and *χ* are functions of *n*_{ph}, *θ* and *ϕ* (see Supplementary Information for exact forms).

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}_{{\mathcal{T}}}/({n}_{{\mathcal{T}}}+{n}_{{\mathcal{S}}})\) (where \({n}_{{\mathcal{T}}({\mathcal{S}})}\) is the number of bitstrings in \({\mathcal{T}}({\mathcal{S}})\)), which reflects contributions from both internal unitary dynamics as well as external decoherence (Fig. 2e). In the absence of dephasing, \({n}_{{\mathcal{T}}}\) should reach a steady-state value after the observed initial drop. However, we observe a slow decay, which we attribute to the dephasing of the qubits, as the data are postselected 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^{26}. We measure the energy of the bound states by comparing their accumulated phase over time relative to the vacuum state \({\left|0\right\rangle }^{\otimes {N}_{{\rm{Q}}}}\). This is achieved by preparing *n*_{ph} adjacent qubits in the \(\left|+X\right\rangle \) state and measuring the following *n*_{ph} body correlator that couples the bound states with the vacuum state:

$$\langle {C}_{j,{n}_{{\rm{ph}}}}\rangle =\langle {\Pi }_{i=j}^{j+{n}_{{\rm{ph}}}-1}{\sigma }_{i}^{+}\rangle =\langle {\Pi }_{i=j}^{j+{n}_{{\rm{ph}}}-1}({X}_{i}+i{Y}_{i})\rangle $$

(4)

for all sets of *n*_{ph} adjacent qubits (Fig. 3a). This protocol is based on measuring the Green’s function of the system. Whereas the correlator above is not Hermitian, it can be reconstructed by measuring its constituent terms (for example, 〈*X*_{j}*X*_{j+1}〉 − 〈*Y*_{j}*Y*_{j+1}〉 + *i*〈*X*_{j}*Y*_{j+1}〉 + *i*〈*Y*_{j}*X*_{j+1}〉 for *n*_{ph} = 2) and summing these with the proper complex prefactors. We note that as \({C}_{j,{n}_{{\rm{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 \(\left|000..00\right\rangle +\left|00..0110..00\right\rangle \). By expanding these states in the momentum basis (*k*-space), it becomes evident that \(\langle {C}_{j,{n}_{{\rm{ph}}}}\rangle \) contains the phase information needed to evaluate the dispersion relation of the *n*_{ph} bound states:

$$\begin{array}{rcl}\left|\psi (t)\right\rangle & = & \frac{1}{\sqrt{2}}\left({\left|0\right\rangle }^{\otimes {N}_{{\rm{Q}}}}+{\sum }_{k}{\alpha }_{k}{{\rm{e}}}^{-i\omega (k)t}\left|k\right\rangle \right)\\ & \to & \langle {C}_{j,{n}_{{\rm{ph}}}}\rangle =1/(2\sqrt{{N}_{{\rm{Q}}}}){\sum }_{k}{\alpha }_{k}^{* }{{\rm{e}}}^{i(\omega (k)t-kj)},\end{array}$$

(5)

where \(\left|k\right\rangle \) and *α*_{k} are bound *n*_{ph} photon momentum states and their coefficients, respectively.

Figure 3b shows the real and imaginary parts of the correlator for the case of two photons. Whereas the real-space data display a rather intricate pattern (Fig. 3b), conversion to the energy and momentum domain through a two-dimensional (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π per cycle due to the discrete time translation symmetry of the Floquet circuit, are in agreement with the predictions of equation (3), as illustrated in coloured 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.

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*θ* (ref. ^{25}). 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 behaviour from the more strongly interacting case studied in Fig. 3c: although no band is observed for most momenta, a clear state emerges near *k* = ±π per site. Interestingly, this observation of a bound state in the weakly interacting regime can be attributed to destructive interference of the decay products of the bound state: a two-photon bound state \(\left|..0110..\right\rangle \) can separate into two possible states, \(\left|..1010..\right\rangle \) and \(\left|..0101..\right\rangle \), which are shifted relative to each other by one lattice site. Hence, they destructively interfere when the momentum is near *k* = ±π per site, which prevents separation. (See the Supplementary Information for band structures of additional fSim angles.)

External magnetic fields can shift the energy bands and reveal the electric pseudo-charge of the quasiparticles constituting the band. We produce a synthetic magnetic flux **Φ** that threads the ring of qubits by performing *Z* rotations 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 (ref. ^{27}). As a consequence, the eigenstates are expected to attain a phase *j*(*n*_{ph}*β*), effectively shifting their quasi-momentum by *n*_{ph}*β*. Figure 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^{25}. 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 relies on integrability or symmetries^{28,29}. 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 Fig. 3 demonstrate an experimental verification of this remarkable theoretical prediction in the integrable limit and constitute 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 in which 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 as used in Fig. 2: *θ* = π/6 and *ϕ* = 2π/3. In the third layer we also use \({\phi }^{{\prime} }=2{\rm{\pi }}/3\) but vary the swap angle \({\theta }^{{\prime} }\) to tune the strength of the integrability breaking perturbation.

Figure 4c shows the probability of measuring three-photon \({\mathcal{T}}\)-bitstrings as a function of time for various \({\theta }^{{\prime} }\) angles. In the limit of small \({\theta }^{{\prime} }\), for which the integrability breaking is weak, the system shows a slowly decaying probability, similar to the unperturbed (integrable, \({\theta }^{{\prime} }=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 three photons (\({\mathcal{T}}\)+\({\mathcal{S}}\)). However, the surprising finding is that even up to \({\theta }^{{\prime} }={\rm{\pi }}/6\), which corresponds to perturbation gates identical to the gates on the main ring, that is, a strong perturbation, there is very little decay in \({n}_{{\mathcal{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 \({\theta }^{{\prime} }\) dependence of the half-width of the band (Fig. 4f). Indeed, we find that the bandwidth is insensitive to \({\theta }^{{\prime} }\) 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, in which certain initial states exhibit periodic revivals and do not thermalize^{30,31}. Moreover, in the case of weak integrability breaking, robustness to perturbations can result from quasi-conserved or hidden conserved quantities^{32,33}. However, the resilience observed here extends well beyond the weak integrability breaking regime typically considered in such scenarios^{34}. Alternatively, the presence of highly incommensurate energy 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^{35,36}. 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 \({\rm{fSim}}\) gates to change the ratio of kinetic to interaction energy, as well as directly measuring multibody correlators \(\langle {C}_{j,{n}_{{\rm{ph}}}}\rangle \), both of which are hard to access in conventional solid-state and atomic physics experiments. Aided by these capabilities, we observed the formation of multiphoton 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.