Dissipative stabilization of squeezing beyond 3 dB in a microwave mode

While a propagating state of light can be generated with arbitrary squeezing by pumping a parametric resonator, the intra-resonator state is limited to 3 dB of squeezing. Here, we implement a reservoir engineering method to surpass this limit using superconducting circuits. Two-tone pumping of a three-wave-mixing element implements an effective coupling to a squeezed bath which stabilizes a squeezed state inside the resonator. Using an ancillary superconducting qubit as a probe allows us to perform a direct Wigner tomography of the intra-resonator state. The raw measurement provides a lower bound on the squeezing at about (6.7± 0.2) dB below the zero-point level. Further, we show how to correct for resonator evolution during the Wigner tomography and obtain a squeezing as high as (8.2± 0.8) dB. Moreover, this level of squeezing is achieved with a purity of (−0.4± 0.4) dB.

While a propagating state of light can be generated with arbitrary squeezing by pumping a parametric resonator, the intra-resonator state is limited to 3 dB of squeezing. Here, we implement a reservoir engineering method to surpass this limit using superconducting circuits. Two-tone pumping of a three-wave-mixing element implements an effective coupling to a squeezed bath which stabilizes a squeezed state inside the resonator. Using an ancillary superconducting qubit as a probe allows us to perform a direct Wigner tomography of the intra-resonator state. The raw measurement provides a lower bound on the squeezing at about (6.7 ± 0.2) dB below the zero-point level. Further, we show how to correct for resonator evolution during the Wigner tomography and obtain a squeezing as high as (8.2 ± 0.8) dB. Moreover, this level of squeezing is achieved with a purity of (−0.4 ± 0.4) dB.

I. INTRODUCTION
One of the most striking predictions of quantum mechanics is that even in the ground state of an harmonic oscillator, any quadrature measurement is noisy. Zero point fluctuations can however be engineered and lowered for one quadrature of the field at the expense of the other. These squeezed states have become a central resource for quantum information processing. They can be used to boost the sensitivity of many measurements including gravitational wave detection [1][2][3][4], perform quantum secure communication [5,6] and used for measurement-based continuous-variable quantum computing [6,7]. Squeezing is usually generated by parametrically pumping a resonator. This process generates squeezing of both the intra-resonator and outgoing fields. While any amount of squeezing can theoretically be obtained for the outgoing field, the steady-state intraresonator squeezing is limited to 3 dB below the zero point fluctuations.
Intra-resonator squeezing beyond 3 dB can in principle be attained by injecting squeezed light into the resonator input using an external source of squeezed radiation [8][9][10]. In practice however, the achievable squeezing in such schemes is limited by losses associated with transporting and injecting the extremely fragile squeezed state into the resonator. A more attractive approach is to use reservoirengineering techniques [11], where tailored driving results in the cavity being coupled to effective squeezed dissipation [12][13][14]. These methods can also surpass the 3 dB limit, and do not involve transporting an externallyprepared squeezed state. Reservoir-engineering intracavity squeezing beyond 3 dB has recently been achieved for mechanical modes, both in optomechanical systems [15][16][17][18] as well as in a trapped ion platform [19].
In this work, we experimentally demonstrate that reservoir-engineering squeezing beyond 3 dB can also be achieved for purely electromagnetic intracavity modes, namely a microwave-frequency mode in a superconducting quantum circuit. Using the well developed circuit-QED toolbox, we also perform a direct tomography of the intra-resonator squeezed state instead of inferring the resonator state from the measured output mode. This is achieved through the use of an ancillary superconducting qubit, which enables in-situ Wigner tomography of the squeezed intracavity microwave mode. The intracavity squeezing factor reaches at least (−6.7 ± 0.2) dB, going well beyond the 3 dB limit. We also probe the non-classicality of the squeezed state by investigating its photon number statistics [20], and use our tomographic method to carefully study the full dynamics of the dissipative generation of squeezing. This work thus presents an interesting platform to stabilize, manipulate and characterize Gaussian states in-situ. Our stabilization technique could also be extended beyond simple squeezed states to other continuous variable states such as cat or grid states [21][22][23][24][25][26] by taking advantage of the large nonlinearities that can be engineered in circuit-QED.

