Four-wave-cooling to the single phonon level in Kerr optomechanics

The field of cavity optomechanics has achieved groundbreaking photonic control and detection of mechanical oscillators, based on their coupling to linear electromagnetic modes. Lately, however, there is an uprising interest in exploring cavity nonlinearities as a powerful new resource in radiation-pressure interacting systems. Here, we present a flux-mediated optomechanical device combining a nonlinear Josephson-based superconducting quantum interference cavity with a mechanical nanobeam. We demonstrate how the intrinsic Kerr nonlinearity of the microwave circuit can be used for a counter-intuitive blue-detuned sideband-cooling scheme based on multi-tone cavity driving and intracavity four-wave-mixing. Based on the large single-photon coupling rate of the system of up to $g_0 = 2\pi\cdot 3.6\,$kHz and a high mechanical quality factor $Q_\mathrm{m} \approx 4\cdot 10^{5}$, we achieve an effective four-wave cooperativity of $\mathcal{C}_\mathrm{fw}>100$ and demonstrate four-wave cooling of the mechanical oscillator close to its quantum groundstate, achieving a final occupancy of $n_\mathrm{m} \sim 1.6$. Our results significantly advance the recently developed platform of flux-mediated optomechanics and demonstrate how cavity Kerr nonlinearities can be utilized for novel control schemes in cavity optomechanics.


INTRODUCTION
Cavity optomechanical systems are the leading platform for the detection and manipulation of mechanical oscillators with electromagnetic fields from the nano-to the macro-scale 1 . Displacement detection with an imprecision below the standard quantum limit 2,3 , sidebandcooling to the motional quantum groundstate 4,5 , the preparation of nonclassical states of motion [6][7][8][9] , quantum entanglement of distinct mechanical oscillators 10,11 , topological energy transfer using exceptional points 12 and microwave-to-optical frequency transducers 13,14 are just some of the highlights that have been reported during the last decade. Essentially all of these impressing results have been achieved with linear cavities and linear mechanical oscillators, but the exploration of intrinsic cavity nonlinearities, often considered undesired and parasitic in optomechanics as they impose limitations on the maximally achievable multi-photon coupling rate 15 , has attracted increasing interest lately [16][17][18][19][20][21][22] . An exciting new scheme to couple a mechanical oscillator to microwave photons in a superconducting LC circuit has very recently been realized: Flux-mediated optomechanical coupling [23][24][25][26] . In this approach, the displacement of a mechanical oscillator is transduced to magnetic flux threading a superconducting quantum interference device (SQUID) embedded in a microwave LC circuit as flux-dependent inductance 16,27,28 . Due to the scaling of the optomechanical single-photon coupling rate g 0 with the external magnetic transduction field in flux-mediated optomechanics 23,28 , record singlephoton coupling rates for the microwave domain have been reported [24][25][26] . In future devices, the optomechanical single-photon regime 29,30 or even the ultrastrong coupling to superconducting qubits 31 seem feasible. In addition to being a flux-tunable inductor, a SQUID simul-taneously constitutes a flexible and highly controllable Kerr nonlinearity, which is widely utilized in superconducting qubits 32 , Josephson parametric amplifiers 33 and four-wave-mixing based bosonic code quantum information processing 34 . Therefore, flux-mediated optomechanics is also an ideal platform for realizing and studying Kerr optomechanics and for the development of new detection and control schemes of mechanical motion.
Here, we implement a flux-mediated optomechanical device with a large single-photon couling rate of up to g 0 ≈ 2π ·3.6 kHz and demonstrate sideband cooling of the mechanical oscillator close to its quantum groundstate by intracavity four-wave mixing (FWM). By using a strong parametric cavity drive, we activate the emergence of two Kerr quasi-modes in the SQUID circuit and realize an optomechanical coupling of these quasi-modes to the mechanical oscillator by an additional optomechanical sideband pump field. The drive-activated Kerr-modes show enhanced properties such as a reduced effective linewidth compared to the undriven circuit and we achieve effective single-photon cooperativities C 0 10. Strikingly, we find that blue-detuned optomechanical sideband-pumping on one of the Kerr-modes leads to dynamical backaction with the characteristics of red-sideband pumping in a standard optomechanical system, in particular to positive optical damping. We use this FWM based bluedetuned optical damping to cool the mechanical oscillator extremely close to its quantum groundstate with a residual occupation of n m ∼ 1.6. Our results demonstrate how cavity Kerr nonlinearities can be used in optomechanics to achieve both, enhanced device performance and new control schemes for mechanical oscillators. At the same time they reveal the potential of flux-mediated optomechanics regarding low-power groundstate-cooling of mechanical oscillators and the future preparation of FIG. 1. A superconducting quantum interference cavity parametrically coupled to a mechanical nanobeam. a Optical micrograph of the circuit. Bright parts are Aluminum, dark parts are Silicon substrate. The LC circuit combines two interdigitated capacitors C with two linear inductors L, connected through a superconducting quantum interference device (SQUID) with total Josephson inductance LS = LJ/2. The circuit is capacitively coupled to a Z0 = 50 Ω coplanar waveguide feedline (top of image) with a coupling capacitor Cc and surrounded by ground-plane. Scale bar corresponds to 50 µm. b Scanning electron micrograph of the constriction-type SQUID, showing the two Josephson junctions and the mechanical oscillator as part of the loop released from the substrate. Inset shows a zoom-in to one of the nano-bridge Josephson junctions. c Circuit equivalent of the device. For the experiments, two magnetic fields can be applied. The field B ⊥ is oriented perpendicular to the chip plane and is used to set the flux bias working point of the SQUID Φ ⊥ . The parallel field B transduces mechanical displacement of the out-of-plane mode to additional flux ∆Φ = B l∆x threading the SQUID loop. d shows the sample integrated into a printed circuit board with two microwave connectors and mounted into a 2D vector magnet. The large split coil is used to generate B , a small single coil behind the chip generates B ⊥ . e Transmission response |S21| of the cavity at B = 25 mT and B ⊥ = 0. From a fit to the data, we extract the resonance frequency ω0 = 2π · 5.2673 GHz, the total linewidth κ = 2π · 380 kHz, and the external linewidth κe = 2π · 110 kHz. Data are shown as circles, fit as black line. e Resonance frequency ω0 vs magnetic flux Φ ⊥ , normalized to one flux quantum Φ0 at B = 25 mT. Circles are data, line is a fit. The two operation points for this paper are marked with stars and denoted "I" for ω0 ≈ 2π · 5.22 GHz and "II" for ω0 ≈ 2π · 5.17 GHz. Details on measurements and fits can be found in Supplementary Note 4. quantum states of motion.

The device
Our device combines a superconducting quantum interference LC circuit with a mechanical nanobeam oscillator embedded into the loop of the SQUID, cf. Fig. 1. Details on device fabrication are given in Supplementary Note 1. At the core of the circuit, the SQUID acts as a magneticflux-dependent inductance L S (Φ), where Φ is the total magnetic flux threading the 21×3 µm 2 large loop. For the tunable optomechanical coupling between the displacement of the mechanical nanobeam and the microwave circuit, two distinct external magnetic fields are required.
First, a magnetic field perpendicular to the chip surface B ⊥ is used to change the magnetic flux bias Φ ⊥ through the SQUID loop, allowing to tune the circuit resonance frequency ω 0 and flux responsivity F = ∂ω 0 /∂Φ. Secondly, a magnetic in-plane field B is used to transduce the out-of-plane displacement ∆x of the mechanical oscillator to additional flux ∆Φ = B l m ∆x, where l m = 18 µm is the length of the mechanical beam. To apply these two fields, the chip is mounted into a homemade 2D vector magnet, consisting of a large split coil for B and an additional small coil mounted below the chip for the generation of B ⊥ , cf. Fig. 1d. The whole configuration is placed in a cryoperm magnetic shielding and attached to the mK plate of a dilution refrigerator with a base temperature T b ≈ 15 mK. More details on the measurement setup are given in Supplementary Note 2.
We perform the experiments presented here at in-plane fields of B = 21 mT and B = 25 mT. Figure 1e shows the transmission response of the cavity for B = 25 mT and at the bias-flux sweetspot. It has a resonance frequency ω 0 = 2π · 5.2673 GHz, a total linewidth κ = 2π · 380 kHz and an external linewidth κ e = 2π · 110 kHz. Figure 1f shows how the cavity resonance frequency can be tuned by ∼ 150 MHz by changing the applied flux bias Φ ⊥ threading the SQUID loop. The curves and cavity parameters at B = 21 mT only deviate slightly from the ones given here, the corresponding additional data can be found in Supplementary Note 4. Due to an improved SQUID design and fabrication, the cavity flux responsivity F is increased by one order of magnitude compared to our previous results 23 , which leads to a significantly enhanced single-photon coupling rate where x zpf is the mechanical zero-point fluctuation amplitude.
The mechanical nanobeam, visible in Fig. 1b and released from the substrate in an isotropic reactive ion etching process using SF 6 plasma 35 , is 500 nm wide and 70 nm thick. From its total mass of m ≈ 1.9 pg and the resonance frequency of the out-of-plane mode Ω m ≈ 2π · 5.32 MHz, we get x zpf = 2mΩm ≈ 30 fm. For an inplane field of B = 25 mT, and the two flux-bias points Φ I and Φ II , cf. Fig. 1f, we obtain single-photon coupling rates g 0,I = 2π · 2.2 kHz and g 0,II = 2π · 3.6 kHz with F I = 2π · 300 MHz/Φ 0 and F II = 2π · 520 MHz/Φ 0 . For the smaller in-plane field of B = 21 mT, the g 0 -values are scaled accordingly, cf. Supplementary Note 5.
The final important parameter of the device is its Kerr nonlinearity, which at the flux sweetspot is K/2π = −30 kHz. For the two flux bias operation points I and II we obtain K I /2π = −40 kHz and K II /2π = −55 kHz, respectively. More details on the determination of the circuit parameters and their flux dependence can be found in Supplementary Note 4.

Driven Kerr-modes and dynamical Kerr backaction
Owing to the Kerr anharmonicity K, the application of a strong microwave drive tone close to the cavity resonance frequency ω 0 significantly modifies the cavity response to an additional probe field. In Fig. 2, we discuss this modified response in the presence of a parametric drive tone with a fixed frequency ω d , when the cavity is tuned to cross this drive tone by means of the bias field B ⊥ . For large detunings between cavity and drive, the circuit response S 21 exhibits a standard single-mode resonance lineshape. However, as the detuning ∆ d = ω d − ω 0 is reduced, the driven cavity susceptibility χ g (Ω) =χ p (Ω) 1 − K 2 n 2 dχ p (Ω)χ * p (−Ω) (2) deviates considerably from a single linear cavity, leading to the regime of parametric amplification and degenerate four-wave mixing, which is experimentally identified by the appearance of a second mode. Here, Ω = ω − ω d denotes the detuning between the probe field at ω and the parametric drive and .
The two Kerr quasi-modes, which we denote as signal and idler resonance, appear symmetrically around the drive with complex resonance frequencies where n d is the parametric drive intracavity photon number. These Kerr-modes have been observed and discussed also in the context of optical cavities and mechanical oscillators [36][37][38] . The signal mode can be identified by the shifted and significantly deepened cavity absorption dip and the idler mode by the resonance peak, indicating net transmission gain by Josephson parametric amplification.
With the activation of the quasi-mode state, we also obtain a highly stabilized effective resonance frequency and linewidth, while the bare cavity suffers from considerable frequency fluctuations due to flux noise. Due to the reduction of frequency fluctuations in combination with a saturation of two-level system losses by the parametric drive (cf. Supplementary Note 6), the effective cavity linewidth is reduced from the flux-noise broadened κ off ∼ 2π · 1.5 MHz to the driven κ on ≈ 2π · 340 kHz. An analysis of the signal mode resonance frequency and linewidth in the presence of the parametric drive is provided in Figs. 2c and d. Within a small region of flux bias values, the drive-tone induced Kerr shift compensates for the flux-noise induced frequency shifts by means of an internal feedback loop. Strikingly, this mechanism yields a stabilization of the driven resonance, which thereby becomes the natural choice of operation regime during the following experiments. In an optomechanical system, any intracavity field also acts back on the mechanical oscillator by altering its resonance frequency and decay rate, an effect known as dynamical backaction 2,39 . Therefore, the effect of the parametric drive to the mechanical oscillator also requires some careful consideration. From the linearized equations of motion for the mechanical amplitude fieldb and the intracavity fluctuation fieldâ in a single-tone driven Kerr cavitẏ with multi-photon coupling rate g α = √ n d g 0 and input fieldsζ,ξ i andξ e , the effective mechanical susceptibility FIG. 2. Activating the driven Kerr quasi-mode state and single-tone dynamical Kerr backaction. a displays color-coded the magnitude of the SQUID cavity response S21 for varying SQUID flux bias Φ ⊥ /Φ0 in the presence of a strong drive placed at ω d . The flux bias range corresponds to a small variation of Φ ⊥ around operation point I and the in-plane field is B = 21 mT. When the flux-tunable resonance frequency ω0, indicated as dashed line and labelled "Drive off", is far detuned from the drive tone, the cavity response exhibits a single broad absorption resonance. As the detuning between cavity and drive ∆ d = ω d − ω0 is reduced, the cavity response is significantly modified and the original resonance is developing into a double-mode structure. The appearance of these driven Kerr quasi-modes indicates the onset of parametric amplification and degenerate FWM in the SQUID circuit. We denote the two modes as signal and idler resonance with the resonance frequencies ωs and ωi, respectively. Arrow on the left indicates the position of the linescan shown in panel b. In addition to the linescan from a (red circles) and the result of the analytical response calculation (solid black line), we show the equivalent linescan without parametric drive (blue circles) and its corresponding theoretical response (dashed black line). The curves without parametric drive are offset by +5 dB for clarity. Panels c and d show the extracted resonance frequency ωs and effective linewidth κ of the signal resonance vs flux bias. Lines show the result of modelling the effective quantities with the driven Kerr cavity equations and taking into account flux-noise broadening and two-level systems. The regime of operation for the experiments reported below is indicated by dashed lines and shaded areas. In this regime, the linewidth is nearly constant with κ /2π ≈ 340 kHz. The width of the operation range corresponds to the flux noise standard deviation, which we estimate to be σΦ ∼ 5 mΦ0. Panel e illustrates the contributions to the dynamical Kerr backaction of the intracavity drive fields to the nanobeam. Optomechanical (OM) up-and downscattering induces cooling and heating/amplification to the mechanical mode, respectively, where gα is the multiphoton coupling rate and χg is the probe susceptibility of the driven Kerr oscillator. In addition, interference between up-and downscattered fields due to degenerate FWM has to be taken into account. f and g show the calculated optical spring and optical damping due to dynamical Kerr backaction. The two blue/red lines and shaded area correspond to g0/2π = (1.78 ± 0.1) kHz. The detuning range ∆ d is slightly increased compared to a-d. In the additional range, the backaction is plotted in gray. The device operation range is indicated by the shaded area in between the vertical dashed lines.
can be derived as for the weak-coupling and high-Q m limit, which is safely fulfilled for our mechanical oscillator with a linewidth of Γ m ≈ 2π · 13 Hz. The single-tone dynamical Kerr backaction with χ g = χ g (Ω m ) and χ g = χ * g (−Ω m ) has almost the same form as in linear optomechanics, but with a modi-fied cavity susceptibility χ g . A striking difference, however, is found in the terms A = −iKn dχp (Ω m ) and A = iKn dχ the drive is located around one mechanical frequency detuned from the cavity |∆ d | ≈ Ω m = 2π · 5.32 MHz, the backaction looks very similar to that of a linear cavity. However, when the drive and the cavity are nearresonant, the backaction is strongly dominated by the intracavity photon number and a Duffing-like behaviour can be observed with a sudden transition from high-to low-amplitude state at ∆ d ≈ −2π · 3 MHz. In the operation regime for the experiments described here, the drive-induced backaction for operation point I is small with Γ opt /2π ∼ −1 Hz and δΩ m /2π ∼ −5 Hz. Using the bare mechanical linewidth Γ m ∼ 2π · 13 Hz, the corresponding phonon occupation is therefore increased by about 10%, a detailed calculation and discussion of the resulting mechanical mode occupation is given in Supplementary Note 7. Due to the considerable cavity flux noise outside of the driven quasi-mode regime, we unfortunateley cannot experimentally access the dynamical Kerr backaction for the detuning range shown in Fig. 2. Nevertheless, with a larger single-photon coupling rate g 0 at operation point II and a stronger drive tone, we observe regimes of mechanical instability induced by the dynamical Kerr backaction, which are in excellent agreement with the prediction from the theory. The corresponding data and analysis are explained in detail in Supplementary Note 7. The presented formalism for the dynamical Kerr backaction can also directly be applied to the sideband-unresolved regime and explain the experimental findings of a recent experiment with a similar SQUID cavity optomechanical device 24 , cf. also Supplementary Note 7.

