Emergent quantum phase transition of a Josephson junction coupled to a high-impedance multimode resonator

The physics of a single Josephson junction coupled to a resistive environment is a long-standing fundamental problem at the center of an intense debate, strongly revived by the advent of superconducting platforms with high-impedance multimode resonators. Here we investigate the emergent criticality of a junction coupled to a multimode resonator when the number of modes is increased. We demonstrate how the multimode environment renormalizes the Josephson and capacitive energies of the junction so that in the thermodynamic limit the charging energy dominates when the impedance is larger than the resistance quantum and is negligible otherwise, independently from the bare ratio between the two energy scales and the compact or extended nature of the phase of the junction. Via exact diagonalization, we find that the transition surprisingly stems from a level anticrossing involving not the ground state, but the first excited state, whose energy gap vanishes in the thermodynamic limit. We clarify the nature of the two phases by pointing at a different behavior of the ground and excited states and we show that at the transition point the spectrum displays universality not only at low frequencies. In agreement with recent experiments, we reveal striking spectral signatures of the phase transition.


Introduction -
The physics of a quantum system coupled to a bath with many degrees of freedom is a fundamental problem of many-body physics.One of the most famous examples is the spin-boson model for a two-level system coupled to a bath of harmonic oscillators [1,2], exhibiting a localization-delocalization phase transition.A related challenging problem is the resistively shunted Josephson junction [3].Schmid [4] and Bulgadaev [5] first addressed this problem by considering the junction as a fictitious particle in a periodic cosine potential, and the resistor as an infinite bath of harmonic oscillators.The junction was predicted to display a superconducting behavior (phase localization) for shunting resistances R < R q and an insulating behavior (phase delocalization) for R > R q , where R q = h/(2e) 2 is the resistance quantum [6], h the Planck constant and e the electron charge.This transition was predicted not to depend on the ratio between the Josephson energy E J and the charging energy E C .Following works using perturbative renormalization group [7][8][9][10][11], as well as Monte Carlo methods [12][13][14][15], further characterized the transition and confirmed its independence on E J /E C .
Experimentally, a few works investigated the predicted insulating phase [16][17][18][19], but were disputed in a recent experimental work [20] reporting the ac linear response of a resistively shunted junction and claiming no insulating behavior.This work revived the interest in this fundamental problem and sparked a lively debate.Yet, new theoretical analyses did not reach a consensus, with claims of absence [20,21], existence [22,23] and non-trivial modification of the predicted phase diagram [24] (although approaches based on the renormalization group are still a matter of debate [25][26][27]).Four key questions are outstanding: (i) which are the circuit Hamiltonian terms that are irrelevant in the thermodynamic limit?(ii) does a transition actually happen and how should one probe the so-called insulating phase?(iii) does the compact versus extended nature [28,29] of the Josephson junction phase play a role?(iv) Is the transition independent on E J /E C ?
The interest in this problem has also been catalyzed by impressive experimental advances in circuit quantum electrodynamics (QED) [30], particularly in the study of Josephson junctions coupled to high-impedance environments [31][32][33][34][35][36][37][38][39].A recent experiment [40] took a different approach: multimode transmission line resonators were used instead of resistors.Importantly, instead of measuring the transport properties of the junction, this experiment focused on the modifications of the linear-response spectral properties.The resonators used in these experiments are finite-size systems, in particular they have a finite free spectral range (mode frequency spacing).In this context, the thermodynamic limit corresponds to vanishing free spectral range.
In this work, we investigate theoretically the emergence of the transition for a Josephson junction coupled to transmission-line resonator when the free spectral range (number of modes) goes to zero (infinity), which, to the best of our knowledge, has never been explored.By approaching finite-size systems with a combination of analytical approaches and of numerical exact diagonalization, we show that for a homogeneous transmission line the transition emerges at R q /Z = 1, where Z is the line characteristic impedance, with the main features being present already for a small number of modes.At the transition point the system displays a universality independent on the ratio E J /E C , that we confirm with exact solutions not only at small E J .Moreover, unlike a wide class of phase transitions [41], we show that the emergent criticality involves a level anticrossing affecting the first excited state and not the ground state.We also show how the occurrence of the phase transition creates distinctive spectral signatures even at high frequencies, and how the ground state behaves differently from the excited states.
Quantum circuit model -Let us consider the multimode circuit QED system in Fig. 1(a), namely a Josephson junction coupled to a lumped-element transmission line.As detailed in the Methods, the quantum Hamiltonian has been derived from the circuit Lagrangian depending on C and Γ, respectively the capacitance and inductance matrices.The system degrees of freedom are described by the vector φ, containing independent flux variables [42], indicated in Fig. 1(a).To get the final form, we have performed a simultaneous diagonalization of C and Γ via the basis change φ = P ϕ.Through a Legendre transform and identifying as the junction phase φ = k P 0k ϕ k , we have obtained the classical Hamiltonian (for an analogous procedure see [43]).Finally, by quantizing the degrees of freedom, we get the quantum Hamiltonian in the charge gauge where ĤJ = 4E C N 2 −E J cos φ and â † k is the bosonic creation operator for the kth mode and g k = 2 √ ℏω k E C P 0k its coupling to the junction.Here we used the following sum rule, due to the normalization identity P T CP = I, ( The frequencies and couplings in (1) depend on the capacitive and inductive matrices.Illustrative values are Nm , where N m is the number of modes.The phase transition is expected to emerge in the thermodynamic limit N m → ∞, that is equivalent to ∆ → 0.
Renormalized Josephson and charging energies -It is insightful to consider the basis diagonalizing the last two terms of (1), namely [44] for the kth mode with α N,k = iN g k ℏω k , and n is a vector with the number of photons in each mode.The only nontrivial term in (1) is the junction Hamiltonian: the dressing of the junction by the line can be determined by inspecting its matrix elements in such a basis.As detailed in the Methods, the coupling to the line effectively produces a renormalized Josephson energy independent of the number of photons in the modes.The remaining effect is a renormalized capacitive energy P 00 vanishes in the thermodynamic limit and has been taken to be exactly zero in a recent study of the phase transition [24].However, the way it approaches zero in the thermodynamic limit is important and cannot be neglected.Indeed, also the renormalized Josephson energy E J tends to zero for ∆ → 0. The properties of the infinite system will hence be determined by the asymptotic behavior of their ratio.We computed such analytical expression by injecting the numerically evaluated values of g k , ω k and P 00 for different values of the impedance Z and decreasing values of ∆.The results are summarized in Figure 2. One can see that both E C and E J decrease for all values of R q /Z while ∆ is decreased (Fig. 2(a)-(b)).However, while E C monotonically decreases while increasing R q /Z, E J first decreases and then increases.The ratio of the two quantities E J / E C instead displays a striking behaviour (Fig. 2(c)): for small enough ∆ it sharply decreases from 1 to very small values and then grows again.Remarkably, the curves for different values of ∆ all cross for Z = R q , with the curves for larger sizes (smaller ∆) tending to smaller values for R q /Z < 1, and to larger values for R q /Z > 1.The inset Fig. 2(d) shows that the slope of the curves at Z = R q grows logarithmically with ∆ −1 .We can conclude that in the thermodynamic limit the ratio E J / E C tends to 0 for R q /Z < 1 (apart from R q /Z = 0, for which no renormalization can occur) and to infinity for R q /Z > 1.
This points at the two expected phases: for R q /Z < 1 the charging energy dominates (so-called insulating behaviour), for R q /Z > 1 the charging energy is negligible (so-called superconducting behaviour).We emphasize that our theoretical analysis predicts a singular point for Z = R q , independently of the bare ratio E J /E C , and that the compact or extended nature of the junction phase does not play a role here.
Exact diagonalization -To complete our investigation and fully characterize the emergence of the two phases we have applied exact diagonalization techniques [45,46] to our Hamiltonian (1) for finite-size systems.To extend the reachable sizes, we have taken full advantage of Josephson potential periodicity and of Bloch theorem.Indeed, we can introduce the eigenstates of the junction Hamiltonian |ν, s⟩ = e iνφ u s ν (φ), where ν is the quasi-charge and s is the band index.In our charge gauge representation, the full Hamiltonian is block diagonal and for each ν we need to diagonalize the following matrix: where ε s ν are the eigenenergies of the bare junction with ν ∈ [−0.5, 0.5] (Brillouin zone).Note that the spectrum is even with respect to ν, so we can restrict to positive values only.
Illustrative energy bands for the bare junction are reported in Fig. 1(b).An example of the obtained energy bands for the interacting system are shown in Fig. 3 versus Z. Different colors correspond to different values of the free spectral range for the same ω p (i.e. to different numbers of modes), and the energies are rescaled by ℏ∆.For R q /Z → 0, the junction and the modes are decoupled: the spectrum is given by the bare junction bands with replicas due to a finite number of bosons in the modes.For finite R q /Z instead, the interaction manifests in energy anticrossings.To investigate a quantum phase transition, the behavior of the ground state and of the low-lying excited states is crucial.While increasing R q /Z, the first energy band becomes narrower, corresponding to a reduced effective charging energy (see Fig. 2(a)).However, the curvature of the first band cannot change, since the sum rule (2) implies that E C ≥ 0. This means that the ground state is always at ν = 0.The first excited band instead changes its shape while varying Z.When the system is not coupled (R q /Z = 0), it is simply a one-photon replica of the first band.Instead, when R q /Z is increased, the slope changes (see last two panels of Fig. 3).The low-energy spectrum of the full system at low impedances is reminiscent of the bare junction spectrum, although at much smaller energies.This qualitative behavior of the first few bands occurs for different sizes, although at smaller energies for smaller ∆ (different colors in Fig. 3).Interestingly, the rescaled bands overlap in the best way for Z = R q , confirming the universality observed in Fig. 2(c).Away from this point, for both larger and smaller values of Z, the rescaled spectra do not overlap any longer.This is reminiscent of the finite-size scaling of continuous phase transitions.
Compact versus extended Josephson phase-In the exact diagonalization studies, we considered an extended phase and applied Bloch's theorem.Now, since the full Hamiltonian for our finite-size system does not couple different quasi-charges ν, the results for a compact phase are independent from the phase being fundamentally extended or compact.Indeed, we found that the ground state always occurs for ν = 0, which corresponds to a periodic wavefunction [47].This independence of the phase transition on the extended/compact nature of the junction phase was also discussed in [23].
Behavior of ground and excited states -In Fig. 4(a) we plot the ground state energy versus R q /Z: surprisingly, it has a weak dependence on the system size, and no precursor of a critical behavior is apparent.Fig. 4(b) instead reports the energies of the first few excited states, separated from the ground state by an energy of the order of ℏ∆ for all the values of Z.In the thermodynamic limit, the system is expected to be gapless for all values of Z: hence, the phase transition needs to emerge with a different mechanism.Indeed, around Z = R q the first three excited levels exhibit an anticrossing.This corresponds to the observation we already made for Fig. 3 that the first excited band passes from being a one-photon replica to a dressed Josephson band.For a compact phase, this occurs around Z = R q , with the splitting of the resulting anticrossing (but not the position) depending on E J , as displayed in Fig. 5(a).As in Fig. 3, while decreasing ∆ all the rescaled spectra exhibit an universal behavior at Z = R q .Hence, when approaching the thermodynamic limit, for a fixed E J the anticrossing becomes narrower and narrower, meaning that the signature of the phase transition at first emerges in the first excited states.At the same time, the spectrum is gapless in the ∆ → 0 limit, so that a singular behaviour in the first excited state necessarily reflects on the ground state.
While a compact phase is the only relevant quantity for the ground state of the system, the response of the system to external perturbations can depend on the full bands.To better highlight the behaviour of the low-lying bands, in Fig. 4(d) we show the behaviour of the gap between the first and second bands at the edge of the Brillouin zone.The ratio of the gap and the FSR displays a behaviour analogous to the renormalized parameters of Fig. 2, with the gap closing for larger systems for R q /Z < 1.This means that the low energy spectrum is approaching a capacitive "free particle" behaviour.The opposite is true for R q /Z > 1.As shown in Fig. 5(b), this behaviour extends to larger values of E J .
Importantly however, the ground state does not behave in the same way, as can be seen for example by computing observables for the junction degrees of freedom.As an example, we show in Fig. 4(c) the charge fluctuations σ 2 N = tr(ρ J N 2 ) − tr(ρ J N ) 2 , with ρJ the reduced density matrix for the junction (analogous results are obtained, for example, for tr(ρ J cos φ)).For R q /Z > 1 the charge fluctuations have a strong dependence on the system size.For R q /Z < 1 instead, the size dependence is much smaller and the fluctuations never fall below the Cooper pair box value that is obtained for R q /Z → 0 (horizontal dashed line) [48].This how- ever does not mean that the transition does not affect the ground state, as signalled by the different size dependence of the charge fluctuations on the two sides of R q /Z = 1.
The emergence of the singularity in the ground state in the thermodynamic limit is however slow and we can only see a precursor.
We believe that this peculiar difference of behaviour between the ground state and the excited states is at the origin of much of the recent controversy over this phase transition, and specifically over the characterization of the so-called insulating state.In particular, in the thermodynamic limit, the ground state is not approaching a purely capacitive behaviour.However, the response of the system to a gate charge (i.e. to a change in the quasi-charge) changes sharply at the transition point.

Effect on higher frequency modes -
To understand what happens to the higher frequency modes at the transition point, we focus on the compact-phase Hamiltonian, given by equation ( 5) for ν = 0.For E J = 0, the excited bands energies are ε n 0 (E J = 0) = 4 E C n, with n a positive integer.In this limit, the compact-phase excited levels can be understood by simply plotting ε n 0 (E J = 0) and the Z-dependent mode frequencies.We show these two energy branches in Fig. 6(a).For energies small enough with respect to ℏω p , the nth photon mode energy crosses the nth junction renormalized band at Z = R q .This degeneracy is lifted for finite E J , resulting in the anticrossings of Fig. 4(b) and Fig. 5(a).Hence, around the transition point all the modes will hybridize with the corresponding dressed band.In the limit ω p → ∞ (a non-dispersive line) this would hold for all the energies, while for finite ω p there is eventually a shift at high energies.
Spectroscopy of high modes -Recent experiments have observed a striking signature of the transition from the dispersion of high-frequency modes [40].With the results of our exact diagonalization, we can calculate the linear- response photonic spectral function ) where E n (|E n ⟩) is the nth excited eigenenergy (eigenstate) of the full system, E G (|G⟩) the ground energy (state) and γ is a phenomenological broadening.Here, as in the experiment, we have considered a SQUID that is equivalent to a junction with E J that can be tuned by a magnetic field.Illustrative spectra versus the external magnetic flux Φ are shown in Fig. 7.We observe a clear change of sign of the photon frequency dispersion across the transition, as observed in [40].This effect was linked to a capacitive/inductive behavior of the junction and here we see that it originates from the fact that at the transition point the nth dressed junction band crosses from above the nth photonic mode frequency [49].In Fig. 7, one can see this dressed-junction state acquiring a finite photonic component near the crossing point for Φ/Φ 0 ̸ = ±0.5.For larger and smaller impedances it is dark, i.e. it does not contribute to this spectral function.
Conclusions -In conclusion, we have revealed through an original analytical approach and exact diagonalization the emergence of the phase transition for a Josephson junction induced by the coupling to a high-impedance environment represented by a multimode transmission line.Our exact approach confirms hence the existence of the much disputed Schmid-Bulgadaev phase transition.Our study has evidenced unconventional properties of such phase transition and distinctive signatures on a high-frequency spectroscopy.Notice that the subtle connection with the superconducting/insulating transport properties of the present phases is not so obvious since this requires a coupling to another circuit, a problem that needs to be addressed in the future.However, the current disagreement over the so-called insulating state may be explained by the different behaviour of the ground state and the excited states we evidenced.The present work also paves the way to investigate emergent quantum phase transitions in complex environments with arbitrary transmission lines that can be engineered in a controlled way.

