A tunable Josephson platform to explore many-body quantum optics in circuit-QED

Coupling an isolated emitter to a single mode of the electromagnetic field is now routinely achieved and well understood. Current efforts aim to explore the coherent dynamics of emitters coupled to several electromagnetic modes (EM). freedom. Recently, ultrastrong coupling to a transmission line has been achieved where the emitter resonance broadens to a significant fraction of its frequency. In this work we gain significantly improved control over this regime. We do so by combining the simplicity of a transmon qubit and a bespoke EM environment with a high density of discrete modes, hosted inside a superconducting metamaterial. This produces a unique device in which the hybridisation between the qubit and up to 10 environmental modes can be monitored directly. Moreover the frequency and broadening of the qubit resonance can be tuned independently of each other in situ. We experimentally demonstrate that our device combines this tunability with ultrastrong coupling and a qubit nonlinearity comparable to the other relevant energy scales in the system. We also develop a quantitative theoretical description that does not contain any phenomenological parameters and that accurately takes into account vacuum fluctuations of our large scale quantum circuit in the regime of ultrastrong coupling and intermediate non-linearity. The demonstration of this new platform combined with a quantitative modelling brings closer the prospect of experimentally studying many-body effects in quantum optics. A limitation of the current device is the intermediate nonlinearity of the qubit. Pushing it further will induce fully developed many-body effects, such as a giant Lamb shift or nonclassical states of multimode optical fields. Observing such effects would establish interesting links between quantum optics and the physics of quantum impurities.

The interaction between light and matter remains a central topic in modern physics despite decades of intensive research. Coupling an isolated emitter to a single mode of the electromagnetic field is now routinely achieved in the laboratory 1 , and standard quantum optics provides a complete toolbox for describing such a setup. Current efforts aim to go further and explore the coherent dynamics of systems containing an emitter coupled to several electromagnetic degrees of freedom 2-6 . Recently, ultrastrong coupling to a transmission line has been achieved where the emitter resonance broadens to a significant fraction of its frequency, and hybridizes with a continuum of electromagnetic (EM) modes 7,8 . In this work we gain significantly improved control over this regime. We do so by combining the simplicity and robustness of a transmon qubit and a bespoke EM environment 6,9,10 with a high density of discrete modes, hosted inside a superconducting metamaterial 11 . This produces a unique device in which the hybridisation between the qubit and many modes (up to ten in the current device) of its environment can be monitored directly. Moreover the frequency and broadening of the qubit resonance can be tuned independently of each other in situ. We experimentally demonstrate that our device combines this tunability with ultrastrong coupling 12-15 and a qubit nonlinearity comparable to the other relevant energy scales in the system. We also develop a quantitative theoretical description that does not contain any phenomenological parameters and that accurately takes into account vacuum fluctuations of our large scale quantum circuit in the regime of ultrastrong coupling and intermediate non-linearity. The demonstration of this new platform combined with a quantitative modelling brings closer the prospect of experimentally studying many-body effects in quantum optics. A limitation of the current device is the intermediate nonlinearity of the qubit. Pushing it Introduction Due to strong interactions between elementary constituants, correlated solids 23 and trapped cold atoms 24 host fascinating many-body phenomena. Attempts to produce similar effects in purely optical systems are hampered by the obvious fact that photons do not naturally interact with each other. If this obstacle can be overcome, there is the tantalizing prospect of probing the many-body problem using the contents of the quantum optics toolbox, such as single photon sources and detectors, high-order correlations in time-resolved measurements, entanglement measures, and phase space tomographies to name a few.
One route to building a many-body quantum optical system is to rely on arrays of strongly non-linear cavities or resonators 22,25,26 , but minimising disorder in such architectures is a formidable challenge. Another route that circumvents these difficulties involves coupling a single well-controlled non-linear element to a disorder free harmonic environment. If the difficult experimental challenge of engineering an ultra-strong coupling can be overcome, thus exceeding the boundaries of the standard single photon regime in quantum optics, this approach could pave the way to bosonic realizations of electronic impurity systems such as the famous Kondo and Anderson models 21 . Our goal here is to achieve a large coupling between a sufficiently non-linear qubit and a quantum coherent environment containing many harmonic degrees of freedom.
When coupling an impurity to a finite size electromagnetic environment, five important frequency scales have to be considered. The impurity is characterized by its qubit frequency ω qubit , i.e. the excitation frequency between its two lowest internal states. Real impurities always possess more than two levels. The anharmonicity α, defined as the difference between ω qubit and the frequency for excitation from the second to third internal impurity state, characterises the departure of an impurity from a trivial harmonic oscillator (α → 0) or a pure twolevel system (α → ∞). The coupling between the impurity and environment is characterized by the spontaneous emission rate Γ at which the impurity exchanges energy with its environment. The environment itself is characterized by its free spectral range δω, which measures the typical frequency spacing between environmental modes, and the spectral broadening κ of these modes due to their coupling to uncontrolled degrees of freedom. The soughtafter multi-mode regime is obtained when Γ is larger than δω so that the impurity is always coupled to several discrete environmental modes, producing a cluster of hybridized qubit-environment resonances 27 . There are several requirements for reaching the many-body regime. First, Γ must be a significant fraction of ω qubit (ultrastrong coupling). This is a prerequisite for multiparticle decay 17,20 . If the coupling is too weak, the system becomes trivial, since only number-conserving processes are relevant (Equivalently Markov and rotating wave approximations apply.) A second requirement is α Γ. If this condition is not met, the non-linearity of the impurity is swamped by the broadening of the impurity levels, and the same frequency will drive transitions between several impurity levels, so that the system as a whole behaves more like a set of coupled harmonic oscillators than like a two-level system coupled to an environment 28 . Within the many-body regime, two limits can be distinguished. In the case of a finite-size environment that we address here (namely, δω > κ), each mode of the system can be addressed and controlled individually, while in the limit of a thermodynamically large environment (δω/κ → 0) one recovers a smooth dissipation-broadened qubit resonance.
The many-body ultra-strong coupling regime (defined by the first two conditions above) is hard to reach in quantum optics experiments because the coupling to three-dimensional vacuum fluctuations arises at order [α QED ] 3 , with α QED 1/137 the fine structure constant 1 . However, for superconducting qubits coupled to transmission lines, the scaling is much more favorable 18,29-31 than in a vacuum. Indeed the ratio Γ/ω qubit can essentially be made arbitrarily large, provided the impedance of the environment matches that of the qubit (see Sec. E of the Supplementary Information). Building on this ability of superconducting circuits to reach very large couplings, several experiments demonstrated the ultra-strong coupling regime in coupled qubit/cavity systems [12][13][14]32 . The rich physics associated to this coupling regime has also been evidenced using quantum simulation 33,34 . The condition Γ > δω has also been fulfilled by coupling superconducting qubits to open transmission lines [35][36][37] or engineered resonators 2,4 . However, it is only recently that the necessary conditions for the many-body regime were demonstrated concurrently 7,8 . The device of Refs. 7 and 8 consists of a flux qubit coupled to the continuum provided by a superconducting transmission line, which realises the thermodynamic limit (δω/κ → 0). Limitations of such setups include the lack of a microscopic model (since it is hard to characterize the waveguide properties of a transmission line outside the relatively narrow 4-8 GHz band where microwave transmission experiments can comfortably be performed), and importantly, that the transmission line is not an in-situ tunable environment.

