Microwave photon detection at parametric criticality

The detection of microwave fields at single-photon power levels is a much sought-after technology, with practical applications in nanoelectronics and quantum information science. Here we demonstrate a simple yet powerful criticality-enhanced method of microwave photon detection by operating a magnetic-field tunable Kerr Josephson parametric amplifier near a first-order quantum phase transition. We obtain a 73% efficiency and a dark-count rate of 167 kHz, corresponding to a responsivity of $1.3 \times 10^{17}~\mathrm{W}^{-1}$ and noise-equivalent power of 3.28 zW/$\sqrt{\rm Hz}$. We verify the single-photon operation by extracting the Poissonian statistics of a coherent probe signal.


I. INTRODUCTION
Similarly to the way in which entanglement is regarded as a resource for enhancing metrological precision, systems displaying critical behavior -with small variations of physical parameters leading to macroscopic changes [1,2] -can be employed for sensitive detection.For example, superconductors near the critical temperature can become metallic even if only tiny amounts of energy are absorbed, and the ensuing change in electrical resistance can be easily measured.This principle underlies the design of transition-edge sensors widely used in astronomy [3].However, classical phase transitions respond only to changes in thermodynamic variables, which restricts the range of applications.In contrast, quantum phase transitions [4] are sensitive to variations in the parameters of a Hamiltonian.In recent years a significant number of theoretical proposals have been advanced toward the investigation of quantum phase transitions for fundamental metrological research [5][6][7][8], for the preparation of quantum resources such as highly-squeezed spin states and extractable work in quantum batteries [9,10], and for applications in the sensing of electromagnetic or mechanical degrees of freedom [11][12][13][14][15][16][17].
The parametrically driven quantum oscillator with Kerr nonlinearity and dissipation encapsulates in a quintessential way the physics of quantum phase transitions.Several paradigmatic models map to it: the Lipkin-Meshkov-Glick model via a Holstein-Primakoff transformation [18], the quantum Rabi model via a Schrieffer-Wolff transformation [19], the Dicke model in the limit of suppressed spin-fluctuations [20], and the PT symmetry-breaking model in the mean-field description [21].However, investigating these models experimentally is often challenging.For example, resolving the structure of the ground state by extracting virtual photons in the Rabi model requires sophisticated protocols [22], while for the parametrically driven harmonic oscillator this is straightforward.Moreover, using finite-component criticality instead of a many-body system presents multiple advantages such as simplicity of the device and Hamiltonian description, avoidance of non-local operations, and direct applicability of quantum control methods.
Here we employ a Kerr Josephson parametric amplifier (K-JPA) realized as a quarter-wavelength superconduct-ing stripe terminated by a SQUID to map in detail the phase diagram and to detect single-photon microwave radiation.The effective electrical length can be tuned fast by using an external magnetic field, allowing us to parametrically pump the device at a frequency close to double of the resonant frequency, inducing 3-wave mixing.Similar devices have been studied in other contexts: amplification, fundamental studies such as the dynamical Casimir effect [23], generation of entanglement and cluster states [24,25], qubit-readout protocols [26], and unconventional forms of quantum computing by using qubits encoded in continuous variables [27][28][29][30][31].However, while the behavior of this system below the threshold is well understood and verified experimentally, the physics near and above the threshold is still uncharted territory from the perspective of applications in sensing.In our experiments, we operate the device as a criticality-enhanced detector by exploiting the crossing of a first-order phase transition triggered by incident quanta.
Our experiments demonstrate conclusively a new concept for a microwave single-photon detector (SPD) -a much sought-after device in mesoscopic physics and nanoelectronics [32].In the optical domain, the energy carried by a photon is as low as 10 −19 J (∼ 1 eV), yet it is still possible to realize single-photon detectors (SPDs) based on photomultiplier tubes, avalanche photodiodes, quantum dots, and superconducting materials including transition edge sensors and nanowires [33][34][35].However, for microwave signals with four-to-five orders of magnitude lower energy, single photon detection becomes extremely challenging [36].Indeed, the standard photodiode principle cannot be applied due to the low photonto-electron conversion efficiency at these frequencies as well as due to the lack of semiconductors with very small gaps.Yet, progress has been made by employing semiconducting nanowire quantum dots coupled resonantly to a cavity, albeit the efficiency obtained so far is only 6% [37].Also recent efforts in hybrid superconductor-graphene-superconductor bolometers has lead to energy resolution as low as 30 GHz [38,39].Other modern microwave single-photon detection schemes and proposals are utilizing the switching of a current-biased Josephson junction from zero to a finite voltage state through absorbing a photon [40][41][42][43][44], the cross-Kerr interaction in a Josephson array embedded in a cavity [45], and various mappings of the photon onto the state of a superconducting qubit [46][47][48][49][50].These detection schemes may also need to be heralded.While some of these schemes have been successfully demonstrated in the laboratory, the design and fabrication is often complicated, with considerable operational overhead and less reliable parameters (e.g.frequency and decoherence times can change over hours for superconducting qubits).In contrast, our detector is simple to fabricate, requires only an easy, bolometer-style operation protocol, and it is non-heralded.Since the design does not use qubit-type elements, the detector is not affected by quasiparticle poissoning and can be operated in a stable way over large periods of time.
The results presented here validate previous theoretical work on phase transitions in similar systems, and in particular the description utilizing a slow variable and its effective potential [51][52][53][54][55][56].This proof-of-concept detector for criticality-enhanced sensing demonstrates that experimental control over this effective potential is possible and may lead to a new class of devices and applications.Indeed, in mesoscopic physics, the formation of effective potentials, typically by circuit engineering, is a key strategy for realizing devices such as superconducting qubits.
The realization of an efficient and robust single microwave photon detector will open up new opportunities in quantum information processing, metrology, and sensing.For example, the SPD is the key ingredient in oneway quantum computing [57] and in quantum illumination with microwave photons [58]; it can be used for qubit readout [16,59] and it would enable dark-matter axion detection with sensitivity below the standard quantum limit [60][61][62][63].Our criticality-enhanced concept lends itself to further theoretical and experimental extensionsfor example machine-learning optimization of the pulse sequence and shapes, investigation of other types of perturbations (e.g.small magnetic fields), and research into new types of transitions obtained by coupling other degrees of freedom (qubits, another oscillator, etc.) to the slow variable.
The paper is organized as follows.We first introduce the device and obtain the quantum phase diagram under parametric operation.We identify the operational point near the border of the phase transition.Further, we describe the pulse sequence, we characterize the device as a detector, and we verify the Poissonian statistics of photons in a probe field.

II. DESCRIPTION OF THE DEVICE
The K-JPA device consists of a λ/4 coplanar waveguide 50 Ω resonator terminated with a SQUID loop, see Fig. 1.The tunnel junctions and wiring layers were fabricated at VTT Technical Research Centre of Finland by Nb and Al sputtering on a high-resistivity silicon substrate [64].The value of the critical current for the Nb/Al-Al 2 O 3 /Nb 1×1 µm 2 JJs is I c ∼ 8 µA.The mea-FIG. 1.
Schematic of the sample and the lowtemperature configuration of the experimental setup.The flux-modulated K-JPA is pumped at twice the frequency of the probe signal.The probe signal propagates through the circulator into the parametric cavity.Bias-tees allow to combine the RF-pump and the DC-control of the detuning.Inset: optical image of the sample with a closeup of the SQUID and feed line.
surements were done at 20 mK in a BlueFors dry dilution cryostat, where the device was mounted in a sample holder and protected from external magnetic noise with a Cryoperm shield, see Supplementary Information Section I.The frequency of the resonator can be tuned by applying a magnetic field, ω 0 (Φ DC ) ≂ ω λ/4 /[1 + χ(Φ DC )] [65,66], where ω λ/4 /(2π) = 6.179GHz is the bare resonance frequency in the absence of the SQUID, Φ DC is the applied magnetic flux, and χ(Φ DC ) is the participation ratio of the Josephson inductance.The range of frequencies span from a maximum of 6.117 GHz at zero field to below 5.5 GHz.At the chosen operational point we have Φ/Φ 0 ≈ 0.362, which results in χ 0 = 0.0227 and ω 0 /2π = 6.042GHz.The K-JPA is pumped with a microwave tone applied also via the SQUID at a frequency ω P close to 2ω 0 , resulting in three-wave mixing [67,68].The DC flux bias and the RF pump share a common on-chip flux line, and the signals are combined by an external bias-tee.The inner conductor of the line is separated from the common ground of the sample and cryostat (see Fig. 1), in order to avoid additional heating by ground currents.The input signals are applied via a circulator with a 4 − 8 GHz frequency band.The pulses are phase locked and controlled by trigger events with high accuracy.The readout is done by using a standard heterodyne scheme, employing a local oscillator tone at LO = 6.028GHz which is modulated at an intermediate frequency IF = 14 MHz.
In a frame rotating at half the pump frequency and in the rotating wave approximation, the Hamiltonian of the system reads are the dimensionless quadratures of the cavity field a, α is the pump amplitude, and K is the Kerr nonlinearity.The detuning ∆ = ω 0 − ω P /2 is defined with respect to half the pump frequency ω P .The coupling constant of the probe signal into the device is denoted by κ and the internal dissipation is characterized by a rate γ.From characterization measurements of the spectrum and gain coefficient we extract κ/2π = 4.44 MHz, γ/2π = 2.30 MHz, and we get K/(2π) ≈ −0.2 kHz by fitting with a simulation of the quantum Langevin equation and by calculations based on circuit parameters.The next-order term is sextic, and from circuit analysis we obtain a value of a few orders less than the Kerr constant, Λ ∼ 10 −6 |K|, therefore we will neglect it.Overall, the device is remarkably robust under cryogenic cycling and aging: over a period of three years, we have not observed significant changes in bandwidth, gain, or flux-dependency of the spectrum.
For certain values of the parameters (∆, α) the system with the Hamiltonian Eq. ( 1), may become unstable, resulting in phase transitions to oscillatory states, as described in the next section.

