Extensible quantum simulation architecture based on atom-photon bound states in an array of high-impedance resonators

Engineering the electromagnetic environment of a quantum emitter gives rise to a plethora of exotic light-matter interactions. In particular, photonic lattices can seed long-lived atom-photon bound states inside photonic band gaps. Here we report on the concept and implementation of a novel microwave architecture consisting of an array of compact, high-impedance superconducting resonators forming a 1 GHz-wide pass band, in which we have embedded two frequency-tuneable artificial atoms. We study the atom-field interaction and access previously unexplored coupling regimes, in both the single- and double-excitation subspace. In addition, we demonstrate coherent interactions between two atom-photon bound states, in both resonant and dispersive regimes, that are suitable for the implementation of SWAP and CZ two-qubit gates. The presented architecture holds promise for quantum simulation with tuneable-range interactions and photon transport experiments in nonlinear regime.


I. INTRODUCTION
Quantum emitters coupled to structured photonic environments constitute both an emerging paradigm of quantum optics [1,2] and a promising platform for quantum information processing [3][4][5][6] and quantum simulation of many-body physics [7][8][9][10].One-dimensional photonic lattices modify the electromagnetic environment, leading to the appearance of finite bands and band gaps in the energy spectrum.A major phenomenon in these systems is the formation of atom-photon bound states in the photonic band gaps [11][12][13][14][15][16][17][18][19][20][21].These states originate from the dressing of the atom with a photonic cloud that remains exponentially localized in its vicinity, thus inhibiting a full atomic decay.In addition, when multiple atoms are coupled to the same photonic lattice to form atom-photon bound states, their interaction is mediated by the overlap of their photonic wavefunctions.Since the photonic localization length can be controlled by varying either the frequency of the atom or the strength of its coupling to the lattice, this architecture supports tuneablerange interactions, opening opportunities for quantum simulation of exotic spin models [22][23][24][25] and quantum computing architectures with connectivity beyond nearest neighbor [26].
Atom-photon bound states have been observed in different systems, ranging from cold atoms coupled to photonic crystal waveguides [27], to optical lattices [28,29], to superconducting circuits [30,31].A seminal experiment in superconducting circuits [30] relied on a microwave photonic crystal consisting of a coplanar waveguide with periodically modulated impedance.In a photonic crystal, the lattice periodicity must be of the same order as the wavelength of the radiation in the band gap.At microwave frequencies, this constraint results in a large footprint, which limits the number of unit cells that can be accommodated on a chip to a dozen and hinders the integration of additional measurement and control circuitry.A new avenue was opened by the introduction of superconducting metamaterials with a deep subwavelength lattice constant, consisting of a set of lumpedelement resonators periodically loading a transmission line, or arranged in a linear chain to form a coupled-cavity array [32][33][34].In these embodiments, the lattice footprint is drastically reduced (nearly two orders of magnitude) and the photonic cloud is more strongly confined.In addition, by staggering the hopping amplitudes between neighboring resonators, lattices with non-trivial topology have been realized, which host topologically protected bound states [34].
Here we introduce a novel circuit quantum electrodynamics (QED) implementation of photonic lattices coupled to quantum emitters, which employs arrays of highimpedance resonators and transmon qubits.Compared to previous work, our implementation features a reduced footprint, an enhanced interaction strength, and a higher extensibility.These benefits are obtained by utilizing arrays of Josephson junctions as compact inductors in the resonator array.We present spectroscopy and timedomain measurements of a proof-of-principle device comprising an array of 21 high impedance resonators and 2 transmon qubits with dedicated measurement and control circuitry.We characterize the mode structure of the array, the emergence of atom-photon bound states, their anharmonicity, the exchange and cross-Kerr interaction between two atom-photon bound states.Our results demonstrate that the presented architecture is endowed with all essential building blocks to carry out quantum simulation and quantum information processing tasks, and to access nonlinear regimes of quantum optics.The foreseen possiblity to accommodate multiple emitters with a limited increase in physical footprint looks particularly promising for carrying out quantum simulations of large many-body spin Hamiltonians [22].

A. Sample and experimental setup
We implement the structured photonic environment as a transmission line made out of 21 high-impedance microwave resonators forming a coupled-cavity array (Fig. 1).Each resonator consists of an array of 10 Josephson junctions of total inductance L r =8.87 nH shunted by a capacitor C r = 91.3fF, resulting in a bare resonant frequency ω r = 5.593 GHz and a characteristic impedance Z r = 312 Ω. Nearest-neighbor resonators are capacitively coupled to form a linear chain and each edge resonator is coupled to a 50 Ω coplanar waveguide.The artificial atoms (Q 1 and Q 2 ), are implemented as superconducting flux tunable transmons [35] capacitively coupled to sites 10 and 12 of the array.Each transmon is additionally coupled to a charge line (XY control), a flux line (Z control), and a readout resonator.The device is realized in aluminum on a silicon substrate with a standard lithographic process [36].The sample is wire bonded to a copper sample holder thermally anchored to the mixing chamber stage of a dilution refrigerator at 10 mK.A summary of the relevant sample parameters and further experimental details are presented in App. A.

Bare coupled cavity array
To characterize the bare coupled cavity array, we tune the resonant frequencies of Q 1 and Q 2 far away from the transmission band and measure the transmission coefficient through the array (Fig. 2).At low driving powers (average photon number in the array mode n ≈ 1), a structure of N = 21 modes forming a pass band emerges in the transmission spectrum, from the first mode at ω 1 /2π = 5.219 GHz to the last one at ω 21 /2π = 6.215GHz.These modes present rather uniform spacing and linewidths in the center of the band, while they concentrate and become narrower in linewidth at the band edges.These features are captured by a tight binding model for the array with only three parameters: the single-cavity bare resonance frequency, ω r , the cavitycavity nearest-neighbor coupling, J, and the coupling of the edge cavities to the input and output transmission lines, κ (see App. B).In the limit κ J, which applies here, the resonance frequencies of the modes are given by the linewidth of each being Eq. ( 1) predicts a pass band of width 4J centered around the bare cavity frequency ω r .From Eq. ( 1), 4J ≈ ω 21 −ω 1 and we extract J/2π = 249 MHz.Comparing it to previous realizations of microwave coupled cavity arrays [33,34], we achieve a larger value for J with a smaller coupling capacitance, thanks to the higher impedance of the resonator (J = 1 2 C J ω 2 r Z r ).From electrostatic simulations we estimate the resonator impedance to be Z r = 1/ω r C r ≈ 312 Ω.This implies a sixfold gain on the capacitive coupling strength compared with 50 Ω resonators.
Compared to the predictions of this simple model, the measured traces present deviations in the frequency distribution of the modes, that we attribute to a 3% scatter in the resistance of the fabricated Josephson junctions [37], resulting in an estimated standard deviation of δω r /2π =25 MHz in the frequencies of the bare (uncoupled) resonators in the array.Importantly, thanks to the small ratio δω r /J = 1/10, this frequency disorder does not significantly affect the properties of the atom-photon bound states, as calculated by numerical diagonalization of the system Hamiltonian for various realizations of the disorder, and directly verified in our experiments below.
In addition, we measure a nonzero transmission outside the photonic pass band, and an alternating transmission background in between resonant modes, which we ascribe to direct cross-talk between the input and output ports of our sample box.For a detailed treatment of these experimental imperfections, see App. C. At higher powers, corresponding to n ∼ 10 2 , each individual mode exhibits the typical phenomenology of a Kerr resonator, due to the nonlinearity inherited by the arrays of Josephson junctions [38].In our design, we estimate the Kerr coefficient for a single mode to be K/2π = 100 kHz, much smaller than the mode linewidth.In fact, the frequency shift produced by more than 100 photons is still smaller than a linewidth as visible in Fig. 2(b).Importantly, this value can be determined by design, both choosing the number of junctions in a single resonator, or the number of unit cells in the resonator array.The extraction of this parameter, and the relation between the Kerr coefficient of individual resonators and that of the array modes, are discussed in App.A 3. In the remainder of this work, we will only consider the linear regime of the coupled cavity array.

