Interacting Qubit-Photon Bound States with Superconducting Circuits

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

nentially localized photonic mode at the qubit position, which together with the qubit forms a qubit-photon dressed bound state [1][2][3][4][5].Photonic crystals are natural avenues to realize these bound states due to their intrinsically tailorable band structure, and characteristic Bloch mode electric field distribution [6] which enables access to strong coupling with qubits [7][8][9][10][11].Bound states in multi-qubit photonic crystal devices are an ideal platform to study many-body quantum optics in one-dimensional systems [5,[12][13][14][15][16][17].Unlike many qubits coupled to a common cavity mode but similar to the case of some optical multimode cavities [18,19], coupling to a band edge creates bound states that intrinsically preserve the spatial organization of qubits, offering the ability to create one-dimensional 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 [15].In this paper, we demonstrate and characterize the underlying, fundamental tunable on-site and inter-bound state interactions in a two-qubit, superconducting microwave photonic crystal device.
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 [11].That work characterized the dependence of localization length on detuning between the qubit and the band edge and further confirmed the existence of the localized state in the bandgap when the bare qubit is in the passband -an unmistakable signature of non-Markovianity (see Supplement).State localization is tunable in situ with frequency through a range determined by device parameters, including qubit-waveguide coupling and band curvature.Compared with previous work, we attain increased localization in this device (Fig. 1b) due mainly to a flatter band dispersion, realized by tailoring the unit cell of the photonic crystal (see Supplement).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 [5], 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 [5].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 [5,11].In our finite system, we observe ∼ 250 MHz separation between the bound state and the band edge with bound state linewidth of 4 MHz when a Bound State Frequency (GHz) Linewidth (MHz) FIG. 1.A platform for interacting dressed bound states a A 16-site microwave photonic crystal is realized by alternating sections of high and low impedance coplanar waveguide.Two transmon qubits (multi-level, anharmonic energy ladder) are in neighboring unit cells in the middle of the device, centered in the high impedance sections for maximal coupling to the band edge at 7.8 GHz (all values presented in units of (2π)Hz.i.e. -ω BE = 7.8 (2π)GHz).For this experiment, the pass band (band gap) refers to states above (below) the band-edge frequency.Each qubit is individually tunable in frequency via a local flux bias line.b Bound state linewidth, an indirect measure of localization, varies with qubit frequency.The wide range over which photon localization can be tuned indicates the feasibility of realizing a chain of strongly interacting bound states.
Experimentally measured and simulated linewidths are shown in red and black, respectively.c The interaction between bound states will be determined by overlap of their localized photonic envelopes with the qubits.d In a larger system, this localization-length-dependent interaction of the bound states would preserve the importance of the spatial organization of qubits in determining the many-body structure of the interactions.
qubit is resonant with the band edge, thus firmly reaching the strong coupling condition (see Fig. 1b and Supplemental Fig. S4a).By fabricating two qubits in the photonic crystal (Fig. 1a), we realize multiple, spectrally-resolvable bound states and can study inter-bound state interaction.
The nature of inter-bound state interaction makes this platform intrinsically well-suited for investigating one-dimensional chains of bound states (see Fig. 1c).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 near identical coupling to the band edge.As the strength of inter-bound 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 chain-like interaction pattern.
Controlling photon-mediated interaction between superconducting qubits has been demonstrated in other one-dimensional systems -such as two qubits in a cavity [20,21] or in a linear waveguide [22].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.2a and 2b, we characterize the dependence of the transition frequencies between the three lowest levels of the bound state on bare qubit frequency, and observe the steady reduction in bound state anharmonicity from over 350 MHz to 0 MHz as the qubit is tuned from deep in the bandgap to the passband.This is dramatically more than the ∼ 10% modification of qubit anharmonicity with frequency expected when a qubit is strongly coupled to a cavity mode [23].
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 qubit [11], it is important to note that this effective qubit differs in frequency and anharmonicity from the bare (multi-level transmon) qubit.Defining the three lowest bare qubit levels as |0 q , |1 q , and |2 q , here the two-excitation bound state is largely due to the coupling of the second qubit transition (|1 q ↔ |2 q ) with the band edge rather than other multi-photon effects [13] (see Supplement).
Numerical simulations, modeling the photonic crystal as a coupled cavity array with free parameters fit to match the band curvature from the dispersion relation [5,24], show similar dependence of anharmonicity on detuning (see Fig. 2a inset and Fig. 2b).Unlike the transfer matrix method [25][26][27], 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 qubit is modeled as a three-level ladder with negative anharmonicity, and with the |0 q ↔ |1 q and |1 q ↔ |2 q transitions coupled with amplitudes g and g √ 2, respectively, to its local cavity-site.It is critical to include level |2 q 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 (Fig. 2c), 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 [28], consistent with single photon transport of a two-level system and resonant pump (see supplement) [29][30][31].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 Transitions between all three eigenstates result in six sidebands [32].These six sidebands are visible in Fig. 2c, though emission intensity varies greatly among them due to eigenstate population.A seventh transition is evident in the data (7.25 GHz in Fig. 2c).This additional transition is due to the fourth effective qubit level (|3 ) while its curvature is reproduced by including a fifth effective qubit level (|4 ) in our numerical simulations (see Supplement).
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 [33].
We observe the flip-flop interaction between the two spatially separated bound states by measuring the avoided crossing in transmission when the bound states are tuned into resonance.As these qubits 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. 3a, 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.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.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 interactionfourth 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.b Numerical simulation for fixed bare qubit frequency of 7.27 GHz, with the one (two) excitation manifold in black (red).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. 3e and 3f, the qubits are maintained on resonance with one another while being simultaneously tuned through the bandgap.
Theoretical modeling (Fig. 3f) shows experimental data to be consistent with localized photonic states and with interaction via wavefunction overlap.In the limit where the qubit is deep in the gap, the Markovian approximation holds.Here the localized mode and flipflop interaction both have the same distance dependence e −x/L where L = a α δ is the localization length, δ 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 Supplement).The corresponding flip-flop interaction Hamiltonian is H ∝ j,l S + i S − j (−1) |x i −x j |/a e −|x i −x j |/L .To further study tunable on-site interaction, we probe the interacting bound states beyond the one-excitation manifold using spectroscopic measurements (see Fig. 4a).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 higher dressed 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 ) and |11 [21], remarkably we measure the two-photon virtual interaction between |20 and |02 , despite the fact that this process is 4th order in coupling g (see Supplement).This two-photon interaction shows consistent dependence on detuning: increasing in strength (from 0 to over 10 MHz) as the bound states shift towards the band edge and the states become more delocalized.Numerical simulations (Fig. 4b) 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 inter-bound state coupling possible via overlap alone.
The widely tunable on-site and inter-bound state interactions demonstrated with this device and consistent theoretical simulations are promising benchmarks towards realizing larger, more complex systems of bound states that can mediate both local and long-range interactions.Beyond stepped impedance coplanar waveguides, there are undoubtably 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 [11].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 [15].
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 [34][35][36][37].Another potential pathway would use an auxiliary microwave field through the crystal itself.Here the qubits could be maintained on resonance deep in the bandgap 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 re-dressing 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 [38]) 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 [15].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.
Supplemental Information: Interacting Qubit-Photon Bound States with Superconducting Circuits