Results/Discussion
In this work, we circumvent the above limitations, by designing circuits that provide independent tunability of both a qubit and a finite size but very large environ-ment, while allowing high-precision spectroscopic measurements of the environment itself (δω > κ) and first principle modeling. Our device, shown in Figure 1b, consists of a transmon qubit, which is relatively insensitive to both charge and flux noise, capacitively coupled to a long one-dimensional Josephson metamaterial, comprising 4700 SQUIDs. Such chains have been studied since the early 90's in the context of the superconductorinsulator transition 38-41 or to explore dual of the Josephson effect 42-44 . Our setup differs in two ways from these previous works. First, we took great care to produce a chain in the linear regime, far from the onset of nonlinear effects such as quantum phase slips 45 , so that one of the basic benifits of quantum optics, i.e. the elimination of non-linearities where they are not wanted, is realized. Second, we performed AC microwave spectroscopy of our device, instead of DC transport measurements. This allows us to characterize the electromagnetic degrees of freedom, also called Mooij-Schön plasma modes 28,46-50 , microscopically. We managed to resolve as many as 50 individual low frequency electromagnetic modes of this non-dissipative and fully tunable environment (See Figure 2c). An essential property of the Josephson metamaterial is its high characteristic impedance Z c = L J /C g 1590 Ω, which being of the same order of magnitude as the effective impedance of our transmon qubit, 760 Ω (here at zero flux), allow us to reach multi-mode ultrastrong coupling. The simplicity of the transmon architecture enables us to either compute from first principles or to extract all the parameters necessary to construct a microscopic model of the full system, without dropping the so-called 'A 2 -terms', a routine approximation in optics 51-55 , that however breaks down at ultra-strong coupling.
Our measurements are based on the frequency-resolved microwave transmission through the whole device. Figure 2b shows the amplitude of the transmitted field at a fixed value of the external magnetic field, and at low probe power. This spectrum reveals a set of resonances in our device, displaying a narrow spectral broadening κ/2π = 20MHz (for the non hybridized modes of the chain) caused by the coupling to the 50 Ω contacts. As the external magnetic field is varied, two modulation periods are seen in the resonance spectrum ( Figure 2a). The short and long periods correspond respectively to a one quantum Φ 0 = h/2e increase of the flux through the large transmon or through the small chain SQUID loops. This feature allows us to adjust independently the flux threading the transmon SQUID loop (Φ T ) from the one threading the chain SQUID loops (Φ C ). The former controls the qubit frequency, while the latter controls the impedance of the environment, and hence the qubit-environment coupling strength. Before studying the hybridization between the transmon and chain, we characterize the chain on its own by setting Φ T = Φ 0 /2, so that the qubit decouples from it (See Fig. 2c and Sec. B of the Supplementary Information). We obtain a good A Josephson platform for waveguide quantum electrodynamics. a Lumped element model of the circuit including the nodes used in the calculations. b Optical microscope image of the sample. The two zoom-in are SEM pictures of the SQUID of the qubit (red square) and the SQUIDs in the chain (blue square). The qubit is capacitively coupled to the chain and to a 50Ω measurement line, via large interdigital contacts. Only a small portion of the Josephson chain, which comprises 4700 SQUIDs in total, is shown. fit between the dispersion relation predicted by a microscopic model and the one extracted from the measured resonances. This allows us to extract all of the parameters necessary to characterize the chain modes.
The transmon qubit becomes active when Φ T = Φ 0 /2, and is expected to hybridize with the chain. Figure 3a shows a low-power spectroscopy of the system as a function of Φ T , keeping Φ C nearly constant. When a chain mode is not hybridized with the qubit, the corresponding spectral line runs nearly horizontally. When Φ T is varied, the qubit frequency sweeps across the resonances of the chain modes and creates a clear pattern of several avoided level crossings (See Figure 3a). We note that at fixed Φ T , several chain modes in the vicinity of the transmon resonance are visibly displaced. Thus, the transmon simultaneously hybridizes with many modes. This is a signature of multi-mode ultra-strong coupling, a topic that will be further addressed below. Evidence that the transmon behaves as a qubit is provided by its saturation spectrum (Figure 3b). Here the fluxes Φ T and Φ C are kept constant, while the transmission through the system is recorded at increasing probe power. For a harmonic system, the resonance positions are independent of driving power. In a anharmonically oscillating classical system, a gradual dependence on the driving power would appear. Experimentally, we observe that below a driving power ∼ −10 dBm, the resonance positions in the transmission spectrum are independent of the driv-ing probe power. As the driving power increases beyond −10 dBm, the peak around 4.6 GHz disappears while the other peaks assume the positions they have when the qubit is inactive. (The horizontal lines correspond to the low power spectrum when Φ T = Φ 0 /2.) This is evidence of saturation, a clear qubit signature 16,35,36 . The fact that several chain modes are shifted when the transmon is saturated constitutes an additional proof that the transmon hybridizes with many modes at once.
To quantify the hybridization of the transmon mode with the chain modes, we compare the normal mode spectrum of the full system at Φ T = Φ 0 /2 to the spectrum at Φ T = Φ 0 /2. As mentioned above, in the latter case, the transmon decouples from the bare modes of the chain. The system with Φ T = Φ 0 /2 therefore has one extra mode in the vicinity of the transmon frequency. We define the relative frequency shift δφ n (Φ T , Φ C ) as the difference in frequency between the nth mode of the coupled and uncoupled chains, normalized to the free spectral range of the chain δω n (Φ T , Φ C ) = ω n (Φ 0 /2, Φ C )−ω n−1 (Φ 0 /2, Φ C ), and including a π factor for later convenience, i.e.
where ω n (Φ T , Φ C ) is the frequency of the nth lowest nonzero mode for a given flux in the transmon and in the chain. This frequency shift is readily extracted from the Remarkably (see Sec. H of the Supplementary Information for a derivation), δφ n (Φ T , Φ C ) in Eq. (1) equals the phase shift experienced by mode n due to the presence of the nearby transmon mode: where φ n (Φ T , Φ C ) is the phase shift of mode n of the full system at transmon flux Φ T and chain flux Φ C . From Eq. (1) it follows that the phase shift equals 0 (resp. π) for modes far below (resp. far above) the renormalized transmon frequency. For hybridized modes in the vicinity of the transmon line, δφ n (Φ T , Φ C ) lies between 0 and π. This behavior is clearly observed in Fig. 4a where the measured relative frequency shifts are reported for a chain flux Φ C = 0 and various transmon fluxes Φ T . The wide frequency dispersion of intermediate δφ n (Φ T , Φ C ) provides direct evidence for a hybridization with up to ten chain modes. In the thermodynamic limit of an infinite chain with perfect impedance matching to the measurement ports, the transmon-induced phase shift δφ n (Φ T , Φ C ) becomes a continuous function δφ(ω, Φ T , Φ C ) of the mode frequency ω. Moreover, it can be shown that the frequency derivative of δφ(ω, Φ T , Φ C ) matches very precisely the theoretically expected lineshape of the dissipative response of the transmon coupled to an infinite environment. (See Sec. I of the Supplementary Information.) This constitutes a central finding of our work: the renormalized transmon frequency ω T and linewidth Γ T can be directly inferred from a measurement of the phase shifts of the individual modes in the finite bath. In terms of measurement protocol however, there is a sharp difference between the chain mode phase shifts and the qubit response functions. Usually, the qubit response is obtained by observing the transmon, and its environment can be viewed as a black box that combines unmonitored decoherence channels as well as the physical ports used for measurement. This procedure constitutes the usual paradigm in the study of open quantum systems. Our protocol is unusual because information about an open quantum system is obtained by monitoring the discrete modes that constitute its dominant environment. We finally turn to a quantitative analysis of our data, including a comparison to the predictions of a microscopic model, in order to determine if the requirements for reaching the many-body regime has been met. By extracting the maximum renormalized transmon frequency ω T,max extracted from the phase shift data of Fig. 4a at chain flux Φ T an integer multiples of Φ 0 , we are able to infer the only remaining unknown system parameter, namely the maximum transmon Josehpson energy E J,T,max . This allows us to estimate the anharmonicity of the transmon α, which ranges from 0.36 GHz at Φ T = 0 to 0.44 GHz at Φ T = 0.3Φ 0 . We emphasize that the condition α Γ T for anharmonic many-body behav-ior is thus fulfilled (see Fig. 4b for the extracted transmon linewidth Γ T , which lies in the range 0.2-0.4GHz). Using the extracted parameters to calculate δφ(ω, Φ T , Φ C ) according to Eq. (11) we find the predicted theoretical lines in Fig. 4a. The excellent agreement between theory and experiment seen here for six different values of Φ T persists for each of the hundreds of (Φ T ,Φ C ) combinations where we have made the comparison. (See Sec. C of the Supplementary Material for a further selection of results.) We stress that this agreement is obtained after all model parameters have been fixed, so that there is no fitting involved in comparing the predicted and measured phase shifts. The quantitative modeling of such a large quantum circuit clearly is an important landmark in the field of open quantum systems. In Fig. 4b we examine the transmon linewidth Γ T that we extracted from the phase shift data, as a function of chain flux Φ C for fixed transmon flux Φ T = 0. Very good agreement (with no fitting parameters) is again obtained with the prediction of our model. These results demonstrate that we can tune the qubit-environment coupling independently from ω T using the flux in the chain, and that we achieved the ultra-strong coupling in our waveguide, i.e. coupling to a large number (here 10) of modes with a sizeable linewidth Γ T /ω T 0.1. A hallmark of ultra-strong coupling is the failure of the rotating wave approximation (RWA), as previously discussed in coupled qubit and cavity systems 12 . We have examined the consequences of the RWA on our microscopic model (see Sec. J of the Supplementary Information), and found a discrepancy of 100MHz in the transmon frequency ω T , showing the quantitative importance of non-RWA terms. We would like to stress that demonstrating the relevance of these terms is much more than obtaining a good data-theory agreement. With counter-rotating contributions in the few percent range, we expect a finite rate for parametric processes in which photon-number is not conserved. In future we plan to use the current platform to observe these interesting many-body effects directly.
In conclusion, this work provides the first demonstration of many-body ultra-strong coupling between a transmon qubit and a large and tunable bath. To obtain full control over the environment, a superconducting metamaterial, comprising 4700 SQUIDs, was employed. Although this quantum circuit contains a huge number of degrees of freedom, we were able to characterise all its parameters in-situ. This allows us to demonstrate unambiguously that our systems meets the three conditions required to reach the many-body regime, namely Γ δω, Γ 0.1ω qubit and α Γ. A novel experimental methodology was implemented to analyze the qubit properties by means of the extraction of the phase shifts of the environmental modes. Despite the large size of our quantum circuit, we succeeded in providing a fully microscopic model which accounts for the transmon response without any fitting parameters. We also found that the qubit linewidth for the long chain agreed with results in the thermodynamic limit, showing that the Extraction of qubit properties from the measurement of its controlled environment a Phase shift δφn of the discrete chain modes as a function of mode frequency ωn for different transmon fluxes ΦT, fixing ΦC = 0. The solid line is a fit using Eq. (2). The inset shows the chosen transmon fluxes ΦT as line cuts in the transmission measurement (with the same color code). b The transmon width ΓT for different fluxes in the chain ΦC, showing control of the coupling to the large but discrete environment. The experimental points (dots) are obtained from an arctangent fit of the data in panel a. For better visibility, only the flux values where the transmon frequency is maximum are included. The blue shaded area represents the theoretical expectation for ΓT, within a confidence interval given by the error in the capacitances in Table I. finite environment has the same influence on the qubit as a truly macroscopic bath. The further possibility to tune the coupling to the environment in-situ, demonstrated by a 50% flux-modulation of the qubit linewidth, opens the way to controlled quantum optics experiments where many-body effects are fully-developped 16,17,19,20,56,57 , as well as more advanced environmental engineering for superconducting qubits 25 .

