Apparent nonlinear damping triggered by quantum fluctuations

Nonlinear damping, the change in damping rate with the amplitude of oscillations 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 resonator, we show that from the interplay of quantum fluctuations and the nonlinearity of a Josephson junction emerges a power-dependence in the resonator response which closely resembles nonlinear damping. The phenomenon can be understood and visualized through the flow of quasi-probability in phase space where it reveals itself as dephasing. Crucially, the effect is not restricted to superconducting circuits: we expect that quantum fluctuations or other sources of noise give rise to apparent nonlinear damping in systems with a similar conservative nonlinearity, such as nano-mechanical oscillators or even macroscopic systems.

(Dated: November 28, 2023) Nonlinear damping, the change in damping rate with the amplitude of oscillations 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 resonator, we show that from the interplay of quantum fluctuations and the nonlinearity of a Josephson junction emerges a power-dependence in the resonator response which closely resembles nonlinear damping.The phenomenon can be understood and visualized through the flow of quasi-probability in phase space where it reveals itself as dephasing.Crucially, the effect is not restricted to superconducting circuits: we expect that quantum fluctuations or other sources of noise give rise to apparent nonlinear damping in systems with a similar conservative nonlinearity, such as nanomechanical oscillators or even macroscopic systems.
Amplitude-dependent (nonlinear) damping 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][17] made for example from carbon nanotubes, graphene [4,18] or superconducting metal [19].In some of these systems the nonlinearity is well explained [20][21][22][23].Most notably the saturation of two-level systems in the environment can cause negative nonlinear damping: the damping rate decreases as the power injected into the system increases [11][12][13]19].But the origin of an increase in damping with power in certain NEMS [4,[15][16][17] 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 [24,25] in superconducting quantum computers [8].
We study a dephasing effect in superconducting circuits [26], which phenomenologically appears as nonlinear damping when measuring the resonant response of a resonator.Central to the observed physics is the 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 [27].For this reason, the phenomena discussed here are applicable to all systems featuring a similar nonlinearity in their resonance frequency, for example the carbon-nanotubes mentioned above, or even a macroscopic mechanical pendulum.We focus on the regime where this nonlinearity is small, as in Josephson parametric amplifiers [25], rather than the single-photon nonlinear regime used to construct artificial atoms in circuit quantum electrodynamics [27].Because of their small nonlinearity, such systems are often thought to be completely described by the classical Kerr oscillator [5,28].
Here however, we report on an effect that is not expected from the classical Kerr oscillator.More specifically, we present a phenomenon triggered by the interplay between the quantum noise and Kerr anharmonicity of the oscillator, which closely resembles nonlinear damping in the steady-state response of the oscillator.The apparent nonlinear damping is first experimentally characterized by probing the frequency response of the resonant circuit.Our observations are then accurately described by a quantum theory of a damped driven Kerr oscillator devoid of ad-hoc nonlinear damping, but which takes into account the effect of quantum noise.Moreover, focusing on an oscillator steady-state below its bistability threshold, a Gaussian state approximation [29] allows us to demonstrate that, in a close vicinity of the resonance, the expected amplitude of oscillations is akin to that of a driven classical Kerr oscillator 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 (Fig. 1) is constructed 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 [30].The junction induces an anharmonicity of strength K = 2π × 80 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 [31] Ĥ 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π × 186 kHz.This is typically due to losses in the different dielectric materials traversed by the electric fields [32].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 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 [27].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 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 Sec.S7).
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 average 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. 2a).
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. 2a), 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 as the drive power increases.At the highest drive power for which data is displayed (P in = −124 dBm), the internal damping rises to 2π × 255 kHz.We note that this power lies below the bistability threshold (see Supplementary Information Secs.S3 and S4F).Such nonlinear damping can be included in the model of Eq. ( 2) through We fit a solution of the resulting equation to the data (see Supplementary Information Sec.S2), observing good agreement (Fig. 2) for γ = 2π × 5.02 kHz.
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 of the internal damping rate, whilst we observe the opposite.Here, we show that in our system this nonlinear damping behavior can be explained purely by dephasing triggered by the joint action of the intrinsic quantum noise and the Kerr anharmonicity 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 numerically 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 Sec.S2), revealing excellent agreement to the data (Fig. 3).Note that we recently became aware of an analytical solution to this Lindblad equation [33,34], which may have simplified our approach.
In Fig. 3, we compare this quantum model to the classical model: the solution to Eq. ( 2), which features neither nonlinear damping nor quantum noise.The only difference between the quantum model -which predicts the increase in Min|S 21 | -and that of Eq. ( 2) -which predicts a constant Min|S 21 | -lies in the value of the commutator [â, â † ].In fact, by taking the trace Tr(â • ) of Eq. ( 3), and assuming the amplitude to be a complex Phase space picture of apparent damping.Here phase space operators are defined by â = (x + ip) / √ 2, see Supplementary Information Sec.S5A for further details.a, Wigner distribution of the steady-state in absence of Kerr nonlinearity (driven at ωr with Pin = −124 dBm).The balance between quantum noise, damping and drive is shown by vectors corresponding to Wigner currents.b, Growth of phase uncertainty of a coherent state under Kerr nonlinearity.The amplitudedependent 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 resonant amplitude's magnitude |⟨â⟩|.c, Wigner distribution of the steady-state with Kerr nonlinearity (at minimum |S21| with Pin = −124 dBm).The Kerr effect is eventually balanced by the damping, quantum noise and drive.Since the drive now opposes both damping and Kerr effect, 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 resonant amplitude's magnitude |⟨â⟩|.The center of mass (⟨â⟩) (white dot) is compared to the classical steady-state (white cross).Since the Wigner current of the Kerr effect 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 Sec.S5B).
number â → a such that [a, a * ] = 0, we arrive at Eq. (2).Quantum noise can therefore lead to the entirety of the change in Min|S 21 |.Thermal noise could lead to a similar effect, but is expected to be negligible in our experiment (see Supplementary Information Sec.S7).
Beyond describing the data, this model can lead to a nonlinear damping equation for the expectation value ⟨â⟩.In the Supplementary Information (Sec.S4), we derive an analytical formula that captures the behavior of the steady-state response that is in good agreement with the numerical simulations.We find that in a resonance scenario, whenever the Kerr effect and thermal noise have only a perturbative effect on the system, the corresponding steady-state expectation value for the amplitude ⟨â⟩ matches that of a classical non-linearly damped driven classical Kerr oscillator.That is, the steady-state amplitude is ruled by an equation resembling Eq. ( 2), but where the quantum and thermal noise lead to a nonlinear Here the familiar + 1 2 stems from quantum noise, which has the same effect as half a quantum of thermal noise.The fact that the nonlinear damping model and the quantum model are both able to describe our measurements is therefore not coincidental: whilst there is no microscopic process leading to nonlinear damping (i.e.loss of energy), there is apparent nonlinear damping in the equation for ⟨â⟩ when accounting for the presence of quantum noise.Similar results were derived for the classical [35] and quantum [36] spectrum of undriven oscillators, and also in work studying the spectrum of a probe field in the presence of a strong pump field [37].The difference here is that we are instead interested in the power-dependence of the scattering parameter S 21 .
We now provide an intuitive explanation as to why there is a decrease in the amplitude of oscillations ⟨â⟩ -leading to an increase in the minimum of |⟨ Ŝ21 ⟩| = |1 − κ ext ⟨â⟩/(2ϵ)| -when quantum noise is considered.Because of the Kerr nonlinearity, the uncertainty in the 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 [38][39][40][41].We introduce x = (â + â † )/ √ 2 and p = −i(â − â † )/ √ 2, 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. 4a 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. 4b, 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 damped driven quantum 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 effect 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. 4c, the evolution of the Kerr effect 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 misaligned with the drive compared to Fig. 4a.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 evolution of the Kerr effect.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 (Sec.S5B), we elaborate on why this decrease in amplitude is nonlinear with driving power.
Using a Gaussian state approximation (see Supplementary Information Sec.S4) , we are able to weigh the influence of effect A and effect B in reducing the value of the resonant amplitude ⟨â⟩.We rely on an analytical comparison of the corresponding amplitude's magnitude |⟨â⟩| and photon number ⟨â † â⟩, and the fact that effect A does not affect the photon number, whereas effect B reduces the amplitude by reducing the photon number from its expected value for a coherent state ⟨â † â⟩ = |⟨â⟩|.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 reduction in photon number is responsible for only half of the observed effect, indicating that half of the increase in Min|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 ).
Whilst we have focused on the case n th ≪ α 2 , where the damping seems to increase with the amplitude of oscillations, the opposite regime n th ≫ α 2 has already been explored experimentally [42,43] and bears some common features with this work.When thermal fluctuations dominate, the state of the oscillator is well described as a statistical mixture of oscillatory amplitudes, each shifting the resonance frequency by a different amount given by the Duffing nonlinearity.This results in a broadening of the resonance line-shape when the oscillator is probed, which has been phenomenologically interpreted as an increase in damping, for example in carbon nanotubes [42].This picture even extends to the case n th ≪ α 2 where a residual broadening persists due to quantum heating of the oscillator by the driving field [44].
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 θ ≃ −0.11π,where the uncertainty ∆u is 83 % of ∆x for a coherent state.
In conclusion, we have shown how the combination of Kerr nonlinearity and noise, and in particular quantum noise, leads to a dephasing that can manifest in the same way as nonlinear damping.Crucially, our findings are not limited to the case of superconducting 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 nature of nonlinear damping effects 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 are available in Zenodo with the DOI identifier 10.5281/zenodo.4565179.