I. MOTIVATION FOR USING PHOTONIC CRYSTALS
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 stopband and a passband that are "slow-light" due to reduced group velocity.Bandgaps 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 [S1].
Due to 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 bandgaps ("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 [S2].Engineering band edges via finite-size photonic crystals has been demonstrated in many systems, including nanophotonic structures [S3, S4] and superconducting microwave coplanar waveguides [S5, S6], and tremendous progress towards integration with ultracold atoms and superconducting qubits, respectively, has been reported [S6, S7].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.

A. 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 cos(
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 bandgap, 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.S1a and b show the simulated dispersion for unit cell parameters in Liu et al. [S6] 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 Transfer matrices are a convenient and accurate method for bare crystal simulation that incorporates both the exact number of cells and boundary conditions [S8].A convenient metric for comparison is transmission across the device, S 21 .In Figs.S1c and d

B. Realistic Parameters and Experimental Aside
The device is fabricated on a ∼ 430µ m thick c-plane sapphire substrate, on which we sputter approximately 200nm 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 photo-lithography, we are limited to impedances between approximately 25 Ω to 125 Ω.
For this device, one unit cell consisted of a high impedance section (Z hi = 124Ω, L hi = 7.8 mm) and a low impedance section (Z lo = 25Ω, L lo = 1.2 mm).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 x 10 mm sapphire chip.While more unit cells could fit on the chip, we saw experimentally that we must wire-bond extensively on-chip (to connect ground planes) to create clean bandgaps 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 GHz and 8 GHz.The frequency of the band edge was designed to be at 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 [S1] 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 resorted to two-tone spectroscopy over direct transmission measurement to detect the bound state when the qubit was deep in the bandgap 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. [S6]); however, this is a trade-off 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.

II. ADDING IN QUBITS
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 125kV e-beam writer, with bridge-free junctions, and are made of evaporated aluminum.The qubits were designed to have a target charging energy (from electro-magnetic simulation) of approximately 450 MHz, and emphasis was 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 uncertainty/error.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 (9mm 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 [S2].
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.

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

A. Hamiltonian
The Hamiltonian for the one-dimensional photonic crystal is given by 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.I A. Here, we have ignored the other bands of the photonic crystal as the qubits couple predominantly to the second band.Fourier transforming (a † k = a L N j=1 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, H c = i,j J i,j a † i a j , where Here, k is a dimensionless integration variable, and we have used the fact that J i,j = J j,i since ω k = ω −k .In Eq. (S2) we have taken N → ∞, as we only know ω k in that limit (see Sec.I A).To model our finite system, which has open boundary conditions and 16 unit cells, we use a 16 site hopping model with open boundary conditions with hopping strengths determined by Eq. (S2).Using the photonic crystal parameters in Sec.I B, we find (by numerical integration) J 0 = 9.3272 GHz, J 1 = 0.7288 GHz, J 2 = −0.0344GHz, J 3 = 0.0178 GHz, J 4 = −0.0034GHz, and J 5 = 0.0014 GHz.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 have calculated the hopping parameters for a given set of photonic crystal parameters.A different choice of photonic cyrstal 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 qubits is given by Here, i labels the qubit, n labels the level of the qubit and ω 0n;i are the bare energy levels of the qubits.In our simulation, the number of qubit levels is truncated at five (i.e.|0 through |4 ) .For our experiment, ω 00;i = 0, ω 02;i = 2ω 01;i − ∆, ω 03;i = 3ω 01;i − 3∆ and ω 04;i = 4ω 01;i − 6∆, where ∆ is the bare anharmonicity of the qubits, which is taken to be the same for both qubits.
We now turn to the coupling of the qubits to the photonic crystal.To an excellent approximation, the coupling between the qubit and the photonic crystal takes place within a single unit cell.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 qubits.In our system, the qubits are on neighboring unit cells, i.e. z 1 = N 2 and z 2 = N 2 + 1, and the coupling for each 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.

B. 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 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 halfwidth, κ 0 is a uniform contribution to photonic halfwidth, and κ is a decay rate on the first and last sites.While there are certainly other forms of loss (such as non-uniform loss on each site), our goal is to reproduce the key features (for example, 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 We have omitted 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 [S9].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.( S5)] [S9].More explicitly, the transmission for a given frequency ω is given by where Here, V l,r n = √ 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.

C. 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. S1d) is 1 GHz, although any κ in the range of .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 κ 0 in the range of 3 MHz to 5 MHz also fits the data well).Fig. S1d 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 [S10].Transmission is then calculated as a function of driving frequency.We find that g 2 = .55GHz and ω 01;2 = 7.9875 GHz match the experimental data well when the bound state is at 7.605 GHz, as seen in Fig. S2a.Calculating transmission when the first qubit frequency is near the pass band and the second qubit frequency is detuned and comparing it to experimental data, we find g 1 = .505[S10].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. S2b, 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 found that κ q = .5MHz accomplishes this task for ω 01;2 = 7.9875.To make sure we have chosen a suitable qubit halfwidth, we again tune the second bare qubit frequency (while keeping the detuned bare qubit frequency fixed).If we have chosen a reasonable qubit halfwidth, we should be able to accurately estimate the linewidth of the bound state for different bare qubit frequencies (while keeping all other parameters fixed).Fig. 1b of the main text shows that our estimate of κ q is reasonable.
Before moving on, we briefly comment on the second transmission method.Fig. S2c shows theoretical data from both simulations for g 2 = .55GHz 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 parameter, ∆.We first diagonalize H tot in the twoexcitation sector for fixed ω 01;i (with the other qubit detuned far away).The theoretical prediction for the dressed anharmonicity 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 ∆ until the theoretically predicted dressed anharmonicity 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 ∆ = .365GHz for both qubits.We then vary the bare qubit frequency and make sure the theoretically predicted dressed anharmonicity is still consis-tent 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.2a and 2b in the main text.
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.331 GHz, J 1 = 0.7308 GHz, J 2 = −0.0345GHz, J 3 = 0.0179 GHz, J 4 = −0.0035GHz, and J 5 = 0.0014 GHz.We see that the ratio of these hopping parameters to the previous set is not less than .97for any term, so we expect the other parameters in our model will not differ by more than 5 percent 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 = .555GHz and that the same range of decay parameters fit the data well, consistent with our expectation.

IV. BOUND STATE FUNDAMENTALS
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. S3).While the bound state is always within the gap, the frequency of the bound state changes with qubit frequency (see Fig. S4).
Additionally, the state becomes less localized as the bare qubit is tuned closer to the band edge (see Fig. S3a,b).
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 QED (see Fig. S3c,d).
In Fig. S4 we measure 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 wavefunction overlap as a change in bound state linewidth -where linewidth increases with localization length (see Fig. 1c in Ref. [S6]).As discussed in Sec.I, the localization of the photonic wavefunction is determined predominantly by the strength of the coupling, properties of the band edge, and the frequency detuning between the atom . Visualizing bound states -A qubit (pink circle) is coupled to one site of a 1D photonic crystal (gray line of alternating width).Coupling a qubit and band edge produces a photonic envelope (blue/purple) that is spatially centered at the qubit location.a and b The localization of the photonic component is determined by the detuning of the qubit from the band edge.a is more localized than b as the detuning is larger in the former.The overlap of the photonic wavefunction with the ends of the crystal (not shown) determines the linewidth of the bound state measured in transmission.The strength of the interaction between the bound states, when qubits are resonant with one another, can be understood in c and d as depending on the localization of the bound states.
and the band edge.The ability to tune the localization with small sized crystals shows the versatility of this platform.