II. SYSTEM AND MODEL
Our device consists in a Josephson Ring Modulator (JRM) [30] coupling one mode (the cavity) which we would like to stabilize in a squeezed state, and a second auxiliary mode strongly coupled to a transmission line (the dump). The cavity and dump have resonant frequencies ω c /2π = 3.741 55 GHz and ω d /2π = 11.382 GHz and decay rates κ c /2π = 40 kHz and κ d /2π = 8 MHz. Our setup also has an ancillary transmon qubit coupled to the cavity; its only role is to perform intra-resonator Wigner tomography ( Fig. 1.a).
When applying a pump at frequency ω − = ω d −ω c , and within the rotating-wave approximation (RWA) and stiff pump condition, the JRM leads to a beam-splitter interaction HamiltonianĤ − / = g −d †ĉ + g * −dĉ † , where the pump amplitude controls the coupling strength g − between the cavity and dump modes described by bosonic operatorsĉ andd. It mediates coherent exchange of photons between the cavity and the dump and thus lossless The dump mode is strongly coupled to a cold transmission line through which the JRM is pumped at both frequencies ω+ = ωc+ω d (two-mode squeezing) and ω− = ω d −ωc (photon conversion). A squeezed vacuum state is stabilized into the cavity as a result. An ancillary qubit with an ancillary readout resonator is used as a Wigner tomograph. The contours of the Wigner functions of each mode are shown as colored regions in the quadrature phase space, while a dashed circle represents the vacuum state. b) Frequencies of the involved modes and drives. c) Pulse sequence. The sum pump at ω+ with amplitude g+ and the difference pump at ω− with amplitude g− are applied for a time ts. After a waiting time tw, the Wigner function of the cavity W (α) is measured using a cavity displacement by −α followed by a parity measurement [27][28][29].
frequency conversion [31,32]. In contrast, a pump applied at frequency ω + = ω d + ω c mediates a parametric down conversion process involving cavity and dump, The pump amplitude controls the coupling strength g + . On its own, this kind of pumping leads to phase-preserving amplification [30,33] and generation of two-mode squeezed states [34]. Note that in order to avoid parasitic nonlinear effects, we operate the JRM at a flux point which maximizes these threewave mixing terms while cancelling the four-wave mixing terms [35,36].
Simultaneously pumping at these two frequencies enables various interesting phenomena such as effective ultrastrong coupling [37,38] or directional amplification [39][40][41]. Here, using a long-lived cavity mode, we show that this double pumping scheme can stabilize a squeezed state [12,13]. Indeed, in the rotating frame, and setting the phase references such that g ± are positive, the total Hamiltonian readŝ In the case where g + < g − , this Hamiltonian can be reinterpreted as a beam splitter interaction between the dump mode and a Bogoliubov modeβ = cosh(r)ĉ + sinh(r)ĉ † with r = tanh −1 (g + /g − ). It readŝ where the coupling strength is In the ideal case where the coupling rate κ d of the dump mode to a reservoir at zero temperature is much larger than any other rates, and where the cavity lifetime κ −1 c is unlimited, the Hamiltonian leads to the relaxation of the Bogoliubov mode into its ground state. In that state, the cavity mode is a vacuum squeezed state with squeezing parameter r = tanh −1 (g + /g − ).
The signature of this squeezing is best seen in the quadrature phase space of the cavity mode. We denote X − and X + the quadratures of the cavity mode that have the smallest and largest variances in a given state. In the vacuum state of the cavity (r = 0), the variance of the quadratures corresponds to the zero point fluctuations X 2 ± |0 = X 2 0 . The squeezing factor one can generate in the ground state of the Bogoliubov mode is simply a scaling of the variances by the factor S ± = X 2 ± /X 2 0 = e ±2r . In the general case, where the Bogoliubov mode is not cooled down to its ground state, these factors become [13] S ± = e ±2r (β ∓ β † ) 2 . ( We thus see that in principle, the 3 dB squeezing limit can be surpassed arbitrarily by having g + approach g − from below (as this causes the squeezing parameter r to diverge). However, in this limit the effective coupling rate G of the Bogoliubov mode to the dump goes down to zero. As a result, the competition between this engineered decay channel and the intrinsic cavity loss (rate κ c ) prevents the Bogoliubov mode from reaching its ground state. This both degrades the effective squeezing of the steady state, as well as its purity. Thus, for any value of g − there exists an optimum value of g + that minimizes the variance X 2 − . This minimum increases with the value of g − and is finally expected to saturate to a level set by the damping rates S − ≥ κ c /(κ c + κ d ), which reflects the fact that the damping rate of the dump κ d sets an upper limit to the coupling of the bosonic mode to the effective squeezed reservoir.
We thus see that a prerequisite for achieving squeezing well beyond 3 dB is to engineer a large ratio κ d /κ c . For our sample parameters, we have κ d 200κ c , leading to a lower bound of S − ≥ −23 dB [13]. into account the thermal equilibrium occupancies n th c and n th d of the cavity and dump modes, and in the limit of the experiment where κ c , G κ d , Eq. (3) leads to (see Appendix C or [18] for formula without approximation) where we introduce κ eff = κ c + 4G 2 /κ d and Γ ± = 4(g − ± g + ) 2 /κ d . Eq. (4) also makes it clear that the intrinsic loss rate κ c and non-zero environmental temperatures also lower the purity of the steady state P = Tr ρ 2 below 1, where ρ is the steady state cavity density matrix. This follows from the fact that P = 1/ S − S + for a Gaussian state.
To measure the squeezing and anti-squeezing factors S ± , we perform a full in-situ Wigner tomography [27][28][29] using an ancillary transmon qubit at frequency ω q /2π = 4.327 31 GHz ( Fig. 1.a). It couples dispersively to the cavity with a dispersive shift χ/2π = −3.28 MHz. A third resonator, at frequency ω r /2π = 6.293 GHz, is used to perform single-shot readout of the qubit state with a fidelity of 96 % in a 380 ns integration time. From the Wigner function, we compute the covariance matrix of the cavity mode quadratures and diagonalize it to extract the minimum and maximum cavity quadrature variances X 2 ± . Due to its coupling to the qubit, the cavity acquires an induced parasitic self-Kerr nonlinearity −Kĉ † 2ĉ2 and a qubit-state-dependent self-Kerr −K eĉ † 2ĉ2 |e e| with K/2π = 20 kHz and K e /2π = 70 kHz (measured in a previous run of the experiment). These non-linearities distort the squeezed state and thus reduce the effective squeezing factor, similarly to what occurs for Josephson parametric amplifiers (JPA) [42]. While no analytical solution taking into account the Kerr effects exists, Eq. (3) and Eq. (4) still provide a good description when G K. In the future, these nonlinearities could be harnessed as a resource to stabilize more complex non-Gaussian states [21,22,43,44].

III. STEADY-STATE SQUEEZING
The key advantage of reservoir engineering is that the desired target state is prepared in the steady state, independent of the initial cavity state: one can simply turn on the pumps and wait. We thus turn on g + and g − for a duration t s = 4 µs (cf Fig. 1.c) that is long enough to establish a steady state, and immediately afterwards measure the Wigner function W (α). To perform the measurement at each amplitude α, we start by applying a calibrated displacement D(−α) to the cavity state using a cavity drive at ω c with a pulse shape chosen to be a 13 ns wide hyperbolic secant and whose complex amplitude is proportional to −α. We then measure the cavity parity operator by reading out the qubit state after performing two π/2 unconditional pulses on the qubit separated by a waiting time π/χ = 152 ns. We perform phase-cycling, running each sequence twice with an opposite phase for the second π/2 pulse, so as to remove most of the parasitic contribution of higher order Kerr effects [45]. The Wigner function is probed on a discretized phase space using a rectangular grid of 25 pixels × 25 pixels approximately aligned to the squeezing axis. Due to the finite window size (cf Appendix E), we could only resolve antisqueezing up to 11 dB (dotted dash line in Fig. 2.b). Each Wigner tomogram is averaged over 5000 realizations. To increase the repetition rate and limit the low-frequency drifts, the cavity is first emptied by applying a difference pump g − , cooling it down to a thermal vacuum state with residual population n th c = n th d = 0.017 ± 0.003 (see Appendix B 3). The qubit is also reset to its ground state using measurement-based feedback. Furthermore, to minimize the low-frequency noise as much as possible, we interleave pump-on-measurements with pump-offmeasurement. We thus obtain experimental squeezing factors S ± = X 2 ± / X 2 off · X 2 off /X 2 0 by first normalizing the measured variances X 2 ± with the measured pumps-off variances X 2 off and then correcting for the thermal occupancy X 2 off /X 2 0 = (0.15 ± 0.03) dB. The pumping strengths g + and g − are calibrated using independent measurements (Appendix B 1). We estimate a statistical uncertainty of ±0.2 dB on the variances extracted from the measured Wigner functions.
The obtained steady-state squeezing and antisqueezing factors are displayed in Fig. 2.a) as a function of g + and g − . We observe a maximum squeezing of (−6.7 ± 0.2) dB well below the −3 dB limit, which we believe to be the highest squeezing factor observed in an intracavity microwave mode. Correspondingly, we extract an anti-squeezing of (7.7 ± 0.2) dB and thus a state purity of (−0.5 ± 0.2) dB. For each value of the rate g + , the largest squeezing we observe occurs for g − close to g + . This trend is expected from the analytical Kerr-free model Eq. (4) of the system, which predicts the working point of largest squeezing as a function of g + (green dashed line in Fig. 2.a) [13]. However, contrary to the expected monotonic increase of the optimal squeezing factor S − with g + (Kerr-free prediction in bottom panels of Fig. 2.a), we find a global maximum squeezing at a finite value of (g − , g + ). Note that for g + > g − the system becomes unstable: the qubit gets ionized [46], preventing us from measuring the Wigner functions (grey shade area).
In Fig. 2.b-c), we show the squeezing and antisqueezing as a function of g + /g − as well as some measured Wigner tomograms, for g − /2π = 1.85 MHz. For g + < 0.7g − , the measured variances are well captured by Eq. (4) with exponentially increasing squeezing factors. For g + > 0.7g − , the measured variances start deviating from the theory (solid lines in Fig. 2.b). As can be seen in Fig. 2.c), the squeezed states are not Gaussian anymore in this parameter region. The Wigner functions develop an S-shape, a typical signature of the cavity self-Kerr. We attribute this effect to higher order terms we have so far neglected: the self-Kerr rate induced by the qubit on the cavity, as well as a residual four-wave mixing term in the JRM Hamiltonian (see Appendix D 2).
While the raw measurement of the Wigner function provides a good estimate of the steady-state squeezing parameter (circles in Fig. 2.b), the finite measurement time needed to perform tomography leads to a systematic error. During this finite measurement time, the pump tones g ± are off, implying that the cavity is no longer coupled to an effective squeezed reservoir. The squeezed state thus degrades due to the intrinsic cavity loss. A further error is caused by evolution under the cavity-self Kerr nonlinearity during this time. Both these effects cause our Wigner function methods to underestimate the true value of the steady-state squeezing.
It is possible to correct for this measurement error and retro-predict via numerical simulation the squeezing fac-tors S ± associated with the state prepared at the end of the stabilization period [47]. To that end, we consider a series of input model Gaussian states for which we numerically implement our experimental Wigner tomography measurement. At the end of these simulations, we obtain a mapping from Gaussian states to measured squeezing and anti-squeezing factors that we are able to invert in order to retro-predict the stabilized state (Appendix D 1). Using this correction improves the best squeezing estimate to (−8.2 ± 0.8) dB (dots with error bars in Fig. 2 It is interesting to compare our stabilization technique to other intra-resonator microwave squeezing generation schemes. One possibility consists in driving a cavity with a squeezed input state that is externally generated by a Josephson parametric amplifier (JPA) [10,48]. High squeezing factors [42,[48][49][50] ( −10 dB) can be achieved in the amplifier output field. However, transferring this state into a cavity is challenging as it is extremely sensitive to microwave losses, resulting in degraded squeezing and purity. For comparison, we consider a resonator driven by a pure squeezing source (in practice a JPA). To achieve the same intracavity squeezing and purity as our setup (S − = −8.2 dB and P = −0.4 dB respectively), the source would need to generate an output squeezing better than ∼−9.1 dB and the losses between the source and the resonator would need to be kept below 0.15 dB. This level of loss is smaller than the typical insertion loss of common microwave components. It is hard to achieve, even if all elements are fabricated in a single-chip architecture. We can also compare against another approach for generating (but not stabilizing) a squeezed state, based on the use of arbitrary state preparation techniques (e.g. the SNAP gate protocol [51]). In Ref. [52], the authors used such an approach to obtain a squeezing factor of −5.71 dB for a purity of −0.86 dB with a non-deterministic success rate of 15 %.

IV. NONCLASSICAL PHOTON DISTRIBUTION
One of the hallmarks of vacuum squeezed states is that they are quantum superpositions involving only evennumber photon Fock states. Such ideal states have the Our device allows us to directly verify this unique, non-classical aspect of the squeezed states we stabilize in our cavity [20]. This is because our coupling to the ancilla qubit is strong enough to place us in photon-number-resolved regime where distinct cavity photon numbers can be resolved by measuring the effective qubit frequency, i.e. χ Γ 2 with Γ 2 = (11 µs) −1 the qubit coherence rate.
After preparing the squeezed state, we perform spectroscopy of the qubit using a narrow-bandwidth π-pulse at a varying probe frequency ω followed by qubit readout (green curve in Fig. 3.a). The observed peak heights at each frequency ω − ω q ≈ nχ allow us to determine the cavity photon-number distribution P(n) [53]. We correct this dataset for the qubit residual thermal population (1 %), the finite fidelity of the π-pulse and readout errors. Interestingly, the peaks are not evenly spaced in frequency due to the higher nonlinear term −K eĉ † 2ĉ2 |e e| with K e /2π = 70 kHz. The photon number distribution P(n) are then obtained from the qubit excitation probability at ω q − n(χ + 2K e n) (vertical lines). For comparison, we also measure the photon number distribution P(n) for two other cavity states: a thermal equilibrium state when no pumps are applied (blue in Fig. 3) and a thermal state that we create by only applying a sum pump g + /2π = 0.43 MHz (orange curve). Note that this thermal state is obtained by tracing out the dump mode for the vacuum two mode squeezed state that is stabilized between cavity and dump [54].
For the squeezed state (green curve with g − /2π = 2.2 MHz, g + /2π = 1.42 MHz), we observe a nonmonotonic behavior: the weight of even photon numbers is enhanced, whereas that of odd photon numbers is suppressed (note the log scale here). The non-zero but small population of odd Fock states indicates a deviation from an ideal squeezed vacuum state. The measured data closely fits to our numerical simulation (dots in Fig. 3.a). The measurement done with the pumps off gives the thermal population of the cavity n th c = 0.017 (blue dots) and also indicates the measurement noise floor. For the thermal state, we observe a Bose-Einstein distribution with a population n th c = 1.5 (orange dots). Even though their Wigner function is always positive, squeezed states are typically regarded as non-classical states as they cannot be represented as statistical mixtures of coherent states. Formally, this means that they do not have well-behaved Glauber-Sudarshan P representations [55]. Equivalently, it also manifests itself in the behaviour of so-called Klyshko numbers K n = (n + 1)P(n − 1)P(n + 1)/nP(n) 2 . A state is non-classical if for one or more integers n, K n < 1 (as this implies that the P function cannot be well behaved). For example, a perfect squeezed vacuum state, as it only includes even photon numbers, exhibits infinite odd Klyshko numbers and zero even Klyshko numbers. Ref. [20] computed the Klyshko numbers for a squeezed state generated by an external JPA and observed a Klyshko number smaller than 1, even though they only observed monotonic behavior in the photon distribution P(n).
For a thermal state, the Klyshko number are given by K th n = (n + 1)/n independently of temperature (orange bars in Fig. 3.b). We observe this universal relation with the prepared thermal state (orange points with errorbar). Interestingly, it is a striking demonstration of the fact that a two mode squeezed state generates a thermal distribution when tracing out one of the modes. It is expected from the maximally entangled state at a given average energy. We do not show the Klyshko numbers when the pumps are off because P(n) is below the noise floor.
For the squeezed state, we observe ample oscillations in the Klyshko numbers (notice the log scale again). We measure K 2 = 0.23 and K 4 = 0.5 that are well below one (green points with errorbars). Similarly to the cavity population P(n), our numerical model (green bars) reproduces the observed Klyshko numbers.

V. STABILIZATION DYNAMICS AND DECAY OF SQUEEZING
Our measurements establish that, as expected, the reservoir engineering scheme we implement is able to stabilize a squeezed state in the cavity. In addition to characterizing the steady state, it is also interesting to ask how long the scheme takes to prepare the steady state. For the ideal (Kerr-free) system, and in the limit of a large dump-mode damping, one can use adiabatic elimination to show that this preparation timescale is κ d /G 2 [13].
We can directly test this prediction in our experiment. The measured squeezing and anti-squeezing factors are shown in Fig. 4.a) as a function of the time t s during which the pumps are turned on for g + /2π = 1.16 MHz and for various values of g − . By normalizing the squeezing and anti-squeezing factors S ± by their steady-state values S ss ± , we observe, as expected, that the steady-state is reached in a typical time of κ d /G 2 that decreases with g − (rectangles in Fig. 4.a). As a consequence, the stabilization time increases with squeezing when considering a fixed g + value, as long as G dominates both the cavity loss rate κ c and its self-Kerr rate K. This is well- understood from the cooling dynamics of the Bogoliubov mode: larger squeezing parameters r = tanh −1 (g + /g − ) are obtained for smaller values of G = g 2 − − g 2 + but they lead to a longer relaxation time. The evolution of the squeezing factors is reproduced using numerical simulation of the master equation (solid lines in Fig. 4.a).
It is also interesting to examine experimentally the time-evolution of the full cavity Wigner functions. In Fig. 4.b), the evolution at g + /2π = 1.16 MHz and g − /2π = 1.85 MHz (global minimum of the squeezing factor) shows how the squeezing establishes with some rotation and distortion of the Gaussian distribution due to Kerr effect as the average number of photons gets larger.
The steady-state is thus reached in about κ d /G 2 but how fast does it disappear once the pumps are turned off? Operating at g + /2π = 1.42 MHz and g − /2π = 2.6 MHz, we perform a Wigner tomography and compute the squeezing and anti-squeezing after a waiting time t w (Fig. 5.a). A fast decrease of the squeezing factor is observed in a characteristic time shorter than the cavity relaxation time κ −1 c . We attribute this deviation from the behavior expected of a perfectly harmonic oscillator (dashed lines) to the self-Kerr effect induced by the trans-mon qubit onto the cavity. The corresponding predicted evolution of squeezing factors is shown with K/κ c = 0.5 as solid lines in Fig. 5.a). In numerical simulations, we observe a transition from over-damped to under-damped oscillations of the squeezing factor S − as K/κ c increases beyond about 1 (Appendix D 3). Since K κ c in the experiment, we are close to a critical damping regime.