Sample fabrication and parameters
The sample was fabricated on a highly resistive silicon substrate, using a microstrip geometry. The ground is defined as the back-side of the wafer which was gold-plated, ensuring a good electrical conductivity. Interdigital capacitances were chosen to connect the transmon and the metamaterial. They do not provide the lowest surface participation factor 58,59 but they allow us to maximize the coupling capacitances Cc of the transmon to the chain, while minimizing the capacitances to ground C g,T and C g,T2 . (See panel a of Fig. 1 for the definitions of the these capacitances.) This system is probed via two 50 Ω transmission lines, one of which is capacitively coupled to the transmon, while the other is galvanically coupled to the chain. The whole device (Josephson junctions, capacitances and transmission lines) was fabricated in a single electron-beam lithography step, using a bridge-free fabrication technique 60 . The Josephson elements of the chain are tailored to be deep in the linear regime (E J /E C = 8400), where E J and E C are the respectively the Josephson energy at Φ C = 0 and the charging energy of a chain element), leaving the transmon as the main source of non-linearity in the system. All parameters of the system are listed in Table I. Here, L J,min is the minimum inductance of a chain SQUID loop, (which occurs when Φ C = 0). There is a slight asymmetry between the two Josephson junctions that constitute a single chain SQUID loop, which is quantified by the asymmetry parameter d. The flux-dependent inductance L J (Φ C ) of a chain SQUID loop is given by and the corresponding Josephson energy by E J (Φ C ) = ϕ 2 0 /L J (Φ C ) with ϕ 0 = /2e the reduced flux quantum. The transmon Josephson junctions are symmetric, and hence the flux dependent transmon Josephson energy is given by The transmon charging energy is not an independent parameter (see Eq. (9)), but is listed in the table, due to the prominent role it plays in what follows. The meanings of the remaining parameters in Table I  The only exceptions are C J , that is obtained from knowledge of the junction areas via an empirical formula, and E J,T,max that we extract via a procedure described in the last subsection below, which uses the data at Φ T equal to multiples of Φ 0 in panel a of Fig. 2.