S1. 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 [31,45].The Hamiltonian writes where ω r satisfies Y (ω r ) = 0 and the Kerr constant is given by 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 [26], 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.As described in Sec.S7, thermal noise is heavily attenuated such that the resonator has a thermal occupation n th < 0.05, a much smaller source of noise than quantum fluctuations.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 in terms ofa Lindblad master equation [46] governing the density matrix of the system ρ, where

S2. 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).Transmission coefficient of the superconducting resonant circuit.Measurements of the transmission coefficient S21, processed following Sec.S2, are used in a curve fitting routine that seeks, for each driving power in the dataset and n th = 0, a numerical solution of the steady-state equation (3) of the damped driven quantum Kerr oscillator with the set of oscillator parameters that fits best to the measured data in accordance to a minimization of the error in: a, the magnitude |S21| and b, phase Arg(S21) of S21 as a function of driving frequency, as well as in c, the minimum Min|S21| of |S21| and d, the driving frequency ωmin at which this minimum occurs as a function of drive power.Data-points are rendered with blue dots and curve fits with solid lines.Each trace in a and b corresponds to a different drive power.At the top trace the drive power is Pin = −135.23dBm, then this increases 1 dBm per trace, reaching Pin = −118.23dBm at the bottom trace.We offset every trace in a and b by -10 dB and -2 rad respectively to ease examination.
The transformation between measured response S 21,meas and fitted response S 21,fit is thus given by (A + Bω)e C+Dω S 21,fit (ω) = S 21,meas , (S10) 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 data-set, the Python library QuTiP [47,48] was used to solve the Lindblad equation of Eq. (S8) with n th = 0.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 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 [49,50], with five free-parameters: ω r , κ int , κ ext , K, and the attenuation that the signal outputted at the VNA experiences before reaching the device.The attenuation is found to be 118.3dB, 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π × 186 kHz , κ ext = 2π × 2.12 MHz, K = 2π × 80 kHz.Excellent agreement between this numerical computation based on the Lindblad equation and the experimental data is found for the complex transmission at all drive powers, as well as for the drive power dependence of both, the minimum value of |S 21 | and the driving frequency at which this minimum is reached, see Fig. S1.
The fitted parameters are confirmed by the measurement of a reference oscillator [51], 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-power resonance frequency measured in Fig. 2. The circuit parameters obtained from this analysis are L = 2.93 nH, C = 288 fF and L J = 0.35 nH.The resonance frequency and Kerr constant calculated from these parameters matches those determined by fitting the data with small deviations of 0.1% and 2% respectively.
Finally, we would like to acknowledge that further measurements could have been performed to validate both the system parameters and the physical mechanism of apparent non-linear damping.The Kerr constant could be verified by strongly driving the system off-resonance and measuring the resulting Stark shift with a low-power probe.The effective non-linear damping parameter γ could be verified through the change in S 11 expected from a varied thermal occupation n th .Here thermal occupation could be varied by injecting a known white noise or by increasing the base temperature of the dilution refrigerator.Lastly, the SQUID tunability could be used to demonstrate an understanding of the apparent non-linear damping with a different Kerr constant and resonator frequency.Unfortunately, such measurements can no longer be carried out in a practical timescale: the device was fabricated for a different purpose and measured long before the analysis presented this work was performed and is no longer in working condition.