Dressed coupled cavity array and atom-photon bound states
We characterize the interaction of each qubit with the coupled cavity array by tuning its frequency across the pass band, while keeping the other qubit detuned [Fig.3(a)].Here we discuss the results for Q 2 ; the results for Q 1 are comparable (see App.A 2).
The low-power transmission coefficient of the array is affected by the presence of the qubit [Fig.3(b)].We observe a minimum in transmission within the band in correspondence with the bare qubit frequency (red dashed line in Fig. 3(b)).In fact, in the low excitation regime all the incoming field is coherently scattered back [39].This effect is particularly visible when the bare qubit frequency is resonant with one of the coupled-cavity array modes.In Sec.II C we show an atom-cavity interaction strength is approximately g i /2π ∼ 300 MHz, five times larger than the frequency spacing of the modes, implying a multimode interaction.In particular, the absence of an avoided crossing between a single mode and the artificial atom, and the monotonic dispersive shift of a mode, indicate the interaction of the artificial atom with multiple modes.Moreover each mode presents a definite standing-wave spatial profile, which sets its effective interaction strength with the artificial atom.For example, mode m = 10, indicated by the white arrow in Fig. 3(b) and (c), is completely decoupled from the qubit due to a corresponding node on the site of the artificial atom.
Comparing our system with a previous realization of multimode strong coupling [40], the vanishing group velocity, v g , at the band edges produces a nontrivial density of photonic states, proportional to 1/v g [41].As the atom frequency approaches the band, it hybridizes with band-edge photons having zero velocity [17,30].This seeds an evanescent field, exponentially localized around the atom (shaded area around the atoms in Fig. 1(a)).In Fig. 3(d) we observe this additional photon-like mode outside the passband, which asymptotically approaches the band edge when the bare qubit frequency is moved towards the center of the band.
To model the transmission spectra, we introduce the Hamiltonian where we introduce the photon annihilation operator a x for the x-th cavity, the transition frequencies ω qi of qubits Q i , their anharmonicities β i , and finally their couplings with the x i -th cavity, g i .We calculate the transmission coefficient in the limit of linear response from inputoutput theory and find a qualitative agreement with our measurements [Fig.3(c); see App.B 6 for details].
When the frequency of the qubit is tuned towards the low frequency edge of the band, the bound state at higher frequency (upper bound state) completely loses its atomic nature, becoming the mode at the high frequency edge of the band.Conversely, the mode at the low frequency edge starts to get dressed with an atomic component (lower bound state, barely visible in Fig. 3(b)), as predicted by the Hamiltonian (3) for a single qubit.Their frequencies ω BSi (i = 1, 2 for Q 1 and Q 2 , respectively), are given by the solutions of the equation with ω BSi −ω r > 2J (see App. B 3 for more details).The upper bound state will be our main focus for the rest of the paper, and we will use its frequency as the independent variable to describe our measurements.As pictorially shown in Fig. 1(b), this state is highly localized and atom-like in nature for atomic bare frequency deep in the gap while it is more photon-like and delocalized for atomic frequency close to the pass-band.These features are directly exploited for exciting the bound state: the finite overlap of the photonic cloud with the resonators at the edges of the array allows for the detection of this state in transmission (see Fig. 3(b)-(d)).
The decay rate of the bound state, extracted from the linewidth of the transmission spectroscopy as function of the bound state frequency, is of the order of γ/2π ≈ 300 kHz, close to the decay rate of the array modes [Fig.3(e), orange dots], due to the large and delocalized photonic component of the bound state [Fig.3(e1)].When the bound state is far from the band, we extract its decay rate from a measurement of its atomic component, via the readout resonators.In this case, the losses are much smaller [Fig.3(e)] due to the mainly atomic component [Fig.3(e2)], and we measure a decay rate γ/2π ≈ 50 kHz.4).(e) Total decay rate of the bound state as function of its frequency.The red data are obtained from the transmission spectroscopy while the blue from the atom spectroscopy (see Sec.II C).The dashed black line represents the theoretically expected decay rate.Insets: calculated distribution of the excitation over the array (green dots) and between the two qubits (blue and red respectively) for two distinct bound state frequencies: ωBS2 = 6.226GHz (e1) and ωBS2 = 6.510GHz (e2).

C. State preparation of a single atom-photon bound state
1. Single-excitation subspace In the limit in which the atom-photon bound state is localized, and thereby not accessible by transmission spectroscopy of the array, we excite it by sending microwave pulses to the qubit via its charge line, and detect it by performing dispersive qubit readout using the readout resonators [see Fig. 1(c1) and pulse scheme in Fig. 4a].A 90 ns Gaussian pulse with frequency ω p and calibrated amplitude to be π-pulse when Q 2 is at its maximum frequency, ω q2,m , is sent to the XY-Q 2 control line (red).At the same time, the frequency of Q 2 , ω q2 , is set with a 140 ns square pulse (with 2 ns rise and 2 ns fall time) on the magnetic flux control, Z-Q 2 with varying amplitude Φ q2 (grey line and red shaded area).The frequency of Q 1 is dynamically tuned to its lowest value with a flux pulse of half flux quantum Φ 0 /2.After 50 ns a 2 µs square readout pulse sent to the multiplexed resonators, is used to read out the average qubit population, P 2 .The excitation spectrum of the bound state is obtained by measuring P 2 as a function of pulse frequency and magnetic flux applied to Q 2 [Fig.4(b)].When the pulse on XY-Q 2 is resonant with the BS transition frequency, the qubit is efficiently excited.The bound state frequency strongly differs from the bare qubit frequency [red dashed line in Fig. 4(b)], and asymptotically approaches the band edge for qubit frequencies in the band [Fig.4(c)].The bound-  state frequency is also in this case well described by the continuum theory, as shown by the solid black line in Fig. 4(b), representing the best fit of Eq. ( 4) to the data.We notice that, keeping the pulse amplitude constant, the bound state population decreases as its frequency approaches the band edge.Compared to the drive rate Ω R,0 for an isolated qubit excited via its charge line by a pulse of given amplitude, the bound state is subject to a reduced drive rate Ω R = Ω R,0 cos(θ), where θ is the mixing angle between the atomic and photonic components.This relation allows us to extract the mixing angle from the measurements in Fig. 4(b), which we find to be in good agreement with our theoretical prediction [Fig.4(d)].