VI. CONCLUSION
Using dissipation engineering, we have shown the stabilization of a squeezed state in a microwave resonator with a squeezing factor greatly exceeding the standard 3 dB limit for coherent in-situ parametric pumping. We directly measure the squeezing factor by performing a direct Wigner tomography using an ancillary qubit. Correcting for state evolution during measurement, we infer that we achieve a squeezing factor of (−8.2 ± 0.8) dB. While reservoir-engineered squeezing of mechanical modes has previously been demonstrated, this is the first demonstration of this method (to our knowledge) in an electromagnetic system. The reservoir engineering technique used here thus extends the state-of-the-art for intra-resonator microwave squeezing. Moreover, the produced squeezed state is close to a pure state with purity of (−0.4 ± 0.4) dB. A displaced vacuum squeezed state could also be stabilized in our system by adding a coherent drive on the dump.
Beyond the stabilization of Gaussian squeezed states, the techniques presented here could be useful for the stabilization of far more complex states. As discussed, Kerr nonlinearities already play an appreciable role in our experiment. Future work could use this nonlinearity directly as a resource for non-Gaussian state preparation. Recent work has demonstrated that the combination of squeezing-via-parametric driving with Kerr interactions can be used to generate cat states [21,43] and even entangled cat states [44]. The combination of dissipative squeezing (as realized here) with Kerr interactions could similarly yield complex cat-like states. Our techniques could also be used to generate squeezed Fock states [56], squeezed Schrödinger's cat states [57] or for the preparation of grid states without the need for measurement [23][24][25][26]. These engineered squeezed states could find many applications. Indeed, used to erase which-path information, they can increase gate fidelity [58]; used to increase distinguishability, they can improve qubit state readout [9,14,59,60]. Squeezing can also be used in spin detection to enhance the light-matter coupling [61,62]. Finally, dissipative squeezing techniques employed on a single site of a lattice of microwave resonators (see e.g. Ref. [63]) can serve as a shortcut for effectively generating highly-entangled many-body states [64,65].

