{'fraction_non_alphanumeric': 0.039935080135457796, 'fraction_numerical': 0.015949062875513037, 'mean_word_length': 4.448979591836735, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 23, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Interacting Qubit-Photon Bound States with Superconducting Circuits

Qubits strongly coupled to a photonic crystal give rise to many exotic physical scenarios, beginning with single and multi-excitation qubit-photon dressed bound states comprising induced spatially localized photonic modes, centered around the qubits, and the qubits themselves. The localization of these states changes with qubit detuning from the band-edge, offering an avenue of in situ control of bound state interaction. Here, we present experimental results from a device with two qubits coupled to a superconducting microwave photonic crystal and realize tunable on-site and inter-bound state interactions. We observe a fourth-order two photon virtual process between bound states indicating strong coupling between the photonic crystal and qubits. Due to their localization-dependent interaction, these states offer the ability to create one-dimensional chains of bound states with tunable and potentially long-range interactions that preserve the qubits' spatial organization, a key criterion for realization of certain quantum many-body models. The widely tunable, strong and robust interactions demonstrated with this system are promising benchmarks towards realizing larger, more complex systems of bound states.

realize these bound states due to their intrinsically tailorable band structure, and characteristic Bloch mode electric field distribution [8] which enables access to strong coupling with qubits [9][10][11][12][13]. Bound states in multiqubit photonic crystal devices are an ideal platform to study many-body quantum optics in one-dimensional systems [6,7,[14][15][16][17][18][19]. Unlike many qubits coupled to a common cavity mode but similar to the case of some optical multimode cavities [20,21], coupling to a band edge creates bound states that intrinsically preserve the spatial organization of qubits, offering the ability to create onedimensional chains of bound states with tunable and potentially long-range interactions. The promise of engineering interaction profiles beyond the intrinsic flip-flop with additional microwave drive tones further opens the possibility of simulating a wide range of quantum spin models in future devices [17]. In this paper, we demonstrate and characterize the underlying, fundamental tunable on-site and interbound-state interactions in a superconducting microwave photonic crystal device coupled to two transmon qubits.
A single dressed bound state, seeded by a single qubit in a crystal, is itself a unique avenue of study. Liu et al. first directly detected such a bound state in a stepped-impedance microwave crystal coupled to a single transmon qubit [13]. That work characterized the dependence of localization length on detuning between the transmon qubit and the band edge and further confirmed the existence of the localized state in the band gap when the bare transmon qubit is in the passband-an unmistakable signature of non-Markovianity, as seen in Figs. 1(e), 1(f), and Appendix D. We note that the transmission dip observed in Figs. 1(e) and 1(f) is due to the reflection from the transmon qubit [22][23][24]. State localization is tunable in situ with frequency through a range determined by device parameters, including transmon qubit-waveguide coupling and band curvature. Compared with previous work, we attain increased localization in this device [ Fig. 1(b)] due mainly to a flatter band dispersion, realized by tailoring the unit cell of the photonic crystal (see Appendix A for a detailed discussion of the experimental parameters of our system). The bound-state localization length in this device is still widely tunable, which is critical for realizing strong, tunable interaction between spatially separated bound states. As the different coupling regimes translate to dramatically altered system behavior [7], it is important to determine which domain our system falls under. In systems such as the one presented here, qubit emission into the waveguide being larger than the other decay rates (coherent atom-photon interaction rates larger than decay rates) is the minimal coupling criterion, upon which the dressed bound state within the gap can be spectrally resolved [7]. The strong coupling criterion corresponds to the situation where a bare qubit resonant with the band edge gives rise to a bound state that is shifted from the band edge by more than the bound state's linewidth [7,13]. In our finite system, we observe an approximately 250-MHz separation between the bound state and the band edge with bound-state linewidth of 4 MHz when a qubit is resonant with the band edge, thus firmly reaching the strong coupling condition [see Figs. 1(b), 1(e), and 1(f)]. By fabricating two transmon qubits in the photonic crystal [see Fig. 1(a) and Appendix B for a discussion on coupling transmons to photonic crystals], we realize multiple, spectrally resolvable bound states and can study interbound-state interaction.
The nature of interbound-state interaction makes this platform intrinsically well suited for investigating one-dimensional chains of bound states [see Fig. 1(c)]. Realizing sizable chains is possible by increasing the number of unit cells-a property that does not impact the Bloch mode distribution or band dispersion. Thus, qubits can be in separate unit cells but realize nearly identical coupling to the band edge. As the strength of interbound-state interaction depends on the spatial overlap of the photonic wave functions with the qubits, the distance separating qubits (set by device design) is directly mapped into the interactions of the system, maintaining the chainlike interaction pattern. Furthermore, in the investigation of bound states, the finite size of the crystal is a practical advantage: the overlaps of bound states with ends of the crystal lead to only quasi-bound states, allowing for detection through transmission measurements (see Appendix C 2 for a proof). Direct detection of such a state in this way was first demonstrated in Liu et al. [13].
Controlling photon-mediated interaction between superconducting qubits has been demonstrated in other one-dimensional systems-such as two qubits in a cavity [25,26] or in a linear waveguide [27]. However, in these cases the distance between the qubits was effectively eliminated (i.e., standing-wave interaction in a cavity) or otherwise reduced (modulo wavelength in a linear dispersion waveguide). Thus, photonic crystals and tunable bound states offer a fundamentally distinct form of interaction.
In addition to determining localization length, the frequency of the bound state also determines on-site interaction strength. In Figs. 2(a) and 2(b), we characterize the dependence of the transition frequencies between the three lowest levels of the bound state on bare transmon qubit frequency, and observe the steady reduction in bound-state anharmonicity from over 350 MHz to 0 MHz as the transmon qubit is tuned from deep in the band gap to the passband. Here, we have defined the bound-state anharmonicity as Δ = 2ω 01 − ω 12 , where ω 01 and ω 12 are the dressed bound-state frequencies [see level diagram in Fig. 2(a)]. This is dramatically more than the approximately 10% modification of transmon qubit anharmonicity with frequency expected when a transmon qubit is strongly coupled to a cavity mode [28].
Therefore, while we may treat the one-excitation and two-excitation bound states as first ( | 1 ) and second ( | 2 ) excited states of a new effective anharmonic transmon qubit [13], it is important to note that this effective transmon qubit differs in frequency and anharmonicity from the bare multilevel transmon qubit. Defining the three lowest bare transmon qubit levels as |0〉 |1〉, and |2〉, here the two-excitation bound is largely due to the coupling of the second transmon qubit transition (|1〉 ↔ |2〉 with the band edge rather other the multiphoton effects [7,15] (see Appendix F 1).
Numerical simulations, modeling the photonic crystal as a coupled cavity array with free parameters fit to match the band curvature from the dispersion relation [7,29] (see Appendix C for details), show similar dependence of anharmonicity on detuning [see Fig. 2(a) inset and Fig. 2(b)]. Unlike the transfer matrix method [30][31][32], this approach can extend beyond the single-excitation manifold to capture the higher levels of the bound state, as well as the Lamb shift of the qubit frequency, observed when including next-nearest-neighbor hopping between coupled cavities. Each transmon qubit is modeled as a three-level ladder (unless otherwise mentioned) with negative anharmonicity, and with the |0〉 ↔ |1〉 and |1〉 ↔ |2〉 transitions coupled with amplitudes g and g 2, respectively, to its local cavity site. It is critical to include level |2〉 to accurately reproduce the two-excitation manifold observed in experiment.
The tunable level structure also emerges in the emission spectrum of a continuously driven bound state [Figs. 2(c) and 2(d)], induced by a single qubit with bare frequency above the band edge. At low drive amplitude or Rabi frequency, transmission across the crystal via the bound state exhibits antibunching [see Fig. 2(c), inset] [33], consistent with single-photon transport of a two-level system and resonant pump (see Appendix D 4) [34][35][36]. As Rabi frequency is increased, we see a Mollow triplet emission spectrum [ Fig. 2(c)], characteristic of a driven two-level system. When the Rabi frequency is on the same order as the anharmonicity, the bound state can no longer be approximated as a two-level system. In this domain, the steady state will be a mixture of the three eigenstates obtained by diagonalizing the drive Hamiltonian in the Hilbert space spanned by |0〉 |1〉, and |2〉. Transitions between all eigenstates result in six sidebands [37]. These six sidebands are visible in Fig. 2(d), though emission intensity varies greatly among them due to eigenstate population. A seventh transition is evident in the data [7.25 GHz in Fig. 2(d)]. This additional transition is due to the fourth effective transmon qubit level ( | 3 ) while its curvature is reproduced by including a fifth effective transmon qubit level ( | 4 ) in our numerical simulations (see Appendix E).
Crucial to reproducing this transition in our theoretical simulations is taking into account that the bound-state level structure cannot be defined by a single anharmonicity, i.e., given the anharmonicity Δ = ω 12 − ω 01 of the bound state, the frequency of the fourth level of the bound state is not simply given by 3ω 01 − 3Δ as is expected for a transmon [38] [see transmon level diagram in Fig. 2(a)].
We observe the flip-flop interaction (by a flip-flop interaction between two qubits with basis states |0〉, |1〉, we mean a Hamiltonian proportional to |01〉〈10| + |10〉〈01|) between the two the spatially bound states measuring the avoided crossing in transmission when the bound states are tuned into resonance. As these qubits (while we really have a multilevel systems coupled to the photonic crystal, we refer to them as qubits when one can ignore higher levels) are a fixed distance apart (9 mm) and there is negligible direct capacitive coupling, the strength of the flip-flop interaction will be entirely determined by the overlap of the localized photonic mode of one qubit with the other qubit, controllable here via the qubit frequencies.
In Fig. 3(a), the frequency of one qubit is held constant while the other is tuned through resonance. Measuring transmission at the single-photon level reveals an avoided crossing between the |01 and |10 levels of the coupled dressed bound states. Transmission amplitude of a bound state dims when the bound state and bare qubit are near resonance [see Fig. 3(a) inset]. From this plot, we can extract a resonant bound-state-bound-state interaction of 120 MHz for a 7.73-GHz bare qubit frequency. In comparison, Figs. 3(b)-3(d) show reduced interaction strength when both qubits are further detuned from the band edge, 6.125 GHz, 6.75 GHz, and 7.625 GHz, respectively, for the fixed qubit frequency.
To characterize this aspect of the two bound-state interaction, we map the magnitude of the avoided crossing as a function of detuning. In Fig. 3(g), the qubits are maintained on resonance with one another while being simultaneously tuned through the band gap. Theoretical modeling [Figs. 3(g) and 8(e)] shows experimental data to be consistent with localized photonic states and with interaction via wave-function overlap. In the limit where the qubit is deep in the gap, the Markovian approximation holds as is evident in the probability that the bare qubits are in the excited state (obtained from the hopping model) in Fig. 3(g) [see Fig. 8(c) for the photonic probability of the bound states]. Here, the localized mode and flip-flop interaction both have the same distance dependence e −x/L where L = a α/δ, is the localization δ is the detuning between the bare qubit and the band edge, a is the unit cell size, and the band-edge dispersion is ω k = ω 0 + αa 2 (k − k 0 ) 2 (see Appendix D 2). The corresponding flip-flop interaction Hamiltonian is When the bare qubit frequencies are near the band edge the probability that the bare qubits are in the excited state decreases [ Fig. 3(g)] and we have non-Markovian behavior (see Appendix D 3 for a detailed on the breakdown of Markovian behavior). While our experiment studies steady states, in other settings, where an initial state evolves in the absence of an input field, non-Markovianity can lead to the preservation of entanglement during time evolution [6,[39][40][41]. It would be interesting to study in future work whether there is also entanglement in steady state and how to best measure it in our system.
We now turn to the qubit nature of the bound states when the bare qubit frequencies are resonant. The qubit part of the wave function of these two bound states (obtained from the effective Hamiltonian of the system) is (approximately) either symmetric (|0〉|1〉 + |1〉|0〉) or antisymmetric (|0〉|1〉 − |1〉|0〉) We theoretically predict that the higher (lower) frequency bound state at resonance in Fig. 3(a) is symmetric (antisymmetric) (see Appendix D 2). This turns out to affect transmission through the system. Intuitively, the antisymmetric state is expected to be dimmed as the exponential parts of the localized photonic states cancel each other out; and hence the linewidth, which is proportional to the wave function at the end of the photonic crystal, is smaller. However, because the band edge is not at zero momentum in our system, it turns out the symmetric state is actually dimmed and has a smaller linewidth, as we prove in Appendix D 2. In Fig. 3(e), we see that the bound states at the same transmission frequency (with different bare qubit frequencies) have drastically different linewidths with the higher-frequency bound state having a smaller linewidth, consistent with our numerical simulations [Figs. 3(f) and 8(d)]. This provides some indirect experimental evidence that the qubit part of the higher (lower) frequency bound-state wave function is indeed symmetric (antisymmetric).
To further study tunable on-site interaction, we probe the interacting bound states beyond the one-excitation manifold using spectroscopic measurements [see Fig. 4(a)]. Similar to spectroscopy of qubits in cavities, we can use transmission at the band edge to help detect bound-state transitions, a technique that provides sharper contrast compared to transmission measurement for the more highly localized bound states and allows detection of higherdressed transitions, such as the transition between |0〉 and |2〉 driven by two photons of frequency ω 02 =2.
With this technique we detect interaction between |02〉, |20〉, and |11〉 of the coupled bound states, observed as avoided level crossings. In addition to the single-photon exchange interaction between |02〉 |20〉 |11〉 and [26], remarkably we measure the two-photon virtual interaction between |20〉 and |02〉, despite the fact that this process is fourth order in coupling g (see Appendix F 2). This two-photon interaction shows consistent dependence on detuning: increasing in strength (from 0 MHz to over 10 MHz) as the bound states shift towards the band edge and the states become more delocalized [see inset of Fig. 4(a)]. Numerical simulations [ Fig. 4(b)] are consistent with experimental data and capture the relative magnitudes of interaction between levels as well as frequency dependence on coupling strengths. Observation of this small interaction highlights the overall strength of interbound-state coupling possible via overlap alone.
The widely tunable on-site and interbound-state interactions demonstrated with this device and consistent theoretical simulations are promising benchmarks towards realizing larger, more complex systems of bound states. Examples of these systems include spin models with both local and long-range interactions, which arise when the bare qubit frequencies are deep in the band gap [17], and complicated multiqubit-photon bound states, which arise when the bare qubit frequencies are in the passband (see, e.g., Ref. [7]). Beyond stepped impedance coplanar waveguides, there are undoubtedly numerous ways to realize superconducting microwave photonic crystals, including lumped element or Josephson junction-based designs, that are equally compatible with superconducting qubits. Regardless of the platform, behavior of bound states due to qubit-band edge coupling will mirror the behavior characterized in this work-elevating this platform above any single experimental design choice.
While the bound states were centered in neighboring unit cells in this device, this is not a limitation or requirement for future experiments as the range of localization can be accordingly set via the basic crystal parameters, as seen by comparing bound-state linewidths measured here with those reported previously [13]. Therefore, one can realize a one-dimensional chain of bound states in a moderately sized photonic crystal, where individual control over the qubits would allow dialing up or down long-distance interaction between sets of qubits.
In this work the interactions were in situ tunable via qubit frequency (DC flux), a static quantity on the timescale of the bound-state lifetime. Dynamically controllable interactions would introduce an additional tool for designing and manipulating spin Hamiltonians [17]. One method for realizing this type of fast timescale control is flux pumping, a technique involving microwave frequency modulation of the qubit frequency along the flux bias line [42][43][44][45]. Another potential pathway would use an auxiliary microwave field through the crystal itself. Here, the qubits could be maintained on resonance deep in the band gap such that there is minimal interaction via bound-state overlap. A single rf control tone can be turned on to induce a transition close to the passband, thus redressing the bound states into new, effective bound states with interaction strength depending on properties of the microwave drive. The addition of several drives or precisely shaped microwave pulses (made possible by commercial high-speed arbitrary microwave waveform generators [46]) promise not only changing interaction strength but also modifying the shape of the interaction itself -from an exponential to a sum of exponentials-leading to a wide range of possibilities including power-law-decaying interactions [17]. These supplementary forms of tunable control would expand the ability of the qubit-photonic crystal platform to realize a broader class of tunable spin models.

APPENDIX A:: MOTIVATION FOR USING PHOTONIC CRYSTALS
In this Appendix, we discuss our motivation for using photonic crystals and their advantages over other devices that might realize bound states. To realize a dressed bound state between a qubit and a band edge, we must couple a qubit to a photonic band edge, where a band edge consists of photonic states corresponding to the transition between a stop band and a passband that are "slow light" due to reduced group velocity. Band gaps and passbands are not unique to photonic crystals-we come across them in numerous other structures such as waveguides (near cutoff frequency) or aperiodic filters.
To motivate the benefit of photonic crystals, let us consider using an aperiodic filter instead. Filters are ubiquitous in the microwave domain with many established design methods that trade off optimizing various parameters such as roll-off or passband ripple. Stepped impedance filters are a standard model for implementing filters, where each impedance step is chosen precisely to meet filter design constraints with no periodicity requirement [30]. Because of this, aperiodic filters may also be more sensitive to fabrication errors than periodic filters. The next requirement is to couple a qubit. However, the lack of periodicity transfers also to the electric field distribution (no Bloch modes), and so we must numerically calculate the optimal location to place the qubit and recalculate for every modification of filter design. While this seems feasible for a single qubit, the lack of Bloch modes is highly problematic from a scalability standpoint as it is not guaranteed that we can couple multiple qubits (or even just two), with near identical coupling strengths, to that band edge.
Furthermore, from a theoretical perspective, periodic crystals are simpler, cleaner structures that are described in the infinite limit by dispersion relations in momentum space, providing useful insight for predicting system behavior. Therefore, while it may be possible to realize dressed bound states with a qubit in an aperiodic structure, the benefits from using a periodic structure far outweigh the likely larger device footprint.
A photonic crystal is an electromagnetic structure formed by periodic modulation of the dielectric constant. This results in dispersion relations characterized by energy bandsalternating band gaps ("forbidden energies") and passbands (continuous density of states). The electric field distribution is characterized by Bloch modes, allowing for the position of qubits at locations that optimize coupling to the band edge [8]. Engineering band edges via finite-size photonic crystals has been demonstrated in many systems, including nanophotonic structures [47,48] and superconducting microwave coplanar waveguides [13,49], and tremendous progress towards integration with ultracold atoms and superconducting qubits, respectively, has been reported [12,13]. There are other ways to create superconducting microwave crystals, including using Josephson junction arrays or lumped element circuits.
Our approach to creating an effective 1D microwave photonic crystal is periodically alternating sections of high and low impedance coplanar waveguides (CPWs). With CPWs, the impedance can be easily changed via the center-pin width and the center-pin to ground plane (the gap) distance.
We define a unit cell as a high impedance section of length L hi and impedance Z hi sandwiched between two sections of low impedance of lengths L lo /2 and impedances Z lo (for symmetry purposes). With a periodic modulation, there are naturally many gaps in the band structure. We chose to more strongly couple the qubit to the second band edge, rather than the first, because it has a smoother passband while still having a steep roll-off.

Crystal simulation and parameters for implementation
We can use the unit cell to calculate the band structure or dispersion relation for an infinite crystal. While we can never make an infinite crystal, calculating the dispersion relation is a very useful starting point and gives us insight into crystal parameters such as band curvature. To a good approximation, the phase velocities in the high and low impedance CPW sections are effectively equal (v p;high ≈ v p;low ≈ v p ≈ 1.248 × 10 8 m/s). This yields where k is the momentum and ω k is the dispersion. We determine the band structure dependence on L lo and L hi for Z lo ≈ 25 Ω and Z hi ≈ 125 Ω.
We look to optimizing the trade-offs across four parameters: the frequency of the band edge, the width of the band gap, the width of the passband, and the curvature of the band. From these simulations, we see that small changes in unit cell impedances do not lead to significant changes in the band dispersion.
For comparison, Figs. 5(a) and 5(b) show the simulated dispersion for unit cell parameters in Liu et al. [13] and the present paper, respectively. The unit cell for the present paper was chosen to have a flatter band dispersion (analogous to effective mass), α, so as to realize more localized bound states.
While the dispersion relation assumes an infinite system, crystals of small, finite length have been shown to realize well-resolved gaps in dispersion where transmission is suppressed and bands where transmission is unimpeded. From an experimental perspective, we use transmission-based measurements to probe the states.
Transfer matrices are a convenient and accurate method for bare crystal simulation that incorporates both the exact number of cells and boundary conditions [31]. A convenient metric for comparison is transmission across the device, S 21 . In Figs. 5(c) and 5(d), comparison of the measured S 21 of the bare crystal (taken at powers high enough to saturate qubit effects) and the calculated S 21 from transfer matrices shows agreement [13].

Realistic parameters and experimental aside
The device is fabricated on an approximately 430-μm-thick c-plane sapphire substrate, on which we sputter approximately 200 nm of niobium. The photonic crystal is patterned using a direct-write laser writer, followed by a dry SF 6 reactive ion etch to transfer the pattern. To ensure our waveguides are reasonably sized and can be reliably fabricated using photolithography, we are limited to impedances between approximately 25 and 125 Ω.
For this device, one unit cell consists of a high impedance section (Z hi = 124 Ω, L hi = 7.8 mm) and a low impedance section (Z lo = 25 Ω, L lo =1.2mm). Impedance estimates are obtained by fitting the measured spectrum with transfer matrix simulations; the dispersion is robust to small impedance deviations in the fabricated sample. We fit 16 unit cells on a 10 mm × 10 mm sapphire chip. While more unit cells could fit on the chip, we see experimentally that we must wire-bond extensively on-chip (to connect ground planes) to create clean band gaps and passbands and included this requirement into the design. By switching to air or dielectric-supported bridges in the future we would be able to shrink the device footprint or include more unit cells in the same-size device.
To integrate with the standard measurement setup and qubit parameters, we choose unit cell lengths such that the second band edge falls between 7 and 8 GHz. The frequency of the band edge was designed to be at 7.8 GHz to match with the traveling-wave parametric amplifier (TWPA, dispersion feature around 8.3 GHz) to provide good amplification for frequencies in the vicinity of the band edge as well as deep in the band gap. The width of the band gap is from 4.75 to 7.8 GHz and needs to be large enough such that we have a wide range to tune the qubit frequency over, which would in turn result in a wide range of accessible localization lengths. The width of the (second) passband needs to be sufficiently large such that we can ignore the third band edge and we can make the approximation that the curvature of the band (near the band edge) is quadratic. Finally, as the curvature of the band plays a role in determining the localization length, to make the states more localized, we chose a shallower band curvature [compare red plots in Figs. 5(a) and 5(b)].
In future iterations, one may consider altering boundary conditions specific to the desired application. For example, methods for impedance matching such as tapers or quarter-wave transformers [30] are straightforward additions that will modify impedance matching at specific frequencies (such as at a band edge). These options were not pursued in this work as we wanted to study bound-state properties across a range of frequencies. If one were interested in only specific frequencies or ranges, then this is a promising improvement. The device is symmetric, so one end is chosen arbitrarily to serve as the input. As detection of the bound state is due to scattering, one may consider modifying the symmetry of the device or detecting signal from both ends of the crystal to improve collection efficiency.
We resort to two-tone spectroscopy over direct transmission measurement to detect the bound state when the qubit is deep in the band gap due to strong localization and poor signal-to-noise ratio (SNR). Choosing a less shallow band curvature will make it possible to see the bound state all the way through the gap (see Liu et al. [13]); however, this is a tradeoff as the bound-state linewidths will increase accordingly. Thus, choosing curvature and crystal parameters such that the linewidth remains as narrow as possible without internal Q or loss effects is essential.

APPENDIX B:: ADDING IN QUBITS
In this Appendix, we discuss the experimental details of coupling the transmons to the photonic crystal. In this device, we capacitively couple a transmon qubit to each of the two central unit cells of the 16-cell crystal. A transmon's anharmonic level structure is due to the nonlinear inductance of Josephson junctions and allows for selective addressing of energy level transitions. However, it is important to emphasize that our realization of a qubit does have higher energy levels set by transmon geometry, unlike the standard theoretical qubit which is synonymous with a two-level system. These higher levels are also coupled to the band edge and therefore accounting for these levels becomes important when looking beyond the first excitation sector.
Our qubits are fully patterned with a 125-kV e-beam writer, with bridge-free junctions, and are made of evaporated aluminum. The qubits are designed to have a target charging energy (from electromagnetic simulation) of approximately 450 MHz, and emphasis is placed on maximizing coupling between the qubit and photonic crystal without significantly altering the unit cell. Finally, although identical in design, in reality the qubits will differ due to fabrication variability. However, as coupling to the waveguide is designed to be the dominant decay channel for the qubits, bound states are robust to small variation in other parameters.
We place qubits in adjacent unit cells (9-mm separation) at the center of the crystal such that we can see significant change in interaction strength with detuning. Each qubit has a local flux bias line for independent control-DC cross talk is corrected for through standard calibration. These lines are low pass filtered; however, as the dominant decay of the qubit is via the crystal, this is not expected to be a limiting factor in coherence.
To determine where to place the qubits within the unit cell such that they maximally couple to the desired band, we must calculate the electric field distribution within the unit cell, distribution determined by the Bloch modes for the crystal [8].
For these crystal parameters, maximally coupling to the second band edge (and minimally to the first band edge) corresponds to placing the qubit in the center of long high-impedance section. Other locations within a unit cell change coupling to each band edge, which is a potentially interesting regime for future experiments. However, here we are interested in reducing the effect of the other band as much as possible. We cannot completely eliminate this coupling-experimentally we can still detect a drop in transmission in the lower passband when the qubit is resonant but this change is orders of magnitude smaller than in the second band.

APPENDIX C:: MODEL, TRANSMISSION, AND FITTING OF PARAMETERS
In this Appendix, we introduce the effective Hamiltonian for our system and discuss the fitting of various parameters of the Hamiltonian in detail.

Hamiltonian
The Hamiltonian for the one-dimensional photonic crystal is given by H c = ∑ k ω k a k † a k , where a k † creates a bosonic excitation with momentum k, and ω k is the dispersion relation of the second band of the photonic crystal, which is discussed in Sec. A 1. Here, we ignore the other bands of the photonic crystal as the qubits couple predominantly to the second band.
Fourier transforming a k † = a/L∑ j = 1 N a j † e ikx j , where L is the system size, a is the unit cell size, N is the number of unit cells, x j = aj, and k = [(2π)/L]n, where n is an integer) gives a hopping model with periodic boundary conditions, Here, k is a dimensionless integration variable, and we use the fact that J i,j = J j,i since ω k = ω −k . In Eq. (C1) we take N → ∞, as we only know ω k in that limit (see Sec. A 1). To model 0014GHz. Unfortunately, we are unaware of an exact analytical solution for J i,j for our system. In our numerical simulations, we keep hopping terms up to J 5 . We calculate the hopping parameters for a given set of photonic crystal parameters. A different choice of photonic crystal parameters would have given a different set of hopping parameters. As such, these parameters should only be understood as estimates. We briefly comment on the change in theory parameters that arises from using different photonic crystal parameters at the end of this section.
The Hamiltonian for the isolated transmon qubits is given by Here, i labels the transmon qubit, n labels the level of the transmon qubit and ω 0n;i are the bare energy levels of the transmon qubits. In our simulation, the number of transmon qubit levels is truncated at five (i.e., |0〉 through |4〉). For our experiment, ω 00;i =0, ω 02;i =2ω 01;i −Δ i , ω 03;i =3ω 01;i −3Δ i , and ω 04;i = 4ω 01;i − 6Δ i , where Δ i is the bare anharmonicity of the ith transmon qubit.
We now turn to the coupling of the transmon qubits to the photonic crystal. To an excellent approximation, the coupling between the transmon qubit and the photonic crystal takes place within a single unit cell (see Appendix D 2 for justification of this statement). Thus, in the rotating-wave approximation, we can write the coupling term of the Hamiltonian as where z i labels the position of the two transmon qubits. In our system, the transmon qubits are on neighboring unit cells, i.e., z 1 = (N/2) and z 2 = (N/2) + 1, and the coupling for each transmon qubit, g i , is different due to small experimental imperfections. The total Hamiltonian of the system is then Finally, we note that this Hamiltonian conserves total excitation number.

Transmission methods
We now discuss the two theoretical methods we use to calculate transmission in the linear drive regime. Neither method relies on a weak coupling approximation between the transmon qubits and the photonic crystal. The first method involves treating the system as an open quantum system, with loss on each site (that is site dependent), subject to a weak drive on the first site. The largest loss terms are at the ends of the one-dimensional photonic crystal. The system can then be described by the following effective non-Hermitian Hamiltonian (in the rotating frame) with driving frequency ω d and strength ϵ, where κ q is the qubit half-width, κ 0 is a uniform contribution to photonic half-width, and κ is a decay rate on the first and last sites. While there are certainly other forms of loss (such as nonuniform loss on each site), our goal is to reproduce the key features (e.g., the locations of the bound state and of the transmission dip, as well as the linewidth of the bound state) using as few parameters as possible. The equations of motion for the quantum operators are ∂a y ∂t = − i ∑ j J y, j − ω d δ y, j a j + ϵδ y, 1 + ∑ i = 1, 2 δ y, z i g i |0 i 1| i − κ 0 + κ δ y, 1 + δ y, N a y , We omit vacuum Langevin noise from the equations of motion as this noise does not affect our calculations. Solving for the steady state of a N in the weak drive limit (〈|1〉 i 〈1|i〉ss = 0) yields the transmission. More specifically, S 21 ∝ 〈a N 〉 ss /ϵ.
The second method we use was introduced by Biondi et al. [29]. Here, we treat the system (which is taken to be the photonic crystal and qubits) as being connected to waveguides with linear dispersions, with velocity v g , at the ends of the photonic crystal. In this method, we take κ q = κ 0 = κ = 0, so that the single-excitation transmission through the system can be expressed in terms of eigenvalues and eigenvectors of the system described by H tot [Eq.
(C4)] [29]. More explicitly, the transmission for a given frequency ω is given by |T(ω)| 2 , where Here, V n l, r = v g g w 0 a 1, N n , Ωn and |n〉 are the eigenvalues and eigenvectors of H tot in the single-excitation sector, and g w is the coupling between the waveguide and the photonic crystal. Intuitively, transmission occurs when the single-excitation eigenstates have the probability of the photon being on both ends of the photonic crystal. Near the bound-state energy E 1B , we can write the transmission as Here, we have assumed that the wave function is symmetric, i.e., V n l = V n r = V . We see that the linewidth (in the limit) is proportional to the wave function at the edge of the system V, as expected.

Fitting of parameters
In this section, we discuss how we fit various parameters of the total Hamiltonian. The unknown parameters include g i and Δ, and, for the first method, also κ 0 , κ q , and κ. Furthermore, ω 01;i is tunable but its value is unknown, and the transmission dip (visible when the bare qubit frequency is in the passband) does not, in general, occur at the bare qubit frequency unless hopping amplitudes J i,j beyond nearest neighbor are negligible.
The first parameters we determine are κ and κ 0 (using the first transmission method). We turn off the coupling of the qubits to the photonic crystal (in the experiment, this is accomplished by saturating the qubits). We set κ 0 and κ q to zero and fit κ. Given that the largest losses occur at the ends of the photonic crystal, fitting κ first is reasonable. κ controls the linewidth of the photonic modes and the transmission amplitude difference between the transmission dips and peaks in the passband. We find that a reasonable estimate for κ [for the experimental data in Fig. 5(d)] is 1 GHz, although any κ in the range of 0.5 to 1.5 GHz also gives a reasonable fit. We then turn on κ 0 , which further reduces the transmission amplitude difference between the transmission dips and peaks in the passband and lowers the transmission peak of the lowest photonic mode. We estimate that κ 0 = 4 MHz (although any κ in the range of 3 to 5 MHz also fits the data well). Figure 5(d) shows that simulated transmission is in good agreement with the experimental data. Given that there are other losses in the system that we have not included, these numbers should be understood as estimates.
We now turn to determining g i . While κ q is set to zero for now, we find that varying it or the other loss parameters does not noticeably change the frequency of the bound state or transmission dip, thus making our estimate of the coupling strength independent of the decay parameters. To begin, we detune the qubit at site N/2 far away from the passband (in the experiment, the detuned qubit is at 4.5 GHz) and then fix the other (i.e., second) bare qubit frequency [50]. Transmission is then calculated as a function of driving frequency. We find that g 2 = 0.55 GHz and ω 01;2 = 7.9875 GHz match the experimental data well when the bound state is at 7.605 GHz, as seen in Fig. 6(a). Calculating transmission when the first qubit frequency is near the passband and the second qubit frequency is detuned and comparing it to experimental data, we find g 1 = 0.505 GHz [50]. To make sure we have chosen suitable coupling strengths, we tune the bare qubit frequency (the detuned bare qubit frequency is kept fixed). If we have chosen the correct parameters, we should accurately capture the locations of the bound state and the transmission dip for different bare qubit frequencies, while keeping the other parameters fixed. Indeed, as seen in Fig. 6(b), we find this is the case for the chosen coupling strengths.
Our next goal is to obtain an estimate for κ q . To do so, we increase κ q , which increases the linewidth of both the transmission dip and the bound state, until the linewidth of the bound state matches the experimental value well for a fixed bare qubit frequency (we note increasing κ 0 also increases the linewidth of the bound state, however κ 0 is already fixed). We find that κ q = 0.5 MHz accomplishes this task for ω 01;2 = 7.9875 GHz. To make sure we have chosen a suitable qubit half-width, we again tune the second bare qubit frequency (while keeping the detuned bare qubit frequency fixed). If we have chosen a reasonable qubit half-width, we should be able to accurately estimate the linewidth of the bound state for different bare qubit frequencies (while keeping all other parameters fixed). Figure 1(b) of the main text shows that our estimate of κ q is reasonable. Before moving on, we briefly comment on the second transmission method. Figure 6(c) shows theoretical data from both simulations for g 2 = 0.55 GHz and ω 01;2 = 7.9875 GHz.
The locations of the bound state and transmission dip occur at the same spot for both methods. The key difference is that the second method does not accurately predict the magnitude of the linewidth of the bound state as it assumes κ q = κ 0 = κ = 0 (but instead has coherent coupling of the crystal to the waveguide). In these simulations, we have taken g w = 2 GHz, as that fits the data when the qubits are saturated well (not shown). We note that any value of g w in the range of 1.5 to 2.5 GHz also gives a reasonable fit to the saturated qubit data and that the frequencies of the bound state and the transmission dip are not sensitive to g w . Simulated transmission data presented in the main text is from method one.
We now fit the last remaining parameters, Δ 1 and Δ 2 , which are the bare transmon anharmonicities. We first diagonalize H tot in the two-excitation sector for fixed ω 01;i (with the qubit on site N/2 detuned far away). The theoretical prediction for the dressed anharmonicity of the qubit on site N/2 + 1 is found by taking the lowest eigenvalue of H tot in the two-excitation sector and subtracting two times the lowest single-excitation eigenvalue. We vary Δ 2 until the theoretically predicted dressed anharmonicity of the second qubit matches the experimentally measured dressed anharmonicity for a given bare qubit frequency (we choose our bare qubit frequency such that the single-particle bound state is at 7 GHz). In doing so, we find, to a good approximation, that Δ 2 = 0.365 GHz. We then vary the bare qubit frequency and make sure the theoretically predicted dressed anharmonicity is still consistent with the experimentally measured value for different bare qubit frequencies.
We find excellent agreement for a wide range of bare qubit values as shown in Figs. 2(a) and 2(b) in the main text. Finally, we choose the value of Δ 1 to give the best fit to the magnitude of the two-photon avoided crossing [see inset of Fig. 4(a) in the main text]. In doing so, we find that Δ 1 = 0.365 GHz fits the experimental data reasonably well.
Before we close this section, we estimate errors in our parameters. To begin, we estimate what change in hopping parameters (we call these different hopping parameters J′) we would get if we choose Z high = 123.5 Ω instead of 124 Ω. Both of these choices fit the experimental data well in transmission matrix simulations. This choice of Z high gives J ′ 0 =9.331GHz, J′ 1 =0.7308GHz, J′ 2 − .0345GHz, J′ 3 =0179 J′ 4 =−0035GHz, and J ′ 5 =0.0014GHz. We see that the ratio of these hopping parameters to the previous set is not less than 0.97 for any term, so we expect the other parameters in our model will not differ by more than 5% from their given predictions. We also expect this to hold true if one uses any reasonable set of photonic crystal parameters. To test this, we estimate g 2 and the range of decay parameters for the second set of hopping parameters. We find g 2 = 0.555 GHz and that the same range of decay parameters fit the data well, consistent with our expectation.

APPENDIX D:: BOUND-STATE FUNDAMENTALS
In this Appendix, we discuss the theory of bound states. Strong light-matter interaction between atoms and slow-light structures, ones with vanishing or significantly reduced group velocity such as photonic band edges in photonic crystals, has been an area of ongoing interest both in theory and recent experiment. The principal interest behind this study is the localized bound photonic state that forms around the atom. This bound state has an exponentially decaying photonic envelope that tunes with detuning of the qubit transition from the band edge (see Fig. 7 In a finite-size system, these localized states overlap with the ends of the crystal, thus facilitating single-photon transport across the crystal at the bound-state frequency and providing an avenue to probe these states. This tunable photonic interaction mechanism provides a platform for simulation of many-body quantum optics in one-dimensional systems, distinct from cavity or waveguide quantum electrodynamics [see Figs. 7(c) and 7(d)].
In Figs. 1(e) and 1(f), we measure (simulate) S 21 at low power to track the bound state as a function of qubit frequency. The bound state shows up as a Lorentzian peak in transmission. We detect the change in wave-function overlap as a change in bound-state linewidth-where linewidth increases with localization length [see Fig. 1(c) in [Ref. 13]]. As discussed in Sec. A, the localization of the photonic wave function is determined predominantly by the strength of the coupling, properties of the band edge, and the frequency detuning between the atom and the band edge. The ability to tune the localization with small-sized crystals shows the versatility of this platform.

Bound-state linewidth dependence on detuning
The linewidth of the bound state is set by the amplitude of the exponentially decaying photonic wave function at the ends of the crystal (ignoring other forms of loss). The envelope amplitude decays as ~e −x/L , where x is the distance from the qubit location and L is the localization length. For simplicity, we consider the case of a single qubit at the center of the crystal (total length d) such that the envelope is symmetric. This results in a linewidth proportional to e −d/2L . Of course, this approximation is only valid in the limit where e −d/2L ≪ 1, meaning the bound state is sufficiently localized compared to the finite length of the crystal.

Single-photon bound states: Exact solution
In this section, we discuss the theory of atom-photon bound states in the single-excitation sector for an infinite photonic crystal coupled to two qubits. While our system is of course finite, these results provide an intuitive understanding of our system. We note that a similar calculation for two qubits was first done, to the best of our knowledge, in Ref. [6] and then later in Ref. [7] for the case when the two qubits have equal coupling strengths and equal qubit frequencies. The results below generalize that work to unequal coupling strengths and unequal qubit frequencies. To begin, we first write the Fourier transformed Hamiltonian from a more microscopic point of view (ignoring decay and ignoring terms that do not affect the single-excitation bound states), To make analytical progress, we assume that the dispersion relation is ω k = ω 0 + αa 2 [k ∓ (π/a)] 2 which is valid around k = ±π/a. While we chose a quadratic dispersion, these results are qualitatively similar for a cosine dispersion. Here, we also restore the momentum dependence of the coupling strength g k, i = ω k /2ℏd i u k r i , where d i is the matrix element between |0〉 and |1〉 states for the ith qubit and u k (r i ) is the spatial function for the kth mode evaluated at the ith qubit position, r i (see Ref. [6]). As ω 0 is much larger than the bandwidth, we can approximate ω k as ω 0 . Furthermore, for our system, u k (r i ) for the second band is nearly independent of momentum as shown in Ref. [13]. Thus, to an excellent approximation, g k,i is constant, which we define as g k,i = g i , and hence the coupling is local in real space. Finally, we note that u k (r i ) for the first band vanishes when the qubit is in the middle of the high-impedance section, thus we can ignore coupling to the first band to an excellent approximation. The most general single-excitation wave function is |ψ 1B = sin ϕ a 1 | 1 1 | 0 2 | 0 + a 2 | 0 1 | 1 2 | 0 Here, the basis states are labeled as |qubit one〉 1 |qubit two〉 2 |photon〉, and the conditions, The sums are evaluated as follows: Shifting the integrals by π/a, making the integrals dimensionless, and extending the limits to infinity gives Here, we have assumed that E 1B < ω 0 . Following the same steps for the remaining integral Here, we use the fact that (z 1 − z 2 ) is an integer. We then have E 1B a 1 = ω 01; 1 a 1 − 1 2 α ω 0 − E 1B × a 1 g 1 2 + a 2 g 1 g 2 cos π z 1 − z 2 e − ω 0 − E 1B /α z 1 − z 2 . (D10) Repeating these steps for a 2 gives the following equation for the bound-state energy, where we use the fact that (z 1 − z 2 ) is an integer again. We note that when the qubits are infinitely far away from each other, we recover the well-known bound-state energy for a single qubit (see, e.g., Ref. [7] or Ref. [13]), Generically, Eq. (D11) yields one or two bound states depending on the coupling strength, qubit frequencies, distance between qubits, and α [7].
To illustrate this point, we consider the case when g 1 = g 2 = g and ω 01;1 = ω 01;2 = ω 01 [which is relevant to the case in Fig. 3(g) of the main text]. While the coupling strengths are not exactly equal in our experimental system, it is a decent approximation to consider them equal. In this case, we expect a symmetric and an antisymmetric solution, i.e., where + is for the symmetric state and − is for the antisymmetric state and Σ ± δE 1B is the self-energy. The condition for the presence of a bound state, as derived in Ref. [7] is − ω 01 − ω 0 > Σ ± (0).
We explicitly consider the experimentally relevant case when |z 1 − z 2 | is odd. In this case, we have Σ + (0) = − g 2 /2α z 1 − z 2 and Σ − (0) = − ∞. For the antisymmetric state, the condition is always satisfied, while for the symmetric state, we only have a bound state if g > 2α ω 01 − ω 0 / z 1 − z 2 . We now apply this formalism to our experimental system. For our experimental system, α = 1.155 GHz and |z 1 − z 2 | = 1. We also take g to be the average of the two coupling strengths determined in the previous section, i.e., g = (g 1 + g 2 )/2 = 0.5275 GHz. Using these numbers, we estimate that we have two bound states for ω 01 − ω 0 < 120 MHz. We remind the reader that this result is only an estimate, as our experimental system is finite and has unequal coupling strengths.
Finally, we investigate the wave function. The ratio of a 1 /a 2 , which we define to be η, is given by a 1 a 2 = η = − 1 2 α ω 0 − E 1B × g 1 g 2 cos π z 1 − z 2 e − ω 0 − E 1B /α z 1 − z 2 E 1B − ω 01; 1 + We see that the photon is exponentially localized around the qubits with localization length a α/ ω 0 − E 1B as in the case of a single qubit (see Fig. 7). We remind the reader that we can have two different bound-state energies (one for the symmetric bound state and one for the antisymmetric one), thus two different localization lengths. In other words, the linewidths of the bound states will, in general, be different, with the bound state closer to the band edge having a larger linewidth. In our system, the symmetric bound state, whenever it exists, is closer to the band edge and thus has a larger linewidth.
Finally, we remark that when g 1 ≈g 2 =g and ω 01;1 ≈ω 01;2 =ω 01 , the symmetric and antisymmetric nature of the bound states can affect transport. Consider the case when the qubits are on the same lattice site. In this limit, the photonic contribution of the antisymmetric wave function will vanish and hence cannot be detected via transport [see Eq.
(D17)]. When the qubits are not on the same site, there will be reduced transmission depending on the separation between the qubits and whether or not the state is symmetric or antisymmetric as can be seen explicitly from Eq. (D17). The fact that this effect depends on the separation between qubits can be traced back to the fact that the band edge is at k = ±π (if the band edge was at k = 0, only transport for the antisymmetric state would be reduced).
To obtain an estimate for how much transmission will be reduced for the symmetric state, we plot the ratio of linewidths, i.e., |c N , for experimentally relevant parameters (|z 1 − z 2 | = 1, N − z 1 = 7, and N − z 2 = 8) in Fig. 8(d). Indeed, we see that the symmetric state has a much smaller linewidth and hence less transmission. This is consistent with the experimental results presented in Fig. 3(e).

Single-photon bound states: Born-Markov solution
We now briefly compare these exact results to the Born-Markov (BM) approximation. For simplicity, we restrict ourselves to the case when the qubit frequencies and coupling strengths are the same. Using a second-order Schrieffer-Wolff transformation to eliminate the high-energy subspace gives the following effective Hamiltonian for the two qubits, In the limit where the qubits are infinitely far apart, we recover the standard expression for the dressed qubit frequency, ω 01 ′ = ω 01 − g 2 2 α ω 0 − ω 01 . (D21) In Fig. 8, we compare the results obtained in the Born-Markov approximation to the exact analytical results. Figure 8(a) shows these results for a single qubit and Fig. 8(b) shows these results for two qubits with |z 1 − z 2 | = 1. As expected, the Born-Markov approximation is a good approximation when the qubit frequency is away from the band edge. Closer to the band edge, the Born-Markov approximation fails, particularly for the lower-energy (i.e., antisymmetric) state. Furthermore, by comparing the blue curve in Fig. 8(a) to the experimentally measured bound state in Fig. 1(e), we clearly see that the experiment is not well described by the Born-Markov approximation. We note that the higher energy red line in Fig. 8(b) ends abruptly at ω 01 ≈ 7.920 GHz as there is only one bound state for ω 01 > ω 0 + 0.120 GHz.
We now analytically show that the Born-Markov approximation is excellent for one of the dressed states close to the band edge. We begin by expanding Eq. (D13) (for the symmetric case, when |z 1 − z 2 | = 1) in the limit that ω 0 − E 1B ≪ α, i.e., when the bound-state energy is close the band edge. This gives Now comparing to the Born-Markov solution for E 1 in Eq. (D19) around ω 0 ≈ ω 01 , we have E 1 ≈ ω 01 − (g 2 =2α). We thus see that one of the dressed states (the symmetric one) is well captured by the Born-Markov approximation, while the other is not as seen in Fig. 8(b).

Single-photon transport via the bound state
As discussed, a bound state mediates transport across the crystal, at the otherwise forbidden frequencies inside the band gap, via the overlap of the photonic mode with the ends of the crystal. However, unlike a cavity mode which accommodates many photons (of the same frequency) due to its harmonic nature, the bound state inherits an anharmonic level structure from the qubit. This will result in single-photon, blockaded transport. For a definitive confirmation, we measure the second-order autocorrelation of the transmitted component of a weak, resonant, continuous drive.
In Fig. 2(c), we plot the emission spectrum of the resonantly driven bound state for low drive amplitudes. Here, we see the familiar Mollow triplet structure featuring sidebands that are linearly displaced from the center peak with increasing drive amplitude. We measure second-order autocorrelation [ Fig. 2(c) inset] for a drive amplitude below the threshold for incoherent triplet emission such that the qubit is not saturated by the drive. This measurement (see [36,[51][52][53][54] for concept) was made possible by a TWPA (MIT Lincoln Labs) to improve SNR and a GPU (CUDA-Matlab) for significant computational speed-up.

APPENDIX E:: EMISSION THEORY
In this Appendix, we discuss the theoretical modeling of the emission spectrum of a resonantly driven bound state. Unfortunately, due to the large dimension of the Hilbert space, we found that a direct approach of calculating the emission spectrum using the master equation for our sixteen-unit-cell system is not numerically feasible. Instead, we diagonalize the Hamiltonian and investigate energy differences. While this approach does not predict the widths and the driving-strength-dependent intensities of the sidebands, it does predict the frequencies of the sidebands. The Hamiltonian of a single transmon qubit with a driving frequency ω d equal to the bound-state frequency is given by (in the rotating frame) [ where Ω is the bare Rabi frequency of the drive and is the only unknown parameter. We stress that ω d is not at the bare qubit frequency but at the frequency of the bound state. In the presence of a drive, the excitation number is no longer conserved, and thus, to make progress, one must implement a cut-off. In our numerical simulations, we have implemented a cut-off of five transmon qubit levels and three photons. Diagonalizing the system, we take the differences of the eigenvalues of states corresponding to large occupation of atomic states with no photons. We now compare the results obtained by diagonalizing Eq. (E1) to the results obtained by driving a dressed-transmon qubit (without explicitly including the photonic crystal). The latter approach was used in Ref. [13]. Here, ω 0n are the dressed-transmon qubit frequencies, ω d = ω 01 , and Ω is the Rabi frequency seen by the dressed-transmon qubit. If the exact wave function for the bound state is |ψ = cos θ | 1 | 0 + sin θ∑ k c k a k † | 0 | 0 with ∑ k c k 2 = 1 (here, our basis states are labeled as | qubit〉|photon〉), we have Ω ≈ Ω cos θ, where θ is given by which, for an infinite system, simplifies to When the bound state is approximately at 7.59 GHz, we find cos θ ≈ 0.68 for our finite system. Figure 9 compares the results obtained from Eqs. (E1) and (E2). We see that, upon taking into account the reduction of the matrix element due to the dressing, the two methods agree. We also see [ Fig. 9(b)] that, as expected, as the bare qubit frequency moves deeper into the gap, θ approaches zero and Ω approaches Ω. In Fig. 9, we have used the anharmonicity value predicted by theory. In Fig. 2(c) of the main text, we use the experimental values of anharmonicity, which is given by ω 02 − 2ω d = − 0.11 GHz.
We assume that the feature around 7.22 GHz in Fig. 2 to illustrate which terms in Eq. (F1) are important for a given bare qubit frequency. For example, when the bare qubit frequency is deep in the gap, the ground state of the twoexcitation sector is mostly in the |2〉|0〉 state as seen in Fig. 10(c). Upon increasing the bare qubit frequency (while still in the band gap), the population of |1〉 increases, while the population of |0〉 stays relatively small. This is because two photons must be exchanged to couple the dominate |2〉|0〉 state and any two-photon state and there are no terms in H tot that directly exchange two photons, thus making it a higher-order process. On the other hand, coupling |2〉|0〉 and a i † | 1 | 0 only requires exchanging one photon. When the bare qubit frequency is at or near the band edge, each transmon qubit level in Eq. (F1) contributes significantly to the bound-state wave function.
Unfortunately, we are unaware of an exact solution similar to the one in Sec. D 2. To make analytical progress, we assume f i,j = 0 which is a valid approximation when the bare qubit frequency is deep in the band gap [as seen in Fig. 10(c)]. Solving the eigenvalue equation, H ψ 2B = E 2B ψ 2B , for the following wave function ansatz (in momentum space), We see that the photon is localized around the transmon qubit, consistent with the exact finite-size numerical results seen in Fig. 10(a). We note this ansatz breaks down when the bare qubit frequency is near the passband as the states become more photonic, in which case we can no longer neglect f i,j . Here, we have assumed that 2ω 0 > ω 02 and ω 0 > ω 12 . These conditions are satisfied in the regime where we experimentally observe the two-photon avoided crossing. We are interested in determining how the interaction decays as a function of distance. Unfortunately, we are unaware of how to analytically evaluate the first two terms (the third term can be evaluated exactly). However, the integrand decreases exponentially as a function of k for the first two terms. Thus, it is reasonable to take the integrand to be a constant value (i.e., the integrand value at k = 0) over a small window, δk, about k = 0 and zero otherwise. This small window is taken to be the momentum value for which the argument in the exponential equals one, i.e., δk = 2 1/ z 2 − z 1 2 − 2ω 0 − ω 02 /α . Evaluating the remaining integral gives H |2, 0 | 0, 2 ≈ − 2g 4 |2, 0 0, 2| ω 0 − ω 12 δk 2π e − z 2 − z 1 2ω 0 − ω 02 /α α 2ω 0 − ω 02 1 ω 0 − ω 12 − 1 ω 0 + ω 12 − ω 02 + g 4 | 2, 0 0, 2 | 1 α ω 02 − 2ω 12 e − z 2 − z 1 ω 0 − ω 12 /α ω 0 − ω 12 − e − z 2 − z 1 ω 0 + ω 12 − ω 02 /α ω 0 + ω 12 − ω 02 e − z 2 − z 1 / α ω 0 − ω 12 α ω 0 − ω 12 + H.c. bound-state avoided level crossings for a fixed qubit whose bare frequency is circa 6.125, 6.75, and 7.625 GHz. As qubits are detuned further from the band edge, bound states are more tightly localized, reducing overlap and thus reducing interaction. (e),(f) Transmission when the qubits are on resonance across a range of qubit frequencies in the experiment and the simulation, respectively. The uneven linewidths of the two bound states when they occur at the same frequency suggest they are symmetric (higher-frequency bound state) and antisymmetric (lower-frequency bound state) states (see main text). (g) Bound-state avoided crossing and qubit population (from simulation) as a function of average bound-state frequency. A steady reduction in interaction strength occurs with increasing detuning from the band edge (moving deeper into the band gap) due to increasing localization of the bound states. Hopping model simulation (black) captures this detuning-dependent behavior observed in experiment (red). Near the band edge, both bound states (blue and cyan) have a significant photonic contribution.

FIG. 4. Interaction between two-excitation levels of two bound states.-
(a) Spectroscopic measurement while tuning one bound state through the other (qubit fixed at 7.2 GHz), reveals survival of strong interaction into the two-excitation manifold. Crossed and dotted lines are guides to the eye to discern the levels belonging to the first (|01〉 and | 10〉) and second (|02〉, |20〉, and |11〉) excitation manifolds, respectively. Here, we have labeled our states as |transmon one; transmon two〉. In addition to the |02〉(|20〉) and |11〉 avoided level crossing, we also detect a two-photon virtual interaction between |02〉 and |20〉 white box). This interaction-fourth order in coupling g-manifests itself in avoided level crossings up to and exceeding 20 MHz. For comparison, the |02〉(|20〉)-|11〉 and |01〉-|10〉 interactions are second order in g and thus both are significantly stronger. Inset: two-photon avoided crossing strength versus average bound-state frequency. (b) Numerical simulation for fixed bare qubit frequency of 7.27 GHz, with the one (two) excitation manifold in black (red).