A. Bound State Linewidth Dependence on Detuning
The linewidth of the bound state is set by the amplitude of the exponentially decaying photonic wavefunction 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.

B. 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 done in Ref. [S11] 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 (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 have chosen a quadratic dispersion, these results are qualitatively similar for a cosine dispersion.The most general single-excitation wavefunction is Here, the basis states are labeled as |qubit one 1 |qubit two 2 |photon .Solving the eigenvalue equation, H|ψ 1B = E 1B |ψ 1B , yields the following coupled equations, Solving for c k from Eq. ( S14) and inserting the result into Eq.( S12) yields 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 gives Repeating these steps for a 2 , gives the following equation for the bound state energy, where we used the fact that (z 1 − z 2 )/a 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, for example, Ref. [S11] or Ref. [S6]), Generically, Eq. (S20) yields one or two bound states depending on the coupling strength, qubit frequencies, distance between qubits and α [S11].
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. 3e 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. a e 1 = a e 2 or a o 1 = −a o 2 .The difference of the bound-state energy and the energy of the band edge, E 1B − ω 0 = δE 1B < 0, is then given by 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. [S11] is We explicitly consider the experimentally relevant case when 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.155GHz 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 = .5275GHz.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, one can investigate the real-space wave function.Fourier transforming c k gives 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. S3).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.