ACKNOWLEDGMENTS
We are grateful to Olivier Arcizet and Alexandre Blais for discussions. This work was initiated during a discussion that happened during Les Houches Summer School in July 2019. We acknowledge IARPA and Lincoln Labs for providing a Josephson Traveling-Wave Parametric Amplifier. The device was fabricated in the cleanrooms of Collège de France, ENS Paris, CEA Saclay, and Observatoire de Paris. This work is part of a project that has received funding from the European Union's Horizon 2020 research and innovation program under grant agreement No 820505. AC acknowledges support from the Air Force Office of Scientific Research MURI program, under Grant No. FA9550-19-1-0399.

Appendix A: Steady-state Wigners tomograms
The Wigners tomograms of all the points of Fig. 2 are available on [66].

Appendix B: Sample and setup
The sample is the same as in Ref. [36] albeit for a different cool-down. The measurement setup is also similar with the addition of the pump at the sum frequency (Fig. 6). The two local oscillators for the pumps are generated by mixing the output of the two microwave sources that are used to generate the dump and cavity drives. Intermediate frequency (IF) signals -tens of MHz -generated by the Quantum Machines' OPX hardware are upconverted by these local oscillators. Finally, we combine and amplify the two pumps before combining them to the dump port inside of the dilution refrigerator.
To successfully stabilize and measure a squeezed state on a well-defined squeezing axis (g ± real), a good phase coherence is required between the pumps and cavity drives. Our setup ensures this condition by deriving the pumps from the dump and cavity local oscillators. One difficulty of our experiment is the large power required for the pumps to reach maximal squeezing factor. This requires the use of a room-temperature amplifier after the mixers (Fig. 6). This amplifier has a slow temperatureinduced drift in gain, leading to a relative error of 2 % on the pump amplitudes (corresponding to the horizontal errorbars in Fig. 2.b).