Multi-tone dynamical four-wave backaction
An interesting question arising now is how the Kerr quasimodes couple to the mechanical nanobeam, when an additional optomechanical pump tone is applied to one of the Kerr-mode sidebands. One might expect that the coupling to the mechanical oscillator is suppressed in this state, similar to the reduced impact of flux noise, as the Kerr-mode frequencies ω s and ω i display only a very weak dependence on flux through the SQUID. Fluctuations of the bare resonance frequency, however, lead to modulations of α d and parametric gain, and therefore will impact the mechanical oscillator by inducing changes in the radiation-pressure force. A straightforward way to investigate this setting experimentally is to apply an additional optomechanical pump tone on the red sideband of the signal resonance, i.e., with a pump frequency ω p ≈ ω s − Ω m . Once in this configuration, a weak probe signal around ω ≈ ω p +Ω m can be used to detect optomechanically induced transparency (OMIT) 41 and thereby characterize the optomechanical interaction. A detailed theoretical description as well as a discussion of the experimental findings for this red-sideband pumping setup is given in Supplementary Notes 12-15. A conceptually less straightforward and more exciting possibility is to pump the idler resonance on its blue sideband ω p ≈ ω i + Ω m , cf. Fig. 3a. A blue-detuned pump is commonly associated with amplification/heating due to the favoured Stokes-scattering to lower energy photons. The Kerr-mode susceptibility χ g close to the idler resonance, however, resembles that of an "inverted" mode. Any small intracavity field in the driven Kerr cavity experiences in addition a mirroring effect due to degenerate four-wave mixing with the parametric drive tone. The presence of the blue-sideband pump field enriches this situation even further. Then the Kerr cavity is effectively oscillating with Ω dp = ω d − ω p due to the presence of two strong fields, and effects arising from non-degenerate four-wave mixing can impact probe fields and mechanical sideband fields and finally also the OMIT response and the backaction to the mechanical oscillator. A clear signature of the parametric state and four-wave mixing is the appearance of optomechanically induced transparency in the probe response of the signal resonance, when the idler Kerr-mode is pumped on its blue sideband. Corresponding data are shown in Fig. 3b and c. Here and in stark contrast to the usual OMIT protocol, the frequency detuning between the idler bluesideband pump and the probe tone is not even close to the mechanical resonance frequency but given by Ω = ω − ω p ≈ 2Ω dp − Ω m . To first order, the observation of this transparency can be understood by considering the intracavity generated tones in addition to the ones that are sent externally. The parametric drive generates an intracavity field with amplitude α d at ω d , and the optomechanical pump at ω p generates an intracavity field with amplitude γ − . Just by this doubly-driven configuration, a third intracavity "pump" field is generated by degenerate FWM at ω + = ω p + 2Ω dp and we denote its amplitude as γ + . Therefore, when ω p = ω i + Ω m , the γ + -field is located at the red sideband of the signal resonance ω + = ω s − Ω m . The beating between a probe field at ω ≈ ω s and the γ + -field is then near-resonant with the mechanical oscillator and will drive it into coherent motion. A second beating component, which is driving the mechanical oscillator, originates from the beating of the γ − -field and the idler field of the weak probe itself, cf. Fig. 3a. These two are also near-resonant with the mechanical oscillator. Once in coherent motion, the mechanical oscillator generates sidebands to all intracavity field Fourier components, some of which interfere with the original probe tone causing the observed appearance of four-wave OMIT. To characterize the dynamical backaction imprinted by the intracavity fields on the mechanical oscillator in the presence of the α d , γ − and γ + fields, we measure the optomechanical transparency response for varying detuning δ p between the γ − -field and the idler-mode blue sideband, cf. Fig. 3. For each detuning, we determine the effective mechanical resonance frequency Ω eff and effective mechanical linewidth Γ eff from a fit to the transparency signal and subtract the intrinsic values Ω m and Γ m . The remaining contributions to the resonance fre-FIG. 3. Four-wave-OMIT and four-wave-backaction for optomechanical blue-sideband pumping of the idler quasi-mode. a shows the experimental protocol. The SQUID cavity is prepared in the quasi-mode state by a strong parametric drive (PD). In addition, we apply an optomechanical (OM) pump tone on the blue sideband of the idler resonance (IR) ωp = ωi + Ωm + δp. Finally, we use a weak probe tone around the signal resonance (SR) to detect optomechanically induced transparency. We repeat this scheme for varying detunings δp. b explains how this protocol to first order leads to coherent driving of the mechanical oscillator. By PD-induced intracavity 4WM, the OM pump (probe tone) gets an idler field on the opposite side of the drive, which has the right detuning to the probe tone (pump) ∼ Ωm to coherently drive the mechanical oscillator. c shows the signal resonance transmission S21 measured with the weak probe field (OM pump off). Circles are data, line is a fit. Vertical bars labelled with A, B, and C indicate zoom regions for the corresponding panels shown in c and ∆s = ω − ωs denotes the detuning between probe field and SR. d probe tone response (OM pump on) in three narrow frequency windows around ω ≈ 2ω d − ωp + Ωm for three different pump detunings δp, cf. panel a. Note that the frequency difference between OM pump and probe field is Ω ≈ Ωm − 2Ω dp , which implies that when the pump field frequency is reduced, the probe field frequency is increasing. Each probe tone response displays a narrow-band resonance, indicating optomechanically induced transparency (OMIT) via excitation of the mechanical oscillator. For each δp, we fit the OMIT response (lines in c) and extract the effective mechanical resonance frequency Ω eff = Ωm + δΩm and the effective mechanical linewidth Γ eff = Γm + Γopt. The contributions δΩm and Γopt, induced by dynamical backaction of all intracavity fields, are plotted in panels e and f as circles vs δp. The result of analytical calculations is shown as two solid lines with shaded area, where the range described by the lines captures uncertainties in the device parameters, cf. Supplementary Note 11. The dashed line shows the result of equivalent calculations without cross-mixing (non-degenerate 4WM) terms. f illustrates schematically one four-wave cross-mixing term that leads to the observed dynamical backaction. Hereby, two mechanical sidebands with frequency difference Ω dp = ω d − ωp and both, the PD and the OM pump, contribute to the interaction. quency and linewidth δΩ m and Γ opt , respectively, correspond to the optical spring and optical damping by the microwave fields. The result, shown in Fig. 3d and e, is quite surprising. Even though the optomechanical pump field is bluedetuned to all cavity resonances ω 0 , ω s and ω i , we observe dynamical backaction with characteristics resem-bling red-sideband pumping in linear optomechanical systems. Most strikingly, we find a positive optical damping, which is usually a clear signature for red-sideband physics and the basis for sideband-cooling of the mechanical mode 4 . We use a linearized, optomechanical multi-tone Kerr cavity model, and implement the hierarchy from the experiment α d γ ∓ â to reveal which interactions are responsible for the observed behaviour, cf. Supplementary Notes 8-10. The resulting effective mechanical susceptibility has still the same form as for a standard optomechanical system, and all the FWM contributions can be captured in J -factors in the dynamical four-wave backaction with g − = γ − g 0 , g + = γ + g 0 , χ g,− = χ g (Ω m ), χ g,α = χ g (Ω m + Ω dp ) and χ g,+ = χ g (Ω m + 2Ω dp ). Closedform expressions for the J are given in Supplementary Note 9. We identify non-degenerate four-wave mixing terms in the J -factors as the dominant origin of the observed backaction. These terms have contributions from the drive field α d , from one of the γ ± fields and couple any two distinct mechanical sidebands which have the frequency difference ±Ω dp , cf. Fig. 3f for a schematic of one of these terms. Hence, these terms correspond to intracavity cross-mixing based on α d and γ ± fields. Using independently determined system parameters, we find excellent agreement between the experimental data and the analytical model when we take these cross-mixing terms into account, cf. solid lines in Fig. 3d and e. If we take only the degenerate FWM terms into account, which are induced by the presence of α d , we find a small and nearly constant backaction for all δ p , cf. dashed lines.
Blue-detuned four-wave cooling close to the groundstate Positive optical damping is commonly related to cooling of the mechanical mode. Therefore, the blue-detuned pumping scheme described in Fig. 3 seems feasible to be utilized as a counter-intuitive, yet innovative, method to eliminate the residual thermal excitations in the mechanical resonator.
To characterize the mechanical mode temperature, we detect the upconverted thermal displacement fluctuations in the signal resonance output field with a spectrum analyzer. For this measurement, the SQUID cavity in the quasi-mode state is pumped with an optomechanical tone on the blue sideband of the idler mode. Using a probe tone, we then measure the signal mode response S 21 in a wide frequency range and the OMIT response in a narrow range. Finally, we detect the output spectrum in the same frequency window where the OMIT is observed. A collection of spectra for varying optomechanical pump power P p is presented in Fig. 4b. From a careful analysis of the combined data sets, cf. Supplementary Notes 8-13, the equilibrium phonon occupation of the mechanical oscillator as well as the phonon occupation resulting from four-wave-cooling can be inferred.
The mechanical oscillator is well thermalized to the mixing chamber base temperature and its residual phonon occupation at the lowest operation temperature T b = 15 mK is about n th m ≈ 70 − 90 phonons. With increasing optical damping caused by the blue-detuned pump tone, we observe a corresponding reduction of the initial thermal occupation and the cooling factor is determined by Γ opt , very similar to usual optomechanical sidebandcooling. The observed four-wave cooling is also very robust with respect to pump and drive strengths and we achieve at both flux bias operation points a final fourwave-cooled occupation extremely close to the quantum groundstate n m ∼ 1.6. Due to the high single-photon coupling rates, it requires only a small amount of effective sideband photons n γ = |γ − | 2 + |γ + | 2 10 to achieve these low occupations. The fact that we use strongly driven Kerr quasi-modes as cold bath, however, modifies the minimally achievable occupation. Due to Josephson parametric amplification of quantum noise in the quasi-mode state, the cavity will acquire an effective temperature, even if the bare cavity is in the quantum groundstate. This drive-induced cavity heating defines the cooling limit for the mechanical resonator. In the state we are operating here, the Josephson gain is small and the effective thermal occupation of the cavity is still considerably below 1. We estimate the current cooling limit due to amplified quantum noise to be ∼ 0.6, where the exact value depends on the drive strength n d and on the bias-flux operation point. With higher bias flux stability the cavity could be stabilized at a point where the Josephson gain is small enough to enable n lim m < 0.3. Achieving the lowest occupation in the current device requires a careful balancing of drive and pump strength and for the highest pump powers, we observe the onset of additional cavity shifts and line broadening, possibly related to drive depletion or higher-order nonlinear effects. With slightly optimized device parameters regarding K and g 0 , we should therefore be able to cool to n m < 1. We emphasize though, that the blue-detuned cooling scheme allowed to achieve a significantly lower phonon occupation than signal-mode red-sideband pumping. With a pump on the red signal-mode sideband, a second cavity bifurcation instability occurs at moderately high pump powers, as the red sideband pump is attracting the cavity, while the blue-detuned pump is repelling it. The related jump to a high-amplitude state with a different signal resonance frequency, prevents us from cooling below n red m ∼ 4. The corresponding red-sideband cooling data and analysis can be found in Supplementary Note 16.