C. Single-Photon Bound States: Born-Markov Solution
We now briefly compare these exact results to the Born-Markov 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 highenergy subspace gives the following effective Hamiltonian for the two qubits, We stress that this formula is only valid for ω 01 < ω 0 .Diagonalizing H BM , we have the following eigenvalues, In the limit where the qubits are infinitely far apart, we recover the standard expression for the dressed qubit frequency, In Fig. S5, we compare the results obtained in the Born-Markov approximation to the exact analytical results.Fig. S5a shows these results for a single qubit and Fig. S5b 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. S5a to the experimentally measured bound state in Fig. S4, 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. S5b ends abruptly at ω 01 ≈ 7.920 GHz as there is only one bound state for ω 01 > ω 0 + .120GHz.
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. (S22) (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. (S25) around ω 0 ≈ ω 01 , we have 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. S5b.

D. Single-Photon Transport via the Bound State
As discussed, a bound state mediates transport across the crystal, at the otherwise forbidden frequencies inside the bandgap, 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. S6, 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. S6 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 [S12, S13, S14, S15, S16] for concept) was made possible by a TWPA (MIT Lincoln Labs) to improve SNR and a GPU (CUDA-Matlab) for significant computational speed-up.

V. EMISSION THEORY
In this section, 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 equa- Here, the blue squares are numerical results obtained from the total Hamiltonian Eq. (S29), the green squares are from the dressed-qubit Hamiltonian Eq. (S30), and the red squares are from Eq. (S30) but with cos θ = 1. a Here, ω 01 = 7.97 GHz, and the bound-state frequency is 7.591 GHz.Using the appropriately reduced matrix element (cos θ < 1, blue squares) is crucial in obtaining accurate results using Eq.(S30).b Here, the bare qubit frequency of 7 GHz and the bound-state frequency of 6.847 GHz are both far from the band edge, so that cos θ ≈ 1 and the reduction in the matrix element can be neglected, so all three data sets lie on top of each other.
When the bound state is approximately at 7.59 GHz, we find cos θ ≈ .68 for our finite system.Fig. S7 compares the results obtained from Eqs. (S29) and (S30).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. S7b) that, as expected, as the bare qubit frequency moves deeper into the gap, θ approaches zero and Ω approaches Ω.In Fig. S7, we have used the anharmonicity value predicted by theory.In Fig. 2c of the main text, we use the experimental values of anharmonicity, which is given by ω02 − 2ω d = −.11GHz.
We assume that the feature around 7.22 GHz in Fig. 2c of the main text is approximately In this section, we discuss the two-photon bound state.For a single qubit, the most general two-excitation wave function is Here the basis states are labeled as |qubit |photon .For i < j, it is convenient to define f i,j = f j,i .In Fig. S8, we plot |d i | 2 and |f i,j | 2 , which are obtained via exact diagonalization of the two-excitation sector for 16 sites [S18].We observe that the photons are localized around the qubit.In Fig. S8c, we plot the populations of qubit levels in the two-excitation ground state of H tot , which are given by to illustrate which terms in Eq. (S33) 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. S8c.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 qubit level in Eq. (S33) contributes significantly to the bound-state wavefunction.
Unfortunately, we are unaware of an exact solution similar to the one in Sec.IV B. 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. S8c These are similar to equations for the single photon case.Thus, we have We see that the photon is localized around the qubit, consistent with the exact finite-size numerical results seen in Fig. S8a.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 .