Calibration of the pumps
This section shows how to relate the IF amplitudes A − and A + to the rates g − and g + .
To calibrate g − , we measure the mean photon number in the cavity after applying the pump when the cavity is initially populated with a coherent state α = √ 6. Depending on the amplitude A − and duration 4σ of the pump pulse, the rate at which the cavity coherently exchanges excitations with the dump varies (Fig. 7). Due to the large dissipation rate of the dump, the oscillations of the cavity mean photon number n are damped. By fitting the oscillations using a master equation, we find, as expected, a linear dependence of g − as a function of A − that we use as calibration.
To calibrate g + , we measure the mean photon number in the cavity n th c after applying a square pulse with amplitude A + for 100 ns when the cavity is initially in vacuum. The mean photon number is measured via cavityinduced Ramsey oscillations [36]. The only difference with the former reference is that the distribution of photon numbers is thermal instead of Poissonian. Hence, the phase acquired by the qubit during the waiting time of the Ramsey sequence differs and leads to a final qubit excitation probability of P e (t) = n th c (1 − cos χt) + 1 2(1 − cos χt)(n th c + 1)n th c + 1 Using a time dependent master equation, the mean photon number is converted into a two-mode squeezing rate g + . The curve g + as function of A + is non linear (Fig. 8), likely due to higher order non-linearities in the Hamiltonian. The calibration g + (A + ) is then obtained by interpolating the measurement.