Derivation of the circuit QED Hamiltonian
We detail here the intermediate steps in the derivation of the Hamiltonian corresponding to the circuit of Fig. 1(a), that can be described in terms of the independent flux node variables [42] included in the vector φ.The procedure is analogous to the one used in [43].The first step is the circuit Lagrangian where C and Γ are the capacitive and inductive matrices.Note that the pre-factors have been introduced to make C and φ dimensionless.Their expressions read: (S.8) The Josephson junction charging energy is given in terms of its capacitance C J by E C = e 2 /2C J .
The second step is to re-express such Lagrangian in terms of the normal modes.To get this, we perform a simultaneous diagonalization of the C and Γ matrices, given by a change of basis φ = P ϕ that satisfies ΓP = CP ω 2 , with ω a diagonal matrix.The Lagrangian can therefore be written as (S.9) Introducing conjugate momenta (charges) N k = ∂L/∂ φk , the third step is to perform a Legendre transformation to get to the Hamiltonian form.It is convenient to change variables, so to identify the argument of the cosine as the junction phase φ = k P 0k ϕ k .This requires a redefinition of the junction charge N = N 0 /P 00 and of the mode charges n k = N k − P 0k N .After this fourth step, we get the Hamiltonian (S.10) By performing a canonical quantization of the Hamiltonian variables, by introducing bosonic creation and annihilation operators for the modes, and using the sum rule (2), one finally obtains the quantum Hamiltonian (1).

