Revealing the finite-frequency response of a bosonic quantum impurity

Quantum impurities are ubiquitous in condensed matter physics and constitute the most stripped-down realization of many-body problems. While measuring their finite-frequency response could give access to key characteristics such as excitations spectra or dynamical properties, this goal has remained elusive despite over two decades of studies in nanoelectronic quantum dots. Conflicting experimental constraints of very strong coupling and large measurement bandwidths must be met simultaneously. We get around this problem using cQED tools, and build a precisely characterized quantum simulator of the boundary sine-Gordon model, a non-trivial bosonic impurity problem. We succeeded to fully map out the finite frequency linear response of this system. Its reactive part evidences a strong renormalisation of the nonlinearity at the boundary in agreement with non-perturbative calculations. Its dissipative part reveals a dramatic many-body broadening caused by multi-photon conversion. The experimental results are matched quantitatively to a resummed diagrammatic calculation based on a microscopically calibrated model. Furthermore, we push the device into a regime where diagrammatic calculations break down, which calls for more advanced theoretical tools to model many-body quantum circuits. We also critically examine the technological limitations of cQED platforms to reach universal scaling laws. This work opens exciting perspectives for the future such as quantifying quantum entanglement in the vicinity of a quantum critical point or accessing the dynamical properties of non-trivial many-body problems.


