A bright on-demand source of anti-bunched microwave photons based on inelastic Cooper pair tunneling

The ability to generate single photons is not only an ubiquitous tool for scientific exploration with applications ranging from spectroscopy and metrology to quantum computing, but also an important proof of the underlying quantum nature of a physical process. In the microwave regime, emission of anti-bunched radiation has so far relied on coherent control of Josephson qubits, where precisely calibrated microwave pulses are needed, and the achievable bandwidth is limited by the anharmonicity of the qubit. Here, we demonstrate the operation of a bright on-demand source of quantum microwave radiation capable of emitting anti-bunched photons based on inelastic Cooper pair tunneling and driven by a simple DC voltage bias. It is characterized by its normalized second order correlation function of $g^{(2)}(0)\approx0.43$ corresponding to anti-bunching in the single photon regime. Our source can be triggered and its emission rate is tunable in situ exceeding rates obtained with current microwave single photon sources by more than one order of magnitude.

Josephson photonics has recently emerged as a way to directly generate and manipulate microwave frequency signals at milliKelvin temperatures without the need for complicated microwave control drives [9][10][11][12][13][14][15][16] . It relies on the phenomenon of inelastic Cooper pair tunneling where a DC voltage biased Josephson tunnel junction at zero temperature can admit a finite direct current, even if the applied bias voltage is smaller than its gap voltage 17 . Although the finite voltage bias makes it impossible for Cooper pairs to tunnel elastically in this regime, they can tunnel inelastically while dissipating their surplus energy into the electromagnetic environment of the junction. The resulting excitations of the environment are photons at microwave frequencies 18,19 . This effect can be harnessed and the properties of the emitted radiation can be controlled by presenting the junction with a specifically tailored electromagnetic environment.
Specifically, a Cooper pair of charge 2e can tunnel inelastically through a Josephson junction biased at voltage V b if the electromagnetic environment of the junction has modes that can absorb the energy difference of 2eV b = hν J = i n i hf i , where ν J = 2eV b /h is the Josephson frequency and n i is the number of photons emitted into a given mode of frequency f i . The power spectral density in units of photons of the emitted radiation at frequency f is then given by the expression 19 : Here, I c is the critical current of the junction, Re Z(f ) is the real part of the impedance describing the electromagnetic environment and P (hf ) is the probability density for a tunneling Cooper pair to exchange an energy hf with this environment 17,20 . The latter quantity is di-rectly proportional to the Cooper pair tunneling rate Γ: While early experiments focused on the rate Γ and the associated direct current 18,[21][22][23] , recent advances in low-noise high-frequency measurements have made it possible to investigate the photonic side of this energy transfer [14][15][16]19,24 . In all previous implementations, the tunneling events were happening either independently 18,19 or through stimulated emission [14][15][16] , in both cases leading to classical Poisson statistics of the emitted photons. In the present work, we demonstrate a device where the electromagnetic environment of the junction is engineered to create anti-bunching in the photon emission, thus showing quantum statistics generated through inelastic Cooper pair tunneling.
A schematic circuit-diagram of our device is shown in the bottom part of Fig. 1a. Its main element is a superconducting quantum interference device (SQUID) made of two parallel 0.02 µm 2 NbN-MgO-NbN Josephson junctions 25 (Fig. 1d). We used a fast flux line to tune its critical current and thus the tunneling rate in situ. One side of the SQUID connects to a quarter-wave transmission line resonator (blue in Fig. 1a), with a fundamental frequency f 0 = 6 GHz, which is followed by an on-chip bias tee and beamsplitter (Fig. 1a, b and Supplementary Information). The other side of the SQUID is grounded through an on-chip parallel RC-circuit (red, for details see Supplementary Information), with resistance R, capacitance C, and a small spurious inductance L p .
We extracted the calibrated and time-resolved autoand cross-correlations (G (1) , G (2) ) 13,26,27 as well as the power spectral density (PSD) of the microwave radiation The SQUID (yellow) is connected to an RC-element (red, with parasitic inductance Lp) and a transmission line resonator with a fundamental mode at f0 (blue). Voltage and flux biases are applied to the SQUID through a voltage divider and an on-chip fast flux line respectively. Photons escape from the resonator through an on-chip bias tee and beamsplitter into two measurement chains each containing an RF switch for calibration (marked "cal") as well as filters and isolators to protect the sample from amplifier noise (Supplementary Information). Power spectral density (PSD) and the photon correlation functions G (1) and G (2) are calculated numerically after down-mixing and digitization ( Supplementary Information). b, Optical image overview of the entire sample as indicated in the framed part of panel a. SQUID and RC-circuit are inside the red frame. The rest of the circuit implements resonator, on-chip bias tee and beamsplitter (Supplementary Information) branching out into two high frequency (RF) outputs, and a DC input. c, Scanning electron micrograph (SEM) of the SQUID (small yellow frame) and RC-circuit (red). d, SEM of the SQUID loop consisting of two vertical NbN-MgO-NbN Josephson junctions. The fast flux line is visible on three sides. Where it is underneath the main circuit, it implements the parallel plate capacitor of the RC. e, Schematic depiction of the inelastic tunneling process. Horizontal black lines represent the Cooper-pair condensates in the superconductors on either side of the insulating layer (gray) of the voltage biased SQUID. The choice of voltage bias (2eV b = hf0 + EC) is such that tunneling with simultaneous photon emission is possible when the capacitor is not already charged. f, Immediately after a tunneling event the energy necessary to charge the capacitor with another Cooper pair is 3EC, blocking further Cooper pair tunneling for a time ≈ RC. emitted by our device using the measurement setup depicted schematically in Fig. 1a and described in detail in the Supplementary Information. The PSD divided by hf has units of photons and corresponds to the quantity computed by equation (1).
The underlying principle of our source is that one tunneling event necessarily acts back onto the next one in order to create anti-correlations in the Cooper-pair current, leading to anti-bunching in the photon emission. On timescales shorter than RC, a tunneling Cooper-pair has to charge the island formed between capacitor and SQUID. To emit a photon, the voltage bias must be chosen such that hν J = hf 0 + E C , where E C = (2e) 2 /2C is the charging energy of the island (Fig. 1e). However, immediately after a first tunneling event the energy necessary to add another Cooper-pair to the island is (2 · 2e) 2 /2C − (2e) 2 /2C = 3E C . This makes the voltage bias insufficient for further tunneling (Fig. 1f). A second Cooper pair can tunnel and emit a photon only after a time RC, when the island is discharged, leading to the desired anti-bunching. This picture is valid if the charging energy associated with one Cooper pair is large compared to the thermal energy fluctuations, i.e. E C k B T eff , as well as the lifetime broadening of the charging energy, i.e. E C /(2RC), or equivalently R h/(4e 2 ). Here, T eff is the effective temperature of E C hf 0 hf 0 hf 1 Figure 2. Power spectral density a, Power spectral density in units of photons measured close to minimum Ic(Φ) (SQUID flux bias of Φ ≈ 0.5 Φ0, with Φ0 = h/2e the flux quantum) as a function of frequency and voltage bias (expressed as the Josephson frequency). The dashed line corresponds to the process hνJ = EC + hf , where EC ≈ 1.5 GHz is the energy required to charge the capacitor with a Cooper pair. The frequency axis is centered around f0 ≈ 6 GHz, the resonance frequency of the transmission line resonator. b, Real part of the impedance as a function of the voltage bias directly extracted from the data shown in panel a (blue dots) and plotted from a fit of our circuit model to the data (solid green line). c, Power spectral density in units of photons at 6 GHz as a function of voltage bias and flux bias in units of flux quanta in the SQUID loop. The horizontal lines are drawn using the values of EC and f0 extracted from panel a and match the observed features. The underlying emission process is sketched for each line with the dashed one corresponding to the first order process hνJ = EC + hf0 and the dotted ones corresponding to second order processes hνJ = EC + 2hf0 and hνJ = EC + hf0 + hf1, where f1 is the frequency of the next higher resonator mode. The turquoise dotted line delimits the region in which no voltage bias can be applied because the SQUID is latching to its current branch (see text). Two filled white circles mark the starting and end-point of the flux pulse giving the results of Fig. 4 and our estimate of the actual trajectory is sketched in white.
the low-frequency electromagnetic environment.
We now verify these conditions and extract the device parameters. To this end, we analyze the PSD in Fig. 2a taken at a flux bias Φ ≈ 0.5Φ 0 , where Φ 0 = h 2e is the superconducting flux quantum. There, I c (Φ) and thus emission rates are low so that tunneling events can be considered independent and equations (1) and (2) are valid 17 . For this flux bias, photon emission is maximal at ν J = 7.5 GHz and f = f 0 = 6 GHz. The difference between the energy hν J provided by a tunneling Cooper pair and the energy hf of the detected photons is the charging energy E C ≈ h × 1.5 GHz ≈ k B × 72 mK (C ≈ 51 fF). Further analysis of this PSD (Supplementary information) yields the real part of the resonator impedance (blue dots in Fig. 2b) and the effective temperature of the low-frequency electromagnetic environment T eff ≈ 21 mK. The resistance R ≈ 32.1 kΩ was determined from an independent DC measurement (Supplementary Information). Additionally, we verify the accuracy of the electrical model presented in Fig. 1a by using it to fit equation (1) to the PSD. The real part of the resonator impedance given by the resulting circuit parameters is shown as a solid green line in Fig. 2b (for detailed parameters see Supplementary Information). It has a full width at half maximum of 575 MHz and agrees well with the measured data. The ≈ 200 MHz-periodic impedance modulations are due to reflections in our output lines, which are not included in the model.
In Fig. 2c we now explore the behavior of the device when the critical current is increased to values relevant for operation, where Cooper pair tunneling events cannot be considered independent any more. To do so, we measure the power spectral density at its maximal value in frequency (at f 0 ) as a function of ν J and the flux bias Φ. The brightest features appear around ν J = 7.5 GHz, corresponding to the desired process. At slightly lower flux bias additional features appear at ν J = 13.5 GHz. They originate from the emission of two photons per tunneling Cooper pair into the mode at f 0 . When the critical current is maximal (Φ = 0), another feature becomes visible at ν J = 25.4 GHz. It corresponds to processes where one Cooper pair emits one photon into mode f 0 and one into the next higher mode of the quarter-wave transmission line resonator at f 1 ≈ 3f 0 . Strikingly, the processes at 7.5 GHz and 13.5 GHz disappear around Φ = 0 where one would expect the rates to be highest. We instead observe a dark zone delimited by the dotted parabola. In this region the critical current is high enough for the SQUID to get trapped in a Bloch oscillation regime [28][29][30][31] , a rudiment of the zero-voltage state observed in larger Josephson junctions. In this state the voltage drops mostly over the RC-element, decreasing the voltage difference across the SQUID below the threshold for photon emission at f 0 . This interpretation is confirmed by an independent measurement of the resonator frequency showing flux-tunability in the region delimited by the dotted line ( Supplementary Information).
We now focus on the key question of this work and investigate the statistics of the radiation emitted by the device. To do so, we measure its normalized second-order correlation function g (2) (t, τ ) given by 32 In the above expression, G (2) (t, τ ) is the unnormalized second-order correlation function dependent on a time t with respect to a reference and on the time delay τ , defined as: The operatorsâ out andâ † out are the annihilation and creation operators of the outgoing field in the transmission line. The denominator of the right-hand side of equation (3) is a normalization factor dependent on the firstorder correlation function G (1) (t, 0) gives the photon emission rate of our device. In the absence of a well defined time reference, an average over t is performed on G (1) (t, τ ) and G (2) (t, τ ) yielding G (2) (τ ), G (1) (τ ) and g (2) (τ ) 32 . Fig. 3a shows the total photon emission rate G (1) (0) of our sample in a region around the one-photon peak visible in Fig. 2c. It is measured by integrating the PSD over frequency between 4.22 GHz and 8 GHz (Supplementary Information). We have evaluated g (2) (τ ) at different points along two lines of constant voltage (flux) bias as indicated by plus (cross) symbols, corresponding to the curves shown in Fig. 3b (Fig. 3c).
Close to Φ ≈ 0.5 Φ 0 , the g (2) (τ ) function shows a marked dip down to approximately 0.5 at τ = 0 and is close to the expected value of 1 elsewhere, indicating that our sample indeed produces strongly anti-bunched photons. Fig. 3b shows that, as we approach lower flux biases to increase the critical current, the sharp dip close to τ = 0 remains, but a broad bunching peak develops around it. We attribute this broad peak to random jumps between the bright voltage state of the junction and the dark zero-voltage state discussed above. This is consistent with the reduction in emission rate towards the lower Voltage bias left corner of Fig. 3a. The persistence of the sharp dip at τ = 0 demonstrates that the intended blocking mechanism remains functional, despite this effect. Fig. 3c shows the sensitivity of the anti-bunching to precise adjustment of the bias voltage. Here we chose the flux bias where we achieve the highest photon rate while bunching due to jumps to the dark state remains negligible: Φ ≈ 0.462 Φ 0 . When biased above 8 GHz, our sample emits bunched light (g (2) (0) > 1). This bunching effect can be understood by considering the energy diagrams in Fig. 1e, f at higher voltage biases: Even though the voltage is initially too high for the resonance condition to be fulfilled, thermal fluctuations still occasionally allow Cooper pairs to tunnel. When a first Cooper pair does tunnel, the RC-circuit is charged, and a second one can more easily follow, leading to bunching. As the bias decreases, so does the value of the second order correlation function, down to g (2) (0) ≈ 0.5±0.1(±3σ). The fact that the lowest value of g (2) (0) is reached below the maximum of the peak is expected: residual double emission events caused by thermal fluctuations of the bias voltage are further suppressed. The blue curve is taken at the maximum emission rate of 77 MHz (ν J = 7.4 GHz) along this cut. At this point we measure g (2) (0) ≈ 0.65 ± 0.06. We conclude that anti-bunching is robust as long as the bias voltage is kept at the nominal resonant value or below, but is rapidly transformed into bunching above.
The solid lines in Fig. 3c are numerical calculations up to fourth order in the critical current 33,34 using the electrical model of Fig. 1a with the extracted device parameters and an effective temperature of T eff = 40 mK. They reproduce well the observed anti-bunching signatures at the lowest bias voltages, including bumps in g (2) at approximately 2 ns but fail to fully explain the bunching signatures for the highest bias voltages. In addition, T eff used here to reproduce the data is significantly higher than the temperature extracted from the data in Fig. 2a at low critical current and at low bias voltage. This discrepancy could indicate that the resistor of the RC element heats up more than expected at the higher photon fluxes and higher voltages used here. Another explanation could be that correlations between more than two Cooper pairs are relevant and that calculations have to be performed beyond fourth order in the critical current, which significantly increases the computational effort and is left for future work.
So far we have focused on the free-running mode of operation, where latching to the dark state prevents us from reaching higher emission rates. However, as we now show, we can make use of this effect to produce photons on demand through the flux-pulsing scheme indicated in white on Fig. 2c. For this, the voltage bias is set to its nominal value (ν J = 7.4 GHz). We start out at a flux bias well in the dark region (left white dot), where no photon emission occurs in the stationary regime. We then apply a flux pulse which frustrates the SQUID (right white dot), suppressing its critical current and resetting it to the voltage branch. At this bias point photon emission is unlikely. Returning to the initial point in the dark region, removes the frustration and allows a photon to be emitted (and the capacitor to be charged). Then, the system is again trapped in the dark state limiting photon emission to at most one per cycle. Fig. 4 shows the data obtained when Gaussian flux bias pulses with FWHM of 1.5 ns are applied every 6 ns. Different pulse durations (within a factor 2) and lower repetition rates do not affect the results significantly. Higher repetition rates, however, cause photon pulses to overlap.
In Fig. 4a we first explore the photon emission rate G (1) (t, 0) as a function of time t with respect to the pulse. One averaged measurement yields the values marked by blue dots, which are separated by the sampling period of 1 ns of our measurement. By shifting the time delay between pulse-generation and the beginning of the sampling window, the photon emission rate can be resolved below the sampling period yielding the orange curve in Fig. 4a. The zero on the time axis matches the time t = 0 at which the measurements in Fig. 4b, c were taken. The width of this peak is due to a combination of the uncertainty in the photon generation time (jitter) and the resonator decay time. g (2) (0, τ ) c Figure 4. Flux-pulsed on-demand operation a, Photon emission rate as a function of time. Blue dots are the rates acquired during one averaged measurement, the orange line is obtained by repeating the measurement with different time delays of the flux-pulse. The shape of the peak is slightly skewed due to the finite time the photon spends in the resonator. The statistical errors in this panel and the next are small compared to the data points (Supplementary Information) b, First order correlation function at t = 0, as a function of the time delay τ (the solid line is a guide for the eye). As the actual time of emission is slightly after t = 0, the peak is centered at slightly positive τ . c, Second order correlation function at t = 0, as a function of the time delay τ . The green dots correspond to the unnormalized second order correlation function (right ordinate) and display periodic peaks every 6 ns given by the period of the applied flux pulse (the line is a guide for the eye). Blue dots with error bars (±3σ) show the normalized correlation function (left ordinate).
The complementary measurement of G (1) (0, τ ) ( Fig. 4b) characterizes the first order coherence of the source. The width of the peak is very close to the peak in G (1) (t, 0) in Fig. 4a, indicating that the broadening of G (1) (t, 0) is not dominated by jitter and that photons generated by our source may be indistinguishable. Also note that G (1) (0, 6 ns) ≈ 0, indicating that there is no coherence between successive pulses as expected for single photons.
In the unnormalized second order correlation function G (2) (0, τ ) this picture is reversed (Fig. 4c). Peaks appear every 6 ns, equal to the period of the applied flux-pulse, proving the on-demand aspect of our source. Note that the peak at τ = 0 is significantly lower than the oth-ers, indicating the anti-bunched character of the emitted radiation.
In order to quantify the anti-bunching we compute the normalized second order correlation function (blue dots) according to equation (3). This is only done for times with high enough emission rate G (1) (0 + τ, 0), in order to avoid divergence of the associated errors. For timedifferences τ = n × 6 ns with n = 0, we obtain g (2) ≈ 1, indicating that photons from two different pulses are independent. At zero time-delay, however, g (2) (0, 0) = 0.43±0.08(±3σ): the photons are strongly anti-bunched, in agreement with the mechanism presented in Fig. 1e, f.
We achieve here stronger anti-bunching than in the free-running case, likely due to the additional blocking effect given by the latching to the zero-voltage state in the dark region. At the same time, the pulsing scheme allows us to maintain very good quantum efficiency and photon flux (0.2 photons per pulse), because of the high emission rates at low flux bias. This makes it likely for a tunnel event to happen during each cycle even for very short flux pulses. We attribute the residual g (2) (0, 0) mainly to the low charging energy of our RC-circuit and its relatively low time constant (RC = 1.64 ns), only slightly larger than the decay time of the resonator τ = 0.28 ns. The latter is limited by practical considerations such as the instantaneous bandwidth of our measurement setup (Supplementary Information). In addition, a parasitic coupling between the flux bias and the RC-circuit, described by the parasitic inductance L p in Fig. 1a, causes the flux pulse to modulate the voltage across the junction. The dashed white line in Fig. 2c indicates the resulting expected trajectory in parameter space, which may not be ideal. Further optimization of this trajectory and modification of the coupling or application of simultaneous voltage and flux pulses could lead to better antibunching. Moreover, this type of device could possibly be optimized as a source of on-demand multi-photon Fock states by addressing the processes appearing at higher bias voltages.
In conclusion, we have demonstrated a Josephson photonics device producing strongly anti-bunched microwave radiation. In doing so we have shown that quantum statistics can be generated from inelastic Cooper pair tunneling. By modulating the effective critical current of the SQUID, using fast flux pulses and locking to a dark state after photon emission, we can produce anti-bunched photons on-demand at very high rates. Increasing the charging energy and the RC time, or replacing the RC-circuit by a high impedance resonant mode, should allow for significant improvements of anti-bunching and quantum efficiency in future iterations of the device. We expect that it can be optimized to be an on-demand single photon source with near unity quantum efficiency. Such a source could then be used for quantum metrology, quantum computation with propagating photons, or quantum measurements in cases where the shot noise of coherent light needs to be avoided.