Gauges and domain of the junction phase
We provide here additional details and remarks about the choice of the circuit gauge and the domain (extended versus compact) of the Josephson junction phase.With the diagonalization procedure detailed in the previous section, the Hamiltonian is in the so-called charge gauge, because the transmission line operators are coupled to the junction via its charge operator.This gauge is analogous to the Coulomb gauge Hamiltonian in cavity QED.Notice however the absence of a diamagnetic term, that has been already implicitly eliminated in the diagonalization procedure.This is equivalent to starting with the diamagnetic term and diagonalizing the modes with a multi-mode Bogoliubov transformation.One could also consider a flux (dipole) gauge Hamiltonian, that can be obtained with a partial diagonalization of the Lagrangian or with a unitary transformation from the charge gauge one.In the flux gauge, the coupling between the junction and the modes is through the junction phase φ.In this work we have chosen to work with the charge gauge Hamiltonian, because it is block diagonal once we apply Bloch theorem.
Notice that a purely charge gauge Hamiltonian is possible if the first eigenvalue coming from the Lagrangian diagonalization is zero.This occurs only if the boundary condition at the end of the array is chosen to be an open one.This is due to the fact that the circuit we consider has de facto a superconducting island in the upper part of the array.If one also had, for example, an extra inductance in parallel to the whole array, the charge in the island would no longer be conserved and additional terms coupling to the flux of the junction would appear in equation (S.10).This subtle point is related to the discussion on the domain over which the junction phase should take values [23,29].In principle, values of the phase differing by 2π should correspond to the same physical state, and hence φ should be a compact variable over one period of the cosine: this corresponds to a quantized conserved charge across the junction.However, once the junction is coupled to another circuit, this compactness is no longer obvious (notice also that we needed to take a linear combination of variables to identify the phase in the Hamiltonian (S.10)).Many treatments, including the ones predicting the Schmid-Bulgadaev transition, consider it to be extended over the whole real axis.In Ref. [20] it was argued that the transition is an incorrect prediction deriving from not considering a compact phase.In our work, we have demonstrated how the phase transition emerges both for a compact and extended phase.