Cavity displacement calibration
The calibration of the displacement of the cavity under a pulsed coherent drive is performed by counting the mean photon number. The method chosen to count the mean photon number is to use the ancillary qubit and readout as a vacuum detector [36]. This method also allows us to extract the cavity decay rate.

Cavity thermal population
The cavity thermal population is extracted from the cavity-induced Ramsey oscillations of the ancillary qubit [36]. With the reset protocol, consisting of a swap pulse (g − ) between cavity and dump modes followed by a measurement-based feedback initialisation of the qubit in its ground state, we measured a mean photon number n th c = 1.7 ± 0.3 · 10 −2 corresponding to an effective temperature of (44 ± 2) mK for the cavity.

Correction and uncertainty on the quadrature variances
We wish to extract the squeezing and anti-squeezing factors by normalizing the measured variances to the  FIG. 6. Schematic of the measurement setup. The rf sources color refers to the frequency of the matching element in the device up to a modulation frequency. Multiple instances of a microwave source with the same color represent a single instrument with split outputs. The sum (blue) and difference (red) pumps are obtained by mixing the cavity and dump rf sources to ensure phase stability. The TWPA [67] was provided by Lincoln Labs. zero-point fluctuations. However, the residual thermal population offsets the measured value of the zero-pointfluctuations by a factor 2n th + 1. This also means that all the measured squeezing factors have to be offset by (0.15 ± 0.03) dB. Due to other sources of uncertainty, such as fluctuations on the cavity displacement pulses, we measure a higher statistical uncertainty for the pump-off variances of ±0.2 dB.
Appendix C: Kerr-free analytical model This derivation, which can be found in Ref. [18], is given here for completeness. When continuously pumping at the difference and sum of the resonance frequencies with rates g − and g + , the Langevin equations reaḋ where the cavity (dump) input field operatorsĉ in (d in ) Solving the Langevin equations for the steady-state, the squeezing S − and anti-squeezing S + factors are given by .