Introduction
The past years have seen many advances in the design of simulators for strongly interacting fermions and bosons, using cold atom lattices [1] and polariton fluids [2]. These platforms are indeed well suited to controllably explore the dynamics of bulk many-body problems [3,4], especially in presence of collective phenomena such as the superfluid to Mott insulator transition [5,6]. The many-body effects triggered by a discrete quantum system at the boundary of a quantum gas (also known as a quantum impurity) have also been investigated in a controlled fashion using various types of electronic quantum dots, leading to the observation of universal scaling laws [7,8], and even more exotic quantum phase transitions [9,10]. With some exceptions [11][12][13], most studies on quantum impurities focused on finite-voltage, zero-frequency DC transport measurements, without the possibility to unveil the finite frequency dynamics of the impurity. In addition, quantum dot systems are strictly limited to fermionic environments, and proposals to realize analogous bosonic impurities in cold atoms [14,15] have not succeeded so far. Nevertheless, bosonic impurity problems have triggered intense theoretical research over the years [16][17][18][19], initially motivated by fundamental aspects of quantum dissipation [20,21]. Bosonic impurities can also be used to describe defects in strongly interacting bulk fermionic systems [22][23][24], using collective degrees of freedom [25]. The scarcity of controlled experiments in the many-body regime of bosonic impurity systems is therefore still a major issue.
It is thus clearly an important technological goal to engineer and characterize truly bosonic quantum impurities, and the simplest advocated path involves the coupling between an ultra-small Josephson junction to a controlled electromagnetic environment. Such devices have been studied theoretically, either in the framework of the spin-boson model [26][27][28][29][30][31][32] in case of a capacitive coupling, or the boundary sine-Gordon (BSG) models [33,34] in case of galvanic coupling, which will be the topic of our study. Indeed, the galvanic coupling limits the sensitivity of the system to fluctuating charges, which is also a nuisance for quantum simulators. The development of circuit quantum electrodynamics (cQED) [35] provides an ideal testbed for the design and the precise measurement of bosonic impurity models, thanks to the recent advances in accessing the finite-frequency response of microwave photons in the quantum limit. While their non-equilibrium dynamical properties open fascinating research directions [36,37], equilibrium spectroscopic studies constitute an important milestone that is necessary to characterize those complex impurity systems, and have been made available only recently [38][39][40][41]. In fact, there are many predictions for the finite-frequency linear response of bosonic boundary models that still await some direct experimental measurement [25], let alone more subtle non-linear phenomena that were more recently considered in the cQED context [27,[29][30][31].
Our main focus here is to unambiguously characterize spectral signatures of nonperturbative effects in a prototypical bosonic quantum impurity problem described by the BSG Hamiltonian. For this purpose, we investigate non-linear effects that are controlled by a single Josephson junction at the edge of a high-impedance superconducting transmission line, combining state-of-the art cQED fabrication techniques and measurements, with fully microscopic many-body simulations. The use of high impedance meta-materials here is crucial to enhance the quantum fluctuations of the superconducting phase variable at the boundary, reaching regimes where linearized theories become invalid and striking many-body effects prevail. The experiment that we designed uses an array of 4250 Josephson junctions, so that our system is close to the thermodynamic limit [40]. As a result, a large number of electromagnetic modes can be resolved, which makes measurements based on phase shift spectroscopy [39] very accurate. In addition, we use here the flux tunability of a SQUID at the terminal junction (the bosonic impurity), which allows us to test for the first time some predictions for the renormalized scale of the BSG Hamiltonian. Our work complements several recent experimental and theoretical studies [41][42][43][44][45] that target non-linear effects in a bosonic impurity model at impedances equal to or larger than the superconducting resistance quantum R Q = h/(2e) 2 ≃ 6.5 kΩ.
Of particular relevance is [45] in which strong non-linear losses were reported. A clear explanation in terms of quantum phase slips was given in the regime α = Z/R Q > 1 whereas the loss in the α < 1 regime could not be quantitatively reproduced with a model based on the non-linearity of the confining potential of the impurity. We target here the exploration of quantum non-linear effects at somewhat lower dimensionless impedances α ≃ 0.3 and at E J /E C < 1, where quartic (Kerr type) and higher order processes at the terminal junction are the dominant source of non-linear effects. In addition, the reactive response of the SQUID is studied in parallel. This allows to measure and explain both sides of the same problem for the first time with this type of system. A bosonic impurity can couple single and multi-photon states of a multi-mode resonator. The resulting appearance of multi-photon resonances in linear response has been demonstrated experimentally [46]. An important physical outcome of our study is a demonstration that these processes can induce a significant many-body dissipation channel in its surrounding transmission line as was also observed in Refs. [40,45]. This effect can be dominant over other known loss mechanisms, either of extrinsic (e.g. loss inside the measurement line) or intrinsic origin (dielectric losses, or magnetic flux noise, see Appendix H). Our microscopic modelling is able to reproduce the measured many-body losses at high frequency in the regime where the Josephson energy of the terminal Junction is a small parameter that can be treated perturbatively. Ref. [45] studies a complementary regime where the E J /E C ratio of the terminal junction is larger than 1. At α ≳ 1, Ref. [45] matches measured many-body losses to theoretical predictions. Viewed together with our work, this demarcates the regime of α ≲ 1 and Josephson energy comparable to charging energy as the frontier for further theoretical work or quantum simulation. We also critically examine scaling predictions from universal models, and we provide a clear path for the future development of superconducting circuits in order to address universal transport signatures, a hallmark of strong correlations.
The manuscript is organized as follows. In Sec. 2, we present the boundary sine-Gordon (BSG) model and our design from a superconducting transmission line terminated by a flux-tunable SQUID. The microwave measurement setup is also introduced, together with the full microscopic model describing AC transport in the device. In Sec. 3, we present our main data and extract both reactive and dissipative responses from the finite frequency spectroscopy. In Sec. 4, we present various numerical calculations of these observables. A self-consistent theory, valid when the SQUID is threaded with magnetic fluxes close to zero, shows a drastic renormalization of the frequency at the bosonic boundary, allowing also to extract the unknown parameters of the device. In addition, a perturbative calculation, valid close to half flux quantum where the Josephson energy becomes a small parameter, is able to describe precisely the dissipative effects due to multi-photon conversion, which are shown to dominate the high frequency response of the junction. We conclude the paper of various perspectives that cQED techniques open for the simulation of strongly interacting bosonic phases of matter.
2 Tailoring the BSG simulator

Design principles
The BSG model describes the quantum dynamics of a resistively shunted Josephson junction, which can be referred to as the weak link, the impurity or the boundary. It has a Lagrangian: The weak link has a critical current 2eE J /ℏ and a shunting capacitance C J that accounts for charging effects when Cooper pairs tunnel through the junction. In an idealized description, φ 0 is viewed as the boundary value φ x=0 of a continuous field φ x , and the environment is described by a continuous relativistic quantum field theory: where c is the phase velocity in the environment. In the language of electronic circuits, this ideal environment is an infinite transmission line with an impedance Z shunting the weak link that is constant at all frequencies, while ℏ∂ t φ 0 /2e is the voltage across the weak link. In principle, the capacitance C J provides an ultraviolet regularization: at high frequencies it shorts the circuit. However, the theoretical analysis of the BSG model often assumes that the environment has a finite plasma frequency ω p , that provides a lower ultraviolet cutoff than E C /ℏ = (2e) 2 /ℏC J . In this universal regime, it is well-established that the BSG model hosts a quantum phase transition. When Z > R Q , environmentally induced zero-point motion delocalizes the phase φ 0 , making it impossible for a dissipationless current to flow through the weak link. When Z < R Q , on the other hand, a dissipationless current can flow. This corresponds to the "superconducting" regime where the phase φ 0 is localized in a minimum of the cosine Josephson potential. We will focus in this work on the "superconducting" phase where Z < R Q , without making a priori assumptions about the shunting capacitance C J , which will turn out to play an important role in the description of our experimental device. Let us first gain a qualitative understanding of the system by expanding the Josephson cosine potential to the second order, which corresponds to replacing the weak link with a harmonic LC oscillator of resonance angular frequency ω J and characteristic impedance Z J , still shunted by an environmental impedance Z env (ω). When the weak link is fully decoupled from its environment, ω J = 1/ √ L J C J and Z J = L J /C J , with L J = (ℏ/2e) 2 /E J . Current biasing the junction means adding a source term ℏI(t)φ 0 /2e to the Lagrangian, contributing to the voltage across the junction ℏ∂ t φ 0 /2e. According to the quantum fluctuation dissipation theorem then, at zero temperature, the phase fluctuations are given by: where is the impedance of the resistively shunted linearized weak link. In the ideal case where Z env (ω) = Z, this gives: When Z J ≪ Z, the weak link itself shorts the resistive shunt, so that the integrand in (5) develops a narrow resonance and φ 2 0 ≃ πZ J /R Q is very small, since we assumed Z < R Q . When Z = Z J /2, φ 2 0 = 4Z/R Q , while if Z J ≫ Z, the ohmic environment shorts the junction over a broad frequency range and φ 2 0 ≃ 4Z ln(Z J /Z)/R Q . These simple electrokinetic considerations dictate the design a non-trivial BSG simulator. Firstly, the observation of non-linear effects requires φ 2 0 ≳ 1, hence Z J ≳ Z. In addition, the Josephson term E J cos(φ 0 ) contains all even powers of φ 0 and thus allows elementary processes in which a single photon disintegrates into any odd number of photons at x = 0. The largest signal of this disintegration is achieved when there is strong hybridization between the weak link and the environment. This requires that Z J is not too much larger than Z, and that ω J is below the plasma frequency ω p of the environment. The ideal situation is therefore to design a device where Z J ≃ Z, with Z on the order of (but smaller) than R Q . A summary of the physical domains of BSG is already given in the right panel of Fig. 1, and will be discussed further in the text below.

The circuit
In our device, which is depicted in Fig. 1, the weak link is actually a SQUID consisting of two nearly identical small physical junctions on opposite sides of a ring. This results in a Figure 1: Left. Schematics of the measured circuit. The Josephson junction array, depicted in blue, is characterized by its lumped element inductance L, capacitance C and ground capacitance C g . The chain is terminated by a nonlinear SQUID, depicted in red and characterized by the flux-tunable Josephson energy E J (Φ) and capacitance C J . a. SEM picture of a small part of the full JJ chain, composed of 4250 sites in total. b. SEM picture of the galvanic coupling between the JJ chain and the nonlinear SQUID. The junction is grounded on its other side. Right. Parameter space of the device. The vertical axis represents the linear response probe frequency ω in units of the array's plasma frequency ω p . The horizontal axis represents the flux-tunable Josephson energy E J (Φ) of the nonlinear SQUID in units of the SQUID's charging energy. The inset shows the two main BSG mechanisms: the reactive linear response is characterized by a renormalized Josephson energy E ⋆ J , and further acquires a dissipative component R ⋆ (ω) due to photon disintegration at the boundary. Region 1 shows the low frequency universal limit of the BSG model, which is only a narrow domain of parameter space. The BSG mechanism produces many-body signatures across region 2, which are strongest in a swathe around the SQUID's resonance frequency shown as grey shaded. In region 3 of moderate phase fluctuations, we constructed a microscopic mean field theory that accurately predicts the renormalization of E ⋆ J . In the high frequency region 4, we developed a perturbative microscopic method that accurately predicts the dissipative response of the device.
Josephson energy [47]: that can be tuned by varying the magnetic flux Φ through the SQUID ring. Here Φ Q = h/2e is the flux quantum. In the above formula, d characterizes the small accidental asymmetry of the SQUID, which is in the 10 −2 to 10 −1 range for the SQUIDS we fabricate. To estimate E J , we constructed several isolated junctions using the same lithographic process as for our full device, and measured their room temperature resistance. Using the Ambegaokar-Baratoff law to extract the critical current, we find E J (0)/h ≃ 25 GHz. The SQUID is further characterized by its internal capacitance C J , predicted to be in the 10 1 fF range from geometrical estimates.
To achieve an environmental impedance close to R Q , we employ a long homogenous array of N Josephson junctions. The junctions are large enough that they behave like linear inductors with inductance L of order 1 nH. Such large junctions possess shunting capacitances C in the 10 2 fF range. The Josephson junctions separate N + 1 superconducting islands, labeled 0 to N , each with a capacitive coupling C g ∼ 10 −1 fF to a back gate. At sufficiently low frequencies, the inductance L shorts the shunting capacitance C, and an infinite array of this type behaves like a transmission line with constant impedance Z = L/C g in the targeted kΩ range. However, above a plasma frequency ω p of the order of the resonance frequency 1/ √ LC of the capacitively shunted inductors, which is in the 10 1 GHz range, electromagnetic excitations are unable to propagate down the chain. We employ an array with N = 4250 junctions. As explained in Appendix A, we were able to perform a sample characterization that yielded the following values for the array parameters: L = 0.54 nH, C = 144 fF, and C g = 0.15 fF, so that ω p = 2π × 18 GHz, and Z/R Q = 0.3.
As seen in Fig. 1, the weak link connects to the rightmost node of the array. In order to perform spectroscopic measurements, the leftmost node is galvanically coupled to a Z tl =50 Ω micro-strip feed-line, in a T-junction geometry. In Appendix B we provide an explicit expression for the environmental impedance Z env of the array coupled to the feedline. Owing to the mismatch between the characteristic array impedance L/C g =1.9 kΩ and that of the 50 Ω transmission lines, Z env has sharp peaks at frequencies corresponding to Fabry-Perot resonances in the array. More than 100 modes can be clearly observed below the plasma cutoff in our device, most of which lie in the ohmic regime, see the top panel of Fig. 2 showing a close-up on five of those modes. These modes are well-resolved, with a maximum free spectral range (level spacing) ∆f FSR ≃ 0.4 GHz that is larger than the line width of each mode. This provides a means to study the system response by spectroscopic analysis, as we detail now.

Measurement protocol
With the feed-line connected to the high impedance array as described in the previous section, we have the two port device with a hanging resonator geometry that is depicted in Fig. 1. We perform spectroscopic measurements on the hanging resonator by sending an AC microwave signal to the input of the feed-line and collecting the signal arriving at its output. The ratio between the complex amplitudes of these two signals defines the transmission of the circuit S 21 . The microwave signal is calibrated to be sufficiently weak that no more than one incident photon on average populates the resonant modes of the circuit. Therefore, the circuit is close to equilibrium, so that standard linear response theory can be used for modelling its transport features. The array terminated by the weak link acts as a side-coupled resonator with a joint impedance to ground Z a+w . Solving the classical scattering problem for electromagnetic waves in the feed-line in the presence of this resonator, yields that S 21 = 1/(1 + Z tl /2Z a+w ), which we can rewrite as Here Z N = (2/Z tl + 1/Z a+w ) −1 is the linear response impedance between superconducting island N of the full device, where the array is connected to the feed-line, and the back gate (or ground). Note that (7) does not involve any semi-classical approximation. All quantum effects are included and Z N (ω) is the linear response impedance obtained from the Kubo formula [27]. We discuss this further in Appendix B, where we provide an explicit formula relating S 21 to the system parameters and the self-energy associated with the weak link.  Fig. 2 for two different values of the magnetic field. By tracking the resonance frequencies as function of flux through the SQUID, we obtain the lower panel of Fig. 2. The flux range here is limited to close to half a flux quantum due to the periodic dependence in magnetic field from Eq. (6). We note that the maximum amplitude of variation of each mode frequency is of the order of the free spectral range. This is due to the change of boundary condition, from nearly closed circuit (large E J ) at low flux to nearly open circuit (small E J ) at half flux. The only component of the device whose electronic properties has such a strong magnetic field dependence is the SQUID at the end of the array. We therefore conclude that our spectroscopic measurements are sensitive to the boundary term in our BSG simulator. Besides this frequency shift of the modes, we further observe a striking flux dependence of the resonance widths. The frequency at which the largest broadening is observed decreases in a manner similar to the expected flux dependence of the SQUID resonance frequency, suggesting that the greatest broadening occurs for modes that are on resonance with the SQUID. (Not shown in Fig. 2, but see Fig. 3 below.) Furthermore, the average broadening of resonances is larger when the SQUID Josephson energy is lowest. These features are not reproduced when the SQUID is modelled as a linear circuit element as in (16), suggesting that the measured transmission contains significant information about inelastic photon processes due to the boundary SQUID, in agreement with what has been observed in [45]. In the rest of the Article, we quantify these inelastic contributions to the measured signal and compare to theoretical predictions.

Shifting of the resonances
The dependence of resonance positions on the external magnetic flux reveals quantitative informations about the reactive effect of the weak link, which we can extract as follows [48]. The resonances we observe are predominantly single-photon in nature, and occur at photon wave vectors k quantized such that (N + 1/2)k l (Φ) + θ l (Φ) = π(l + 1 2 ) where Φ is the external flux, and θ l (Φ) is the phase shift associated with resonance l, which vanishes when the weak link is replaced by an infinite impedance. Provided that the asymmetry of the SQUID is small, we can take its effective inductance as infinite at Φ = Φ Q /2. The relative phase shift then measures the phase shift induced by the Josephson potential of the weak link. If we denote the l th resonance frequency as f l (Φ) and notice that θ l+1 − θ l = O(N −1 ) then it follows that where ∆f l (Φ) and ∆f FSR,l are respectively the frequency shift at Φ with respect to Φ Q /2 and the free spectral range for the mode l (see Fig. 2). As was done in other recent recent experiments [39,40] on dynamic quantum impurities in the context of superconducting circuits, we extracted the phase shift δθ l as a function of external magnetic field and frequency.
To explore the relation between the phase shift and the properties of the weak link, we pose the following question: What would be the phase shift if the non-linear weak link was replaced by an effective linear inductor of inductance L ⋆ J = (ℏ/2e) 2 /E ⋆ J , with E ⋆ J a renormalized Josephson energy ? A formula for the phase shift in terms of E ⋆ J , C J and the array parameters can be derived by finding the wave vectors k where the impedance between node N and ground, due to the array terminated in the weak link, vanishes. (See Appendix B for details.) In Fig. 3 we show the experimentally extracted relative phase shifts δθ l as a function of mode frequency f l , for various fluxes Φ, together with the best fit to the effective linear theory (solid lines). We find excellent agreement between the experimentally extracted phase shifts and the theoretically predicted curve. We use the derived formula (31) in Eq. (8), and fit the extracted phase shifts δθ l at different fluxes Φ, the fitting parameters being C J and L ⋆ J (Φ). (We use the array parameters quoted in Sec. 2.2.) From this, we extract C J = 14.5 ± 0.2 fF as well as the effective weak link as a function of flux. Its discussion is postponed to Sec. 4.

Broadening of resonances
The flux dependence of the broadening of spectroscopic resonances contains quantitative information about dissipation caused by photon disintegration in the weak link. We extract it as follows. Close to a resonance, the array together with the weak link has an impedance Z a+w,l ≃ Z dis,l − iZ reac,l × (ω − ω l )/ω l . Here Z dis,l is real and caused by dissipation internal to the array or weak link. It does not have a significant frequency dependence on the scale of the resonance width. Similarly −iZ reac,l × (ω − ω l )/ω l is the reactive (imaginary) part that vanishes at the resonance frequency ω l and has been expanded in frequency around the resonance. From this follows that resonances in our setup then have the familiar 'hanging resonator' line shape, where Q i,l = Z reac,l /2Z dis,l and Q e,l = Z reac,l /Z tl are internal and external quality factors characterizing respectively dissipation in the weak link plus array, and in the feed-line. As a result, 1 − |S 21 | 2 , as a function of frequency, has a Lorentzian line shape with halfwidth (1/Q i,l + 1/Q e,l )f l in frequency.
To study the internal dissipation due to the weak link, we therefore fit the hanging resonator line shape to individual resonances, and extract the internal broadening γ int,l = f l /Q i,l as a function of the resonance frequency and external magnetic field. In practice, connecting the feed-line to the array adds a small reactive part to the feed-line impedance. This causes a small peak asymmetry, which we include as another fitting parameter. Full details are provided in Appendix B. Besides the nonlinear processes taking place in the weak link, more mundane processes in the array can also contribute to γ int . In recent years, several groups have been investigating the mechanisms that may be responsible for internal losses in superconducting resonators. For most materials [49], including resonators made out of Josephson junctions [50], the main mechanism that induces internal damping in the single-photon regime is the coupling with a bath formed by two-level-systems in dielectrics nearby the resonator. This effect, discussed further in Appendix G, does not depend on the external flux Φ, and can thus be calibrated at Φ = 0, where the non-linear contributions to the internal losses nearly vanishes, since the SQUID phaseφ 0 variable is well localized by the strong Josephson potential. It produces a constant-in-flux contribution γ diel that we subtract from the total internal broadening.
In the results that we present below, we plot the resulting γ J = γ int −γ diel which represents the contribution of the broadening that is due to nonlinear effects due to the weak link only.
We analyze the non-linear damping γ J of the modes due to the boundary junction as follows. In Fig. 4, we plot resonance frequencies between 2 GHz and 9 GHz as a function of flux Φ. Each vertical column of data points in Fig. 4 represents a mode of the chain obtained from the frequency traces shown in Fig. 2. The color of each data point shows in log scale the broadening γ J due to the boundary junction, that we extracted by fitting the resonance line shape to (34), with the estimated dielectric losses in the chain subtracted. A dashed line indicates the effective weak link resonance frequency ω ⋆ as extracted from the phase shift data. We observe internal broadening varying from essentially zero (especially at fluxes close to zero where the system is nearly linear) to values exceeding 100 MHz, with excellent correlation between the effective weak link resonance frequency ω ⋆ J and the maximum internal broadening. As the frequency cuts shown in the right panel of Fig. 4 reveal, γ J (ω) decays exponentially for ω > ω ⋆ J . The individual data sets with Φ/Φ 0 > 0.45 each show γ J decreasing by two decades as the probe frequency is scanned.
We have further estimated contributions to the internal broadening in the weak link itself due to other mechanisms not directly related to BSG physics. Obvious candidates are coupling with normal quasiparticles or dielectric loss [45,51,52] (cf. Appendix H) or inhomogeneous broadening from fluctuations in magnetic flux Φ through the SQUID. The latter, discussed in Appendix I, indeed depends on Φ, but is sufficiently small to be discarded in our setup. It would furthermore produce a Gaussian line shape which is not what we observe. Normal quasi-particle tunnelling or dielectric losses in the SQUID both peak at frequencies close to the weak link resonance frequency. However, even under unrealistically favorable assumptions for these processes, they can contribute at most between 10 0 and 10 1 MHz to broadening (cf. Appendix H). We therefore conclude that the results in Fig. 4 are a clear manifestation of BSG physics in the weak link. The magnitude of the damping can be used to calculate the round-trip decay probability in the circuit for the single-photon excitations via: and is equal to 0.25 for the maximal measured damping, a hallmark of ultra-strong coupling showing the large dissipation induced by the nonlinearity. In Ref. [45], similar round-trip decay probabilities are obtained in the transmon regime (1 < E J /E C < 5) and α ≃ 2, while smaller decay probabilities were obtained in the transmon regime at α ≃ 0.7. We now vindicate these qualitative effects by a microscopic modeling of the device.

Theoretical modelling of the observed many-body physics
Modelling theoretically the many-body effects in our experiment is challenging. An exact treatment is not feasible, due to the huge Hilbert space associated with the large number (a few hundreds) of modes that are involved in the ohmic range of the spectrum. However, we can take advantage of the tunability of the SQUID junction to investigate in a controlled way the regimes of large and small Josephson energies E J (Φ) of the boundary junction. We therefore discuss these two regimes separately.

Reactive effects at large E J
At given α < 1, fluctuations of the boundary phase φ 0 are controlled by the ratio E J /E C .
Here we focus on the regime where E J /E C is sufficiently large that phase fluctuations do not much exceed unity. The weak link has a charging energy E C of around h × 10 GHz. At zero external flux, the weak link Josephson energy is more than twice as large, and the approximation in which the boundary Josephson energy is replaced by that of a linear inductor is adequate to capture the reactive aspects of the dynamics. Moving away from zero flux, a better approximation is to replace the bare value E J by a renormalized one E ⋆ J , also called the self-consistent harmonic approximation (SCHA) [25,53,54]. This mean field theory is known to remain accurate at moderate phase fluctuations, when the phase explores more than the very bottom of the cosine Josephson potential, but does not tunnel out of the potential well. This regime corresponds to the region 3 of Fig. 1. To implement the SCHA for our circuit, we write the Josephson potential aŝ The term in parenthesis is viewed as a perturbation that will be dropped, and E ⋆ J is chosen to make the resulting error as small as possible. This leads to the self-consistency conditions that the expectation value of the perturbation with respect to the ground state of the effective linear system should vanish. This self-consistency condition can be rewritten as where phase fluctuations φ 2 0 are computed using Eqs. (3)(4) and the environmental impedance (28) . For given bare Josephson energy E J (0) at zero flux and SQUID asymmetry d, the self-consistency condition (13) allows us to generate a curve E ⋆ J (Φ), which is expected to be accurate at flux Φ not too close to Φ Q /2. We treat E J (0) and d as free parameters and adjust this curve to the E ⋆ J (Φ) data that we experimentally extracted with the aid of phase shift spectroscopy, see Sec. 3.2. We use data for 0.35 < Φ < 0.47Φ Q . This provides estimated values for the zero-flux bare E J (0)/h = 27.5 GHz and SQUID asymmetry d = 2%. The estimated asymmetry is reasonable for our fabrication process for small junctions given that we aimed for a perfectly symmetric SQUID. The value of E J (0) is in good agreement with the estimate of E J (0)/h = 25±1 GHz obtained from the measurement of the room temperature resistance of isolated test junctions fabricated on the same wafer and at the same time as the full device. This confirms the accuracy of the SCHA at relatively large E J (Φ).
The renormalization of E J is a textbook feature of the BSG model. Deep in the overdamped limit Z J ≫ Z, Eq. (5) yields φ 2 0 ≃ 2α ln E C /(2πα) 2 E ⋆ J , with α = Z/R Q the dimensionless resistance of the perfectly Ohmic environment. This is known as the scaling regime. Solving the self-consistency condition (13) then yields the well known scaling law: showing a strong downward renormalization of the Josephson energy E ⋆ J ≪ E J in the non-perturbative regime 0.1 < α < 1. Note that E ⋆ J cannot exceed the bare value E J , which is why it has been bounded in Eq. 14. Note also that this scaling law predicts a superconducting to insulating Schmid transition at the critical value α = 1, where the Josephson energy E ⋆ J renormalizes to zero due to a divergence of the phase fluctuations φ 2 0 . At frequencies sufficiently below the plasma frequency, the weak link in our device sees an effective environmental impedance Z c ≃ L/C g = 1.9 kΩ so that α = 0.3. It is interesting to ask how the renormalization of E ⋆ J that we observe in our device compares to the renormalization predicted for an idealized system in the scaling regime.
In Fig. 5, we plot the observed renormalization E ⋆ J (Φ) of the boundary Josephson energy in our BSG device, as a function of the bare scale E J (Φ). The solid line shows the result that the fully microscopic SCHA calculation predicts for our device, and the dashed line shows the scaling law (14) for an idealized system with the same E C and α as in our device, but with infinite plasma frequency ω p → ∞. We observe that E ⋆ J (Φ) starts with a weak renormalization E ⋆ J /E J ≃ 0.9 at low flux (large bare E J ). Close to half flux quantum (small bare E J ), the measured renormalized scale becomes as small as E ⋆ J /E J ≃ 0.2 ≪ 1. Except for this low flux regime, where it becomes invalid, the full SCHA provides an excellent description of the data. This underscores the fact that a detailed characterization of the environment is necessary in order to achieve agreement between theory and experiment in cQED simulators [39,55]. By nature of the universal regime E J ≪ E C , the analytical scaling law (14) should meet the full SCHA result at small E J , which is what is seen indeed for E J (Φ) < 1 GHz. However, this is already the domain where the phase fluctuates very strongly, and both the SCHA and the scaling law are inapplicable.
Indeed, from the observed E ⋆ J (Φ), we can estimate phase fluctuations using the selfconsistency condition (13) as ⟨φ 2 0 ⟩ = 2 ln(E J /E ⋆ J ). We plot the estimated ⟨φ 2 0 ⟩ as a function of the renormalized impedance Z ⋆

Dissipative effects at small E J
For α < 1, the BSG model is known to flow to a Kondo-like strong coupling fixed point in the limit of zero temperature and for frequencies below a small emergent scale that characterizes the low-frequency inductive response of the weak link. The response of the system at these low frequencies are beyond the reach of a perturbative treatment [25].
Here we denote that scale E ⋆ J because in the universal regime Z J ≫ Z, it has the same scaling as in Eq. 14. Note however that the device we are modelling is not in the universal regime. Nonetheless we may expect E ⋆ J < E J . In order to tackle the strong non-linear regime of small Josephson energy, we develop perturbative theory that is controlled in the high frequency domain ℏω ≫ E ⋆ J , using E J as a small parameter (compared to E C and ℏω p ), which corresponds to the region 4 of Fig. 1. For simplicity, we outline here the zero-temperature calculation based on time-ordered Green's functions. Experiments on our device are performed at a temperature T ≃ 30 mK≃ 0.6 GHz (in units of h/k B ), that is of the same order as E J (Φ Q /2) and we therefore have to include finite temperature in our numerical calculations. The generalization to finite temperature is discussed in Appendix D.
If we set E J to zero, the impedance between the node zero of the array and ground is An explicit expression (28) for the impedance Z env that shunts the weak link is provided in Appendix B. At zero temperature, this impedance is related to the time-ordered Green's function (The superscript 0 of the Green's function indicates that it is calculated at E J = 0.) As discussed in Appendix B, the weak link self-energy, given by Dyson equation, Σ(ω) = fully captures the effect of the nonlinear Josephson potential on the linear response of the system at zero temperature. Indeed, R Q ℏω/2πiΣ(ω) enters the linear response functions we eventually wish to calculate as the impedance of a circuit element connecting node 0 of the array to ground.
At first sight, it seems that a straightforward expansion in powers of E J (Φ) would allow us to calculate Σ(ω) for Φ in the vicinity of Φ Q where E J (Φ) is small. Indeed, dissipative effects show up at second order in E J . However, because the system's response is nonperturbative at frequencies below an emergent scale E ⋆ J /h, a regularization procedure is required in order to extract the response at frequencies in the measurement window between 2.5 and 11 GHz (well above E ⋆ J /h) perturbatively. We first discuss the formal perturbative expansion of the self-energy and subsequently the regularization procedure.
To second order in E J , we find the self-energy: The vertex Josephson energy, is given by In Appendix C we present two independent derivations of this result.
In principle, a strict perturbative calculation would use the bare Green's function Eq. (16). Quite generally, formula (17) implies that the self-energy introduces linear response resonances associated with a single incoming photon disintegrating into multiple photons at the weak link. At zero temperature, these multi-photon resonances occur at frequencies that are sums of single-photon resonance frequencies. Owing to the approximately linear dispersion relation of the array ω = vk in the ohmic regime, single-photon resonance frequencies are almost equally spaced. As a result there is a large near-degeneracy in these multi-photon resonances. For instance, if we denote the lowest bare resonance frequency by ω 1 , then there are 16 multi-photon resonances, each involving an odd number of photons, at frequency 10 ω 1 , which corresponds to a single photon resonance in the middle of the experimentally accessible frequency window. This leads to a highly singular behaviour of the self-energy in the vicinity of these degenerate clusters of multi-photon resonances, when it is built on bare Green's functions. In our device this is not mitigated appreciably by the slight curvature of the photon dispersion or by geometric irregularity [46]. This singular behaviour is however spurious as it does not take into account the significant many-body level repulsion between multi-photon states coupled directly or indirectly by the highly non-linear terminal junction. We therefore self-consistently dressed all propagators with self-energy insertions to obtain what is also called a skeleton diagram expansion, or self-consistent Born approximation, which introduces many-body level repulsion and smoothens the self-energy. This is why we used the full interacting Green's function G φ 0 φ 0 (ω) in Eq. (17), which is determined self-consistently together with Σ(ω) from Dyson equation: with G 0 φ 0 φ 0 (ω) given by Eq. (16). Let us now discuss the regularization procedure. When naively expanding in E J around zero, the Debye-Waller factor E v J /E J = exp − φ 2 0 /2 is zero, due to a logarithmic divergence in φ 2 0 . As a result, an unphysical answer Σ = 0 is obtained, so we do need to regularize the self-energy at low frequencies by introducing a counter-term E cutoff , where E cutoff must be larger than the true renormalized scale E ⋆ J . The intuitive picture behind this regularization procedure is as follows. We imagine adding an extra linear inductor L cutoff = (ℏ/2e) 2 /E cutoff in parallel to the weak link to our model. It only adds a parabolic potential E cutoff φ 2 0 /2 that remains flat for φ 0 = O(1) to the Hamiltonian. At very low frequencies, this inductor shorts the weak link, thus providing an infrared regularization, but at frequencies in the measurement window, it hardly carries any current, and thus should not affect results. We have taken E cutoff = 0.05 E J (Φ Q /2) ≈ 2πℏ × 0.02 GHz, and have checked that other choices of the same order of magnitude give the same results in the high frequency regime ℏω ≫ E ⋆ J where the calculation is controlled. Using the expansion presented above, and considering again the bare Josephson energy E J and SQUID asymmetry d as free parameters, we compare in Fig. 6 the measured internal linewidth (dots) to the theoretical predictions (lines), for three values of the flux Φ/Φ q = 0.48, 0.49, 0.5. Note that two of those three curves are sufficient to fully determine E J and d, so that the theoretical curve at Φ/Φ q = 0.48 contains no fitting parameter. The estimated values of the fitting parameters are reported in Tab. 1 for the various measurements that have been performed (including the room temperature critical current and the fit of the renormalized scale E ⋆ J ), which give all very consistent results. For smaller flux values, corresponding to the region 2 of Fig. 1, the junction frequency ω J enters the measurement windows, and our theory surprisingly still describes qualitatively the maximum observed in the loss function γ J (ω k ) for ω l ≃ ω ⋆ J . However, the magnitude of γ J is largely underestimated in the calculation, see Appendix F. This dis-crepancy is due to a breakdown of the expansion in powers of E J , as we confirmed by computing all the order E 3 J perturbative terms. At Φ < 0.4Φ Q we find that the E 3 J contributions are of the same order as the E 2 J contributions, while they remain negligible for the larger flux values of Fig. 6. At smaller fluxes, the superconducting phase is trapped near minima of the periodic Josephson potential, and non-perturbative 2π-phase slip processes between minima provide the dominant contribution to the damping process, which are not taken into account in our perturbative treatment. Deviations from our model thus gives an estimate of these phase slip processes at α ≲ 1. These have also be investigated theoretically and experimentally at α ≳ 1 [43][44][45].
Finally, we stress from figure 6 that the universal scaling law (See Appendix E) controlling the junction damping, γ J (ω) ∼ ω 2α−2 , is not obeyed in our measurement, rather an exponential decay is observed instead. This is expected because the scaling laws of the BSG model should be manifest on dynamical quantities only if E ⋆ J ≪ ℏω ≪ E C , corresponding to the region 1 of Fig. 1. However, the charging energy E C of the junction is too small to fullfill both constraints together. The observed exponential decay can be explained qualitatively from the influence of the high energy cutoff on the photon conversion processes. When increasing the probe frequency, the number of available photonic states at higher frequency decreases exponentially (due to the reduction in combinatorics), drastically reducing the possibility of recombination of a single photon into various multi-photon states. One solution to observe the power law mentioned above would be to optimize the design of the boundary junction in order to increase E C , but also to push the plasma frequency ω p to higher values by increasing the transparency of the Josephson junctions in the chain, or by replacing them by a disordered superconductor. This would require important technological advances in the field of cQED.