III. PHASE TRANSITION DIAGRAM
The existence of phase transitions between different phases (or states) of a parametrically pumped nonlinear oscillator has been predicted theoretically [51][52][53][54][55][56].In order to analyze the phase transitions and map the full phase diagram, we capture the quadrature voltage amplitudes X and Y by heterodyning the amplified output of the device, which we express hereafter in arbitrary units.For every value of the pump amplitude α and the flux bias Φ/Φ 0 we acquire the raw X, Y data over a relatively large time duration of 0.1 s.We sweep α and Φ/Φ 0 in the range from 0.352 to 0.37 (corresponding to detunings ∆/(2π) from -6 MHz to 6 MHz).Note that increasing the flux causes a decrease of the resonator frequency ω 0 and therefore a decrease of detuning, while the pump/readout frequency is fixed in experiment.Since the phase of the complex signal X + iY is arbitrary, we rotate it for convenient visualization such that the stretching and displacement are oriented along the X axis.We calculate the first and the second order statistical moments of the X, Y distributions and characterize the states in steadystate measurements with respect to the average values ⟨X⟩ , ⟨Y ⟩ and the variances In Table I we show the resulting classification of these states: vacuum, squeezed vacuum, double-displaced state and coherent oscillation (also called "period-2 vibrations" 51).
The phase diagram obtained in this way is presented in Fig. 2. The phase diagram agrees very well with the theoretical predictions based on identifying a slow variable Q and its associated effective potential U(Q) (see Appendix B and Supplementary Information Section II).By changing the pump amplitude and the flux biasing we can control the state of the K-JPA near the critical threshold as follows: • A second-order transition exists at negative detunings ∆ < 0. The system changes smoothly from the squeezed vacuum to the double-displaced state, see Fig. 2a, point 3 .
• A first-order phase transition is present at positive detunings ∆ > 0 (see for example point 5 in Fig. 2a); here the mean value has a jump 0 → ⟨X⟩.This transition border is definitely sharp and suitable for defining the operational point for photon detection.
For the detection of weak signals approaching the level of single photons, the first-order transition border at positive detunings ∆ = ω r − ω P /2 and critical pump amplitude value α c (∆) = ∆ 2 + (κ + γ) 2 /4 is preferable owing to its sharpness and reliable discrimination of switching events due to the large average quadrature ⟨X⟩.Three local minima exist in the effective potential when approaching from the right of the transition borderline, similar to the Landau theory of first-order transitions (see also Supplementary Information Section II).In contrast, the second-order transition appearing at negative detunings is characterized by a smooth change from the zero state to states with relatively small positive or negative quadratures.
To observe transitions from squeezed vacuum to coherent oscillation experimentally, we monitor the in-phase quadrature X for small powers of the incoming probe field, and observe sudden bursts to large values.Several operational points were probed at different positive detunings; the suitable ones should have a reasonable small switching probability in the absence of probe excitation (dark count rate) and also should present a measurable change in these probabilities when the probe field is turned on.Finally, we decided for an operational point close to the triple-state point 6 , at a flux Φ/Φ 0 ≈ 0.3618 (positive detuning ∆/(2π) = 0.7 MHz)   and pump amplitude α/(κ + γ) ≈ 0.51, marked with a red cross in Fig. 2.
Next, the detection protocol can be devised as follows.In the central well of the effective potential we prepare the state characterized by zero mean amplitude ⟨X⟩ = 0 (squeezed vacuum state).In the absence of a probe signal, the state of device can change by activation over the barrier due to quantum and thermal fluctuations [51,69]; these events are called dark counts.If a microwave probe pulse is applied, the effective potential gets tilted and the activation is easier (see Appendix D and Supplementary Information Section II).

