{'fraction_non_alphanumeric': 0.04152620933760324, 'fraction_numerical': 0.013336423394572729, 'mean_word_length': 4.406652238295933, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Emergence of nonlinear friction from quantum fluctuations

Nonlinear damping, a force of friction that depends on the amplitude of motion, plays an important role in many electrical, mechanical and even biological oscillators. In novel technologies such as carbon nanotubes, graphene membranes or superconducting resonators, the origin of nonlinear damping is sometimes unclear. This presents a problem, as the damping rate is a key figure of merit in the application of these systems to extremely precise sensors or quantum computers. Through measurements of a superconducting circuit, we show that nonlinear damping can emerge as a direct consequence of quantum fluctuations and the conservative nonlinearity of a Josephson junction. The phenomenon can be understood and visualized through the flow of quasi-probability in phase space, and accurately describes our experimental observations. Crucially, the effect is not restricted to superconducting circuits: we expect that quantum fluctuations or other sources of noise give rise to nonlinear damping in other systems with a similar conservative nonlinearity, such as nano-mechanical oscillators or even macroscopic systems.

Nonlinear damping, a force of friction that depends on the amplitude of motion, plays an important role in many electrical [1], mechanical [2] and even biological [3] oscillators. In novel technologies such as carbon nanotubes, graphene membranes [4] or superconducting resonators [5], the origin of nonlinear damping is sometimes unclear. This presents a problem, as the damping rate is a key figure of merit in the application of these systems to extremely precise sensors [6,7] or quantum computers [8]. Through measurements of a superconducting circuit, we show that nonlinear damping can emerge as a direct consequence of quantum fluctuations and the conservative nonlinearity of a Josephson junction. The phenomenon can be understood and visualized through the flow of quasi-probability in phase space, and accurately describes our experimental observations. Crucially, the effect is not restricted to superconducting circuits: we expect that quantum fluctuations or other sources of noise give rise to nonlinear damping in other systems with a similar conservative nonlinearity, such as nano-mechanical oscillators or even macroscopic systems.
Nonlinear damping, the change in damping rate with the amplitude of oscillations, is ubiquitous in nature. It was famously described mathematically by van der Pol [1] in the context of his work on vacuum tube circuits [9]. Now, it is used to describe the physics of a diverse set of systems, such as the rolling of ships in waves [2] or the nervous system [3]. It has attracted recent interest due to its appearance in novel experimental platforms such as nanoscale ferromagnets [10], superconducting circuits [5,[11][12][13] and nanoelectromechanical systems (NEMS) [14][15][16] made for example from carbon nanotubes, graphene [4,17] or superconducting metal [18]. In some of these systems the nonlinearity is well explained [19][20][21][22]. Most notably the saturation of two-level systems in the environment can cause the damping rate to decrease, as the power injected into the system increases [11][12][13]18]. But the origin of an increase in damping in certain NEMS [4,[14][15][16] or superconducting resonators [5] remains speculative. Understanding the origin of nonlinear damping in some of these systems is critical due to the importance of their energy damping rates in applications such as NEMS based mass sensing [6] or spectrometry [7], as well as quantum-limited amplification [23,24] in superconducting quantum computers [8].
We study nonlinear damping through the experimental platform of superconducting circuits [25]. Central to the observed physics is the conservative nonlinearity induced by a Josephson junction: that the resonance frequency varies with the oscillation amplitude, which can be further approximated as a Duffing or Kerr nonlinearity [26]. For this reason, the phenomena discussed here are applicable to all systems featuring a similar nonlinearity in their resonance frequency, ranging from the carbonnanotubes mentioned above, to the simple mechanical pendulum. We focus on the regime where this nonlinearity is small, as in Josephson parametric amplifiers [24], rather than the single-photon nonlinear regime used to construct artificial atoms in circuit quantum electrodynamics [26]. Because of their small nonlinearity, such systems are often thought to be completely described by the classical Kerr oscillator [5,27].
Here however, we report on a novel form of nonlinear damping originating purely from the quantum noise of a Kerr oscillator. The nonlinearity is experimentally characterized by probing the oscillator's frequency response. Our observations are accurately described by a model devoid of ad-hoc nonlinear damping, but which take into account the effect of quantum noise. Perturbative calculations based on this model demonstrate that the expected amplitude of oscillations obeys a classical Kerr equation with nonlinear damping. Finally, we provide an intuitive picture in which the phenomenon can be understood as the oscillator experiencing dephasing induced by its own photon shot noise.
The circuit used in this experiment ( Superconducting Kerr oscillator circuit to explore nonlinear damping. a, The Kerr oscillator is constructed from an inductor, capacitor and Josephson junction, and is side-coupled to a transmission line with a coupling rate κext. The circuit undergoes internal damping at a rate κint. b, Optical micrograph of the device, where light gray corresponds to superconducting molybdenum-rhenium, and dark gray to the insulating silicon substrate. An interdigitated capacitor on the right is connected to a meandering inductor on the left. The circuit couples to a transmission line (coplanar waveguide) at the top. c, Scanning electron micrograph of the SQUID: two aluminum/aluminum-oxide Josephson junctions connected in parallel. As the flux threading the SQUID is fixed, it behaves in this context as a single junction.
from an inductor, capacitor and superconducting quantum interference device, or SQUID. The SQUID is fluxbiased to its sweet spot (integer flux quantum), and behaves as a single Josephson junction [28]. The junction induces an anharmonicity K = 2π × 76 kHz five orders of magnitude smaller than the resonance frequency ω r = 2π × 5.17 GHz. The cosine potential of the junction is accurately described in this limit K ω r by the Kerr effect in the Hamiltonian [29] whereâ is the annihilation operator for photons in the circuit. Intuitively, the junction is acting as an inductor, with an inductance which increases with the number of photonsâ †â in the circuit. As a consequence, the resonance frequency of the circuit is lowered with each added photon, labeled as the Kerr term in Eq. (1).
The circuit undergoes internal damping, losing energy at a rate κ int = 2π × 190 kHz. This is typically due to losses in the different dielectric materials traversed by the electric fields [30]. Additionally, the circuit is coupled to a transmission line, through which we drive the circuit with a microwave signal. Conversely, the transmission line leads to energy leaking out of the circuit, which is Observation of a resonator response suggesting nonlinear damping. a, Measured transmission |S21| (dots) for different drive powers. Whilst the shift in resonance frequency is expected from a classical analysis of the Kerr oscillator using Eq. (2), Min|S21| is expected to remain constant (dashed line). b, Measured Min|S21| (dots) as a function of drive power. Eq. (2) yields Min|S21| = |1 − κext/(κint + κext)|, suggesting a damping rate which increases with power κint → κ nl int (|a|) from κ nl int = 2π × 190 kHz to κ nl int = 2π × 340 kHz. Indeed, adding nonlinear damping κ nl int (|a|) to Eq. (2) leads to theoretical predictions (solid lines in a) in good agreement with the data. characterized by an external damping rate κ ext = 2π×2.1 MHz. As a consequence, the total damping rate and spectral linewidth κ = κ int + κ ext is much larger than the shift in resonance frequency K due to an added photon: κ K. The circuit is thus far from the regime of superconducting qubits [26]. We will call it a Kerr oscillator and first attempt to describe its behavior following the classical equation for the steady-state amplitude of its oscillations a Here ∆ = ω r −ω d is the detuning of the driving frequency ω d to the resonance frequency ω r , and the strength of the drive = κ ext P in /2 ω r is given by P in the power of the drive impinging on the device.
The circuit is made by patterning a thin film of sputtered molybdenum-rhenium alloy on silicon, and subsequently fabricating the aluminum/aluminum-oxide tunnel junctions (see Supplementary Information). The device is thermally anchored to the ∼ 20 milliKelvin stage of a dilution refrigerator, and the input (output) microwave wiring is attenuated (isolated) to lower its microwave mode temperature, such that the number of photons n th excited by thermal energy is negligible. The transmission coefficient S 21 = 1 − κ ext a/(2 ) is then measured using a vector network analyzer (VNA) for varying microwave power (see Fig. 2(a)).
We note an increase in both the detuning ∆ min which minimizes transmission, and the value of the minimum Min|S 21 |. The classical prediction ∆ min = Kα 2 resulting from Eq. (2) -where α = 2 /κ is the expected maximum amplitude -accurately matches the shift of the resonance. However, by plugging the maximum amplitude α into the expression for S 21 , we obtain a constant value for Min|S 21 | = |1 − κ ext /κ| (dashed line in Fig. 2(a)), which disagrees with the measurement.
In a classical approach to the problem, a powerdependence of the internal damping rate therefore has to create this change. Since κ ext is determined by the geometry of the circuit, it should remain unchanged by the power of the drive. For κ ext /κ to vary and produce the observed change in Min|S 21 | = |1−κ ext /(κ ext +κ int )|, the internal damping should increase to 2π × 295 kHz as the power increases. Such nonlinear damping can be included in the model of Eq. (2) through κ int → κ nl int = κ int +γ|a| 2 . We fit a solution of the resulting equation to the data (see Supplementary Information), observing good agreement (Fig. 2).
Whilst providing an accurate model for our observations, adding ad-hoc nonlinear damping offers no explanation as to the physical mechanism underlying the effect. Usually, the most prominent source of nonlinear damping in superconducting circuits is the saturation of two-level systems (TLSs) in the environment [11][12][13]. However, with increasing driving power, the saturation of TLSs will result in a decrease the internal damping rate, whilst we observe the opposite. Here, we show that in our system nonlinear damping emerges purely from the internal quantum noise of the oscillator.
We first show that approaching the problem quantum mechanically, without adding nonlinear damping, perfectly describes our measurements. The effect of quantum noise is included in the model through the steadystate Lindblad equation whereρ is the density matrix describing the steady-state of the oscillator. By solving this equation for varying drive strengths and frequencies, the resulting amplitude â = Tr(âρ) is used to obtain S 21 . With only the circuit parameters as free variables, and notably a constant value for the internal damping, this model is fitted to all S 21 traces (see Supplementary Information) , revealing excellent agreement to the data (Fig. 3). Note that we recently became aware of an analytical solution to this Lindblad equation [31,32], which may have simplified our approach.
In Fig. 3 Beyond describing the data, this model analytically leads to a nonlinear damping equation for the expectation value â . As detailed in the Supplementary Information, we determine a perturbative solution to the quantum model. We find that when the drive detuning, Kerr effect and thermal noise have only a perturbative effect on the system, quantum and thermal noise lead to a nonlinear damping coefficient in an equation for â where the damping is given by κ + γ| â | 2 . The term [â,â † ] = 1 is the contribution of quantum noise and notably has the same effect as half a quantum of thermal noise n th = 1/2. The fact that the nonlinear damping model and the quantum model are The amplitude-dependent resonance frequency (Kerr effect) translates to a radius-dependent rotation around the origin. The center of mass of the distribution rotates at a frequency Kα 2 . In a frame rotating at that frequency, the effect of the Kerr nonlinearity is to increase the uncertainty in phase (this is the frame adopted in panel c). The larger the uncertainty in phase (the extreme case being a ring around the origin), the closer the center of mass of the distribution gets to the origin (i.e. | â | → 0). This is the first contribution (effect A) to a reduced amplitude | â |. c, Wigner distribution of the steady-state with Kerr nonlinearity (at minimum S21 with Pin − 122 dBm). The Kerr effect is eventually balanced by the damping, quantum noise and drive. Since the drive now opposes both damping and Kerr, it is less effective at opposing the damping and driving the state away from the origin (compared to panel a). This brings the distribution closer to the origin, and constitutes the second contribution (effect B ) to a reduced amplitude | â |. The center of mass ( â ) (white dot) is compared to the classical steady-state (white cross). Since the Wigner current of the Kerr grows with the amplitude squared | â | 2 ∝ 2 and the drive and dissipation currents grow with | â | and respectively, the reduction in | â | does not linearly follow the driving strength (see Supplementary Information).
both able to describe our measurements is therefore not coincidental: nonlinear damping in the equation for â is a direct consequence of quantum noise.
We now provide an intuitive explanation as to why there is a decrease in the amplitude of oscillations â -leading to an increase in | Ŝ 21 | = |1 − κ ext â /(2 )| -when quantum noise is considered. Because of the Kerr nonlinearity, the uncertainty in the cavity photon number operatorâ †â , translates to uncertainty in the resonance frequency of the oscillator ω r − Kâ †â /2 (see Hamiltonian of Eq. (1)). This has two consequences. Effect A: the signal leaking out of the oscillator into the transmission line inherits the frequency fluctuations of the oscillator. Since we are measuring a single frequency component with our VNA, we will measure a signal of smaller amplitude. Effect B : the driving is less effective at exciting the oscillator because the resonance frequency of the oscillator is fluctuating and no driving frequency will lead to resonant driving. The average number of photons in the oscillator will decrease. This interpretation can be more thoroughly explored in phase space, by making use of the Wigner distribution and Wigner current [33][34][35][36]. We introducê such that the amplitude â is given by the center of mass of the distribution through â = √ 2 dxdp(x+ip)W . The Wigner current J, governs the dynamics of the Wigner function W through the continuity equation ∂ t W + ∇ J = 0. It provides an intuitive visualization of the flow of quasiprobability in phase space.
As a pedagogical starting point, we show in Fig. 4(a) the distribution and different contributions to the current for a coherent state of amplitude α. This state corresponds to the steady-state that would be reached in our resonantly driven system without Kerr nonlinearity. The damping tends to bring each point of the distribution back to the origin. The drive however, is sensitive to phase and acts in a single direction. These two currents are balanced by the quantum noise, which creates a diffusion of the quasi-probability.
In Fig. 4(b), we look at how the Kerr effect deforms the same coherent state, with the damping, driving, and noise temporarily inactive. We see the consequence of the amplitude-dependent resonance frequency of a Kerr oscillator. In phase space, the resonance frequency sets the rate at which a point rotates around the origin. And the amplitude is given by the distance to the origin. The resulting deformation of the coherent state does not bring any point in phase space closer to the origin (total energy, or photon number, remains constant). The center of mass, however, will move closer to the origin. To be convinced of the latter, one can imagine the extreme case of the Kerr deforming the coherent state into a ring circling the origin, so that center of mass would be the origin, and â = 0. This mechanism for reducing â corresponds to effect A previously discussed.
In the steady-state of our experiment, simulated in Fig. 4(c), the Kerr evolution is eventually balanced by the other currents. Due to the large spread of the state in phase, the diffusion induced by quantum noise is weaker, and the damping current further mis-aligned with the drive compared to Fig. 4(a). Since the drive is not parallel to the combined currents of damping and noise, it is less effective at countering them, so less effective at driving the system. Or in other words, in addition to countering the damping, the drive also has to counter the Kerr evolution. As a consequence, the average photon-number tends to decrease , which is the second contribution to a lower amplitude (effect B). In the Supplementary information, we elaborate on why this decrease in amplitude is nonlinear with driving power.
Using a perturbative approach to the quantum model (see Supplementary Information) , we are able to weigh the influence of the two effects A and B in reducing the amplitude â . We rely on an analytical comparison of the amplitude | â | and photon number â †â . For a coherent state, we always have | â | = â †â |. With respect to a coherent state of amplitude α, the reduction in â †â corresponds to half the reduction in | â | in our system (without thermal noise). This means that a loss of photon number is responsible for only half of the observed effect, indicating that half of the increase in |S 21 | can be attributed to effect A, and half to effect B. The same conclusions can be drawn with thermal noise (assuming n th α 2 ).
Finally, we note that for all driving strengths featured in our measurements, quadrature squeezing occurs along an axis u = cos(θ)x + sin(θ)p, rotated by an angle θ with respect the the x−axis. At the highest driving power (Figs. 3,4), the most highly squeezed quadrature is characterized by θ −π/10, where the uncertainty ∆u is 86 % of ∆x for a coherent state.
In conclusion, we have shown how the combination of Kerr nonlinearity and noise, and in particular quantum noise, can lead to observations of nonlinear damping. Crucially, our findings are not limited to the case of su-perconducting resonators. Indeed, preliminary calculations based on our analytical model indicate that this effect has the correct order of magnitude to play a role in the nonlinear damping observed in NEMS systems [4] -however driven by thermal rather than quantum noise. We are therefore confident that this phenomenon can play a valuable role in identifying the origin of nonlinear damping in a broader class of systems, such as NEMS or other Josephson circuits, which will be critical to their use in emerging technologies ranging from carbon nanotube sensors to superconducting quantum computing.
Data and code availability The experimental data and software used to generate the figures in the main text and the supplementary information is available in Zenodo with the DOI identifier 10.5281/zenodo.4565179.

S1. DEVICE FABRICATION AND SETUP
The device shown in Fig. 1 is fabricated in two steps [37]. First, we fabricate the input/output waveguide structures, meandering inductor and capacitor. On a chip of high-resistivity silicon, cleaned in solutions of RCA-1, Piranha, and buffered hydrofluoric acid (BHF), we sputter 60 nm of molybdenum-rhenium (MoRe). A three layer mask (S1813/W(tungsten)/PMMA-950) is then patterned using electron-beam lithography, and is used in etching the MoRe by SF6/He plasma. The mask is finally stripped using PRS 3000.
Secondly, we fabricate the Josephson junctions using the Dolan bridge technique [38]. We first pattern a Methylmethacrylate (MMA) /Polymethyl-methacrylate (PMMA) resist stack with e-beam lithography. After development of the resist, and to ensure a good contact between the aluminum of the junctions and the MoRe, we clean the sample with an oxygen plasma and BHF. Evaporation of two aluminum layers (30 nm and then 50 nm thick) under two angles (± 11 degrees), interposed by an oxidisation of the first aluminum layer, forms the junctions. Removal of the resist mask in N-Methyl-2-pyrrolidone (NMP) at 80 degrees Celsius completes the sample fabrication.
The edges of the chip are wire-bonded to a printed circuit board (PCB) and the chip is placed in a copper box thermally anchored to the 20 milliKelvin stage of a dilution refrigerator. The input to the device is wired through the PCB to the output of a room-temperature vector network analyzer (VNA) and the outputted signal is attenuated at each plate of the dilution refrigerator. The output of the device is wired to the input of the VNA after being amplified with a high-electron-mobility transistor (HEMT) amplifier at ∼4 Kelvin, and a room-temperature amplifier.

S2. THEORETICAL DESCRIPTION
A Hamiltonian describing the series assembly of a inductor L, capacitor C and junction with Josephson inductance L J can be entirely determined with only the knowledge of the admittance Y (ω) = 1/Z(ω) across the Josephson junction, if we replace the latter by a linear inductor L J . This approach assumes weak anharmonicity and damping of the circuit, and is often referred to as black-box quantization [29,39]. The Hamiltonian writeŝ where ω r satisfies Y (ω r ) = 0 and the anharmonicity is given by (S3) The weak anharmonicity assumption which leads to this Hamiltonian writes K ω r .
This circuit loses energy through resistive losses at a rate κ int , and can exchange energy with a transmission line at a rate κ ext . The total rate at which the circuit loses energy is then κ = κ int + κ ext . On one end of the transmission line, we feed a coherent signal with power P in oscillating at ω d . Following quantum input-output theory [25], the dynamics ofâ(t), in a frame rotating at ω d , is given by where ∆ = ω r − K − ω d , and the strength of the drive is characterized by = κ ext P in /2 ω r . Note the factor 2 in the denominator which corresponds to the fact that there are two directions of propagation in the feedline and that only one is occupied by the driving signal. The termŝ corresponds to both thermal noise and quantum vacuum noise. We assume it to be well described by quantum white noise, a stationary random process which is characterized by its 0 mean ŝ = 0 and the correlation functions Here n th corresponds to the average number of excitations induced by the thermal environment at a temperature T with k B corresponding to Boltzmann's constant, to which is added "+1" corresponding to the quantum vacuum fluctuations. Thermalized at 20 mK, and with a frequency ω r 2π × 5.172 GHz, the experimental device has a thermal occupation n th = 4 × 10 −6 . We may thus safely make the approximation n th 0 when describing the experiment. The S 21 parameter is obtained from â(t) as the factor 2 again reflecting that only half of the signal emitted by the circuit will travel towards the receiver of the VNA. As an alternative to the Langevin equation, one may also formulate the problem as a Lindblad master equation [40] ∂ρ ∂t governing the density matrix of the systemρ, where

S3. DATA PROCESSING AND FITTING
Even at lowest driving power, the response of the device does not perfectly fit to a Lorentzian curve, indicating the presence of additional resonances in the measurement chain which could not be calibrated out experimentally. To eliminate these, as well as the change in phase length of the cabling with frequency, we subtract (divide) an affine function of frequency to the measured phase (amplitude). The transformation between measured response S 21,meas and fitted response S 21,fit is thus given by where A, B, C and D are determined through a fit of a low-power response, where the nonlinearity does not come into play, and the response of the device alone S 21,fit (ω) is assumed to be We then reduce the amount of noise as well as the superfluous number of frequency points in the data-set by replacing blocks of 10 successive frequency data-points by their average. The treated data is notably used in the construction of Figs. 2,3. The reduction in number of data-points also facilitates the fitting. We were able to numerically compute S 21 over the 500 frequency points of the data-set in a minimization routine. For each driving power of the dataset, the Python library QuTiP [41,42] was used to solve the Lindblad equation of Eq. (S8). For each power, the absolute difference between the 500 (complex) numerical and experimental points constitute a first contribution to the minimized cost function. The difference in the minimum amplitude of S 21 , and the frequency at which S 21 is minimized, are also added to the cost-function each with a weight of 200 points. The function is minimized using a modified Powell algorithm [43,44], with five free-parameters: ω r , κ int , κ ext , K, and the attenuation that the signal outputted at the VNA will experiences before reaching the device. After convergence of the minimization routine, we performed a slight manual adjustment to the parameters which visually improve the fit. The attenuation is found to be 118.3 dB, which is consistent with the physical attenuation installed at room-temperature and at the different stages of the dilution refrigerator. The device parameters converge to ω r = 2π × 5.172 GHz, κ int = 2π × 189 kHz , κ ext = 2π × 2.12 MHz, K = 2π × 76 kHz. Excellent agreement between theory and experiment is found for the complex transmission at all powers (see Fig. S1,3(a)), the minimum values of |S 21 | as a function of power, and the driving frequency at which this minimum is reached (Fig. S2,3(b)).
The fitted parameters are confirmed by the measurement of a reference oscillator [37], built using the same geometry as the device in Fig. 1 (and in the same fabrication run) but where the SQUID is replaced by a short circuit. We simulate the resonance frequency of the reference oscillator using the finite-element software Sonnet, and compare it to the experimentally measured value. The discrepancy between these allows us to determine the kinetic inductance of the 60 nm sheet of MoRe. Using a sheet inductance of 1.575 pH/sq in Sonnet, the simulated resonance frequency matches the measured value. We then add a lumped element inductor L J at the location of the SQUID in the simulation, and vary it to determine the value of the lumped element inductor L and capacitor C, by fitting simulated resonance frequencies to 1/ (L + L J )C. The Josephson inductance L J is found when this simulated resonance matches the low- (S12) In the absence of (thermal and quantum) noiseŝ = 0, using the notation â = a, one can simply write â †ââ = |a| 2 a. This yields the classical steady-state equation (S13) Note that ∆ = ω r − ω d in the classical case, as the extra −K in the definition of ∆ in the quantum Langevin equation (Eq. (S4)) is a consequence of the quantum fluctuations similar to the Lamb shift [45]. One can see this by rewriting the Kerr termâ †â †ââ using the commutation relations, which would yield a different expression for ∆ in the quantum equation, but would not yield a different Kerr term in the classical steady-state equation.
To simulate nonlinear damping, we consider that the internal damping rate can depend on power κ int → κ int + γ|a| 2 yielding We solve this equation by computing the magnitude |a| and phase ϕ of the amplitude a = |a|e iϕ separately. An equation for |a| is obtained by multiplying Eq. (S13) by its conjugate, yielding We solve this equation by computing the eigenvalues of the polynomial's companion matrix [46] using Python [43]. The phase is then given by ϕ = arctan 2(∆ − K|a| 2 ) κ + γ|a| 2 (S16) Beyond a threshold driving strength , Eq. (S15) exhibits more than one real solution (bifurcating regime). All the data and analysis presented in this work concerns driving powers lower than this bifurcating point. We define the bifurcation power as the first power for which a drive frequency used in the experimental dataset yields more than one real solution to Eq. (S15) with γ = 0.
We may search for the largest attainable amplitude for a given driving strength in Eq. (S15) by assuming a small deviation from the largest attainable amplitude without nonlinearities |a| = α(1 + δ) , |a| n α n (1 + nδ) (S17) with δ 1. The amplitude α is the solution to − κ 2 a + = 0 , (S18) the resonantly driven system, without nonlinearities. Assuming a small deviation from α means we are considering the nonlinearities K, γ as well as the detuning ∆ to be perturbations, such that Since ∆ Kα 2 , we note that these approximations are only valid at the lower powers of our experimental data. Injecting the perturbed expression for |a| into Eq. (S15), we obtain an equation for δ γ 2 /4 + K 2 α 6 (1 + 6δ) + (γκ/2 − 2∆K) α 4 (1 + 4δ) + κ 2 /4 + ∆ 2 α 2 (1 + 2δ) = 2 . (S20) Expanding the solution to second order in K, γ, ∆ through the approximation of Eq. (S19) yields which is minimized for with a maximum (to leading order in K, γ) In order to demonstrate that a classical nonlinear damping model also provides an adequate description of the data, we fit the data to a solution to Eq. (S14). The data processing, and the construction of the cost function, is identical to that described in Sec. S3. The free-parameters, however, are different: we fix the internal and external damping, the oscillator frequency to ω r − K and the attenuation to the values fitted with the quantum model. The remaining free parameters of the fit are the nonlinear damping rate γ, and the Kerr nonlinearity K. After convergence of the same modified Powell algorithm [43,44], and slight manual adjustment, we obtain γ = 2π × 4.58 kHz and K = 2π × 78.3 kHz (2 kHz higher than the value obtained by fitting the quantum model). In Figs. S2,S1(e,f), we show that this model provides a reasonable fit to the data.

S5. QUANTUM AND NOISY SOLUTION (PERTURBATION)
Here we derive a perturbative solution to the Langevin equation of Eq. (S4) in the presence of quantum and/or thermal noise. The aim is to solve the equation for â in the steady-state. Contrary to the noise-less case studied in the previous section, â †ââ = |a| 2 a.

A. Assumptions and small quantities
We assume that the thermal fluctuations, as well as the nonlinearity, slightly move the quantum state away from a coherent state of amplitude â . This approximation comes in two forms. First we will assume the thermal fluctuations entering through the noise termŝ to be small through the condition Secondly, we will assume the anharmonicity to be small, through the condition In this regime, it is practical to express the equations using the deviationd from the steady-state expectation values of the operators involved in the dynamics. Forâ, this means we will writê a = â +d .
We will taked to be a small quantity through the assumption and we will neglect d †ndm with respect to â † n â m when n + m > 2. Note that for K = 0, the steady-state would be a displaced thermal state, andd would correspond to thermal and quantum noise with (d †d ) m scaling with n m th . Assuming that the thermal fluctuations are small n th | â | 2 is therefore a necessary condition to assumed to be small.

B. Preliminary results
We first derive expressions for useful expectation values. Injecting Eq. (S28) in â †â for example yields where we have used the definition ofd to obtain d = d † = 0. We proceed similarly to obtain and by invoking the approximation of Eq. (S29), we also have C. Reformulation of the problem By injecting Eq. (S34) in Eq. (S25), we rewrite the equation for â as A steady state solution for â thus requires knowledge of d †d and d2 . These expectation values will be determined by deriving the equation of motion ofd †d andd 2

D. Equation of motion for d †d
From Eq. (S32), we have d †d = â †â − â * â such that We should thus search for the equation of motion for â †â . Utilizing the Langevin equation of Eq. (S4), we first obtain an equation forâ †â When taking the expectation value of this equation, we follow Ref. [47] to treat the noise terms. Forŝ given by quantum noise, with ŝ = 0 and Eqs. (S5), then givenÂ an arbitrary system operator, we have such that â †ŝ +ŝ †â = − √ κn th . Using Eq. (S32) again to rewrite â †â we have Using the equation of motion for â of Eq. (S36), we finally obtain (S44)

E. Equation of motion for d 2
From Eq. (S33) we have d2 = â 2 − â 2 , such that the equation of motion for d2 writes We should thus search for the equation of motion for â 2 . We start by writing the equation of motion forâ 2 as If quantum noise is taken into account, then we have {â † ,â} = 1 + 2â †â rather than {a * , a} = 2a * a without. To keep track of the quantum noise, we have colored this term in blue throughout the rest of the derivation. When taking the expectation values, we treat the noise terms inŝ following Eq. (S40), and an additional result from Ref. [47] ŝ(t)Â(t) = (n th + 1) resulting in ŝâ +âŝ = 0. By rewriting â 2 using Eq. (S33) and â †â3 using Eq. (S35), we get Using this equation as well as the equation of motion for â of Eq. (S36), we finally obtain

F. Steady-state solution
We obtain a steady state solution for the equations of motion (S36,S44,S50) by equating all time derivatives to zero. We first have where Ω = ∆ − 2K[| â | 2 + 1/4]. Re-injected into Eq. (S44), we get We now invoke the approximations of small anharmonicity K κ Eq. (S27) and small thermal population n th | â | 2 Eq. (S26), and assume the detuning to be close to satisfying the resonance condition derived classically in Eq. (S23) to write Ω −K| â | 2 . This yields and (S55) damped, following where ∆ = (∆ − 2Kn th ) and the nonlinear damping, induced by both thermal and quantum noise, is given by At the detuning ∆ which cancels the imaginary part of this equation, Cardano's formula provides a solution for | â | given by in good agreement with the data (see Fig. S3) By making a final approximation K| â | 2 κ (which is only valid at lower driving powers), and neglecting K 2 | â | 4 /κ 2 , we obtain an equation of a form identical to the classical nonlinear damping equation Eq. (S14) Following the analysis of this equation carried out in Sec. S4, and notably equations (S23,S24), the absolute value of the amplitude is maximized for There are two corrections to the un-perturbed amplitude α: the terms in blue are due to the quantum commutation relations, and the terms in red are due to thermal effects. We note that other authors have already derived more precise analytical equations for the amplitude (in absence of thermal noise) [31,48,49].

G. Discussion
The fact that an inclusion of noise leads to the same equations as nonlinear damping model is confirmed by our fit of the latter model to data. As covered in Sec. S4, the nonlinear damping is found to be 2π × 4.58 kHz, which is very close to the value expected from this calculation. Indeed, in absence of thermal noise n th = 0, Eq. (S57) gives γ = 2K 2 /κ = 2π × 5 kHz.
These calculations also provide a value for the average photon number (relative to the thermal photon number). We derive this quantity close to resonance, where Eq. (S53) holds. Given the expression for â †â of Eq. (S32), â †â = | â | 2 +d †d , and using our expressions for | â | (Eq. (S61)) andd †d (S54), we obtain to leading order in  S4. Wigner current in a resonantly-driven harmonic oscillator We overlay the Wigner function of the steadystate (a coherent state), with the Wigner current associated with the different terms of the Lindblad equation. The sum of these vector fields J plays the role of a current in a continuity equation ∂W/∂t = ∇ J equivalent to the Lindblad equation. a, On resonance ωr = ω d , there is no harmonic evolution in the rotating frame. b, The driving moves the Wigner function towards increasing x. c, The damping moves the Wigner function towards the origin. d, The diffusion term (quantum noise) expands the Wigner function in all directions. e, Since the drive acts along the x−axis and the damping acts more radially, these two terms do not counteract each other. Rather, they move the Wigner function towards a point in phase space, the classical solution to the problem. This is counteracted by the diffusion (quantum noise). Together, these effects give the coherent state its finite size in phase space. f, Another way of showing this is to add the effect of the damping and diffusion, which yields a current acting in the x-axis rather than radially, which counteracts the driving.
yielding, in the limit given by our assumptions (S27) and (S26) the following result The interplay of noise and Kerr nonlinearity leads to a reduction in â †â equal to half the reduction in | â | (Eq. (S63)). The steady-state of the driven-dissipative Kerr system differs thus from a classical or coherent state, for which â †â = | â |, or a thermal state, for which â †â | â | + n th /2| â | if n th | â |. We can thus attribute half of the reduction of amplitude | â | to a reduction of photon number or energy of the oscillator.

A. Mathematical details
We follow Ref. [33] to write the evolution of the Wigner function W as a phase-space continuity equation where ∇ = ∂ x ∂ p and J = J x J p denotes the Wigner current.
This method has already been used to study the Kerr oscillator [34,35], however in the absence of a driving force.
The slightly different Hamiltonian of the Duffing oscillator, this time with both a driving force and damping has also  been studied, however not in the rotating frame [36]. Here we derive an expression for the Wigner current of the driven-dissipative Kerr oscillator in the rotating frame, characterized by the Hamiltonian The slight difference in the form of nonlinearity with respect to Eq. (S1) (allowed by the commutation relations) makes for a more favorable expression when Wigner-transforming the Hamiltonian to phase space coordinates x, p.
We introduce the phase-space operators asx,p througĥ such that [x,p] = i, yieldingĤ omitting constant contributions. To compute the Wigner current we first need the Wigner transform of the Hamiltonian, also known as the inverse of the Weyl transform, defined by [33] H(x, p) = 2 which yields omitting constant contributions. It is easier to prove thatĤ is the Weyl transform of H(x, p) rather than the fact that H(x, p) is the Wigner transform ofĤ. To do so, we may use the McCoy formula [50] p m x n −→ 1 2 n n r=0 n r   Fig. 4(e). e, Zoom-in showing both environmental and driving contributions. Whilst the drive acts on the x−axis, the environment acts more radially. They only partially compensate each other, and the result of summing these currents is shown in f. Together, they tend to diminish the spread in phase of the Wigner function. g, This is compensated by the anharmonicity which tends to increase the spread in phase. To make this clear, we have divided the sum of the harmonic and anharmonic currents into two contributions. In g-left we have shown the part of the current which preserves the spread in position ∂ x 2 /∂t = ∂ p 2 /∂t = 0. This corresponds to the total Wigner current. In g-right we have shown the part of the current which increases the spread in position, and exactly compensates the combined effect of the drive and environment. The latter is plotted in Fig. 4(c).
where −→ designates a Weyl transformation. A consequence of this formula is that The two first relations prove the correspondence between the harmonic and driving terms. Utilizing all three relations, we can demonstrate the correspondance between the Kerr terms inĤ and H(x, p), up to a constant factor which plays no role in successive manipulations of H(x, p).
The unitary evolution of the Wigner function, equivalent to the evolution of the state vector dictated by Schrödinger's equation, is given by [33] where {{H, W }} is called the Moyal bracket. The arrows above the partial derivatives indicates whether the term on the right or left should be differentiated. For example: and Evolution of the driven-dissipative Kerr oscillator starting in a coherent state. We compute the time-evolution of the density matrix following Eq. (S8) in QuTiP [41,42]. The starting state is taken to be a coherent state of amplitude α (see Fig. 4(a)), the amplitude corresponding to the classical solution to the problem. The simulation is carried out with a power Pin = −122, until the steady-state shown in Fig. 4(c) is reached. The time evolution of various observables are plotted over time in order to further illustrate the amplitude reduction mechanism shown in Fig. 4. We see in b that the Kerr causes | â | to reduce over time (negative time-derivative) whereas the drive and environment counteracts this effect (positive time-derivative) until equilibrium is reached at the steady-state (the total time-derivative is 0). This causes the reduction of | â | shown in a. The dashed line corresponds to the steady-state value for | â | in a and c. In c, we plot the square of the average photon number, and see that a reduction in photon number accounts for half of the total reduction in | â |. From d, we find that the photon number is conserved under the Kerr evolution (derivative is 0), whereas the interaction of the drive and environment is at the origin of the decrease. The phase-space interpretation of this effect is that through an increase in the phase variance, the environmental Wigner current is no long parallel to the driving current (see Fig. 4(c)), creating a net current of the probability towards the origin (i.e. a reduction in photon number). Since the amplitude of the damping current is reduced closer the origin (i.e. is smaller for lower photon numbers), whereas the drive is not, equilibrium is found when the probability gets closer to the origin. In e, we plot the evolution of the variance in phase ∆φ. Its increase is shown to be due to the Kerr effect in f, and countered by the interaction of drive and environment.
Since the Hamiltonian only contains terms x n p m with n + m ≤ 4, we may write in this context The Moyal bracket for the harmonic part of the Hamiltonian writes  Fig. 4(b) illustrating this effect. The Wigner currents are drawn with arbitrary lengths for pedagogical purposes and do not reflect their true values. (a -1) A point in phase space which has acquired a phase offset ∆ϕ will be subject to the Kerr and harmonic evolution scaling with the amplitude squared, or in other words 2 .
(2) When the damping and driving is active, steady-state is reached when the Kerr current is compensated by the tangential component of the drive which scales with sin(∆ϕ). (3) For an equilibrium to occur, and assuming the dephasing is small sin(∆ϕ) ∆ϕ, the dephasing should scale linearly with . In the inset we show a simulation of the phase uncertainty ∆ϕ = ∆φ − 1/n which is acquired with increasing driving. To compute ∆ϕ, we have subtracted the phase uncertainty of a coherent state 1/n = 1/ â †â , which does not contribute to the amplitude damping, to the total uncertainty in phase ∆φ. As expected, ∆ϕ increases linearly with . (4) By projecting the equilibrium point on the axis of the center of mass of the distribution (COM), we find a contribution to the reduction of the COM proportional to 2 , this scaling is consistent with our experimental and theoretical findings, and helps to understand why effect A leads to a damping which is nonlinear in . b When we look at the balance of currents in the radial direction, we find that the phase offset ∆ϕ causes the contribution of the driving current to be reduced by a factor ( − cos(∆ϕ))/ 2 . This is consistent with the reduction in photon number proportional to the amplitude squared, and helps to understand why effect B leads to nonlinear damping.
For the drive term The first order derivatives of the Moyal bracket applied to the nonlinearity write And the higher order derivatives of the Moyal bracket applied to the nonlinearity yield Given a change in density matrix ∂ρ ∂t , we can compute the change in amplitude from where ∂ ∂t â = ∂ ∂t Tr(âρ) = Tr(â ∂ρ ∂t ) .
Consistently with the explanations surrounding Fig.4, we observe that the Kerr effect causes a decrease of the amplitude ( Fig. S7(b)) (effect A). And this effect is eventually compensated by the combined effect of the drive and environment. We also observe that the average number of photons decreases due to the interaction of drive and damping ( Fig. S7(f)), and this effect seems to follow the increase in phase variance (Fig. S7(f)) (effect B).
Whilst the explanations above and surrounding Fig. 4 help in understanding the physical mechanism behind the damping, they do not elucidate its nonlinear nature. To understand the latter, we turn to estimating the equilibrium of Wigner currents (see Fig. S8 for visual point of reference). We consider the state of Fig. 4(b), after the Kerr effect has acted on the coherent state. Assuming we are driving at the new resonance frequency ∆ = Kα 2 , we determine the scaling of currents for a point which has undergone a characteristic amount of dephasing ∆ϕ. We consider a point in phase space which has a distance from the origin given by x + ∆x, where ∆x = 1/ √ 2 is the uncertainty in position of a coherent state and x = /( √ 2κ). We find that the Kerr and harmonic evolution Wigner current J Kerr,1 + J harmonic scales with 2 , and points in a perpendicular direction to the damping. Our numerics show that in practice the second contribution to the Kerr current J Kerr,2 is negligibly small, and we see in Fig. 4(c) that the influence of quantum noise at large dephasing is also negligible. In this discussion we will simplify the steady-state condition as requiring the orthogonal Wigner currents of Kerr/harmonic evolution and damping to be matched by the drive which is oriented horizontally. As illustrated in Fig. S8 this condition leads to an understanding of the nonlinear evolution of the damping, which decreases the total amplitude α by a factor 2 , resulting in the nonlinear scaling of amplitude α 2 ∝ α 3 found both experimentally and theoretically.

S7. QUADRATURE SQUEEZING
We determine the degree of squeezing of the steady-state solutions to Eq. (S8) for different powers. We consider squeezing along the orthonormal coordinate set u, v with u = cos(θ)x + sin(θ)p. The coordinate u corresponds to the operatorû = e iθâ † + e −iθâ for which we compute the uncertainty ∆u for varying θ. As shown in Fig. S9, for each measured power, there is a angle θ for which the uncertainty is below the uncertainty of a coherent state, demonstrating quadrature squeezing.