DISCUSSION
The results we presented here demonstrate clearly that the young field of flux-mediated optomechanics is quickly advancing towards an exciting and competitive optomechanical platform, which intrinsically allows for novel FIG. 4. Blue-detuned four-wave-cooling of a mechanical oscillator close to its quantum groundstate. a Schematic representation of the experiment. A parametric drive is used to activate the quasi-mode state and an OM pump is sent to the blue sideband of the idler resonance ωp ≈ ωi + Ωm. The signal resonance output power spectral density is measured using a spectrum analyzer around ω = ωp + 2Ω dp + Ωm ≈ ωs. b Power spectral densities normalized to the optomechanical pump input power Pp for various pump powers. Frequency axis is given with respect to ω = ωp + 2Ω dp + Ωm. With increasing pump power, the linewidth of the upconverted mechanical noise spectrum is increasing, indicating four-wave dynamical backaction damping. Simultaneously, the area of the normalized signal decreases, indicating cooling of the mode. From fits (lines and shaded areas) to the data (points), we determine the resulting phonon occupation nm. In c we show the cooled phonon number vs Γ eff /Γm in a collection of several different datasets. Intracavity drive photon numbers vary between different points in the range 40 < n d < 100. Circles correspond to data from measurements at operation point I and squares to data from operation point II. Stars show the points that correspond to the data shown in b, taken at operation point I. All measurements have been taken at B = 25 mT. Inset shows the result of a thermal calibration measurement, indicating that the mechanical oscillator mode equilibrates with the fridge base temperature and the residual thermal occupation at T b = 15 mK is n th m ≈ 70 − 90. Dashed lines and shaded area display the theoretically calculated range of four-wave-cooled phonon occupation, taking into account a possible range of 60 ≤ n th m ≤ 100 and 45 ≤ n d ≤ 90. Parametric amplification of cavity quantum noise limits the minimally achievable phonon occupation in our parameter regime to n lim m ∼ 0.6. For the highest powers, we exceed this theoretical limit by only a factor ∼ 3. d shows the effective effective mechanical linewidth vs intracavity sideband photon number nγ = |γ−| 2 + |γ+| 2 for points from c, which have nearly constant n d ≈ 60 ± 10, demonstrating that we achieve significant cooling with a small number of photons. Line corresponds to theory with Γm = 2π · 15 Hz e shows an OMIT scan at the point of largest cooling with an effective linewidth Γ eff ≈ 2π · 1.5 kHz, which corresponds to an effective four-wave cooperativity of C fw 100. f shows the corresponding power spectral density in units of quanta with noise squashing due to a small, but finite effective temperature of the cavity by amplified quantum noise. Error bars in c consider uncertainties in the fitting procedure and in the bare mechanical linewidth, for details see Supplementary Note 14. ways of manipulating mechanical motion. Our device provides a large single-photon coupling rate of up to g 0 = 2π · 3.6 kHz and achieves large cooperativities of up to C fw > 100 for small numbers of intracavity photons. By using strong parametric driving, we show how the intrinsic Josephson-based Kerr nonlinearity can be utilized as a resource for improved sideband-resolution and frequency-stability and for the implementation of a novel four-wave-mixing-based phonon control scheme. In combination, these properties enabled us to use four-wavecooling in a Kerr cavity to prepare a MHz mechanical nanobeam resonator close to its quantum groundstate.
Future device improvements can be achieved by reducing the SQUID loop inductance further in order to increase the flux responsivity and the single-photon coupling rate. One order of magnitude is a feasible goal in this direction, as related platforms have already demonstrated such high responsivities 24,25 . This improvement alone would bring the device to a cooperativity of 10 4 and to the onset of the strong-coupling regime with g ∼ 2π · 150 kHz ∼ κ/2. With increased in-plane fields, up to ∼ 1 T with e.g. Niobium or granular Aluminum, those numbers could be improved by another order of magnitude. In the current device, however, the main limiting factor to achieve higher coupling rates and cooling the mechanical oscillator into the groundstate was external flux noise coupling into the SQUID in large in-plane fields. We suspect that the origin of this flux noise is in the vector magnet leads and the used current sources, respectively, or in parasitic out-of-plane components that lead to flux instabilities, vortex avalanches and microwave-triggered vortex motion in proximity to the SQUID. Flux noise in the leads and current sources could potentially be reduced by using a superconducting magnet in persistent current mode. And although our current setup can locally cancel parasitic out-of-plane fields, it cannot do so over the complete chip simultaneously due to the geometry of the small coil. A global compensation might be necessary, however, to completely avoid any flux instabilities arising from the out-of-plane fields, which can cause flux fluctuations also in large distances from their occurrence. Using intrinsic Kerr nonlinearities as a resource in optomechanical systems has just begun. Further interesting directions in Kerr optomechanics might involve intracavity squeezing, intracavity Josephson parametric amplification, intracavity cat-state generation, groundstate cooling in the sideband unresolved regime or enhanced quantum transduction. Significantly larger Kerr nonlinearites than the ones presented here, implemented in superconducting transmon qubits, have also been discussed recently for mechanical quantum state preparation 31,45,46 . Similar schemes investigating and exploiting the Kerr nonlinearity of SQUID circuits could furthermore be implemented naturally in the platform of photon-pressure coupled circuits 42-44 .  Here we present a step-by-step description of the device fabrication. The individual steps are schematically shown in Supplementary Fig. S1, where we omitted step 0, the patterning of the electron beam lithography (EBL) alignment markers, as well as the wafer dicing steps and the final device mounting. • Step 0: Marker patterning The fabrication of the device starts by the patterning of alignment markers on top of a 2 inch silicon wafer using electron beam lithography (EBL). The marker structures are patterned using a CSAR62.13 resist mask and a sputter deposition of 50 nm Molybdenum-Rhenium alloy. After undergoing a lift-off process, the only remaining structures on the wafer are the markers. The complete 2 inch wafer is then diced into individual 14 × 14 mm 2 chips, which are used individually for the subsequent fabrication steps. On each of these fabrication chips, we structure 2 device chips with dimensions of 5×10 mm 2 , each of which contain one coplanar waveguide microwave feedline and seven quantum interference LC circuits.

• Step 1: Junctions patterning
As first real step of the device fabrication we pattern two nanobridges (the later Josephson junctions) for each LC circuit using CSAR62.09, cf. Supplementary Fig. S1a. The two bridges of each pair of nanobridges forming one superconducting quantum interference device (SQUID) are hereby always identical. All bridges have a length of ∼ 100 nm but vary in width between 30 and 60 nm for different SQUIDs in order to compensate for small variations and uncertainties in final structure size and select the most suitable device during the experiment. The nanobridges also have two 700×1150 nm 2 large pads for achieving good galvanic contact to the rest of the circuit, which is patterned in fabrication step 2. After the EBL exposure, the sample is developed in Pentylacetate for 60 seconds followed by a 1:1 solution of MIBK:IPA (Methyl IsoButyl Ketone:IsoPropyl Alcohol) for another 60 seconds and finally rinsed in IPA. Once the resist is developed, the chip is loaded into a sputtering machine where a 15 nm think layer of Aluminum (1% Silicon) is deposited. After the deposition, the sample is placed horizontally at the bottom of a glass beaker containing a small amount of room-temperature Anisole and left in an ultrasonic bath for a few minutes. During this time, the remaining resist is dissolved and the Aluminum layer sitting on top is lifted off, the result is schematically shown in Supplementary Fig. S1b.

• Step 2: Microwave cavity patterning
After the junctions are patterned, we once again spin-coat the sample with CSAR62.13 and pattern the SQUID arms together with all the remaining superconducting structures. After the EBL exposure, the sample is developed as for the previous fabrication step and afterwards loaded into a sputtering machine. Hereby, the nanobridges themselves are covered and protected by resist, cf. Supplementary Fig. S1c. At this point and prior to the deposition of the second Aluminum layer, an Argon milling process is perfomed in-situ in order eliminate any oxide present on top of the contact pads. This measure is necessary to generate good electrical contact between the two layers. After the sputtering process of the second, 70 nm thick Aluminum (1% Silicon) layer, the sample undergoes an ultrasonic lift-off process similar to the one in Step 1, the result is shown schematically in Supplementary Fig. S1d. • Step 3: Nanobridge thinning by Ar ion milling In order to reduce the cross-section and the critical current of the nanobridges even further, we apply a short ion milling step to the SQUID at this point. To do so, we pattern and develop another layer of CSAR62.13 on top of the device as described in Steps 1 and 2, which protects the whole chip except for rectangular windows around the SQUIDs themselves, cf. Supplementary Fig. S1e. From test measurements, we observe that if we do not protect the rest of the circuit from the milling in this step, we obtain a significant reduction of the circuit quality factor, which we think might be due to ion implanation into the substrate. Note that with the milling parameters we use for this step, we do not get a directional milling, but mainly a narrowing of the nanobridges from the sides. This is also the reason why we need the contact pads in the first place. If we work with bare nanobridges in Step 1, they are milled away completely during the essential in-situ native oxide removal in Step 2.

• Step 4: Dicing
Right before the final release of the mechanical oscillator, the sample is once again diced to two smaller 5×10 mm 2 sized chips in order to fit into the sample mountings and the microwave PCB (Printed Circuit Board). The remaining 2 mm at each edge of the original 14 × 14 mm 2 large chip is only a margin for the fabrication and is disposed of.

• Step 5: Mechanical beam release
For the final EBL step, a CSAR62.13 resist was once again used as mask and the development of the pattern was done in a similar way as for the first two layers. Once the etch mask, consisting of a small window close to the outer side of the SQUID loop (cf. Supplementary Fig. S1f ), is patterned, the sample undergoes an isotropic, reactive ion etching process in SF 6 at a sample temperature of ∼ −10 • C for two minutes S1 . During this time the Silicon substrate under the SQUID arm/the mechanical beam is etched without attacking the aluminum layer forming the cavity and the mechanical beam. Once the beam is released, we proceeded with an O 2 plasma ashing step in order to remove the remaining resist from the sample. At this point the fabrication is completed, the result is shown schematically in Supplementary Fig. S1g.

• Step 6: Device mounting
After the fabrication, the sample is glued into a microwave printed circuit board (PCB) using GE varnish and wirebonded both to ground and to 50 Ω connector lines. An optical image of the chip and the PCB, both mounted into a magnet, is shown in Fig. 1 of the main paper.

Setup configuration
The experiments reported in this paper were performed in a dilution refrigerator with a base temperature T b ≈ 15 mK. Within the outer vacuum can of the system, a mu-metal shield is installed to provide basic magnetic shielding for the whole sample space from the 3 K plate to the mK plate. A schematic diagram of the experimental setup and of the external measurement configuration used in the reported experiments can be seen in Supplementary Fig. S2. The PCB, onto which the fabricated sample was glued and wirebonded, is mounted into a 2D vector magnet casing and connected to two coaxial lines. The complete configuration including the vector magnet is placed in a magnetic cryoperm shield. The vector magnet combines two distinct superconducting magnets, a small one for the generation of an out-of-plane field and a larger split coil for the in-plane field. The coils are used to independently generate a magnetic field in the two different directions by providing a DC current to the corresponding coil. A more detailed information about the design and setup of the vector magnet is provided in the following subsection.
Since the optomechanical circuit that we present in this paper was designed in a side-coupled geometry, the input and output signals were sent/received through separate coaxial lines in order to measure the transmission spectrum of the feedline to which the system is coupled. The input line is heavily attenuated in order to balance the thermal radiation from the line to the base temperature of the fridge and the output line contains a cryogenic HEMT (High-Electron-Mobility Transistor) amplifier working in a range from 4 to 8 GHz and two isolators to block the thermal radiation from the HEMT to reach the sample.
Outside of the refrigerator, we used a single measurement scheme for all the different experiments. The VNA was used to measure the response spectrum S 21 of the electromechanical system, one microwave generator sends a coherent signal at ω d as parametric drive for the SQUID cavity and the second microwave generator sends a tone at ω p as optomechanical pump for the parametrically driven cavity. Finally, a spectrum analyzer was used to record the output power spectrum around the cavity resonance.
For all experiments, the microwave sources and vector network analyzers (VNA) as well as the spectrum analyzer used a single reference clock of one of the devices. Vector magnet design Figure 1 of the main paper shows photographs of the sample mounted on the PCB and fixed in the vector magnet bobbin. The two large parallel coils on each side of the sample are wound from a single wire (niobium-titanium in copper-nickel matrix) and in the same orientation and therefore form a Helmholtz-like split coil (the distance between the coils is slightly larger than their effective radius), which creates a nearly homogeneous in-plane magnetic field at the location of the device. At room temperature the coil has a resistance of R ≈ 6 kΩ, which approximately corresponds to 2000 windings of superconducting wire on each side. From the coil geometry and the number of windings, we estimate the current-to-field conversion factor to be 70 mT/A. On the backside of the sample/PCB platform within the magnet bobbin is a second small coil mounted for providing the out-of-plane magnetic field used to tune the SQUID flux bias point, cf. main paper Fig. 1. This out-of-plane coil can also be used to compensate for a parasitic out-of-plane component of the in-plane field due to misalignments of the sample/PCB with respect to the in-plane field axis (estimated to be around 2 • − 3 • from the SQUID flux response). For in-plane fields B 25 mT, however, the compensation is not yet critical. For larger in-plane fields, vortices start to penetrate the film and there is a dramatic reduction in the cavity quality factor observable. The room-temperature resistance of the out-of-plane coil is R ⊥ ≈ 120 Ω which corresponds to approximately 400 turns of superconducting wire and to a conversion factor of 1 mT/A. The superconducting wires leading to each of the coils from the 3 K plate are twisted in pairs, in order to reduced the amount of captured flux noise. Furthermore, since the critical temperature of the wire is about ∼ 12 K, the wires can go unbroken until the 3 K stage. Above this plate, the wires are no longer superconducting and therefore a transition to normal conducting wires is required. For this, we connected each of the superconducting in-plane coil wires to 9 wires of a 24-line copper loom provided by Bluefors and each of the out-of-plane coil wires to 3 wires of the loom. From the 3 K stage until room temperature the current flows in parallel through the respective loom wires, decreasing the additional heat load on the plate. With this approach we are able to send I ∼ 0.5 A through the in-plane coil without any considerable heat added to any of the plates and maintaining the fridge base temperature. At room temperature we are left with 4 cables, two for each coil, which are used with individual directed current (DC) sources to independently generate the magnetic fields.