IV. PHOTON DETECTION
We can now demonstrate single-photon detection at the first-order phase transition.In this case, the perturbation triggering the switching is provided by the absorbed photon.
a. Pulse sequence.The protocol uses a specific sequence of pulses for the pump, the probe, and for the readout window, as presented schematically in Fig. 3.The detection procedure starts with the unpumped device in the ground (vacuum) state, after which we bring it at the operational point (α/(κ + γ) ≂ 0.51, Φ/Φ 0 ≂ 0.3618, ∆/(2π) = 0.7MHz) by ramping up the pump.We send a rectangular-shaped calibrated probe pulse containing n photons on average to the input port of the K-JPA and we record the output signal.During the readout window of duration τ = 1.5 µs the recorded trace is averaged to define the detection result R τ = (1/τ ) τ 0 dt|X(t) + iY (t)|.In the final step of the measurement we switch off the pump, creating the dead time necessary for the cavity to relax back to the ground state.This is needed because the strong oscillations, once established, will not relax by themselves and the detector will not be able to register another photon.This dead time is set to 3 µs, from the falling edge of pump pulse to the rising edge of next pump pulse, sufficient to bring back the parametric cavity to its ground state.The time latency, dead time, and the readout window value define the minimum time before the next pulse sequence, resulting in a maximum repetition rate of 0.2 Mcounts/s, see Fig. 3.The calibration of the mean photon number n is based on the determination of the gain coefficient in the output line and on a reflection measurement of the attenuation in the input line, (see Appendix A and Supplementary Information Section I).We collect sufficient quadrature data (N = 2 • 10 4 pulses) at every value of n to obtain reliable statistics and to construct histograms normalized over the number of pulses N in the sequence, see Fig. 4.
From these probability distributions one can see that the detection is photon non-resolvable, due to the overlap of distributions for different n > 1.Thus, the device is inherently a threshold detector.
b. SPD basic characterization.A threshold detector yields one bit of information corresponding to either of the two states "switched/click" or "non-switched/nonclick", depending on weather the response of the device is below or above a certain threshold set by the experimentalist [35].Here we suggestively denote these two states as 0 and 1 + , defined by R τ < R th and R τ > R th respectively.Here R th is the response threshold.Ideally, the detector should react to the presence of a probe coherent field at the input by switching with probability where η is the single-photon efficiency (see also Appendix C and Supplementary Information Section II).In reality, the detector may switch also when no probe signal is present, producing a false positive (dark count); we denote the probability of these events by p dark .As a result, the probabilities of recording 0 or 1 + are not the ideal p 1 + and p 0 .Instead, we have the true-positive probability P 1 + and its complement P 0 with a structure obtained by the composition rules for probabilities P 1 + ≡ p 1 + ∨dark = p 1 + + p dark − p 1 + p dark and P 0 ≡ p 0∧¬dark = (1 − p dark )p 0 , with normalization From these definitions we can relate the efficiency to the probabilities P 1 + and p dark that are measured directly For our detector, the efficiency depends weakly on n (see Supplementary Information Section II); this effect becomes visible at large values of n.We can include this dependence by the replacement η → η ϵ (n), with η ϵ (n) defined phenomenologically as η ϵ (n) = η 1 1+ϵ(n−1) where η is the efficiency at n ≃ 1 and ϵ ≪ 1 is a subunit constant.
We start by identifying a suitable operational point and then optimizing the threshold value R th .This is done for different pulse durations τ p = {0.25;0.5; 1; 2}µs, all with n = 1 photons.The protocol [33,40] consists of the following sequence of measurements (see also Fig. 3).Firstly, we do not send any probe pulses (PROBE-OFF experiment), obtaining the false positive ratio FPR ≡ p dark ; secondly, we send N probe pulses containing n photons each (PROBE-ON experiment) obtaining the true positive ratio TPR ≡ P 1 + .Then we build the receiver operation characteristic (ROC) with floating threshold value R th , see Fig. 5 a).We choose the operational point parameters (α, ∆) such that the area-under-curve (AOC) is maximal.Next, we find the optimal value of the threshold by constructing the TPR(1 − FPR) = P 1 + (1 − p dark ) curves (see Fig. 5b).The threshold R th ≈ 2.45 corresponds to the maximum of the τ p = 1µs trace and it produces a near-maximum for the others as well.This threshold value coincides to the R th producing the farthest deviation of the traces from the diagonal reference line.
Next, by using Eq. ( 22) we can extract the efficiency η.The results, with the operational point and threshold optimized for each value of τ p = {0.25;0.5; 1; 2} µs are presented in Fig. 5c).We can see that the detector can be thus calibrated to photon detection with similar efficiencies η and dark count probabilities p dark for pulses of different durations.
A more detailed analysis is shown in Fig. 5 d).Here we present the measured P 1 + and p dark for a wide span of n photons together with numerical simulations using a Fokker-Planck equation describing the dynamics of switching from the effective potential shown in Fig. 3 (see Appendix D and Supplementary Information Section II).We extract the ideal probability p 1 + = (P 1 + −p dark )/(1−p dark ) and we fit it with 1−exp(−ηn) and with 1 − exp(−η ϵ n) (here ϵ = 0.1), obtaining η = 0.73.Furthermore, we observe that at low n ≪ 1 a coherent state yields approximately the same probabilities as a state (1 − n)|0⟩⟨0| + n|1⟩⟨1|, which represents a source emitting single photons randomly with probability n.From the POVM formalism we obtain in both cases p 0 | n≪1 ≈ 1 − nη and p 1 + | n≪1 ≈ nη.This relation can be readily verified from Fig. 5d), and it has the meaning that the detection probability is the product between the emission probability and the single-photon detection efficiency η.In other words, from Fig. 5d) we can see that the dark count probability is the intercept of P 1 + with the n = 0 vertical axis, while the efficiency is the slope of p 1 + at the same point, η = (dp 1 + /dn)| n=0 .
c. Noise equivalent power.Seen as a power sensor, the sensitivity to the microwave field present at the input can characterized by the power responsivity (susceptibility) at low photon numbers We obtain χ P = 1.3 × 10 17 W −1 .The noise equivalent power (NEP) (see also Appendix C) can be obtained from [70][71][72] For Γ dark = 0.167 MHz and ω p /(2π) = 6.042GHz we get NEP ϵ = Observation of Poissonian statistics.In general, the effect of a subunit efficiency can be understood as a reduction of the number of detected particles by η (see Appendix C).This can be modeled by the use of an attenuator with transmission coefficient η, which reduces the average number of particles, followed by detection with unit efficiency [73][74][75].Moreover, it is known that the Poissonian character of a coherent state is preserved
Pulse sequence for the probe pulse turned OFF or ON, with the illustrations on the right showing the effective potential U representation as a function of the slow quadrature Q (see Appendix B and Supplementary Information Section II).Each measurement starts by sending a pump pulse at ωP, with the system initiated in the ground state (i), in order to bring it at the operational point (ii).After a small delay time the probe pulse is applied, resulting in excitations out of the central well of the effective potential (iii).The detection session ends by switching off the pump in order to bring the cavity back to the ground state (i).The measurement yields a detection result Rτ defined as the time-average of |X + iY | during the readout window width τ .The role of the PROBE-OFF experiments is to measure the FPR contributed by dark counts and technical noise in the output part of the setup.The optimal threshold value, responsible for the "click-no click" decision, is eventually set up according to ROC and TPR(1-FPR) optimization, see Fig. 5a)b).Then, in the PROBE-ON experiments, the decision about detection events is taken by comparison of Rτ with the chosen threshold value R th .
under attenuation [76], as it is apparent from the formulas p 1 + = 1 − exp(−ηn) and p 0 = exp(−ηn).Although apparently innocuous, it is not at all guaranteed that this ideal-detector exponential dependence reproduces correctly the response of the real detector, even when corrected for the dark counts.
Here we demonstrate the Poissonian statistics of an input coherent field over a wide range of values of n and thresholds R th by comparing the experimentallymeasured P 0 with the values obtained by factoring the measured probability of no dark events and the idealdetector no-event probability exp(−η ϵ n), in other words Here η ϵ is extracted at each value of the threshold from the measured probability at just one point, n = 1, and ϵ = 0.1.Fig. 6 shows this comparison for the case τ p = 1 µs over the entire range n ∈ {0, 10}.In view of the high degree of agreement obtained, we conclude that the results verify convincingly the Poissonian statistics of the incoming photons.
We can further corroborate the single-photon operation by using two binomial test whereby the underlying probabilities entering in the binomial distribution are obtained from the Poissonian statistics, see Supplementary Information Section IV.

V. CONCLUSIONS
We have realized a detector of microwave photons by employing the sensitivity to perturbations of a Kerr Josephson parametric amplifier operated near the critical threshold of a first-order phase transition.We have mapped the full phase diagram and showed that the system is described by an effective potential in the slow quadrature.The operational point is chosen such that the system is in a metastable state near a first-order transition, and consequently even single-photon microwave pulses can trigger the system to switch from the squeezed state to an oscillatory state.These oscillations are further amplified to values above the thermal noise and eventually detected using a homodyne scheme at room temperature.Our results include measurement and optimization of threshold via the ROC and the TPR(1-FPR) curve, and we further extract the efficiency and the dark count rate.We also show how this detector can be used to observe the Poissonian statistics of photons in a coherent probe field.Our results demonstrates that the effective potential can be controlled and the phase transition can be harnessed to producing a useful device.
For future experiments, the dark count rate may be further lowered by reducing technical noise and by increasing the sharpness of first order-transition.This can be done by engineering the effective potential for example by using SNAIL (Superconducting Nonlinear Asymmetric Inductive eLement) structures to achieve a smaller Kerr coefficient.Our results open the way to using criticality for sensitive detection of electromagnetic fields in superconducting-circuit experimental platforms.

DATA AVAILABILITY
The calibration is done in two steps: measurement of the system gain with a spectrum analyzer and measurement of attenuation in the input line with a vector network analyzer, at the operational frequency.
System gain calibration procedure consists of measuring the power spectrum density of the Johnson-Nyquist noise for different temperatures of a 50 Ω resistor.We use the formula for thermal noise power with equal source-load impedances P T = ∆f k B T , from which the received power per unit of bandwidth can be written in terms of noise temperatures of source T and of the amplifiers T [RT] , T [HEMT] , multiplied by the system gain The term G [HEMT] = 5.2 ± 0.1K is the combined noise temperature of amplifiers, where the largest contribution comes from the HEMT amplifier.Next, we extract the total gain G [Σ] = 82.4± 0.2dB by linear interpolation.
To obtain the attenuation A = S 21 −G [Σ] , we replace in the input line the noise generated by the resistor with the actual probe field and we use a vector network analyzed to find the transmission coefficient S 21 .Knowing the attenuation we can then find the power P p that reaches the sample, and therefore the average number of photons We estimate the calibration errors to a maximum of 0.35 dB at every value of n, including calibration of the power after the mixer using the SPA, calibration of the system gain and attenuation, and errors due to reconnection, see also Supplementary Information Section II.This value is plotted as the error in n in Fig. 5.