Conclusion
In this work, we demonstrate a two-fold interplay between the boundary and the bulk from finite frequency measurements. On one hand, the bosonic environment induces a reactive response on the boundary degree of freedom, strongly renormalizing its resonance frequency. This effect is captured in the regime of large Josephson energy compared to its charging energy at the boundary junction, so that fluctuations of the superconducting phase variable remain moderate, and an effective linear model can apply (region 3 in Fig. 1). On the other hand, the boundary is also able to induce a dramatic dissipative response onto its environment, due to efficient frequency conversion into multi-photon states, which were shown to dominate over known sources of photonic losses, in accordance with what was reported in [45]. When the Josephson energy is small enough, it can be used as an expansion parameter, leading to a perturbative theory which accounts well for the measured high frequency response (region 4 in Fig. 1). Both approaches led to consistent estimates of the unknown parameters at the boundary junction. To compare our measurements using quantum many-body theory, we developed a fully microscopic model of the circuit. We found excellent agreement in regimes where the non-linear effects could be controlled. Interestingly, these two extreme regimes border a large domain of parameters where non-perturbative phenomena fully develop, and our circuit challenges all theoretical approaches we are aware of (region 3 in Fig. 1). We also evidenced that the use of universal scaling laws have to be taken with a grain of salt in superconducting circuits, due to the limited measurement bandwidth and relatively low ultraviolet cutoff set by the junction charging energy (a few GHz) and the plasma frequency (about 18 GHz). This scaling regime, where power laws in various response functions should develop, corresponds indeed to a parameter space that cannot be easily explored (region 1 in Fig. 1).
Future experimental developments of bosonic impurities in cQED could lead to the observation of more dramatic many-body phenomena, such as quantum criticality [16], for instance the Schmid or spin-boson transitions that are predicted to occur at larger dissipation. Nevertheless, our work demonstrates that precursor effects of quantum phase transitions are worth investigating, because they exacerbate many-body behavior. The direct detection of the down-converted photons in cQED remains also a topic of interest, not only from the point of view of many-body physics [27,31,46], but also because they could be used as a potential quantum information resource.