B. Two-Photon Avoided Crossing
In this section, we discuss the two-photon avoided crossing seen in Fig. 4 of the main text.As in previous sections, we assume the two qubits have equal coupling strengths and equal bare frequencies.To begin, we first write the total Hamiltonian in the two-excitation sector (in the rotating frame), using a different notation,

FIG. 2 .
FIG. 2. Probing the bound state energy levels aThe anharmonicity of the bound state is dependent on bare-qubit frequency, demonstrating a tunable on-site interaction strength.In blue (red), the first (second) transition of the bound state is measured across a range of bare qubit frequencies (inset -simulation).b Decreasing anharmonicity with increasing bound state frequency shown in red for experimental data and black for simulation.c Emission spectrum of a resonantly driven (∼ 7.59 GHz) bound state (induced by a qubit at 7.9 GHz, which is above the band edge located at 7.8 GHz) as a function of drive power.At low drive power, only the Mollow triplet is observed.With increasing power we see four additional sidebands, two on either side of the original Rabi sidebands, which together are the transitions between the three lowest levels of the anharmonic bound state (|0 , |1 , |2 ).The white crosses are from numerical simulations (see Supplement).We have included five qubit levels in our simulation.See text for the discussion of the seventh sideband around 7.25 GHz.

Frequency
FIG. S1.First two bands of simulated dispersion for Z lo = 25Ω, Z hi = 124Ω and a L lo = 0.45 mm, L hi = 8 mm (parameters from [S6]), and b L lo = 1.2 mm, L hi = 7.8 mm (parameters from present paper).c Overlay of simulated S 21 from the transfer matrix method (blue; same parameters as b) and measured high-power S 21 (black) shows good agreement in bare crystal characteristics.d Overlay of simulated S 21 from the transfer matrix method (blue; same parameters as b) and from the hopping model (red; with κ = 1 GHz and κ 0 = 4 MHz) and measured high-power S 21 (black) near the band edge.Drop in transmission circa 8.3 GHz is due to TWPA dispersion, not device defect.
, 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 [S6].

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 bandgap.The width of the bandgap is from 4.75 GHz to 7.8 GHz and needed 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 Fig S1aand S1b).

