One hundred second bit-flip time in a two-photon dissipative oscillator

Current implementations of quantum bits (qubits) continue to undergo too many errors to be scaled into useful quantum machines. An emerging strategy is to encode quantum information in the two meta-stable pointer states of an oscillator exchanging pairs of photons with its environment, a mechanism shown to provide stability without inducing decoherence. Adding photons in these states increases their separation, and macroscopic bit-flip times are expected even for a handful of photons, a range suitable to implement a qubit. However, previous experimental realizations have saturated in the millisecond range. In this work, we aim for the maximum bit-flip time we could achieve in a two-photon dissipative oscillator. To this end, we design a Josephson circuit in a regime that circumvents all suspected dynamical instabilities, and employ a minimally invasive fluorescence detection tool, at the cost of a two-photon exchange rate dominated by single-photon loss. We attain bit-flip times of the order of 100 seconds for states pinned by two-photon dissipation and containing about 40 photons. This experiment lays a solid foundation from which the two-photon exchange rate can be gradually increased, thus gaining access to the preparation and measurement of quantum superposition states, and pursuing the route towards a logical qubit with built-in bit-flip protection.


I. INTRODUCTION
The quest for a physical system suitable for quantum information processing is intensifying at the dawn of quantum computing. Microscopic entities such as single ions [1] and electrons [2] came forward as initial promising candidates. More recently, human-made superconducting circuits encode quantum bits in the form of single excitations of electromagnetic modes [3]. Despite impressive progress [4], these systems continue to be individually plagued by too many errors, blocking their deployment into scalable quantum machines.
A qubit with computational states |0 and |1 undergoes errors spanned by two processes: bit-flips that randomly swap |0 and |1 , and phase-flips that scramble the phase of quantum superpositions of |0 and |1 . Unlike phase-flips, bit-flips have a clear classical analogue. Interestingly, classical bits in a typical static random access memory have bit-flip times in the range of 10 15 seconds [5], 17 orders of magnitude larger than their quantum counterpart [6]. This observation sparked interest in a qubit whose computational states would be stable over macroscopic time-scales. A nonlinear dynamical system that is open to its environment through a carefully tailored interaction may exhibit a rich dynamical phase space hosting multiple stable steady states [7,8]. The switching between these states in dissipative Kerr oscillators has long been used for amplification [9,10] and ultra-low power classical logic [11,12]. Their stability is ensured by regular energy damping, * zaki.leghtas@ens.fr that, at the same time, decoheres their quantum superpositions. Surprisingly, there exists a mechanism, known as two-photon dissipation, that provides stability without inducing decoherence. A recent qubit, coined the cat-qubit [13], is embedded in the cavity field of a superconducting resonator that exchanges pairs of photons with its environment [14,15]. This defies the common intuition that a qubit must be well isolated from its environment. Instead, the subtle interplay of drive and dissipation pins down the cavity field on two coherent steady states with complex amplitudes denoted ±α without affecting quantum superpositions of the two.
Increasing the number of photonsn = |α| 2 in these two steady states has two opposing effects [16]. On the one hand, their distinguishability by an inevitably coupled uncontrolled environment increases linearly with n. This results in a linear increase of the phase-flip error rate. Therefore, for this system to be suitable for quantum information processing, it must operate at low photon number. On the other hand, as soon as their separation exceeds their vacuum fluctuations, that is |α − (−α)| 2 = 4n > 1, their wave-function overlap rapidly decreases, reducing random tunneling between them and hence exponentially increasing the bit-flip time. It is remarkable that, at least in principle, it is possible to reach macroscopic bit-flip times with computational states pinned by two-photon dissipation containing only a handful of photonic excitations.
Previous experiments have succeeded in implementing a two-photon exchange mechanism to observe the squeezing of a Schrödinger cat state out of vacuum [15], the dynamics of a quantum gate [17], the exponential suppression of bit-flips and linear increase of phase-flips [16]. A recent insight has been to define a cat-qubit out of the A cavity is endowed with a special mechanism (dashed left mirror) that exchanges pairs of photons (blue double waves) at variable intensity (control knob) with a cold bath. Two meta-stable pointer states emerge, represented by the blue distributions centered around amplitudes ±α. A fraction of the cavity field (blue waves) escapes through the weakly transmissive mirror and is collected by our heterodyne detector. By monitoring the signal over time, we track individual trajectories undergoing bit-flips (blue time trace). (Bottom) False-color optical micrograph of the experimental superconducting circuit in a coplanar waveguide geometry. The cavity is a λ/2 resonator (blue) that radiates a field aout through a weakly coupled 50 Ω port. It also couples to a two-photon exchange device composed of a buffer mode (red) shunted to ground through an asymmetrically threaded SQUID (ATS) as emphasized in the insets. DC currents enter through on-chip bias tees (green) and impose phase biases ϕL,R. A differential pump (purple arrows) and a common buffer drive (red arrows) are channeled through filtered transmission lines (orange).
interplay of Kerr non-linearity and single mode squeezing [18,19]. However, in all these implementations, the bit-flip time saturated in the millisecond range, limited by errors impinging from the cat-qubit tomography apparatus [16], and possible dynamical instabilities [20,21]. In this work, we aim for the maximum bit-flip time we could achieve in a two-photon dissipative oscillator. To reach this goal, we first design a circuit with the objective of removing all suspected sources of dynamical instabilities and ancillary systems that could propagate uncorrectable errors. We fabricate a two-photon exchange dipole element close to the regime where its energy landscape exhibits a single global minimum at any operating point, a possible requirement for stability [21,22]. Second, we entirely remove the tomography apparatus: our design does not contain a transmon and readout mode. Instead, we directly measure the field radiated by the cavity through a travelling wave parametric amplifier (TWPA) [23], thereby accessing individual oscillator state trajectories. We measure a bit-flip time exceeding 100 seconds for computational states pinned by two-photon dissipation and containing about 40 photons. Our design choices came at the cost of a two-photon exchange rate dominated by single-photon loss, hence losing our ability to prepare quantum superposition states and hence measuring the phase-flip rate. Guided by this benchmark, future experiments can then gradually enter the regime suitable to implement a qubit where twophoton loss is the dominant dissipation mechanism.