A Determining the array parameters
In order to find the chain parameters C, C g and L, we measure its dispersion relation using standard two tone spectroscopy, which allows us to accurately measure resonance positions from below 1GHz up to the plasma frequency at 19GHz The results as a function of wave number are shown in Fig. 7. For k ≪ 1 the dispersion relation (21) is linear: ω k = k/ LC g . For k > C g /C on the other hand, the dispersion relation saturates to the plasma frequency ω p = 1/ L(C + C g /4). These asymptotic behaviors allow us to fit C g and L once C is known. We determine C from knowledge of the area of the junctions composing the chain: C = 45fF × area[µm −2 ] to set C = 144 fF. From the fitting to the measured dispersion relation: we then estimate L = 0.54 nH and C g = 0.15 fF. Hence, we find the array impedance Z c = 1.9 kΩ and the plasma frequency ω p = 18 GHz.

B Green's functions and impedance
The phase-phase correlation function and the linear response impedance play a crucial role in the analysis performed in this work. Here we elucidate their connection, and work out the various ingredients that are relevant for modelling our device. At zero temperature, it suffices to study time-ordered Green's functions, which is what we will discuss here for the sake of simplicity. At finite temperature, we employ (equilibrium) Keldysh Green's functions in our numerical computations. The necessary generalizations are discussed in Appendix D. Associated with the phase variables φ n on each island n = 0, 1, . . . , N of the array, we define the Green's function has the property G φm,φn (−ω) = G φn,φm (ω). At positive frequencies, G φm,φn (ω) = G R φm,φn (ω), i.e. the time-ordered Green's functions, convenient for diagramatic expansions, are equivalent to retarded Green's functions that describe the system linear response. The identification between retarded Green's functions and impedance embodied in Eq. (16) in the main text extends to all superconducting islands in the array.
We view G φm,φn (ω) as an (N +1)×(N +1) matrix G(ω). The corresponding impedance matrix describes a N + 1-port system obtained by associating a port with each superconducting island in the array, with one node of the port connected to the island, and the other to the back gate. See Figure 8. The operator corresponding to a current bias at port n is −ℏI(t)φ n /2e while the voltage across port m is ℏ∂ tφm /(2e). Hence, at positive frequencies, iωR Q G(ω)/2π is the impedance matrix of the N + 1-port system. The Dyson equation for G reads [G 0 (ω)] −1 − Σ(ω)/ℏ G(ω) = 1. Here G 0 (ω) is the matrix Green's function when the weak link Josephson energy E J is set to zero. The self-energy Σ(ω) incorporates the effect of the weak link. Since its energy E J (1 − cos φ 0 ) only involves the phase on island n = 0, Σ(ω) m,n = δ m,0 δ n,0 Σ(ω).
We thus identify at ω > 0 as the N + 1-port circuit admittance matrix in the presence of the weak link, while [iωR Q G 0 (ω)/2π] −1 is the same, in the absence of the weak link. The fact that Σ(ω) contributes additively to the admittance and has the form δ m,0 δ n,0 Σ(ω) then implies that as far as linear response is concerned, the effect of the weak link cosine potential is exactly equivalent to that of connecting a circuit element with impedance R q ℏω/2πiΣ(ω) across the nodes of port n = 0. Given the important role that impedance plays in our modelling, we now microscopically characterize the impedance of our device. We start by considering the array on its own, which we can view as an impedance network with a port at either end. One node of either port is connected to respectively the first or last superconducting island of the array and the other node of either port is connected to the back gate. The array 2 × 2 Figure 9: Diagrams that show the circuits as well as the positioning of current sources and voltage probes used to define the impedances Z a , Z ab , Z env and Z a+w . In each case with k given by Eq. (21). See Fig. 9 for the definition of Z a and Z ab . For completeness we mention that when N → ∞, the resulting single-port element has impedance which indeed reduces to L/C g when ω ≪ ω p . The total environmental impedance that shunts the weak link is then that of the finite array connected to the feed-line at the far end (see Fig. 9), Another impedance of special significance is Z N (ω), the impedance between array island N and the back gate, which through Eq. (7) determines the measured transmission in the feed-line. Since at island N , the array is shunted by the transmission lines that carry the input and output signal, where Z a+w is the impedance due to the array terminated in the weak link. In analogy to (28), it is given by which defines Z w , the impedance of the weak link (see Fig. 9). If we model the weak link Josephson term as an effective linear inductor L ⋆ J (Φ) = (ℏ/2e) 2 /E ⋆ J (Φ), as within the self-consistent harmonic approximation (SCHA), then The relationship between the effective inductance L ⋆ J (Φ) and the phase shift θ is obtained by setting k = [π(l + 1 2 ) − θ(Φ)]/(N + 1/2) solving Z a+w = 0 for θ. This yields Hence, we have an analytical expression for the relative phase shift defined in Eq. (8) where we approximate L ⋆ J to be infinite for Φ = Φ Q /2. If, on the other hand, we wish to include dissipation, we need to calculate the selfenergy. The expansion is performed around the limit L ⋆ J → ∞. In this case the weak link linear response impedance can be related to the self-energy using Eqs. (25), (26): Let us now investigate the line-shape of S 21 = 2Z N /Z tl in more detail, including the small reactive contribution to the feed-line impedance, which we have ignored up to now. In the vicinity of a resonance, one can write Here (Z tl − iX)/2 is the impedance due to the 50 Ω transmission lines carrying the input and output signals. Ideally this impedance would be purely real and equal to Z tl /2 = 25 Ω. In practice, the contact between the chain and the feed-lines contributes a small impedance −iX/2 in series, which is typically inductive (X > 0) and has a smooth frequency dependence on the scale of the free spectral range. Similarly, Z dis − iZ reac (ω − ω l )/ω l is the impedance of the array that terminates in the weak link, expanded to first order in both frequency around the resonance and real Z dis , the purely dissipative response at the resonance contained in Σ(ω). Here −iZ reac (ω − ω l )/ω l represents the reactive response in the vicinity of the resonance, i.e Z reac is real. The parameters X, Z reac and Z dis can be taken as frequency-independent in the vicinity of a resonance. The resonances thus have the line shape where Q i = Z reac /2Z dis and Q e = Z reac /Z tl are internal and external quality factors characterizing respectively dissipation in the weak link plus array, and in the external environment. This is the line-shape that we fit to the measured transmission resonances in order to extract the internal broadening γ int = ω n /Q i /(2π) and the precise resonance frequencies ω n /2π.