Double-excitation subspace
The physics of the bound states explored so far was restricted to a single excitation.When higher-excitation subspaces are considered, the nonlinear nature of the transmon starts to play a role and leads to deviations from the linear regime.To explore these nonlinear features, we "climb up" to the two-excitation subspace using a sequence of two pulses [Fig.5(a)].We first send a resonant π-pulse to excite the long-lived, single-excitation bound state.We subsequently send another pulse of fixed amplitude and varying frequency ω p2 to search for the first-to-second-excited state transition of the bound state.Finally, because our readout pulses are optimized so that the ground state is the most distinguishable one, we find it convenient to add a third pulse identical to the first, so that if the second pulse does not affect the bound state, then the system is brought back to the ground state.Such a scheme allows us not only to perform a complete spectroscopy of the second excitation subspace but also, compared to previous approaches [31], to efficiently prepare the two excitation bound state.
In order to quantify the deviation from the linear case we define the dressed anharmonicity parameter as BS2 where the superscripts (1) and ( 2) stand for the number of excitations.Experimentally, this quantity is determined as the difference in frequency between the second and first pulses when the second transition is resonantly excited.The measured β dress is always negative; its magnitude is maximum when the bound state frequency is the farthest from the band edge, and monotonically decreases towards zero as the bound state frequency approaches the band edge [Fig.5(b), dots].To theoretically capture the measured nonlinearity, we describe the transmon as a nonlinear resonator with bare anharmonicity β/2π = −257 MHz (see Eq. ( 3)).This quantity is estimated by performing the same measurement as in Fig. 5(a) for ω q2 /2π = 3.329 GHz, well below the band and the read-out resonator frequency.We use the parameters J, ω r and g 1 extracted from the fit of the single excitation bound state spectroscopy, to diagonalize the Hamiltonian in the single-and double-excitation subspace (see App.The observed anharmonicity for the bound state is intermediate between that of a fully linear emitter, for which no anharmonicity would be observed, and that  of the most nonlinear emitter, a genuine two-level atom (β → −∞), plotted for comparison as a dashed line in Fig. 5(b).A detailed numerical study, presented in App.B 5, indicates that double-photonic, double-atomic, and hybrid excitations are all present across the frequency range considered.In addition, we find that the localization length of the doubly excited bound state is also renormalized according to the nonlinearity of the emitter, as originally discussed in [17].

A. Two atom-bound states level splitting
To investigate the interaction between bound states, we bring their frequencies close to each other, send an excitation pulse to one of the qubits, and perform joint readout of the two qubits [Fig.6(a)].Sweeping the frequency of one qubit while keeping the other one fixed, we observe characteristic avoided crossings in the excitation spectrum, detected as peaks in the measured populations of both qubits as a function of the qubit and probe frequencies [Fig.6(b-g)].Interestingly, the excitation pulse on Q 2 , excites a population fraction in both bound states.In fact, when the two atom-photon bound states resonantly interact, they hybridize forming even (+) and odd (−) states (see App.B 4).In Fig. 6(b-g) we report the measurement results for three different frequencies of Q 1 : when it is tuned well inside the band ω q1 /2π ≈ 5.9 GHz, (b) and (e), at the band edge, ω q1 /2π ≈ 6.2 GHz, (c) and (f), and finally at its largest bare frequency possible, ω q1 /2π = 6.332GHz (d) and (g).
Tuning the bound state frequencies closer to the band edge, results in a larger avoided crossing, corresponding to a larger interaction strength.It is important to notice that the two qubits do not present a direct coupling, in fact their interaction derives from the mutual overlap of their photonic clouds which becomes progressively more extended and populated closer to the band edge [see Fig. 3(e1)-(e2) where we calculated the photonic cloud distribution of the bound state depending on the distance from the band edge].
When the frequencies of the interacting bound states approach the band edge, we observe a vanishing population of the (-) bound state, when its frequency crosses the band edge [Fig.6(b,e)].This phenomenon, referred to as "melting" of one bound state into the modes of the band [17], only occurs in the presence of two interacting bound states, (+) and (-).The even state is characterized by a bonding behavior and its energy is raised compared to the one of the individual bound state, leading to an increased localization of its photonic cloud.The odd state instead presents an anti-bonding behavior, with its energy pushed towards the band edge making the state more extended.This progressive delocalization leads to a disappearance of this bound state solution.
We define the bound states interaction strength U = [ω BS+ (ω q ) − ω BS− (ω q )]/2 as half of the level splitting when the bare qubit frequencies are on resonance, ω q = ω q1 = ω q2 , where the ± stands for even and odd dressed bound states, respectively.This quantity is displayed in Fig. 6(h) (red dots) as a function of the frequency at the midpoint of the level splitting, ωBS = ω BS− + U/2, which gives the approximate frequency of the individual (single artificial atom) bound state.We measure the bound states interaction strength ranging from 27 to 52 MHz in a frequency range of 250 MHz.This tunability is directly related to the variation of the overlap between the photonic clouds.Based on numerical simulation, we predict that the interaction strength could be tuned by two orders of magnitude through a straightforward optimization of the device parameters (see App. C 2).Fig. 6(h) also shows that the coupling does not present a monotonic behavior.In fact, it increases for bound state approaching the band edge, but it exhibits a maximum close to the melting condition.This behavior relies on our definition of the BSs coupling and on the melting condition.Indeed, when the odd bound state merges with the continuum in the band, the splitting is measured between the even bound state and the band edge.
Pulse frequency, Next-nearest neighbor Nearest neighbor Separately, we find that the expected interaction strength calculated with the eigenvalue equation derived by the system Hamiltonian Eq. (3) (see App. B 4 for more details) systematically underestimates the measured coupling strength [Fig.6(h), dotted line].As shown in App.C 1, parasitic capacitances are responsible for a nonnegligible coupling between each qubit and their nextnearest array sites.In particular, the bound-state interaction is affected by a direct coupling between the qubits and the cavity in-between them (x = 11).After including these additional couplings in the Hamiltonian, we obtain a much better agreement with the data [Fig.6(h), solid  line].This behavior cannot be explained, in contrast, by adding just a direct qubit-qubit coupling to the model.

B. Time resolved excitation exchange
After establishing a static interaction between the bound states, we now exploit it in a dynamic setting to realize an excitation swap between two bound states [Fig.7(a)].When both qubits are at their largest frequencies, we excite the bound state of Q 2 with a π-pulse and we tune it in resonance with the bound state of Q 1 , by applying a flux pulse on each individual flux line.After an interaction time τ , we bring the two qubits back to their initial frequencies and we read-out the populations P 1 and P 2 , in Q 1 and Q 2 , respectively.Fig. 7(b) and (c) show P 1 and P 2 as a function of the bound states interaction time, τ , and Q 2 frequency ω q2 , while the Q 1 bound state frequency is kept at ω q1 = ω q1,m in (b)-(c), and ω q1 /2π = 6.15GHz in (d)-(g).These are the same parameters chosen for measuring the avoided crossing in Fig. 6(d) and (g).The population of the two qubits as a function of the interaction time, restricted to resonant bound states, is presented in Fig. 7(f) and (g).These data, corresponding to a slice of the chevron pattern (dashed blue and red line in Fig. 7 (b)-(e)), can be fitted to a dumped sinusoidal function (solid red and blue line), from which we extract a complete swap time of 18 ns.This value is in good agreement with the measured interaction strength.The excitation swap between the two bound states is realized by non adiabatically tuning the two bare bound states in resonance.The adiabatic threshold for this process is related to the interaction strength [42], with a relative time scale ∼ 1/U ∼ 5 ns.The flux pulse we implement has a rise (and fall) time of t raise = 1 ns, thus we fulfill this constraint.Nevertheless the two bound states interact during the rise time, and we take this into account with a time shift of the measured data in Fig. 7(f).
The Chevron pattern of the population P 1 and P 2 measured for the bare frequency ω q1 /2π ≈ 6.02 GHz shown in Fig. 7(d)-(e) highlights a faster excitation swap of 13 ns, as expected from the larger interaction strength closer to the band edge.Focusing only on the resonant case (see Fig. 7(g)) we notice that after the pulse protocol we observe a total population of only P 1 + P 2 ∼ 0.7.The population loss originates from the weakly adiabatic regime in which we operate, i.e. t raise 1/∆ ωBS ∼ 0.5 ns, where ∆ ωBS is the gap between the average individual bound state frequency ωBS and the band edge (see Sec. III A).
When the bare qubit frequency is tuned in the band faster than the latter adiabaticity threshold, the bare qubit state is not an eigenstate of the system.In this case, part of the atomic population redistributes among the resonators according to the projection of the bare qubit state on the bound state.In this process, a fraction of the atomic population is converted into itinerant photons and released into the waveguide.In particular, the released population is at most P released sin 2 (θ) with θ being the mixing angle θ defined in Sec.II C. In fact, where the bound states are mainly qubit-like, the total population P 1 + P 2 ∼ 0.95 (see Fig. 7(f)), while a larger photonic dressing preserves a smaller fraction of population, P 1 + P 2 ∼ 0.7, as shown in Fig. 7(g).This last value could be improved by slowing down the protocol but still keeping t raise 1/U in order to induce the excitation swap in the first place.
The lifetime of the interaction is determined by the lifetime of the average individual bound state (dashed black lines in Fig. 7(f)-(g)), which we extract from the fit to be 1.2 µs and 950 ns, respectively.

C. Two atom bound states ZZ interaction
When the two atom photon bound states are close in frequency to the band edge, but detuned from each other, exciting one of them results in a shift in the transition frequency of the other.This ZZ (or cross-Kerr) interaction originates from dispersive and resonant interactions between energy levels in the double-excitation manifold.The interaction mechanism relies on the overlap of a photonic cloud with finite localization, as it happens in the single excitation subspace.
We investigate the ZZ interaction with the pulse scheme displayed in Fig. 8(a).After setting the bare qubit frequencies with a flux pulse on each flux line, we excite Q 1 bound state with a π-pulse.We apply a second pulse to Q 2 with frequency ω p2 , and finally read out the state of the qubits.When the pulse on Q 2 is such that BS2 , the bound state of Q 2 is excited.We repeat the same pulse sequence by tuning the bare frequency of Q 1 , mapping the transition frequencies ω BS2 as the function of ω q1 .
Fig. 8(b) shows the ZZ interaction between two excited bound states for the bound state of Q 2 tuned at the constant frequency ω BS/2 /2π = 6.653GHz, while the frequency of the Q 1 bound state is tuned from band edge toward its largest frequency 6.55 GHz.Following the level |φ 11 on the right panel of Fig. 8(c), we notice that this level crosses |φ 02 relative to the two-excitation bound state of Q 2 .In Fig. 8(b), at lower frequencies, we measure a residual ZZ interaction up to 49 MHz.The dashed line shows the calculated ZZ interaction numerical evaluated with the complete Hamiltonian Eq. (C3).
Note that in [31] the avoided crossing between |φ 02 and |φ 11 was infer using high power coherent spectroscopy.With this approach, second order transitions may induce Stark shift in the measured frequencies, while the measurement of the ZZ interaction with a pulse scheme is directly relevant for gate implementation.
Fig. 8(c) shows the energy structure of bare qubits (left panel) and the same energy levels for bound states as function of the bare frequency of Q 1 .The color scale highlights the population of the two atoms going from blue, for Q 1 , to red for Q 2 , where the yellow color instead stands for excitation equally distributed on the two.Fig. 8(c) is a simplified version of the complete energy structure of our system.In fact, the three photonic bands that originate from the coupled cavity array, are not included in the figure (see App. B 5 for more details).

IV. DISCUSSION AND CONCLUSION
In summary, we introduced a new implementation of a finite band waveguide made out of an array of high impedance resonators to which we coupled two artificial atoms.We measured the transmission spectrum of the coupled resonators array, highlighting the underlying multimode interactions with the artificial atom, and observed the formation of the atom-photon bound states in the band gap.We demonstrated full control in accessing and preparing the atom-photon bound states in both the single and double-excitation subspace.We characterized the resonant interaction between two bound states in the static and dynamic regime, measuring an effective coupling strength up to 52 MHz, and an excitation swap time down to 13 ns.Finally, we investigated the cross-Kerr (ZZ) interaction between two detuned bound states, reaching a value up to 49 MHz.As shown in simulations, we expect a straightforward optimization of the device parameters to significantly improve our ability to control as well as suppress these interactions.
Comparing to previous implementations of superconducting qubits coupled to gapped waveguides, our approach based on high impedance resonators makes it possible to reach atom-cavity coupling strengths of a few hundred MHz, while maintaining the resonators footprint comparable to the one of the artificial atoms.Small footprint and strong interactions translate in a higher extensibility of our platform, with the foreseen possibility of adding more qubits as well as anchoring points to move towards two-dimensional lattices [43].
The tuneable-range interactions between atom-photon bound states available in this platform find application well to the quantum simulation of spin models [22].At the same time, the possibility to implement fast and high-contrast SWAP and CZ gates using the array as a quantum bus could be further investigated in the context of quantum computing.From the perspective of quantum optics, this platform is amenable to studies of correlated nonlinear photon transport [14,44,45] and quantum nonlinear optics protocols [15].By varying the coupling strength between neighboring sites, it is possible to engineer the band structure of the array and to endow it with nontrivial topological properties [34].The intrinsic nonlinearity of the array, whose strength can be adjusted by design, could be utilized to implement recent theory proposals describing exotic light-matter interaction effects [46][47][48].Finally, from the viewpoint of materials, replacing Josephson junction arrays with highkinetic inductance superinductors [49,50] may lead to a more robust fabrication process and reduced disorder in the array.

ACKNOWLEDGMENT
The authors are grateful to Peter Rabl, Angelo Carollo, and Christopher Warren for useful discussions.MS, PD and SG wish to express their gratitude to Lars Jönsson for making the sample holder.They acknowledge financial support from the Swedish Research council and the Knut and Alice Wallenberg Foundation.DEC acknowledges support from the European Union's Horizon 2020 research and innovation programme, under European Research Council grant agreement No 101002107 (NEWSPIN); the Government of Spain (Europa Excelencia program EUR2020-112155, Severo Ochoa program CEX2019-000910-S, and MICINN Plan Nacional Grant PGC2018-096844-B-I00); Generalitat de Catalunya through the CERCA program, Fundació Privada Cellex, and Fundació Mir-Puig.The device was fabricated at Myfab Chalmers.

Appendix A: Experimental details 1. Experimental Setup
The complete experimental setup is shown in Fig. 9.The sample is wire-bonded in a nonmagnetic oxygen-free copper sample box (see App. C 3), mounted to the mixing chamber stage of a dilution refrigerator operating at 10 mK, and shielded by an additional copper and µmetal cans.Two additional shields (in order, copper and µ-metal) are placed around the sample in a light-tight fashion.The signal to the Q 1 and Q 2 charge lines (XY-Q 1 and XY-Q 2 respectively), to the waveguide (WG in ) and to read-out resonators (RO in ) is delivered by an highly attenuated coaxial lines (nominal total attenuation of -60 dB).The flux lines (Z-Q 1 and Z-Q 2 ) do not have any attenuation at the last two stages (total nominal attenuation -33 dB) but they are equipped with a 4 GHz low pass filter for reducing flux noise from higher stages.The transmitted signals (through the couple cavity array and through the readout feedline) is amplified with high electron mobility transistor (HEMT) amplifier at 3 K, further amplified at room temperature.The microwave control for Q 1 (blue area) and Q 2 (red area), is achieved up-converting the inphase (I) and the quadrature (Q) components of a low frequency pulse generated by an arbitrary waveform generator (AWG), while the flux control is obtained injecting current generated directly by two AWG channels.The qubits readout pulse is up and down-converted with the same local oscillator (LO) (yellow area).Finally, the waveguide coherent spectroscopy is set up with a vector network analyzer (VNA) (green area).

State preparation of Q 1
Using the same pulse scheme as the one shown in Fig. 10(a) applied on Q 1 we measure its population as a function of the bare qubit frequency (dashed white line) and driving pulse frequency, ω p .The result, reported in Fig. 10(a), shows Q 1 largest frequency, ω q1,m = 6.332GHz, smaller than Q 2 , but compatible with the fabrication yields.The fit of the bound state frequency given by Eq. ( 4) is in good agreement with the data (solid black line).The plot of the extracted ω BS,1 as a function of the bare qubit frequency is shown in Fig. 10(b), and the atomic population fraction is reported in Fig. 10(c), where the solid line is the calculated value with the parameters extracted from the fit.

Photon number estimation
We calibrate the photon number in the array modes using the mode Kerr nonlinearity inherited by the Josephson junctions array resonators.The self-Kerr coefficient of the bare resonator scales with the inverse square of the number of junctions, Moreover the nonlinearity is further diluted in each mode as with an effective Kerr K = K r /N = 2π × 100 kHz; in fact two excitations in one mode, in the momentum space, are spatially distributed to all the resonator of the array [51].We can experimentally confirm this estimation investigating the power response of the array modes.
Fig. 11 shows the phase and magnitude of a coherent microwave tone in the range of mode number 8 at 5.552 GHz.Adapting the treatment in [52] to a trans- mission configuration, and assuming a total input attenuation at the sample equal to 70 dB, we can globally fit the data to where κ is the external coupling rate, κ tot is the total loss rate, δ = (ω − ω 0 )/κ tot is the relative detuning, ξ = (P in / ω)κK/κ 3 tot is the (adimensional) input power, K is the mode Kerr and finally ñ is given by the solution of the algebraic equation The global fit (solid black lines) reproduces the data well.

System parameters
The complete list of the system parameters measured in our sample is shown in Table I.

Appendix B: Theoretical model 1. Circuit model
The ideal circuit model of the system is sketched in Fig. 12 and shows N LC-resonators, with capacitance C r and inductance L r , capacitively coupled in series via a capacitance C J .Two sites of the array, x 1 = 10 and x 2 = 12, are coupled to two transmon qubits, represented by a Cooper pair box with capacitance C q and Josephson energy E J , via the coupling C gi .
To start, we consider the bare coupled cavity array waveguide, whose unit cell has a lattice constant d = 200 µm.For probing frequency ω that respects the homogeneity condition d < λ g /4 for the electromagnetic waves, we can treat the coupled-cavity array as a composite right-left-handed transmission line.Defining the lattice unit impedance Z l (ω) = 1/iωC J and admittance Y l (ω) = (iωC r + 1/iωL r ), we can follow standard procedures [53] to obtain the dispersion relation: This equation can be recast in terms of the resonator frequency, ω r = 1/ Cr L r , with Cr = C r + 2C J , in the limit of C J Cr and reads where Z r = L r / Cr is the characteristic impedance of the resonator.From the dispersion relation Eq. (B1), we calculate the group velocity at the center of the band, which gives us a left handed dispersion with negative phase velocity [53].Finally, we observe that the characteristic impedance of the transmission line is equal to: After having characterized the waveguide properties we can take into consideration the full circuit model including the qubits.The Lagrangian of the circuit can be written in terms of the resonator and qubit fluxes, Φ x and Φ qi respectively [54], and reads where ΦT = ( Φq1 , Φq2 , Φ1 , ..., ΦN ) and is the potential energy term with Φ 0 = /(2e) being the flux quantum.The first term of (B5) is the kinetic term that is governed by the capacitance matrix (B7) where Cqi = C qi + C gi .The Hamiltonian of the system is obtained from the Lagrangian (B5) via the usual Legendre transformation [54] In the matrix notation the Hamiltonian reads where now the coupling between the charges are governed by the inverse of the capacitance matrix C −1 .This matrix when inverted connects all the elements of the circuit.Nevertheless, we can use the fact that the coupling capacitances are much smaller than the resonator ones, C J , C gi Cr , Cqi , to obtain the following simplified ex- with the definitions . In order to quantize the Hamiltonian Eq. (B8) we first express the charge and the flux of the resonators in terms of the annihilation and creation operators: In this way the canonical commutation relations The qubits can be described by nonlinear resonators with a Kerr nonlinearity β i and annihilation (creation) operators b (b † ).The qubit frequency is given by the frequency difference between the two lowest energy states, |0 and |1 , of the transmon and it is a function of the flux on the qubit, i.e. ω qi = ω qi (Φ qi ).The transition between the first two transmon levels is determined by the dipole moment D qi = 1|Q qi |0 and defines the qubit-cavity couplings g i = D qi Cr ω r /(2 C2 g ).Finally, by defining the cavity-cavity hopping J = Cr ω r /(2 CJ ) we obtain the system Hamiltonian given in Eq. (3).Note that in this derivation we neglected direct parasitic capacitive couplings, disorder in the circuit elements and we made some simplifications on the inverse capacitance matrix.We discuss the effect of these approximations in detail in App. C.

Bare coupled-cavity array in the tight binding picture
Here we re-discuss the bare array of coupled resonators starting from the standard tight-binding model with uniform nearest-neighbor couplings with ladder operators a x fulfilling usual bosonic commutation rules, [a x , a † x ] = δ x,x .It is worth stressing that here N must be finite and the array subject to open boundary conditions, in contrast to most treatments in the literature of bound states which focus on the thermodynamic limit and periodic boundary conditions.The array's free Hamiltonian (B11) is diagonalized as (see e.g.Ref. [55]) with the normal modes a k given by and the normal frequencies ω k by Note that the discrete spectrum (B14) obtained by a tight binding array with open boundary conditions coincides with the one obtained by the circuit model in Eq. (B2).

One-atom bound states
Consider the case that only one atom is effectively coupled to the coupled-cavity array (the other being fardetuned from the photonic band).Then the total Hamiltonian Eq. ( 3) in terms of normal modes (B13) and in a frame rotating at frequency ω r , reads (we omitted the label i in this case) where δ = ω q − ω r is the detuning of the qubit from the bare frequency of each resonator and x q the cavity to which the atom is directly coupled.Note the sine-shaped atom-mode interaction strength (colored coupling), stemming from the open boundary conditions which the array is subject to [56].As shown next, such a colored coupling generally could lead to results slightly different from the usual white coupling under periodic boundary conditions [16][17][18]56].Within the single-excitation subspace defined by n = 1, an atom-photon bound state |φ can be worked out as an eigenstate of the total Hamiltonian whose corresponding energy lies outside the photonic band, that is such that H|φ = ω|φ with |ω| > 2J.Expanding |φ in the basis {b † |g, 0 , {a † k |g, 0 }} as and then inserting this into the Schrödinger equation yields that the frequency ω must fulfills the equation where the self-energy Σ(ω) is given by The solutions of Eq.(B19) with |ω| > 2J correspond to the atom-photon bound states.In the standard case of periodic boundary conditions, two bound states -one with energy above and one below the band -always exist [17,18].In our case, however, the presence of array edges may affect the existence of bound states.By substituting in the eigenvalue Eq. (B19) the limiting values of the selfenergy Σ(ω) for ω → ±2J, the condition for existence of bound states is obtained as where − (+) indicates a solution with energy above (below) the band.In our experimental setup, g ∼ J, N = 21 and x q = 10, 12 so that Eq. (B22) is always satisfied and the self-energy is well approximated by performing the thermodynamic limit In practice, this means that the atom in our setup is sufficiently far from the edges and the number of resonators large enough that the bound state can be calculated as if the array were infinitely long (in line with standard treatments).Following Ref. [17], the bound state corresponding to a solution ω BS can thus be worked out in the form where θ is given by defines a bosonic ladder operator.Here, λ = λ(ω BS ) as given by (B21), while s = sgn(ω BS ) so that s = +1 when ω BS lies above the band and s = −1 when it falls below.The mixing angle θ measures the degree of hybridization: the dressed state is fully atomic for θ = 0 and fully photonic for θ = π/2.Eq. (B26) fully defines the (normalized) photonic component, showing that the spatial mode is exponentially localized around the atom's position x q .Accordingly, parameter λ represents the localization length of the photonic cloud surrounding the atom.

