Quantum bath engineering of a high impedance microwave mode through quasiparticle tunneling

In microwave quantum optics, dissipation usually corresponds to quantum jumps, where photons are lost one by one. Here we demonstrate a new approach to dissipation engineering. By coupling a high impedance microwave resonator to a tunnel junction, we use the photoassisted tunneling of quasiparticles as a tunable dissipative process. We are able to adjust the minimum number of lost photons per tunneling event to be one, two or more, through a dc voltage. Consequently, different Fock states of the resonator experience different loss processes. Causality then implies that each state experiences a different energy (Lamb) shift, as confirmed experimentally. This photoassisted tunneling process is analogous to a photoelectric effect, which requires a quantum description of light to be quantitatively understood. This work opens up new possibilities for quantum state manipulation in superconducting circuits, which do not rely on the Josephson effect.

In order to evaluate the properties of the resonant modes of the system, we simulate the sample using the Sonnet software. The resonance frequency and the characteristic impedance of each mode are calculated from the admittance Y (ω) seen by the junction. We first numerically simulate the admittance matrix Y ij (ω) of the two port network that is obtained by placing a first port at the position of the junction and a second one on the 50 Ω measurement line. From the admittance matrix, we then compute the admittance seen by the junction by eliminating the second port and taking into account the capacitance of the junction C J . We arrive at where Y 0 = 1/50 Ω −1 is the admittance of the measurement line. The resonant frequencies ω n of the modes are obtained by looking for the zeros of Im[Y (ω)], while their characteristic impedance is given by [1] Z n = 2i ω n ∂Y (ω) ∂ω ω=ωn −1 .

(S.2)
In order to estimate C J and the precise value of the grAl sheet inductance, we compare the measured frequencies of the modes to the simulated ones. We find that a value of C J = 1.75 fF and a sheet inductance of 0.56 nH/□ explains well our results. It is compatible with the 150 × 150 nm 2 area of our junction. The measured values of the resonance frequencies are obtained from photo-assisted current measurements. The circuit is dc-biased below the gap and the current through the junction is measured as a function of the frequency of an incoming microwave signal. We observe sharp resonances that we identify to the modes coupled to the junction. In table S1, we report the measured resonance frequencies as well as the simulated properties of the modes. The mode with the largest impedance probed in the experiment corresponds to the fundamental mode of the quarter wavelength resonator, which is labeled with the index n = 1. Because of the reduced length of the second section of the DBR with respect to the Bragg condition, there exists a lower frequency mode, which resonates at 1.9 GHz. This mode has a relatively large impedance, even though it is not fully localized into the last high impedance section of the sample. The modes with even n = 2, 4, 6 have a voltage minimum at the position of the junction and therefore a small impedance Z n as expected for a quarter wavelength resonator. After n ≥ 7, the mode structure is affected by the DBR and becomes more complicated, we only list the modes that have a significant impedance. Note that the mode impedance as defined in (S.2) is not an intrinsic property of the mode, as in the usual sense of a mode impedance. Here, it depends on the position of the junction and quantifies the vacuum quantum voltage fluctuation as seen by the junction. This definition coincides with the usual one when the junction is located at the maximum of voltage.
Finally, we estimate the Kerr coefficient of the n = 1 mode due to the non-linearity of the kinetic inductance. Test wires with the same geometry as the resonator and evaporated in the same conditions exhibit a critical current density of approximately 1 mA µm −2 . Following [2], we estimate from this measurement that the Kerr coefficient is approximately 2π × 200 Hz, which is negligible.

II. SAMPLE FABRICATION
The sample is micro-fabricated on a 300 µm thick, thermally oxidized, silicon substrate through the following steps.

Granular Aluminum resonator and DBR structure patterning
The design is written on a 1.3 µm thick trilayer resist stack (PMMA 495/PMMA 950/PMMA 950) using electronic beam lithography. After revelation, 10 nm of grAl are deposited by e-beam evaporation at 0 • angle ( fig. S2a). The oxygen pressure during the evaporation is adjusted to 1.1 × 10 −5 mbar and the evaporation rate is 1.5Å/s. In a following evaporation, without exposing the sample to air, we deposit 30 nm of pure Al with an angle of 45 • in the direction perpendicular to the one of the wires ( fig. S2b). Because of the the resist thickness and the evaporation angle, only regions with a width exceeding 1.3 µm are covered during this second evaporation. This corresponds, from left to right, to the measurement line, the second DBR section and a final 2.5x2.5 µm 2 square contact that is used to connect the junction.
In a second step, we pattern the remaining part of the 50 Ω measurement line and a central ground plane, which is used to connect the other side of the junction (fig. S2c). These structures are realized through optical lithography  and liftoff of a 60 nm thick Al layer. In order to obtain a good dc contact with the previously deposited measurement line, an argon ion beam etching step is used before the evaporation. At the end of the process a 30 nm layer of Nb is deposited at the back of the sample for the ground plane.