II. THE TWO-PHOTON DISSIPATIVE OSCILLATOR
An oscillator exchanging pairs of photons with its environment in addition to usual energy relaxation (see Fig. 1) is modeled by the following Hamiltonian and loss operators: where a is the annihilation operator of the mode referred to as the memory, 2 is the two-photon injection rate, κ 2 the two-photon loss rate, and κ a is the energy damping rate. When the two-photon injection rate overcomes the damping rate, two meta-stable pointer states emerge: where |α is a coherent state with complex amplitude α, verifying α 2 = 2 κ2 ( 2 − κ a /4) if 2 > κ a /4, and α 2 = 0 otherwise.
The two-photon dissipation mechanism is engineered by implementing a two-to-one photon exchange interaction with a dissipative mode referred to as the buffer [15], modeled by the Hamiltonian where b is the annihilation operator of the buffer, g 2 is the two-to-one photon coupling rate and d the buffer drive amplitude. In the limit where the buffer energy decay rate κ b is larger than |g 2 |, we recover Eq. (1) with κ 2 = 4|g 2 | 2 /κ b and 2 = 2g 2 d /κ b [15].

III. CIRCUIT DESIGN
Our two-photon dissipative oscillator is implemented in a circuit quantum electrodynamics coplanar waveguide architecture (see Fig.1). The memory mode is the fundamental mode of a λ/2 resonator. We measure coupling and internal loss rates κ c a /2π = 40 kHz and κ i a /2π = 18 kHz. In order to minimize dielectric losses [24], we target the relatively low frequency of ω a /2π = 4.0457 GHz. A thermal population of about 1% was measured on a twin sample using a transmon (see [25] Sec. S4 C). An undesired side effect of coupling the memory to a lossy mode -the buffer -is to increase the decay rate of the memory due to the Purcell effect. To prevent this, we designed a stop-band filter centered at the memory frequency, consisting of three λ/4 sections [26] on both routes linking the memory to its cold bath. The buffer mode consists of a metallic plate of charging energy E C /h = 73 MHz shunted to ground through an Asymmetrically Threaded SQUID (ATS) [20]. The ATS is formed by two Josephson junctions in a loopeach of Josephson energy E J /h = 37 GHz -split in its center by an inductance made of five junctions of total inductive energy E L /h = 62 GHz. This layout defines two loops that can be biased in DC flux ϕ L,R . We can hence independently control the common and differential flux through the ATS: ϕ Σ = 1 2 (ϕ L + ϕ R ) and ϕ ∆ = 1 2 (ϕ L − ϕ R ). Radio-frequency (RF) signals are routed to the ATS through a 180°hybrid coupler. The buffer drive propagates in phase through both arms of the two photon exchange apparatus. When reaching the ATS, these waves combine, inducing currents in the inductance and thereby driving the buffer mode. On the other hand, the pump propagates with opposite phases, inducing common flux in the ATS.
In the process of choosing the ATS parameters, we were guided by the intuition that dynamical instabilities would be avoided in a system with 2E J /E L 1 [21,22]. However, this criterion needs to be balanced with the requirement of large g 2 (see [25], Sec. S2). In this experiment, we favoured stability over coupling strength and chose: 2E J /E L = 1.2, a factor of 3.3 smaller than our previous implementation [20]. Moreover, we engineered a weak hybridization between the memory and buffer mode in order to minimize undesired nonlinear couplings such as the Kerr effect, with a rate estimated below 1 Hz.
Previous experiments constructed the Wigner distribution of the memory field through a non-linear coupling to a transmon qubit and its readout resonator [20]. However, the finite thermal occupation of the transmon was suspected to limit the bit-flip time to the millisecond range. Instead, we monitor our memory through a minimally invasive detection tool: a weakly coupled transmission line connected to a TWPA. This added leakage channel slightly decreases the total quality factor but has the advantage of not inducing any additional non-linear couplings to a lossy ancillary system.