C Self-energy
Here we derive the self-energy expression (17) used in Sec. 4.2 to model the dissipative response of the BSG model. We perform the same calculation twice, using two equivalent approaches. In both cases, we perform Gaussian averaging of exponents whose arguments are linear in bosonic creation and annihilation operators. In the first derivation, we perform the required normal ordering by hand using Wick's theorem. In the second derivation, we represent the Wick contractions by Feynman diagrams. This is not as compact, but has the virtue of showing all multi-photon decay channels explicitly. The formal structure of our expansion is similar to that encountered for the bulk cosine nonlinearity in the Sine Gordon model so that the correctness of our result can be checked against results obtained in that context [58]. Subsequently, we discuss the self-consistent Born approximation.

C.1 First derivation
Since the perturbation contains cos φ 0 , a diagramatic representation of the expansion to second order will already contain an infinite number of diagrams. Fortunately the amputation process can be automated as follows. In the path-integral language, and in the time domain, the G φ 0 ,φ 0 (t 2 − t 1 ) Green's function, with external legs amputated, can be calculated from where /ℏ is the action associated with the weak link cosine perturbation, and ⟨. . .⟩ 0 denotes a Gaussian path integral over the field φ 0 (t), such that −i ⟨φ 0 (t)φ 0 (0)⟩ 0 = G 0 φ 0 ,φ 0 (t). Without the functional derivatives, the right-hand side of Eq. 35 would sum over all connected diagrams without external legs. The functional derivatives cut one bare φ 0 propagator of each such diagram, to produce two stubs, one at 0 and one at t, where external legs can be grafted. We expand (35) to second order in E J and go over from the path integral to the operator description where ⟨f f being any functional of the field φ 0 (t) and the right-hand side being the time-ordered expectation value of interaction-picture operators with respect to the zero-order ground state. We obtain Because the zero-order problem is harmonic, the field operatorφ 0 (t) is linear in boson creation and annihilation operators. One can expand the sin and cos functions of the field operators into exponentials. Under time-ordering, the product of exponentials of field operators equals the exponential of the sum of the operators. One is thus left with evaluating ⟨T exp X⟩ where X is linear in boson creation and annihilation operators. It is well known that the result is ⟨T exp X⟩ = exp T X 2 /2 . Thus one straightforwardly obtains for instance where is the "tree-level" vertex energy. Calculating the remaining expectation values in (36) in a similar manner, we find At second order in E J , the self-energy is related to the amputated Green's function through where Σ (n) (t) is the n'th order in E J contribution to Σ(t). We thus see that the linear in G 0 φ 0 φ 0 (t) part of sin G 0 φ 0 φ 0 (t) corresponds to the second term on the right-hand side of (40) and that