Full Model
The circuit diagram for the lumped-element model is shown in Fig. 1. It consists of N + 2 nodes, where N is the number of SQUIDs in the chain. To describe the circuit, we use the Cooper pair number operator n j , which gives the number of Cooper pairs in node j, and the superconducting phase operator ϕ j , which gives the superconducting phase at node j. They satisfy the canonical commutation relations [ n j , ϕ l ] = iδ jl 61 . Here j, l ∈ [L, R, 1, 2, 3, . . . , N ] with L and R referring to the left and right transmon nodes. As explained before, the SQUIDs of the chain are linear inductors, to a very good approximation. We define n T = ( n L , n R , n 1 , . . . , n N ) and ϕ T = ( ϕ L , ϕ R , ϕ 1 , . . . , ϕ N ). In this notation, the Hamiltonian of the circuit is given by C is the capacitance matrix, such that elements [ C] jl = [ C] lj equal the capacitive coupling between the charges on islands j and l.
In the same way, elements [ J] jl = [ J] lj of matrix J contains the Josephson energy that couples the superconducting phase on island j and island l. Both matrices are (N + 2) × (N + 2). Their explicit forms are given below in Eqs. (6) and (7). In both matrices, the boundary conditions that determine entries (L, L) and (N, N ) are obtained by assuming that the nodes to the left of node L and to the right of node N are grounded.
The elements in the capacitance matrix are given by We define the operators n T = ( n R − n L )/2 + constant and ϕ T = ( ϕ R − ϕ L )/2 associated with the transmon dynamics. Introducing these operators, and noting that the total transmon charge n R + n L is concerved, we can rewrite the Hamiltonian as where we defined ϕ N +1 ≡ 0. The transmon charging energy E C,T is given by The coupling of n T to the charge on island j is given by Chain modes phase shift in the thermodynamic limit In Fig. 4a we compare the measured relative frequency shift δφn(Φ T , Φ C ) to the theoretically predicted transmon phase shift δφ(ω, Φ T , Φ C ) with which it is expected to agree in the thermodynamic limit. Here we provide the analytical formula for the phase shift φ(ω, Φ T , Φ C ) of a mode with frequency ω. (See Sec. G of the Supplementary Information for the derivation.) It reads In this expression ωp(Φ C ) = 1/ L J (Φ C ) (C J + Cg/4) is the plasma frequency of the chain, and . (12) has dimensions of capacitance. Finally, E S (Φ T ) is an effective linear inductor energy associated with the Josephson junctions in the transmon, which nonetheless incorporates the transmon nonlinearity, and is given by (See Sec. F of the Supplementary Material for further detail.) The theoretical δφ(ω, Φ T , Φ C ) curves plotted in panel a of Fig. 4 were obtained from φ(ω, Φ T , Φ C ) similarly to Eq.
(2) as the difference Analysis of the experimental data To extract the relative frequency shift δφn(Φ T , Φ C ) from the data presented in Fig. 2 of the main text, we go about as follows. At a fixed value of the magnetic field that determines Φ T and Φ C , we fit each of the peaks in the transmission spectrum individually with a Lorentzian. This gives the center frequency of the peaks. From these peak positions, we obtain δφn(Φ T , Φ C ) experimentally using Eq.
(1) at a particular Φ T and Φ C . Next we extract the transmon frequency ω T for the transmon coupled to the chain from the experimentally determined δφn(Φ T , Φ C ). The details are as follows. Empirically, we find that the experimentally determined δφn(Φ T , Φ C ) vs. ωn data-points fit an arctangent lineshape very well, for suitable choices of the parameters A, ω T , and Γ T . (In the parameter regime where our device operates, the theoretically predicted phase shift δφ(ω, Φ T , Φ C ) also closely approximates this line shape.) We therefore fit the measured δφn vs. ωn at fixed Φ T and Φ C to Eq. (15), interpreting ω T as the frequency and Γ T as the resonance width of the transmon when it is coupled to the chain. Before we can quantitatively compare the experimental results for δφn(Φ T , Φ C ) to the theoretically predicted δφ(ω, Φ T , Φ C ) (Eq. (11)), one final model parameter, namely the maximum transmon Josephson energy E J,T,max must be determined from the experimental data. The general procedure is as follows. Our theoretical model predicts that in the regime where the actual device operates, this transmon frequency is very nearly equal to the isolated (L J → ∞) transmon frequency, i.e. the chain only slightly renormalizes the transmon frequency, and indeed, we see little Φ C dependence in the extracted ω T . At Φ T = nΦ 0 , n = 0, ±1, ±2, . . ., where the transmon Josephson energy is maximal, we therefore use the isolated transmon result (see Sec. F of the Supplementary Information): Taking the average over n of the experimentally determined ω T (Φ T = nΦ 0 ), and using our first principle estimate for E C,T in Table I, we obtain The theoretical curves in Fig. 4a were then obtained using the system parameters in Table I in Eqs. (11), (14), (13), and (4). The full data set covers many transmon periods. Within a given transmon period, we generally analyze data at several values of Φ T in the interval from -0.3 φ 0 to 0.3 φ 0 . Each transmon period is measured at different Φ C . We also take into account the small variation in Φ C as the transmon flux sweeps through one flux quantum. The experimental points in Fig. 4b are obtained as the transmon width Γ T closest to Φ T = 0. The error bars come from the least square fit using Eq. (15). The theoretical width is obtained from a fit of the phase shift δφ(ω, Φ T , Φ C ) with the arctangent of Eq. (15).

Data Availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.   Supplementary information for "A tunable Josephson platform to explore many-body quantum optics in circuit-QED"

A. Experimental setup
The full measurement setup is shown in Fig. S1. The device was placed in a dilution refrigerator at a base temperature of 20 mK, and the transmission measurements were performed using a Vector Network Analyzer (VNA). An additional microwave source was used for two-tone measurements, while a global magnetic field was applied via an external superconducting coil. Both the coil and the sample were held inside a mu-metal magnetic shield which is coated on the inside with a light absorber made out of epoxy loaded with silicon and carbon powder. The output line included two isolators at 20 mK, a HEMT amplifier at 4 K and a room temperature amplifier. The input line is attenuated at various stages, including a home-made filter that prevents stray-radiations from reaching the sample. We adopted a coaxial geometry with a dissipative dielectric (reference RS-4050 from resin systems company). The bandwidth of the measurement setup goes from 2.5 GHz to 13 GHz.