APPENDIX B: THEORETICAL MODEL FOR PHASE TRANSITIONS
To understand the phase diagram presented in Fig. 2, we examine the dynamics generated by the general Hamiltonian for the parametric oscillator with a Kerr nonlinearity K and with the pump strength α = |α| exp(iθ P ), (7) By using the following definition for the quadratures Here the last terms ξ Qin and ξ Pin are total quadrature input noises, while α c (0) = (κ + γ)/2 is the pump amplitude corresponding to the critical point at ∆ = 0 observed in the phase diagram.Indeed, considering the linearization of Eqs.(11,12) and with the notation is the determinant whose zeroes yield the separation between phases in Fig. 2. For values of |α| and ∆ not too far from the transition, we see from the equations above that the time scales of the dynamics separate: Q becomes a slow variable and P a fast variable.As a result, P will follow Q adiabatically, and by imposing Ṗ = 0 we can obtain this adiabatic value of P from Eq. ( 12) and use it in Eq. ( 11).After neglecting some relatively small terms, we obtain the equation of motion of Q in terms of an effective potential U(Q) with This effective potential explains fully the phases observed experimentally, see Fig.
In addition, for ∆ > 0 and if |α| > α c (0) and |α| < α c (∆) we get two maxima at ±Q max , where Thus, depending on the detuning and pump strength, the system can be localized in either the central well (socalled zero-amplitude state) or in one of the wells at Q min (so-called oscillatory states), see Fig. 2.

APPENDIX C: QUANTUM EFFICIENCY
Here we give a brief overview of the general quantum theory of non-number-resolving detectors.These detectors -also called bucket detectors -are common in quantum optics, and they have two states, denoted here by 0 and 1 + .In this framework, the detection can be regarded as a form of hypothesis testing [77,78].The associated POVM operators are [79] resulting in probabilities p 0 = Tr{Π 0 ρ} = ∞ n=0 (1 − η) n P n and p 1 + = Tr{Π 1 + ρ} = 1 − p 0 , where ρ is the state at the input of the detector with photon number probabilities P n = Tr{|n⟩⟨n|ρ} [35,79,80].The efficiency η is identical to the probability p 1 + of detection when ρ = |1⟩⟨1| is a Fock state consisting of a single photon.
m! |n⟩⟨m| with n average number of photons is present at the input, we have p 0 = exp(−ηn) and p 1 + = 1 − p 0 .This reproduces exactly the Poissonian statistics with number-of-photon probabilities P n = nn n! e −n of the coherent field.Also the single-photon detection efficiency can be related to the coherent-state detection probability at n = 1; we have The next step is to introduce the dark count events, which occur with probability p dark .We can obtain the probability that the detector does not click p 0∧¬dark as the product of the probability p 0 and the probability p ¬dark = 1 − p dark of no dark counts, p 0∧¬dark = (1 − p dark )p 0 , assuming that the two processes are independent of each other [81].
Then, the probability of a detection event (resulting from either an input photon or a dark click) is On the other hand, for a coherent state with n average number of photons, we have p 0 = exp(−ηn) and p 1 + = 1 − p 0 , or and as a result we find Another way to obtain this relation is by combining p 0∧¬dark = p 0 p ¬dark , and p 0 = exp(−η n n).
Analogies with detector models from quantum optics We can get some further insights by considering a minimalist model for switching, often used in quantum optics.We treat the detector as a box containing a single mode a with internal dissipation rate γ and having a port with coupling κ -which allows us to input a coherent state with amplitude b and provides an additional dissipation channel.Note that the average number of photons in the input field in a time interval t is n = |b| 2 t.We also assume that a switching mechanism is present in the box, acting at a rate Γ, and that the switched state provides a measurable detection signal.In our experiments, this signature is provided by high-amplitude oscillations.If the information about the switched state is not recorded, the switching is seen by the mode a as another dissipative channel.Overall, we have where D[a] is the Lindblad operator defined as and the last term results from the Hamiltonian evolution under the continuous input field b.The solution of this equation is a stationary coherent state Next, we assume that we are able to register the switching events, for example if the switched configuration corresponds to a state that has an easy observable experimental signature.We can now choose to record or not these events by setting a threshold value; let us call η th the fraction of events recorded in this way.From Eq. ( 24) and we unravel only the last part.This results in , where ρ a is the (unnormalized) density matrix conditioned on nonswitching [82].The probability associated with recording a count in a time dt is dp 1 + = η th ΓTr(a † aρ a )dt, resulting from the action of the measurement operator √ η th Γa.
Putting everything together, we arrive at The solution of this equation (see also [50]) can be obtained as Since on the other hand n = |b| 2 t and dp 1 + = η|b| 2 exp(−η|b| 2 t)dt, we can identify the efficiency as From this relation one can clearly see that η th can be used to modify η from zero to a maximum value 4κΓ/(Γ+ κ + γ) 2 ≤ 1.The latter quantity can reach a maximum value of κ/(κ + γ) if the rate Γ can be adjusted such that Γ = κ + γ.Finally, if the internal losses are negligible we can tune Γ such that Γ = κ, achieving an exact balance between the rate of input photons and their recording by switching.In this situation η = η th and therefore in this optimal case the setting of the threshold fully determines the efficiency.
Noise equivalent power A convenient measure of sensitivity typically used in single-photon detection and bolometry is the noise-equivalent power [70][71][72].Consider a switching process as above, where the events happen randomly such that the probability that there is a switching event in a short time interval dt is Γdt.We define a generic random signal associated with switching events at times t j as r(t) = j δ(t − t j ) and its fluctuations around the ensemble average r n (t) = r(t) − ⟨r(t)⟩.The double-sided spectral density of r n at a frequency f is in this case S rnrn (f ) = Γ, and the single-sided spectral density is twice this value.From the Wiener-Khinchin theorem we obtain the temporal mean square r 2 n = ∞ 0 df S rnrn (f ) and since the detection is in a finite bandwidth B we obtain r 2 n = BΓ.The noise-equivalent power (NEP) is defined by considering an equivalent input signal that would produce the noise r n .The power of this signal should be P n = 1 η r n hf -higher than just r n hf by a factor of 1/η to compensate for the limited efficiency of the detector.Then, NEP is defined as NEP = P 2 n /B and therefore APPENDIX D: SWITCHING RATES The noise terms in Eqs.(11,12) originate from the thermal fluctuations of the electromagnetic field in the input line, with correlations As a result, the system prepared in the zero-amplitude state at Q = 0 may switch either to the left or to the right wells located around ±Q min .This manifests in the measurements as dark counts.
The presence of a coherent field b = |b| exp(−iφ) at the input can be modeled by adding a probe Hamiltonian The field contains on average n = |b| 2 τ p in the time interval τ p .We use again the same rotating frame U defined by half the pump frequency and we apply the rotating wave approximation, obtaining a new set of Heisenberg-Langevin equations The equations are invariant under φ → φ + 2π.The maximum switching rate is obtained by maximizing the effect of |b| in the slow (and larger) variable Q, therefore the phase that achieves this is φ = θ P /2 + π/4.By following the same procedure as before, we can obtain a modified effective potential for the slow variable as From here we see that the effect of the probe field is to tilt the potential U(Q), which is reminiscent of the way in which a current bias applied to a Josephson junction produces a washboard potential.The minima around Q = 0 is only mildly affected by this (see also Supplementary Information Section II).Most of the effect of the probe field on the potential amounts to a lowering of the left barrier to approximately U(Q max ) − √ 2κ|b|Q max and raising the right barrier to approximately Next, we construct the associated Fokker-Planck equation for the probability density where the diffusion is D = κ+γ 2 (n T + 1/2).We simulate numerically this equation with initial conditions corresponding to the starting point SP (α = 0).The probability that no switching event has occurred until time t is obtained as (33) where Q th is the threshold value.The dark count probability is obtained by simply setting b = 0 in the formula above, The probabilities obtained from these equations can be directly compared to the experimental values.The efficiency can be obtained immediately by using p 0∧¬dark = p 0 p ¬dark , where p 0 = exp(−η n n), as explained also in the main text.To model our system, which is essentially a quarter-wavelength coplanar waveguide cavity terminated with a SQUID, we adopt the approach used in Refs.[1][2][3][4].The full classical Lagrangian is where Here C l and L l are the capacitance respectively inductance per unit length of the coplanar waveguide of total length d -therefore the total capacitance is C cav = dC l and the total inductance is L cav = dL l .Since the two JJs in the SQUID are approximately identical, the SQUID Lagrangian L SQUID simplifies to where we denote C S = 2C J , ϕ S = (ϕ J,1 + ϕ J,2 ) /2, and E S (Φ) = 2E J cos(πΦ/Φ 0 ), and the term containing Φ can be dropped if the flux variation is adiabatic [2].