IV. EXPERIMENT CALIBRATION
The ATS contributes a non-linear potential energy of the form U ϕΣ,ϕ∆ (ϕ) = 1 2 E L ϕ 2 −2E J cos(ϕ Σ ) cos(ϕ+ϕ ∆ ), where ϕ is the phase drop across the central inductance [20]. The buffer and memory modes hybridize through their capacitive coupling, and hence to the lowest order of their hybridization strength υ, we have ϕ = ϕ b (b+b † + υ(a + a † )), where the buffer zero-point phase fluctuations verifies ϕ b = (2E C /E L ) 1/4 [27]. We operate the ATS at ϕ Σ = −π/2 + ε p cos(ω p t), and ϕ ∆ = π/2. At this operating point the buffer resonates at ω b /2π = 6.1273 GHz with an energy decay rate κ b /2π = 16 MHz. By tuning the pump frequency at ω p = 2ω a − ω b and driving the buffer mode at ω d = ω b , we synthesize the Hamiltonian of Eq. (2), with g 2 = − 1 2 E J ε p υ 2 ϕ 3 b . We start the experiment by measuring the buffer mode frequency map as a function of the two DC currents. Conveniently, the desired operating point (ϕ Σ , ϕ ∆ ) = (−π/2, π/2) is easily recognisable since it corresponds to a saddle point of this map (see [25], Sec. S2). Consequently, the buffer and memory modes are first order insensitive to flux noise. The next step is to activate the RF pump and buffer drive. We pick the largest pump power that does not deteriorate the buffer and memory modes spectra. The modes frequencies are Stark shifted in the presence of this strong pump. Therefore, a precise calibration of the pump and drive frequencies is required to rigorously verify the frequency matching conditions: ω p = 2ω a − ω d , and ω d = ω b . To this end, we acquire the memory mode fluorescence as a function of detunings from these matching conditions. As the drive amplitude d is increased, the region over which the drive and pump combine to populate the memory expands around the frequency matching point (see Fig. 2a). For the remaining of the experiment, we place ourselves at the center of these regions (colloquially referred to as diamonds).
We calibrate the number of photons in the memory by measuring I 2 + Q 2 (see [25], Sec. S4) as a function of the drive amplitude, where I and Q are the in-phase and out-of-phase quadratures of the radiated field acquired over an integration time T m . (see Fig. 2b). Three notable features are apparent. First, in the limit of strong drives, the radiated energy scales linearly with the drive amplitude, a signature of the conversion of 1 buffer photon to 2 memory photons. This is in stark contrast with the common quadratic scaling for a driven harmonic oscillator. Moreover, the offset of this asymptote from the origin excludes the Kerr effect as the underlying process (see [25], Sec. S2 D). Second, the output power is close to zero up until a critical drive amplitude, reminiscent of a non-linear dissipative phase transition [28]. This transition occurs when the two-photon injection rate overcomes the memory losses, setting the scale for the product d g 2 to its value at the critical point d g 2 = κ a κ b /8. Furthermore, we compute the classical dependence ofng 2 2 on d g 2 , demonstrating that it is invariant under a change of g 2 (see [25], Sec. S4). Importantly, quantum fluctuations blur the transition out of vacuum resulting in a non-scale invariant curvature, a striking deviation from the sharp transition expected in the classical regime. A full quantum model is necessary to capture this third notable feature, from which we extract the key parameter g 2 /2π = 39 kHz, and deduce κ 2 /2π = 370 Hz andn for every d . This places our experiment in the regime where κ a /κ 2 = 150 1. In the future we will increase the hybridization factor υ to enter the regime suitable for a qubit implementation: where the two-photon exchange rate largely dominates the cavity losses. Finally, with the photon number calibration in hand, the measurement records I and Q are rescaled to respectively coincide with a measurement of (a + a † )/2 and (a − a † )/2i.
In fact, the observed phase transition corresponds to a spontaneous symmetry breaking, where the cavity field adopts two opposite phases (or any quantum superposition of the two in the absence of single photon loss). We observe the emergence of these two phases by continu- ously acquiring, for each drive amplitude, 10000 times the I-quadrature of the radiated field integrated over T m = 1 ms, for a total measurement duration of 10 s (see Fig 2c). For the lowest drive amplitudes ( d /2π 2.7 MHz), the cavity state remains in the vacuum, as signaled by the Gaussian distribution centered at I = 0. This distribution then broadens around the critical point (2.7 MHz d /2π 3.5 MHz), due to the significant overlap of the distributions of states |±α at small α and the multiple flips in between during the acquisition time T m = 1 ms. For (3.5 MHz d /2π 4 MHz), the two states are well resolved, and their approximately equal weights hint towards a bit-flip time larger than the acquisition time of 1 ms and smaller than the full experiment duration of 10 s. For d /2π 4 MHz, the field stays pinned to one of the two computational states, hinting towards bit-flip times exceeding 10 s.

V. TRAJECTORIES AND BIT-FLIP TIMES
We access the dynamics of the memory by tracking individual trajectories over time (see Fig. 3). For each trajectory, we set the drive amplitude at a fixed value d , and record the I-quadrature of the radiated field. In order to resolve quantum jumps, we set the integration time T m to be simultaneously smaller than the bit-flip time and sufficiently large to average out the heterodyne detection noise. To capture the statistical properties of each trajectory, we plot the cumulative distribution function of the interval between two consecutive jumps, denoted τ jump . It shows approximately an exponential law, revealing an underlying Poisson process. The average of τ jump , that defines the bit-flip time, undergoes a spectac-ular increase from 1 ms to 0.3 s to 206 s for an increase of photon numbern from 11 to 28 to 43. With respect to the bare cavity lifetime of 2.7µs, this represents a 10 8 increase of the bit-flip time, and (although inaccessible with our measurement scheme) an estimated 2 × 43 = 86 fold decrease of the phase-flip time.
We quantify the scaling of the bit-flip time with the photon number by repeating the trajectory acquisition procedure for multiple drive amplitudes. From each trajectory we extract the bit-flip time and the corresponding photon number, and display them in Fig. 4. We observe two distinct regimes. Initially, the bit-flip time rises exponentially multiplying by a factor of about 1.4 for every added photon. In theory, in the limit where κ a /κ 2 1, this factor would approach e 2 ∼ 7.4 [13]. In this experiment we favoured stability over coupling strength, placing ourselves in the opposite regime κ a /κ 2 ∼ 150, which is expected to decrease this factor as confirmed by numerical simulations (see [25], Sec. S5). Forn 40 photons, the bit-flip time saturates in the 100 second range. Although the origin of this saturation is yet to be established, its timescale is compatible with the measured rate of highly correlated errors in a large array of qubits [29], possibly due to high energy particle impacts [30,31].

VI. CONCLUSION
In conclusion, we have measured timescales of order 100 seconds for bit-flips between pointer states of a twophoton dissipative oscillator containing about 40 photons. To reach these macroscopic bit-flip times with mesoscopic photon numbers, we designed a two-photon exchange circuit in a regime expected to circumvent dy- namical instabilities, and employed a minimally invasive detection tool that collects the oscillator's radiated field.
Our experiment thus puts a scale on the bit-flip times attainable with two-photon dissipation, a necessary mechanism for quantum information processing with cat-qubits [32]. Future work could be to uncover the phenomena causing these bit-flip events [30,31,33] by monitoring oscillator trajectories over timescales of days or weeks. Also, gradually increasing the two-photon ex-change rate so that it supersedes all loss mechanisms would lead to the regime suitable for the cat-qubit where quantum superposition states can be prepared and measured, thereby paving the way towards fully protected chains of cat-qubits [34,35].