S3. CLASSICAL, NOISE-LESS SOLUTION AND BIFURCATION POINT
Taking the expectation value of Eq. (S4), in the steady-state (d⟨â⟩/dt = 0), yields 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 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 [52].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 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 [53] using Python [49].
The phase is then given by Beyond a critical drive power there are drive frequencies for which Eq. (S15) possesses not only one but three real solutions.In that case, the steady-state response of the oscillator will exhibit bistability as the drive frequency or detuning is varied.The point (detuning) at which the onset of bistability takes place is known as a critical point (detuning).The same applies to Eq. (S15) with γ = 0, in which case the ensuing critical power P c = 2 √ 3ℏω r κ 3 /(9|K|κ ext ) ≃ −122.2 dBm is obtained through a simple stability analysis, see e.g.Appendix C of Ref. [54].In what follows, we shall then use P c as an estimate of the actual critical power that follows with γ ̸ = 0.An examination of Figs.S1(a,b), suggests that this is a reasonable approximation.
Most of our measurements employ drive powers lower than P c .Furthermore, we observe that the phenomenon we are concerned with, namely, the reduction in the oscillator's amplitude as the drive power increases, manifests clearly in all those measurements.Either a fully classical description of a steady-state in a nonlinear quantum system (the superconducting resonant circuit in our case) or that based on a linearization of quantum fluctuations around such a steady-state, are accurate far from a critical point.However, both descriptions are known to fail in the vicinity of a critical point, see Ref. [29] for more details on this regard.As a consequence, we focus all the theoretical analysis we present in this work on a frequency response of the resonant circuit describable by a steady-state oscillator well below the bistability threshold.To ensure this condition, we only consider drive powers P in ≤ −124 dBm, further motivated by Sec.S4.
Bearing in mind the above, we now 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 with δ ≪ 1.The amplitude α is the solution to 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 δ 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 γ) 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.S2.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 [49,50], we obtain γ = 2π × 5.02 kHz and K = 2π × 78.3 kHz (1.7 kHz lower than the value obtained by fitting the quantum model).In Figs.S2,S3, we show that this model provides a reasonable fit to the data we address in our theoretical description (below the bistability threshold of the steady-state of our model oscillator).

S4. QUANTUM AND NOISY SOLUTION IN THE GAUSSIAN REGIME
Here we derive an approximate 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 the ensemble average ⟨â⟩ in the steady-state and in a small neighborhood of the resonance.In order to lighten our notation, we drop the time dependency in the equations and, unless it is explicitly stated otherwise, identify the ensuing dynamical variables with those that define the steady-state of the system.Contrary to the noise-less case studied in the previous section, here ⟨â † ââ⟩ ̸ = |a| 2 a.This nonlinear term in Eq. (S25) above leads to a Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) like hierarchy of coupled differential equations, in which, e.g, the dynamics of the moment ⟨â n ⟩ of order n requires knowledge of the higher order moment â † a n+1 , with n ∈ N + [55].To perform a tractable analytical study of such infinite hierarchy of equations we shall truncate it appropriately.We do so using a Gaussian state approximation which allows for closing the hierarchy, thus guaranteeing a physically meaningful steady-state solution [29,55,56].The resulting closed system of coupled differential equations involves only first order statistical moments (i.e.ensemble averages) and second order statistical moments (covariances) of the system's canonical operators.This approach, different to assuming a coherent state ansatz for the oscillator's quantum state, does not prescribe an operator's ensemble average to be entirely ruled by the classical dynamics of the system.Instead, in this case an operator's ensemble average is treated as an unknown variable to be determined self-consistently with its covariances, thus including some information of the quantum dynamics of the system into that of the operator's ensemble averages.Next, we present the conditions under which we may apply such an approach, followed by the description of a method for obtaining the ensemble average ⟨â⟩ resulting from our Gaussian state ansatz.