B. Chain dispersion relation
In this section we explain how we experimentally obtained the dispersion relation of the chain, and how we used it to determine the model parameters L J,min , C g and C J that characterize the chain. (See Methods -Sample fabrication and parameters in the main text.) The experimental data in Fig. S2a (reproduced from Fig. 2b in main text) was obtained by first tuning the external magnetic field to a point where Φ T = Φ 0 /2 so that E J,T (Φ T ) = 0. (Green dashed line in the inset to Fig. S2a.) This leads to vanishingly low transmon frequency. As a result, the transmon does not contribute any degrees of freedom that can hybridize with the bare modes of the chain. In order to realize a good fit, one needs to measure the spectrum in a wide frequency band (0.1 GHz to 20 GHz). One is however limited by the bandwidth of the setup (2.5 GHz to 13 GHz). This difficulty can be overcome by performing a two-tone measurement 1-3 , taking advantage of the fact that the array is not perfectly linear. As a consequence, when applying a microwave tone at a given resonance of the chain, the other resonant frequencies are shifted by the cross Kerr effect. With the Vector Network Analyzer (VNA), we proceed by measuring the transmission of the system at a fixed frequency ω VNA = ω 1 where ω 1 matches a given resonance frequency of the circuit. Then with a microwave source we apply a second tone at a variable frequency ω MW . Whenever ω MW equals any other resonance frequency of the circuit, ω 1 shifts toω 1 due to the cross Kerr effect, so that ω VNA =ω 1 , which leads to a dip in transmission. The value of ω MW at these dips provides all the resonances of the system. Because we measure at a constant frequency ω VNA inside the setup bandwidth, we are not limited by the frequency range of our measurement setup anymore. A typical two-tone measurement is shown in Fig. S2b, where we fit each dip separately with a Lorentzian. The center frequencies obtained from these fits are the experimental points in Fig. S2a for the eigenmodes of the chain.
In order to fit the experimental data for the chain modes, we assume that the left end of the chain is open when Φ T = 0.5 Φ 0 (E J,T = 0). We also take the right end of the chain to be grounded. Given that the chain SQUIDS are designed to have a Josephson energy several thousand times their charging energy (E J /E C = 8400), we can model the Josephson junctions in the chain as linear inductors with inductance L J (Φ C ) = ϕ 2 0 /E J (Φ C ), with ϕ 0 = /2e the reduced flux quantum. The theoretical dispersion relation can be obtained applying Kirchoff's laws to one chain cell  Fig. 1a in the main text.) Denoting the flux at node j as Φ j , we obtain Now if we use as ansatz a plane waves Φ j = A exp i (ωt − κja) + B exp i (−ωt + κja) and solve for ω we obtain the dispersion relation for a bare chain The boundary conditions at site 0 (vacuum) and at site N (grounded) read which restricts the values of κa to κa = n − 1 2 π N n = 1, 2, . . . , N.
Since the areas of chain SQUID loops are much smaller than that of the transmon SQUID loop, at Φ T = Φ 0 /2 we can tune the flux Φ C through each chain SQUID to a multiple of Φ 0 , so that the chain SQUID inductance is minimal, i.e. L J (Φ C ) = L J,min , without appreciably changing the transmon flux Φ T from its value Φ 0 /2.  The error for C J is just the error we obtain for the measurement of the Josephson junction's area using a Scanning Electron Microscope. The error for L J,min and C g are the values where the deviation between the experiment and the fit was below 5 %.

C. Additional phase shift data
In this section we present a further selection of relative phase shift data δφ n obtained for various (Φ T ,Φ C ) combination. This is only a small subset of the full data set, and the agreement between theory and experiment exhibited here is representative of the full data set. Results presented here complement Fig. 4a of the main text. The parameters used to obtain the theory curves are the ones of Table S1.

D. Transmon qubit capacitances estimation
In order to obtain the capacitances listed in Table S1 (reproduced from the main text) we use EM simulation software (Sonnet). This software solves Maxwell's equations in three dimensions for the specified design of our device and gives the scattering parameters of the system as a function of frequency. We simulate two parts of the design independently, the interdigital capacitors and the SQUID of the transmon qubit.

Interdigital capacitors
Since we are only interested in modelling the capacitors of the transmon, we remove the chain from the simulation and replace the Josephson junction of the transmon by a linear inductor, L test . We place two ports at both ends of the design, and set the characteristic impedance of the port on the left to Z left and for the port on the right to Z right .
From the EM simulation we obtain the transmission of the system, S 21 , as a function of frequency. We fit the prediction of the linear model of the qubit to this. This model consists of the capacitance network shown in Fig. S4 in red. The transmission of this system is given by Eq. (S7) 5 where A, B, C and D are the ABCD matrix elements 6 for the capacitance network plus the linear inductor.
In theory L test , Z left , Z right do not affect the obtained capacitances and can be chosen arbitrarily. However, due to the fact that the lumped element model is an idealization, we observed a small shift of the capacitances as a function of L test . (This shift was not observed as a function of Z left or Z right ). To minimize the effect of this shift, we set L test = 22 nH which gives a resonance frequency close to ω T .
We perform two simulations. In the first one we set Z left = Z right = 50 Ω. Due to the low impedance of the ports, we can neglect C g,T2 and we therefore fit only C c , C g,T,0 and C sh,0 . Then we perform a second simulation with Z left = 50 Ω and Z right = 3000 Ω. Now we fit only C g,T2 keeping the other capacitances constant. In this way we obtain all the capacitances in Fig. S4 independently. Note that the self-capacitance of the junction C J,T cannot be simulated and is therefore obtained from Eq. (S6). The errors are obtained as the maximum range where the difference between simulation and model is smaller than 10 %. The two fits are shown in Fig. S5.