C.2 Second derivation
As an alternative to the above calculation, the self-energy can equivalently be represented diagramatically as follows. Expanding to second order, we construct all amputated diagrams with up to two vertices, that cannot be split in two by cutting an internal propagator: −i(Σ (1) + Σ (2) ) = + + + . . .
We have to organize this list. We focus on the effect of tadpoles. First, we single out a vertex. It can be dressed with any number of tadpoles, i.e. loops with a single propagator. We draw the rest of the diagram as a box, with any even number of lines between it and the singled out vertex. The sum over tadpoles reads where the whole sum as been absorbed into a new vertex, depicted as a grey disk. The factorization above worked because symmetry factors are multiplicative: if s is the symmetry factor of a diagram, the same diagram with n more tadpoles on some vertex will have symmetry s × 2 n n!. Thus, the dressed vertex equals E V0 J of Eq. 38. Using this vertex dressing, we reduced the list of diagrams to, −i(Σ (1) + Σ (2) ) = + + + + . . .
The symmetry factors in the second line and third lines work out such that the sums become respectively the cosine and sine of the propagator, with the leading term removed. Thus the second and third lines above exactly correspond to the second and third lines in Eq. 41. Because the Hamiltonian is even in the phases, and in particular in φ 0 , photon number parity is conserved, that is, one photon can decay into an odd number of photons only. This selection rule is fully respected by our calculation, and the diagrammatic representation of the self-energy is especially convenient to see this fact. Namely, non-linear vertices can have only an even number of legs, and in the sin(G) − G self-energy term which is responsible for the decay, each vertex (0 or t) has one and only one external leg, so the two vertices must be connected by an odd number of photonic lines.