Frequency
FIG.S2.Dependence of transmission S 21 on drive frequency, dependence used for determining the coupling strength of the second qubit, i.e. the one at site N/2+1 = 9. a Solid blue line is theoretical data for ω 01;2 = 7.9875 GHz while blue are experimental data.We choose parameters of the Hamiltonian such that the bound state frequency, the transmission dip frequency, and the linewidth of the bound state match the experimental data.b Solid blue line is theoretical data for ω 01;2 = 7.941 GHz and blue dots are experimental data.Solid red line is theoretical data for ω 01;2 = 8.038 GHz and red dots are experimental data.c Comparison of transmission methods for ω 01;2 = 7.9875 GHz.The solid blue line is from method one, while the solid red line is from method two.The bound state and transmission dip occur at the same frequencies for both methods.
FIG. S4. a Experimental data and b hopping model simulation for S 21 vs. single-qubit frequency and probe frequency.The bare band edge is at 7.797GHz.The bright peak in the bandgap is the dressed qubit-photon bound state.The bound state always exists within the bandgap for qubit frequencies both above and below the band edge -a clear signature of non-Markovianity.In this figure, the other qubit is far detuned and has negligible effect.

FIG. S5 .
FIG. S5.Exact solution versus Born-Markov Approximation.Here, the red lines are the exact solution while the blue lines are the Born-Markov approximation.The dashed black lines mark the band edge and the thick black line is the bare qubit frequency.Here we take, ω 0 = 7.8 GHz, α = 1.155GHz and g = .5275GHz. a Single qubit.b Two qubits for |z 1 − z 2 | = 1.