arXiv:2308.07084v2 [quant-ph] 9 Dec 2023
For our purpose, it is convenient to decompose the time-dependent parameter Φ into a static bias part plus a small flux pumping part, i.e., Φ = Φ DC + Φ P (t), which means Under DC biasing we can make a linear approximation, i.e., sin [ϕ (d, t)] ≈ ϕ (d, t), for the boundary condition Eq. ( 4) and plug the ansatz into it.We obtain a constraint on the eigenmodes kd: where E L,cav = ℏ 2e 2 /L cav and E S (Φ DC ) = 2E J cos (πΦ DC /Φ 0 ).For the fundamental mode k 0 , the capacitive term on the RHS −C S (k 0 d) 2 /C cav , can be neglected [4].Taking a series expansion for LHS, we can get the approximate solution in which L S (Φ DC ) ≡ Φ0 4πIc cos(πΦDC/Φ0) .Let us return to the Lagrangian and only consider a single mode, ϕ(x, t) = ϕ k (t) cos(kx).We therefore have ϕ(d, t) = ϕ S (t) = ϕ k cos(kd).By employing the linearized DC constraint condition Eq. ( 7), the Lagrangian L sys takes the form of a nonlinear parametric oscillator where Also, for the fundamental mode k 0 , we have M 0 ≈ 1 + χ(Φ DC ).After introducing the conjugate momentum n k = (1/ℏ)∂L sys /∂ φk = [ℏC k /(2e) 2 ] φk , by expanding cos ϕ S to the sixth order and taking advantage of small pump and field amplitudes, i.e., πΦ P /Φ 0 ≪ 1 and ϕ S ≪ 1, we arrive at the Hamiltonian for the parametric oscillator The quantized Hamiltonian can be rewritten in terms of bosonic creation and annihilation operators (with commutation relations [ϕ k , n k ] = i and a, a † = 1), defined by where Anticipating the pump frequency ω P ≈ 2ω k , we may write Φ P (t) = Φ P e iωPt + Φ * P e −iωPt and obtain in which The value of kd should be determined from the constraint condition Eq. ( 7).Hereafter we will only consider the fundamental mode k 0 and hence will drop the subscripts of α 0 , K 0 and Λ 0 for convenience.Still, for χ ≪ 1, α can be approximated as The values of the non-linearity coefficients, K and Λ, are calculated as follows.
First, the Josephson coupling energy E J and the SQUID capacitance C S , determined from the critical current I c ≈ 8 µA of a single JJ and its plasma frequency ω pl /(2π) ≈ 80 GHz, are E J ≈ 2.633 × 10 −21 J and C S ≈ 192.4 fF, respectively.At the operational point for detection, the DC flux applied is Φ DC /Φ 0 ≈ 0.3618.
Next, we extracted L cav ≈ 2.023 nH (or equivalently, the bare resonance frequency ω λ/4 /(2π) ≈ 6.179 GHz) and L S (Φ DC ) ≈ 45.94 pH from the fitting of the resonant frequency under DC-flux bias sweep, using the formula derived from Eq. ( 8) (and see also Ref. [5]) To better fit the data, we have adopted the method used in [6], which takes the finite loop inductance (L loop = 4.846 pH, obtained through Ansys HFSS simulation) of the SQUID into account, such that L S (Φ DC ) is corrected to where ϕ min − (Φ DC ) is determined by minimizing the dc SQUID potential [6, 7] with negligible ∝ I term, ϕ − ≡ (ϕ J,1 − ϕ J,2 ) /2 and the screening parameter β L ≡ 2L loop I c /Φ 0 ≈ 0.0375.Figure S1 shows the fitting result.From the fitting curve we extract the resonance frequency at zero flux bias 6.117 GHz, and the inductive participation ratio at the operational point χ 0 ≈ 0.0227.Alternatively, the cavity capacitance and inductance can be determined based on the design, as L cav = L l d ≈ 1.917 nH and C cav = C l d ≈ 744.6 fF.

Measurement setup
The full experimental setup, including the instruments at room temperature, is presented in Fig. S2.The K-JPA chip is thermally anchored to the mixing chamber (≈ 20 mK) of a dry dilution cryostat.The chip has three ports: one is used for the signal and the two ports are used together for DC biasing and flux pumping.A 10-MHz Rubidium frequency standard serves as the common reference clock.For the probe signal, the intermediate frequency (IF=14 MHz) is combined with the local oscillator (LO = 6.028GHz) using an IQ mixer, while for the pump the frequencies are doubled.To obtain a clean signal we use the technique of image band suppression with a 90-degree hybrid coupler.The same is achieved on the pump line by strong bandpass filtering using a combination of homemade bandpass cavities for 2LO and image band rejection after upconversion.The probe and pump are attenuated at different temperature stages (-66 and -46 dB respectively), and several circulators are used to provide isolation between input and reflected signals and to reduce back action.The reflected signal is amplified by a cryogenic high-electron-mobility transistor (HEMT) amplifier and a room-temperature amplifier (RTA), successively.The detection is done by a heterodyne arrangement, where the signal is downconverted to IF = 14 MHz and digitized by the data acquisition card; after this it is multiplied by exp(−i2π * IF * t) in order to downcovert it to zero frequency.
For the calibration of the mean photon number we use a resistor realized as a 50 Ω termination in one of the circulators in Fig. S2).The total system gain is extracted from the linear interpolation with temperature, as shown in Fig. S3c).It also includes frequency transition losses/amplification in the spectrum analyzer circuit part.We obtain G [Σ] = 82.4± 0.2 dB, where the error is associated to a single calibration.With the vector network analyzer (VNA) we measure the transmission coefficient S 21 in continuous-wave operation at a frequency 6.042 GHz, where the cavity is far-detuned.The total attenuation provided by all elements in the probe line (from the output represented as "-" of the hybrid coupler on Fig. S2 to cavity) can be calculated as A = S 21 − G [Σ] .Such calibration scheme assumes lossless connections to the sample and it neglects the losses in the circulators; however, these can be estimated and included as a small systematic error correction.We obtain A = −97.2± 0.2 dB.This figure includes also the room temperature −20 dB attenuator shown in Fig. S2.Room-temperature and cryogenic cables as well as connectors amount for an additional 11 dB of attenuation.Finally, the probe power P p reaching the sample is obtained by multiplying the attenuation A with the input power at room temperature.The latest is calibrated with the SPA, using the full setup including mixers, LO cancellation and filters.

Operation as detector
To operate the device as a detector, a sequence of pulses and operations is applied, see Fig. S4.First, the system is brought from the starting point (SP) to the operational point (OP) by applying the pump of amplitude α for a time τ P .Then the probe pulse is applied, with a fixed delay of 0.228 µs with respect to the pump.This delay is chosen such that the probe and the pulse arrive at the sample with the phase difference that maximizes the switching and also to ensure that the ringing is quenched.The average number of photons in the probe pulse can be obtained from the power P p obtained in the calibration described above, by using the relation between power and input field The probe field adds the following Hamiltonian to the system Hamiltonian In our experiments we have ω p /(2π) = 6.042GHz and P p = −174 dBW = −144 dBm ≈ 4.0 × 10 −18 W, corresponding to |b| ≈ 1 × 10 3 √ Hz.Finally, the readout operation is performed during a time-window τ .
Finally, all the parameters presented above are subjected to experimental errors.We have done a thorough estimation of these errors and their origin, see Table II.
FIG. S3.Calibration protocol.(a) Simplified schematic of the circuits used for measuring the gain in the measurement chain.
A 50 Ω load is gradually heated and its thermal noise power is amplified and measured at room temperature with a spectrum analyzer (SPA).(b) Simplified schematic of the experimental setup for obtaining the attenuation coefficient A from the "-" port hybrid coupler (see Fig. S2).This is done by measuring the S21 scattering parameter with a vector network analyzer (VNA).

Discussion: error contribution to n
The average photon number in a pulse n is calculated (see also Methods) from the probe power P p = AP reaching the sample as To obtain the point-to-point error in n we consider three sources of error: errors in P, errors in A, and errors caused by reconnection cables from calibration to measurement.We estimate each of these errors to a maximum of σ(dB) ≃ 0.2 dB, where we use the convention that the dB error in a power-like quantity A is σ A (dB) = 10 log(1 ± σ A /A).For small errors the ± sign does not make a significant difference.Then the accumulated error can be calculated from the standard formula of error propagation