A. Assumptions and small quantities
We assume that the thermal fluctuations, as well as the nonlinearity, have a small effect on the system.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 deviation d from the steady-state expectation values of the operators involved in the dynamics.For â, this means we will write We will take d to be a small quantity through the assumption such that we will neglect ⟨ d †n dm ⟩ with respect to ⟨â †n âm ⟩ − ⟨ d †n dm ⟩ and thus use ⟨ d †n dm ⟩ ≃ 0 whenever n + m > 2. This second order ladder approximation amounts to considering that the steady-state is well described by a Gaussian state.We note that for a Gaussian state ⟨( d † d) m ⟩ scales with n m th .Therefore, considering that the thermal fluctuations are small n th ≪ |⟨â⟩| 2 is necessary for Eq.(S29) to hold.
The assumption (S29) enables a linearization of the steady-state equations of the ensemble averages and covariances of the oscillator's canonical operators that fully characterize the resulting Gaussian steady-state.As shown for a similar paradigm in Ref. [56], given a limit of weak anharmonicity, c.f. Eq. (S27), and far from a multistability threshold, this linearization leads to an analytical steady-state solution in good agreement with numerical calculations.This approach is then applicable for our system below the bistability threshold.Although lacking the precision of other methods [33,34], it provides a faithful description of our experimental observations, with the additional benefit that quantities of interest (such as averages and covariances of canonical operators) can be expressed more simply.

B. Preliminary results
We first derive expressions for useful expectation values.Injecting Eq. (S28) in ⟨â † â⟩ for example yields where we have used the definition of d 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 of d † d and d2 .