Stray capacitances from the transmon SQUID
The transmon qubit has a SQUID with a large loop (∼ 55 µm × 1.2 µm). Due to its large size, the capacitances associated to this SQUID are not negligible. In Fig. S6a the SQUID design with the different capacitances is given. The lumped element model used for simulating the system is shown in Fig. S6b.
We follow the same procedure as before. Given the small number of fitting parameters (C sh,S and C g,S ) we can perform a single fit with Z left = 50 Ω and Z right = 3000 Ω. The result of the fit is shown in Fig. S7 with the obtained capacitance values. The SQUID increases both the shunting capacitance and the ground capacitance of the transmon qubit. Before quantitatively modeling the system, in this section we try to gain a qualitative understanding of how the transmon decay rate Γ T depends on the characteristic impedance of the chain. Since we are only aiming for a qualitative description, we treat the transmon SQUID loop as an LC circuit, ignoring its non-linearity. We retain the capacitive couplings C c , that couple the transmon to the chain and the 50 Ω transmission line. We drop the ground capacitances C g,T and C g,T2 that shunt the chain at high frequencies, thus idealizing to the situation where the chain produces an optimal broadening of the transmon resonance. We consider an infinite chain. Since we are not interested here in modeling frequency dependent transport through the system, but only in the effect the chain has on the transmon, we replace the complicated frequency dependent impedance of the chain Z chain (ω) with its constant characteristic impedance R = L J /C g . We replace the 50 Ω (low impedance) transmission lines by ground connections. These assumptions produce the simple linear circuit depicted in Fig. S8.
Within this linear model, the resonance frequency ω T and decay rate Γ T of the transmon are obtained as respectively the real and imaginary parts of the relevant pole of the frequency dependent impedance between A and B in the circuit diagram. This impedance is given by . Toy model of the whole circuit. Simplified lumped element model used to qualitatively understand the link between the transmon decay rate and the impedance of the environment.
We have to note here that we have oversimplified the model, which now predicts an overdamped regime at small C sh .
In the real device, overdamping is prevented by the sizable capacitances C g,T and C g,T2 , which we have dropped. A quick fix, is to use a value for C sh that is comparable to C g,T and C g,T2 (several tens of fF), rather than its actual value of 4.4 fF. The behavior of the resonance frequency is easy to understand. At small R, one effectively has an LC circuit with capacitance C sh + C c /2 and resonance frequency ω T = 1/ L J,T (C sh + C c /2), while at large R, one has an isolated SQUID loop with resonance frequency ω T = 1/ L J,T C sh . In Fig. S9 we present the behavior of Γ T vs. R for C c = 119 fF (its actual value) and C sh = 80 fF, chosen to give rough quantitative agreement with the experimental results we present in the Fig. 4 of the main text, although the qualitative behavior does not change if we change C sh moderately. At small R, the behavior of the decay rate Γ T is (proportional to R) while at large R, it is given by (proportional to 1/R). A good estimate of the value R * that maximizes Γ T is obtained by equating the small and large R asymptotic expressions for Γ T . This yields where Z T,simp = L J,T /C sh is the characteristic impedance of the transmon, in the simplified circuit of Fig. S8. Thus the largest coupling (as measured by Γ T ) is obtained when the characteristic impedances of the transmon and chain match up to factors of order one, a result that is familiar in microwave engineering. At this optimal chain impedance, the decay rate Γ T is proportional to ω T with a proportionality constant that is a function of C c /C sh . This constant can reach values of order one, implying ultra-strong coupling is attainable. In our actual device, we find Γ T to be a decreasing function of L J , in the L J window to which we have access, suggesting that the lowest chain impedance that we can reach, is larger than the optimal value R * . When we compare the characteristic impedances Z T = /(2e) 2 2E C,T /E J,T 760 Ω (transmon) and Z C = L J /C g 1590 Ω (chain) of the actual device, we see that indeed Z C > Z T . Note that here we took a realistic estimate for the transmon impedance, that includes the effect of the capacitances C g,T and C g,T2 which were dropped in our qualitative model. Had we naively taken Z T = L J,T /C sh we would have obtained Z T = 1990 Ω, which though still close to the chain impedance, might have lead us to expect to observe a maximal value for Γ T as we sweep the L J window to which we have access.

F. Dealing with the transmon nonlinearity
The results in Fig. 3 in the main text confirm that the transmon qubit is a non-linear quantum circuit element that is strongly coupled to the chain. Here we review a standard way to deal with this anharmonicity. 7 The method is known as the self-consistent harmonic approximation (SCHA) because the anharmonic term is replaced by a harmonic one whose magnitude is determined self-consistently, via the variational principle. We also determine the regime of validity of the approximations we introduce. Let us consider the complete Hamiltonian of the device, neglecting only the weak non-linearity in the chain elements: Here we found it convenient to define an operator ϕ N +1 ≡ 0 which is not an extra degree of freedom, but simply the zero operator. To shorten notation we don't indicate the Φ T dependence of E J,T or the Φ C dependence of E J explicitly here. Were it not for the term −E J,T cos( ϕ T ), the quantum system described by the Hamiltonian in Eq. (S12) (Eq. 8 in the main text) would have been equivalent to a set of coupled harmonic oscillators, and therefore straightforward to solve. The term −E J,T cos( ϕ T ), not being quadratic in ϕ T , produces an interacting many-body problem. The strategy will be to replace the transmon terms in H with more tractable, yet accurate counterparts. For this purpose our starting point is to consider the Hamiltonian H in the limit where the inductances L J between chain nodes go to infinity (E J → 0), so that the charge on each chain island is conserved, and we can treat n j , j ∈ {1; . . . ; N } as ordinary numbers. Since the transmon then does not couple to any dynamical degrees of freedom, we refer it as isolated. The conserved chain charges contribute to the offset charge for n T . We will abuse notation slightly and still denote the transmon's charge degree of freedom, which now incorporates this additional offset, by n T . The isolated transmon Hamiltonian, in which reference to the conserved charges n α , α ∈ {1; . . . ; N } has been eliminated, reads Due to charge being quantized in units of 2e, the state space of H T is restricted to states |ψ for which where the offset charge n T (an ordinary number) contains contributions from the total transmon charge, the charge on each chain island, and from gate charges. We denote the eigenbasis of n T by |ν , i.e.
Owing to (S14), ν is quantized such that The matrix elements of H T in the |ν basis read The problem is equivalent to that of a charge e particle of mass E −1 C,T confined to a ring of circumference 2π that is threaded by a flux of n T times the flux quantum Φ 0 = h/2e. The phase observable ϕ T plays the role of the position coordinate, and the particle has an electrostatic potential energy E J,T (1 − cos ϕ T ). This system is easily solved numerically. In Figure S10 the low energy spectrum is plotted as a function of the offset charge n T , for three E C,T /E J,T ratios. States with energies sufficiently less than the height 2E J,T of the cosine well are insensitive to the offset charge n T . This is easy to understand in the equivalent picture of the particle confined to a ring: Sensitivity to the flux inside the ring requires the interference of paths with different winding numbers around the ring. However, at energies below 2E J,T , paths with non-zero winding number are exponentially suppressed by the tunneling amplitude to go through the cosine barrier. States with energy 2E J,T on the other hand are sensitive to the offset charge n T . As n T varies form 0 to 1/2 (half a Cooper pair), each of these energies sweep through an interval (or band) of width comparable to the spacing between levels. In the equivalent picture of a particle on a ring, this is a manifestation of the Aharonov Bohm effect. In a real experiment, the offset charge n T is subject to environmental noise. Performing spectroscopy on the levels sensitive to n T will therefore produce a noisy signal in which the extracted level energy "jumps around" inside the band through which the energy sweeps as n T is varied. When we reduce L J from infinity to its actual value, thus coupling the transmon to dynamical degrees of freedom in the chain, a numerically exact solution is no longer possible, given the large size of the Hilbert space. An obvious approximation scheme for states with energies below 2E J,T is the following. For these states, the phase observable ϕ T is unlikely to make excursions over the top of the cosine barrier (phase slips) at ϕ T = ±π. For such states it should therefore be permissible to replace the ring to which the particle is confined with the whole real line, and the cosine potential with a parabola. Note that this approximation ignores the restriction (S14), responsible for charge quantization. For states with energies E J,T , the phase is confined very close to the minimum at ϕ = 0 of the cosine potential, and the replacement E J,T (1 − cos ϕ 0 ) → E J,T ϕ 2 0 /2 is legitimate. This leads to a harmonic spectrum ω n = E C,T E J,T (n+1/2). However, the quadratic approximation can be improved to have a larger regime of validity, in the following way. We approximate the eigenstates of the isolated transmon as those of the parent Hamiltonian where the parameter E S is optimized according to some criterium in order to give the best possible agreement with the exact solution. Here we use the criterion that the energy E S should be chosen to minimize H T where the expectation value is taken with respect to the ground state of H P . For given E S , the eigenstates and energies of H P are |n = (B † ) n √ n! |0 , E (0) n = ω T (n + 1/2), n = 0, 1, 2, . . . , where Expressed in terms of the bosonic opertors B and B † , and manipulated into normal ordered form, the full transmon Hamiltonian reads The expectation value 0| H T |0 evaluates to The minimal value for 0| H T |0 is produced by E S satisfying the equation which can also be written as an equation determining the transmon frequency ω T . We note that in principle, the approximation can be further improved by treating H P with the optimized value (S23) for E S as the zero'th order approximation and treating H T − H P as a small perturbation. In general, the leading corrections in such a perturbation expansion are of first order in λ 2 . However, for the ground and first excited states, it is one order higher, i.e. λ 4 . Thus, the approximation H T H P is particularly accurate for the ground and first excited states of the transmon, which given the probe power and plasma frequency of the chain, are the ones we are interested in, in the experiment. The assumed smallness of λ allows us to solve approximately the self-consistency equation (S23) to get which gives an excitation energy In Figure S11 we compare the results (S24) and (S26) to the exact excitation energy of the isolated transmon, and find that in the regime where sensitivity to the offset charge n T is weak, i.e. E J,T /E C,T 1, the approximation (S26) is indistinguishable from the more sophisticated (S24). Now we turn to the case E J > 0 where the transmon is coupled to the dynamical degrees of freedom in the chain. We take the same approach as before. In principle, it is possible to generalize the previous calculation in the following way. When E J > 0, we may choose E S so that it minimizes 0| H |0 , where H is the full Hamiltonian in (Eq. 4 in the main text) and |0 is the ground state of the parent Hamiltonian that is obtained from H by making the replacement E J,T (1 − cos ϕ T ) → E S ϕ 2 T /2. We have implemented this approach, but find that it yields an insignificant improvement upon the simpler approach of simply taking E S = E J,T − E J,T E C,T /4, at least in the one photon sector, and for the parameters of the current device.