author contributions
AG and MH designed the sample and built the initial experiment based on an initial idea by MH and discussions with FP and EDF. AG, SJ, DH, FG and JLT brought up the fabrication process. AG, FG and JLT built the sample. AG performed initial measurements and dataanalysis. FB and RA improved the measurement setup, datataking and calibration routines, took final data and performed final data analysis. JL provided theory support. AG wrote the paper with figures from FB and input from all authors.

I. SAMPLE FABRICATION
The samples were fabricated in a multi-layered process on a Si(500 µm):SiO 2 (500 nm) substrate. We first sputtered an MgO(20 nm) buffer layer and an NbN(80 nm)-MgO(4 nm)-NbN(200 nm) tri-layer, into which we subsequently etched steps defined by a combined optical (OL) and electron-beam lithography (EBL) process. In this way we could simultaneously define small features (Josephson junctions) and large features (coplanar waveguides) for many samples on a 4 inch wafer. After this, a 500 nm thick dielectric layer (Si 3 N 4 ) was deposited conformally by chemical vapor deposition (CVD) and then, following an OL step, etched directionally to leave spacers isolating the tri-layer sidewalls. A counter electrode NbN layer (300 nm) was then deposited and, after another combined OL and EBL step, etched where needed to define junctions. More details on this process and an in-depth characterization of the fabricated junctions and the different superconducting thin films can be found in reference 25. On the device at hand the SQUID is formed by two 0.02 µm 2 Josephson junctions and has a total room temperature resistance of ≈ 230 kΩ. The gap voltage extracted from a measurement of the current-voltage characteristic (not shown) is V gap = 4.8 mV.
The Si 3 N 4 layer also serves as the dielectric in a parallel plate capacitor between the tri-layer and the counterelectrode, which implements the capacitive part of the RC element. It is visible in Fig. 1c and d. The top NbN electrode passes over the grounded flux-bias line (implemented using the tri-layer) and contacts the resistor on the other side.
We fabricated the resistive part of the RC using chromium, which was chosen for its relatively high resistivity allowing us to limit stray capacitances by keeping the resistor short. We implemented it in a final fabrication step directly on the MgO buffer layer in an area where all other layers had been etched away. For reasons detailed in Sec. II D this element alternates thin resistive lines (15 nm) with much thicker cooling pads (100 nm). Both were deposited in the same fabrication step by using angle evaporation and liftoff. For this, an EBL step (ZENON-ZEP520A and AR-300-70) first defined 300 nm wide trenches connecting much larger rectangles (see Fig. 1c and Fig. 5). The thin lines were deposited by evaporating perpendicularly to the wafer plane into the narrow trenches. Afterwards, the direction of evaporation was tilted so that its angle with respect to the plane of the wafer was 35 • and the angle between its projection onto the plane of the wafer and the trenches was 90 • . In this orientation 174 nm of chromium (vertical thickness: ≈ 100 nm) were deposited. Since the resist layer was > 300 nm thick, the Cr from this evaporation step did not reach the bottom of the trenches and was removed during lift-off.