D. Equation of motion for ⟨ d † d⟩
From Eq. (S32), we have 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. [57] 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 E. Equation of motion for d2 ⟩ 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. [57] ⟨ŝ(t) 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 solve for the covariances to get where Ω = ∆ − K(2|⟨â⟩| 2 + 1/2).By injecting our expressions for ⟨ d † d⟩, ⟨ d2 ⟩ (c.f.Eqs.(S51) and (S52), respectively) in Eq. (S36) and after some algebra, we find that, in the steady-state, the amplitude is ruled by the following equation Equation (S53) shows that the interplay of the oscillator's nonlinearity and the noise gives rise, among other effects, to a nonlinear broadening of the oscillator's steady-state response, the strength of which is given by These effects resulting from the interplay of the noise and the oscillator's nonlinearity manifest most prominently as we approach to a resonance scenario, for if, on the contrary |∆/κ| → ∞, the distribution 1/(4Ω 2 /κ 2 − 4K 2 |⟨â⟩| 4 /κ 2 + 1) → 0.Moreover, for a sufficiently weak anharmonicity, in a narrow frequency window enclosing the resonance and on resonance itself, we shall see next that 1/(4Ω 2 /κ 2 − 4K 2 |⟨â⟩| 4 /κ 2 + 1) ≃ 1.In that case, a nonlinear broadening of the resonance is the dominant effect that stems from the presence of noise in the system.Indeed, in order to be able to approximate the value of the aforementioned distribution by 1 it is necessary that |4Ω 2 /κ 2 − 4K 2 |⟨â⟩| 4 /κ 2 | ≪ 1.The solutions to this inequality depend on the value of K|⟨â⟩| 2 /κ.That is, the solution to the inequality is given by the sets (S55) and by the set Using the oscillator's parameters measured in our experiment and considering either n th = 0, n th = 0.05 (the upper bound for the average number of thermal photons in our experiment as estimated in Sec.S7 B) or n th = 1/2, a numerical solution of the master equation (S8) reveals that our anharmonicity is weak enough so as to guarantee K|⟨â⟩| 2 /κ ≲ 0.48 for every detuning, −4MHz ≲ ∆/(2π) ≲ 2.92MHz, and all the drive powers P in ≤ −124 dBm addressed in our analysis.Drive powers above P in = −124 dBm but still close to the critical power (at about P c ≃ −122.2 dBm) that sets the onset of a bistable steady-state, not only break down the condition K|⟨â⟩| 2 /κ < 1/2, but lead to a solution that starts satisfying more narrowly the necessary assumption (S29).Let us then focus on this latter regime for which 0 < K|⟨â⟩| 2 /κ < 1/2 is fulfilled.For a given driving strength resonance is attained with the detuning ∆ ⋆ for which the magnitude |⟨â⟩| of the steady-state amplitude reaches its maximum |⟨â ⋆ ⟩|.FIG.S4.Conditions leading to a standard nonlinear broadening, proportional to |⟨â⟩⋆| 2 , of the resonance of a damped driven quantum Kerr oscillator close to a Gaussian state.In a regime complying with 0 < K|⟨â⟩| 2 /κ < 1/2 for every applied input power, (a) the on resonance detuning, ∆⋆, lies in between the detuning bounds ±∆ b that derive from the frequency bounds ±Ω b delimiting the on resonance equivalent of the set defined in Eq. (S56) and, at the same time (b) both Ω⋆ = ∆⋆ − K(2|⟨â⟩⋆| 2 + 1/2) and the magnitude of the oscillator's amplitude maximum |⟨â⟩⋆| are such that the quantity All the plotted quantities are obtained from a numerical solution of the Lindblad eqation (S8) with parameter values as extracted from the experimental measurements and shows the dependence on P in of this resonance detuning ∆ ⋆ and the detuning bounds to the on resonance version of the set defined in Eq. (S56).Note that knowledge of a maximum |⟨â⟩ ⋆ | as computed through a numerical solution of the master equation (S8) enables us to determine both, its corresponding detuning ∆ ⋆ as well as the ensuing detuning bounds ±∆ b .That is how we obtain each detuning value we plot in Fig. S4(a), using in the numerical computation n th = 1/2 and the same parameter values we extract from the experimental data.Clearly, we observe that −∆ b < ∆ ⋆ < ∆ b , i.e., that (−∆ b , ∆ b ) defines an interval of detunings comprising the resonance.In Fig. S4(b), we check the validity of the assumption −∆ b ≪ ∆ ⋆ ≪ ∆ b or, more explicitly, the relative error of the approximation 1/(4Ω The curve shows that the deviation of the distribution from the target value 1 remains lower than a 2% for the entire range of input powers used in our theoretical analysis.The same curve as computed using either n th = 0 or n th = 0.05 shows an error that stays below a 5% for input powers P in ≤ −126 dBm, after which the error increases up to a 10% for P in = −124 dBm. As anticipated above, we may then conclude that Eq. ( S56) is rather well fulfilled on resonance.Eq. (S56) can be equally satisfied for an arbitrarily small range of detunings ∆ ⋆ − ϑ ≤ ∆ ≤ ∆ ⋆ + ϑ around the resonance, with ϑ a real constant such that 0 < ϑ ≪ |∆ ⋆ |.Thus, in order to evaluate the oscillator's steady-state amplitude near resonance, we constrain our analysis to such range of detunings and set 1/(4Ω 2 /κ 2 − 4K 2 |⟨â⟩| 4 /κ 2 + 1) ≃ 1.This allows us to approximate the equation for the steady-state amplitude as We may now associate the resonance scenario with the approximated value of the maximum of the steady-state amplitude's magnitude that results from this Eq.(S57).We achieve this with the detuning in good agreement with the data (see Fig. S5).In panel a we use the same regime of parameters as the measured device, notably the same Kerr parameter K, with both no thermal population n th = 0 and half a quantum of thermal population n th = 1/2.In panel b, we show that the accuracy of the analytical model increases as we decrease K by a factor 10. For the calculations with thermal population, good agreement is found in the regime n th ≪ |⟨â⟩| 2 covered by our assumptions.
Similarly to Sec.S3, to ease the subsequent analysis and acquire a clearer insight of the behaviour of the oscillator's resonant amplitude, we derive next a simpler approximated solution to Eq. (S57) based on perturbation theory.Given that K|⟨â⟩| 2 /κ < 1/2, if we further take into account our original assumptions of weak anharmonicity and strength of thermal fluctuations, c. f. (S27) and (S26), respectively, we find that ε ≪ 1.The solution of the oscillator's steady-state amplitude that we seek comes from a series expansion ⟨â⟩ = l=0 ε l ⟨â⟩ l in powers of the small parameter ε = γ/κ.For our simple description of the resonance, we content ourselves with the regime for which the first order correction in ε is dominant.This leaves us with the expansions ⟨â⟩ ≃ ⟨â⟩ 0 + ε⟨â⟩ 1 and . Likewise, we expand the detuning in the form ∆ = ∆ 0 + ε∆ 1 .By replacing these expansions into Eq.(S57), neglecting terms of higher than ε 1 and grouping terms of the same order in ε, we obtain the following set of equations In the resonance scenario described by the set of equations above, the detuning must be such that it cancels each of the brackets multiplied by the imaginary unit in the left hand side of Eqs.(S59) and (S60).Solving the resulting equations for the corresponding resonant amplitude yields ⟨â⟩ ⋆ ≃ ⟨â⟩ ⋆,0 + ε⟨â⟩ ⋆,1 with ⟨â⟩ ⋆,0 = α and ⟨â⟩ ⋆,1 ≃ −|⟨â⟩ 0 | 2 ⟨â⟩ 0 = −α 3 .The simple approximation for the maximum of the absolute value of the amplitude is maximized for where ∆ ⋆,0 = −K + Kα 2 + 2K(n th + 1/2), and reads There are two corrections to the unperturbed 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) [33,58,59].

G. Discussion
The fact that an inclusion of noise leads to the same equations as the nonlinear damping model is confirmed by our fit of the latter model to data.As covered in Sec.S3, the nonlinear damping is found to be 2π × 5.02 kHz, which is close to the value expected from this calculation.Indeed, in absence of thermal noise n th = 0, Eq. (S54) gives γ = 2K 2 /κ = 2π × 5.6 kHz.
We can also use the approximate description of the resonance scenario embodied in Eqs.(S61) and (S62) above to provide a value for the corresponding average photon number.Given its equation of motion (S42), in the steady-state and on resonance, we have ⟨â † â⟩ ⋆ = n th + α(⟨â⟩ ⋆ + ⟨ â † ⟩ ⋆ )/2.Based on our previous perturbation series expansion and hence 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.(S64)).The steady-state of the driven-dissipative Kerr system differs thus from a coherent state, for which We can thus attribute half of the reduction of amplitude |⟨â⟩ ⋆ | to a reduction of photon number or energy of the oscillator.
The subscript ⋆ is omitted in the rest of the manuscript where the resonant condition is implied through context.