II. THEORY OF PARAMETRIC CRITICALITY
Here we give an overview of the theoretical model of a parametrically pumped harmonic oscillator with dissipation, and we show that the phase diagram observed experimentally is reproduced by this model.Starting with the Hamiltonian we introduce the quadratures and rewrite the Hamiltonian up to a constant ℏω 0 /2 in the form where α = |α| exp(iθ P ).Note also that (P 2 + Q 2 )/2 = a † a + 1/2, in obvious analogy with the quantum harmonic oscillator.
Note that Q and P are the quadratures of the field inside the resonator, while the quadratures X and Y utilized in the main text denote the recorded quadratures of the heterodyned signal (after amplification and digital downconversion to zero frequency), discussed also later in this subsection.
It is also instructive to write the Hamiltonian in terms of the usual coordinate q and momentum p of the harmonic oscillator, obtaining Next, we introduce U = exp[−i(ω P t/2 + θ/2)a † a] and perform the rotations These are the quadratures in the rotating frame.A more familiar expression can be written in terms of the annihilation operator a, using e iξa † a ae −iξa † a = e −iξ a, namely U † aU = e − i 2 (ωPt+θ) a and its adjoint U † a † U = e i 2 (ωPt+θ) a † .Acting on H sys with U † and U we obtain the new effective Hamiltonian U † H sys U − iℏU † U and further, applying the rotating wave approximation (RWA), we find where ∆ ≡ ω 0 − ω P /2.
In terms of the annihilation and creation operator, this Hamiltonian reads which can be obtained also directly from Eq. ( 22) by applying U and the RWA.Next, by using [Q, P] = i we obtain the Heisenberg-Langevin equations of motion The terms ξ Qin and ξ Pin are noises in the Q and P quadratures and will be discussed in more detail later.The determinant of the linear truncation of this system and without the noise terms is and its zero value defines the stability line in the phase diagram in the absence of nonlinearity.

Mean-field theory and pump-induced shifts
By reconstructing the creation and annihilation operators from Eqs. (31,32) we find To have a complete analysis of the instability threshold, we discuss briefly another effect that manifests as a shift of resonance frequency [5], see Fig. S6.To estimate this pump-induced shift explicitly, we proceed as follows.Recalling the equation of motion for the SQUID at the boundary Eq. ( 4), we use the decomposition Eq. ( 6) and expand into [πΦ P (t)/Φ 0 ] 2 .The linearized boundary constraint Eq. ( 7) will get modified, that is, after averaging over time, E S (Φ DC ) should be replaced by E S (Φ DC ) (1 − |πΦ P /Φ 0 | 2 ).In particular, for the lowest mode, ω 0 (Φ P )/ω 0 (0) = k 0 (Φ P )/k 0 (0) ≈ 1 − χ |πΦ p /Φ 0 | 2 , or equivalently, using Eq. ( 15) and identifying ω 0 (Φ P ) with ω 0 (α) we have Accordingly, the modified critical threshold is determined by and the result is where µ = (κ+γ)/2 χω0 tan 2 (πΦDC/Φ0) , and the other solution with opposite sign in the front of 1 + 4µ ∆ (κ+γ)/2 − µ can be neglected with our parameter µ ≈ 5.24 × 10 −3 .Note that Eq. ( 39) reduces to Eq. ( 37) in the limit µ ≪ 1, as can be seen by performing a Taylor expansion We plot Eq. ( 39) in Fig. S6, obtaining a slight asymmetry.We can see that this effect is rather small.In the simulations, this pump-induced shift in |α c | is assimilated to experimental errors in determining the detuning.

Slow variable and effective potential
Following the ideas from Ref. [10], we notice that if θ is chosen such that θ − θ P ≈ π/2, then P becomes a fast variable following the slower variable Q. Therefore we can set Ṗ = 0, neglect certain small nonlinear terms, and extract P from Eq. ( 32) and insert it into Eq.( 31).This leads to where with U(Q)is an effective potential for the slow variable Q.This potential has very interesting properties, see Fig. S5.We note first that in the absence of nonlinearity the potential is a parabola if |α| < [(κ + γ)/2] 2 + ∆ 2 (harmonic oscillator, stable solution) and an inverted parabola otherwise (unstable solution).Let us now analyze what happens if K ̸ = 0. To start with, the value Q 0 = 0 is always a solution of the equation ∂ Q U(Q) = 0, and can be either a local maximum or a global minimum, as we will see below.
Then, if |α| < (κ + γ)/2, it follows that Q 0 = 0 is the only solution and it corresponds to a minimum.Next, we recall that for our system K < 0 so we analyze the solutions accordingly.If |α| > (κ + γ)/2 and ∆ > 0, or if |α| 2 − ((κ + γ)/2) 2 > |∆| 2 for ∆ > 0 then we obtain two minima at From Eq. ( 41) we obtain after some algebra Lastly, if |α| > (κ + γ)/2, ∆ > 0, and |α| 2 − ((κ + γ)/2) 2 < |∆| 2 , in addition to the Q min solutions above we obtain two more solutions corresponding to local maxima Note also that the expressions for Q max and Q min are consistent with the mean-field result Eq. (36).Indeed, a † a = (Q 2 + P 2 )/2 − 1/2 and since Q is the larger variable, we have approximately N + ≈ Q 2 min /2 and N − ≈ Q 2 max /2.By direct insertion into the expression Eq. ( 41) we find or explicitly A further simplification can be obtained by using the approximation ∆ , which is the case for our operational point, therefore we write the second square bracket as 3 2 ) and neglect the small term, obtaining A quantity that will become important when calculating the switching rates is the curvature around the extreme points.By direct calculation from the formula of the effective potential we obtain and Note that the curvatures do not depend on the nonlinearity K.Alternatively, we can write Fig. S7 shows the effective potential U(Q) at various points (with locations) in the phase diagram and Tab.III lists some numerical values of the potentials.We have now developed all theoretical tools to understand the phase diagram.The effective potential can be calculated at various points, as indicated in Fig. S7.At the SP, the shape of the effective potential is very simple, and it can be approximated by the second-order derivative U ′′ (Q) = 0.5216(κ + γ), therefore the stationary distribution at the SP is The square root of the variance of this distribution is and the partition function ensuring normalization as above, Z SP = √ 2πσ SP .Numerically, σ SP ≈ nT + 1/2.The instant the pump is turned on, the shape of the effective potential changes to that corresponding to the OP, and the system starts to evolve with W SP (Q) as the initial state.To have an estimate of the timescales involved let us consider what would happen in a flat potential.In this situation, we would apply Fick's diffusion law to find that the distribution takes the form σ t Fick = σ 2 SP + 4πDt Fick .An estimate can be obtained by the condition Q max = σ t Fick , which yields a value of t Fick = 0.7 µs.
In Fig. S8 we plot the effective potential U as a function of Q, observing a three-well structure.We evaluate the probability of the particle being in the central well as which is the same as the probability of not obtaining dark counts.The limits of integration are taken as ±Q th , typically close to ±Q th , since a detection event with in-between these values most likely correspond to the particle being still confined in the central well.The results of these simulation are shown in Fig. S9.In Fig. S10 we plot the natural logarithm of this quantity over a time of 4 µs (blue line).We observe that after around 0.3 µs this curve becomes nearly linear, meaning that the decay is exponential.On the other hand, from the experiments described in the main text we have for example for a τ p = 2 µs pulse that p dark = 0.256, therefore from 1 − p dark = exp[Γ dark (τ P − 0.228 µs)] we find Γ dark = 0.167 MHz.An exponentially decaying function with this rate approximates extremely well the numerical results for t delay = 0.35µs (not very different from the experimental 0.228 µs delay)as shown by the yellow line in Fig. S10.