VII. ACKNOWLEDGMENTS
We thank Lincoln Labs for providing a Josephson Traveling-Wave Parametric Amplifier. The devices were fabricated within the consortium Salle Blanche Paris Centre. This work was supported by the QuantERA grant QuCOS, by ANR 19-QUAN-0006-04. Z.L. acknowledges support from ANR project ENDURANCE, and EMERGENCES grant ENDURANCE of Ville de Paris. This work has been supported by the ParisÎle

S1. DEVICE FABRICATION AND WIRING
a. Wafer preparation We sputter 120 nm of Nb on a 2-inch intrinsic silicon wafer, with a 280 µm thickness and a resistivity larger than 10 kΩcm. We fabricate twelve 10×11 mm chips on the same wafer. We dice the individual chips at the end of the fabrication process and select the sample that is best suited for the experiment.
b. Circuit patterning We pattern the large features of the circuit (> 5 µm) using laser lithography. We spin positive resist (S1805), expose the pattern, then develop in MF319 for 1 min and rinse in deionized (DI) water for 1 min. The wafer is then etched in a reactive ion etching machine with a SF6 plasma and a 10 s overetch. A 30 min lift-off step follows in a 50°C acetone bath with sonication. Finally, the sample is rinsed in isopropyl alcohol (IPA) for 1 min, blow dried and cleaned in an O 2 plasma for 20s, thus stripping residual organic contaminants.
c. Junction patterning Our Josephson junctions are fabricated from Dolan bridges patterned with electron beam (e-beam) lithography. We spin two layers of resist: first, methacrylic acid/ methyl methacrylate (MAA EL13) baked for 3 min at 185°C and second, poly(methyl methacrylate) (PMMA A3) baked for 30 min at 185°C. Once the e-beam patterning completed, we develop in a IPA:H 2 O (3:1) bath at 6°C for 2 min, rinse for 10 s in IPA and blow dry. Finally, residual organic contaminants below the bridges are stripped by an O 2 plasma for 10 s.

d. Junction deposition
The wafer is then introduced in an e-beam evaporator. We start with a 2 min argon milling step at an angle of ± 30°to prepare for a good electrical contact with the Nb layer. We deposit two layers of aluminium (35 nm then 70 nm thick) at an angle of ± 30°, separated by a static oxidation in a pure O 2 atmosphere at 10 mbar for 10 min. Before venting to air, the chamber is filled with 300 mbar of O 2 for 5 min. We lift-off in a 50°C acetone bath for 1 h, transfer the wafer to a new acetone bath for 5 min and sonicate for 10 s, then rinse in IPA and blow dry. Images of the fabricated junctions are displayed in Fig. S1.  The inductive energy of the ATS writes [16] U ϕΣ,ϕ∆ (ϕ) = 1 2 where E L is the inductive energy of the shunt inductance, E J ± ∆E J are the Josephson energies of the left and right Josephson junctions respectively, ϕ is the superconducting phase difference across the ATS and ϕ Σ,∆ = (ϕ L ± ϕ R )/2 are related to the common and differential flux threading the ATS with ϕ L,R threading the left and right loop of the ATS respectively. a. Symmetries The ATS potential has the following translational symmetries and an inversion symmetry center at (ϕ Σ , ϕ ∆ ) = (π/2, π/2) such that U π/2+ϕΣ,π/2+ϕ∆ (ϕ) = U π/2−ϕΣ,π/2−ϕ∆ (−ϕ) (S3) S3. Lumped element model of the circuit. The buffer (red) with bare frequency ω b,0 /2π is constituted with an ATS (inductive energy EL, mean Josephson energy EJ , asymmetry ∆EJ ), and a capacitor with charging energy EC . The ATS loops are threaded with fluxes ϕL, ϕR (green) and is connected to the memory (blue), with bare frequency ωa,0/2π. The phase ϕ is indicated with an arrow (black).
Combining these three symmetries gives rise to a second inversion symmetry center located at (ϕ Σ , ϕ ∆ ) = (π/2, −π/2). Hence, all the information about the system is contained in the region ϕ Σ ∈ [0, π], ϕ ∆ ∈ [−π/2, π/2]. Note that provided ∆E J = 0, we have additional symmetry axes ϕ Σ = 0 and ϕ ∆ = 0 such that b. Saddle points Let us study the potential around the inversion symmetry points (ϕ Σ , ϕ ∆ ) = (π/2 + , ±π/2 + δ) For small and δ For = δ = 0, the potential reaches its minimum at ϕ min = 0. At , δ = 0, we search for a first order perturbation of ϕ min . Solving for ∂ ∂ϕ U (ϕ min , , δ) = 0, we get Around the minimum ϕ min , the second derivative of the potential with respect to ϕ, i.e. the inductive energy of the ATS writes The ATS inductive energy has no linear terms in or δ so the points (ϕ Σ , ϕ ∆ ) = (π/2, ±π/2) are critical points of the inductive map of the ATS as a function of and δ. Its quadratic dependence around the critical point has the following matrix representation the determinant of which writes The determinant is negative (provided ∆E J < E J ) hence the critical point is a saddle point. This property is used to tune the DC working point experimentally (see Fig. S5). When ∆E J = 0, the two points (ϕ Σ , ϕ ∆ ) = (π/2, ±π/2) are non equivalent saddle points of the ATS with inductive energy E L ∓ 2∆E J .