S5. WIGNER CURRENT A. Mathematical details
We follow Ref. [38] to write the evolution of the Wigner function W as a phase-space continuity equation where This method has already been used to study the Kerr oscillator, however in the absence of a driving force [39,40] or in a non-rotating frame [60].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 [41].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 as x, p through 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 [38] H 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 [61] 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 function, equivalent to the evolution of the state vector dictated by Schrödinger's equation, is given by [38] ∂ 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 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 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  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.
which in our regime of parameters has a negligibly small contribution to the total current.The non-unitary evolution of the Wigner function, equivalent to the Lindblad equation of Eq. (S8), is given by [41] These expressions are utilized in the discussion surrounding Fig. 4. We supplement Fig. 4 by the discussion below, where we provide details on the construction of the figure and further arguments in favor of its interpretation.We also provide in Fig. S8 a detailed plot of the Wigner currents for the steady-state shown in Fig. 4c.Finally, we plot in Fig. S6 and Fig. S7 the Wigner currents for the pedagogical cases of a resonantly and off-resonantly driven harmonic oscillators respectively.
B. Supplementary information for Fig. 4 In Fig. 4b, we show a coherent state, of amplitude α, and the state which results from evolving that coherent state under the Kerr effect Ĥ/ℏ = λKâ † â − (K/2)â † â † ââ for a time t = 1/(45K).This unitary evolution is computed using QuTiP [47,48].To ensure that the state is centered around the x axis in phase space, we add a harmonic evolution characterized by the coefficient λ = 15.After evolution of the coherent state due to the Kerr effect, the average photon number remains unchanged.This is theoretically expected since the Hamiltonian dictating this evolution Ĥ = ℏλKâ † â − ℏ(K/2)â † â † ââ commutes with the photon number operator â † â.S7.Wigner current in an off-resonantly driven harmonic oscillator.We overlay the harmonic a, drive b, damping c, and diffusion (quantum noise) d Wigner currents on the Wigner distribution of the steady-state.e, Contrary to the resonant case, the sum of diffusion and damping no longer cancel the drive.Instead, these currents point away from the direction of the harmonic evolution, such that the net total current is 0. f, The remaining current moves points of the Wigner distribution along lines of equal probability such that ∂W/∂t = 0 .
The amplitude |⟨â⟩|, or the distance between the center of mass of the distribution and the origin in phase space, decreases.In Fig. 4b, the decrease in amplitude is of 0.34 percent, which translates to a considerable 0.33 dB difference in |S 21 |.We have called this effect A.
As a result of the deformation of the Wigner distribution, the drive is less affective at countering the environmental effects (noise and damping).If we call ⃗ J env = ⃗ J damping + ⃗ J diffusion , and integrate the absolute value of this current, and the driving current over one of the distribution of Fig. 4b, we find Whilst this is true for the coherent state and the state deformed by the Kerr effect, we obtain a different result if we project the driving current in the same direction as the environmental effects.For a coherent state, we have So the drive is aligned with the environmental effects and exactly counters them.However, for the deformed state | ⃗ J env |dxdp = 6.10 MHz, and Jenv| |dxdp = 6.07 MHz.So after effect of the Kerr, whilst the drive matches the environmental effects in absolute magnitude, in the direction of the environmental current the driving current is smaller.The Wigner distribution thus tends to move towards the origin in phase space, hence reducing the total photon number.We have called this effect B. This is further illustrated if we track the evolution of various observables as a coherent state of amplitude α evolves to the steady-state following Eq.(S8), see Fig. S9.We look at the amplitude |⟨â⟩| = |Tr(âρ)|, the square root of the number of photons ⟨â † â⟩ = Tr(â † â ρ) and the variance of the phase ∆φ = Tr( φ2 ρ) − Tr( φρ) 2 .Additionally, we look at the time derivative of these observables induced by various terms of Eq. (S8).
We distinguish the influence of different terms by computing the change in density matrix induced by drive and damping  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. 4c.
and by the harmonic and Kerr evolution ) .

(S87)
For the photon number 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 [47,48].The starting state is taken to be a coherent state of amplitude α (see Fig. 4a), the amplitude corresponding to the classical solution to the problem.The simulation is carried out with a power Pin = −124, until the steady-state shown in Fig. 4c 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 longer parallel to the driving current (see Fig. 4c), 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.
Origin of the nonlinear nature of the damping.Following a unitary evolution of a coherent state due to the Kerr effect, the state acquires an uncertainty in phase characterized by ∆φ.We show in both panels the Wigner function of Fig. 4b 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 are active, the 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 an apparent 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 an apparent nonlinear damping.
Consistently with the explanations surrounding Fig. 4, we observe that the Kerr effect causes a decrease of the amplitude (Fig. S9b) (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. S9f ), and this effect seems to follow the increase in phase variance (Fig. S9f ) (effect B).
Whilst the explanations above and surrounding Fig. 4 help in understanding the physical mechanism behind the dephasing which manifests in the same way as damping, they do not elucidate its nonlinear nature.To understand the latter, we turn to estimating the equilibrium of Wigner currents (see Fig. S10 for a visual point of reference).We consider the state of Fig. 4b, after the Kerr effect has acted on the coherent state.Assuming we are driving at the new resonance frequency[62] ∆ = 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. 4c 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. S10 this condition leads to an understanding of the nonlinear evolution of the apparent damping, which decreases the total amplitude α by a factor ϵ 2 , resulting in the nonlinear scaling of amplitude αϵ 2 ∝ α 3 found both experimentally and theoretically.