C.3 Self-consistent Born Approximation
We can sum over a larger subset of diagrams by dressing the zero-order propagators appearing in the above result by all possible self-energy insertions. This takes into account that when a photon disintegrates at the weak link, it does not disintegrate into bare photon modes of the harmonic system, but into modes that are themselves hybridized with the weak link. If we ignore this dressing of propagators in the self-energy, the approximate interacting Green's function contains resonances when an incoming photon has a frequency equal to the sum of any (odd) n single-photon resonances of the harmonic zero-order problem. Given the nearly linear dispersion relation at low frequencies, this incorrectly predicts dense clusters of nearly degenerate n-photon resonances.
The dressing of propagators in the self-energy incorporates the level-repulsion between these resonances, which spreads them out over the free spectral range, thus giving a smooth background, rather than pronounced many-body peaks. The dressing of propagators in the self-energy leads to the replacement G 0 (41). However, if this is done blindly, there will be double-counting of some diagrams. For instance, because This happens because the full propagators used in the self-energy, built on a self-energy expansion to first order, already incorporates in an exact manner any terms in the perturbation that are quadratic in φ 0 . To cure the double-counting, we must therefore remove the second-order in G part of the cosine term in the dressed self-energy. Thus we arrive at the dressed self-energy of Eq. (17) in the main text, in which the dressed Green's function must be found self-consistently from the Dyson equation (19) in the main text.

D Keldysh technique
Finite temperature results are often obtained from diagramatic calculations that employ imaginary time Green's functions. To obtain the retarded Green's function, one has to perform an analytic continuation from imaginary to real time. This step is hard to perform numerically at the desired spectral resolution for our system with its many sharp resonances. We therefore rather use the Keldysh formalism in equilibrium, which does not involve such analytic continuation, to compute retarded Green's functions at finite temperature.
Instead of the time-ordered Green's function we employed previously, we have to use the contour-ordered Green's function, Here T c orders operators along a time contour with a forward branch σ = + from −∞ to ∞ and a backward branch σ = −, from ∞ to −∞. If (σ, σ ′ ) = (+, +) the largest time of t and t ′ is to the left. If (σ, σ ′ ) = (−, −), the largest time of t and t ′ is to the right. If (σ, σ ′ ) = (−, +), t is to the left, while if (σ, σ ′ ) = (+, −), t ′ is to the left. The expectation value refers to a thermal average. Thanks to the close similarities between contour ordering and time-ordering, the self-energy for G can be calculated using the same machinery as in Appendix C. In the path integral language there are independent fields φ 0± (t) associated with the forward and backward branches of the time contour. The weak link action is in analogy to the calculation in Appendix C. Summing the same class of diagrams as in Appendix C, one obtains The vertex energy is (Any component of G σσ ′ can be used in the vertex energy, since they are all equal at coinciding times.) In equilibrium, Σ σσ ′ (t, t ′ ) and G σσ ′ (t, t ′ ) only depend on the time difference t − t ′ , and can be Fourier-transformed from time-difference to frequency. The self-energy Σ and the Green's function G, viewed as 2 × 2 matrices with entries arranged according to ++ +− −+ −− , must be found self-consistently from Eq. (46) together with the Dyson equation in matrix form: In equilibrium, the contour-ordered Green's function is related to the retarded Green's function through so that After the self-energy Σ of the contour-ordered Green's function has been calculated, the retarded self-energy (that enters the admittance matrix) can be extracted through This self-energy is then substituted into Eq. (32) for the weak link impedance when the measured transmission signal is calculated.

E Scaling laws
In the main text, we remarked that the experimentally observed internal broadening of the chain modes decays exponentially as a function of frequency above ω ⋆ J . This is also what our microscopic theory predicts. However, in the theoretical literature, the BSG model is more usually associated with power-law dissipation. Here we review how the power law comes about, and explain why experimental realizations in cQED will have a hard time to exhibit such behavior.
The starting point is to assume that fluctuations of φ 0 are determined by the environmental impedance over a broad frequency range. In other words Z J ≫ Z env or ℏω p ≪ (2e) 2 /C J . Approximating the environment as an Ohmic impedance Z env (ω) = αR Q , yields a zero order (time-ordered, zero-temperature) Green's function where E ⋆ J is the effective weak link Josephson energy. At short times then where the omitted terms remain finite at small times. Using the bare propagator in the self-energy expression (17), the logarithmic divergence at small t then gives for ℏω ≫ 2παE ⋆ J and α < 1/2, where C is a constant with a positive real part. This power law would give the broadening of resonances a power-law frequency dependence ∼ ω 1−2α above the weak link resonance frequency. However, the above analysis ignores the finite ultraviolet cutoff of the physical system, typically set by the charging energy E C of the boundary junction. To see the predicted power law in an experimental realization would require a scale separation of several decades between 2παE ⋆ J /ℏ and the ultraviolet cutoff, as well as performing linear response measurements at frequencies that are orders of magnitude smaller than the ultraviolet cutoff, that is typically in the 10 1 GHz range, as illustrated by Fig. 10. Figure 11: Nonlinear damping γ J as a function of the frequency for Φ/Φ Q = 0.40 and 0.42. The dots correspond to the experimental data, the full lines are calculations from the perturbative treatment at second order, while the dashed ones correspond to the third order. At these lower magnetic fluxes, the circuit is in a regime where ℏω ∼ E ⋆ J . Hence, the perturbative approach fails, as the discrepancy between the second and third order shows. However, the theory is still able to reproduce the maximum of damping when ω ∼ ω ⋆ J F Perturbative breakdown at intermediate E ⋆ J Data corresponding to the low frequency range ℏω < E ⋆ J cannot be correctly described by the diagrammatic theory, as it is only valid for ℏω ≫ E ⋆ J . While the calculations matches quantatively the experimental data for the large flux values where E J is small enough (see Fig. 6), we see in Fig. 11 that the theoretical predictions (full lines) underestimate the measured losses by an order of magnitude for two smaller flux values. In order to confirm the non-perturbative nature of this discrepancy, we pushed the perturbative expansion to third order in E J .
In describing as simply as possible this higher class of diagrams, we only draw in what follows the diagrams with the lowest number of intermediate lines between each vertex, but have summed over all possible numbers of lines. This results in making the following replacements: G φ 0 φ 0 (t) → sin(G φ 0 φ 0 (t)) and G φ 0 φ 0 (t) 2 /2 → 1 − cos G φ 0 φ 0 (t). We also write here as double lines the full propagators of the skeleton expansion. The third order class of diagrams thus reads: Note that we did not include any nested diagram, since they are already generated by the skeleton expansion.
In Fig. 11 we compare the second order (full lines) and third order (dashed lines) diagrammatic results for the nonlinear damping rate at fluxes Φ = 0.40Φ Q and Φ = 0.42Φ Q with the corresponding experimental data (dots). There is a significant difference between the second and third order perturbative expansion implying that the expansion is not converged, although higher order corrections at least go in the right direction. At the larger fluxes of Fig. 6 where E J is markedly smaller, we find that the third order contribution is small compared to the second order one (provided ℏω ≫ E ⋆ J ), validating our theory. The failure of the diagrammatic apporach comes from the fact that, when ℏω ∼ E ⋆ J , the phase is partially trapped in the Josephson potential. Therefore, another source of damping, the phase slip between different minima of the Josephson potential, must be taken into account [44]. Since, the SCHA consists of replacing the cosine potential by an effective quadratic potential, these phase slip phenomena cannot be caught. However, although our theory is not quantitative in this regime, it correctly predicts that the maximum of γ J occurs at ℏω ∼ ω ⋆ J .