Two-atom bound states
Let us now study the case of bound states in the singleexcitation sector when two atoms are coupled to resonators x 1 and x 2 respectively.The total Hamiltonian is the natural generalization of (B16) and will now feature the detuning of each qubit δ i = ω qi − ω r and the corresponding coupling strength g i with i = 1, 2.
In this case, an effective interaction between bound states arises when the photonic clouds of the individual bound states overlap.This interaction changes the bound-state energies, which are now given by the real solutions with |ω| > 2J of the transcendental equation [17] (ω Here, Σ i is the self-energy of the ith qubit (in absence of the other qubit) as given by Eq. (B23) for g = g i .
Eq. (B27) admits up to four real solutions (and as many bound states).Based on the performed measurements in our experiment, we focus on the case of boundstate energies above the band.Then the pair of bound states with energies ω ± are given by where we defined the dressed ladder operators and the mixing angle with .
The parameter ξ = ξ(ω ± ) gives the amount of the hybridization between the two atoms and reads In the regime of equal coupling g 1 = g 2 (quite close to our experimental realization) and in correspondence with the avoided crossing, δ 1 = δ 2 , the two bound states get completely hybridized, ξ(ω ± ) = 1, and the labels +(−) signify the bound states with even (odd) symmetry with respect to the atoms' midpoint.Under certain conditions, the interaction can push the odd-parity bound state frequency into the propagating band (in which case the state just no longer exists).Similarly to (B22), an existence condition for this antisymmetric state can be worked out by taking the limit ω → 2J + in (B27), which yields (for The melting of this bound state into the photonic band reported in Fig. 6 indeed occurs for parameters such that (B32) is not satisfied.The interaction between one-atom bound states can cause also a (coherent) excitation transfer between the two atoms described by an effective spin Hamiltonian (the field is adiabatically eliminated).This occurs in the dispersive regime where the qubits are far-detuned from the band edge, i.e., ω qi − ω r − 2J g i , and cos θ ± ≈ 1 such that the bound states are mostly atomic.The effective spin Hamiltonian then reads [22,57] where we assumed that the atoms are tuned on resonance with one another, i.e., ω q1 = ω q2 = ω q , and set δ e = ω q − ω r −2J (detuning from the band edge).This Hamiltonian captures the excitation exchange dynamics between the qubits in Fig. 7.The coherence time is just given by the single atom bound state decay into the waveguide and into the other dissipation channels as discussed in the main text.

Two-photon bound states
We next study the energy of two-excitation bound states.In the ideal case β → −∞ (ideal two-level atom case), two-photon bound states are known to occur entailing strong nonlinear effects [17,18].Yet, even for finite β, the nonlinear energy spacing of levels above the first two can still affect two-excitation bound states causing measurable deviations from the fully linear case β = 0 (as discussed in the main text).To show this in more detail, we first note that any state in the two-excitation sector can be written as Here, c ij is the probability amplitude of having one excitation on transmon i and one on transmon j (including the case i = j), while c i (x) is the probability amplitude corresponding to one excitation on the i-th transmon and another one in the waveguide.Finally, u(x, y) = u(y, x) is the symmetric wavefunction of the two-photon bound state component.By plugging Eq. (B34) in the Schrödinger equation generated by Hamiltonian Eq. ( 3), we numerically solve the resulting set of coupled equations.
In the following, we separately address the one-and two-atom case.(one above and one below the continuous bands).For comparison, we also plot the bound state energies in the limiting cases of a linear resonator (black solid line) and a two-level atom (black dashed line).In Fig. 5 of the main text we defined the dressed-state anharmonicity as the difference in the bound state energy between the linear and nonlinear case.While Fig. 13(a) shows occurrence of a two-excitation bound state, it is natural to wonder how hybridized such a state is, and additionally, if it features a significant two-photon component (or alternatively if it is mostly populated by the excitation of the second transmon level).To clarify this point, in Fig. 13(b) we plot the population distribution of the two-excitation bound state, which clearly shows the photon dressed nature of the bound state in the considered parameter regime.Notably, the two-photon wavefunction exhibits a different localization length depending on the qubit frequency, in this respect similarly to the single-excitation case [see Fig. 13(c)-(d)].

b. Two atoms
As discussed in the main text, the two-excitation spectrum in the case of two atoms becomes quite involved, as shown in Fig. 14.Besides the two-photon unbound states with energies ω ∈ [2(ω r − 2J), 2(ω r + 2J)] (green box), two sidebands occur (above the main band) with energies ω ∈ [ω + − 2J, ω + + 2J] and ω ∈ [ω − − 2J, ω − + 2J], where the ± sign refers to the two-atom single-excitation bound state discussed in Section B 4. Above the bands, there appear three bound states stemming from the hybridization of the bare qubit states |20 , |02 and |11 .Once coupled to the array, the transmons get highly hybridized both with one another and the array field, resulting in the dressed states |φ 20 , |φ 02 and |φ 11 (the labels reflect the resemblance of each dressed state to the bare state of corresponding indices away from the avoided crossings).In particular, the color scale used in Fig. 14 reflects the population distribution among the two qubits, ranging from blue (excitation on Q 1 ) to red (Q 2 )

Array Transmission
We calculate the transmission through the array applying the input-output theory relation to the Heisenberg equation [58].The intra-cavities field a x for the x-th resonator, and b m for the m-th qubit, with a driving field a in (t) = a in e −iωt applied on resonator where κ = δ x1 κ r + δ x21 κ r + κ nr is the sum of radiative (in case of the edges cavities) and nonradiative decay rate for the resonators, γ is the decay rate for the qubits, and H is the Hamiltonian in Eq. (3).In the steady state, a x (t) = e −iωt a x and the equations for the resonators and qubits fields become where we defined the tuning ∆ r = ω r −ω and ∆ m = ω qi − ω.Solving the algebraic system and applying the inputoutput relation a out + a in = √ κ nr a x we calculate the transmission coefficient from cavity 1 to 21.

Appendix C: Experimental imperfections 1. Parasitic capacitance
As observed in the main text, the frequency distribution of the coupled-cavity array modes and the BSs interaction strength, cannot be quantitatively reproduced neglecting capacitive couplings beyond the nearest neighbor.
Solving Poisson's equation with a finite element method solver (Comsol Multiphysics, electrostatic package) we estimate a parasitic capacitance between next nearest neighbor resonators C (2) J ≈ 0.52 fF and the one between qubits and resonator C (2) g ≈ 0.73 fF.Adding this contribution to the capacitance matrix C, we obtain   where we introduced the next neighbor couplings g and J (2) = ω r Cr /(2 C(2) J ).

Interaction strength optimization
In Fig. 15 we report the interaction strength of different design based on the values of the device studied in this work.In particular, we calculate the expected bound state interaction for the same parameters of the current sample, but increasing the free sites between the coupling points of the qubits.In the new configuration the qubits are coupled to cavity 9 and 13, instead of 10 and 12.The same calculation is repeated, but reducing the qubitcavity coupling to 50 MHz.With the latter configuration, we expect an on-off ratio of 1000 within 1 GHz range.

Electromagnetic cross-talk
Four features of the transmission spectrum across the coupled-cavity array cannot be quantitatively explained by the ideal case.The transmitted signal reported in Fig. 16(c) (compare it with Fig. 2 plotted with a linear scale), highlights the non vanishing transmission outside the band, the "pairing" effect of modes, the transmission minima within the band and finally the asymmetric mode distribution.These non-idealities are due to the non-negligible cross-talk between the input and output ports of the sample box and the next nearest neighbor interaction discussed in the section C 1.
Fig. 1(a) reports a micrograph of the sample bonded to its sample holder.The bond wires connecting the sam-ple box ports to the signal launcher on chip are approximately 1.5 mm long.We model their electromagnetic cross-talk as two power dividers, the first of which separates the incoming signal S in into a part that reaches the second power divider after an electrical delay e iθ and a part through the coupled-cavity array, S in .In our model the two signals are then recombined by the second power divider.
The solid line in Fig. 16c shows the prediction for our model with the fitting parameters = 0.22 and θ = 0.34π.In comparison with a pure tight binding model represented by the dashed black line, we can reproduce the main features of the transmission spectrum.Unfortunately, the mode distribution is still not completely described by the model.We attribute this behavior to fabrication imperfection, mainly on the smallest resonator feature represented by the junctions in the resonators.The design of the JJ resonator intrinsically mitigates the possible variance by a factor of √ 10.Moreover, although the individual modes are affected by this imperfection, collective behavior of the coupled-cavity array as a waveguide is still close to the ideal because of the large resonator-resonator coupling J.
We would like to stress that the cross-talk affects only the measured transmitted field and not the intrinsic mode structure.

Magnetic cross-talk
The flux lines Z-Q 1 and Z-Q 2 have a linear magnetic cross-talk that has been calibrated and compensated.The flux reported here is therefore the net flux in each SQUID loop, and not the one produced by the aforementioned flux lines.
The relation between the flux in each SQUID, Φ q1 and Φ q2 , and the room temperature voltages applied to the coaxial cable, V 1 and V 2 respectively, is expressed by: where we introduced the inductance matrix for the flux lines L fl .From the measurement of the periodicity of the readout resonator of each qubit at different flux we can extract the coefficients of the inductance matrix and the flux offest in each SQUID loop.In order to decouple the two flux lines, we can redefine the voltages in eache of them as: For our sample we measure an offset Φ off 1 /Φ 0 = 0.091 and Φ off 2 /Φ 0 = 0.084, while the relative magnetic crosstalk for Q 1 is L 12 /L 11 = 0.041 and for Q 2 is L 21 /L 22 = 0.063.The model of the cross talk consists in two power splitter / combiner: a small portion of the signal bypass the waveguide and is recollected at the output, constructively or destructively interfering with the signal through the waveguide, and producing a non-zero transmission outside the band.(c) The measured transmission through the sample (green dots) is not well reproduced considering only the ideal case (input-output theory at low power, red dashed line).The sum on the signal through the array and the cross talk (black solid line) explain the non zero transmission and the "pairing" effect.

FIG. 2 .
FIG. 2. Single coherent tone spectroscopy of the coupled cavity array.(a) Frequency landscape of the resonant modes of the systems: Q 1 (blue), Q 2 (red) (tuned away from the band in this measurement), readout resonators (yellow), and array modes (green).b) Transmitted amplitude, |S21| vs frequency, ω, for different input powers, Pin.The traces are vertically offset in steps of 1.1 and filled to their own baseline for clarity.

