Optimized heat transfer at exceptional points in quantum circuits

Superconducting quantum circuits are potential candidates to realize a large-scale quantum computer. The envisioned large density of integrated components, however, requires a proper thermal management and control of dissipation. To this end, it is advantageous to utilize tunable dissipation channels and to exploit the optimized heat flow at exceptional points (EPs). Here, we experimentally realize an EP in a superconducting microwave circuit consisting of two resonators. The EP is a singularity point of the Hamiltonian, and corresponds to the most efficient heat transfer between the resonators without oscillation of energy. We observe a crossover from underdamped to overdamped coupling via the EP by utilizing photon-assisted tunneling as an \emph{in situ} tunable dissipative element in one of the resonators. The methods studied here can be applied to different circuits to obtain fast dissipation, for example, for initializing qubits to their ground states. In addition, these results pave the way towards thorough investigation of parity--time ($\mathcal{PT}$) symmetric systems and the spontaneous symmetry breaking in superconducting microwave circuits operating at the level of single energy quanta.

Systems with effective non-Hermitian Hamiltonians have been actively studied in various setups in recent years [1][2][3][4][5][6][7]. They show many intriguing properties such as singularities in their energy spectra [8][9][10][11][12][13]. A squareroot singularity point in the parameter space of a non-Hermitian matrix is called an exceptional point, EP, if the eigenvalues coalesce [12,13]. Previously, EPs have been shown to emerge, for example, in non-superconducting microwave circuits, laser physics, quantum phase transitions, and atomic and molecular physics [12,13]. The fascinating effects of EPs include the disappearance of the beating Rabi oscillations [14], chiral states in microwave systems [15], and spontaneous symmetry breaking in systems with parity-and time-reversal (PT ) symmetry [16][17][18][19]. In the quantum regime, PT -symmetric systems may show features that are different from the semiclassical predictions, such as new phases owing to quantum fluctuations [19,20]. Despite the active research on EPs, they have not been thoroughly investigated in superconducting microwave circuits to date [21].
Superconducting microwave circuits provide an ideal platform to realize various quantum devices, such as ultrasensitive photon detectors and counters [22][23][24][25], and potentially even a large-scale quantum computer [26,27], or a quantum simulator [28] in the framework of circuit quantum electrodynamics [29,30]. Notably, superconducting qubits have been shown to approach the required coherence times [31,32] for quantum error correction [33,34]. However, despite the tremendous interest in superconducting microwave circuits in recent years [35][36][37], there are still many issues to be solved before a fully functional quantum computer is possible. For example, the precise engineering of energy flows between different parts of the circuit in scalable architectures is of utmost importance since unwanted heat is a typical source of decoherence in qubits [38,39]. In many error correction codes, qubits are repeatedly initialized, which requires fast and efficient cooling schemes [40][41][42]. One promising method for absorbing energy and initializing qubits to their ground state is based on resonators with tunable dissipation [43].
The recently developed quantum-circuit refrigerator (QCR) [44,45] provides great potential for both qubit initialization and thermal management since it enables tunability of energy dissipation rates over several orders of magnitude in a superconducting microwave resonator [46].
Operation of the QCR relies on inelastic tunneling of electrons through a normal-metalinsulator-superconductor (NIS) junction [47]. The tunneling electrons can absorb or emit photons to a resonator which allows to control the coupling strength to a low-temperature bath in situ. This tunable coupling strength has also been shown to induce a broadband Lamb shift [46]. Furthermore, elastic tunneling can be utilized for temperature control of the normal-metal electrons [48][49][50], and for precise thermometry down to millikelvin temperatures [50,51]. Recently, NIS junctions have also been utilized in a realization of a quantum heat valve [52], and phase-coherent caloritronics [53].
In this work, we combine the advantages of tunable dissipation and EPs to optimize the heat flow in a superconducting microwave circuit. To this end, we investigate a circuit consisting of two coupled resonators, one of which is equipped with NIS junctions and a fluxtunable resonance frequency (Fig. 1). We denote the NIS junctions and the normal-metal island that is capac-itively coupled to the resonator as a QCR. Thanks to the voltage-tunable dissipation within the QCR and the fluxtunable resonance frequency, an EP arises in the Hamiltonian that describes the modes of the coupled resonator system. We investigate the emergence of the EP using frequency and dissipation as control parameters (Fig. 2) and verify its properties experimentally by measuring the microwave transmission coefficient (Fig. 3). The optimal heat flow given by the coupling strength can be reached at the EP (Fig. 4). Different types of tunable resonators have been studied in recent years [42,[54][55][56][57][58][59][60][61] but not with voltage-tunable dissipation. Our work demonstrates a platform to control the local heat transport between neighboring nodes in a quantum electrical circuit. In addition to thermal management within superconducting multi-qubit systems, these methods may be applicable to thermally assisted quantum annealing [62] and to studies of the eigenstate thermalization hypothesis in many-body quantum problems [63]. Furthermore, our work is an important step towards the investigation of PT -symmetric systems at the quantum level that can be realized with circuit quantum electrodynamics architectures [21,64].

Experimental samples
Our samples consist of two coplanar waveguide resonators, R1 and R2, which are capacitively coupled to each other, as depicted in Fig. 1(a), (b) (see also Supplementary Fig. S1). The resonator R1 has a fixed fundamental frequency at 2.6 GHz. This mode does not couple strongly to the resonator R2 owing to the voltage node in the middle of the resonator R1 where the coupling capacitor C C is located. Therefore, we focus on the first excited mode of R1, with frequency f 1 = ω 1 /(2π) = 5.2 GHz, which has a voltage antinode at the coupling point of the resonators. The resulting capacitive coupling between the resonators has a strength g/(2π) = 7.2 MHz. In contrast to R1, the resonator R2 has a flux-tunable resonance frequency ω 2 (Φ) owing to a superconducting quantum interference device (SQUID), and a voltage tunable loss rate κ 2 (V b ) owing to the QCR. Here, Φ and V b are the magnetic flux applied to the SQUID loop and the voltage bias of the QCR, respectively. The inductance of the SQUID and, hence, also the resonance frequency of R2 are periodic in flux with a period of the flux quantum Φ 0 = e/(2h). Consequently, due to the coupling of the resonators, R1 also shows flux-dependent features. We show the QCR in Fig. 1(c) and schematically present its operation principle in Fig. 1(d). The difference in photon absorption and emission rates originates from the gap of 2∆ in the density of states of the superconductor, and the difference can be utilized to cool down quantum circuits [44,47]. The sample fabrication  Supplementary Table S1.

Exceptional points
To study exceptional points, we utilize two control parameters in the effective Hamiltonian of the system: we use the voltage-tunable dissipation rate κ 2 (V b ) and the flux-tunable detuning between the resonators δ = ω 2 (Φ) − ω 1 . We study the system consisting of the two resonators in a rotating frame with a frequency corresponding to the uncoupled mode frequency of R1, similarly as in Ref. [3]. Thus, the excitations of the system can be described with the effective non-Hermitian Hamiltonian in matrix form in the basis ψ = (A, B) T where A and B are field amplitudes in R1 and R2, respectively, (see Methods and Ref. [60]) where κ 1 is the decay rate of the resonator R1. The eigenvalues of H can be written as and the corresponding eigenvectors are where Thus, the eigenvalues and eigenvectors coalesce when the square-root term s vanishes resulting in an EP. Consequently, there is only a single eigenvalue and, importantly, there is also only a single eigenvector. The EP occurs at |κ 2 − κ 1 | = 4g, and δ = 0. In the following, we assume that κ 1 κ 2 , which is valid for our samples, as discussed below. Consequently, the condition for the EP simplifies to κ 2 = 4g.
To visualize the system singularity, i.e., the EP, we show the real and imaginary parts of λ ± in Fig. 2. The eigenvalues form a self-intersecting Riemann surface in the parameter space of κ 2 and δ. The imaginary part corresponds to mode decay, and real part to mode frequency deviation from the uncoupled mode frequency of R1. Our system consisting of the resonators R1 and R2 can be considered as a single damped harmonic oscillator, where the energy oscillates between the two resonators. In the underdamped case, κ 2 < 4g, the modes have an equal decay rate at zero detuning, and there is an anti-crossing of the mode frequencies. In contrast, in the overdamped case, κ 2 > 4g, there is an anticrossing in the mode decay rates as a function of the detuning, and the mode frequencies are equal at zero detuning. Consequently, one of the modes remains lossy whereas the other one has a low decay rate at different detunings.
Let us connect the meaning of this critical point to the efficiency of energy transfer between the two resonators. In terms of coupled dissipative systems, the EP separates the system between the overdamped and underdamped regime being the point of critical coupling. It follows from the dynamics of the coupled system that at this point the energy is transfered between the two resonators optimally fast without back and forth oscillation [42,60]. In particular, the rate of heat transfer at zero detuning is given by Here, the branch with the upper signs corresponds to a mode located predominantly in the primary resonator R1, and the branch with the lower ones in the dissipative resonator R2 ( Supplementary Fig. S2). Consequently, by reaching the EP at κ 2 = 4g, we operate our sample at a point of optimally efficient and nonreciprocal heat transfer out of R1.

Experimental observations
To explore the dissipative dynamics of the two coupled resonators, we measure the flux-and frequencydependent scattering parameter S 21 describing the transmission from Port 1 to Port 2 for different bias voltages using a vector network analyzer. We tune the magnetic flux in a range where the frequency of R2 crosses that of R1. As shown in Fig. 3(a), we observe a transition from anticrossing to a single mode already indicating the presence of an EP in-between. A broader range of bias voltages is shown in Supplementary Figs. S4 and S5. To generate a quantitative description of our system, which is required for the investigation of EPs, we numerically simulate the scattering coefficient as shown in Fig. 3(b). Here, we model the SQUID as a flux-tunable inductor, and the QCR as an effective resistance R eff (see Methods and Supplementary Fig. S3). We extract R eff by fitting the circuit model to the experimental results, and use R eff to obtain the damping rates of the dissipative resonator R2 as a function of the bias voltage (Methods). The extracted resistances are given in Supplementary Fig. S6. The model in Fig. 3 is in very good agreement with the experimental results. In addition to R eff , we also extract the coupling capacitance and the critical current of the SQUID from the simulation. The coupling capacitance is found to be 3.8 fF which agrees well with the finite-element-method simulation that yields approximately 5 fF. In Fig. 3(a) the crossing of the modes as a function of the flux shifts slightly towards lower flux values. We attribute this shift to heating of the SQUID, which results in a reduced critical current of the SQUID, and hence a larger inductance and lower resonance frequency of R2. In principle, we vary also the Lamb shift [46], which, how- The EP is obtained at the intersection of κ2 and the critical coupling 4g. Furthermore, the figure shows the theoretical coupling strength of the QCR, κQCR, without taking dephasing and other voltage-dependent losses into account. The uncertainty of the data points is of the same order as the marker size. (b) Effective decay rates of the coupled system κ eff = −2Im(λ±) calculated from κ2 at zero detuning using Eq. (2). The two branches at high voltages correspond to λ+ and λ− with the modes located predominantly in one of the resonators as indicated. The maximum decay rate for R1 is obtained at the EP. The damping rates of the modes are equal at eV b /(2∆) < 1 due to hybridization. ever, causes only a minor effect since the resonance of R2 is very broad, and therefore, we neglect it in our model. To verify our above assumption κ 1 κ 2 , we measure the internal quality factor of the primary resonator R1 with R2 far detuned at Φ = Φ 0 /2. From the quality factor, we extract a loss rate κ 1 /(2π) 260 kHz for both samples which is substantially lower than the extracted value of κ 2 . We measure a slight temperature and power dependence of the quality factor of R1 as expected in the case of two-level fluctuators dominating the losses [65] ( Supplementary Fig. S8).
To demonstrate the presence of an EP, we show the extracted damping rates (Methods) of the bare resonator R2 in the absence of R1 in Fig. 4(a) as functions of the bias voltage for Samples A and B. Both samples have similar damping rates κ 1 and κ 2 . The value of κ 2 can be tuned by approximately two orders of magnitude. At low bias voltages eV b /(2∆) < 1, the rate κ 2 is below 4g, and at eV b /(2∆) > 1 the damping rate exceeds 4g. We describe the origin of the tunable damping rates using a model that contains the photon absorption and emission at the QCR given by the rate κ QCR , as well as constant internal losses κ int,2 , and voltage-dependent residual losses κ r,2 . We have designed our sample in such a way that κ QCR covers the critical damping rate, and thus the losses originating from the QCR are sufficient to realize the EP. The damping rate κ QCR shown in Fig. 4 is calculated using the measured electron temperatures of the normal-metal island (see Methods and Supplementary Fig. S7). At higher voltages, the major contribution in κ 2 is given by the residual damping coefficient, as discussed in Methods. This damping coefficient includes dephasing, quasiparticle losses, and resistive losses, and we extract its value based on the experimental rates. The dephasing has a similar effect on the measured coherent photon population as photon absorption due to other loss mechanisms although pure dephasing does not reduce the total photon number in the resonators.
In Fig. 4(b) we show the damping rates of the coupled circuit calculated from κ 2 using Eq. (2) as κ eff = −2Im(λ ± ). The maximum energy decay rate for the mode in R1 is obtained at the EP as discussed above, and it is given by Thus, in the optimal case, the decay rate is limited by the coupling strength between the resonators. DISCUSSION We have experimentally realized an exceptional point, EP, in a superconducting microwave circuit. We study the presence of the EP by observing a transition from an avoided crossing to single modulating resonance frequency. At this point, we achieve a maximum heat transfer between the two resonators without back and forth oscillation of the energy. The effective dissipation rate of the resonator R1 is bounded from above by the coupling strength in the optimal case, κ max,1 = 2g = 14 MHz. In the far detuned case with zero bias voltage, the effective rate of R1 is reduced down to κ min,1 = κ 1 = 260 kHz. Thus, the effective damping rate can be tuned by a factor of approximately 50 in our samples. Larger tuning can be obtained by minimizing the internal losses [66]. The measurement results are in very good agreement with our model. The circuit is based on a QCR, which enables the investigation of the crossover from an underdamped to critically damped and further to overdamped circuit. In addition to the realization of an EP, the circuit also behaves as a frequency-and voltage-tunable heat sink for quantum electric circuits that can be applied, for example, in quantum information processing for initializing qubits to their ground state by absorbing energy [43].
The tunability of the damping rate enables one to obtain the fastest possible photon absorption allowed by a given coupling coefficient.
In the future, it is interesting to further investigate the EP by modifying the circuit design. By introducing tunnel junctions to both resonators, one obtains a continuous line of EPs instead of an isolated singularity point. Incorporating qubits also enables the investigation and utilization of the EP with single energy quanta. Furthermore, one can investigate a circuit consisting of resonators with tunable damping realized with a QCR and a tunable coupling realized with a SQUID or a qubit [56,58]. The use of several microwave resonators will result in a more versatile parameter space [67], and hence yields an interesting platform for studying fundamental physics. Dynamic encircling of the EP with topological energy transfer [3,9] can be realized with superconducting resonators in a straightforward manner using standard microwave techniques. It requires fast tuning of the magnetic field, which can be realized by fabricating a flux bias line on the chip. Topological energy transfer with microwave pulses may provide an asset for applications in quantum information processing and other quantum technological devices. In addition, EPs are suitable for investigating PT symmetry on the level of single microwave photons. Here, superconducting circuits provide an attractive architecture owing to the ability to design system parameters yielding, for example, ultra-strong-and deep-strongcoupling regimes [21]. Furthermore, we see EPs as candidates to realize nonreciprocal signal routing beneficial for active quantum circuits [64].

Quantum-circuit refrigerator
We use a QCR to absorb and emit photons in the resonator R2. The resonator transition rate from the occupation number m to m can be written as [47] (5) where V = V b /2, R T is the tunneling resistance, M mm is the corresponding matrix element, R K = h/e 2 ≈ 25.8 kΩ is the von Klitzing constant, and the normalized rate for forward tunneling is given by The matrix element describing the transition can be written in terms of the generalized Laguerre polynomials L l n (ρ) as [47] where ρ = πα 2 /(ω 2 C l x 2 R K ) is a environmental parameter, where C l is the capacitance per unit length of the coplanar waveguide, 2x 2 is the length of the resonator R2, and the capacitance fraction α is given in terms of the capacitance between the normal-metal island and the center conductor C N , and junction capacitance C j as α = C N /(C N + 4C j ). In the equations above, we have neglected the effects owing to the charging of the normal-metal island since the capacitance of the island is relatively large. Furthermore, the rates for single-photon transitions can be expressed as [47] Γ m,m−1 = κ QCR (N + 1)m, where κ QCR denotes the coupling strength of the QCR, and the Bose-Einstein distribution at the effective temperature of the electron tunneling, T QCR , is given by where These equations are derived by defining Elastic tunneling in normal-metal-insulator-superconductor junctions Typically, the elastic tunneling is the dominating tunneling process. The electric current through a single NIS junction can be written as [50,68] (13) where T N denotes the normal-metal temperature, and V is the voltage across the junction. For a symmetric SINIS structure, we apply a voltage V b = 2V . Importantly, this equation has a monotonic dependence on the temperature of the normal metal but only a very weak dependence on the temperature of the superconductor. Thus, we may use NIS junctions as thermometers measuring the electron temperature of the normal metal.
The tunneling electrons transfer heat through the insulating barrier. The average power is given by [50] Based on this equation, we can reduce and increase the temperature of the normal metal. The applied voltage at the SINIS junction produces a total Joule heating power P = V b I, which is unequally divided between the N and S electrodes.

Quantum mechanical model
We analyze the temporal evolution of the coupled resonators following Ref. [60]. The Hamiltonian can be written in the rotating wave approximation aŝ The first term describes the energy of the primary resonator R1 with annihilation operatorâ, the second term the energy of R2 with annihilation operatorb, and the third term describes the coupling between the resonators. Here, we have neglected driving. Furthermore, this equation is valid only for a linear resonator. The effects owing to nonlinearity are discussed below. The dynamics of the system can be obtained from the Lindblad master equation for the density matrix of the coupled system,ρ, as where the Lindblad superoperator is given by L[x]ρ = xρx † − 1 2 {x †x ,ρ}. We can write the resulting equations of motion as [60] We define the resonator fields as â = A exp(−iω 1 t), b = B exp(−iω 1 t). Consequently, the equations assume the form where δ = ω 2 − ω 1 . These equations can be written in a matrix form as a time-dependent Schrödinger equation where ψ = (A, B) T , and as given in Eq. (1). Here, H is a non-Hermitian Hamiltonian scaled with . Equations (19) and (20) can also be written as a second-order differential equation, (23) When the resonators are tuned into resonance, δ = 0, we can express Eq. (23) as This equation describes a damped harmonic oscillator, where the energy is transferred between the resonators R1 and R2 at an angular frequency g 2 + κ 1 κ 2 /4. Due to the asymmetric damping rates in the two resonators, the total dissipation rate of the system is time-dependent and reaches its maximum value when the excitations are in R2. The damping ratio is given by Here, κ 2 is a function of voltage V b , which allows us to examine the transition from an underdamped system, ξ < 1, through critical damping, ξ = 1, to an overdamped system, ξ > 1. Critical damping is obtained when |κ 2 − κ 1 | = 4g. The total damping rate of R1 is given by κ 1 = κ int,1 + κ ext , and of R2 by , where κ int,1/2 denote the internal losses, κ ext the losses to the external measurement circuit, κ QCR (V b ) the photon-assisted tunneling in Eq. (12), and κ r,2 (V b ) the residual voltagedependent losses in R2. In our samples κ 2 (V b ) κ 1 , and g κ int,1 κ ext , as discussed below. Therefore, we obtain an approximate condition for the critical damping as The critical damping, which corresponds to the EP, is obtained at eV b /(2∆) ≈ 1 where the photon number remains low, and therefore, the slight nonlinearity caused by the SQUID is of negligible importance. However, at eV b /(2∆) > 1, the QCR generates thermal photons that result in photon-number-dependent losses, as discussed below.

Sample parameters
The main parameters for the samples are summarized in Supplementary Table S1. The coupling strength between the resonators can be estimated as [69] where the voltages are given by V i = ω 1 /(2x i C l ), i = 1, 2, the angular frequency of the second mode of the resonator R1 is ω 1 /(2π) ≈ 5.223 GHz for Samples A and B, and C l is the capacitance per unit length. Consequently, the critical damping is obtained with κ 2 = 4g ≈ 2π×29 MHz. The external quality factor corresponding to the leakage from the resonator R1 to the transmission line through the capacitances C TL is given by [70] Q ext = 2x 1 C l /(4Z L ω 1 C 2 TL ) ≈ 9 × 10 5 . Consequently, the corresponding damping rate is κ ext = ω 1 /Q ext ≈ 2π × 6 kHz. The loaded quality factor of the second mode of R1 is approximately Q L = 2 × 10 4 , when the resonators are far detuned. Thus, the internal losses in R1 dominate over the losses to the transmission line, Q int ≈ Q L . Furthermore, we obtain the damping rate κ 1 = ω 1 /Q L ≈ 2π × 300 kHz. The real part of the complex wave propagation coefficient, γ = α + iβ, describes the damping in the waveguide, and it can be calculated as [70] where n m is the mode number with n m = 2 denoting the first excited mode. The internal losses without the photon-assisted tunneling in the QCR are somewhat higher in the resonator R2 than in R1 since the design and fabrication of the QCR and the SQUID have not been optimized for low loss rates. The internal loss rate for R2 can be extracted at zero detuning and zero bias voltage from the saturation level of extracted κ 2 values since κ int,1 κ 2 (0), and hence the losses in R2 dominate over those in R1. We obtain from the circuit model the internal loss rate for R2 as κ int,2 ≈ κ 2 (0) = 2π × 16 MHz.
The photon number inside R1 when R2 is far detuned can be estimated as [71] n = 4Ω 2 d /κ 2 1 ≈ 10, where the driving strength is given by Ω d = C TL V in V 1 / , the input voltage is obtained from the input power as V in = √ P in Z L , and the input power is P in ≈ −115 dBm. The input power is −100 dBm for Sample A in Fig. 3 and Supplementary Fig. S4, and −115 dBm for Sample B in Supplementary Fig. S5. We also measure the resonators at different power levels. When R1 and R2 are in resonance at V b = 0, the total photon number is approximately equally divided between the resonators if κ 1 ≈ κ 2 . However, in our samples κ 1 < κ 2 , especially at V b > 2∆/e, and therefore the number of coherent photons is lower in R2 than in R1. When the Q factor of the resonator is reduced to 200, which is of the order of the critical damping, photon numbers close to unity are obtained with an input power P in ≈ −85 dBm.

Residual losses in the resonator R2
We attribute the residual voltage-dependent losses to dephasing, and to dissipation sources such as quasiparticle generation in the superconductors and resistive losses in the normal metal. Firstly, the resonator R2 is slightly nonlinear owing to the SQUID, and hence, an increasing incoherent photon number results in dephasing. Dephasing can be added in Eq. (16) with a term κ φ L[b †b ]ρ, where the dephasing rate κ φ depends on the number of thermal photons in the resonator. Similarly, in the case of superconducting qubits, the dephasing can be written as κ φ L[σ z ]ρ, whereσ z is a Pauli operator. The factor κ φ causes a similar effect as κ 2 in Eqs. (16)-(26) although it does not decrease the total photon number in the resonators. The photon number variance for a thermal state is of the form [38,39] n(n + 1), and therefore, thermal photons cause more dephasing than the coherent photons with a variance of n, where n is the average photon number. Consequently, we assume that κ φ = ω φ n(n + 1), where ω φ is a proportionality coefficient. Furthermore, as discussed above, the number of the coherent photons is low in R2 due to the relatively high loss rate. The steady-state photon number in the resonator can be estimated as [47] where we assume that the photon number of the effective bath, to which R2 is coupled through κ int,2 , vanishes owing to the very low cryostat temperatures of approximately 10 mK. The photon number depends linearly on the bias voltage at voltages above the superconductor energy gap, as shown in Supplementary Fig. S6.
Secondly, we take the quasiparticle losses into account. The critical temperature of Nb is approximately 9 K, and therefore, the quasiparticle density remains low in it. However, the critical temperature of Al approximately 1.2 K, which enables higher quasiparticle density than in Nb. We observe a decrease in the critical current of the SQUID, which indicates increased temperature in the Al leads of the SQUID, and hence heat dissipation. The quasiparticle loss rate [72,73] κ qp ∝ n qp ∝ √ P , where P is the absorbed power. The Al leads at the NIS junctions receive half of the Joule power P = IV b at high voltages, whereas the other half is absorbed to the normal metal. Thus, the power is quadratic in voltage, which is linear in the estimated photon number. Therefore, the expected quasiparticle losses are linear in photon number, κ qp = ω qp n, where ω qp is a proportionality coefficient. The dc power dissipated in the junctions is substantially higher than the microwave input power. At eV b /(2∆) = 2, the dc power is approximately 30 pW compared to a microwave power of −100 dBm = 0.1 pW. The normal metal in the QCR acts as an effective quasiparticle trap [73] minimizing the quasiparticle losses. Some fraction of the power dissipated at the QCR leaks to the SQUID.
There is an approximately 10-µm-long section of normal metal between the actual Nb resonator and the NIS junctions, which may cause some losses. The loss rate at the resistor depends on the current profile of the microwave mode, which can depend on the voltage V b . Nevertheless, we assume these losses to be small due to the QCR being at the end of the resonator. Furthermore, there is a layer of superconducting Al below the normal metal due to the shadow evaporation technique, which decreases the current in the resistor, and hence also the resistive losses. The very weak resistive losses are quadratic in the voltage amplitude of the microwave resonator which is linear in photon number. Thus, it can be approximated as κ res = ω res n with a proportionality coefficient ω res .
Consequently, the total voltage-dependent losses in R2 including the dephasing, quasiparticle losses in the superconductors and the resistive losses are given by κ r,2 = κ φ + κ qp + κ res = κ φ n(n + 1) + κ qp n + κ res n. (28) The quasiparticle and resistive losses are expected to be very weak, as described above, but a small contribution cannot be excluded. Nevertheless, we expect the dephasing to dominate over the quasiparticle and resistive losses. Therefore, in the numerical analysis, we take the photonnumber-dependent losses into account as κ r,2 = ω r,tot n(n + 1), (29) with only one fitting parameter ω r,tot effectively describing the different loss methods discussed above. From the experimental damping rates of the dissipative resonator R2, we extract the the coefficient ω r,tot ≈ 2π × 22 MHz. The good agreement with the experimental damping rate κ 2 and the model with the quadratic residual losses κ r,2 in Fig. 4(a) gives further support for the approximation in Eq. (29). We do not take this loss rate into account in Eq. (27) for simplicity, and also due to the fact that pure dephasing does not decrease the photon number. The odd modes of R1 do not show flux dependence as expected due to the voltage node at the coupling capacitor. However, they do show some dependence on the voltage V b . Similar dependence can be observed also for the even modes at Φ/Φ 0 = 0.5 where the inductance of the SQUID ideally vanishes and thus decouples the QCR from the resonator R1. We attribute this observation to unintentional asymmetry of the sample. Furthermore, the QCR may be weakly coupled to the input and output microwave fields through some spurious mode of the sample holder. The very broad resonance at high bias voltages enables the coupling to the spurious modes. We note that the spurious modes may be partially responsible for the κ r,2 . However, we do not quantitatively model these losses. Instead, they are effectively included in the parameter ω r,tot in Eq. (29).

Full model for κ2 and κ eff
The parameters κ 2 and κ eff are obtained as follows. First, we extract the effective resistance corresponding to the QCR by fitting the classical circuit model to the experimentally obtained scattering parameter S 21 using a least-squares algorithm. Second, we calculate the quality factor of the resonator R2, Q R2 , for the obtained effective resistance, as discussed below. The coupling rate is related to the quality factor as κ 2 = ω 2 /Q R2 . The full model denoted by the line in Fig. 4(a) is obtained by fitting to the experimental transition rates according to Eqs. (12), (27), and (29). Here, we use ω r,tot as the only fitting parameter since we fix κ int,2 to the saturation value at zero bias, as discussed above. Subsequently, we may proceed to the effective damping rates κ eff = −2Im(λ ± ), which can be obtained from κ 2 with the help of Eq. (2). The damping rates above the critical damping, κ 2 > 4g, result in the two branches at bias voltages V b 2∆/e. The lines in Fig. 4(b) are obtained similarly from Eq. (2).

Classical circuit model
To simulate the scattering parameter S 21 , we use a classical circuit model similar to the one presented in Ref. [61]. We analyze the samples using standard microwave circuit analysis [74]. The input impedance of the resonator R2 is , (31) where the impedance of the SQUID and the capacitors between the SQUID and the center conductor is given by Z S = iωL S + 2/(iωC S ), the impedance of the coupling capacitor between the resonators by Z C = 1/(iωC C ), γ is the complex propagation coefficient discussed above, and the terminating impedance consisting of the effective resistance of the NIS junctions and the capacitor between the normal-metal island and the center conductor is modeled as an effective resistor with resistance R eff . The inductance of the SQUID is calculated as where the maximum supercurrent through the SQUID is I 0 , and the flux quantum is Φ 0 = h/(2e).
The scattering parameter S 21 describing the voltage transmission from Port 1 to Port 2 can be calculated using the transmission matrix method as [74] where Z L is the characteristic impedance of the external measurement cables, and with We analyze the losses in the resonator R2 also in the absence of coupling to R1. In particular, we omit matrices M 1 and M 2 from Eq. (33). The resonator R2 causes a dip in the amplitude of the transmission coefficient S 21 , whereas R1 causes a peak. The quality factor can be estimated directly from the ratio of the center frequency and the width of the peak or dip. Alternatively, more advanced methods can be used [75].

Sample fabrication
The samples are fabricated on a Si wafer with a thickness of 500 µm and a diameter of 100 mm. First, a 300nm-thick layer of SiO 2 is thermally grown on the wafer with resistivity ρ > 10 kΩ cm. Subsequently, a 200-nmthick layer of Nb is sputtered on top of the oxide. The resonators are patterned on the Nb layer with optical lithography and reactive ion etching. We cover the complete wafer with a 40-nm-thick layer of Al 2 O 3 fabricated using atomic layer deposition. This oxide layer serves as an insulating barrier in the parallel plate capacitors and separates the QCR lines from the ground plane. The nanostructures are defined using electron beam lithography and two-angle shadow evaporation followed by a lift-off process. The SQUID consists of two Al layers with thicknesses of 40 nm each. The first Al layer is oxidized in situ in the evaporation chamber at 1.0 mbar for 5 min. The SINIS junctions consist of Al (40 nm) and Cu (40 nm), and the Al layer is similarly oxidized as in the SQUID. The shadow evaporation technique results in overlapping metal layers.

Measurement setup
The measurement setup is schematically presented in Fig. S3. The samples are measured in a commercial drydilution refrigerator with a base temperature of approximately 10 mK. The scattering parameters are measured with a vector network analyzer which contains both the microwave source and the detector. The microwave signal is attenuated at different temperature stages to avoid heat leakage from higher temperatures to the sample. We employ amplifiers at 4 K and at room temperature. The NIS junctions are controlled by applying a bias voltage or current through continuous thermocoax cables from room temperature down to the base temperature. Magnetic flux for the SQUID is produced using a superconducting coil with a bias current.

Normalization of scattering parameters
All measured scattering parameters S 21 are normalized. Initially, we normalize the phase winding originating from the electrical delay τ ≈ 50 ns in the measurement setup outside the sample by multiplying with exp(iωτ ). Consequently, the resonance produces a circle on the complex plane as the frequency is increased over the resonance. We transform this circle to its canonical position where max |S 21 | is on the positive real axis and the circle intersects the origin. Finally, we normalize the amplitude to unity by dividing with max |S 21 |.
Supplementary Table S1. Parameters. The parameters for Sample B that differ from those for Sample A are given in parenthesis. See Methods and Supplementary Fig. S3(b) for details. The resonance frequency of the first excited mode of the resonator R1, f1, is a measured value, the characteristic impedances of the transmission lines in the resonator, Z0, and in the external measurement circuit, ZL, are nominal values. The lengths of the resonator sections x1 and x2 are design values, and the effective resistance is calculated as [70] √ ε eff = c/(2f1x1), where c is the speed of light in vacuum. We obtain the values for the capacitance per unit length C l , and the capacitance CTL from finite element method (FEM) simulations. The capacitance CC is obtained by fitting the circuit model to the measured scattering parameter |S21| in good agreement with FEM simulations, and the capacitances CN, CS and Cj are calculated using a parallel-plate model. The coupling strength g is obtained from CC. The loaded quality factor of the first exited mode of R1 Qint,1 and the tunneling resistance RT are measured values, and the Dynes parameter γD is estimated as the ratio of the asymptotic resistance and the resistance in the superconductor gap. The critical current at zero-bias Ic,0 is given by the flux corresponding to the crossing of the modes in the circuit model in good agreement with a control sample with slightly smaller junction area and a critical current of approximately 200 nA. The damping rate κ1 is given by the ratio ω1/Qint,1, and the damping rate κint,2 is extracted from the saturation value of κ2 at zero bias. The proportionality coefficient for the residual losses ωr,tot is a fitted value.
Parameter      The theoretical curve is calculated using Eq. (13). We extract the electron temperatures of the normal-metal island using a linear voltage-to-temperature conversion below 300 mK, and above that we extract the temperatures from the voltages corresponding to the different experimental bath temperature points. At high temperatures, the low sensitivity reduces the reliability of the extracted island temperatures. (c) Electron temperature of the normal-metal island as a function of bias voltage at different bath temperatures.