Tunnel junction patterning
The Al/AlOx/Al junction is realized through a standard double angle evaporation using a PMMA bilayer resist without any suspended bridge ("Manhattan technique", see fig. S2d). The width of the junction arms is 150 nm and the oxidation is performed at 10 mbar oxygen pressure (static) for 20 minutes. In a final step, we connect the junction contacts to both the resonator and the ground plane using a last step of optical lithography, followed by Ar ion beam etching and the evaporation of 80 nm of Al ( fig. S2e).

III. EXPERIMENTAL SETUP
The sample is glued on a gold-plated copper sample holder and bonded to a microwave PCB. The sample holder is fixed on the 10 mK stage of a Cryoconcept dry dilution refrigerator. The experimental setup used to measure the sample is shown in figure S4. All the microwave measurements presented in the manuscript are performed with a vector network analyzer connected to the RF input and output ports. The microwave drive is sent to the sample through an attenuated coaxial line and its reflection is collected by a different line using a circulator. The reflected signal is amplified by a cryogenic HEMT amplifier at 4 K followed by two room temperature amplifiers. In order to dc-bias the junction, we use a voltage divider followed by a filtering capacitor and a bias tee. A 10 kΩ resistance is mounted in series with the sample to measure the current through the junction.

IV. JUNCTION ADMITTANCE
When the impedance of the mode coupled to the junction is small compared to the quantum of resistance (λ ≪ 1), the resonator field can be treated classically. In this case, the effect of charge tunneling can be modeled by treating the junction as a frequency and voltage dependent admittance Y J [3]. This admittance also varies with the microwave amplitude seen by the junction. Supposing a weak microwave amplitude, the real and imaginary parts of the admittance can be expanded to first order in the microwave amplitude as [3,4]: Here, I is the quasiparticle current that is defined as where n R (n L ) is the dimensionless density of states in the right (left) contact, which tend to one at voltages far below and far above the gap, and f (V ) is the Fermi distribution. quantities are given by These expressions can be derived using equations (2) and (3) of the main text in the limit of small λ. In a quantum approach, they correspond to the usual single photon loss and frequency shift that arise when an harmonic oscillator mode is linearly coupled to a bath. In the strong impedance regime, where λ ≈ 1, this approach is not valid, and a quantum treatment of the field, and in particular of the displacement operator appearing in the tunnel Hamiltonian, is necessary.