FIG. 3 .
FIG. 3. Single coherent tone spectroscopy as a function of Q2 frequency.(a) Frequency landscape: Q 1 is tuned to its lowest frequency.The Q 2 frequency is tuned from its largest frequency at 6.6 GHz to 4.1 GHz.b) Low-power coherent tone spectroscopy through the resonators array as a function of the Q 2 frequency (dashed red line).(c) Calculated transmission spectrum via input-output theory using the single-excitation Hamiltonian Eq. (3).(d) Detail of the measurement in (a) showing how the atom-photon bound state approaches the band and eventually becoming the last array modes.The dashed black line shows a data fit via Eq.(4).(e) Total decay rate of the bound state as function of its frequency.The red data are obtained from the transmission spectroscopy while the blue from the atom spectroscopy (see Sec.II C).The dashed black line represents the theoretically expected decay rate.Insets: calculated distribution of the excitation over the array (green dots) and between the two qubits (blue and red respectively) for two distinct bound state frequencies: ωBS2 = 6.226GHz (e1) and ωBS2 = 6.510GHz (e2).

FIG. 4 .
FIG. 4. Isolated bound state in single excitation subspace.(a) Pulse scheme.A Gaussian pulse (red line) with frequency ωp drives Q 2 , while Q 1 is far detuned by flux pulse (blue shaded area).The readout pulse (in yellow), reads the population in both qubits.On the side, Q 1 and Q 2 (blue and red dots, respectively), the photonic band (green dots) and the readout resonators (yellow dots) frequencies are depicted in relation to the flux pulses.(b) Population of Q 2 as a function of flux pulse amplitude, Φq2, and driving pulse frequency, ωp.The black line shows the fit of the bound state frequency given in Eq. (4) as a function of the expected bare Q 2 frequency (red dashed line).(c) Bound state frequency extracted from panel (b).(d) Atomic fraction of the excitation.The black line shows the expected value (Eq.(B25)) using the parameters extracted from the fit in (b)-(c).
B 5) and to calculate the dressed anharmonicity, finding a good agreement with the measured data [Fig.5(b), solid line].