SUPPLEMENTARY NOTE 3: POWER CALIBRATION
In order to estimate the input power on the on-chip feedline of the device, we use the thermal noise of the HEMT (High-Electron-Mobility Transistor) amplifier as calibration method. The cryogenic HEMT amplifier thermal noise power is given by where k B is the Boltzmann constant, T HEMT is the noise temperature of the amplifier, which, according to the specification datasheet, is approximately 2 K, and ∆f = 2000 Hz is the measurement IF bandwidth. The calculated noise power is P HEMT = −162.6 dBm, or as noise RMS voltage ∆V = 1.66 nV. Taking into account the room temperature attenuators of 60 dB as well as additional 3 dB of room-temperature cable losses between the VNA output and the directional couplers for the pump tones and assuming an attenuation between the sample and the HEMT of 2 dB we extract a frequency-dependent input attenuation for the pump tones as shown in Supplementary Fig. S3. In addition and for confirmation, we perform a fixed-frequency measurement of the signalto-noise ratio using the pump signal generator itself and a spectrum analyzer for selected frequency points around 5.22 and 5.17 GHz. We observe agreement between the two methods better than 0.5 dB.

Circuit model
A simplified circuit equivalent of the SQUID cavity used in this experiment is shown in Supplementary Fig. S4a. We model it as a simple parallel RLC circuit capacitively coupled by a coupling capacitance C c to a microwave feedline with characteristic impedance Z 0 as shown in b, cf. also Ref. S2 . The resistance in this model captures all intracavity losses. The resonance frequency, external and internal linewidth of the circuit shown in b are given by respectively. Each of the two physical capacitors in the main circuit, cf. main paper Fig. 1, is an interdigitated capacitor (IDC) with N = 148 fingers, each 100 µm long and 1 µm wide. With the gap between two fingers of also 1 µm and the relative permittivity of the Silicon substrate r = 11.7, we obtain for each of the IDCs C ≈ 824 fF using the analytical . For each frequency point, we determine from the 501 traces the signal-to-noise ratio and with the assumption of a frequency-independent HEMT noise temperature and 2 dB losses between the sample and the HEMT, we get the input line attenuation as plotted. The gray area shows where the cavity was during the calibration. Due to its presence, the attenuation in this range can not be considered a reliable value. Our experiments, however, mainly take place around 5.22 GHz and 5.17 GHz (labelled with I and II, respectively) and therefore the presence of the cavity at around 5.26 GHz does not lead to any calibration problems. We also note, that we observe almost identical amplitude oscillations in the transmitted signal, indicating that we are indeed dealing with strong cable resonances. expressions provided in Ref. S3 . The total capacitance is then approximately given by C tot = 2C + C c ≈ 1.65 pF, where we included also the (mostly negligible) coupling capacitance C c ≈ 6.5 fF. The value for C c was obtained via the external cavity linewidth of κ e ≈ 2π · 110 kHz, the feedline characteristic impedance Z 0 = 50 Ω and the resonance frequency ω 0 = 2π · 5.267 GHz. Using the resonance frequency, we can also estimate the total inductance as L tot = 1 In the linear regime, a capacitively side-coupled LC circuit is described by the S 21 response function with detuning of the probe tone from the resonance frequency and the internal and external linewidths κ i and κ e , respectively. Implicitly, we assume symmetric coupling to the left and right feedline part in this relation. Due to considerable cable resonances in our setup, however, this assumption might be not strictly valid. We also observe, that for a consistent modeling of all our datasets, small adjustments to κ e in different experimental situations are leading to higher agreement between data and theory. The different microwave components in the setup (cables, attenuators, directional couplers, isolators etc) affect the ideal cavity transmission spectrum by amplitude and phase modulations, and we consider a modification in the response function by introducing a frequency-dependent complex-valued microwave background. The modified cavity response is written as where we consider a frequency-dependent complex background and an additional, possible interference rotation of the resonance circle around its anchor point with the phase factor e iθ . In our fitting routine the background is extracted by first excluding the cavity resonance from the response and fitting the remaining data with Eq. (S6). After complex division of the data with the background model, the remaining cavity response is fitted independently. As final step the original data are fitted with the full function for S 21 including the background again using the obtained fit values from the first two independent fits as starting values for the full fit. From the final fit, we remove the background of the full dataset by complex division for the resonance data shown this paper. Also, we correct for the additional rotation factor e iθ . In Supplementary Fig. S5, we show an exemplary fit of the cavity response around resonance as raw data and as background-corrected data in both, the complex plane and in the magnitude of S 21 . From the fit to the data, taken at B = 25 mT and B ⊥ = 0 (the sweetspot), we obtain ω 0 = 2π · 5.2672 GHz, κ i = 2π · 269 kHz and κ e ≈ 2π · 111 kHz.

The SQUID Josephson inductance
The total flux Φ in a superconducting quantum interference device (SQUID) with non-negligible loop inductance L loop is given by where Φ b is the bias flux by external magnetic fields, J is the screening current circulating in the SQUID loop and Φ 0 = h 2e = 2.07 · 10 −15 Tm 2 is the flux quantum. Note that L loop contains both, the geometric and the kinetic inductance contribution to the inductance of the SQUID loop. In the absence of a bias current and for identical Josephson junctions with a sinusoidal current-phase relation, the circulating current is given by with the zero-flux-bias of a single junction I c . Using the screening parameter β L = 2L loop Ic Φ0 , the relation for the total flux can be written as We use this equation to numerically calculate the total flux in the SQUID for a given external flux.
With the total flux in the SQUID known, the Josephson inductance of a single junction (S10) and the total Josephson inductance of the SQUID can be determined.

Cavity field dependence
Using the flux-dependence of the SQUID Josephson inductance and our simplified circuit model, the resonance frequency of the cavity as function of the perpendicular bias flux Φ ⊥ can be written as with the linear inductance participation ratio Λ = L L + L J0 (S13) and the total flux in the SQUIDΦ The zero-bias junction inductance is hereby given as L J0 = L J (Φ = 0). The first experimental step to fit the flux-dependence of the cavity resonance frequency and to determine Josephson inductance L J and screening parameter β L is a calibration of the bias flux axis and to find the current-to-flux conversion for the small coil generating Φ ⊥ , respectively. Supplementary Fig. S6a shows as circles the experimentally obtained resonance frequencies at B = 0 for a sweep of the bias flux Φ ⊥ . The dataset combines the data points obtained during a bias flux upsweep and a downsweep. This is necessary as the SQUID has a non-negligible loop inductance, which leads to a hysteretic flux response S2,S4,S5 . The distance between two neighboring flux archs corresponds to one flux quantum Φ 0 and via this procedure the current-to-flux conversion is obtained. Subsequently, the flux-dependence of ω 0 can be fitted using Eqs. (S12) and (S9). From the fits, we obtain the zero-bias junction critical current I c and the screening parameter β L , the corresponding fit curves are shown as lines in Supplementary Fig. S6a. From the fit at zero in-plane field we get I c ≈ 2.6 µA and a screening parameter β L ≈ 0.7. Using L J0 = Φ0 2πIc , we get for the inductance of a single Josephson junction L J0 ≈ 127 pH, which corresponds to a linear inductance participation ratio Λ ≈ 0.89 and a total SQUID loop inductance L loop ≈ 278 pH. Enabling the optomechanical coupling between the nanobeam and the SQUID cavity requires an additional in-plane magnetic transduction field, and therefore we also record the resonance frequency flux-dependence at the in-plane fields of B = 21 mT and B = 25 mT, where we operate the device for the optomechanical experiments. The result is shown in Supplementary Fig. S6b as circles. From the data, we observe a small decrease of the sweetspot resonance frequency with increasing B . In addition, we observe a slight widening of the flux arch with increasing B , indicating a nonlinear increase of the kinetic contribution to the SQUID loop inductance and a consequently increased β L . From the fits, we get for both in-plane fields a slightly reduced critical junction current I c ≈ 2.2 µA and slightly increased screening parameters β L,21 ≈ 0.79 and β L,25 ≈ 0.82. These value correspond to Λ ≈ 0.865, L loop,21 ≈ 371 pH and L loop,25 ≈ 385 pH. We observe that the loop inductance seems to increase by more than both, the Josephson inductance and the linear circuit inductance due to the in-plane field. Our suspicion is that this effect is caused by a modification of the nanobridge current-phase relation in the in-plane field, but for a final conclusion more experiments would have to be conducted.
For the optomechanical multi-photon interaction two more quantities of the SQUID cavity and their flux dependence are highly important. The first is the flux responsivity F = ∂ω 0 /∂Φ b , i.e., the change of resonance frequency with change of bias flux through the SQUID loop. It is directly proportional to the optomechanical single-photon coupling rate g 0 , cf. Supplementary Note 5. The responsivity is identical to the slope of the flux tuning curve shown in Supplementary Fig. S6b and the numerically obtained results for both, experimental data and the fit curve, are shown in Supplementary Fig. S7a. The bias-flux operation points relevant for this paper are marked with a dotted and dashed line, respectively, and labelled as "I" and "II". The corresponding flux responsivities are F I ≈ 2π ·300 MHz/Φ 0 and F II ≈ 2π · 520 MHz/Φ 0 , respectively, and nearly identical to each other for the two chosen in-plane fields. The second important quantity is the Kerr anharmonicity related to the nonlinear Josephson inductance of the SQUID. It is given by and depends in addition on the in-plane field via the in-plane dependence of the nanobridge critical current or Josephson inductance, respectively. The result of this calculation, based on the flux arch fits of Supplementary Fig. S6b is shown in Supplementary Fig. S7b. The dependence of the anharmonicity on flux bias Φ ⊥ shows a very similar trend for all in-plane fields with different starting values at the sweetspot Φ ⊥ = 0. The completely unbiased cavity has K ≈ −2π·17.5 kHz. For an in-plane field of B = 21 mT, we obtain K I,21 ≈ −2π·39 kHz and K II,21 ≈ −2π·54 kHz at the operation points "I" and "II", respectively, and for B = 25 mT, we find K I,25 ≈ −2π·41 kHz and K II,25 ≈ −2π·56 kHz. As the difference between the two in-plane field is small and subject to uncertainties due to uncertainties in the circuit parameters, we will work with the same approximate anharmonicities for both in-plane fields of K I ≈ −2π · 40 kHz and K II ≈ −2π · 55 kHz.

SUPPLEMENTARY NOTE 5: THE OPTOMECHANICAL SINGLE-PHOTON COUPLING RATE
The optomechanical single-photon coupling rate in flux-mediated optomechanics is given by S2,S6 g 0 = γFB l m x zpf (S16) where F is the cavity frequency flux responsivity, B is the in-plane magnetic field, l m is the length of the mechanical nanobeam and γ is a scaling factor on the order of unity taking into account the mode shape of the beam. The zero-point fluctuation amplitude of the mechanical displacement is given by where m is the effective mass of the beam and Ω m its resonance frequency. For our nanobeam, cf. main paper Fig. 1b, we get from a scanning electron micrograph the length as approximately l m = 18 µm. Using the beam film thickness of ∼ 70 nm, its width of ∼ 500 nm and the density of Aluminum ρ Al = 2700 kg m −3 , we get a total mass of m = 1.7 · 10 −15 kg. In the experiment, we observe a mechanical resonance frequency Ω m ≈ 2π · 5.32 MHz and therefore, we find a zero-point fluctuation amplitude of x zpf ≈ 30 · 10 −15 m = 30 fm.
With the flux responsivities shown in Supplementary Fig. S7a and assuming γ ≈ 1, we calculate the corresponding single-photon coupling rates g 0 , the result is shown in Supplementary Fig. S8. For the different operation points, we obtain single-photon coupling rates g 0 between g min 0 ≈ 2π · 1.78 kHz (point I) and g max 0 ≈ 2π · 3.57 kHz (point II). For all presented results, we will add the corresponding coupling rates either in the figure legend or in the caption.
where the intracavity field α is normalized such that |α| 2 = n c is the intracavity photon number and |S in | 2 corresponds to the input field photon flux.