G. Analytical formula for the Scattering Phase Shift
In the section Methods -Analytical formula for the Scattering Phase Shift we presented a formula (Eq. 11) for the phase shift φ(ω, Φ T , Φ C ) of a mode with frequency ω of the full system at transmon flux Φ T and chain flux Φ C : Here we give the derivation. We work in the thermodynamic limit where N → ∞. This means that both matrices C and J (Eqs. (S32) and (S33)) become semi-infinite. Furthermore, the approximation discussed in the previous section eliminates the anharmonic term in H, and introduces four new non-zero matrix elements in J, namely Making the flux-dependence of E S explicit, we have from Eq. (S25) with E J,T (Φ T ) = E J,T,max |cos (πΦ T /Φ 0 )|. Explicitly the matrices C and J then read (see Eqs. (6) and (7) in the main text) We now have a fully linear system. The phase shifts in Eq. (S28) are obtained by solving the classical equations of motion. We start by defining a vector with the superconducting phases in each island π T = (ϕ L , ϕ R , ϕ 1 , ϕ 2 , · · · , ϕ N ).
The equations of motion for the mode at frequency ω are given by The solution to the part of Eq. (S34) involving degrees of freedom in the chain can be taken as Here a is the length of the unit cell of the array, N (ω) a frequency dependent amplitude and κ(ω, Φ C ) is the wave number of a wave that propagates in the chain with angular frequency ω. It can be obtained from the dispersion relation in Eq. (S2), where ω p (Φ C ) = 1/ L J (Φ C ) (C J + C g /4) is the plasma frequency of the chain. The phase shift φ(ω, Φ T , Φ C ) in Eq. (S35) is determined by the components of Eq. (S34) involving the transmon islands. They read Using Eq. (S37) and Eq. (S38) we eliminate ϕ L and solve for ϕ R in terms of ϕ 1 to obtain Substituting this into Eq. (S39) and using Eq. (S36) for κ(ω, Φ C ) we obtain with Using the mode definition in Eq. (S35) for ϕ 1 and ϕ 2 in Eq. (S41) leads to Finally, using again the expression in Eq. (S36) for κ(ω, Φ C ) we obtain the expression for the phase shift φ as a function of the system parameters.
In the main text we defined the relative frequency shift δφ n in terms of the discrete mode frequencies of the full system (transmon plus finite chain of N nodes) at respective transmon fluxes Φ T and Φ 0 /2. We repeat the definition here: Here we relate this to the relative scattering phase shift δφ(ω, Φ T , Φ C ) for the infinite system (see Eq. 14 in the Methods section of the main text) by showing that For conveniece we assume that island N + 1 is grounded. (The precise boundary condition becomes immaterial in the N → ∞ limit.) Had the chain been open to the left of node 1, the eigenmodes would have been ϕ j ∝ cos κ 0 n aj , j ∈ 1, 2, 3, . . . , N with wave numbers given by κ 0 n = (n − 1/2)π/N a with n ∈ 1, 2, 3, . . . , N . In the presence of the transmon to the left of chain node 1, the eigenmodes inside the chain change to cos (κ n aj − φ n ). Here φ n is the additional phase introduced by the transmon. Now the κ n depend on φ n too. Assuming the boundary conditions that the nodes to the left of transmon island L and to the right of chain island N are grounded, they are given by The modes of the system follow a dispersion relation ω n (Φ T , Φ C ) = ω(Φ C , κ n ). The notation must be understood as follows: ω n (x, y) and ω(x, y) denote distinct functions. The two arguments of the former refer to respectively the flux in a transmon and in a chain SQUID, and for given fluxes, the function assumes the value of the frequency of system mode n. The first argument of the latter function ω(x, y) refers to the flux in a chain SQUID, while the second argument refers to the unquantized wave number of a mode in the infinite chain. The function evaluates to the frequency corresponding to the given wave number (which does not depend on the transmon flux). For sufficiently large N we can expand the dispersion relation around κ 0 n a with corrections of order N −2 . The dependence on the transmon and chain fluxes is included. Similarly, we can expand ω n−1 (Φ T , Φ C ) around κ 0 n a to obtain Here we made use of the fact that φ n − φ n−1 = O N −1 . Therefore we can set φ n−1 = φ n introducing an error of O N −2 which can be ignored for large N . We can now obtain the terms in Eq. (S45), Substituting this into Eq. (S45), and noting that φ n (Φ T , Φ C ) = φ (ω n (Φ T , Φ C ), Φ T , Φ C ), we obtain Eq. (S46).