A. Coplanar waveguides
The basic coplanar waveguide (CPW) structures (center conductor and ground planes) were defined in the first optical lithography step and then etched into the tri-layer. We implemented ground-bridges over the CPWs by opening up vias through the Si 3 N 4 above the ground planes on either side of the center conductor and connecting them with the counter-electrode NbN-layer deposited in a later fabrication step (Fig. 6).
In order to determine the required geometries for CPWs of different characteristic impedances we used a Pythonbased 2D boundary element simulation program to compute capacitance and inductance per unit length. We took into account the kinetic inductance of the different NbN thin films, which we estimated from independent measurements 13,25 . It was necessary to correct for the presence of the ground-bridges, since the capacitance to ground of these sections is greatly enhanced. We performed a separate simulation for this specific geometry and then averaged the two types over the length of the line using effective widths for the bridges that take into account their lateral capacitance. This was done by approximating them to microstrip lines across the central conductor and comparing the resulting capacitance to that of a parallel plate capacitor. Table I shows the simulated parameters for different types of CPW elements on the device at hand. )  146  64  1361  26  2  49  1000  90  114  898  33  3  8.5  100  70  140  719  29  5  7.5  100  50  204  512  27  10  5  100   Table I. Some relevant parameters of the CPW geometries used in this work including the characteristic impedance (Z0), the capacitance (C) and inductance (L) per unit length and the percentage of the total inductance due to kinetic inductance (L kin ). S is the width of the center conductor and W is the width of the gap separating the center conductor from the ground plane on either side. The ground bridges have a cross-sectional width of 2 µm and the distance between the centers of two consecutive bridges along the CPW is given by d bridge .