Single-tone solution
If the cavity is driven with a single tone at frequency ω d , the input field is given by and for the intracavity field we make the Ansatz with real-valued S in and α d . Any phase difference between the drive and the intracavity field is captured in the drive phase φ d . Inserting drive and intracavity field Ansatz into the equation of motion gives where ∆ d = ω d − ω 0 is the detuning between drive and undriven cavity resonance frequency. Multiplication of this equation with its complex conjugate leads to the determination polynomial for the drive intracavity photon number In general, this polynomial has three roots for n d . The real-valued solutions correspond to physical states and for certain parameters all three solutions are real-valued. This regime corresponds to the bifurcation regime, where two of the three oscillator states are stable, one low-and one high-amplitude solution. The middle solution is unstable and irrelevant for the experiments described here. The phase difference between drive and intracavity oscillations can be found via Nonlinear cavity response modeling We measure the power dependence of the SQUID cavity response S 21 by means of a vector network analyzer at the bias flux sweetspot and at B = 25 mT. The result is shown and discussed in Supplementary Fig. S9, where in a the magnitude and in b the phase of the complex S 21 is plotted. The data shown have been background-corrected as described above using the background fit of the first line. To model the response, we employ the single-tone model described in the previous section and calculate the response via where P d is the drive power on the on-chip microwave feedline. For the model, we use the Kerr anharmonicity obtained from the independent cavity modeling K = −2π · 30 kHz, the resonance frequency ω 0 = 2π · 5.2672 GHz and the corresponding external linewidth κ e = 2π · 106.5 kHz as obtained from the lowest-power response of the current dataset.
As apparent from the increase in resonance absorption dip depth and the reduction of total linewidth with increasing drive power, we also have to consider nonlinear dissipation, that decreases with increasing power. As first step in the nonlinear resonance analysis, we fit each of the nonlinear response curves using the Kerr polynomial and using a single decay rate for each power as fit parameter. The result of this procedure for κ is shown in Supplementary  Fig. S9c, where we plot the fit linewidth vs the photon number n d on resonance. We model this decrease of linewidth with the functional dependence for two-level systems where n crit describes the characteristic photon number for two-level saturation. From a fit to the linewidth data we obtain κ 0 = 2π · 209 kHz, κ 1 = 2π · 145 kHz, and n crit ≈ 0.26, cf. line in Supplementary Fig. S9c. As next and final step, we implement this analytical function into the Kerr polynomial and solve for the final photon number at each drive power and detuning. To find convergence in the solution for the photon number n d due to the power dependent κ we have to iterate the polynomial solution multiple times in this approach for each frequency point, feeding back in each iteration the κ(n d ) from the previous iteration. The result is added as lines in Supplementary Fig. S9a and b and shows good agreement between experimental data and theoretical modeling. The only free parameter used in the end is the in-fridge attenuation between the VNA output and the sample, and we find best agreement for choosing G STK = −89.2 dB. This value is very close to the independently estimated probe attenuation of G SNR = −89.5 dB at the corresponding frequency, cf. Supplementary Fig. S3 with consideration that the probe attenuation is 3 dB larger than the shown pump attenuation.

Linearized two-tone solution
If the Kerr cavity is driven by a strong drive tone S d and a weaker second tone S p with frequency ω p , the total input field is given by where we again choose S d to be real-valued. As Ansatz for the intracavity field, we then use where γ − and γ + are complex-valued and γ + corresponds to the idler field of γ − , generated by degenerate four-wavemixing due to the Kerr nonlinearity. Note that with this choice of Ansatz, we omit all higher order Fourier components to the total intracavity field, as in the operation regime of our device, they can be neglected to first order. Inserting drive and intracavity field Ansatz into the equation of motion yields γ + e i2Ω dp t + iK n d + |γ − | 2 + |γ + | 2 α d e iΩ dp t + γ − + γ + e i2Ω dp t + iKα d e iΩ dp t γ * − + γ * + e −i2Ω dp t α d e iΩ dp t + γ − + γ + e i2Ω dp t + iKα d e −iΩ dp t γ − + γ + e i2Ω dp t α d e iΩ dp t + γ − + γ + e i2Ω dp t + iK γ − γ * + e −i2Ω dp t + γ * − γ + e i2Ω dp t α d e iΩ dp t + γ − + γ + e i2Ω dp t With n d + |γ − | 2 + |γ + | 2 ≈ n d and after omitting all terms not linear in the small quantities γ − , γ + , we obtain where ∆ d = ω d − ω 0 and ∆ p = ω p − ω 0 describe the detunings of the two input field frequencies from the undriven cavity resonance. Sorting for frequency contributions leaves us with three equations The first of these three equations is the steady-state equation for the single-tone case and can be solved for n d using the approach presented above. To solve the other two coupled equations, we use ∆ p = ∆ d − Ω dp and define the susceptibilities and and get The driven Kerr-modes and their response To find the resonance frequency of the quasi-mode with susceptibility χ g , we solve for the complex frequencyω = ω p , for which χ −1 g = 0. Therefore, the condition is which is solved byω where the real part corresponds to the resonance frequency and the imaginary part corresponds to half the mode linewidth. So, as a consequence of the presence of the strong drive, the system has two resonances, which are split symmetrically with respect to the drive frequency ω d . The two Kerr-modes correspond to the cases when the cavity field γ − is resonant or when its idler field γ + is resonant and they have been discussed also in the context of nonlinear optical cavities S7,S8 and mechanical oscillators S9 . For the experiments described here, the argument of the square root will always be positive and hence, we get as resonance frequency of signal and idler mode and both modes are having a constant linewidth of κ. The most relevant regime for our experiment is given by ∆ d < 0 and ∆ d < Kn d . Then, the signal resonance is given by and the idler resonance by Hence, if S p is a probe tone scanning the driven SQUID cavity, we expect the response to be given by and to observe two resonances symmetrically positioned around the drive, that correspond to the two Kerr-modes of the driven Kerr cavity.

Kerr cavity two-tone response modeling
To observe and model the two-tone response of the driven Kerr cavity, we set a drive tone with fixed power P d and fixed frequency ω d to a point of non-zero flux bias and sweep the bare resonance frequency of the cavity ω 0 through the drive frequency by sweeping the bias flux Φ ⊥ . For each flux bias value during the sweep, we record a response trace S 21 of the SQUID cavity using the VNA. In Supplementary Fig. S10 we show in comparison the result of such a measurement with and without parametric drive. When ω 0 ∼ ω d , the flux-dependence of the cavity is modified considerably with respect to the case without a drive. The resonance frequency of the observed mode becomes nearly constant with flux and on the opposite side of the drive, a second resonance line appears. This second mode on the right side of the drive tone with net transmission gain indicates that we are entering the quasi-mode regime and have Josephson parametric amplification in both Kerr-modes. Main paper Fig. 2 shows and discusses a similar dataset in more detail including linescans and also analyzing the linewidths and the Kerr-mode resonance frequency and comparing the experimentally obtained values with theoretical calculations. The overall behaviour of cavity and Kerr-modes is in excellent agreement with the two-tone model of a Kerr oscillator as demonstrated by the high-level agreement between theory lines and experimental data for the resonance frequencies. The only subtlety we have to consider additionally is an effective linewidth broadening due to low-frequency flux noise out of the two-mode regime. From the data, it is obvious that the lineshape is not just broadened but distorted also on timescales comparable to the measurement time. Although due to this noise modulation the cavity response is not anymore described by its ideal response with an increased linewidth S11 , we fit it using Eq. (S5). The apparent κ we obtain from this procedure is however still a good measure for the effective linewidth. These flux-noise broadened linewidths are considerably larger than the intrinsic energy decay rate of the cavity and to take the noise-broadening and the two-level systems simultaneously into account we model the noise-free linewidth using GHz. The dashed line shows the theoretical resonance frequency without drive. When ω d ≈ ω0, the response of the cavity deviates significantly from the undriven response. Two modes are visible, one on the left side of the drive with increased absorption and reduced linewidth compared to the undriven case, and a second on the right side of the drive as a peak. The increased depth and reduced linewidth of the signal mode absorption dip reflects both, reduced impact of flux noise on the cavity line and Josephson parametric amplification. The Josephson amplification is also apparent in the idler mode on the right side of the drive, which shows net gain of the input signal. Due to the underlying four-wave mixing process, the resonance frequencies of the two Kerr-modes with respect to the parametric drive are equal in magnitude and opposite in sign.
where the first term κ 0 is the bare power and flux-independent decay rate, and the second term describes the generalized two-level system impact for detuning between drive and cavity/probe with the two-level system dephasing rate Γ 2 S10 . Finally, we phenomenologically take into account a broadening of the linewidth using where F is the flux responsivity and σ Φ is the rms value of the flux fluctuations through the SQUID. This relation seems to resemble closely what has been found numerically for tunable and frequency-fluctuating cavities S11 . The line in main paper Fig. 2 is a modeling of the experimentally obtained linewidths with Eqs. (S44) and (S45). The fit parameters are κ 0 = 2π · 200 kHz, κ 1 = 2π · 145 kHz, n crit = 0.26, Γ 2 = 2π · 300 kHz and σ Φ = 5.1 · 10 −3 Φ 0 . Note that the values used and obtained here are also in good agreement with the numbers we extracted from the modeling of the nonlinear single-tone response.
For the flux operation point II, we can use almost exactly the same parameters, except for a slightly increased κ 0 = 2π · 210 kHz.

Classical equations of motion
To model the optomechanical system with a Kerr nonlinearity, we use the classical equations of motion (EOM) with the mechanical displacement x, the mechanical resonance frequency Ω m , the mechanical damping rate Γ m , and the cavity pull parameter These equations are identical to the EOMs of linear classical optomechanics S12 , except for the additional Kerr term K|α| 2 in the equation for the intracavity field.

Single-drive solution
For a single cavity drive field we make the Ansatz and look for the steady-state solutionẍ =ẋ = 0. For the equilibrium offset displacement x 0 , we obtain where we used For the intracavity field amplitude α d , we find with the modified Kerr anharmonicity As for our device K > 2π · 10 4 Hz and Ωm < 2π · 10 Hz, we can assume in good approximation K ≈ K. From here it is straightforward to calculate the intracavity drive photon number n d and the phase φ d using the third order polynomial as for the bare Kerr cavity.