V. QUANTUM MASTER EQUATION
In order to describe the quantum dynamics of the resonator field in the high impedance regime, we use a master equation formalism and treat the junction as a bath and the resonator as an open quantum system. Following the standard quantum optics approach based on the Born-Markov and on the secular approximation [5], one obtains a master equation in the Lindblad form for the density matrix of the resonator modes as detailed, for example, in [6]. The master equation is written in terms of jump operatorsÂ l operators that come from the expansion of the displacement operator exp[iλ(â +â † )] in the Fock state basis. Their effect is to create l photons when l > 0 and annihilate −l photons when l < 0. They are defined, for positives l, aŝ A l>0 = ∞ n=0 |n + l⟩ ⟨n + l| e iλ(â+â † ) |n⟩ ⟨n| . (S.8) The operators with l < 0 are defined throughÂ −l = (−1) lÂ † l . For l ≥ 0, the matrix elements ofÂ l can be obtained from the following identity n are the generalized Laguerre polynomials. In the main text, we define α nl = | ⟨n + l|e iλ(â+â † ) |n⟩ | 2 , which can be rewritten in terms ofÂ l operators as α nl = | ⟨n + l|A l |n⟩ | 2 for l ≥ 0. When the non-linearity is strong (λ ≈ 1), quantum jumps, generated by theÂ † l operators, lead to transitions from |n⟩ to |n − l⟩. This is in stark contrast to the usual quantum optics situation, where the two possible jump operators areâ andâ † , which only couples adjacent Fock states. In the case where many modes are coupled to the junction, jump operators can be defined for each mode using the value of λ obtained from numerical simulations of the mode impedance seen by the junction (see I). In order to obtain more compact expressions, we definê A l0,l1,... =Â l0 ⊗Â l1 ⊗ . . . , (S.10) whereÂ ln is the operator defined in (S.8) associated to mode n. The other ingredients entering the master equation are the real and imaginary part of the single sided Fourier transform of the bath correlation function. In the case of a tunnel junction, one obtains With these definitions, the master equation for the density matrix of the modes is (S.14) The first term of equation (S.13) describes the Lamb shift of the resonator energy levels with associated Hamiltonian H LS , while the second term corresponds to photo-assisted tunneling events, during which one electron tunnels and photons are absorbed and/or created in the resonator modes. As required by gauge invariance, the master equation is unchanged by V → −V . While deriving this equation, we assume that the mode frequencies are not commensurate when performing the secular approximation. The bath functions γ(V ) and ϵ(V ) can be related to the I(V ) characteristic of the junction by noticing that The last approximations are valid for positive bias and low temperature. With these approximations, the master equation can be rewritten as This master equation is only valid for positive bias and as long as the energy shift n l n ℏω n remains small compared to the gap 2∆, which is the case in our experiment. This multimode master equation is in general difficult to solve. In the quantum Zeno experiment presented in the main text, only the 6 GHz mode is pumped. We therefore assume that the other modes do not participate to the dynamics and remain in vacuum. When tracing out these modes, we have to compute The modified current characteristicĨ takes into account the dynamical Coulomb blockade (DCB) due to the other modes. It is defined asĨ The multi-mode DCB factor is the product of the blockade factors | ⟨0|exp iλ n (â n +â † n ) |0⟩ | 2 = e −λ 2 n for each mode, except the 6 GHz mode. As mentioned in the main text, this is equivalent to an increase of the tunnel resistance by a factor Π n̸ =1 e λ 2 n . The master equation can be further simplified by keeping only the l n = 0 term in the sum over the other modes. This assumes that the effect of these modes is simply to renormalize the tunnel resistance. We then obtain For voltages below the gap, the simplification of the jump terms only neglects the contribution of the 1.9 GHz mode. In particular, we neglect processes where photons in the 6 GHz are converted into photons at 1.9 GHz and the remaining energy is taken by the tunneling electron (see figure S6). The higher frequency modes have no influence because of energy conservation. Computing the loss rate for state |n⟩ using (S.26) leads to equation (2) of the main text. The expression ofĤ LS is more affected than the jump terms, because all the modes contribute as discussed in more details below (see VII).