B. Circuit Hamiltonian
The dynamics of the circuit displayed in Fig. 1 is well captured by a reduced lumped element model (see Fig. S3) with the following Hamiltonian [16]: where a, b are the memory and buffer annihilation operators. The buffer's angular frequency verifies ω b,0 = √ 8E L E C / , where E L , E C are the energies associated to the buffer's inductive and capacitive shunt respectively. The angular frequency of the memory is denoted ω a,0 . We denote 2E J the sum of the Josephson energies of the single junctions composing the SQUID loop, and 2∆E J their difference. During fabrication we aim for the smallest possible junction asymmetry, however in practice we are left with ∆E J /E J ≈ 0.5% (see [25], Sec. S3 A) which leads to spurious Kerr and cross-Kerr effects. We neglect ∆E J in the rest of the analysis. The ATS is threaded with a common and differential flux ϕ Σ,∆ = 1 2 (ϕ L ± ϕ R ), where ϕ L,R is the flux threading the left and right loop of the ATS. In the limit where the hybridization factor υ between the buffer and memory is much smaller than 1, the phase across the ATS denoted ϕ, verifies By flux pumping the ATS around a well chosen DC working point [16] ϕ Σ = π 2 + p cos(ω p t) we engineer a two-to-one photon exchange Hamiltonian between the memory and the buffer, provided the pump frequency ω p is close to the matching condition ω p = 2ω a − ω b . This two-to-one photon exchange Hamiltonian converts the strong single photon losses of the buffer into an effective two-photon loss channel for the memory. Likewise, a microwave drive at frequency ω d close to the buffer frequency, is converted into an effective twophoton drive of the memory (or squeezing) at frequency (ω d + ω p )/2. By definition, this frequency is close to the memory frequency.
For the memory, the combination of the two-photon loss and two-photon drive, stabilizes two coherent states with frequency (ω p + ω d )/2 of equal amplitude and opposite phase. The heterodyne demodulation frequency ω dm for the memory is constrained accordingly (S13) By going in the frame rotating at frequency ω dm for the memory and ω d for the buffer and performing first order rotating wave approximation (RWA), the Hamiltonian (S11) writes [16] and ω a and ω b are respectively the memory and buffer frequency accounting for the AC-stark shift due to the pump [16]. Incorporating the buffer drive and the dissipation of the two modes, the dynamics of the system is governed by where d is the buffer drive strength and κ a and κ b are the single photon loss rate of the memory and the buffer respectively.
We gain further insight on the dynamics of the system by performing the adiabatic elimination of the buffer. This is justified provided g 2 κ b . Following the method of [36], we find that the reduced dynamics of the memory is given by At ∆ a = ∆ b = 0, we recover eq. (1)

D. Steady-state photon number
In this paragraph, we derive the stationary mean photon number in the memory using a semi-classical approximation In the interaction picture, the dynamics arising from (S15) writes We perform a mean-field approximation on mode a and b, and compute the steady-state of the simplified dynamics. The operator a and b are replaced by their mean value, the complex numbers α and β.
This system always admits a solution in which the memory is in vacuum and corresponds to This solution is stable provided it is the only solution of eq. (S19) for a given set parameters. Assuming α = 0, we can write where θ a = arg(α). Solving for β in the first equation and injecting in the second one, we get a. Zero-detuning When ∆ a = ∆ b = 0, eq. (S22) simplifies into leading to We recover the critical point when the two-photon drive overcomes the cavity dissipation.
In the absence of calibrated input or output lines, the power radiated by the memory is defined up to a constant, in particular the quantity |g 2 α| 2 writes which, as a function of d g 2 , has a slope 1 and an xintercept κ a κ b /8 . In terms of the single mode effective quantities, eq. (S24) rewrites for 2 ≥ κ a /4 If instead of two-photon dissipation, Kerr effect of amplitude K was limiting the amplitude of the pointer states [18,37], a similar semi-classical analysis predicts a mean photon number for 2 ≥ κ a /4 which is qualitatively different from what is observed in this experiment. b. General case Since θ a and α are on separate sides of eq. (S22), we can geometrically solve the system in the complex plane. The right-hand side is a circle of radius | d /g * 2 | centered on z. The left-hand side is the positive real axis. In this picture, there can be 0, 1, or 2 intersections between this circle and the real positive axis, giving rise to 1, 3 or 5 solutions for the system, one for vacuum plus two for each intersection since ±α are both valid solutions. Experimentally observed solutions are the ones that give rise to the largest field in the memory. Hence |α| 2 , the mean photon number in steady-state, writes: otherwise. (S26) In the ∆ a , ∆ b coordinates, the region where α 2 is nonzero forms a diamond shape. Here after, we provide the equation for the borders of this feature colloquially referred to as a diamond.
The ∆ a , ∆ b -plane is divided in two domains depending on the sign of the quantity ∆ a ∆ b − κ a κ b /4 .
The top-right and bottom-left edges of the diamonds are located in the positive domain and given by Note that the slope of this edge is given by −κ b /κ a . The bottom-right and top-left edges are located in the negative domain and given by Hence the border of the diamond only depends on the product d g 2 and does not carry information on g 2 nor d independently. Moreover, we cannot determine g 2 and d even with the full diamond information (not only the edges). Indeed, from the measurement of the rates κ a , κ b and the knowledge of the applied detunings ∆ a , ∆ b , one has independently access to the quantity In the absence of photon number calibration, we only learn from eq. S26 that the power radiated by the memory is proportional to Thus, we only have access to the product d g 2 but not g 2 or d independently.

S3. TUNING THE EXPERIMENT
Our experiment requires DC and RF powering. For optimal operation, two DC currents and two RF powers and frequencies need to be fine tuned. In this section we describe the sequence of calibration experiments we perform to set the working point.