Single-drive Kerr backaction
If we allow also for fluctuations of the displacement and the intracavity field, we get the Ansatz and the equation of motion for the mechanical oscillator becomes For the intracavity field, we find For the linearization, we omit now all terms not linear in the small quantities δα and δx, we apply K ≈ K and remove the steady state solution. The remaining equations are which can be Fourier-transformed to Using the convention δα = δα * (−Ω), the observation that δx(Ω) = δx * (−Ω) and the definitions we write these equations as and solve for Note that here we used the earlier definition of the two-tone Kerr susceptibility We calculate the cavity response and the dynamical Kerr backaction for a sideband-unresolved optomechanical system with parameters close to the device discussed in Ref. S13 . In a, we plot the magnitude of the transmission matrix element S21 at a side-coupled cavity with a resonance frequency ω0 = 2π·8.167 GHz, a total linewidth κ = 2π·2.8 MHz, an external linewidth κe = 2π·1.4 MHz and an anharmonicity K = −2π·2.5 kHz for varying drive powers. We chose the drive powers such that the characteristic Duffing-like tilting of the resonance line to the left is clearly visible and keep the highest drive power below bifurcation. Using the obtained intracavity photon numbers n d , a mechanical resonance frequency Ωm = 2π · 274.4 kHz, an optomechanical single-photon coupling rate g0 = 2π · 57 Hz and our model presented in the text, we subsequently calculate the optical damping Γopt and the optical spring δΩm induced by the drive in a Kerr cavity. The result is plotted in panels b and c, respectively, and seems to agree well with the experimental results reported in S13 . It is furthermore interesting to compare the obtained dynamical Kerr backaction with the dynamical backaction for a completely identical system but without anharmonicity. The corresponding calculations for the linear system K = 0 are shown as insets. The most striking and exciting difference is the eightfold enhancement of the optical damping on the "red" side of the Kerr cavity. This enhancement can also be understood by the increased slope of the intracavity field S14 , cf. a, and the subsequently enhanced asymmetry for cavity photon up-scattering and down-scattering compared to the linear case. Different lines in each panel correspond to different drive powers and the detuning ∆ labelling the x-axes corresponds to the detuning from the zero-power resonance frequency ω0.
Inserting everything into the equation of motion for the mechanical oscillator, we obtain for a (real-valued) external driving force F ex (Ω) with the effective mechanical Kerr susceptibility describes the dynamical backaction of a single-tone driven Kerr cavity to the mechanical oscillator with the multiphoton coupling rate g α = α d g 0 = √ n d g 0 . We note, that the expression for Σ k is formally equivalent to the dynamical backaction in a linear cavity. The first term in square brackets describes the quasi-mode susceptibility χ g for the blue motional sideband field (cooling), while the second term -its conjugate at the opposite side of the drive tone χ g -is responsible for the red motional sideband field (amplification). The additional factors in parentheses 1 − iKn d χ p and (1 + iKn d χ p ) take into account that the blue and red motional sidebands are interfering with each other due to the Kerr-drive induced four-wave mixing. The blue motional sideband coincides with the idler field of the red sideband and vice versa and their interference will contribute and modify the simple picture of dynamical backaction in a linear cavity.
For a high quality factor mechanical resonator and in the weak-coupling regime, we can approximate Ω ≈ Ω m and get Note that χ eff m and χ eff 0 do not have the same dimension and 2iΩ m χ eff m ≈ χ eff 0 for Ω ≈ Ω m . Nevertheless we will call both a susceptibility for simplicity. The optical spring δΩ m and optical damping Γ opt are then given by Very recently, there has been an experimental report on dynamical backaction with a SQUID Kerr cavity in the sideband-unresolved regime S13 and we demonstrate in Supplementary Fig. S11 that our expressions lead to very similar results to the ones reported in Ref. S13 using parameters comparable with the reported ones. Using our expressions for the Kerr backaction, we can also model with high accuracy the regimes for mechanical self-oscillation induced by the strong parametric cavity drive in the device presented here. To do this, we perform a similar experiment to the one discussed in main paper Fig. 2, but at bias flux operation point II. At this operation point, we have a larger g 0 ∼ 2π · 3.56 kHz as well as a larger Kerr nonlinearity K ≈ −2π · 55 kHz. In addition, we use a slightly higher drive power. The result of the probe tone transmission S 21 for a bias flux sweep of the cavity resonance through the parametric drive tone is shown and discussed in Supplementary Fig. S12. The red dashed line shows how the bare cavity resonance frequency would be moving with bias flux Φ ⊥ in absence of a drive with frequency ω d . We can discriminate between five different regimes in the displayed data set.
• Regime 1: For the highest bias flux values, the cavity follows the undriven behaviour and its resonance frequency increases with reduced flux. The cavity linewidth and -shape is significantly distrorted by low-frequency flux noise. Panel c shows a linecut in this regime with a very broad and noisy single cavity absorption dip.
• Regime 2: At approximately Φ ⊥ /Φ 0 ∼ 0.476, the drive is positioned close to the blue sideband of the cavity ω d ≈ ω 0 +Ω m and multiple resonance lines appear in S 21 , four absorption modes and four gain modes are visible. The frequency distance between two neighboring absorption modes or two neighboring gain modes is always the mechanical frequency ∆ω ≈ Ω m . The appearance of a multiple-modes response is characteristic for a cavity with a strongly oscillating resonance frequency. In our device, the behaviour in this regime is generated and explained by the optomechanical instability and mechanical self-oscillations induced by the parametric drive being at the same time a very strong optomechanical blue-sideband pump tone. In panel b, we plot the calculated optical Kerr damping based on our equations for the dynamical backaction and on the independently determined device parameters. It is clearly visible that regime 2 corresponds to negative optomechanical damping, which exceeds the intrinsic mechanical linewidth Γ m ≈ 2π · 10 Hz and therefore the mechanical oscillator is in the instability regime of self-oscillations. In panel d, a linecut of regime 2 is shown. The original cavity is labelled with "s" and its oscillation-induced replicas with "s1", "s2", and "s3". In addition, we observe 4 versions of the idler Kerr-mode as small peaks, where the original mode is labelled with "i" and its replicas with "i1", "i2", and "i3".
• Regime 3: For 0.462 Φ ⊥ /Φ 0 0.468 the observed resonances return to a single absorption dip on the left side of the drive and a single small gain mode on the right side of the pump, indicating that the cavity frequency is not oscillating anymore. A linecut in this regime is shown in panel e.
• Regime 4: For 0.453 Φ ⊥ /Φ 0 0.462 the negative optical damping once again exceeds the intrinsic mechanical linewidth and a second regime of instability is entered. A linecut in this regime is shown in panel f, where three signal-modes and three idler-modes are visible and labelled as in regime 2. As the frequency difference between signal and idler Kerr-mode in this regime is close to the mechanical frequency, each mode almost overlaps with one replica of the corresponding mirror-mode and they form dip-peak pairs.
• Regime 5: For Φ ⊥ /Φ 0 0.453 the cavity jumps to the low-amplitude oscillation branch and the impact of both, the parametric drive and the dynamical backaction on the cavity lineshape become negligible. The cavity continues to shift with flux just as the undriven cavity would do. A linecut in this regime is not explicitly shown.
We note again that for the calculation of the Kerr backaction and the instability regimes as shown in panel b, the only free parameter was the line attenuation for the drive tone, which was adjusted to −81.5 dB for a good agreement between theory and experiment, a value very close to what was obtained from the line calibration in Supplementary Note 3 for operation point II. The other parameters used here are ω 0 , K, F, Ω m , Γ m , κ e , κ 0 , κ 1 , n crit , Γ 2 , σ Φ and they were all obtained from independent measurements or taken from theoretical estimates in the case of g 0 .
FIG. S12. Observing and modeling Kerr-backaction induced mechanical self-oscillations. a Magnitude of the cavity response S21 vs SQUID flux bias around operation point II in presence of a strong drive with ω d = 2π · 5.1721 GHz. The red dashed line shows the expected resonance frequency in absence of a drive. In the response, five regimes can be discriminated, which are described in detail in the text. Black horizontal dashed lines show the boundaries of the different regimes, which are closely related to the dynamical Kerr backaction onto the mechanical resonator. b shows the calculated optical damping for the flux range shown in a. The two solid red lines are the result for g0 = 2π · 3.4 kHz and g0 = 2π · 3.7kHz, the red-shaded area captures the range in between these values. For two particular flux ranges (which correspond to two different ranges of detunings between the bare cavity and the parametric drive and therefore to different intracavity drive photon numbers) the total damping rate of the mechanical oscillator with a intrinsic damping rate of Γm ≈ 2π · 10 Hz becomes negative. These ranges are indicated by gray areas and the threshold Γtot < 0 is indicated by a vertical black line. Here, the mechanical oscillator will become unstable and undergoes self-oscillations. These self-oscillations induce oscillations of the SQUID cavity resonance frequency due to the optomechanical interaction, which will in turn lead to the observation of multiple replicas of the cavity and idler modes in regime 2 and regime 4. c-f show individual linecuts of the cavity response shown in a, one for each regime from 1 to 4, where regime 1 corresponds to the highest bias flux values. In the regimes 2 and 4, shown in panels d and f, respectively, these Kerr-mode replicas are visible and labelled with sX and iX with X = 1, 2, 3. The original signal and idler Kerr-modes are labelled with s and i.

Phonon population with Kerr backaction
To calculate the equilibrium phonon population in the mechanical oscillator with Kerr backaction, we will use the linearized equations of motion for the quantum fieldsâ,â † andb,b † representing the classical intracavity fluctuations δα, δα * and β, β * with δx = x zpf (β + β * ), respectively. For the input noise, we useζ m for the mechanical oscillator andξ i ,ξ e for the internal and external cavity baths, respectively. We denote the input noise operators of the cavity at different frequencies with subscripts "+" and "−". The equations of motion becomė or in Fourier space and with the equations for the creation operators toô Note that for the external bath of the cavity, we assumed a single port, although strictly speaking the temperature on the left side of the transmission line could be different from the right side. As first step, we shorten the expressions by usinĝ and decouple the equations forâ andâ † . We get We can go even more compactâ As next step, we substituteb and obtainâ with which we can finally expressb only as a function of the input noise operatorŝ where the effective mechanical susceptibiliy is given by .

(S99)
The last approximation is valid for a high-Q m mechanical oscillator. Using the common relations for the expectation values of the noise correlators, we can calculate the phonon power spectral density in the mechanical resonator under dynamical Kerr backaction as and the total phonon number via This integration is performed numerically. We also note that we assumed a constant effective cavity bath occupation n c = κ e n e + κ i n i κ (S104) in the relevant frequency range. In Supplementary Fig. S13 we discuss the phonon occupation that results from a parametrically driven Kerr cavity for the relevant operation points I and II. Related discussions on optomechanical cooling using a Kerr cavity can be found in Refs. S14-S16 . Due to the dynamical backaction induced by the parametric drive, the phonon occupation of the mechanical oscilltor will deviate from the steady-state equilibrium value. In Fig. 2 of the main paper, we show the dynamical backaction of the drive vs flux bias in the SQUID. In a, we show the corresponding effect on the phonon occupation, normalized to the (maximum) case without parametric drive n th m = 80. For Φ ⊥ /Φ0 ∼ 0.36, the drive is located on the blue sideband of the cavity and amplification/heating is observed. For Φ ⊥ /Φ0 ∼ 0.327, the drive is located on the red sideband and cooling of the mechanical mode is visible. In the regime of interest, where our experiments take place (gray area), the effective occuption is increased by about 10%. The dashed and solid blue lines correspond to g0 = 2π · 1.8 kHz and g0 = 2π · 2 kHz, respectively. All other parameters are taken from independent measurements. In b, we show the equivalent occupation at operation point II with g0 = 2π · 3.4 kHz and g0 = 2π · 3.7 kHz for the dashed and solid red lines, respectively. The striped areas correspond to the instability regime Γm + Γopt < 0 in the case g0 = 2π · 3.7 kHz. All parameters are identical to the ones used for the backaction calculation shown in Fig. S12b, except for the drive power which is 1 dB reduced as this is the regime for our cooling experiments. In the experimentally relevant flux regime, marked by solid gray shading, the phonon occupation is increased by about a factor of ∼ 2.5 due to dynamical backaction.

Three-tone linearization
Finally, we discuss the linearized equations of motion with three input fields, i.e., where S 0 (t) is a third, weak probe field. We choose as Ansatz with real-valued and time-independent α d , complex-valued and time-independent γ − , γ + and complex-valued and time-dependent δα.
For the mechanical oscillator, we get with this (S108) In our experiments presented in the main paper we choose both Ω dp and 2Ω dp to be very far-detuned from the mechanical resonance frequency Ω m and therefore we can neglect the pure driving force terms with ±Ω dp and ±2Ω dp . After omitting the steady-state solution we get Gα d m δα * e iΩ dp t + δαe −iΩ dp t + G m γ − δα * + γ * − δα + G m γ + δα * e i2Ω dp t + γ * + δαe −i2Ω dp t .

SUPPLEMENTARY NOTE 9: THREE-TONE DYNAMICAL KERR BACKACTION
To calculate the dynamical backaction induced by the doubly-driven Kerr cavity, we omit any probe drives S 0,j for now. Also, we only keep terms linear in γ − , γ + . Finally, we omit δx j for j = 0, 1, 2, as these will not contribute to first order to the dynamical backaction. Under these conditions, we can write down the remaining terms in the four next-iteration field components contained in Eq. (S122).
With the unique replacements we can write these equations shorter as Note that the indices on the B and D terms are not describing a frequency shift like in the δα j , χ g,j and A j terms. The notation with the overline also does not refer to negative frequencies or complex conjugation here. Instead, these rather indicate their unique definition given in Eqs. (S128-S135).
Next, we inject the field relations back into the original equation for δα 0 (without the probe input) and obtain We perform a final variable substitution now where the nondegenerate four-wave mixing terms are still coded in blue and we obtain Inserting this result into the equation of motion for the mechanical oscillator and omitting higher order displacement terms leads to the mechanical susceptibility in the weak-coupling and high-Q m regime with the four-wave backaction α χ g (Ω m + Ω dp )J α (Ω m + Ω dp ) − χ * g (−Ω m + Ω dp )J * α (−Ω m + Ω dp ) +|g + | 2 χ g (Ω m + 2Ω dp )J + (Ω m + 2Ω dp ) − χ * g (−Ω m + 2Ω dp )J * + (−Ω m + 2Ω dp ) (S146) This multi-tone dynamical Kerr backaction has a very similar shape as a standard linear multi-tone backaction expression for several pump tones whose frequency difference is far detuned from the mechanical frequency. The main difference, besides the modified susceptibility χ g , is the J -factors. These J -factors take into account that the intracavity field one mechanical resonance frequency detuned from each of the pump tones α d , γ − , γ + also has contributions from the other fields due to four-wave-mixing. Without any Kerr nonlinearity, all As would be zero and all J s would be 1. With exclusively degenerate FWM, only the black terms in Eqs. (S141)-(S143) would survive. These terms describe the interference between the red and blue sidebands of α d , which are the idler fields of each other, but they also describe the interference between the red (blue) sideband of the γ − field with the blue (red) sideband of the γ + field. Also these form two pairs of signal and idler fields. In addition to these terms, there are the (in the equations blue-colored) non-degenerate FWM contributions. These modify the total backaction significantly, as can be seen in the main paper Fig. 3 or in Supplementary Fig. S14. Their origin and impact can be understood in two different ways. The first way is to consider that the cavity resonance frequency is permanently oscillating with the frequency Ω dp due to the beating of the α d -field with the γ ± -fields and by taking into account the dependence of the cavity susceptibility on the total intracavity field intensity via the Kerr nonlinearity. In this scenario, when a sideband of one of the tones is generated by mechanical motion at ±Ω m , higher-order sidebands of the scattered field will be generated by the oscillating susceptibility of the modulated cavity. These higher-order sidebands are detuned by the cavity oscillation frequency Ω dp from the original field and they will fall on top of other first order mechanical sidebands at ±Ω m ± Ω dp . The second way to understand this effect is that there are four-photon processes occurring, which involve one photon at the γ − -frequency ω p or at the γ + -frequency 2ω d − ω p , one photon at ω d and two sideband photons at different mechanical sideband frequencies.
The result is that by these processes the intracavity fluctuation field at Ω also gets contributions from five other fluctuation frequencies, as can be clearly seen in the Eq. (114). Note that in the general picture, an infinite number of fields will contribute to the field at Ω, but in our experimental situation, we can restrict the Fourier components to the most dominant ones.

SUPPLEMENTARY NOTE 10: FOUR-WAVE OMIT
If we take into account input probe fields at the frequencies of the relevant field components, the relations become Keeping all these terms, we get for the intracavity field To calculate the cavity response around the probe frequency, we will only have to keep a single term of these probe fields later, the one proportional to S 0,0 . To express the total driving force to the mechanical oscillator though, we have to keep them all for now. The four-wave mixing will generate fields also at frequencies that beat with the α d and the γ + field and therefore drive the mechanical oscillator. Nevertheless, the equations can be simplified according to the experimental situation.

Signal resonance red-sideband pumping
If the optomechanical pump field γ − is around the red sideband of the signal resonance and we probe around one mechanical frequency detuned from this pump, we have Ω ≈ Ω m . Due to the high quality factor of the mechanical oscillator, mechanical motion with Ω − Ω dp or Ω − 2Ω dp will be suppressed and we can neglect these terms in the equation for the field. As equation of motion for the mechanical oscillator under these conditions, we obtain which can also be written as Injecting this back into the equation for the intracavity field, we get δα(Ω) χ g (Ω) = i 1 − g − g * − χ g (Ω)P − (Ω) + g + χ * g (−Ω + 2Ω dp )A(Ω)P + (Ω) J − (Ω)χ eff 0 (Ω) κ e 2 S 0 (Ω) (S154)