B. Resonator
The resonator is formed by a λ/4 segment of CPW with a characteristic impedance of Z R = 146 Ω (design value) connected on one side to the junction (point 0 in Fig. 7a). When the other end is short-circuited, this CPW segment forms a resonator with resonances at f n = (2n + 1) c 4l where c is the phase velocity in the CPW and l its length. These resonances have a characteristic impedance of Z n = Figure 7. a, Circuit of the on-chip bias tee and beam splitter. Each simple purple line represents a CPW. Two parallel lines stand for coupled CPWs. Characteristic impedances, ports and λ/4 segments are indicated. The SQUID is marked by a cross and is grounded through the RC circuit. b, Simplified circuit at low frequency. The two RF ports are grounded and the SQUID is connected to the DC bias port. c, Simple circuit at resonance frequency. The RF ports are dynamically connected to the SQUID and the DC port is dynamically grounded.
slightly lowers f n and Z n and the fact that the other end of the CPW is not connected to a short circuit but to a beamsplitter (see below) at point C in Fig. 7a gives rise to a finite quality factor. Close to resonance this beamsplitter loads the resonator CPW with an effective impedance of Z C = 12 Ω leading to a quality factor of the resonator of 2n+1 , corresponding to a width of approximately 600 MHz.

C. On-chip beamsplitter and bias tee
A key element of the device is a 4 port on-chip network of quarter-wave CPW segments acting as bias tee and beamsplitter. It connects the junction to the DC port at low frequency. At the operating frequency it connects it equally to two RF ports, while isolating them from each other. Figure 7 gives a schematic representation of the device and Fig. 8 shows an optical micrograph of its central area. The different ports are also indicated in Fig. 1a and b.
Intuitively, the working principle of the device can be understood by considering separately what happens at the resonance frequency f 0 of the λ/4 segments and at zero frequency. In the latter case, the only port connected to the junction is the DC bias port, while the RF ports are grounded (Fig. 7b). In contrast to that, at f 0 the open end of the stub at A transforms to a short at B and to an open again at C. This means that the RF ports and the SQUID are isolated from the DC port and that the impedance seen by the SQUID will be dominated by the RF ports. At the same time, this acts as a filter reflecting noise leaking down the DC measurement setup.
The double lines on either side of point C are capacitively coupled CPWs. While at low frequency they connect the RF lines to ground and disconnect them from the SQUID, on resonance they act as a transformer connecting the SQUID at point C to a transformed RF port impedance of Z coupler = (C C /C T ) 2 Z L ≈ 24 Ω. Here Z L = 50 Ω is the impedance of the RF ports 2 and 3, C C ≈ 160 pF m −1 is the coupling capacitance between the two lines of the coupler and C T ≈ 230 pF m −1 the total capacitance of the line connected to the SQUID (i.e. the sum of its capacitance to ground and to the other line). On resonance, the two couplers together yield a total impedance at point C of Z C = Z coupler /2 = 12 Ω. Over the whole relevant frequency range, we can roughly approximate the effective impedance generated by this bias tee and beamsplitter at point C by a λ/4 transformer with characteristic impedance of Z TL = √ Z C Z L = 24 Ω connected to a single port of impedance Z L = 50 Ω. A quantitative prediction of the device parameters and frequency response was obtained using a Python-based simulation program for linear networks. The simulator has been written in our group and computes the admittance matrix of a network defined by nodes, ports and circuit elements connecting them. The latter can be either lumped elements, simple transmission lines or coupled transmission lines as the ones used in the beam splitter.
Quantities of interest are the isolation between the DC and the RF ports and the filtering of high frequency signals coming down the DC bias line (Fig. 9). We perform this analysis in terms of the amplitudes of the device S-parameters using the definition S ij = 20 log(|V − i /V + j |), where V − i and V + j are the incoming and outgoing complex voltage amplitudes at ports i and j respectively.
We first focus on the reflection of incoming signals at the DC port S DC,DC (solid blue line) in Fig. 9. It is almost unity over a wide frequency range of several GHz around resonance, meaning that high frequency noise coming down the DC bias line is reflected before reaching the sample. The solid green and dotted red lines show the transmission parameters between the RF ports and the DC port. Isolation is good (S DC,i < −20 dB over a span of ≈ 1 GHz) and becomes excellent around resonance. Fig. 10 gives an overview of the different S-parameters concerning the RF ports. Both, the reflection of incoming signals at the ports (S 1,1 , solid blue line), and the transmission between the two ports (S 1,0 and S 0,1 , solid green and dotted red lines) exhibit the behaviour of a bandstop filter with a bandwidth corresponding to the full width at half maximum of the resonant structure ≈ 575 MHz (see also Fig. 2b).

D. Heating and cooling effects in the resistor
The resistor connects to superconducting leads on both sides (SQUID and ground). Therefore electrons cannot be cooled via the leads, but the resistor is deposited directly on the MgO buffer layer in order to facilitate thermalization. Electron phonon coupling in the thick pads between the resistive lines helps to cool the electron temperature of the resistor below the smallest energy scale in the system, which is the charging energy of the capacitor (E C /k B = 2e 2 /(k B C) ≈ 70 mK). This approach is based on the following assumptions: 1. The Joule heating in the resistor can be estimated using the average current given by the RC time. For a sample with 32.1 kΩ and 50 fF the RC time is ≈ 1.6 ns and the resulting current is I = 2e/RC ≈ 0.2 nA. This leads to a Joule heating power of I 2 R ≈ 1.2 fW. We also assume the film temperature to be constant in time. This is the case, because the fluctuations in current on the timescale of RC are averaged out by the heat capacity of the electrons. The electronic heat capacity of a metal of volume V and conduction electron density n at a temperature T is 35 : Here k B ≈ 1.38 × 10 −23 J K −1 is the Boltzmann constant and T F ≈ 5 × 10 4 K is the Fermi temperature. For one conduction electron per atom the electron density n evaluates as n = ρN A /A r , with N A ≈ 6.02 × 10 23 mol −1 the Avogadro constant, ρ the volumetric mass density and A r the relative atomic mass. For Cr, ρ ≈ 7.19 × 10 3 kg m −3 and A r ≈ 52 × 10 −3 kg mol −1 . The volume of the resistor is dominated by the cooling pads and is V ≈ 2000 µm 3 . The energy needed to heat the structure up by ∆T = 1 mK is C e ∆T = 28.3 eV. Comparing this to the Joule heating power yields a rise time: ∆t = C e ∆T /P ≈ 4.5 ms. Fluctuations on the order of ns can, therefore, safely be neglected.
2. The second assumption is that the Kapitza thermal resistance originating from the coupling of the phonon populations across the boundaries between the different layers of our substrate (Si(500 µm):SiO 2 (500 nm):MgO(20 nm)) and the chromium pads is negligible. Following the reasoning of 36 we consider that a thin film cannot have an independent phonon population if its thickness is inferior to the wavelength of the most energetic phonons at a given temperature. Assuming that the Si base of the substrate is well thermalized with the copper sample holder, the combined thickness of the remaining films is d = 620 nm. For T T D = 460 K, where T D is the Debye temperature in Cr 35 , only acoustic phonons are relevant, which have an energy E ph = hv c /λ = k B T . Here v c = 5.9 × 10 3 m s −1 is the speed of sound in chromium 37 . We make the above condition more stringent by requiring that even quarter-wave resonances can be excluded. The minimum temperature needed to excite phonons of λ/4 = d is ≈ 114 mK, a value clearly larger than the base temperature of our dilution refrigerator (15 mK).
3. Finally, we assume that there is no heat diffusion from the resistor into the superconducting contacts. This is justified by the large gap of NbN (T c ≈ 15 K). Thermal electrons are far below the gap and the dissipated power has to be evacuated through electron-phonon coupling in the cooling pads.
The electron temperature of a resistor with volume V , taking into account Joule heating and coupling to a phonon bath of temperature T ph , is given by 36,38 : The electron phonon coupling constant Σ has values around 2 × 10 9 W m −3 K −5 for most metals 38 . To work with a lower bound we consider it to be 0.2 × 10 9 W m −3 K −5 in our case 39 . The resistor has 12 cooling pads connected by 12 thin resistive line of length 20 µm and width 0.3 µm each. In the absence of cooling pads and for T ph = 15 mK the electron temperature of the wire would be ≈ 89 mK. Taking into account the additional volume of the pads the average electron temperature of the entire structure is ≈ 20 mK. While most of the heating occurs in the lines, most of the cooling happens in the pads.
In spite of this acceptable average T e , the local electron temperature in the resistive lines could still be too high. Due to their small volume, electron-phonon coupling in the lines is negligible and they are in the hot electron interaction limit 40 . The temperature profile along the wire in the normalized coordinate x is then given by 38 : The resulting temperature curve is very flat, with a difference between the T e on the borders (electron temperature of the pads) and the maximum in the middle of the wire of less than 1 %.