Appendix D: Modeling the Kerr effect
The Kerr effect is not included in the analytical model described in Appendix C. It induces spurious effects, which reduce the maximal squeezing factor and accelerate the relaxation of squeezing. In this section, we show how to take these effects into account. We simulate our system using the QuantumOptics.jl library [68]. The steady-state simulations are run on an Nvidia Geforce 1080Ti GPU, which allows us to reach Hilbert space dimensions of about 1800. All of the other simulations are run on the CPU. Except for the Wigner tomography retro-prediction, the qubit is not simulated but we take into account the Kerr effect it induces on the cavity. In the case of the Wigner tomography retro-prediction, the dump is adiabatically eliminated.

Retro-prediction of the Wigner tomography
In order to correct for the error introduced by the cavity evolution during Wigner tomography, we resort to simulations of the cavity and qubit alone. Indeed, in the absence of pumps, the effect of the JRM on the cavity is negligible. We numerically implement our experimental Wigner tomography pulse sequence on a truncated Hilbert with up to 50 excitations for the cavity and the two qubit states. Starting from a range of initial squeezed states for the cavity, with variances (S i − , S i + ), we simulate This data-set provides a function f W that maps actual variances (S i − , S i + ) of the pre-measured quantum state to the variances (S W − , S W + ) extracted from the measured Wigner tomograms. As this function empirically appears bijective, the retro-prediction is performed by interpolating its inverse f −1 W . The interpolated f −1 W for the initial squeezing and anti-squeezing, as well as the simulated points, are shown in Fig. 9.a and b respectively. The retro-predicted initial squeezing and anti-squeezing corresponding to Fig. 2.a-b are shown in Fig. 9.c-d respectively.
Assuming the ±0.2 dB of uncertainty on the measured S W ± , and retro-predicting the evolution during Wigner tomography, we obtain an uncertainty ∆S ± on the retropredicted squeezing and anti-squeezing factors that depends on the value of the measured squeezing and antisqueezing ( Fig. 9.e-f).