Idler resonance blue-sideband pumping
If on the other hand the pump field γ − is located on the blue sideband of the idler resonance and we probe around one mechanical frequency away from the corresponding γ + field, we have Ω − 2Ω dp ≈ Ω m . In this case, mechanical motion with Ω and Ω − Ω dp will be irrelevant. Then, Injecting this back into the equation for the intracavity field, we get Optomechanical cavity response The response in both cases is given by • We start the experimental cycle with choosing the bias-flux operation point, either point I, and an in-plane magnetic field B . We ramp the in-plane current to its corresponding value.
• A parametric drive tone is sent to the cavity with fixed frequency ω d and power P d to match the chosen operation point.
• The cavity bias flux is adjusted manually to prepare the SQUID cavity in the driven Kerr-mode state.
• The frequency of the optomechanical pump is chosen to be either on the red sideband of the signal resonance or on the blue sideband of the idler resonance. The pump is activated with fixed frequency ω p and power P p .
The measurement • For the actual measurement, we start a python-based control and data acquisition script, which is programmed to wait for a terminal starting command before each data point.
• Prior to running the measurement, we input some fixed parameters to the script such as all values of the attenuators.
• We then manually adjust the probe VNA to a parameter set regarding frequency window, probe power and bandwidth in order to measure a clean OMIT response curve.
• Upon a terminal command, the script begins the acquisition and first catches all relevant information such as powers, frequencies, frequency spans, bandwidths as well as magnet DC currents from all relevant measurement equipment.
• The parameters obtained from the manually adjusted OMIT settings on the VNA are then re-used for all subsequent measurements. Based on the mechanical frequency and cavity frequency, an array of optomechanical pump frequencies is generated, which corresponds to an array of δ p . At the same time a corresponding set of VNA frequency ranges is generated to track the OMIT response for all the different pump frequencies.
• The script performs a narrow-band VNA scan to measure the OMIT response and stores the data in file 1, where all subsequent narrow-band scans for varying pump frequencies are attached as well. • The script performs a wide-band VNA scan to measure the cavity response and stores the data in file 2, where all subsequent wide-band scans for varying pump frequencies are attached as well.
• At this point the script will expect an input via the terminal, which tells whether we want to take the exact same measurement again for identical parameters or if we are going to proceed to the next pump detuning.
• After receiving our choice, the script sets the VNA to the cavity center frequency with a fixed span of 1 kHz and waits upon a terminal command for measuring the two VNA scans of the next point. During this waiting window, we have the possibility to counter possible bias flux drifts by manually adjusting the out-of-plane current, while permanently monitoring the cavity response at the response minimum.
• Both measures described in the latter two bullet points are critical to obtain a consistent set of data, as sometimes the bias flux and cavity starts to drift considerably on a slow timescale (∼seconds). This drift can significantly distort the captured OMIT response, which cannot be measured too fast due to the small mechanical linewidth. Another risk is that the cavity leaves the driven Kerr-state without the manual feedback control loop in between measurement points.
• After the cavity is stabilized and the measurement can proceed, the script repeats the cycle from gathering all relevant parameters from all machines to taking the two VNA traces and the waiting and stabilization time.

Data analysis
• Data analysis starts with a fit of the wide-band signal resonance response S 21 using Eq. (S5). From this fit, we obtain effective parameters for κ , κ e and ω s and a fit of the complex background.
• Using the background fit parameters, we calculate the complex background in the narrow-band frequency window of the corresponding OMIT response scan and divide it off the measured signal.
• We convert the frequency axis of the OMIT response to frequencies relative to γ − in the red-sideband case and relative to γ + in the blue-sideband case and fit the background corrected and frequency-shifted OMIT response using Eq. (S5) as well. As we are only interested in the resonance frequency and the effective linewidth of the OMIT resonance at this point, this procedure is as reliable but significantly simpler than using the full four-wave OMIT expression. As fit parameters we obtain Ω eff m and Γ eff .
• We substract the corresponding bare values Ω m and Γ m to obtain the dynamical backaction contributions δΩ m and Γ opt . We note that to make them match with theory for the red-sideband case and the blue-sideband case simultaneously with a single set of otherwise identical parameters, the bare values slightly differ between the red and blue case. For the mechanical resonance frequency, this blue-red-difference is about 8 Hz and for the linewidth it is about 1 Hz. As the bare mechanical frequency and linewidth depend strongly on temperature (cf. Supplementary Note 14), a small temperature difference of the mechanical oscillator during the two measurements might be the origin of these differences.
• For the range of theoretical values (shaded area between the two solid lines in main paper Fig. 3 and Supplementary Fig. S14), we consider uncertainties in the parameters going into the theoretical calculations. These include variations of the optomechanical single-photon coupling rate g 0 = 1.85 ± 0.05 kHz, of the driven cavity linewidth κ = 349 ± 20 kHz and of the bare cavity resonance frequency ω 0 = 5.2236 ± 0.1 MHz.

SUPPLEMENTARY NOTE 12: DATA FOR OMIT AND DYNAMICAL BACKACTION ON THE RED SIDEBAND OF THE SIGNAL RESONANCE
In this section, we present data on four-wave-OMIT and dynamical four-wave backaction for an optomechanical pump field γ − on the red sideband of the signal resonance with ω p = ω s − Ω m + δ p , where δ p is the detuning between the pump tone and the red sideband. The experimental setting is schematically shown in Supplementary Fig. S14a In this configuration, we follow the usual OMIT protocol in the experiment, i.e., we pump at a frequency ω p around the signal resonance red sideband and probe the cavity response S 21 in a narrow frequency window around ω ≈ ω p + Ω m . In addition, we measure the transmission in a wider range to capture the complete cavity absorption. We repeat this scheme for different detunings δ p , where the range of δ p is chosen to cover more than 2 signal resonance linewidths around the red sideband. One exemplary wide-band cavity transmission is displayed in Supplementary Fig. S14b, in c three corresponding narrow-band measurements are shown for three different δ p , clearly showing the characteristic OMIT window, representing the mechanical oscillator. From a fit to the OMIT response, shown as lines, we extract the mechanical oscillator resonance frequency and the effective mechanical linewidth. After subtracting the bare values, we obtain the dynamical backaction contributions δΩ m and Γ opt , which are plotted in d and e as circles.

SUPPLEMENTARY NOTE 13: MULTI-TONE KERR OPTOMECHANICS WITH NOISE INPUT
Working with quantum formalism for the equations of motion with noise input we obtain for the mechanical oscillator and for the intracavity fluctuation fieldŝ The latter two equations can be combined intô Just as for the classical equations, we need now the next iteration field componentŝ which lead to an expression forâ 0 given bŷ or, in its shortest version and the corresponding equation forâ † Signal resonance red sideband pumping and keep only first order terms to obtain We can find our earlier obtained four-wave dynamical backaction in this relation and writê For a pump on the red sideband of the signal resonance, a high mechanical quality factor and the detection frequency to be Ω ≈ Ω m , we can simplify the relations, i.e., keep only dominant terms and obtain a 0 = −ig − J − (Ω)χ g,0b0 + χ g,0M0+ (S174) (S175) a 1 = −ig α J α (Ω + Ω dp )χ g,1b0 + χ g,1M1+ (S176) a † 1 = ig α J * α (−Ω + Ω dp )χ g,1b0 + χ g,1M † 1− (S177) a 2 = −ig + J + (Ω + 2Ω dp )χ g,2b0 + χ g,2M2+ (S178) For the detection frequency range, we therefore get where we applied χ 0,0 ≈ 0 for Ω ≈ Ω m . Note that the cavity noise is well described for all frequencies in this approximation, but the upconverted mechanical noise is limited to one of the many sidebands. We can resolve and sort now for input noise frequency components, where we only keep cavity input noise terms around the signal and the idler resonances. The result iŝ For the cavity output field on one side of the feedline, we get This can be used just as in the red sideband case to calculate the output field power spectral density in units of phonons.
To calculate the corresponding phonon occupation, we use relations (S191) -(S196) and keep only cavity noise input which gives the mechanical power spectral density The integration of this relation over all frequencies then results in the effective phonon occupation in presence of the optomechanical coupling. in-plane magnetic field B . We ramp the in-plane magnet current to its corresponding value.
• A parametric drive tone is sent to the cavity with fixed frequency ω d and power P d to match the chosen operation point.
• The cavity flux bias is adjusted manually to prepare the SQUID cavity in the driven Kerr mode state.
• The frequency of the optomechanical pump is chosen to be either on the red sideband of the signal resonance or on the blue sideband of the idler resonance. The pump is activated with frequency ω p and power P p .
Thermal calibration measurement • The optomechanical pump is positioned on the red sideband of the signal resonance and the fridge is sitting at base temperature.
• For the actual measurement, we start a python-based control and data acquisition script.
• Prior to running the measurement, we input some fixed parameters to the script such as the values of all external attenuators in the input lines and the number of room-temperature amplifiers on the output line.
• We manually adjust the probe VNA to the parameter set regarding frequency window, probe power and bandwidth in order to measure a clean OMIT response. We hereby choose a red-sideband pump power, that is large enough to yield a clear OMIT response, but low enough to keep the effective cooperativity in the regime ∼ 1.
• Upon a terminal input command, the script begins data acquisition and first catches all relevant parameters such as powers, frequencies, frequency spans, bandwidths as well as magnet DC currents from all participating electronic devices.
• Afterwards, the script takes three datasets.
• First, it performs a scan of the narrow-band OMIT window using the VNA and the manually adjusted settings.
• Secondly, the center frequency and frequency span are sent to the spectrum analyzer to measure a power spectrum in the OMIT frequency range. During this spectrum analyzer measurement, the VNA frequency is set to a frequency 4 kHz detuned from the detection window of the spectrum analyzer to avoid any interference and the cavity response at a single frequency point is monitored permanently during the spectrum acquisition. This measure enables to control the cavity response by out-of-plane current feedback for the case the cavity is drifting during the spectrum measurement.
• Lastly, a wide-band VNA scan of the complete cavity response is taken.
• Each of the three data traces is stored in a separate file.
• For most temperatures, we repeat the measurement once.
• We adjust the fridge temperature to its new set value and after a temperature settling time of ∼10 minutes, we begin the cycle from the beginning at the new base temperature.
• Parametric drive and optomechanical pump powers P d and P p , respectively, were adjusted during the temperature sweep to keep driven cavity state and effective cooperativity nearly constant.