VI. ZENO DYNAMICS AND PUMP CALIBRATION
We now use the previously derived master equation to explain the quantum Zeno dynamics observed in the experiment. Quantum Zeno dynamics happens when we dc-bias the junction at 363 µV, which lies in the voltage region where the l ≥ 2 processes are the dominant loss process. We measure the resonator reflected spectrum S 11 for different microwave powers (see inset of figure 3b of the main text) and extract the mode amplitude from these measurements. This requires to know the pump amplitude, which is defined as η = κ c P/ℏω, where κ c is the coupling loss rate and P is the incoming microwave power on the sample. The coupling κ c is measured independently. The incoming power P depends on the total attenuation of the experimental setup that needs to be calibrated. We do so by measuring the resonance fwhm and by comparing it to the one of a two level system. In a two level system, the fwhm is equal to √ κ 2 + 4Γ 2 , with κ being the total loss rate and Γ the power broadening, which varies with the pump amplitude η as Γ = √ 2η. By assuming the same power dependence on the measured Γ, we arrive at the attenuation of the line that pumps the resonator.
We extract the loss rate κ and the broadening Γ by fitting the spectra measured at different power using the two-level formula [7] where δ is the frequency detuning between the pump and the 0 → 1 resonance. In figure S7, we plot the spectrum measured at P = −162 dBm together with the fitted S 11 (red curve). This procedure is repeated for every power, and by fitting the slope of Γ with power, we obtain the attenuation and the value of η as used in figure 3 of the main text. Once the pump is calibrated, we extract the resonator field intensity from the measured spectra through a b FIG. S7: a Squared modulus and b phase of the reflection spectrum measured at P = −162 dBm. The red curve shows the fitted spectrum assuming a two-level system model.
A plot of |⟨â⟩| 2 at resonance is shown in figure 3a of the main text. We compare this measured intensity to the predictions of the following master equation This master equation is the simplified single mode master equation obtained in (S.26), where we add a pump term with amplitude η, a detuning term proportional to δ, and also include a single photon loss rate κ in order to take into account single photon loss from various origins. This loss rate is κ = 2π × 4.2 MHz (κ = 2π × 2.7 MHz) for V = 363 µV (V = 284 µV). These values are chosen to reproduce the observed resonance width at low pump power. The expectation value ofâ is then obtained from the steady state solution of the master equation (calculated using the QuTip Python software [8]) and allows us to plot the curve in figure 3a of the main text. We also compute the expected S 11 and apply the same fitting procedure as for the experimental data and obtain the theoretical power broadening Γ (shown in figure 3b of the main text) and the loss rate κ (plotted in figure S8).
In figure S9, we compare the results of simulations including or not the Lamb shift term H LS . The non-linearity introduced by the Lamb shift significantly contributes to the blockade of the 1 → 2 transition in addition to the Zeno effect. The figure of merit of the observed blockade can be quantified by the range of pump amplitude over which the mode behaves as a saturated two level system. Supposing that the mode is subject only to one and two photon loss, without any non-linearity, this range is set by the ratioγ 2 /γ 1 , whereγ n is the effective loss rate of state |n⟩ in this model. We find that, without Lamb shift, the ratioγ 2 /γ 1 is about 15, withγ 2 being 2π × 65 MHz and γ 1 = 2π × 4.2 MHz. When we include the Lamb shiftγ 2 increases approximately to 2π × 100 MHz with the ratiõ γ 2 /γ 1 being about 25.
Finally, we simulate a master equation modeling ideal Zeno dynamics (see figure S9). We suppose that the Hamiltonian is given by equation (S.31) without the Lamb shift term. We then include two jump operators. The first one describes single photon loss with rate κ and the second one is the projector |0⟩⟨0| + |1⟩⟨1|, which measures if the state is in the {|0⟩ , |1⟩} subspace. The projection rate γ P is adjusted to fit the data in the region η < 2π × 3 MHz, we obtain γ P = 2π × 215 MHz, which is comparable to the rate of 2π × 100 MHz found above. FIG. S8: Comparison between the measured (blue data) and simulated (red curve) evolution of the loss rate κ. The curves show that κ remains constant, as in a two level system, as long as the Zeno effect limits the dynamics of the system to the first two Fock states. When the Zeno dynamics break down, the loss rate increases quadratically with the pump amplitude.

VII. LAMB SHIFT
From the expression ofĤ LS derived in V, the Lamb shift δω n of the |n⟩ state of the 6 GHz mode, with all the other modes in vacuum, is given by As discussed above, the main effect of the other modes is to renormalize the tunnel resistance, but the remaining terms modify the V dependence beyond a simple renormalization of R T . Including all the modes presented in I, we use this expression to compute the shift of the |0⟩, |1⟩ and |2⟩ states, from which we deduce the curves plotted in figure 4 of the main text. As previously discussed in V, a natural approximation consists in keeping only the l n = 0 terms in the sum, which leads to δω n = − 1 2e which corresponds to equation (4) of the main text. Keeping only the l = 0 term, we obtain a semi-classical expression, which is equivalent to the prediction of a classical model where the tunnel resistance is renormalized as R T → R T Π n e λ 2 n . We compare these different approximations and the exact expression in to our data in figure S10. When we compare our data to the theoretical predictions, as here or as in figure 4 of the main text, we allow for a small frequency shift of the transition frequency (vertical offset) and also for a small shift of the gap (horizontal offset). With these adjustments, our data are in better agreement with the exact multi-mode calculation (S.33) than with the simplified formula obtained in (S.36). FIG. S10: Comparison of the different predictions for the shift δω 01 of the fundamental resonance together with the experimental data (blue circles). The red curve corresponds to the multi-mode prediction (S.33), the green curve to the single mode approximation (S.36) and the orange curve to the semi-classical model, which is obtained by keeping only the l = 0 term in (S.36).