Probe fields
Next, we will study the effect of adding a probe field at the input of the external mode.This is equivalent to adding a probe Hamiltonian of the form Here |b| is the amplitude of the coherent input field and b(t) = b exp (−iω p t) = |b| exp (−iω p t − iφ) is the free-evolution incident probe field.One recognizes here the standard form of the interaction Hamiltonian from the input-output theory.
In the rotating frame U † H p U can be made time-independent by using ω p = ω P /2 and by applying the rotating wave approximation, obtaining or, in terms of the quadratures Q and P, Alternatively, we may see this as a field added to the input operators of the external modes, or From the general relation for the input power in the external modes P  in (t)⟩ we find that the power added by the coherent input field is P probe = ℏω p |b| 2 .The energy in a time τ p added by this field is therefore ℏω|b| 2 τ p , which at the same time equals ℏω p n. Therefore n = |b| 2 τ p .This allow us to extract the value of |b| directly from the experimental data.
The total Hamiltonian with the probe turned ON is then FIG. S11.The probe field tilts the potential, see Eq. ( 79).
The overall recorded probability for the absence of switching events is then, using the notations in the main text, resulting in a measured dark count probability during the probe pulse time τ p .As done in the main text, it is convenient to separate the switchings due to the dark counts from those caused by the input photons in the form and we can write as in the main text where p 0 = exp(−ηn) = exp(−η|b| 2 τ p ), and therefore we can identify η in a consistent manner.The efficiency is therefore the ratio between the "true" rate of detection (with dark counts eliminated) and the average rate of incoming photons.These definitions are in agreement with those utilized in various other contexts in optics [14,15] and microwave photonics [16].
To illustrate the role played by the detection threshold in determining the efficiency, we measure P 1 + and p dark at n = 1 as a function of R th ranging from 0 to 10. From these data we extract p 1 + and the efficiency η as η = ln 1−pdark 1−P1 + , see Fig. (S15).The optimal threshold value is represented by a vertical dashed line.Note that in principle we can increase the efficiency by reducing the threshold, but this comes at the cost of increasing the dark count rate.Furthermore, at too low values of R th , the efficiency becomes somewhat ill-defined (slightly above 1 and then dropping abruptly to zero).This corresponds to setting Q th too much inside the central well, that is The dynamics of the quadratures can be observed experimentally in time domain.The relevant modes here are the external ones, which are amplified through the measurement chain and eventually recorded though the heterodyning The final recorded quadratures, which we call X and Y , are directly proportional to these amplified values.
where G [HEMT] and G [RT] are the gains of the HEMT and room-temperature amplifiers and the other terms are the respective noises added in each quadrature.In Figs.S16 and S17 we present the time-dependence of the quadratures X(t) and Y (t) for probes τ p = 1 µs, and for τ p = 2 µs respectively.

Kramer's theory
Kramer's theory assumes quasi-equilibrium and a deep enough potential well around Q = 0.These conditions are not strictly satisfied for our experiment at the OP, that is why we have to resort to numerical simulations.However, for different parameters α, ∆ we can easily bring the system to a point where these assumptions are valid.
Then, can obtain the switching rate from the zero-amplitude state at Q = 0 over the barrier with height U(Q max ) to either the left or the right wells located at ±Q min as twice the Kramers rate corresponding to each of these processes, Next, our goal is to derive approximate analytical results for η by using Kramer's theory.To do so, we have to examine the new effective potential U b .
The dominant effect when b ̸ = 0 comes from the change in the height of the barriers at larger values of Q -in the vicinity of Q max .To a good approximation, the new maxima occur at around ±Q max and therefore the height of the . Also, we observe that the curvatures around these maxima are not changed significantly by the presence of the probe field.This can be seen in the following way: the second derivative of U around Q 0 is dominated by the quadratic term in Q, and it is independent on the probe field.
Putting all these together, we get for the case when the probe field |b| ̸ = 0 is applied If the argument of the cosh function is small enough we can expand resulting in an efficiency For example, taking the point 4 in from Table III at a temperature corresponding to half a quantum nT = 0.5 we get Γ dark = 0.08 MHz and η = 0.84.In Fig. S18(a) we plot the efficiency η given by Eq. ( 96).As we approach the transition, the efficiency increases.However, too close to the transition we cannot use Kramer's approximations because the exponent becomes too small.We plot this exponent in Fig. S18(b), and we take as a reference the value 3.5 (see the dark gray line).In the efficiency plot, this line corresponds to η ≈ 0.5 (it does not vary that much since the prefactor of Γ does not vary much either in this range).
A smaller effect is the change in the position of the former local minimum at Q = 0, which now changes to a value Q 0 .In order to estimate this value, we assume that we can neglect terms of the order of Q 4 and higher in U probe and we find the minimum with respect to the quadratic and linear term only.We obtain the position Q 0 < 0 of the new minimum at and which is also negative (the probe field deepens a bit the potential around origin).
In other words Moreover, around ±Q max the second derivative is simply a constant U ′′ (Q max ), since the contribution of the probe field is zero.
Then the total switching rate in the presence of the probe field reads or, with the total switching rate Γ dark for escaping either to the left or to the right defined as we can write This yields at low values of the probe field |b| the following expressions for the rate Again, we can immediately obtain This result demonstrates that for low input powers the efficiency is constant, due to the fact that the switching rate is proportional to |b| 2 .In general, at larger values of n, these approximations are not well satisfied because the switching rate is highly nonlinear in the input field b.From Eq. (102) we can see that for larger values of |b| there will be additional terms in the Taylor expansion.Let us examine the next-order terms.From the exponential function, the next-order term is positive and proportional to |U(Q 0 )| 2 , which is fourth order in |b|.Another positive fourth order in |b| is obtained from the expansion of the hyperbolic cosine function.In conclusion, at higher values of n we expect to see a slight increase in the rate Γ b due to a term proportional to n2 and accordingly an increase in η, seen as a term proportional to n.Since the origin of these terms are exponential functions, we can try a phenomenological fit with the expression η ϵ (n) = η(n = 0) exp [1/(1 + cn)], where c ≪ 1 is a relatively small constant that should be determined empirically and η(n = 0) is the efficiency at n ≈ 0. In the main text we denoted by η the value of efficiency that best approximate the data at around n ≈ 1, therefore η ϵ (n = 1) ≡ η and we have We also assume that this type of dependency extends to the non-equilibrium case, which is confirmed both numerically (simulations of the time-dependent Fokker-Planck equation) and experimentally.We get very good fitting of the data for the entire range n ∈ [0, 50] for c = 0.12.We note that this correction is not large, as c ≪ 1, but still for practical applications requiring precision and for the analysis at large n it is useful.Since η refers to the efficiency near n = 1, a more transparent form of Eq. ( 106) is where ϵ = c/(1 + c) = 0.107.To simplify the notations, in the main text we refer to η ϵ (n) as η ϵ , keeping the n-dependence implicit.

Discussion
Let us review here the final results.We note that all of them depend sensitively on the distance to the transition (in frequency) which we denote briefly as So we have For the dark count rate and for the efficiency we got We can also approximate further so in the end we have These expressions, although valid only in the Kramers quasi-equilibrium approximation, allow us to provide heuristic recipes for future improvements of the device.For example, an important aim would be to reduce the dark count rate.For detectors employing activation over a barrier, a general rule of thumb is that reducing the dark count can be easily achieved by making the well deeper; however, this obviously reduces in a similar way the efficiency.This can be seen in our case from Eq. (112) and Eq.(113) above.However, one can still achieve this feat by exploiting the fact that in Eq. ( 113) the dark count rate Γ dark appears multiplied by Q 2 max .For example, if we could lower K by a factor s (K → sK) and simultaneously change the detuning of the operational point such that ϵ ∆ → √ sϵ ∆ the depth of the well will remain approximately constant U(Q max ) → U(Q max ).Now, due to the change in the prefactor we have Γ dark → √ sΓ dark but because of the multiplication with Q 2 max the efficiency will be the same η → η.Reducing and even nullifying K is possible for example by employing SNAILs (Superconducting Nonlinear Asymmetric Inductive eLements).

III. PHASE-SENSITIVE DETECTION
In continuous-wave operation, the switching rate depends on the phase φ of the probe, with the phase of the pump fixed and taken as reference.We can change φ by introducing a delay between the 14 MHz IF signal used to generate the probe and the 28 MHz IF signal employed for the pump.A delay of 71.4 ns (one period) imprints a phase difference of 2π.This technique is used also to calibrate the delay between the probe and pump in the pulse sequence as illustrated in Fig. (S4).Indeed, the delay of 228 ns utilized there corresponds to a phase difference of approximately 3π + π/4 (with precision limited in practice by the sampling time of 4 ns of the hardware).The reason for adding an extra 3π of delay is to reduce the effect of transients such that the system is stabilized at the operational point when the pulse probe comes.The 3π value is sufficient to achieve stabilization and at the same time is short enough such that the system has not decayed significantly by dark events before the probe pulse arrives.
To put in evidence the role of the phase, we perform the experiment presented in Fig. S19.Here we define the detector response as the slow variable ℜ(X + iY ) = X in order to select switching events only in one direction (either to the left or the right well of the effective potential).In this way we can observe the full 2π periodicity implied by Eqs.(75,76).Here ∆φ is defined as ∆φ = φ − (π/2 + θ P )/2.Indeed, recall that in the main text we have specified directly our choice of θ as θ = θ P + π/2.The effective potential then reads resulting in the ∆φ -dependent tilting and a subsequent oscillatory behavior of switching.In Fig. S19 a) we show two pairs of experimental traces, corresponding to switching to the left (P ← 1 + and p ← dark ) or to the right (P → 1 + and p → dark ).For values around π (mod 2π) in the first case, respectively 0 (mod 2π) in the second case, the probabilities P ← 1 + and P → 1 + drop to nearly zero.The crossing points of p ← 1 + and p → 1 + coincide with the level of p dark , which is a signature of the indistinguishability of detection results at ∆φ (mod π) = π/2.
Note that phase-sensitive detection of coherent state with well-known or controllable phase also allows us to reduce the value of p dark to half, see Fig. S19 b).In this situation we get Γ ← dark = 0.0835 MHz and the resulting left-switching efficiency is η ← = 0.76.Similarly, we can define a left-switching power sensitivity (susceptibility) and we obtain χ ← P = 1.68 × 10 17 W −1 .Similarly, the left-switching noise equivalent power is NEP ← ϵ = 2.22 zW/ √ Hz.