A. Cross correlation measurements
Following the argumentation presented in 26,27 we combine quadrature measurements from the two different measurement channels to compute the correlation functions of the field emitted by our sample while rejecting contributions arising from amplifier noise. Here, we give a brief overview of the guiding principle behind this approach using a very general noise model to show that unwanted residual noise can be subtracted during data treatment through combinations of simple "On/Off" measurements.

Noise model of the measurement chain
The entire measurement process with different noise sources from the output of the sample to the quadrature measurement is summarized in Fig. 11 and contains the following steps: 1. The cavity mode a on our sample is connected to 3 lines (a DC line and RF lines 0 and 1). Input-output theory 4 can straightforwardly be adapted to this situation by including a unitary scattering matrix S describing the coupling between the ports.b Here γ i is the resonator energy decay rate into line i = DC, 1, 2,b i the outgoing field operator on line i andm i the incoming field operator. As the S matrix is unitary, different output modes commute, [b i , b † j ] = δ ij . The coupling to the DC line is designed to be negligible, γ DC = 0 (see Fig. 9) and on resonance the two RF ports are designed to be symmetrically coupled. Therefore the cavity field leaks out into a modê with rate γ ≈ 2γ 1 = 2γ 2 . Any imbalance and residual loss can be accounted for in the gain and noise of the subsequent amplification step.
Note that, according to equation (9), the outgoing modes of the beam splitter are a linear combination of the incoming modes on the transmission lines which are connected to the RF ports. They are terminated by the cold loads on the isolators which are at a temperature T hf 0 /k B such that the input fieldsm i can be considered in the vacuum state and will drop out in all normally ordered expectation values 26 we consider below. Therefore, ideally only the sample output can produce correlations between different modes. However, we additionally fully remove any residual spurious correlation between these modes, independently of whether they arise from hot input modes or crosstalk at subsequent stages (see Sec. III A 2).
2. In the next step the signal is amplified, which necessarily adds noise. With the power gains g i and the noise modesĥ i , we can write the amplified signals on either chain as 41 : We take the noise on the amplifiers to be independent from their inputs giving [b i ,ĥ † i ] = 0 and b iĥ † i = 0. We suppose that the two noise modes commute, but we do not consider them to be uncorrelated: ĥ † 0ĥ 1 = G (1) ×,amp (t, t + τ ). We also do not assume them to be independent from the noise of incoming modes meaning m † iĥ j = 0 and [m i ,ĥ j ] = 0, for instance due to imperfections of the isolators. No assumptions on the form of the noise cross-correlation are made.
3. Lastly, quadrature measurements are performed on both outputs which introduces additional noise. One can define a complex envelope from the two quadratures, which is proportional to the input field and a noise mode 26 : where ı 2 = −1. These complex envelope operators are defined from the classical outputs of the quadrature measurements and fulfill the equality Ŝ i (t) = S i (t) , where S i (t) = X i (t) + ıP i (t) is a complex number.
The modesk i (t) commute with the amplified fields and all the other modes except for the usual relation [k i (t),k † j (t + τ )] = δ(τ )δ i,j , but they can have non-vanishing correlations.
The noise terms intervening after the first beam splitter can be condensed into one model i per channel i. Even though the noise from mixing will most likely be negligible compared to the amplifier noise in an experimental setup, we keep it here for the sake of completeness.
This operator still respects the commutation relations [l † i (t),l j (t + τ )] = δ(τ )δ i,j and the operatorsŜ i take the form:

Correlations between complex envelopes
The complex amplitudes contain the original fieldâ and can be combined in different ways to obtain the correlation function G (1) . If we were just using one measurement channel, the only option would be: The last equality is derived directly from the noise model described above (see: reference 13, Appendix C). The first term in the last line is the desired correlation function. The second term corresponds to the summed direct correlations of all noise sources on this channel and the last term originates from the commutator [l i (t),l i † (t + τ )] = δ(τ ) and can be seen as the contribution of the vacuum fluctuations of model i . In a realistic experiment its divergence at τ = 0 would be smeared out due to the necessarily finite measurement bandwidth 26 .
A better choice to extract the correlation is: The noise cross-correlation G ×,noise (t, t + τ ) is taken between the noise sources on both channels and is intuitively and experimentally (Fig. 14) much smaller than G (1) i,noise (t, t + τ ). Moreover, the delta function does not emerge, sincê l 0 andl 1 commute. The extraction of G (1) (t, t + τ ) is a simple matter of measuring Γ × (t, t + τ ) with the sample in the "On" and "Off" state and then subtracting the results. For the sake of completeness we have carried the quantum treatment of the field to the very end of the measurement chain. In practice, commutators of fields after the first amplification stage are negligible and a classical treatment is sufficient. Nevertheless, the noise-subtraction was performed in the case of the direct measurement, the much larger noise correlations would make the extraction very sensitive to any change in gain or amplifier noise between the "On" and "Off" measurements.
In a similar manner, one can define several combinations of complex envelopes to obtain the second order correlation function and again some include delta functions and direct noise correlations, while others do not. A choice falling into the latter category is: Again, the noise cross-correlation G ×,noise (t, t + τ ) is the only term that remains when the sample is not active. All the other noise terms are known from the determination of G (1) (t, t + τ ) and can be subtracted to find the correlation G (2) (t, t + τ ). This measurement subtracts any system noise or spurious correlation and contains only one calibration constant, the normalization factor g0g1 2 . In the normalized second order correlation function g (2) this factor cancels as well, so that g (2) is fully self-calibrated without the need to know the gain and noise temperature of the measurement chain.

IV. MEASUREMENT SETUP
A schematic representation of the entire measurement chain is given in Fig. 12. Only elements inside the dilution refrigerator are described in detail. On the right side, the temperatures of the different stages are marked. The setup can be divided into the following parts: Sample: The sample under test as described in the previous sections, with two RF outputs and two inputs for the DC voltage bias and the flux (current) bias.

RF side:
The high frequency branch of the chain includes everything beyond the two RF ports.
Calibration: The blue rectangle highlights a Radiall (R591-763-600) six way switch on each channel. Its output is either connected to the sample or to other terminations used for calibration. More detail is given in section IV A 2.
Filtering and Isolation: Next on each line are bandpass filters (Microtronics BPC50403, passband between 4 GHz and 8 GHz) and three circulators (Raditek RADC-4.0-8.0-Cryo-S21-qWR-M2-b) used as isolators by putting a thermalized 50 Ω load on the third port. The lowest circulator on the left side is instead connected to an attenuated RF input line. This source is usually not active during the experiment and is only used to characterize the sample.
Amplification and measurement: The remaining green section of the chain in Fig. 12 starts at the cold amplifiers mounted on the 4 K stage and ends at the analog to digital converter (ADC) of our measurement computer. It is described in section IV A.