Renormalization of the junction
We detail here some intermediate technical steps about the derivation of the renormalized Josephson energy and capacitive energy due to the coupling of the junction to the multi-mode transmission line.Let us consider the matrix elements of the Hamiltonian (1) expressed in the basis diagonalizing the sum of the last two terms of the Hamiltonian (they contain the bosonic operators for the transmission line modes).For each mode, one can introduce a shifted bosonic operator bk = âk + i N g k ℏω k such that (S.13)The overlap of two displaced number states has the following expression, that can be obtained from the matrix elements of the displacement operator between number states [44], for m > n, where L (m) n are the associated Laguerre polynomials.The matrix elements (S.13) can hence be expressed as (S.15) Notice that for N = M the bare junction Hamiltonian matrix element is not modified.Hence the dressing of the junction by the transmission line only takes place for off-diagonal matrix elements in charge space, that is for the Josephson potential term.
Even if this expression is apparently complicated, there is an easily readable important effect on all the matrix elements.Since the Josephson potential can be also written as While the charging kinetic term of the junction is not directly affected by the dressing in (S.15), the last term of (S.12) proportional to N 2 modifies it.This effect is analogous to the polaron dressing increasing the effective mass of a particle coupled to bosonic fields.In the present case, this corresponds to a decrease of the effective charging energy.Given the sum rule (2), the residual charging energy is given by equation (4).
Notice that, while the sum appearing in the shift of the kinetic term is bounded because of the sum rule (2), the sum appearing in the exponential renormalizing the Josephson energy (3) is not, and will diverge in the thermodynamic limit ∆ → 0. This determines an exponential suppression of all the matrix elements (S.15), irrespective of the occupation numbers of the different modes.This is therefore the dominant effect in the thermodynamic limit.