FIG. 5 .
FIG. 5. Isolated bound state in the two-excitation subspace.(a) Pulse scheme.A π-pulse (red line) excites Q 2 , a second Gaussian pulse brings it to the third level when ωp2 is resonant with the ef transition.When ωq2 is not resonant, a π-pulse brings the qubit back to the ground state.(b) Absolute value of the dressed anharmonicity β dress .The arrow indicates the value of the bare anharmonicity.The black solid line shows the calculated β dress with the parameters extracted from the fit in Fig. 4(b).Finally the black dashed line indicates the dressed anharmonicity of an ideal two-level atom.
FIG. 6. Avoided crossing and interaction strength between two bound states.(a) Pulse scheme.A Gaussian pulse (red line) with frequency ωp drives the two bound states tuned in resonance with two flux pulses (red and blue shaded areas).The populations of the qubits are read out (yellow line) at their highest frequencies.(b)-(g) Individual bound state population on Q 1 [P1 (b)-(d)] and on Q 2 [P2 (e)-(g)] as a function of the bare frequency of Q 2 for different frequencies of Q 1 .The dashed black lines are the expected bound state frequencies calculated from Hamiltonian (C3).The green dashed line indicates the photonic band edge.(h) Bound states coupling U , as defined in the main text, extracted from the measurements reported in (b)-(g).The black dashed line represents the expected value calculated from the Hamiltonian Eq. (3) while the solid line is obtained by the extended model given in Eq. (C3), which takes into account the next-nearest neighbor interaction.