DC side:
Diplexer: A Marki DPXN-M50 diplexer connects the sample DC port to the rest of the voltage bias setup at frequencies below its cross-over frequency (50 MHz) and to a cold 50 Ω termination above.
Transformer and amplifier: The current passing through the sample was measured via the voltage drop it caused over a 20 Ω cold resistor. A CMR-direct low temperature transformer (LTT-h) with a winding ratio of 30 : 1 was used to amplify the voltage across the measurement resistor and up-convert its impedance for better matching with the input impedance of the Celians (EPC1-B) amplifier. The latter has a variable gain (40 dB, 60 dB or 80 dB) and a voltage noise of ≈ 0.7 nV √ Hz −1 . It is connected to the transformer by a Thermocoax cable. The shield of the cable is connected to ground on the base-stage and to the inverted input of the amplifier at 300 K to avoid ground loops. The amplifier is followed by a Measurement Computing (USB-1608GX-2AO) analog to digital converter, with a maximum sampling frequency of 500 kHz. The transformer is housed inside a homemade filter-box including several stages of 10 kΩ resistors, 100 pF feed-through capacitances and homemade common mode choke. The entire system acts as a bandpass filter between 300 Hz and 11 kHz and is optimized for lock-in measurements at a frequency of ≈ 1 kHz.
All resistive elements are NiCr resistors (Sumusu RR1220P-XXX-D) and were tested at 4 K displaying a maximal change in their values of < 2 %. This element was only present during the measurement of the resistance of the RC circuit.
Calibration: The blue rectangle highlights a two position switch (Radiall R 572 F33 000) allowing us to connect the DC measurement to a 5 kΩ resistor (Susumu RR1220P, change at 4 K < 1 %) to ground in order to calibrate the current measurement. Like the preceding component it was not present during most measurements.
Bias box: The orange square represents the base temperature part of a filtered voltage divider housed in a thermalized copper box. Its input is grounded through two parallel 50 Ω resistors, followed by a filtering stage (4 µF capacitor and silver epoxy filter).

K filter:
Finally, the line is filtered at 4 K using a home-made element containing three series 10 kΩ resistors, with two 1 nF capacitances to ground between them.
Flux bias line: The flux bias line is visible at the very left side of the circuit diagram. It consists of semi-rigid cupronickel coaxial cables (Coax Co., Ltd SC-219/50-CN-CN). It is attenuated (20 dB) at 4 K and filtered at base temperature (home-made 50 Ω-matched Eccosorb CRS 124 low-pass filter 13,42 ; cut-off frequency ≈ 1 GHz).

A. RF measurement setup
The main part of the high frequency measurement chain after the calibration stage and the isolators is indicated by the green box in Fig. 12. In order to be able to reject amplifier noise it has two independent channels 13,26,43 starting with the first amplification stage consisting of two Low Noise Factory cryogenic amplifiers (LNF-LNC4-8A) at the 4 K. They work in the band between 4 GHz and 8 GHz and their gain and noise temperature are 44 dB and 2 K respectively. The amplifiers are followed by 3 dB attenuators for protection, thermalization and reduction of standing waves. The lines leading up to room temperature from the cold amplifiers are semi-rigid cupronickel coaxial cables (Coax Co., Ltd SC-219/50-CN-CN), while the lines between the base-stage and the amplifiers are niobium titanium (SC-219/50-NbTi-NbTi).
The first room temperature amplification stage uses Miteq (AMF-5F-04000800-07-10P) amplifiers with a gain of 50 dB and a noise temperature of approximately 50 K.
After some additional filtering (Microtronics BPI 17594, bandpass between 4.25 GHz and 7.75 GHz) the signal is down-converted using Marki (MM1-0312SS) mixers and a Rohde Schwarz SMF100A high frequency source as local oscillator, which is split on a AA-MCS power divider (AAMCS-PWD-2W-2G-18G-10W-Sf) to simultaneously act on both channels.

Frequency down-conversion scheme
The correlation functions measurement presented in Sec. III A relies crucially on having access to the time-resolved measurement record of the sample response in a bandwidth containing the entire resonator. Here we explain our particular sampling process. The signal is initially contained in a frequency band between approximately 4 GHz and 8 GHz, corresponding to the amplifier bandwidth. The digitizer sampling rates allow to sample slices with instantaneous bandwidth of up to 1 GHz, given by the Shannon-Nyquist sampling theorem. This can be done in the first Nyquist band going from 0 Hz to 1 GHz or in the second Nyquist band from 1 GHz to 2 GHz. In both cases, the other band has to be filtered out to avoid aliasing. Figure 13a shows an example of a situation where both positive and negative bands contribute to the measurement. Here the signal (shown in red) is down-converted to the position on the frequency axis indicated by the light red shape. Parts of it lie in either band and, since the ADC does not distinguish between positive and negative frequencies, are simply summed up in the final measurement result. This effect renders data extraction difficult and doubles the amplifier noise. Figure 13b shows a much more favourable situation. The signal is mixed down with an LO frequency of about 7.5 GHz and doesn't fall into the positive Nyquist band anymore. The sampled slice corresponds directly to the original signal between 5.5 GHz and 6.5 GHz, which is centered on the emission frequency the device discussed in this work. In a similar way we can use different local oscillator (LO) frequencies to directly and unambiguously sample the entire bandwidth of our cold amplifiers. Note that this would not be possible, if we were using the first Nyquist band. Then, the effect shown in Fig. 13a could not be avoided for some parts of the signal range.
To obtain a complete map of the PSD in the 4 GHz to 8 GHz range we use the following set of LO frequencies: 3.3 GHz, 3.9 GHz, 4.2 GHz, 4.5 GHz, 7.5 GHz, 7.8 GHz, 8.1 GHz, 8.7 GHz and 9.0 GHz.

Calibration
Even though all correlation function measurements are self calibrated because of the normalization and noise subtraction discussed in section III A, we still have to calibrate the power spectral density emitted by our sample in units of photons. To do so we need to know the gain and noise temperature of our output lines.
The calibration part of the RF setup is indicated by the two blue boxes on the high frequency branch in Fig. 12. On each channel it consists of a 6-port switch (Radiall R591763600) connecting the input of the amplifiers to either the sample or the two calibration references. They are two 50 Ω loads, thermally isolated from the switch by short NbTi coax lines and thermalized at base temperature and at the still respectively. The noise power spectral density coming from a 50 Ω load on a matched line at temperature T is: The total signal on channel i after amplification is: Here, g i and N i are the amplifier gain and the combined effective noise photon number of the amplifier and the cables leading up to it from the 50 Ω load.
At base temperature (T base ≈ 15 mK), S (in) (T ) ≈ 1 2 hf , i.e. approximately independent of temperature, so that accurate knowledge of the temperature of the cold loads is not essential. At the still temperature (T still ≈ 0.9 K), however, the emitted noise is S (in) (T ) ≈ k B T . Therefore, we track the temperature of the loads at still temperature with a dedicated germanium thermometer mounted on their thermalisation copper bracket.
By measuring the signal when switched to the still load (S (out) still ) as well as the signal coming from the base stage load (S (out) base ), we can extract both gain and noise photon number of each measurement chain: Since our measurement is spectrally resolved, this can be done for each frequency point and of course for each local oscillator on both channels, thus fully calibrating the system. Figure 14 shows the result of a typical gain calibration. Different colours represent different local oscillators used for the measurement shown in the main text. The lowest coloured curves in the bottom panel show the cross noise photon number between both channels, which was obtained by combining complex envelopes from both channels similar to the correlation function measurements described above. It is 0 within the uncertainty of our measurement. The black dashed-dotted curve on top of them shows the cross noise measured when the sample with no applied voltage bias (in the "off" mode) served as the low temperature reference instead of the cold 50 Ω loads.