S6. 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. S11, for each measured power, there is a angle θ for which the uncertainty is below the uncertainty of a coherent state, demonstrating quadrature squeezing.On the left and on top, we plot the marginal distributions for u and v respectively (blue), overlayed with the distributions of a coherent state (dashed lines) with the same amplitude ⟨â⟩ as the plotted Wigner function.In this case, the uncertainty ∆u for the steady-state is 83 percent of the uncertainty one would obtain for a coherent state, demonstrating quadrature squeezing.

S7. DEVICE FABRICATION AND SETUP
A. Fabrication The device shown in Fig. 1 is fabricated in two steps [51].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 [63].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 oxidization 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.

B. Experimental setup
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), the signal coming from the VNA is attenuated at each plate of the dilution refrigerator.The cryogenic wiring is detailed in Fig. S12.The noise photon occupation n i of a coaxial cable microwave mode, at a plate i, with frequency ω is given by where n i−1 (ω) is the occupation for the plate with the next higher temperature, n 0 = n BE (293 Kelvin, ω), A i is the attenuation at the plate i with temperature T i , such that for an attenuator with value -10 dB, A = 10.
The noise temperature of modes coming towards the device from the input is then n th,in < 0.108, which is only an upper bound as additional attenuation in the cabling and connectors has not been taken into account.For modes coming towards the device from the output wiring, the noise is given by the 50 Ω Johnson noise in the isolator n th,out ≃ 0. The thermal occupancy of the resonator n th is the sum of the occupancy of each thermal bath to which the resonator is coupled to, weighted by the coupling to each bath.Since the device is over-coupled to the feedline, the average occupancies of modes coming towards the device from both directions is a good approximation of the resonator occupancy: n th < 0.05 ≪ 1/2.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.

C. Frequency-dependence of damping
By sweeping the current in a coil located under the device we can change the static magnetic field traversing the device.Through this field, we are able to tune the Josephson inductance of the SQUID and hence the resonance frequency of the circuit.In the data presented in this paper, we operate at the first order magnetic field insensitive point of the SQUID -also known as flux "sweet-spot".In Fig. S13, we show our analysis of the data presented in Ref. [51], where S 21 is recorded as the flux through the SQUID is changed.This data set is acquired at a power of -138 dBm, far below the power at which the apparent nonlinear damping manifests.We find a change in the depth of the resonance peak as a function of detuning from the sweet-spot.However, the detuning required to match the change induced by power far exceeds any change in resonance frequency occurring through the Kerr nonlinearity.Indeed, the shift induced by driving the Kerr oscillator remains under 2 MHz (see Fig. S3) but the minimum changes by 4 dB, whereas the changes in Min|S 21 | shown in Fig. S13 in the detuning range of a few MHz remains smaller than 1 dB, dominated by experimental noise.This means for example that any frequency-dependence in the measurement port admittances as seen from the resonator is not responsible for the non-linear damping discussed in this work.