A. DC currents
The first calibration experiment we perform is to extract the buffer and memory frequencies as a function of the common and differential flux in the ATS loop (see Fig. S4). From these maps we identify the circuit parameters and locate the ATS saddle points.
In the following we describe the measurement protocol to acquire the buffer frequency flux map. We set a tone at frequency f on the buffer port and record its reflected amplitude and phase as a function of the DC voltages V Σ,∆ = (V L ± V R )/2 (see Fig. S2). The physical controls V Σ,∆ are transformed to the common and differential flux basis ϕ Σ,∆ to match the symmetries of the circuit Hamiltonian (S11). A variation in the reflected signal is detected at flux points ϕ Σ,∆ (f ) where the buffer frequency enters the vicinity of f . This sequence is repeated by scanning f in between 5.2 GHz and 9 GHz in steps of 100 MHz. In Fig. S4, the frequency f is encoded in the color of pixels located at ϕ Σ,∆ (f ). We repeat the same protocol on the memory port to extract the memory frequency flux map. The theory plots in Fig S4 are obtained for the numerical diagonalization of the Hamiltonian in Eq. (S11). From the ATS symmetries, we know that there exist two nonequivalent families of saddle points, those generated from (ϕ ∆ , ϕ Σ ) = (−π/2, π/2), and those from (ϕ ∆ , ϕ Σ ) = (π/2, π/2). The junction asymmetry ∆E J lifts the degeneracy of the buffer frequency between these points.
We refine the flux and frequency sweeps around these saddle points in order to precisely pin down their location. In Fig. S5, we directly display the reflected amplitude on the buffer port at different frequencies f . A saddle point is easily identified as the closing of the buffer frequency contour line. Note that the saddle point at (−π/2, π/2) appears at f b1 = 6.00 GHz as shown in the top middle panel. The second one appears at f b2 = 6.04 GHz as shown in the bottom middle panel.
The parameters entering Eq. (S11) are listed in table I. The charging energy E C is extracted from 3D finite elements electromagnetic simulations. The inductive energy E L and the junction asymmetry ∆E J are computed from the buffer frequencies at the saddle points that verify, in the weakly hybridized limit: The Josephson energy E J is extracted from the maximum buffer frequency f bmax measured in the buffer flux map of Fig. S4. This maximum value f bmax = 8.9 GHz is reached for (ϕ ∆ , ϕ Σ ) = (0, 0) and verifies, in the weakly hybridized limit: The memory frequency f a is extracted from the memory frequency flux map in Fig. S4 at the saddle points. Due to the weak hybridization with the buffer, the memory frequency difference at the two saddle points is negligible. Finally we numerically find the hybridization factor υ that produces a memory frequency flux map in agreement with the data in Fig. S4. κ i a , κ c a are the internal and coupling loss rates of the memory, κ b is the loss rate of the buffer, g2 is the two-photon coupling rate, κ2/2π is the two-photon dissipation rate andnsat is the average number of photons at the bit-flip time saturation. The following parameters are given with confidence intervals: κ i a ∈ [15,22]

B. Phase-locking
In the laboratory frame, at any point in time, the phase of the pointer states resulting from the junction mixing process is given by where θ p and θ d are respectively the pump and drive tone phase. The pump tone is directly generated by the microwave signal generator and pulsed via a microwave switch whereas the drive tone is pulsed via an IQ-mixer (Fig. S2).The resulting phases of the tones are where θ LO p and θ LO d are the Local Oscillator (LO) phases of the microwave generator and θ IF d is the Intermediate Frequency (IF) signal phase delivered by the Arbitrary Waveform Generator (AWG) channel to pulse the drive tone. The radiated signal from the memory is demodulated in a frame with phase where θ LO dm is the phase of the demodulation LO and θ IF dm is the phase of the Analog-Digital Converter (ADC). In order to phase-lock the pointer states with the demodulation frame, we should ensure The three LOs are generated with a single four channel Anapico signal generators, and the two IFs with a single Quantum Machines OPX. The accuracy of these instruments ensure that all the LOs share the same time reference and all the IFs share the same time reference. However, given the high frequencies at stake, the instrument sharing the same 50 MHz clock is not sufficient for this two time references to be considered identical. The LO time is referred to as t and the IF time as t . Hence and the phase-locking condition imposes the frequency matching conditions

C. Pump and drive frequencies
Once the DC biases are tuned at one saddle point and the memory and buffer resonance frequencies are determined by direct spectroscopy in reflection on their respective ports, we proceed to tune the pump and drive tones. The pump frequency is determined via two-tone spectroscopy: a weak drive tone is used to perform buffer spectroscopy while sweeping the pump frequency around the frequency matching condition ω p = 2ω a − ω b . As this operation is being performed, we also perform the heterodyne detection of the field radiated by the memory (Fig. S6).
When the two-to-one photon exchange is resonant, a sharp feature is observed within buffer resonance, referred to as a diamond, and the memory starts to radiate power. The discrepancy between the ideal and measured diamond shape is used as a witness for the appearance of higher order processes. The pump amplitude is set so as to maximize g 2 while mitigating detrimental higher order effects.
In order to minimize the amount of data collected and accurately zoom on the diamond feature Fig. S6a is acquired in the following way : • at a fixed pump frequency (x-axis of Fig. S6) the buffer spectroscopy is done by varying the buffer IF frequency with a fixed LO frequency. Due to the frequency constraint of eq. (S36) the heterodyne detection of the radiated memory field is done by varying the memory IF frequency with a fixed LO frequency.
• for the demodulation frequency to remain close to the memory frequency while varying the pump frequency, (y-axis of Fig. S6) the pump and drive LO frequencies are varied in opposite directions. In this way, we have ω p + ω LO b = 2ω LO dm = cst. In these coordinates, ∆ b and 2∆ a are varied along the x-axis and ∆ b is varied along the y-axis. If there was no Stark shift, the buffer resonance would be a diagonal line of slope 1 and the memory line (when the two-photon drive is tuned with the memory mode) a vertical line. In practice, these two lines are distorted and we numerically fit the buffer and memory frequencies as a function of the pump frequency to perform the change of basis leading to the diamond of Fig. S6d in the ∆ a , ∆ b coordinate system.
To perform this change of coordinates, we evaluate the following functions from measurements by linear interpolation: For each data point, we can now compute the actual value of ∆ a and ∆ b . This enables us to display radiated energy by the memory in the basis of Hamiltonian eq. (S15): previously distorted, the diamonds recover their shape.
By construction, diamond center should coincide with the zero detuning point. By exploiting the diamonds inversion symmetry, we can verify that the maximum of the auto-correlation function gives back the zero detuning point, up to a slight discrepancy due to the diamond imperfection. Experimentally, we find the zero detuning point as the convergence point of the diamond feature at vanishingly small drive amplitude (see Fig. S6c). All three centers (construction, auto-correlation, experimental) are shown in Fig. S7 and lie in a small region at the center of the diamond.