IV. BINOMIAL DISTRIBUTION FOR SEQUENCES WITH POISSONIAN PROBABILITY P0
To further substantiate the agreement with the photon statistics of coherent states, we present two tests for the case of τ p = 1 µs pulse duration.
In the first test (Fig. S20) we split the N = 2 • 10 4 pulse sequence into fixed number of blocks, M , with each block containing nN/M photons on average.In this test n is varied.This distribution of sequences leads to a binomial distribution with a certain probability of click/no-click for each pulse.We plot the obtained histograms (based on switching events) for such distribution, and find the expectation value and variance for the binomial distribution with nN/M photons in the block.We also calculate analytically the histograms corresponding to the binomial distribution, the expectation values and variances.Recall that P 1 + = p 1 + ∨dark and P 0 = p 0∧¬dark = p 0 p ¬dark , with P 1 + = 1 − P 0 and p 0 inheriting the Poissonian character from the input beam, p 0 = e −ηϵ n.By applying the binomial distribution we have E[X 2 ] − E[X] 2 = P 1 + P 0 N/M.
In the second binomial test shown in Fig. S21 the number of blocks M is varied, while n ≂ 1 is kept fixed, thus keeping constant the probabilities P 1 + and P 0 .We again obtain very good agreement with the theory.
To summarize, the results in Figs.S20, S21 confirm that we can get the correct binomial statistics when splitting the sequence into blocks, where the intrinsic probabilities P 1 + and P 0 in each block obey the Poisson distribution.

2 .
FIG. 2. Experimentally-obtained phase diagram.The different states of the device can be diagnosed by the statistical properties of the output quadratures X and Y , see also Table I. a) The highly-contrasted variance of X (a,top) shows the double-displaced state, while the highly-contrasted mean (a, bottom) shows the oscillatory state.The region below the border is vacuum or squeezed vacuum.At zero frequency detuning one can notice the triple-phase point where the squeezed vacuum, the double-displaced, and the coherent oscillatory phase coincide.The red crosses mark the operational point for detection.b) Photon output power dependency over pump amplitude and detuning sweep.The diagrams represent the effective potential U(Q) as a function of the slow variable Q in each relevant region.

FIG. 4 .
FIG. 4. Quadratures and thresholds (a) We present a sample of measured quadratures (X, Y ) (arbitrary units) in time-domain.The green bands denote the readout windows of duration τ .(b) Normalized histograms of detection events yielding the results Rτ (discretized values of segments ∆Rτ = 0.2), taken at every value of the average number of photons n.The duration of the pulses was τp = 1 µs and for each n we used a total number N = 2 • 10 4 of pulses.

FIG. 5 .
FIG. 5. SPD's basic characterization results.We present: a) the receiver operating characteristics (ROC), with the dashed lines indicating the maximum separation from the diagonal, corresponding to the optimal threshold; b) TPR(1 − FPR) = P1 + (1 − pdark) as a function of threshold; c) single-photon efficiencies η and dark count probabilities pdark for probe pulses of different duration but with the same average number of photons per pulse n = 1; d) measured probabilities P1 + and pdark as a function of n for probe pulses with the same duration τp = 1µs.The probability p1 + is obtained from these data as p1 + = (P1 + − pdark)/(1 − pdark) and it is fitted with 1 − exp(−ηn) and 1 − exp(−ηϵ n) with η = 0.73 and ϵ = 0.1.The calibration error in the number of photons n from one value to another (depicted as yellow or blue semitransparent areas) is estimated as ±0.35 dB.

FIG. 6 .
FIG. 6. Surface and contour plots confirming the Poisson statistics of a probe signal containing n photons.We present the overall probability to detect zero photons P0 as a function of n in pulse and threshold value.The data shown as surface plot is the experimental P0, whereas the contour plot is the reconstructed P0 = [exp(−ηϵ n)](1 − pdark) using the values pdark and ηϵ obtained experimentally for each value of the threshold R th ..

QTF
Centre of Excellence and InstituteQ, Department of Applied Physics, Aalto University, FI-00076 Aalto, Finland (Dated: December 12, 2023) I. EXPERIMENTAL SETUP System parameters and Hamiltonian FIG. S2.Illustration of experimental setup.The sample is cooled down to 20 mK in a BlueFors LD400 dry dilution refrigerator.See text for detailed information.
FIG.S3.Calibration protocol.(a) Simplified schematic of the circuits used for measuring the gain in the measurement chain.A 50 Ω load is gradually heated and its thermal noise power is amplified and measured at room temperature with a spectrum analyzer (SPA).(b) Simplified schematic of the experimental setup for obtaining the attenuation coefficient A from the "-" port hybrid coupler (see Fig.S2).This is done by measuring the S21 scattering parameter with a vector network analyzer (VNA).(c) Linear fitting of the measured power spectral density divided by kB of the thermal noise at different temperatures of the 50 Ohm termination.

52 FIG
FIG. S6.Critical threshold αc(∆), normalized with respect to the coupling rate κ and γ as a function of the detuning ∆.The blue continuous line represents the case without pump-induced shift, while the red dashed line shows the border of the transition with pump-induced shift (µ = 5.24 × 10 −3 ).

6 FIG
FIG. S8.The effective potential U(Q) at the OP has a central well with a local minimum U(0) = 0 MHz (upper picture presents a zoom-in around the zero) as well as a left and a right well with minima at ±|Qmin|.

FIG. S9 .
FIG. S9.Simulations of the probability density W (Q, t) around the central well (upper plot) and around the left and right wells (lower two plots).The initial probability density corresponds to α = 0 and it is a relatively narrow Gaussian.As soon as α is ramped to the value corresponding to the operational point, the distribution starts to spread, first expanding into the full width of the central well.After the delay time of 0.228 µs, the decay into the left and right wells is nearly exponential.At large times, the variable Q will eventually reside with 50% probability in the left or right well, and with negligible probability in the central well.

FIG. S10 .
FIG. S10.The blue line represents the natural logarithm of the probability of being in the central well, calculated as 1 − pdark(t) = −Q th −Q th W (Q, t)dQ.The dashed orange line is the linear interpolation of ln[1 − pdark(t)] by the linear function −Γdark(t − t delay ) + ln[1 − p th (t delay )] with t delay = 0.35 µs and Γdark = 0.167 MHz.The inset shows pdark(t) at short times, demonstrating that for around 0.3 µs the dark counts are close to zero, after which the system enters a regime where the decay from the central well is approximately exponential.
FIG.S14.Venn diagram for probabilities as used in the main paper.
FIG.S16.Measured quadrature data (in arbitrary units) for the case with no probe signal and for the case with the probe on (pulse duration 1 µs and n = 1 photons in each pulse).The readout windows for each measurement operation in the pulse sequence are given by shaded regions bordered by green dashed lines.The duration of a single sequence is 5 µs.
FIG.S17.Measured quadrature data (in arbitrary units) for the case with no probe signal and for the case with the probe on (pulse duration 2 µs and n = 1 photons in each pulse).The readout windows for each measurement operation in the pulse sequence are highlighted by green dashed lines.The duration of a single sequence is 6 µs.(b)(a)

[ 1 ]
FIG. S19.Characterization of phase sensitivity.(a) We extract the directional switching probabilities P ← 1 + , P → 1 + , p ← dark , and p → dark as a function of ∆φ by varying the phase of the probe field with the phase of the pump as reference.The sum of left-and rightprobabilities yield the non-directional probabilities P1 + = P ← 1 + + P → 1 + and pdark = p ← dark + p → dark .(b) Experimental demonstration of phase-sensitive detection by fixing ∆φ = 0 showing the dependence of left-switching probabilities P ← 1 + and p ← dark as a function of n.The ±0.35 dB error is marked with yellow or blue semitransparent areas at each experimentally recorded value of n.

TABLE I .
State classification by the mean and variance of quadratures with the vacuum state taken as reference.