FIG. 7 . 2
FIG. 7. Time resolved energy swap.(a) Pulse sequence.A πpulse (red line) excites Q 2 and after 30 ns the two bound state tuned in resonance with two flux pulses (blue and red shaded area) for a variable duration τ .The qubits are tuned back to their unbiased frequencies and their populations are read out (yellow line).(b)-(e)Chevron pattern.Measured bound state population of Q 1 (left column) and Q 2 (right column) as a function of the interaction time and of the bare frequency of Q 2 The two rows correspond to two different values of ωq1, at its largest frequency, (b)-(c), and 6.02 GHz, (d)-(e), respectively.(f)-(g) Horizontal line cuts in panels (b)-(e) (blue and red dashed line) displaying the population of Q 1 (blue) and Q 2 (red) as function of the interaction time.In both plots the relaxation time is fitted to an exponential decay (black dashed line).

FIG. 8 .
FIG. 8. Bound state interaction in two excitation subspace.(a) Pulse sequence.After tuning the bound states frequency with flux pulses (red and blue shaded area) they are excited (red and blue lines) and their population readout.(b) Cross Kerr interaction and avoided crossing between the states |02 and |11 for the frequency of the Q 2 bound state ω BS/2 /2π = 6.653GHz.The dashed line represent the expected value calculated with the Hamiltonian Eq. (C3).(c) Energy levels scheme for the bare (left) and dressed (right) states.