FrequencyFIG. S7 .
FIG.S7.Theoretical simulations of the emission spectrum for different pump powers.Here, the blue squares are numerical results obtained from the total Hamiltonian Eq. (S29), the green squares are from the dressed-qubit Hamiltonian Eq. (S30), and the red squares are from Eq. (S30) but with cos θ = 1. a Here, ω 01 = 7.97 GHz, and the bound-state frequency is 7.591 GHz.Using the appropriately reduced matrix element (cos θ < 1, blue squares) is crucial in obtaining accurate results using Eq.(S30).b Here, the bare qubit frequency of 7 GHz and the bound-state frequency of 6.847 GHz are both far from the band edge, so that cos θ ≈ 1 and the reduction in the matrix element can be neglected, so all three data sets lie on top of each other.

ω 23 (
which gives ω03 − 3ω d = −.48GHz).This assumption is in decent agreement (around 50 MHz off) with results obtained by diagonalizing the full system when the bound state is at 7.59 GHz.This 50 MHz disagreement can be traced back to the 20 MHz disagreement in ∆, the anharmonicity, between theory and experiment when the bound state is at 7.59 GHz (if ∆ is off by 20 MHz, level |3 is expected to be off from its value by around 3 times this amount as ω03 ≈ 3ω 01 − 3∆).Finally, we find that ω04 − 4ω d = −1.78GHz (this value is also consistent with results obtained by diagonalizing Eq. (S29)) matches the experimental data well as seen by overlaying the theoretical data from the dressed qubit with the experiential data (Fig.2cof main text).In particular, we have captured the feature that appears around 7.22 GHz and −10 dB.The unique bending of this feature can be traced back to the fact that the level structure of the dressed qubit [Eq.(S30)] does not behave like a normal transmon due to the strong coupling to the photonic crystal.VI.MULTIPHOTON THEORYA.Two-Photon Bound State ). Solving the eigenvalue equation, H| ψ2B = E 2B | ψ2B , for the following wave-function ansatz (in momentum-space),| ψ2B = b|2 |0 + k dk a † k |1 |0 ,(S35)yields the following equations,bE 2B = ω 02 b + √ 2g √ n k e ikaz 2 dk ,(S36)dk E 2B = (ω k + ω 01 ) dk +

e
ika(j−z 2 ) E 2B − ω 01 − ω k = − √ 2 b cos(π(j − z 2 )) 2 α(ω 0 + ω 01 − E 2B ) e − ω 0 +ω 01 −E 2B α |j−z 2 | .(S39) FIG. S8.Visualizing the two photon bound state.Here, the qubit is on site nine.a The blue dots are |d i | 2 versus position for ω 01 = 7.6 GHz.The red dots are |d i | 2 for ω 01 = 8.15 GHz as a function of position.b Plot of |f i,j | 2 as a function of i and j for ω 01 = 8.15 GHz.c The population of qubit levels in the two-excitation ground state of H tot as a function of bare qubit frequency.

e Measuring bound state separation as qubits are simultaneously tuned, maintaining resonance, shows changing bound state interaction strength with frequency. f Bound state avoided crossing 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 bandgap) due to increasing localization of the bound states. Hopping model simulation (black) captures
Interacting Bound States -Interaction between bound states is characterized by the avoided crossing seen in transmission while tuning one qubit (y-axis) through resonance with the other (fixed).a An avoided crossing of 240 MHz is observed when the fixed qubit is at 7.73 GHz.The two points where transmission amplitude of a bound state dims are understood as the bound state peak being resonant with the qubit frequency.a (inset) Hopping model simulation of the one-excitation manifold is consistent with experimental observation.The lamb shift in the hopping model originates from next-nearest neighbor interaction between coupled cavities.b, c, d Tunable bound state interaction strength is illustrated in example 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.
this detuning-dependent behavior.