Details about the exact diagonalization
We report here some technical details about the exact diagonalization.All numerical computations were performed with the Julia programming language [46].We have diagonalized the Hamiltonian (5) at fixed quasicharge ν, taking advantage of the fact that the full Hamiltonian is block diagonal in the quasi-charge.We repre- sent (5) on the basis of eigenstates of the uncoupled system {|ν, s⟩ Nm k=1 |n k ⟩}, where s labels the bands and |n k ⟩ is a Fock state of the kth mode.A finite basis is obtained by taking a fixed number of bands s ∈ [1, N bands ] and by introducing an energy cutoff limiting the photonic states to the ones satisfying j ℏω j n j < E cut .Both E cut and N bands are numerical parameters that we need to take large enough to ensure convergence of the results.The convergence with respect to both parameters is shown in Fig. 8 and Fig. 9.
To determine the Hamiltonian matrix to diagonalize, we have first solved the Schrödinger equation for the bare junction using Bloch's theorem in order to obtain the bare energy spectrum and the matrix elements of the junction charge operator.Then, we have numerically diagonalized the capacitive and inductive matrices to obtain the impedance-dependent frequencies of the transmission line modes and their coupling to the junction.
Since the size of the basis grows exponentially with the system size, in order to push to the limit the exact diagonalizations it is important to have an efficient way of generating it given the energy cutoff constraint and to fill the (sparse) Hamiltonian.We have used numerical alogorithms similar in spirit to those described for example in [45]: in particular, we have generated the basis incrementally by introducing a lexicographic order and filled the Hamiltonian by searching through ordered tags that uniquely identify each state.We have finally diagonalized the resulting sparse matrix by employing libraries that perform iterative solutions for computing the lowest eigenvalues.
Notice that the size of the basis increases quickly when increasing the number of modes and the numerical cutoff.The number of modes for which we performed exact diagonalization is hence moderate with respect to the ones considered in the computation of the renormalized Josephson and charging energies.It is interesting that the universality of the spectrum, that is expected to hold for large enough systems, is already well visible at these moderate sizes.