S4. PHOTON NUMBER CALIBRATION
It is of central importance that our macroscopic bit-flip times were observed for states containing only a few tens of photons. Indeed it is only in the low photon number regime that this system can operate as a coherent qubit. A reliable calibration of the number of photonsn is therefore key to this work. In this section, we describe two calibration methods, and check their consistency. We start by computing the mapping between the cavity field properties and the measured quadratures. Then, we detail the method used in the main text to calibraten from the curvature at the critical point. Finally, we describe a method relying on the measurement of the detection efficiency η.

A. Heterodyne detection
The heterodyne detection of the field radiated by the memory results in two signals that are integrated over a integration time T m to give out (I, Q) pairs time traces where G is the gain of the amplification chain, κ c a is the coupling rate of the memory, η is the quantum detection efficiency, ρ t is the instantaneous state and dW I , dW Q are the noises added to each quadrature that verify dW 2 I = dW 2 Q = dt. The statistics of the distribution of the (I, Q) pairs collected over time gives information about the memory state. In particular, we can verify that in the general case [38, 39] and in the limit of small T m wheren is the mean photon number, I 2 + Q 2 is the statistical average over the (I, Q) pairs collected over time and ρ ∞ is the steady-state density operator of the cavity. In our specific case, we have verified both numerically and experimentally that this limit is practically reached for T m = 10 µs. In eq. (S40), the offset GT m can be calibrated out from the average of I 2 + Q 2 when the cavity is in vacuum which results in the average energy radiated by the cavity over a period T m of For various values of buffer drive amplitude d , we measure the average energy radiated by the memory for a duration T m = 10 µs according to eq. (S40) (see Fig. S8) which is proportional ton. The following paragraphs aim at calibrating this proportionality constant. First, we calibrate the axes such that the only unknown parameter is g 2 using the semi-classical analysis. Then, we use the quantum fluctuation at the critical point to determine g 2 . Finally, we explain how we use this calibration to estimate the memory photon number in a given trajectory.
a. Semi-classical As shown by equation (S24) in the semi-classical approximation, at zero detuning, the critical point appears when | d g 2 | = κ a κ b /8 and from there the mean photon number in the memory increases linearly with the drive amplitude. Using these properties, we calibrate the drive amplitude axis in units of | d g 2 |: the x-intercept of the linear dependence at large photon number (semi-classical regime) is located at κ a κ b /8. The x-axis being calibrated, we linearly stretch the y-axis such that the asymptotic slope is 1 in the strong drive regime. According to eq. (S25), this transformation enforces the y-axis to |αg 2 | 2 and leads to the scaled data of Fig. S8.
This rescaling crucially depends on the values of κ a and κ b which are determined as follows. We measure the reflection coefficient of the memory in the presence of a pump tone slightly detuned from the frequency matching condition. This enables to capture the shift of parameters (frequency, internal losses, coupling losses) due to nonlinear effects arising from the pump while disabling the two-photon losses. From this measurement we numerically fit κ i a /2π ∈ [15,22] kHz and κ c a /2π ∈ [39, 42] kHz. The same protocol fails to determine precisely κ b due to the background induced by the band-stop filters and the strong dependence of the buffer parameters on the pump frequency. Instead, we use the diamond property derived in eq. (S27) that the top-right and bottom-left edge of  dm (x-axis), for increasing drive amplitude. On these panels, the two-photon drive LO frequency is set to ω LO dm = ωa + 100 MHz. When the two-to-one photon exchange transition is resonant, the engineered two-photon drive populates the memory. The average occupation of the memory is determined thanks to an undercoupled port via heterodyne detection. (d) radiated energy from the memory in units of circulating photon number (color) as a function of the pump and drive detuning from the frequency matching condition ∆a (x-axis), and the drive detuning from the buffer ∆ b (y-axis). In these coordinates, the feature takes the shape of a regular diamond.
the diamond have a slope of −κ b /κ a . We find κ b /2π in the range [13,20] MHz.
We later propagate the parameter range found for κ a and κ b on the rest of the calibration to give a robust confidence interval for g 2 andn.
We can independently check the calibration of |αg 2 | 2 by studying the excess internal losses arising from the two-photon dissipation. When the two-photon dissipation becomes resonant, the internal losses of the memory measured by direct spectroscopy increase drastically and become non-linear as a function of the probe power. The effective internal losses of the memory write where κ i a is the bare internal losses of the cavity and |α| 2 is the average circulating photon number due to the spectroscopy tone. The excess losses rewrite 8|αg 2 | 2 /κ b and provides an independent calibration of |αg 2 | 2 that we find in good agreement with the previous method.
b. Quantum signature At the critical point, the semi-classical analysis fails to capture the curvature of the mean photon numbern as a function of the drive amplitude d (Fig. 2). This curvature results from the quantum fluctuations at the dissipative phase transition [28]. Instead, we perform a quantum analysis and compute the average photon number in the steady state ρ ∞ of the Lindblad equation generated by (S16),n = Tr(ρ ∞ a † a) using the steadystate function imported from the QuTiP python package. Once we express |αg 2 | 2 as a function of | d g 2 | (see Fig. S8) the only fitting parameter is g 2 . Given the range of κ a and κ b , we estimate g 2 /2π ∈ [30,46] kHz.
c. Trajectory calibration We analyse the bit-flip time scale over several orders of magnitude, hence we increase the integration time T m to keep manageable amount of data in long bit-flip traces. Thanks to the previous calibration, we can readily get the memory photon number from the (I, Q) statistics of the trace. Indeed, from eq. S40, we have both G from the value of I 2 + Q 2 vac and 2Gκ c a η from the calibration ofn. When given a trace with different integration time T m , we determinen FIG. S8. Two photon coupling calibration. (a) radiated energy from the memory (y-axis) as a function of drive amplitude (x-axis). There are two regimes: when the drive amplitude is small, single photon loss overcome the two photon drive, and the memory stays in the vacuum. Passed the critical point, the memory gets populated by a coherent state with photon number asymptotically proportional to the drive amplitude. The axes units are chosen so that the critical point is at κaκ b /8 and the asymptotic slope is 1. The data correspond to an integration time of 10 µs with 10000 averages (crosses). The semi-classical model (green solid line) captures the position of the critical point but fails to explain the curvature of the experimental data. A full numerical simulation is used to reproduce the data where the only fitting parameter is g2 (red lines). (b) Zoom on the curvature around the critical point (gray rectangle from (a) ), emphasizing the agreement between simulations and experimental data.