G Dielectric losses in the chain
In the remaining sections, we investigate whether more mundane loss mechanisms could provide an alternative explanation of our data. We start by considering the losses generated in the dielectric of the junction capacitances C of the chain, that can be modeled by writing that C = (ϵ ′ + iϵ ′′ ) d where ϵ ′ and ϵ ′′ are respectively the real and imaginary part of the dielectric permittivity while d is a parameter proportional to the length which depends on the capacitance geometry. ϵ ′ gives the capacitive response of C while ϵ ′′ is its dissipative part. Hence, the admittance of C is given by: where tan δ = Im(C)/Re(C) = ϵ ′′ /ϵ ′ . To find the damping induced by the dielectric, we use the dispersion relation 21), replacing C by C(1 + i tan δ) and taking the limit ka ≪ 1 (which is equivalent to the mode number n ≪ N , valid in the frequency windows that we probe). By defining the complex wavenumber κ = ka, and the dimensionless frequency x = ωl c /v φ , with l c the screening length and v φ the velocity of plasma modes, we have: Because of dielectric losses, the wavevector has a complex part: κ = κ ′ + iκ ′′ . We suppose that losses are weak enough so that tan δ ≪ 1 and hence κ ′ /κ ′′ ≫ 1. At first order in these quantities, we find: We then consider the small wavenumber limit |κ| ≪ 1/l c , that is equivalent to n ≪ N/l c ∼ 100, which describes the modes below 10 GHz as seen in Fig. 7 (above this frequency the dispersion relation starts to bend). If |κ| ≪ 1/l c , then x ≪ 1 and the chain behaves as an ideal transmission line sustaining TEM modes (see Eq. (59) with x ≪ 1 ) and the quality factor of the modes are given by [59]: Since 2πQ int = ω/γ int we end up with: Eq. (62) is used to fit simultaneously the internal damping for the magnetic fluxes Φ/Φ Q equal to 0, 0.2 and 0.3. For these three fluxes, the damping of the modes does not vary. Therefore, they do not appear to be caused by the SQUID nonlinearity. It has been noticed for chains of junctions [60] that tan δ has a slight frequency dependence which can be parametrized as: where A is the amplitude and b should be close to unity. The results of the fit where A and b are the free parameters is given in Fig. 12. The good agreement between the model and the data shows that for these magnetic fluxes the damping of the modes are dominated by dielectric losses in C. From that fit, we estimate that A × 2π × 1GHz = (3.4 ± 1.5).10 −4 and b = (0.5 ± 0.2). Hence, tan δ ∼ 10 −4 in the gigahertz range, which is consistent with what is found in similar devices, confirming that the dielectric is a good suspect for the damping observed in this range.

H Losses at the boundary junction
In this section, we will investigate whether another mechanism could explain the observed damping of the chain modes, via the dissipation coming from the capacitive or inductive part of the SQUID at the boundary junction. We saw in the previous section that the dielectric used in the junctions of the chain can generate a damping. The same effect can take place at the level of the boundary junction, that we model by adding a resistance in parallel to the junction capacitance such that: Since the SQUID itself is composed of two junctions, it can be expected that it can trigger damping of the circuit modes. On the other hand, since the circuit is superconducting, Then, any wave propagating into the circuit will be reflected respectively by a coefficient √ R 1 and √ R 2 (in amplitude) at site N and 0, because of the impedance mismatch with the measurement line or with Y ⋆ J . Let E em (0) be the electromagnetic energy stored in the circuit at a time t = 0. After a time t RT = 1/∆f FSR (the round trip time) the energy in the circuit is given by: E em (t RT ) = |R 1 (ω) R 2 (ω) |E em (0).
Hence, the energy decays exponentially with respect to time such that: where τ d is the characteristic damping time in energy. This damping time is inversely proportional to the damping frequency, such that: Because we want to estimate the internal damping frequency, we consider that the reflection is perfect at site N , |R 1 (ω) | = 1, while the reflection at site 0 is given by: whereZ C is the characteristic impedance of the chain, where the plasma frequency is taken into account. Hence,Z C (ω) = Z C / √ 1 − LCω 2 . Therefore the internal damping at the junction site is given by Now that we have an analytical formula relating R J to the damping frequency, we can use it with Eq. (64) and (65) to estimate the corresponding loss γ J . The results are displayed in the middle panel (for dielectric losses) and lower panel (for the quasiparticles) of Fig. 13. The estimates are given for the circuit parameters corresponding to Φ/Φ Q = 0.44, one of the flux for which γ J reaches its maximum.
We estimated the parameter tan δ coming from the dielectric of the chain junctions to be about 10 −4 . There is no physical reason why that of the SQUID should differ from it by several orders of magnitude. However, we see in the middle panel of Fig. 13 that even taking an unrealistic value of 5.10 −2 , γ J is underestimated the measured losses by about an order of magnitude. Finally, regarding the influence of non-equilibrium quasiparticles, the measured values ranges from 10 −8 to 10 −5 . However, we see in the lower panel in Fig. 13 that even a quasiparticle density x qp = 10 −3 gives a γ J that is off by two orders of magnitude. Hence, neither the dielectric nor the quasiparticles can realistically explain the order of magnitude of the measured mode damping γ J .

I Magnetic flux noise
The very convenient feature of a SQUID, its magnetic flux tunability has a price: its Josephson energy is sensitive to the noise of this control parameter. Since we saw that the frequency of the modes depends on the SQUID parameters, we expect the frequency of the modes to be time dependent. Hence, this noise could generate an inhomogeneous broadening of the modes. This broadening would then be more important where the modes are strongly influenced by the SQUID frequency, meaning those close to ω ⋆ J . Furthermore, the broadening is also expected to be larger for larger noise in E J . Hence, using Eq. (6) we see that the broadening should be larger close to Φ/Φ Q = 0.5 where E J depends strongly on magnetic flux (expect at the sweet spot induced by d but since it is very small in our case the sweet spot is also extremely small). These two properties for such an inhomogeneous broadening are qualitatively compatible with our observations. On the other hand, this This integral is obviously ultraviolet divergent. This is because such a noise is observed for a restricted frequency range, given by f ∈ [10 −4 Hz, 10 9 Hz] [65]. However, since we are looking for a frequency broadening, the low frequency cutoff must by, at least, larger than the integration bandwidth used to acquire the S 21 data, otherwise we would be able to measure the time dependence of δf n (Φ). We therefore set the low frequency cutoff to 1Hz. In addition, we set β = 1 for simplicity, as this will not influence future conclusions. With these assumptions we find: Finally using Eq. (72), (74), (76) and (79) together, we can estimate the broadening induced by a magnetic flux noise. The result of this estimation is given in Fig. 14. This study shows that despite a qualitatively good behavior with respect to flux and frequency, the broadening is more than two orders of magnitude smaller than the measured γ J . Hence, we can safely neglect the influence of a magnetic flux noise.