FIG. 1 .
FIG. 1. Superconducting Kerr oscillator circuit.a,The Kerr oscillator is constructed from an inductor, a capacitor and a SQUID (which behaves as and is depicted by a single 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 effectively behaves in this context as a single junction.

FIG. 2 .
FIG. 2.Observation of a resonator steady-state response suggesting nonlinear damping.a, Measured transmission magnitude |S21| (dots) for different drive powers.Whilst the shift in resonance frequency is expected from a classical analysis of the damped driven 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π × 193 kHz to κ nl int = 2π × 255 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.At the three highlighted points, the expectation values of photon number |a| 2 (where the minimum of |S21| is achieved) are 1.1, 6.8 and 13.2.
FIG. 4.Phase space picture of apparent damping.Here phase space operators are defined by â = (x + ip) / √ 2, see Supplementary Information Sec.S5A for further details.a, Wigner distribution of the steady-state in absence of Kerr nonlinearity (driven at ωr with Pin = −124 dBm).The balance between quantum noise, damping and drive is shown by vectors corresponding to Wigner currents.b, Growth of phase uncertainty of a coherent state under Kerr nonlinearity.The amplitudedependent 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 resonant amplitude's magnitude |⟨â⟩|.c, Wigner distribution of the steady-state with Kerr nonlinearity (at minimum |S21| with Pin = −124 dBm).The Kerr effect is eventually balanced by the damping, quantum noise and drive.Since the drive now opposes both damping and Kerr effect, 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 resonant amplitude's magnitude |⟨â⟩|.The center of mass (⟨â⟩) (white dot) is compared to the classical steady-state (white cross).Since the Wigner current of the Kerr effect 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 Sec.S5B).
FIG. S1.Transmission coefficient of the superconducting resonant circuit.Measurements of the transmission coefficient S21, processed following Sec.S2, are used in a curve fitting routine that seeks, for each driving power in the dataset and n th = 0, a numerical solution of the steady-state equation (3) of the damped driven quantum Kerr oscillator with the set of oscillator parameters that fits best to the measured data in accordance to a minimization of the error in: a, the magnitude |S21| and b, phase Arg(S21) of S21 as a function of driving frequency, as well as in c, the minimum Min|S21| of |S21| and d, the driving frequency ωmin at which this minimum occurs as a function of drive power.Data-points are rendered with blue dots and curve fits with solid lines.Each trace in a and b corresponds to a different drive power.At the top trace the drive power is Pin = −135.23dBm, then this increases 1 dBm per trace, reaching Pin = −118.23dBm at the bottom trace.We offset every trace in a and b by -10 dB and -2 rad respectively to ease examination.
FIG.S2.Quantum vs. classical descriptions: transmission coefficient.Data-points (blue dots) of the transmission coefficient S21, processed following Sec.S2, are compared to the best curve fits resulting from both a numerical solution of the Lindblad equation (S8) (solid line), and a numerical solution of the classical steady-state equation (S14) with nonlinear damping (dotted line), following, respectively, the minimization routines described in Sec.S2 and Sec.S3.We also display the solution to the steady-state equation (S13) of the linearly damped driven classical Kerr oscillator (dashed line).a, c, e The magnitude |S21| and b, d, f phase φ[S21] = Arg(S21) of the transmission coefficient are plotted as a function of the drive frequency.Panels (a,b), (c,d) and (e,f) correspond to powers -135 dBm, -127 dBm and -124 dBm respectively; all of them associated with a steady-state oscillator below the bistability threshold.
FIG.S3.Quantum vs. classical descriptions: minimum of |S21|.a, The minimum Min|S21| of transmission |S21| and b, the frequency ωmin which minimizes it, are plotted as a function of the drive power.We obtain excellent agreement between experimental data (dots), the fit based on a numerical solution of the Lindblad master equation (solid line), as described in Sec.S2, and the fit based on a numerical solution of the classical nonlinear damping model (dotted line), as described in Sec.S3.On the contrary, in a the minimum of transmission |1 − κext/κ| that follows from the solution to the classical steady-state equation (S13) with linear damping (dashed line) contrasts sharply with the experimental data.In b, due to the overlap between the solid and dotted lines, the inset shows their difference as a function of drive power.The drive power is bound to Pin ≤ −124 dBm, in which case the steady-state of the oscillator lies below the bistability threshold.

Fig
Drive power P in (dBm) Drive power P in (dBm) D Detuning (MHz)

10 n
FIG. S5.Test of the analytical model for |⟨â⟩⋆|.Numerical calculations of |⟨â⟩⋆| (full curves) are compared to the analytical expressions of Eq. (S58) (dashed curves).In panel a we use the same regime of parameters as the measured device, notably the same Kerr parameter K, with both no thermal population n th = 0 and half a quantum of thermal population n th = 1/2.In panel b, we show that the accuracy of the analytical model increases as we decrease K by a factor 10. For the calculations with thermal population, good agreement is found in the regime n th ≪ |⟨â⟩| 2 covered by our assumptions.

3 FIG
FIG. S6.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.
FIG. S7.Wigner current in an off-resonantly driven harmonic oscillator.We overlay the harmonic a, drive b, damping c, and diffusion (quantum noise) d Wigner currents on the Wigner distribution of the steady-state.e, Contrary to the resonant case, the sum of diffusion and damping no longer cancel the drive.Instead, these currents point away from the direction of the harmonic evolution, such that the net total current is 0. f, The remaining current moves points of the Wigner distribution along lines of equal probability such that ∂W/∂t = 0 .

20 FIG
FIG. S8.Wigner current in the steady-state of the driven Kerr oscillator.We overlay the contributions of the harmonic a, anharmonic b, drive c, and environmental d terms of the Lindblad equation on the Wigner function of Fig.4e.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.4c.

)
Given a change in density matrix ∂ ρ ∂t , we can compute the change in amplitude from ∂|⟨â⟩| ∂t FIG. S9.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[47,48].The starting state is taken to be a coherent state of amplitude α (see Fig.4a), the amplitude corresponding to the classical solution to the problem.The simulation is carried out with a power Pin = −124, until the steady-state shown in Fig.4cis 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 longer parallel to the driving current (see Fig.4c), 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.

W
FIG. S11.Quadrature squeezing.The angle θ defines the orthonormal coordinate set u, v with u = cos(θ)x + sin(θ)p where the uncertainty in û named ∆u is minimized.For each driving power, we plot the minimum ∆u in a and the angle which achieves this minimum in b.The uncertainty ∆u is expressed relative to the quadrature uncertainty of a coherent state shown as a dashed line.c, Wigner function plotted in the u, v basis for Pin = −124.2dBm where the smallest ∆u is achieved.On the left and on top, we plot the marginal distributions for u and v respectively (blue), overlayed with the distributions of a coherent state (dashed lines) with the same amplitude ⟨â⟩ as the plotted Wigner function.In this case, the uncertainty ∆u for the steady-state is 83 percent of the uncertainty one would obtain for a coherent state, demonstrating quadrature squeezing.
FIG.S12.Cryogenic wiring diagram.On the input side, at each temperature stage (corresponding to a plate of our dilution refrigerator) the signal passes through a thermalized attenuator.This reduces the noise photon occupation ni.On the output side, the signal passes through two circulators which isolate the device from noise from the HEMT amplifier at 3 K.