C. Quantum detection efficiency
A different route leading ton is to measure the detection efficiency η. Indeed, for a coherent state containing a number of photonsn, the measured mean I and standard deviation σ(I) of the I-quadrature verifȳ Inversely, from the calibration ofn from the previous section, we estimate η 3%.
In this section, we independently evaluate η by fabricating a device containing a memory mode coupled to a transmon that serves as an in-situ measurement ofn. The fabricated memory mode is identical to the one described in the main text, and the entire two-photon exchange apparatus is replaced with a transmon qubit (see Fig. S10). The chip was mounted in a similar sample holder and measured with an identical wiring as the experiment described in the main text. The characteristics of the chip are listed in Table. III. T1 19.3 µs T2 24.3 µs κ c a /2π 38 kHz κ i a /2π 17 kHz χ/2π 1.75 MHz   TABLE III. Parameters of the device used to calibrate the quantum detection efficiency η. The transmon qubit lifetime and coherence times are denoted T1 and T2. The memory coupling and internal loss rates are denoted κ c a and κ i a , and χ corresponds to the dispersive coupling rate between the transmon and the memory.
Our evaluation of η follows three steps. First, we perform a standard spectroscopy in reflection of the memory mode in order to emulate a measurement signal that is directly proportional to the intra-cavity field amplitude a . Second, for a given amplitude a in , we calibrate the cavity photon numbern = a † a by resolving the photon number splitting of the qubit. Third, for each calibrated photon number we measure the fluctuations of the outgoing field a out and retrieve η by inverting Eq. (S44). We detail each step of this procedure below. FIG. S10. False-color optical micrograph of the detection efficiency chip in Nb (grey) on Si (dark blue). The memory resonator (blue) is capacitively coupled to a transmon (green).The inset is centered on the transmon and its Josephson junction in Al (light grey). We can address the memory and collect the reflected signal (blue waves) via the bottom 50Ω port. The left 50Ω port is dedicated to drive the transmon (green waves). This sample was also used to measure the memory thermal population of about 1%.
a. Memory spectroscopy For various incoming signal amplitudes S in , we perform a spectroscopy measurement recording the reflected signal S out . Using the results of the resonance fit, we can then translate the data in the (I, Q) plane in order to emulate a transmission signal: S t = A a , where A is an unknown scaling factor to be calibrated.
b. Photon number resolved qubit spectroscopy For various resonant signal amplitudes S in we activate a drive on the transmon at a fixed amplitude S q with a varying detuning ∆ q . The data S t (∆ q , S in , S q ) are then fitted to the result of a numerical simulation that we detail in the following. Using the steadystate function of the QuTiP package [40, 41], we solve the following dynamics where a (resp. q) is the memory (resp. qubit) mode annihilation operator, κ 1 = 1 T1 , κ φ = 1 T2 − 1 2T1 , Ω a (resp. Ω q ) is the drive on the memory (resp. qubit) and the remaining parameters are defined in table III. From this simulation we extract a (∆ q , Ω a , Ω q ), that is used to fit the dataset S t (∆ q , S in , S q ), where the fit parameters are the proportionality constants relating S in to Ω a , S q to Ω q and S t to a (see Fig. S11). c. Output field statistics For every drive amplitude S in , the previous fit estimates the intra-cavity field a , and hencen. By acquiring histograms of the output field S t , we now invert Eq. (S44) and retrieve η 7%, a factor two larger than the previously estimated value. This deviation can be attributed to differences in the RF connections of the two samples. These values may be explained by lossy elements (two circulators, one Eccosorb filter and two directional couplers) between the sample and the TWPA.

S5. BIT-FLIP TIME SIMULATIONS
We numerically simulate the dynamics of the memory described in Eq. (1) using the mesolve function imported from the QuTiP python package [40,41]. We run the simulation for three different values of g 2 (or equivalently κ 2 ). For each of these values, we sweep 2 in order to varȳ n = |α| 2 in the range of 4 to 40 photons. We initialize the memory in the coherent state |α , and fit the expectation value of the annihilation operator a to an exponentially decaying function. The extracted decay time corresponds to the bit-flip time. In Fig. S12, we display the computed bit-flip time as a function of the productn × (g 2 /2π) 2 , since it is a well calibrated quantity in our experiment. The data lie in the vicinity of the simulation results for g 2 /2π = 39 kHz, thus confirming our calibration of g 2 .