FIG. 1 .
FIG. 1.(a) Representation of a circuit consisting of a Josephson junction (red) of Josephson energy EJ and capacitive energy EC coupled to a transmission line resonator.The fluxes used to describe the circuit are indicated with black dots.(b) Energy bands of the uncoupled junction in terms of the quasicharge ν for EJ = 0.5EC .(c) Frequencies of the transmission line modes versus the mode number.(d) Junction-line couplings in the charge gauge as a function of the mode frequency for different values of the impedance Z.For these last two panels: plasma frequency ℏωp = 5EC and free spectral range ℏ∆ = 0.2EC (see definition in the text).

FIG. 2 .
FIG. 2. Upper diagram: in the thermodynamic limit the full system can be mapped into a renormalized junction.(a)-(b) Renormalized capacitive and Josephson energies computed from (4) and (3) for different values of ∆.(c) Ratio of the renormalized Josephson and charging energies.(d) Slope of the ratio at Z = Rq showing logarithmic growth in the thermodynamic limit.For these plots EJ = EC and ℏωp = 2EC .

FIG. 3 .
FIG.3.Exact diagonalizaton results for the full system.The different panels correspond to different impedances Z and the different line colors to different free spectral ranges (see legend).The vertical axis is rescaled by ℏ∆.In these plots EJ = 0.5EC and ℏωp = 2EC .

FIG. 4 .
FIG. 4. (a) Ground state energy for different free spectral ranges.(b) Corresponding dependence of the first excited state energies for a compact phase (ν = 0).(c) Junction charge fluctuations in the ground state for different free spectral ranges.(d) Band gap at the edge of the Brillouin zone ν = 0.5.In all the panels EJ = 0.5EC and ℏωp = 2EC .

FIG. 5 .
FIG. 5. (a) Dependence on the Josephson energy EJ of the anticrossing between the first few excited states for a compact phase.Here ℏ∆ = 0.55EC and ℏωp = 2EC .(b) Band gap at the edge of the Brillouin zone showing the universality at Rq/Z = 1 for a wide range of EJ .Solid lines ℏ∆ = 0.55EC , dashed lines ℏ∆ = 0.41.

11 )
The eigenstates of this part of the Hamiltonian are the tensor product of a charge eigenstate for the junction and of number states for the displaced operators bk .Note that these are displaced number states for the original operators âk , namely the states|N, n⟩ = |N ⟩ Nm k=1 |α N,k , n k ⟩already introduced in the main text.The Hamiltonian matrix elements between these states are hence ⟨M, m| Ĥ |N, n⟩ = δ(M − N )δ m,n k ℏω k n k + ⟨M, m| ĤJ |N, n⟩ − δ(M − N )δ m,n of the junction Hamiltonian being ⟨M, m| ĤJ |N, n⟩ = ⟨M | ĤJ |N ⟩ Nm k=1 ⟨α M,k , m k |α N,k , n k ⟩ .

(S. 16 )
only matrix elements with |N − M | = 1 are nonzero.Hence the exponential Franck-Condon-like factor before the product in equation (S.15) can be interpreted as giving the renormalized Josephson energy (3), independently from the number of photons in the different modes.