Drift compensation for on-off measurements
This calibration scheme involves commuting the switches fixed to the base stage of our dilution refrigerator and cannot be repeated too often to avoid heating. In between calibrations, the gain of the amplification chain and phase delay may drift, mainly due to temperature fluctuations. To compensate for these drifts we use the fact that in all measurements (PSD, g (1) and g (2) ) presented in this article we periodically turn the bias voltage on and off, with a period of the order of 1 second. We do so in order to be able to remove all unwanted cross-correlation terms not due to the signal emitted by the sample. In the "off" part of the cycle when the sample does not emit photons, we continuously measure the noise of each amplification chain. This noise is dominated by the HEMT noise which is very insensitive to the actual physical temperature of the amplifier and therefore provides a more stable reference than the gain of the full chain containing many active components. We therefore renormalize the gain as follows: In this expression S (out) off is the measured signal when the sample is switched off, averaged over several cycles. Thus, we can account for slow drifts in the gain of the entire chain as illustrated in the left panel of Fig. 15 for a measurement lasting several days.

Additional drift compensation for G (2) measurements
The correlation function measurements performed according to Eqs. (16) and (17) are complex valued with a trivial complex phase given by the difference in phase of the amplitude gain of the two amplification chains. These measurements often last several days and in order to average them correctly, any phase drift must be removed prior to averaging.
To do so, we use the fact that the G (1) measurement has several orders of magnitude lower statistical noise than the G (2) measurement, so that any changes in G (1) averaged over several minutes are due to slow drifts. Fig. 15 shows such a measurement. Each group of blue dots corresponds to a specific time delay τ , with each point representing an averaged measurement. The group of points with the largest magnitude corresponds to τ = 0. We first compensate for magnitude drifts as described above and then shift the phase of the G (1) measurements by δ i + 2π∆f i t by fitting to the points at τ = 0, ±1 ns, ±2 ns for each averaged measurement, indexed by i. δ i accounts for drifts in the gain phase, ∆f i is the difference between demodulation frequency and resonator frequency. After applying these compensations for magnitude and phase drift we obtain the G (1) curves marked by orange points which align on the real axis.

B. Numerical signal processing
Signals are bandpass filtered between 1 GHz and 2 GHz and recorded by an AlazarTech ATS9373 digitizer on 2 channels at a rate of 2 GSamples/s and a resolution of 12 bits. This data is streamed to PC memory and analyzed in real time on 16 CPU cores working in parallel. First we subtract from each channel the data recorded on the other channel convoluted with an appropriate finite impulse response filter (16 taps) in order to reduce crosstalk occurring on the card from approximately −40 dBc to approximately −60 dBc.

Power spectral density measurements
Power spectral density calculations are performed by calculating fast Fourier transforms f i,n of the data in channel i = 0, 1 for each block of data n. We then calculate the power spectral densities PSD i,j = n f * i,n f j,n . PSD i,i is the power spectral density of channel i. PSD 0,1 = PSD * 1,0 is the cross power spectral density. As all signals have low coherence and low dynamic range, we do not apply any windowing in order to get the best signal to noise ratio.

Correlation function measurements
The correlation functions G (1) and G (2) are calculated from complex amplitudes S i of the output signal of channels i = 0, 1 according to Eqs. (16) and (17). We describe here how these envelopes are calculated and how G (1) and G (2) are computed numerically.
At the output of the amplification chain we have voltage signals V i (t) centered around the resonator frequency f 0 . For the g (1) and g (2) we need the complex time-dependent envelope S i (t) of this signal: We down-convert this signal with a local oscillator at f LO .
Where is the bandpass filter matching the 2nd Nyquist band of our digitizer from 1 2∆T to 1 ∆T with ∆T = 500 ps the sampling interval of our digitizer. We choose the local oscillator frequency such that the down converted signal is at the center of our filter, i.e. .
The sampled voltages V i,n = V IF (n∆T ) are = R i,n/2 (−1) n/2 (n even) 1. The real and imaginary part of the envelope are not measured at the same time: S i at even n and S i at odd n, i.e. they are shifted by one sampling interval ∆T . The center of the kernel is at ∆T /2 so that the filtered quadratures are centered at the same time.
2. The bandpass filter has a time response close to a sinus cardinalis, with strong side lobes which can lead to counter-intuitive long-time correlations features. The finite impulse response filter smoothes this response.
3. The filter is slightly wider than our resonator and therefore accepts unnecessary amplifier noise, degrading signal to noise ratio and dramatically increasing data acquisition times. The finite impulse response filter closely matches the resonator bandwidth.
By combining the quadratures we then obtain the desired complex envelopes S i,m for channels i = 0, 1 sampled at m + 1 2 ∆T . In order to calculate the Γ (1) × and Γ (2) × according to Eqs. (16) and (17) we first calculate the fast Fourier transform of a block of N data points of S 0,m , S 1,m and S 0,m S 1,m padded with an equal amount of 0 (necessary to avoid wraparound), which we call, respectively F 0,n , F 1,n and F ×,n . The Fourier transform F (1) of Γ (1) × is then F (1) i = F * 0,n F 1,n and the Fourier transform F (2) n of Γ (2) × is F (2) n = F ×,−n F ×,n . These Fourier transforms are then averaged over many blocks of size N . The inverse Fourier transform is performed in post processing and a scaling factor (N − |τ /2∆T |) −1 is applied in order to account for reduced overlap in the finite convolution product with increasing τ .

V. ERROR ANALYSIS
The number of averages performed for each measurement presented in Fig. 2, 3, 4 of the main text is given in table II. We estimated the statistical error of our data by saving averages over n b blocks of n avg samples on each channel. The standard deviation of the final result is then σ = σ b / √ n b − 1, where σ b is the standard deviation of the results for each block. We found the error of the first order correlation function measurement to be negligible (e.g. for G (1) (0, 0) = 94.82 MHz in Fig. 3a, b the statistical error is ±3σ = ±0.05 MHz). This also holds for the power spectral density and emission rate data shown in Fig. 2 and Fig.3a. For the second order correlation function g (2) (0, τ ) the error is significant because the number of averages necessary for a given signal to noise ratio scales exponentially with the order of the measured correlation function (see reference 26 and reference 13 Ch.3). In Fig. 16 we show g (2) (0) for each curve in Fig. 3b and Fig. 3c of the main text together with the extracted error bars (±3σ). The point at τ = 0 was specifically chosen because the measured signal is weakest and the amplifier noise correlation highest at this value making the error bars an upper bound for the other data points.  We have performed a complementary measurement of the dark region in Fig. 2c of the main text by measuring the flux-and voltage bias dependent transmission through the two RF ports of the sample with a vector network analyzer (VNA). This was done by sending the stimulus to the sample via the circulator of one of the amplifier chains (visible on the left side of Fig. 12) and then measuring the amplified response on the other channel.
As discussed in section II, we expect the two RF ports to be decoupled on resonance, where the last λ/4 element between the node at point C and the SQUID acts as an open stub 44 . Figure 17a gives a schematic representation. It only shows the parts of the on-chip beam splitter and bias tee relevant to this measurement. The horizontal section corresponds to the two arms going down at 45 • angles in Fig. 7 coupling point C to the RF ports. At the resonance frequency of the device the stub decouples the two RF ports leading to an anti-resonance dip in the transmission measurement (Fig. 10).
From Fig. 17a we can see that the SQUID is at the open end of the resonator, where the voltage is at its maximum, Measured frequency f (GHz) Figure 18. Frequency of the anti-resonance in a VNA transmission measurement between the two RF ports of the sample as a function of flux in units of flux quanta Φ0 and voltage bias. The frequency becomes flux-tunable where the SQUID is on its current branch and has a finite contribution to the inductance of the stub shown in Fig. 17. The dotted curve is a function | cos(πΦ/Φ0)|, rescaled so its apex matches the top of this region. Outside of the boundary the SQUID is on the voltage branch and the frequency of the anti-resonance does not depend on flux.
while the current is zero. This means that it is situated at the point of maximum impedance on resonance. This resonant circuit can be modelled by a series LC-circuit, since it produces a dip in the measured signal. The inductance coming from the SQUID shunts the effective capacitance to ground as shown in Fig. 17b. Here, we neglect the impedance of the RC circuit in series with the SQUID because we suppose the junction impedance to dominate. The anti-resonance in the measurement occurs at the frequency where the total impedance to ground of this effective circuit is zero: Here, C eff and L eff are the effective capacitance and inductance of the resonator and L SQUID is the tunable inductance of the SQUID, which (in the balanced case) is given by In this expression I c is the critical current of the SQUID, Φ is the external applied flux and Φ 0 is the flux quantum.   18 shows the center frequency of the measured anti-resonance dip as a function of flux and voltage bias applied to the SQUID. It clearly shows that the frequency is only flux-tunable in the region where no photon emission occurs in Fig. 2c of the main text. Outside of this region the frequency becomes constant. This can be understood by considering that the series resistance of the RC circuit tilts the load line of the voltage bias at the junction. For bias voltages such that V b < RI c (Φ) the SQUID stays on its current branch undergoing Bloch oscillations. Then it contributes a finite inductance to the resonant structure making it flux-tunable.
Only when the voltage is sufficiently high, or I c sufficiently suppressed by the flux bias, a voltage drops over the SQUID. Then its effective inductance becomes infinite and does not contribute to the frequency of the anti-resonance anymore. The boundary between the two regions is well described by a function | cos(πΦ/Φ 0 )| (dotted line in Fig. 18, rescaled to match the apex of the bright region). This line is also drawn on Fig. 2c for comparison. It delimits the region without photon emission corroborating the interpretation given in the main text.