Four-wave-cooling measurement
• For the four-wave-cooling measurement we start a python-based control and data acquisition script, which is programmed to wait for an input terminal starting command before each data point.
• Prior to running the measurement, we input some fixed parameters to the script such as the values of all external attenuators in the input lines and the number of room-temperature amplifiers on the output lines.
• We set the parametric drive power and the optomechanical pump powers to the desired values P d and P p , respectively, and manually adjust the probe VNA to the parameter set regarding frequency window, probe power and bandwidth in order to measure a clean OMIT response.
• Upon a terminal input command, the script begins data acquisition and first catches all relevant parameters such as powers, frequencies, frequency spans, bandwidths as well as magnet DC currents from all participating electronic devices.
• At this point, the script waits for the final command to measure three traces. Once we observe a stable response in the OMIT window on the VNA screen, the script is continued.
• The first data trace acquired by the script is a narrow-band VNA scan of the OMIT response using the manually adjusted settings.
• Secondly, the center frequency and frequency span of the VNA trace are sent to the spectrum analyzer and an output power spectrum in the OMIT frequency range is acquired. During this spectrum measurement, the probe VNA frequency is set to a value several spectrum analyzer frequency spans detuned from the detection window of the spectrum analyzer window to avoid any interference between the VNA signal and the power spectrum. The VNA is continuously scanning a single frequency point of the cavity response to enable monitoring and feedback control of the bias flux for the case the cavity response is drifting due to flux drifts.
• Lastly, a wide-band VNA trace of the complete signal resonance S 21 is acquired.
• The three data traces are stored in individual files, where also the subsequent equivalent traces for the next settings are appended.
• The VNA setting are reset to the OMIT window and manual VNA control is enabled by the measurement script.
• At this point, we can choose by a terminal input between a repetition of the same measurement or a continuation to the next settings.
• In case of continuation to the next settings, the optomechanical pump power P p is manually set to its new value and the VNA probe settings are adjusted for a the next measurement iteration. The most important parameter is the frequency range for the OMIT and the spectrum analyzer window, which has to be increased with increasing dynamical backaction due to the considerable increase of the mechanical linewidth from the bare value of ∼ 15 Hz to the largest effective linewidth of ∼ 1.5 kHz for the highest achieved effective cooperativity.
• Depending on the optomechanical pump power, we also adjust the parametric drive power P d in some cases to keep the driven cavity response nearly constant. We suspect this measure is necessary as for large optomechanical pump powers, the experiment is at the edge of the linearized regime with respect to γ − , γ + and parametric drive depletion is occurring.
• Once the new parameters are adjusted, the data acquisition cycle starts from the beginning.
Thermal calibration data analysis • The data analysis of the thermal calibration begins with fitting the cavity response S 21 for each temperature using Eq. (S5). This fit provides us with a value for the linewidth κ and a fit function for the complex background.
• We divide off the complex background from both, the cavity response and the OMIT response in their corresponding frequency ranges. In addition we correct for a small phase rotation which is intrinsic to the Kerr cavity susceptibility and therefore not captured by the complex background.
• Next, we use the full Kerr-mode model function Eq. (S43) to fit the background-corrected cavity response once again. As fixed parameter for this second fit, we use κ from the first fit, the parametric drive power P d and the in-line attenuation. We also use the independently obtained κ e = 2π · 120 ± 20 kHz of the undriven cavity and allow for small variations, necessary to match the observed resonances. As fit parameters in this second fit, we get the detuning between drive and undriven cavity ∆ d and by using the Kerr polynomial (S22) we get the corresponding intracavity drive photon number n d .
• Next, we fit the background-corrected and rotated OMIT VNA response with Eq. (S5) and obtain Γ eff and Omega m as fit parameter.
• We fit the OMIT VNA response with the full four-wave OMIT model Eqs. (S154) and (S158). We input as fixed parameters K and g 0 from the flux arch, the parametric drive power P d , the optomechanical pump power P p and the in-line attenuation. Also, we use the fitted κ e from the previous full model cavity fit. As starting values for the remaining model parameters, we use ∆ d , n d and κ from the previous fits of the cavity. In the routine, we allow for a change of the bare cavity resonance frequency, i.e., of ∆ d due to possible flux fluctuations. Based on the characteristic polynomial Eq. (S22), we then dynamically adjust n d to the modified ∆ d . Additionally, we allow for changes in κ of up to ±100 kHz, with a lower limit of κ min = 2π · 320 kHz. As in the experiment between the VNA cavity scan and the VNA OMIT scan several minutes pass, during which the thermal noise spectrum is recorded, these allowed changes in fit parameters reflect possible drifts of the bare cavity resonance due to flux fluctuations in this time span. We calculate the intracavity γ − and γ + fields based on Eqs. (S35) and (S37). Ultimately, we obtain from this OMIT fit with the full model a value for the bare mechanical linewidth Γ m and the mechanical resonance frequency Ω m . As last fit parameter, we allow for a small correction of the OMIT resonance circle, which is necessary due to the uncertainty of the background extraction during the initial cavity fitting routine at the cavity resonance frequency. This is taken into account by allowing for a multiplication of the OMIT response by a small complex scaling factor (1 + x) e iβ with x, β 1.
• Finally, we fit the measured thermal noise spectrum using the full model Eqs. (S184) and (S185). We input as fixed parameters here all relevant quantities as obtained from the previous full model OMIT fit. The only remaining fit parameters at this point are the total detection output gain, converting the PSD in numbers of quanta to an absolute power and the uncalibrated equilibrium occupation of the mechanical oscillator n m . The last used parameter n add is adjusted to match the temperature dependence of the uncalibrated n m in the linear regime to the Bose distribution, which corresponds to a calibration of n m to n th m . As a result we obtain a number for the added photons n add ≈ 14 and the residual thermal phonon occupation shown in main paper Fig. 4. In Supplementary Note 15, we present some additional data on the temperature dependence of Ω m and Γ m as obtained from this procedure.
• Note: Due to slow SQUID cavity resonance frequency fluctuations and drifts, the measurement time of the spectrum acquisition was limited and we took these data with a reasonable compromise between number of data points, frequency span and bandwidth. For the lowest cooperativities however, when the mechanical linewdith is close to the intrinsic linewidth, the resolution bandwidth of the spectrum analyzer and the mechanical linewidth are comparable in size. In this case the effect of the resolution bandwidth is to smoothen and broaden the real mechanical power spectral density Lorentzian. To consider this effect in the fit curve for the power spectral density, we apply a moving-average-filter with the corresponding width in frequency space to the theory curve within the fitting routine itself. In the thermal calibration experiment, the resolution bandwidth of the spectrum analyzer was set to 5 Hz and therefore considerably smaller than the effective mechanical linewidths, cf. Supplementary Fig. S15.
• Error bars in thermal occupation: During the fitting procedure we simultaneously calculate the corresponding error bars for each point. These translate the impact of deviations of the cavity from its operation point. From this we estimate the difference in thermal occupation by fitting the thermal noise spectrum with the cavity parameters κ , ∆ d and n d . This difference is plotted as the error in y-axis of the inset of in Fig. 4c of the main paper.
Four-wave cooling data analysis • The data analysis of the four-wave cooling experiment begins with fitting the cavity response S 21 for each power setting of P d and P p using Eq. (S5). This fit provides us with a value for the linewidth κ and a fit function for the complex background.
• We divide off the complex background from both, the cavity response and the OMIT response in their corresponding frequency ranges. In addition, we correct for a small phase rotation which is intrinsic to the driven Kerr cavity susceptibility, and therefore not captured by the complex background.
• Next, we use the full Kerr-mode model function Eq. (S43) to fit the cavity response once again. As fixed parameter for this second fit, we use κ from the first fit, the parametric drive power P d and the in-line attenuation. We also allow for small variations of the bare external decay rate of each operation point κ I,II e and allow for κ e = 2π · (κ I,II e ± 10) kHz, which is necessary to match the observed resonances for all powers. As fit parameters in this second fit, we get the detuning between parametric drive and undriven cavity ∆ d . Additionally, by using the Kerr polynomial Eq. (S22) within the fit routine we get the corresponding intracavity drive photon number n d .
• We fit the background-corrected and rotated OMIT VNA response with Eq. (S5) and obtain Γ eff and Ω m as fit parameters.
• We fit the OMIT VNA response with the full four-wave OMIT model Eqs. (S154) and (S158) for red-signalsideband pumping or Eqs. (S157) and (S158) for the blue-idler-sideband case, respectively. We input as fixed parameters K and g 0 as determined from their flux dependence, the parametric drive power P d , the optomechanical pump power P p and the input line attenuation. Also, we use the κ e as determined from the full model cavity fit. As starting values for the remaining model parameters, we use ∆ d , n d and κ from the previous fits of the cavity. In the routine, we allow for adjustments of the bare cavity resonance frequency, i.e., of ∆ d up to ±2 MHz due to possible flux fluctuations, cf. operation range in main paper Fig. 2. Based on the characteristic polynomial Eq. (S22), the intracavity drive photon number n d is adjusted correspondingly within the fit routine. Additionally, we allow for adjustments of the total linewidth κ of up to ±100 kHz, but with a lower limit κ min,I = 2π · 320 kHz at operation point I and κ min,II = 2π · 350 kHz at point II. Several minutes pass in = 100 mK. The measurement scheme is detailed in Supplementary Note 14. Points show data, lines and shaded areas show fits to the full four-wave model, where all system parameters have been obtained from a VNA measurement of the cavity and the OMIT response, except for the thermal phonon occupation n th m (T b ) and the total output gain conversion, defining the absolute power scale of the spectra. Frequency-axis is given with respect to the red sideband optomechanical pump frequency ωp. b shows the resonance frequency shift ∆Ωm = Ωm(T ) − Ωm(15 mK) of the mechanical oscillator with dilution refrigerator base temperature. Panels c and d show the effective and intrinsic mechanical linewidth Γ eff (T ) and Γm(T ) vs base temperature T b . The effective linewidth is broadened by dynamical backaction and obtained from a fit to the OMIT response, the intrinsic linewidth is obtained from a fit to the OMIT response using the full four-wave model. The parameters obtained from this temperature dependence are used to calculate the residual thermal phonon occupation vs fridge temperature, the result is shown in main paper Fig. 4, where the values for each temperature have been averaged. Thermal calibration measurements were done at operation point I with an in-plane field of B = 21 mT. Square points in b-d indicate the results for the two datasets shown in a.
the experiment between the VNA scan of the OMIT response and the VNA scan of the signal resonance, during which the thermal noise spectrum is recorded. We allow for adjustments of some of fit parameters between the two scans, which reflects possible drifts of the bare cavity resonance due to flux drifts in this time span. We note that in principle also K and g 0 might experience small drifts due to the effective change in bias flux. As the flux-drift related variations of these two parameters are small within the operation range however, cf. main paper Fig. 2 and Supplementary Note 4, we work with constant average values for the analysis here. Based on Eqs. (S35) and (S37), we calculate also the intracavity fields γ − and γ + . As last fit parameter, we allow for a small correction of the OMIT resonance circle, which is necessary due to the uncertainty of the background extraction during the initial cavity fitting routine at the cavity resonance frequency. This is taken into account by allowing for a multiplication of the OMIT response by a small complex scaling factor (1 + x) e iβ with x, β 1.
• Finally, we fit the measured thermal noise spectrum using the full model Eqs. (S184) and (S185) for the redsignal-sideband case and Eqs. (S201) and (S185) for the blue-idler-sideband case, respectively. Once again, in this fitting procedure we allow for fluctuations of the cavity decay rate δκ = ±2π · 60 kHz and δ∆ d = ±2π · 0.8 MHz which are limited to small deviations around the average values obtained from the cavity and OMIT fitting routines. Note that the OMIT measurement and the cavity scan were taken prior and posterior (respectively) to the thermal noise detection and therefore we account for possible deviations of the cavity state due to the fluctuations in the system. Finally, the only remaining fit parameters are the total detection output gain and the equilibrium occupation of the mechanical oscillator n th m . Note that this last fit parameter was allowed to vary between 70-90 phonons for blue sideband driving and 60-90 phonons for red-sideband driving. Without these restriction, the fit often fails and the corresponding thermal phonon numbers are oscillating unsystematically between 50 and 130 phonons. The corresponding values without restrictions are considered in the error bars though, see below. The number of added photons n add ≈ 14 was determined via the thermal calibration procedure.
• Note: Due to slow SQUID cavity resonance frequency fluctuations and drifts, the measurement time of the spectrum acquisition was limited and we took these data with a reasonable compromise between number of data points, frequency span and bandwidth. For the lowest cooperativities however, when the mechanical linewdith is close to the intrinsic linewidth, the resolution bandwidth of the spectrum analyzer and the mechanical linewidth are comparable in size. In this case the effect of the resolution bandwidth is to smoothen and broaden the mechanical power spectral density Lorentzian. To consider this effect in the fit curve for the power spectral density, we apply a moving-average-filter with the corresponding width in frequency space to the theory curve within the fitting routine itself. An additional broadening effect of the spectrum might arise due to slow mechanical frequency fluctuations induced by a variation of the optical spring during bias flux drifts, cf. main paper Fig. 2.
• Based on the full ensemble of system parameters obtained by this multi-step analysis and fit procedure, we finally infer the resulting cooled phonon number n m of the mechanical oscillator by integrating Eq. (S187) for the red-signal-sideband case and Eq. (S203) for the blue-idler-sideband case, respectively. The results are plotted in main paper Fig. 4 and in Supplementary Fig. S16.
• Error bars for cooled number of phonons: During the fitting procedure we simultaneously calculate the corresponding error bars of each point. These translate the impact of deviations of the cavity from its operation point and of a different thermal occupation in the extraction of cooled photon number. For this we estimate the difference in the cooled phonons by fitting the thermal noise spectrum with the cavity parameters, which were extracted from the OMIT fit, in this case without any restriction to the thermal phonon number. This difference is plotted as the error in the y-axis of Fig. 4 and in Supplementary Fig. S16. Furthermore, we calculate the error in the extraction of the effective mechanical linewidth by computing the difference of Γ eff obtained from the full model noise fit and the one obtained from the fit of the OMIT VNA response with Eq. (S5). In addition, we consider an uncertainty in Γ m of ±1 Hz The sum of these errors is plotted as the error bar in the Γ eff /Γ m direction.

SUPPLEMENTARY NOTE 15: THERMAL CALIBRATION OF THE RESIDUAL MECHANICAL PHONON OCCUPATION
We perform the thermal calibration measurement and data analysis as described in Supplementary Note 14. In Supplementary Fig. S15 we present some of the results obtained from this experiment. In particular, we show the detected thermal noise spectrum in units of quanta for the lowest and highest temperatures T min b = 15 mK and T max b = 100 mK, including the fit from the full model, and we show the obtained temperature dependence of the mechanical parameters Ω m , Γ eff and Γ m . From the data it is clear, that both, the mechanical frequency and the intrinsic mechanical linewidth increase significantly with temperature and have a strong dependence on temperature at the lowest temperatures. Either intrinsic or pump-power-induced small variations of the chip temperature might therefore be a source for the observation that for good agreement of our datasets with the theory, we have to consider variations between the datasets of some 10 Hz in the mechanical resonance frequency and of a few Hz for the mechanical linewidth in the range 10 Hz< Γm 2π < 15 Hz.

SUPPLEMENTARY NOTE 16: FOUR-WAVE-COOLING WITH A PUMP ON THE RED SIDEBAND OF THE SIGNAL RESONANCE
We perform a four-wave cooling experiment with an optomechanical pump positioned on the red sideband of the signal resonance, cf. Supplementary Fig. S16. Data acquisition and data analysis are described in Supplementary Note 14.
FIG. S16. Red-sideband four-wave cooling of the mechanical oscillator. a Schematic representation of the experiment. A parametric drive is used to activate the Kerr quasi-mode state and an optomechanical pump is sent to the red sideband of the signal resonance ωp ≈ ωs − Ωm. From S21 measurements of signal resonance and OMIT, as well as a signal mode output power spectrum, we determine the cooled mechanical phonon occupation. The result for increasing sideband pump power Pp and effective mechanical linewidth, respectively, is shown in panel b. Inset shows that Γ eff ∝ Pp. Circles and stars are data, the lowest achieved occupation is nm ∼ 4, limited by cavity bifurcation instability. Dashed lines and shaded areas display the theoretical value range of nm, taking into account 60 < n th m < 90. Data were taken at operation point I and at B = 21 mT. In addition, the error bars represent possible deviations in the extraction of the plotted values based on their difference from the ones calculated based on the parameters extracted from the corresponding OMIT fit. For more details, see Supplementary Note 14. Panel c shows selected power spectra data for the points in b, which are plotted as stars.