I. Relation between the Phase Shift and the impurity response function
To elaborate the link between the local impurity response function and the phase shift induced by the transmon qubit, we define three spectral densities corresponding to the phase-phase, phase-charge and charge-charge response of the transmon up to constant prefactors. In order to calculate these spectral densities we turn to the quantum mechanical problem (In the previous section, we could compute the phase shift by solving the classical equations of motion). With each mode ω we associate cannonical bosonic operators b ω and b † ω . These are related to the charge operator n j and phase operator ϕ j for each of the islands as where the profile ϕ j (ω) is normalized such that j=L,R,1,2,··· The normalization constant in Eq. (S35) is thus set to From the definition of n T and ϕ T (see the text below Eq. 7 in the Methods section of the main text) follows Explicitly, we find: In the Heisenberg picture b ω (t) = e iωt b ω and b † ω (t) = e −iωt b † ω . Using this to calculate the spectral densities A j (ω) for ω > 0 we obtain Now we compare the three correlation functions with the frequency derivative of δφ(ω, given in Eq. (S28). In Fig. S12 we plot the four curves. We see that the four curves overlap around the transmon frequency ω T . This means that the width and the center frequency obtained from the scattering phase shift are good estimations of the real width Γ T and frequency ω T of the transmon.
FIG. S12. Connection between phase shifts and qubit dissipation. Comparison between the three correlation functions and the energy derivative of the phase shift.

J. Breakdown of the rotating wave approximation
In this section, we investigate the applicability of the rotating wave approximation (RWA), a common technique for analysing the light-matter interaction at sufficiently weak coupling. The regime in which the light-matter coupling is so large that this approximation becomes inaccurate, is referred to as ultra-strong coupling. We will find that indeed, the RWA leads to errors of a few percent for our device.
To set up the RWA, we have to identify the harmonic oscillator basis that diagonalizes the (finite) chain part of the Hamiltonain. For this purpose, it is convenient to define N × N matrices C −1 chain and L −1 chain with entries i.e. C −1 chain and ϕ 2 0 L −1 chain are the lower right N × N blocks of respectively the inverse of the full (N + 2) × (N + 2) capacitance matrix C, and of the full (N + 2) × (N + 2) Josephson matrix J. We then define an N × N matrix Π chain and a positive definite diagonal matrix ω chain such that the columns of Π chain contain the eigenvectors of C −1 chain L −1 chain while the diagonal entries of ( ω chain ) 2 are the corresponding eigenvalues, i.e.
Since the eigenvalues are real, we can and will choose the entries of Π chain to be real as well. This definition determines each column of Π chain up to a normalization constant. We fix these constants by demanding that where [ω chain ] l is the l'th diagonal entry of ω chain . We also define a matrix Ξ chain = ϕ 2 0 L −1 chain Π chain ω −1 chain / .
It is easy to verify that the row l of Ξ T chain contains the left-eigenvector of C −1 chain L −1 chain that is associated with eigenvector ([ω chain ] l ) 2 . As a result Ξ T chain Π chain is guaranteed to be a diagonal matrix. Furthermore, due to the normalization condition (S68) we chose for Π chain , the diagonal entries of Ξ T chain Π chain are all equal to unity. Thus We now define N operators Owing to (S70) and the fact that [ n j , ϕ k ] = iδ j,k , the operators b chain,k , k = 1, 2, . . . , N are bosonic annihilation operators, i.e. [ b chain,j , b chain,k ] = 0 and [ b chain,j , b † chain,k ] = δ j,k . Furthermore, for the chain part of the Hamiltonian we obtain For the term in the Hamiltonian that couples the transmon to the chain, we find A standard way to proceed from here is to truncate the full Hilbert space of the transmon to the subspace spanned by two lowest energy eigenstates |0 T and |1 T of the isolated transmon Hamiltonian (S13). This leads to a Hamiltonian of the Jaynes-Cummings type, ubiquitous in Quantum Optics. At sufficiently weak coupling, the expectation is that this should be accurate for studying the situation where a near-resonant excitation from the chain induces a transition between the ground and first excited states of the transmon. The operator n T is replaced by where we've chosen the overall phases of |0 T and |1 T such that 0 T | n T |1 T is real. Alternatively, a more controlled way to proceed is to make the self-consistent harmonic approximation (SCHA) (see Sec. F), which we have shown to be well-justified, and to express n T in terms of the resulting bosonic transmon operators [see Eq. (S20)], i.e.
After the SCHA, the Hamiltonian becomes quadratic, and no further approximations are required. We will however still consider the effect of making the RWA on this quadratic Hamiltonian, in order to asses whether or not the assumptions underpinning the RWA are valid in our device. The RWA approximation now involves dropping the transmon-chain coupling terms in which the transmon (in the unperturbed basis) is excited while a boson is emitted into the chain, or the transmon is de-exited while a boson is absorbed from the chain. We adopt the standard nomenclature and refer to the dropped terms as "counterrotating" (based on their time-dependence in the Dirac picture). Depending on whether this approximation is made in conjunction with truncating the transmon Hilbert space or with the SCHA, we either obtain an RWA Hamiltonian H RWA,1 = (E 1,T − E 0,T ) |1 T 1 T | + or [ω chain ] n b † chain,n b chain,n + 1 2 For both Hamiltonians, the ground state is trivial: the transmon is in its unperturbed ground state, and there are no bosonic excitations in the chain. We measure energy relative to this ground state. Both Hamiltonians leave invariant the subspace spanned by states in which there are no bosons in the chain, while the transmon is in its unperturbed first excited state, or there is one boson in the chain while the transmon is in its unperturbed ground state. The excited states relevant for spectroscopy at low driving power are found by diagonalizing the RWA Hamiltonians in this subspace. � (���) FIG. S13. Deviations from the microscopic model under the RWA assumption. Curves show the analytical formula for the relative frequency shift (Eq. 14 in the main text). Open circles show the corresponding result calculated using HRWA,1 and closed circles the result calculated using HRWA,2. In all cases the system parameters were taken as in Table 1 of the main text. The flux ΦC through chain SQUIDS was held fixed at zero. Results for six different fluxes ΦT ∈ [0, 0.3Φ0] through the transmon SQUID are shown.
In Figure S13 we compare the relative frequency shift (Eq. 1 in the main text) predicted by H RWA,1 and H RWA,2 to the analytical SCHA formula (Eq. 14 in the main text) derived for an infinite chain. We have also computed SCHA results for the finite chain of 4700 islands, and found that they lie on top of the infinite chain curves. We omit them from the figure to avoid clutter. We note that H RWA,1 and H RWA,2 give very similar results. This is consistent with our claim that at low energies, the SCHA Hamiltonian from which H RWA,2 derives, is a good approximation to the full Hamiltonian from which H RWA,1 is derived. We ascribe the small difference between results for H RWA,1 and H RWA,2 to the truncation by hand in H RWA,1 of the transmon Hilbert space to two states. (No such by-hand truncation was required in H RWA,2 .) If we fit an arctan line shape (Eq. 18 in the main text) to the SCHA curves in the figure, we find that the transmon resonance occurs at a frequency within about 0.01 GHz from h −1 times the energy difference between the ground and first excited states of the isolated transmon. Using the same procedure on the relative frequency shift predicted by either H RWA,1 or H RWA,2 on the other hand, gives a resonance frequency that is ∼ 0.1 GHz higher than h −1 times the ground to first excitation energy of the isolated transmon. We conclude that the RWA approximation is not quantitatively accurate, producing an error of between 2% and 5% for the transmon