VII. RESISTANCE MEASUREMENT
The current through the sample as a function of bias voltage was independently measured (see Sec. IV). The result is given in Fig. 19. Each curve was taken at a different flux bias, going from nearly full frustration (lowest current) to minimum frustration (largest current). The slope of the current branch directly gives the resistance in series with the SQUID, which is entirely dominated by the resistor of the RC circuit. From these curves we extract a value of 32.1 kΩ. Moreover, the measured curves show good agreement with the data presented in Sec. VI.

VIII. EXTRACTION OF SYSTEM PARAMETERS USING P(E)-THEORY
The rate at which Cooper pairs inelastically tunnel through a voltage biased Josephson junction (or equivalently a voltage biased SQUID acting as an effective junction with an adjustable critical current) can be computed with Fermi's Golden Rule considering the junction as a perturbation to the modes of its electromagnetic environment. The result depends on correlations between the phase fluctuations at the junction, which are related to the real part of the environmental impedance seen by the junction via the fluctuation-dissipation theorem. An overview of this calculation (called P (E)-theory) is given in reference 17.
The expression for the Cooper pair tunneling rate into the direction of the voltage bias is given in Eq. 2 of the main text. It depends on the function P (hν J ) giving the probability density for the environment to absorb an energy hν J from a Cooper pair tunneling through the junction. Here, ν J = 2eV b /h is the energy given to the Cooper pair by the voltage bias expressed in frequency units. To simplify the notation we rewrite P (ν J ) = hP (hν J ). This function obeys the following normalization and detailed balance relations: P (−ν J ) = e −βhνJ P (ν J ) Here, β = 1/(k B T ), with k B the Boltzmann constant and T the temperature of the electromagnetic environment. The bias range V from ν J = 0 to 30 GHz of our measurements covers the dominant processes at ν J = E C /h ≈ 1.5 GHz, E C /h + f 0 ≈ 7.5 GHz and E C /h + f 1 ≈ 19.5 GHz, so that we can consider P (ν J ) normalized over the bias range V with good approximation.
The second equation signifies that P (−ν J ) vanishes at zero temperature. This happens, because at negative energies the Cooper pairs have to tunnel against the voltage bias, thus drawing their energy from the thermal excitations of the electromagnetic environment. At finite temperatures P (ν J ) is given by the integral equation 13,17,45 : This form, called the Minnhagen equation, depends only on the real part of the frequency dependent environmental impedance Z R (f ), the temperature, and the superconducting resistance quantum R Q = h/(4e 2 ) ≈ 6.5 kΩ.
It can be shown 13,19 that the rate density of photon emission into the electromagnetic environment at frequency f due to forward tunneling Cooper pairs can be expressed in terms of the critical current of the junction I c , the impedance Z R (f ) and the P(E)-function as: This quantity is proportional to the power spectral density emitted by the sample with a conversion factor given by the photon energy hf and can be directly measured (see Fig. 2a of the main text). By integrating both sides of the above equation over the voltage bias and using the normalization property from equation (33) we find: The integral is performed over the data shown in Fig. 2a of the main text. Evaluating equation (36) at a voltage bias ν J + f results in: This can be combined with equation (37) to yield the expression of P (ν J ) where σ P is a weight function, which we chose large where γ is large. Fig. 20 shows the extracted P (hν J )-function of our experiment. It displays a prominent peak around ν J = E C ≈ 1.5 GHz corresponding to forward tunneling of Cooper pairs without photon emission while only charging the capacitance of the RC circuit. This peak does not appear in Fig. 2c, which only measures the photons emitted around f 0 . The next two peaks at 7.5 GHz and 13.5 GHz are due to the hν J = E C + hf 0 and hν J = E C + 2hf 0 processes described in the main text. The second one of these two peaks is linked to a higher order process (involving two photons at f 0 ) which occurs with much smaller probability than the first one 19 . The last peak at 19.5 GHz has again a higher amplitude, but does not show up in Fig. 2c. It comes from the Cooper pair current linked to the resonance condition hν J = E C + hf 1 . The photons emitted by this process are outside of the measurement bandwidth of our setup. Energy ν E = E /h (GHz)  Figure 20. Extracted P (hνJ) for forward tunneling Cooper pairs as a function of the energy difference a CP has to loose. Negative values on the abscissa correspond to thermally activated tunneling against the voltage bias. The blue dots are extracted from the data with no assumption on the environmental impedance and the orange line is a fit of the data based on our impedance model (see Fig. 1a).
Note that the photon emission rate density measured at a frequency f and voltage bias ν J such that ν J +f < f , gives us access to P (ν J ) for negative arguments. Since this is related to thermally induced tunneling against the voltage bias, we can use equation (34) to compute β and thus the effective temperature of the electromagnetic environment: where the weight function σ β is chosen largest where we get the best signal to noise ratio for β. This temperature together with the P (ν J ) function and the measured photon rate density enables us to extract the critical current.
Integrating the Minnhagen equation (equation (35)) and using the normalization property equation (33) we get with equation (37) can be regrouped and integrated over frequency to give The weight function σ f must be entirely contained in the measurable bandwidth in order to be able to evaluate these integrals from our measurement results. In our case this bandwidth is 4 GHz and is not large enough to fully contain the main peak of the P (ν J ) function. We therefore chose σ V to contain positive and negative values, so that σ f , the convolution product of σ V and P , remains approximately limited to our measurement bandwidth. Through identification of these two equations we then obtain an expression for the critical current depending only on known quantities and fundamental constants: Finally, we can extract the impedance seen by the voltage biased SQUID (shown in Fig. 2c of the main text) using equation (37): From the photon rate density shown in Fig. 2a of the main article we find an effective temperature of 20.9 mK and a critical current of 0.82 nA.
To verify our approach we now introduce a realistic circuit model based on the sketch given in Fig. 1a of the main text. It consists of the RC element with resistance R = 32.1 kΩ (measured independently) and capacitance C (with a parasitic inductance L p in series with the capacitance; see Fig 1) as well as a stepped transmission line resonator with characteristic impedances Z 0 and Z 1 .
We perform a fit of the photon rate density using this impedance model with C, L p , Z 0 , Z 1 , T and I c as free parameters. We find the values Z 0 = 110 Ω, Z 1 = 22 Ω, C = 56.7 fF, L p = 53 pH, I c = 0.85 nA and T = 21 mK in good agreement with our earlier analysis. The impedance Z 0 of the transmission line section close to the SQUID is lower than its design value, which could be explained by the influence of the SQUID capacitance. The shape of Z R (f ) given by our circuit model for these values is plotted in Fig. 2c of the main text and reproduces the extracted curve up to impedance modulations likely due to spurious reflections in our output lines, which are not part of our model. The P (hν J )-function found by this fit is plotted in Fig. 20.