Optimized heat transfer at exceptional points in quantum circuits

TL;DR

This paper demonstrates the experimental realization of an exceptional point in superconducting microwave circuits, enabling optimized heat transfer and tunable dissipation, with implications for quantum computing and PT-symmetric systems.

Contribution

It introduces the first experimental realization of an exceptional point in superconducting circuits, showing tunable heat flow and dissipation control at the quantum level.

Findings

Observation of a crossover from underdamped to overdamped coupling at the EP

Demonstration of in situ tunable dissipation using photon-assisted tunneling

Potential applications in qubit initialization and PT-symmetric system studies

Abstract

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.

Tables

Table 1. 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, f 1 subscript 𝑓 1 f_{1} , is a measured value, the characteristic impedances of the transmission lines in the resonator, Z 0 subscript 𝑍 0 Z_{0} , and in the external measurement circuit, Z L subscript 𝑍 L Z_{\textrm{L}} , are nominal values. The lengths of the resonator sections x 1 subscript 𝑥 1 x_{1} and x 2 subscript 𝑥 2 x_{2} are design values, and the effective resistance is calculated as  Göppl  et al. ( 2008 ) ε eff = c / ( 2 ​ f 1 ​ x 1 ) subscript 𝜀 eff 𝑐 2 subscript 𝑓 1 subscript 𝑥 1 \sqrt{\varepsilon_{\textrm{eff}}}=c/(2f_{1}x_{1}) , where c 𝑐 c is the speed of light in vacuum. We obtain the values for the capacitance per unit length C l subscript 𝐶 l C_{\textrm{l}} , and the capacitance C TL subscript 𝐶 TL C_{\textrm{TL}} from finite element method (FEM) simulations. The capacitance C C subscript 𝐶 C C_{\textrm{C}} is obtained by fitting the circuit model to the measured scattering parameter | S 21 | subscript 𝑆 21 |S_{21}| in good agreement with FEM simulations, and the capacitances C N subscript 𝐶 N C_{\textrm{N}} , C S subscript 𝐶 S C_{\textrm{S}} and C j subscript 𝐶 j C_{\textrm{j}} are calculated using a parallel-plate model. The coupling strength g 𝑔 g is obtained from C C subscript 𝐶 C C_{\textrm{C}} . The loaded quality factor of the first exited mode of R1 Q int,1 subscript 𝑄 int,1 Q_{\textrm{int,1}} and the tunneling resistance R T subscript 𝑅 T R_{\textrm{T}} are measured values, and the Dynes parameter γ D subscript 𝛾 D \gamma_{\textrm{D}} is estimated as the ratio of the asymptotic resistance and the resistance in the superconductor gap. The critical current at zero-bias I c,0 subscript 𝐼 c,0 I_{\textrm{c,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 subscript 𝜅 1 \kappa_{1} is given by the ratio ω 1 / Q int,1 subscript 𝜔 1 subscript 𝑄 int,1 \omega_{1}/Q_{\textrm{int,1}} , and the damping rate κ int,2 subscript 𝜅 int,2 \kappa_{\textrm{int,2}} is extracted from the saturation value of κ 2 subscript 𝜅 2 \kappa_{2} at zero bias. The proportionality coefficient for the residual losses ω r,tot subscript 𝜔 r,tot \omega_{\textrm{r,tot}} is a fitted value.

Parameter	Value
$f_1$	5.223 GHz
$Z_0$	50 $\Omega$
$Z_{\text{L}}$	50 $\Omega$
$x_1$	12 mm
$x_2$	6.0 (6.5) mm
${\epsilon }_{\text{eff}}$	5.73
$C_{\text{l}}$	155 pF/m
$C_{\text{TL}}$	0.8 fF
$C_{\text{C}}$	3.8 fF
$C_{\text{N}}$	98 fF
$C_{\text{S}}$	460 fF
$C_{\text{j}}$	6.2 fF
$g/\left(2\pi \right)$	7.2 MHz
$Q_{\text{int,1}}$	$2.7\times {10}^4$ $\left(2.0\times {10}^4\right)$
$R_{\text{T}}$	8.4 (9.5) k$\Omega$
${\gamma }_{\text{D}}$	$1\times {10}^{-4}$
$I_{\text{c,0}}$	340 (300) nA
${\kappa }_{\text{1}}/\left(2\pi \right)$	190 (260) kHz
${\kappa }_{\text{int,2}}/\left(2\pi \right)$	16 MHz
${\omega }_{\text{r,tot}}/\left(2\pi \right)$	22 MHz

Equations

$$H=\left(\begin{array}{cc}-i\frac{{\kappa }_1}{2}& g\\ g& \delta -i\frac{{\kappa }_2}{2}\end{array}\right),$$

$${\lambda }_{\pm }=\frac{1}{4}\left(2\delta -i{\kappa }_1-i{\kappa }_2\pm s\right),$$

$${\psi }_{\pm }={\left(\frac{-2\delta -i{\kappa }_1+i{\kappa }_2\pm s}{4g},1\right)}^{\text{T}},$$

$$s=\sqrt{4{\delta }^2+16g^2-({\kappa }_2-{\kappa }_1{)}^2-i4\delta ({\kappa }_2-{\kappa }_1)}.$$

$${\Gamma }_{m,m^{\prime }}(V)=M_{m{m}^{\prime }}^2\frac{2R_{\text{K}}}{R_{\text{T}}}\sum _{\tau =\pm 1}\overrightarrow{F}[\tau eV+\hslash {\omega }_2(m-m^{\prime })],$$

$$\overrightarrow{F}(E)=\frac{1}{h}{\int }_{-\infty }^{\infty }\mathrm{d}{E}^{\prime }{n}_{\mathrm{S}}(E^{\prime })[1-f(E^{\prime },T_{\text{S}})]f(E^{\prime }-E,T_{\text{N}}),$$

$$n_{\text{S}}(E)=\left|\text{Re}\left(\frac{E/\Delta +i{\gamma }_{\text{D}}}{\sqrt{(E/\Delta +i{\gamma }_{\text{D}}{)}^2-1}}\right)\right|.$$

$$M_{m,m^{\prime }}^2=\{\begin{array}{ll}{e}^{-\rho }{\rho }^{m-m^{\prime }}\frac{m^{\prime }!}{m!}[L_{m^{\prime }}^{m-m^{\prime }}(\rho ){]}^2,& m\ge m^{\prime },\\ e^{-\rho }{\rho }^{m^{\prime }-m}\frac{m!}{m^{\prime }!}[L_{m^{\prime }}^{m^{\prime }-m}(\rho ){]}^2,& m<m^{\prime },\end{array}$$

$$\Gamma {}_{m,m-1}$$

$$\Gamma {}_{m,m+1}$$

$$N_{\text{QCR}}=\frac{1}{\exp \left(\frac{\hslash {\omega }_2}{k_{\text{B}}{T}_{\text{QCR}}}\right)-1},$$

$$T_{\text{QCR}}=\frac{\hslash {\omega }_2}{k_{\text{B}}}{\left\{\ln \left[\frac{{\sum }_{\tau =\pm 1}\overrightarrow{F}(\tau eV+\hslash {\omega }_2)}{{\sum }_{\tau =\pm 1}\overrightarrow{F}(\tau eV-\hslash {\omega }_2)}\right]\right\}}^{-1}.$$

$${\kappa }_{\text{QCR}}=\frac{\Gamma {}_{m,m-1}}{m}-\frac{{\Gamma }_{m,m+1}}{m+1}.$$

$$I(V)=\frac{1}{e{R}_{\text{T}}}{\int }_0^{\infty }{n}_{\text{S}}(E)[f(E-eV,T_{\text{N}})-f(E+eV,T_{\text{N}})]\text{d}E,$$

$$P=\frac{1}{e^2{R}_{\text{T}}}{\int }_{-\infty }^{\infty }{n}_{\text{S}}(E)(E-eV)[f(E-eV,T_{\text{N}})-f(E,T_{\text{S}})]\text{d}E.$$

$${\hat{H}}_{\text{RWA}}=\hslash {\omega }_1{\hat{a}}^{†}\hat{a}+\hslash {\omega }_2{\hat{b}}^{†}\hat{b}+\hslash g(\hat{a}{\hat{b}}^{†}+{\hat{a}}^{†}\hat{b}).$$

$$\frac{\text{d}\hat{\rho }}{\text{d}t}=-\frac{i}{\hslash }[{\hat{H}}_{\text{RWA}},\hat{\rho }]+{\kappa }_1\mathcal{L}[\hat{a}]\hat{\rho }+{\kappa }_2\mathcal{L}[\hat{b}]\hat{\rho },$$

$$\frac{\text{d}⟨\hat{a}⟩}{\text{d}t}$$

$$\frac{\text{d}⟨\hat{b}⟩}{\text{d}t}$$

$$\frac{\text{d}A}{\text{d}t}$$

$$\frac{\text{d}B}{\text{d}t}$$

$$\frac{\text{d}}{\text{d}t}\psi =-iH\psi ,$$

$$H=\left(\begin{array}{cc}-i\frac{{\kappa }_1}{2}& g\\ g& -i\frac{{\kappa }_2}{2}+\delta \end{array}\right),$$

$$\frac{{\text{d}}^{2}A}{\text{d}{t}^2}+\left(\frac{{\kappa }_1+{\kappa }_2}{2}+i\delta \right)\frac{\text{d}A}{\text{d}t}+\left(g^2+i\delta \frac{{\kappa }_1}{2}+\frac{{\kappa }_1{\kappa }_2}{4}\right)A=0.$$

$$\frac{{\text{d}}^{2}A}{\text{d}{t}^2}+\frac{{\kappa }_1+{\kappa }_2}{2}\frac{\text{d}A}{\text{d}t}+\left(g^2+\frac{{\kappa }_1{\kappa }_2}{4}\right)A=0.$$

$$\xi =\frac{{\kappa }_1+{\kappa }_2}{2\sqrt{4g^2+{\kappa }_1{\kappa }_2}}.$$

$${\kappa }_2=4g.$$

$$n=\frac{{\kappa }_{\text{QCR}}{N}_{\text{QCR}}}{{\kappa }_{\text{QCR}}+{\kappa }_{\text{int,2}}},$$

$${\kappa }_{\text{r,2}}={\kappa }_{\varphi }+{\kappa }_{\text{qp}}+{\kappa }_{\text{res}}={\kappa }_{\varphi }n(n+1)+{\kappa }_{\text{qp}}n+{\kappa }_{\text{res}}n.$$

$${\kappa }_{\text{r,2}}={\omega }_{\text{r,tot}}n(n+1),$$

$${\kappa }_2(V_{\text{b}})={\kappa }_{\text{QCR}}(V_{\text{b}})+{\kappa }_{\text{r,2}}(V_{\text{b}})+{\kappa }_{\text{int,2}}$$

$$Z_{\text{R2}}=Z_{\text{C}}+\frac{Z_0\left\{Z_{\text{S}}+Z_0\tanh (\gamma x_2)+\frac{Z_0\left[R_{\text{eff}}+Z_0\tanh (\gamma x_2)\right]}{Z_0+R_{\text{eff}}\tanh (\gamma x_2)}\right\}}{Z_0+\tanh (\gamma x_2)\left\{Z_{\text{S}}+\frac{Z_0[R_{\text{eff}}+Z_0\tanh (\gamma x_2)]}{Z_0+R_{\text{eff}}\tanh (\gamma x_2)}\right\}},$$

Full text

Systems with effective non-Hermitian Hamiltonians have been actively studied in various setups in recent years Milburn et al. (2015); Xu et al. (2016); Doppler et al. (2016); Mandal and Tewari (2016); Ding et al. (2016, 2018); San-Jose et al. (2016). They show many intriguing properties such as singularities in their energy spectra Kato (1966); Dembowski et al. (2001); Heiss (2004); Berry (2004); Heiss (2012); Cao and Wiersig (2015). A square-root singularity point in the parameter space of a non-Hermitian matrix is called an exceptional point, EP, if the eigenvalues coalesce Heiss (2012); Cao and Wiersig (2015). Previously, EPs have been shown to emerge, for example, in non-superconducting microwave circuits, laser physics, quantum phase transitions, and atomic and molecular physics Heiss (2012); Cao and Wiersig (2015). The fascinating effects of EPs include the disappearance of the beating Rabi oscillations Dietz et al. (2007), chiral states in microwave systems Dembowski et al. (2003), and spontaneous symmetry breaking in systems with parity- and time-reversal ($\mathcal{PT}$) symmetry Bender and Boettcher (1998); Guo et al. (2009); Chtchelkatchev et al. (2012); El-Ganainy et al. (2018). In the quantum regime, $\mathcal{PT}$-symmetric systems may show features that are different from the semiclassical predictions, such as new phases owing to quantum fluctuations Kepesidis et al. (2016); El-Ganainy et al. (2018). Despite the active research on EPs, they have not been thoroughly investigated in superconducting microwave circuits to date Quijandría et al. (2018).

Superconducting microwave circuits provide an ideal platform to realize various quantum devices, such as ultrasensitive photon detectors and counters Inomata et al. (2016); Govenius et al. (2016); Ding et al. (2017); Opremcak et al. (2018), and potentially even a large-scale quantum computer Ladd et al. (2010); Clarke and Wilhelm (2008), or a quantum simulator Georgescu et al. (2014) in the framework of circuit quantum electrodynamics Blais et al. (2004); Wallraff et al. (2004). Notably, superconducting qubits have been shown to approach the required coherence times Barends et al. (2014); Kelly et al. (2015) for quantum error correction Lidar and Brun (2013); Terhal (2015). However, despite the tremendous interest in superconducting microwave circuits in recent years You and Nori (2011); Devoret and Schoelkopf (2013); Wendin (2017), 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 Clerk and Utami (2007); Goetz et al. (2017). In many error correction codes, qubits are repeatedly initialized, which requires fast and efficient cooling schemes Geerlings et al. (2013); Bultink et al. (2016); Wong et al. (2018). One promising method for absorbing energy and initializing qubits to their ground state is based on resonators with tunable dissipation Tuorila et al. (2017).

The recently developed quantum-circuit refrigerator (QCR) Tan et al. (2017); Masuda et al. (2018) 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 Silveri et al. (2018). Operation of the QCR relies on inelastic tunneling of electrons through a normal-metal–insulator–superconductor (NIS) junction Silveri et al. (2017). 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 Silveri et al. (2018). Furthermore, elastic tunneling can be utilized for temperature control of the normal-metal electrons Nahum et al. (1994); Leivo et al. (1996); Giazotto et al. (2006), and for precise thermometry down to millikelvin temperatures Feshchenko et al. (2015); Giazotto et al. (2006). Recently, NIS junctions have also been utilized in a realization of a quantum heat valve Ronzani et al. (2018), and phase-coherent caloritronics Fornieri and Giazotto (2017).

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 flux-tunable resonance frequency (Fig. 1). We denote the NIS junctions and the normal-metal island that is capacitively coupled to the resonator as a QCR. Thanks to the voltage-tunable dissipation within the QCR and the flux-tunable 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 Wang et al. (2013); Pierre et al. (2014); Baust et al. (2015); Vissers et al. (2015); Wulschner et al. (2016); Adamyan et al. (2016); Pierre et al. (2018); Wong et al. (2018); Partanen et al. (2018) 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 Dickson et al. (2013) and to studies of the eigenstate thermalization hypothesis in many-body quantum problems Nandkishore and Huse (2015). Furthermore, our work is an important step towards the investigation of $\mathcal{PT}$-symmetric systems at the quantum level that can be realized with circuit quantum electrodynamics architectures Metelmann and Türeci (2018); Quijandría et al. (2018).

Results

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\text{\,}\mathrm{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_{\textrm{C}}$ is located. Therefore, we focus on the first excited mode of R1, with frequency $f_{1}=\omega_{1}/(2\pi)=$5.2\text{,}\mathrm{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\pi)=$7.2\text{\,}\mathrm{MHz}. In contrast to R1, the resonator R2 has a flux-tunable resonance frequency $\omega_{2}(\Phi)$ owing to a superconducting quantum interference device (SQUID), and a voltage tunable loss rate $\kappa_{2}(V_{\textrm{b}})$ owing to the QCR. Here, $\Phi$ and $V_{\textrm{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 $\Phi_{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\Delta$ in the density of states of the superconductor, and the difference can be utilized to cool down quantum circuits Tan et al. (2017); Silveri et al. (2017). The sample fabrication is described in Methods. We study two samples with different R2 resonator lengths, Sample A ($12\text{\,}\mathrm{m}\mathrm{m}$) and Sample B ($13\text{\,}\mathrm{m}\mathrm{m}$). The R1 resonator has a length of $24\text{\,}\mathrm{mm}$ in both samples. Sample parameters are summarized in 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 $\kappa_{2}(V_{\textrm{b}})$ and the flux-tunable detuning between the resonators $\delta=\omega_{2}(\Phi)-\omega_{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. Doppler et al. (2016). Thus, the excitations of the system can be described with the effective non-Hermitian Hamiltonian in matrix form in the basis $\psi=(A,B)^{\textrm{T}}$ where $A$ and $B$ are field amplitudes in R1 and R2, respectively, (see Methods and Ref. Pierre et al. (2018))

[TABLE]

where $\kappa_{1}$ is the decay rate of the resonator R1. The eigenvalues of $H$ can be written as

[TABLE]

and the corresponding eigenvectors are

[TABLE]

where

[TABLE]

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 $|\kappa_{2}-\kappa_{1}|=4g$, and $\delta=0$. In the following, we assume that $\kappa_{1}\ll\kappa_{2}$, which is valid for our samples, as discussed below. Consequently, the condition for the EP simplifies to $\kappa_{2}=4g$.

To visualize the system singularity, i.e., the EP, we show the real and imaginary parts of $\lambda_{\pm}$ in Fig. 2. The eigenvalues form a self-intersecting Riemann surface in the parameter space of $\kappa_{2}$ and $\delta$. 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, $\kappa_{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, $\kappa_{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 Pierre et al. (2018); Wong et al. (2018). In particular, the rate of heat transfer at zero detuning is given by $\kappa_{\textrm{eff}}=-2\textrm{Im}(\lambda_{\pm})\approx\kappa_{2}[1\mp\textrm{Re}(\sqrt{1-(4g/\kappa_{2})^{2}})]/2$. 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 $\kappa_{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 frequency-dependent 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_{\textrm{eff}}$ (see Methods and Supplementary Fig. S3). We extract $R_{\textrm{eff}}$ by fitting the circuit model to the experimental results, and use $R_{\textrm{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_{\textrm{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\text{\,}\mathrm{fF}$ which agrees well with the finite-element-method simulation that yields approximately $5\text{\,}\mathrm{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 Silveri et al. (2018), which, however, 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 $\kappa_{1}\ll\kappa_{2}$, we measure the internal quality factor of the primary resonator R1 with R2 far detuned at $\Phi=\Phi_{0}/2$. From the quality factor, we extract a loss rate $\kappa_{1}/(2\pi)\lesssim$260\text{,}\mathrm{kHz}$$ for both samples which is substantially lower than the extracted value of $\kappa_{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 Zmuidzinas (2012) (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 $\kappa_{1}$ and $\kappa_{2}$. The value of $\kappa_{2}$ can be tuned by approximately two orders of magnitude. At low bias voltages $eV_{\textrm{b}}/(2\Delta)<1$, the rate $\kappa_{2}$ is below $4g$, and at $eV_{\textrm{b}}/(2\Delta)>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 $\kappa_{\textrm{QCR}}$, as well as constant internal losses $\kappa_{\textrm{int,2}}$, and voltage-dependent residual losses $\kappa_{\textrm{r,2}}$. We have designed our sample in such a way that $\kappa_{\textrm{QCR}}$ covers the critical damping rate, and thus the losses originating from the QCR are sufficient to realize the EP. The damping rate $\kappa_{\textrm{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 $\kappa_{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 $\kappa_{2}$ using Eq. (2) as $\kappa_{\textrm{eff}}=-2\textrm{Im}(\lambda_{\pm})$. The maximum energy decay rate for the mode in R1 is obtained at the EP as discussed above, and it is given by $\kappa_{\textrm{max,1}}=\max[-2\textrm{Im}(\lambda_{+})]\approx\kappa_{2}/2=2g$. 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, $\kappa_{\textrm{max,1}}=2g=$14\text{,}\mathrm{MHz}. In the far detuned case with zero bias voltage, the effective rate of R1 is reduced down to $\kappa_{\textrm{min,1}}=\kappa_{1}=$260\text{\,}\mathrm{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 Megrant et al. (2012). 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 Tuorila et al. (2017). 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 Baust et al. (2015); Wulschner et al. (2016). The use of several microwave resonators will result in a more versatile parameter space Demange and Graefe (2012), and hence yields an interesting platform for studying fundamental physics. Dynamic encircling of the EP with topological energy transfer Dembowski et al. (2001); Doppler et al. (2016) 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 $\mathcal{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-strong-coupling regimes Quijandría et al. (2018). Furthermore, we see EPs as candidates to realize nonreciprocal signal routing beneficial for active quantum circuits Metelmann and Türeci (2018).

Acknowledgements

We thank A. A. Clerk for discussions, and J. Govenius and M. Jenei for assistance. We acknowledge the provision of facilities and technical support by Aalto University at OtaNano - Micronova Nanofabrication Centre. We acknowledge the funding from the European Research Council under Consolidator Grant No. 681311 (QUESS), and Marie Skłodowska-Curie Grant No. 795159, the Academy of Finland through its Centres of Excellence Program (project Nos. 312300, 312059) and grants (Nos. 265675, 305237, 305306, 308161, 312300, 314302, 316551), the Vilho, Yrjö and Kalle Väisälä Foundation, the Technology Industries of Finland Centennial Foundation, the Jane and Aatos Erkko Foundation the Alfred Kordelin Foundation, and the Emil Aaltonen Foundation.

Author contributions

M.P. was responsible for sample design, fabrication, measurements, data analysis, and writing the manuscript. K.Y.T, and J.G. contributed to the sample design. L.G. deposited the Nb layer. J.G., V.S., K.K., and K.Y.T. contributed to measurements and data analysis. M.P., K.Y.T, and R.E.L. developed the fabrication process. V.V., R.K., J.I., and D.H. contributed to the measurements. M.P.S., A.M., and E.H. contributed to the theory. M.M. supervised the project. All authors commented on the manuscript.

Methods

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^{\prime}$ can be written as Silveri et al. (2017)

[TABLE]

where $V=V_{\textrm{b}}/2$, $R_{\textrm{T}}$ is the tunneling resistance, $M_{mm^{\prime}}$ is the corresponding matrix element, $R_{\textrm{K}}=h/e^{2}\approx 25.8$ k$\Omega$ is the von Klitzing constant, and the normalized rate for forward tunneling is given by

[TABLE]

where $f(E,T)=1/\{\exp[E/(k_{\textrm{B}}T)]+1\}$ is the Fermi–Dirac distribution, $k_{\textrm{B}}$ is the Boltzmann constant, and the density of states in a superconductor can be expressed with the help of the Dynes parameter $\gamma_{\textrm{D}}$ as

[TABLE]

The matrix element describing the transition can be written in terms of the generalized Laguerre polynomials $L_{n}^{l}(\rho)$ as Silveri et al. (2017)

[TABLE]

where $\rho=\pi\alpha^{2}/(\omega_{2}C_{\textrm{l}}x_{2}R_{\textrm{K}})$ is a environmental parameter, where $C_{\textrm{l}}$ is the capacitance per unit length of the coplanar waveguide, $2x_{2}$ is the length of the resonator R2, and the capacitance fraction $\alpha$ is given in terms of the capacitance between the normal-metal island and the center conductor $C_{\textrm{N}}$, and junction capacitance $C_{\textrm{j}}$ as $\alpha=C_{\textrm{N}}/(C_{\textrm{N}}+4C_{\textrm{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 Silveri et al. (2017)

[TABLE]

where $\kappa_{\textrm{QCR}}$ denotes the coupling strength of the QCR, and the Bose–Einstein distribution at the effective temperature of the electron tunneling, $T_{\textrm{QCR}}$, is given by

[TABLE]

where

[TABLE]

These equations are derived by defining

[TABLE]

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 Giazotto et al. (2006); Pekola et al. (2010)

[TABLE]

where $T_{\textrm{N}}$ denotes the normal-metal temperature, and $V$ is the voltage across the junction. For a symmetric SINIS structure, we apply a voltage $V_{\textrm{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 Giazotto et al. (2006)

[TABLE]

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_{\textrm{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. Pierre et al. (2018). The Hamiltonian can be written in the rotating wave approximation as

[TABLE]

The first term describes the energy of the primary resonator R1 with annihilation operator $\hat{a}$, the second term the energy of R2 with annihilation operator $\hat{b}$, 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, $\hat{\rho}$, as

[TABLE]

where the Lindblad superoperator is given by $\mathcal{L}[\hat{x}]\hat{\rho}=\hat{x}\hat{\rho}\hat{x}^{\dagger}-\frac{1}{2}\{\hat{x}^{\dagger}\hat{x},\hat{\rho}\}$. We can write the resulting equations of motion as Pierre et al. (2018)

[TABLE]

We define the resonator fields as $\langle\hat{a}\rangle=A\exp(-i\omega_{1}t)$, $\langle\hat{b}\rangle=B\exp(-i\omega_{1}t)$. Consequently, the equations assume the form

[TABLE]

where $\delta=\omega_{2}-\omega_{1}$. These equations can be written in a matrix form as a time-dependent Schrödinger equation

[TABLE]

where $\psi=(A,B)^{\textrm{T}}$, and

[TABLE]

as given in Eq. (1). Here, $H$ is a non-Hermitian Hamiltonian scaled with $\hbar$. Equations (19) and (20) can also be written as a second-order differential equation,

[TABLE]

When the resonators are tuned into resonance, $\delta=0$, we can express Eq. (23) as

[TABLE]

This equation describes a damped harmonic oscillator, where the energy is transferred between the resonators R1 and R2 at an angular frequency $\sqrt{g^{2}+\kappa_{1}\kappa_{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

[TABLE]

Here, $\kappa_{2}$ is a function of voltage $V_{\textrm{b}}$, which allows us to examine the transition from an underdamped system, $\xi<1$, through critical damping, $\xi=1$, to an overdamped system, $\xi>1$. Critical damping is obtained when $|\kappa_{2}-\kappa_{1}|=4g$. The total damping rate of R1 is given by $\kappa_{1}=\kappa_{\textrm{int,1}}+\kappa_{\textrm{ext}}$, and of R2 by $\kappa_{2}(V_{\textrm{b}})=\kappa_{\textrm{int,2}}+\kappa_{\textrm{QCR}}(V_{\textrm{b}})+\kappa_{\textrm{r,2}}(V_{\textrm{b}})$, where $\kappa_{\textrm{int,1/2}}$ denote the internal losses, $\kappa_{\textrm{ext}}$ the losses to the external measurement circuit, $\kappa_{\textrm{QCR}}(V_{\textrm{b}})$ the photon-assisted tunneling in Eq. (12), and $\kappa_{\textrm{r,2}}(V_{\textrm{b}})$ the residual voltage-dependent losses in R2. In our samples $\kappa_{\textrm{2}}(V_{\textrm{b}})\gg\kappa_{\textrm{1}}$, and $g\gg\kappa_{\textrm{int,1}}\gg\kappa_{\textrm{ext}}$, as discussed below. Therefore, we obtain an approximate condition for the critical damping as

[TABLE]

The critical damping, which corresponds to the EP, is obtained at $eV_{\textrm{b}}/(2\Delta)\approx 1$ where the photon number remains low, and therefore, the slight nonlinearity caused by the SQUID is of negligible importance. However, at $eV_{\textrm{b}}/(2\Delta)>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 Jones et al. (2013) $g=C_{\textrm{C}}V_{1}V_{2}/\hbar\approx 2\pi\times$7.2\text{,}\mathrm{MHz}, where the voltages are given by $V_{i}=\sqrt{\hbar\omega_{1}/(2x_{i}C_{\textrm{l}})}$, $i=1,2$, the angular frequency of the second mode of the resonator R1 is $\omega_{1}/(2\pi)\approx$5.223\text{\,}\mathrm{GHz} for Samples A and B, and $C_{\textrm{l}}$ is the capacitance per unit length. Consequently, the critical damping is obtained with $\kappa_{2}=4g\approx 2\pi\times$29\text{,}\mathrm{MHz}. The external quality factor corresponding to the leakage from the resonator R1 to the transmission line through the capacitances $C_{\textrm{TL}}$ is given by Göppl *et al.* ([2008](#bib.bib70)) $Q_{\textrm{ext}}=2x_{1}C_{\textrm{l}}/(4Z_{\textrm{L}}\omega_{1}C_{\textrm{TL}}^{2})\approx 9\times 10^{5}$. Consequently, the corresponding damping rate is $\kappa_{\textrm{ext}}=\omega_{1}/Q_{\textrm{ext}}\approx 2\pi\times$6\text{\,}\mathrm{kHz}. The loaded quality factor of the second mode of R1 is approximately $Q_{\textrm{L}}=2\times 10^{4}$, when the resonators are far detuned. Thus, the internal losses in R1 dominate over the losses to the transmission line, $Q_{\textrm{int}}\approx Q_{\textrm{L}}$. Furthermore, we obtain the damping rate $\kappa_{1}=\omega_{1}/Q_{\textrm{L}}\approx 2\pi\times$300\text{,}\mathrm{kHz}. The real part of the complex wave propagation coefficient, $\gamma=\alpha+i\beta$, describes the damping in the waveguide, and it can be calculated as Göppl *et al.* ([2008](#bib.bib70)) $\alpha=n_{\textrm{m}}\pi/(4x_{1}Q_{\textrm{int}})\approx$7\text{\times}{10}^{-3}\text{\,}{\mathrm{m}}^{-1}, where $n_{\textrm{m}}$ is the mode number with $n_{\textrm{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 $\kappa_{2}$ values since $\kappa_{\textrm{int,1}}\gg\kappa_{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 $\kappa_{\textrm{int,2}}\approx\kappa_{2}(0)=2\pi\times$16\text{,}\mathrm{MHz}$$.

The photon number inside R1 when R2 is far detuned can be estimated as Oelsner et al. (2017) $n=4\Omega_{\textrm{d}}^{2}/\kappa_{1}^{2}\approx 10$, where the driving strength is given by $\Omega_{\textrm{d}}=C_{\textrm{TL}}V_{\textrm{in}}V_{1}/\hbar$, the input voltage is obtained from the input power as $V_{\textrm{in}}=\sqrt{P_{\textrm{in}}Z_{\textrm{L}}}$, and the input power is $P_{\textrm{in}}\approx-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_{\textrm{b}}=0$, the total photon number is approximately equally divided between the resonators if $\kappa_{1}\approx\kappa_{2}$. However, in our samples $\kappa_{1}<\kappa_{2}$, especially at $V_{\textrm{b}}>2\Delta/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_{\textrm{in}}\approx-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 $\kappa_{\phi}\mathcal{L}[\hat{b}^{\dagger}\hat{b}]\hat{\rho}$, where the dephasing rate $\kappa_{\phi}$ depends on the number of thermal photons in the resonator. Similarly, in the case of superconducting qubits, the dephasing can be written as $\kappa_{\phi}\mathcal{L}[\hat{\sigma}_{z}]\hat{\rho}$, where $\hat{\sigma}_{z}$ is a Pauli operator. The factor $\kappa_{\phi}$ causes a similar effect as $\kappa_{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 Clerk and Utami (2007); Goetz et al. (2017) $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 $\kappa_{\phi}=\omega_{\phi}n(n+1)$, where $\omega_{\phi}$ 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 Silveri et al. (2017)

[TABLE]

where we assume that the photon number of the effective bath, to which R2 is coupled through $\kappa_{\textrm{int,2}}$, vanishes owing to the very low cryostat temperatures of approximately $10\text{\,}\mathrm{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\text{\,}\mathrm{K}$, and therefore, the quasiparticle density remains low in it. However, the critical temperature of Al approximately $1.2\text{\,}\mathrm{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 Barends et al. (2011); O’Neil et al. (2012) $\kappa_{\textrm{qp}}\propto n_{\textrm{qp}}\propto\sqrt{P}$, where $P$ is the absorbed power. The Al leads at the NIS junctions receive half of the Joule power $P=IV_{\textrm{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, $\kappa_{\textrm{qp}}=\omega_{\textrm{qp}}n$, where $\omega_{\textrm{qp}}$ is a proportionality coefficient. The dc power dissipated in the junctions is substantially higher than the microwave input power. At $eV_{\textrm{b}}/(2\Delta)=2$, the dc power is approximately $30\text{\,}\mathrm{pW}$ compared to a microwave power of $-100$ $\textrm{dBm}=$0.1\text{,}\mathrm{pW}$$. The normal metal in the QCR acts as an effective quasiparticle trap O’Neil et al. (2012) minimizing the quasiparticle losses. Some fraction of the power dissipated at the QCR leaks to the SQUID.

There is an approximately 10-$\mu$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_{\textrm{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 $\kappa_{\textrm{res}}=\omega_{\textrm{res}}n$ with a proportionality coefficient $\omega_{\textrm{res}}$.

Consequently, the total voltage-dependent losses in R2 including the dephasing, quasiparticle losses in the superconductors and the resistive losses are given by

[TABLE]

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 photon-number-dependent losses into account as

[TABLE]

with only one fitting parameter $\omega_{\textrm{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 $\omega_{\textrm{r,tot}}\approx 2\pi\times$22\text{,}\mathrm{MHz}$$. The good agreement with the experimental damping rate $\kappa_{2}$ and the model with the quadratic residual losses $\kappa_{\textrm{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_{\textrm{b}}$. Similar dependence can be observed also for the even modes at $\Phi/\Phi_{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 $\kappa_{\textrm{r,2}}$. However, we do not quantitatively model these losses. Instead, they are effectively included in the parameter $\omega_{\textrm{r,tot}}$ in Eq. (29).

Full model for $\kappa_{2}$ and $\kappa_{\textrm{eff}}$

The parameters $\kappa_{2}$ and $\kappa_{\textrm{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_{\textrm{R2}}$, for the obtained effective resistance, as discussed below. The coupling rate is related to the quality factor as $\kappa_{2}=\omega_{2}/Q_{\textrm{R2}}$. The full model denoted by the line in Fig. 4(a) is obtained by fitting

[TABLE]

to the experimental transition rates according to Eqs. (12), (27), and (29). Here, we use $\omega_{\textrm{r,tot}}$ as the only fitting parameter since we fix $\kappa_{\textrm{int,2}}$ to the saturation value at zero bias, as discussed above.

Subsequently, we may proceed to the effective damping rates $\kappa_{\textrm{eff}}=-2\textrm{Im}(\lambda_{\pm})$, which can be obtained from $\kappa_{2}$ with the help of Eq. (2). The damping rates above the critical damping, $\kappa_{2}>4g$, result in the two branches at bias voltages $V_{\textrm{b}}\gtrsim 2\Delta/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. Partanen et al. (2018). We analyze the samples using standard microwave circuit analysis Pozar (2011). The input impedance of the resonator R2 is

[TABLE]

where the impedance of the SQUID and the capacitors between the SQUID and the center conductor is given by $Z_{\textrm{S}}=i\omega L_{\textrm{S}}+2/(i\omega C_{\textrm{S}})$, the impedance of the coupling capacitor between the resonators by $Z_{\textrm{C}}=1/(i\omega C_{\textrm{C}})$, $\gamma$ 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_{\textrm{eff}}$. The inductance of the SQUID is calculated as $L_{\textrm{S}}(\Phi)=\Phi_{0}/[2\pi I_{0}|\cos(\pi\Phi/\Phi_{0})|],$ where the maximum supercurrent through the SQUID is $I_{0}$, and the flux quantum is $\Phi_{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 Pozar (2011)

[TABLE]

where $Z_{\textrm{L}}$ is the characteristic impedance of the external measurement cables, and

[TABLE]

with

[TABLE]

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 Petersan and Anlage (1998).

Sample fabrication

The samples are fabricated on a Si wafer with a thickness of $500\text{\,}\mathrm{\SIUnitSymbolMicro m}$ and a diameter of $100\text{\,}\mathrm{m}\mathrm{m}$. First, a 300-nm-thick layer of SiO2 is thermally grown on the wafer with resistivity $\rho>$10\text{,}\mathrm{k\SIUnitSymbolOhm}\text{,}\mathrm{cm}$$. Subsequently, a 200-nm-thick 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 Al2O3 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\text{\,}\mathrm{nm}$ each. The first Al layer is oxidized in situ in the evaporation chamber at $1.0\text{\,}\mathrm{mbar}$ for $5\text{\,}\min$. The SINIS junctions consist of Al ($40\text{\,}\mathrm{nm}$) and Cu ($40\text{\,}\mathrm{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 dry-dilution 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\text{\,}\mathrm{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 $\tau\approx$50\text{,}\mathrm{ns}$$ in the measurement setup outside the sample by multiplying with $\exp(i\omega\tau)$. 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}|$.