3 FIG. 9 .
FIG. 9. Complete wiring diagram and room temperature setup.The coherent tone to the waveguide is generated and measured with a Vector Network Analyzer (VNA).Q 1 and Q 2 driving (blue and red respectively) are made by analog upconversion of the pulsed generated by an arbitrary waveform generator (AWG).The flux lines are controlled by two AWG channels.Finally the qubits read-out (yellow) is controlled by thr up-and down-conversion of a pulsed generated by an AWG, and logged by an analog to digital converter (ADC).

6 Q1
FIG. 10.State preparation of the isolated Q 1 .(a) Population of Q 1 as a function of flux pulse amplitude, Φq1, and driving pulse frequency, ωp.The black line shows the fit of the bound state energy given in Eq. (4) to the measured data.The white dashed line shows the bare Q 1 frequency calculated with the parameters extracted from the fit.(b) Bound state frequency extracted from panel (a) as function of the bare Q 1 frequency.(c) Atomic fraction of the excitation, extracted from the relative population measured in the Qubit.The black line shows the expected value of the mixing angles as predicted from the ideal case theory given in Eq. (B25) using the parameters extracted from the fit in (a)-(b).

FIG. 11 .
FIG. 11.Single coherent tone spectroscopy of a single mode.(a) Magnitude and (b) phase of the transmission of the mode of the coupled-cavity array at 5.552 GHz as a function of power.The self-Kerr from the JJ array resonator is inherited by the modes.The coherent tone addresses only one cavity mode, so the global fit of a Kerr nonlinear resonator (continuous line) is used to extract the self-Kerr, K, and the photon number, n of the array mode.The shade areas indicate the deviation zero transmission (deviation from zero phase).

FIG. 12 .
FIG. 12. Circuit model for N = 21 capacitively coupled LCresonators (in green).The resonators at the edges of the array are capacitively coupled to 50 Ω transmission line.Two transmons qubits are coupled at sites x = 10 (blue) and x = 12 (red).
Fig. 13(a) shows the upper part of the two-excitation spectrum, ω > 2ω r , as a function of the atom frequency.The full spectrum (of which only the upper part is shown) features a band of two-photon unbound states defined by ω ∈ [2(ω r − 2J), 2(ω r + 2J)] plus a pair of sidebands ω ∈ [ω (1) BS − 2J, ω (1) BS + 2J], where ω (1) BS is the energy of the single-excitation bound states.Additionally, there exist a pair of two-photon bound states (see section II C 2) with discrete energies ω (2) BS such that |ω (2) BS − ω r | > 4J(one above and one below the continuous bands).For comparison, we also plot the bound state energies in the limiting cases of a linear resonator (black solid line) and a two-level atom (black dashed line).In Fig.5of the main text we defined the dressed-state anharmonicity as the difference in the bound state energy between the linear and nonlinear case.While Fig.13(a) shows occurrence of a two-excitation bound state, it is natural to wonder how hybridized such a state is, and additionally, if it features a significant two-photon component (or alternatively if it is mostly populated by the excitation of the second transmon level).To clarify this point, in Fig.13(b)we plot the population distribution of the two-excitation bound state, which clearly shows the photon dressed nature of the bound state in the considered parameter regime.Notably, the two-photon wavefunction exhibits a different localization length depending on the qubit frequency, in this respect similarly to the single-excitation case [see Fig.13(c)-(d)].

FIG. 14 .
photonic band C1)The inverse capacitance matrix, retaining only the first order terms in C gi , C J , C J , C possible to write the Hamiltonian:

FIG. 15 .
FIG.15.Interaction strength between two bound states.The four lines represent the calculated coupling between two bound states, for qubits coupled to different array sites(xi) and with different qubit-cavity coupling strengths g.

5 0
FIG.16.Transmission cross-talk.(a) Micrograph of the sample bonded in the copper sample box.The 1.5 mm aluminium bond wires connecting the inner conductor of a coaxial cable to the bond pad on the chips have a direct electromagnetic cross talk.(b) The model of the cross talk consists in two power splitter / combiner: a small portion of the signal bypass the waveguide and is recollected at the output, constructively or destructively interfering with the signal through the waveguide, and producing a non-zero transmission outside the band.(c) The measured transmission through the sample (green dots) is not well reproduced considering only the ideal case (input-output theory at low power, red dashed line).The sum on the signal through the array and the cross talk (black solid line) explain the non zero transmission and the "pairing" effect. GHz]

TABLE I .
Complete list of the system parameters.