Steady-state simulations
As seen in Fig. 2.b, the analytical Kerr-free model fails to quantitatively describe the squeezing factor at g + > 0.7g − (Fig. 10.a). Here, we compute how higher or-der terms in the Hamiltonian may explain this difference. The first term we consider is the Kerr effect −Kc † 2 c 2 induced by the qubit on the cavity (Fig. 10.b). This simulation accurately predicts the optimal g + but still fails to reproduce the measured squeezing factors above g + /g − = 0.7. Experimentally, we aim for a JRM flux bias that maximizes the three-wave mixing term while cancelling the four-wave mixing term. However, small deviations from this sweet spot create four-wave mixing terms between the cavity, dump and pumps. Contrary to the retroprediction simulations which model a situation where the pumps are turned off, these extra terms may have a significant impact on the squeezing factor where the pumps are turned on. In the RWA, the four wave-mixing term leads to three kinds of interactions, a cross-Kerr between cavity and dump K cd c † cd † d, an AC-Stark frequency shift due to the pumps 2(|p − | 2 +|p + | 2 )(K pc c † c+K pd d † d) and parametric squeezing drive due to pump inter-modulation K pc p + p * − c † 2 + K pd p * + p * − d † 2 + h.c.. The JRM also induces a self-Kerr interaction for the dump, but we neglected it as it is one order of magnitude smaller than K pd p * + p * − in our case [69] and much smaller than the dissipation rate κ d anyway. The rates K cd , K pc and K pd are not measured in this run. Realistic values K cd /2π = 250 kHz, K pc |p − | 2 /2π = 172 kHz and K pd |p − | 2 /2π = 172 kHz can b) a) c) FIG . 10. a, b and c) Crosses, retro-predicted squeezing factors for g−/2π = 1.85 MHz as in Fig. 2.b. Shaded areas correspond to different models with 2 % uncertainty on g+ and g−; in a), Kerr-free analytical model, in b), steady-state simulations with K/2π = (20 ± 2) kHz, in c), steady state simulations including in addition to the Kerr effect some JRM four-wave mixing terms, K cd /2π = 250 kHz, Kpc|p−| 2 /2π = (172 ± 4) kHz and K pd |p−| 2 /2π = (172 ± 4) kHz.
change the squeezing factors at the large g + , which comforts the assumption that higher order nonlinearities may explain the deviations we observe between our analytical model and the measured squeezing factors.
To numerically compute the steady-state squeezing and anti-squeezing as a function of g − and g + , we use an iterative method to find the Liouvillian eigenvalues on a truncated Hilbert space comprising up to 60 excitations for the cavity and 30 excitations for the dump.

Simulations of the squeezing dynamics
The dynamics of stabilization and decay of squeezing are computed by solving the master equation on a truncated Hilbert space comprising up to 20 excitations for the cavity and 16 excitations for the dump.
To understand the effect of the self-Kerr term on the squeezing decay, we simulate the evolution of squeezing for varying waiting time t w and various self-Kerr rates K (Fig. 11). We initialize the cavity state at t w = 0 in a Gaussian state with the measured S − = −5.7 dB and S + = 6.2 dB of Fig. 5.a. Dismissing the Kerr effect (K = 0), we observe an exponential damping of squeezing due to the cavity relaxation. For nonzero K, the squeezing factor also oscillates in time. Our experimental value K/2π = 20 kHz is closed to the critically damped regime where the effective decay time is maximally reduced. This observation highlights the crucial role of Kerr effect in the imperfections of our Wigner tomography technique used to estimate the variances.