Strong nanomechanical Duffing nonlinearity and interactions induced through cavity optomechanics

Abstract

Nonlinearity is a key resource in both classical and quantum signal processing. Nonlinear nanomechanical elements have found applications ranging from sensing to computing, while networks of nonlinear resonators, as well as nonlinearly coupled networks of linear resonators, constitute promising platforms for simulating complex dynamics. Here, we experimentally demonstrate an approach to realizing strong mechanical nonlinearity in nanomechanical resonators, fully controlled through optical laser drives. The mechanism exploits the nonlinearity of the radiation-pressure interaction in a cavity optomechanical system, which gives rise to a nonlinear optical spring effect. The resulting Duffing nonlinearity is conveniently tunable in strength via pump laser power, while its sign is controlled by laser detuning. Moreover, we demonstrate that the nonlinear optical spring mediates effective interactions between mechanical modes coupled to a common cavity, inducing tunable nonlinear interactions between them that impact spectral response and dynamics. These results establish cavity optomechanics as a versatile and in-situ reconfigurable platform for engineering nonlinear dynamics in resonators and networks.

Nonlinearity is a pervasive and often decisive feature of nanomechanical resonators [26]. Duffing-type dynamics give rise to phenomena such as frequency shifts, bistability, hysteresis, bifurcations, limit cycles and chaos [6, 15, 12, 20]. It is also widely explored for applications in sensing, signal processing, and mechanical computing approaches [29, 30, 38]. More broadly, exciting opportunities arise from controlled nonlinearity in both single resonators and networks or arrays: it can manipulate noise and fluctuations to yield strong squeezing [13] and thermal engines [42], produce solitons and frequency combs [18, 37], and enable novel computing paradigms including neuromorphic approaches and reservoir computing [44, 45, 31]. In the quantum domain, nonlinearity enables quantum information processing and leads to intriguing correlated phases of matter [10, 28]. But also in classical arrays with nonlinear interactions fascinating phases of matter and collective dynamics have been predicted and observed, including nonlinear and dynamical topological phases [24, 34, 35, 36, 27, 19] and quantized nonlinear Thouless pumping [17, 16]. Specifically, networks with nonlinear interactions have attracted interest, for example for enabling topological solitons that could occur at self-induced travelling domain walls [7, 9, 8].

Nanomechanics provides an excellent testbed for such dynamics, as well as a promising technological application platform. However, even if nanomechanical nonlinearities of geometrical origin are ubiquitous, they are generally weak, requiring significant amplitudes on the order of the dimensions of a nanostructure (e.g. the thickness of a flexural beam). Moreover, as they are rooted in material properties and device geometry, they exhibit very limited in-situ tunability.

Here, we experimentally demonstrate a nanomechanical resonator with an optically tuneble, strong Duffing nonlinearity. Rather than from device geometry or material properties, this nonlinearity originates from cavity optomechanical backaction [2]: it is linked to the intrinsic nonlinearity of the cavity optomechanical interaction and fully induced and controlled through laser light. Previously, light-controlled Duffing nonlinearity was observed in dissipatively coupled optomechanical systems [25, 47] for amplitudes of $\sim$100 nm and also predicted in coupled waveguides [1]. We show that nanomechanical nonlinearity can generally emerge from a standard dispersively coupled driven cavity optomechanical resonator, and can be orders of magnitude larger than intrinsic mechanical nonlinearity. In particular, the nonlinear shape of the cavity response as a function of displacement leads to a nonlinear optical spring effect, which results in tunable nonlinearities of different orders. We test this principle in an on-chip optomechanical device with exceptionally strong optomechanical coupling [22], so that motion explores a large fraction of the cavity response. We observe mechanical bistability for picometer-level vibrations, due to a Duffing nonlinearity whose sign and magnitude is completely controlled through an optical drive. Moreover, we demonstrate a scheme for realizing a pure inter-mode (cross-Kerr) nonlinearity between coupled resonators.

Figure 1: Laser-controlled Duffing nonlinearity. (a) Canonical optomechanical system: a cavity formed by two mirrors, where the movable right mirror experiences mechanical restoring force $F_{\text{m}}$ and radiation pressure $F_{\text{opt}}$. Mirror displacement $x$ linearly shifts the detuning $\Delta$ between cavity and pump laser. We assume the fast-cavity limit $\kappa\gg\Omega_{\text{m}}$, with $\kappa$ the cavity linewidth and $\Omega_{\text{m}}$ the intrinsic mechanical frequency. (b) Cavity photon number $n_{\text{c}}$ depends nonlinearly on $\Delta$ (black), resulting in optical backaction $F_{\text{opt}}$ that is nonlinear in $x$. Tuning the nominal laser detuning $\Delta_{0}$ across the Lorentzian response controls this nonlinearity. At the indicated $\Delta_{0}=-\kappa/(2\sqrt{3})$, a linear approximation (blue) for small displacements $Gx\ll\kappa$ gives the optical spring effect; for larger displacements $Gx\sim\kappa$, a cubic approximation (red) additionally yields Duffing nonlinearity. (c) Scanning electron micrograph of sliced nanobeam device. Flexural modes of the suspended silicon beams couple dispersively to a broadband photonic crystal cavity. (d) Driven response of the mechanical resonance at $\Omega_{\text{m}}/2\pi=17.6~$MHz for positive detuning $\Delta_{+}=\kappa/(2\sqrt{3})$ (left) and negative detuning $\Delta_{-}=-\kappa/(2\sqrt{3})$ (right). Forward (blue) and backward (red) frequency sweeps are shown for increasing spring laser powers $P_{\text{in}}$ (offset for clarity), with theory overlaid (black). (e) Cubic coefficients $\beta$ for increasing $P_{\text{in}}$ and drive amplitude $A_{\text{max}}$. Estimates from the experimental response (circles) match well with Eq. (6) (black line) for both $\Delta_{+}$ (green, $\beta<0$) and $\Delta_{-}$ (orange, $\beta>0$). (f) Ring-down at $P_{\text{in}}=1.0$ mW for $\Delta_{+}$ (green, shift $>0$) and $\Delta_{-}$ (orange, shift $<0$), for varying initial amplitude $A_{0}$. The nonlinearity does not affect the resonator’s exponential amplitude decay (top), while it manifests in the evolving instantaneous frequency $\Omega(t)$ once the mechanical drive is switched off for $t>0$.

We consider a mechanical resonator with intrinsic resonance frequency $\Omega_{\text{m}}$ and mass $m$, whose displacement $x$ linearly shifts the resonance frequency $\omega_{c}=\omega_{0}-Gx$ of an optical cavity with optomechanical coupling strength $G$, nominal frequency $\omega_{0}$ and linewidth $\kappa\gg\Omega_{\text{m}}$ (Fig. 1a). The latter condition—known as the fast-cavity or unresolved-sideband limit—implies that the cavity photon number $n_{\text{c}}$ responds instantaneously on mechanical timescales, such that photon dynamics can be adiabatically eliminated. A pump laser with intensity $P_{\text{in}}$ and frequency $\omega_{\text{L}}$, nominally detuned from the cavity by $\Delta_{0}=\omega_{\text{L}}-\omega_{0}$, drives it to a mean photon number $n_{\text{c}}$. Under the nonlinear optomechanical interaction, the resonator displacement evolves as

$$m\ddot{x}=-kx-m\gamma\dot{x}+F_{\text{opt}},$$ (1)

where $k=\Omega_{\text{m}}^{2}/m$ is the mechanical spring constant, $\gamma$ its energy damping rate and $F_{\text{opt}}=\hbar Gn_{\text{c}}$ the radiation pressure force exerted by the intracavity field.

Importantly, $F_{\text{opt}}$ depends on the displacement $x$ as it shifts the instantaneous spring laser detuning $\Delta=\Delta_{0}+Gx$, which in turn modulates the intracavity photon number $n_{\text{c}}=n_{\text{max}}h(2\Delta/\kappa)$ through the Lorentzian cavity response

$$h(u)=\frac{1}{1+u^{2}}.$$ (2)

Here, we define the maximum photon number $n_{\text{max}}=4P_{\text{in}}/(\hbar\omega_{0}\kappa)$ attained on resonance ($\Delta=0$), and the relative detuning $u=2\Delta/\kappa$. For small displacements $Gx\ll\kappa$, Eq. (2) can be linearised around the average relative detuning $u_{0}=2\Delta_{0}/\kappa$, yielding an optical backaction force

$$F_{\text{opt}} =F_{\text{opt}}^{\text{DC}}+\frac{2\hbar G^{2}n_{\text{max}}h^{\prime}(u_{0})}{\kappa}x=F_{\text{opt}}^{\text{DC}}-k_{\text{opt}}x$$ (3)

that varies linearly with $x$. The static radiation pressure $F_{\text{opt}}^{\text{DC}}=\hbar Gn_{\text{max}}h(u_{0})$ only displaces the resonator’s equilibrium position $x_{0}$, which we account for by redefining $x\to x-x_{0}$. Meanwhile, the optical spring constant $k_{\text{opt}}$ supplements the intrinsic mechanical spring $k$, shifting the mechanical resonance frequency by

$$\delta\Omega=\frac{k_{\text{opt}}}{2m\Omega_{\text{m}}}=-\frac{2g_{0}^{2}n_{\text{max}}h^{\prime}(u_{0})}{\kappa}.$$ (4)

This well-known shift is conveniently expressed using the vacuum optomechanical coupling rate $g_{0}=Gx_{\text{zpf}}$ and the resonator’s zero-point quantum fluctuation amplitude $x_{\text{zpf}}=\sqrt{\hbar/(2m\Omega_{\text{m}})}$, even though Eq. (4) equally applies in the classical domain.

For larger amplitudes $Gx\sim\kappa$, the instantaneous detuning samples a significant fraction of the cavity response, so that the nonlinear shape of $h(u)$ needs to be taken into account. By moving the laser detuning $u_{0}$ across this Lorentzian, we tune the nonlinear terms in the optical force’s Taylor series. To induce a cubic nonlinearity to leading order, we select $u_{0}=u_{\pm}=\pm 1/\sqrt{3}$ (Fig. 1b) where $h^{\prime\prime}(u_{0})=0$. Consequently, $|h^{\prime}(u_{0})|$ is maximal, so that the linear spring shift $|\delta\Omega|$ is largest. The equation of motion

$$\ddot{z}+\Omega^{2}z+\gamma\dot{z}+\beta z^{3}+\mathcal{O}(z^{4})=f(t),$$ (5)

of the resonator’s normalized displacement $z=x/x_{\text{zpf}}$ is then the Duffing equation, featuring a cubic term with coefficient (see Supplementary Information)

$$\beta=-\frac{8\Omega_{\text{m}}n_{\text{max}}h^{\prime\prime\prime}(u_{\pm})g_{0}^{4}}{3\kappa^{3}}=-6\Omega_{\text{m}}\delta\Omega\frac{g_{0}^{2}}{\kappa^{2}}.$$ (6)

In Eq. (5), $\Omega=\Omega_{\text{m}}+\delta\Omega$ is the spring-shifted linear frequency of the mechanical resonator and $f(t)$ represents a driving force. We neglect the higher-order terms $z^{n}$ for $n>3$ in the optical force as their coefficients scale with $(g_{0}/\kappa)^{n}\ll 1$.

Crucially, this optomechanically induced Duffing coefficient $\beta$ is fully tunable in both magnitude and sign. By setting $u_{0}=-1/\sqrt{3}$, the linear spring shift $\delta\Omega<0$ is negative, so that the corresponding $\beta>0$ represents a hardening nonlinearity, while for $u_{0}=1/\sqrt{3}$ a softening nonlinearity $\beta<0$ is obtained. Moreover, the magnitude of $\beta$ scales linearly with $\delta\Omega$ and is, therefore, directly tuned by the pump laser power $P_{\text{in}}$.

In principle, this optically tuneable mechanical nonlinearity occurs in any fast-cavity optomechanical system. Whether it leads to observable effects depends on the vibrational amplitude. The critical phonon number $n_{\text{NL}}$ above which a coherently driven resonator’s response becomes bistable [41] occurs when the nonlinear mechanical frequency shift exceeds the linewidth $\gamma$. It can be expressed as

$$n_{\text{NL}}=\frac{8}{3\sqrt{3}}\frac{\Omega_{\text{m}}\gamma}{|\beta|}=\frac{4}{9\sqrt{3}}\frac{\gamma}{|\delta\Omega|}\frac{\kappa^{2}}{g_{0}^{2}}.$$ (7)

Owing to its large ratio $g_{0}/\kappa$ and fast-cavity operation, a well-suited optomechanical system to observe nonlinearity at low amplitude is the sliced nanobeam photonic crystal cavity [23], in which the effect of the nonlinear cavity response on the optical read-out of mechanical motion has been observed before [22] as well as strong spring shifts [33]. We employ this sliced nanobeam platform to demonstrate a controlled Duffing nonlinearity. Our on-chip device features a suspended silicon beam with a mechanical resonance at $\Omega_{\text{m}}/2\pi=17.6~$MHz (effective mass $m=1$ pg) coupled to a telecom optical mode (linewidth $\kappa/2\pi=313$ GHz) at $g_{0}/2\pi=3.5\pm 0.3~$MHz. For this high-quality resonator ($\gamma/2\pi=4.1$ kHz), we readily achieve spring shifts $|\delta\Omega|>30$ kHz when illuminating the cavity with a focused ‘spring’ laser of up to $\sim$1  mW incident from free space. Here we remark that the cavity coupling efficiency is approximately 3%, so the the coupled laser power is maximally a few tens of $\upmu$W.

In Fig. 1d, we show the resonator’s response to coherent driving $f(t)=f_{0}\cos(\omega_{\text{d}}t)$ at frequency $\omega_{\text{d}}$, for a range of incident spring laser powers $P_{\text{in}}$ and detunings $u_{0}$. The drive force $f(t)$, with fixed amplitude $f_{0}$, is applied optically by modulating an additional weak ‘force’ laser resonant with the cavity, while a far-detuned ‘detect’ laser reads out the resulting motion (see detailed methods description in Supplementary Information). When the spring laser is off, we observe a Lorentzian response that is approximately symmetric in frequency, indicating essentially linear resonator dynamics for this drive force strength, apart from a small residual shift arising from the detection laser and intrinsic nonlinearity (see below). Turning on the spring laser at detuning $\Delta_{0}=u_{0}\kappa/2=\kappa/(2\sqrt{3})$, where $\delta\Omega>0$, induces an observed softening nonlinearity that increases in strength with $P_{\text{in}}$. The driving force is kept the same, so the optomechanically induced nonlinearity clearly exceeds the intrinsic mechanical nonlinearity of the beam. Notably, at higher $P_{\text{in}}\geq 0.6$ mW, the resonator exhibits the bistable ‘shark fin’ response typical for a Duffing resonator, and associated hysteresis depending on the sweep direction of $\omega_{\text{d}}$. Conversely, for red-detuned $\Delta_{0}=u_{0}\kappa/2=-\kappa/(2\sqrt{3})$, where $\delta\Omega<0$, we find a hardening nonlinearity of similar strength but opposite sign. In that case, the characteristics of the response appear reflected in frequency. At both detunings, the experimental frequency response agrees well with Duffing theory (black curves in Fig. 1d, see Supplementary Information) for the coefficient $\beta$ predicted by Eq. (6).

Evidently, the nonlinear resonance shift $\sim\beta A^{2}/\Omega_{\text{m}}$ opposes the linear spring shift $\delta\Omega$. This can be understood intuitively: at $u_{0}=\pm 1/\sqrt{3}$ the derivative $|h^{\prime}(u_{0})|$ is maximal, so that any detuning excursions will reduce the average derivative of the cavity response and, thus, the effective spring shift. Ultimately, for extremely large amplitudes $GA\gg\kappa$, the cavity will be far off-resonant for most of the mechanical cycle, so that the resonator reverts back to its intrinsic resonance frequency (in the absence of other nonlinearities).

The cubic coefficient $\beta$ can be extracted from the experimental data, by inverting the relation $A_{\text{max}}=\sqrt{8(\omega_{\text{d,max}}-\Omega)\Omega/(3\beta)}$ between the optimal drive frequency $\omega_{\text{d,max}}$ and amplitude $A_{\text{max}}$ of the maximum response [41]. As shown in Fig. 1e, experimentally obtained estimates of $\beta$ match well with Eq. (6), showcasing the full tunability of our scheme. At zero pump power, a small residual nonlinearity $\beta<900$ Hz² is estimated from the observed shift in the response peak. This is related to the intrinsic nonlinearity as well as the optomechanical nonlinearity of the detuned detection laser, which remains on even for zero spring laser power. We estimate the intrinsic geometrical nonlinearity contribution at $\beta_{\text{intr}}\approx 300$ Hz² from finite element modelling (Supplementary Information) and the contribution from the detection laser at $\beta_{\text{det}}\approx-30$ Hz².

The largest observed optomechanical nonlinear coefficient $\beta\approx 16000$ Hz² is thus $50$ times larger than the strength of the intrinsic nonlinearity. Appreciable Duffing nonlinearity emerges at critical amplitude $A_{\text{NL}}=\sqrt{n_{\text{NA}}}x_{\text{zpf}}$ as low as $A_{\text{NL}}=17\,000\,x_{\text{zpf}}$, only $20$ times larger than the thermal amplitude of our room-temperature resonator. Moreover, with $x_{\text{zpf}}=22$ fm (estimated from simulations), $A_{\text{NL}}=0.36$ nm represents a lateral displacement of less than a nanometer, much smaller than the critical amplitudes observed through different schemes [25] or typical intrinsic mechanical nonlinearities for flexural nanomechanical resonators of similar dimensions.

We probe the time-domain dynamics of our nonlinear resonator in a ring-down experiment (Fig. 1f). The resonator is first excited to a high amplitude $A_{0}$ by a spring laser modulation, which is subsequently switched off at $t=0$. Initially, for $t<0$, the instantaneous frequency $\Omega(t)$ of the resonator’s motion (see Supplementary Information), shown in Fig. 1f, bottom, is locked to the drive. After the modulation is switched off, the resonator is free to evolve at its natural frequency, which for a Duffing resonator depends on its amplitude [15, 20]. Consequently, at $t=0$, the estimated instantaneous frequency $\Omega(t)$ jumps to a nonlinearly shifted value. As the resonator’s amplitude rings down (top), this nonlinear shift gradually reduces, until the resonator oscillates at its linear (spring-shifted) frequency $\Omega$.

Figure 2: Nonlinear shift in susceptibility to thermal fluctuations. (a) A strong mechanical drive induces large amplitude oscillations, shifting the resonator’s susceptibility $\chi_{\text{NL}}$ to additional (thermal) forces. (b) Thermomechanical spectra of the resonator, coherently driven by modulations of a weak, resonant ‘force’ laser at increasing depth $c_{\text{d}}$. The strong ($P_{\text{in}}=1.0$ mW) spring laser is detuned by $\Delta_{-}=-\kappa/(2\sqrt{3})$ (top) or $\Delta_{+}=\kappa/(2\sqrt{3})$ (bottom). Bright bands indicate thermal fluctations on top of the driven amplitude, while fainter bands result from nonlinear transduction. (c) Fitted center frequencies of the Lorentzian thermal contribution plotted against the square of the driven amplitude, for $\Delta_{-}$ (orange) and $\Delta_{+}$ (green) at two powers $P_{\text{in}}$. Expected susceptibility shifts from Eq. (8) (black) are overlaid in (b) and (c).

We continue by exploring the resonator’s thermal dynamics in the presence of nonlinearity. Figure 2 shows the measured thermomechanical spectrum of the resonator while it is simultaneously driven by a force laser modulation at increasing depth $c_{\text{d}}$. The coherent high-amplitude oscillations appear as a narrow peak with increasing amplitude, saturating the color scale. Notably, the broad spectrum of thermal vibrations now shifts in frequency, with the direction of the shift depending on the spring laser detuning: for $u_{0}=-1/\sqrt{3}$ it moves up in frequency as $c_{\text{d}}$ increases, while for $u_{0}=1/\sqrt{3}$ it moves down. We note that additionally a weak, frequency-reflected ‘ghost’ band appears that moves in the other direction. We attribute this to a combination of nonlinear transduction [22] that mixes the thermal fluctuations with the strong coherent oscillations, or a small magnitude of nonlinearly driven motion at such frequencies, as observed before in nonlinear electromechanical [13] and microwave [5] systems.

Figure 2 demonstrates that in the presence of optomechanical nonlinearity, coherent oscillation alters the resonator’s susceptibility $\chi_{\text{NL}}$ to small thermal forces. Using harmonic balance [21] (Supplementary Information), we find that $\chi_{\text{NL}}$ remains Lorentzian like its linear counterpart $\chi_{\text{m}}$ (Fig. 2a), but with a center frequency shifted by

$$\delta\Omega_{\text{NL}}=\frac{3}{4}\frac{\beta A^{2}}{\Omega_{\text{m}}}=-\frac{9}{2}\delta\Omega\frac{g_{0}^{2}}{\kappa^{2}}A^{2},$$ (8)

where $A$ is the normalized amplitude of the strong coherent oscillations. The thermal peak frequencies fitted from Fig. 2b, plotted in Fig. 2 against the measured coherent amplitude $A$, agree well with the predicted values from Eq. (8).

Figure 3: Optically-mediated nonlinear coupling between mechanical resonators. (a) Multi-mode optomechanical system: both mirrors are free to move, so that photon number $n_{c}$ depends nonlinearly on both displacements $x_{1},x_{2}$ via the Lorentzian cavity response. The resulting radiation pressure $F_{\text{opt}}$ on each mirror induces a nonlinear, amplitude-dependent optical spring $k_{\leftrightarrow}(x_{1},x_{2})$ coupling the two mechanical elements. (b, left) Forward-sweep frequency response of resonator $2$ ($\Omega_{\text{m},2}/2\pi=12.8$ MHz) for increasing drive depth $c_{\text{d}}$, while coupled to resonator $1$ ($\Omega_{\text{m},1}/2\pi=17.6$ MHz) via a cross-resonator optical spring induced by modulating a single spring laser (detuning $\Delta_{-}=-\kappa/(2\sqrt{3})$) at their difference frequency $\Delta\Omega_{12}$. (b, right) Similar, but with an additional static spring laser at detuning $\Delta_{+}=\kappa/(2\sqrt{3})$ that cancels the intra-resonator nonlinearity while preserving the cross-resonator nonlinearity, resulting in reduced peak splitting with increasing amplitude. (c) Response of uncoupled resonator $1$ with both spring lasers active. Cancellation of nonlinearity restores a Lorentzian lineshape and eliminates hysteresis between forward (blue) and backward (red) sweeps (offset for clarity). (d) Response of coupled resonator $2$ at fixed drive depth, for single-laser (left, $c_{\text{d}}=0.08$) and two-laser (right, $c_{\text{d}}=0.045$) driving. Hysteresis between forward (blue) and backward (red) sweeps is observed for single-laser driving, yet largely absent for two-laser driving.

The above results demonstrate how nonlinear optical backaction can control the nonlinear dynamics of just a single mechanical mode. Interestingly, rich multimode nonlinear dynamics can be induced for multiple resonators that are simultaneously coupled to a common cavity. In that case, the cross-resonator optical spring, produced by their mutual backaction, can generate effective nonlinear hopping interactions between the mechanical modes. This has the prospect of creating cross-Duffing interactions, analogous to cross-Kerr effect in optics and photonics. Such interactions are essential for a variety of fascinating phenomena, including nonlinear energy localization [40], frequency pulling in synchronization [32], correlated phases of matter [14] and nonlinear dynamical topological solitons [7, 9].

Here we analyse such coupled nonlinear dynamics. Two cavity-coupled resonators with displacements $x_{1},x_{2}$ (Fig. 3a) modulate the instantaneous cavity detuning $\Delta=\Delta_{0}+G_{1}x_{1}+G_{2}x_{2}$ through their optomechanical coupling rates $G_{j}$. Consequently, the optical force $F_{\text{opt},j}(x_{1},x_{2})$ acting on each resonator $j$ now depends on both their displacements. For small amplitudes, linearizing $F_{\text{opt},j}(x_{1},x_{2})$ yields a cross-resonator optical spring that dynamically couples $x_{1}\leftrightarrow x_{2}$ by a linear spring constant $k_{\leftrightarrow}$ (see Supplementary Information). However, for large displacements $G_{1}x_{1},G_{2}x_{2}\sim\kappa$, the nonlinear cavity response must again be taken into account, resulting in a displacement-dependent cross-resonator spring $k_{\leftrightarrow}(x_{1},x_{2})$.

The sliced photonic crystal nanobeam hosts several flexural mechanical overtones (Supplementary Information) that couple to the cavity, allowing us to study such nonlinear coupling. We focus on a pair of modes that includes the aforementioned $\Omega_{\text{m},1}/2\pi=17.6$ MHz mode and a lower-frequency mode at $\Omega_{\text{m},2}/2\pi=12.8$ MHz with comparable damping rate $\gamma_{2}/2\pi=3.7$ kHz and vacuum optomechanical coupling rate $g_{0,2}/2\pi=3.8\pm 0.1$ MHz. With the spring laser, both modes experience spring shifts up to $|\delta\Omega_{j}|\sim 30~$kHz so that they resonate at $\Omega_{j}=\Omega_{\text{m},j}+\delta\Omega_{j}$. Due to the large frequency separation $\Delta\Omega_{12}=\Omega_{1}-\Omega_{2}\gg\gamma_{j}$, the static cross-resonator optical spring $k_{\leftrightarrow}$ has virtually no effect on their dynamics. That frequency difference can, however, be overcome by modulating the intensity $P_{\text{in}}$ of the spring laser at $\Delta\Omega_{12}$. The resulting time-modulated spring $k_{\leftrightarrow}(t)$ stimulates frequency conversion between the modes, by producing sidebands in the optical force that are mutually resonant. For small amplitudes, this leads to an effective linear hopping interaction that resonantly couples both modes at a rate $J=c_{\text{h}}\sqrt{\delta\Omega_{1}\delta\Omega_{2}}/2$ (Supplementary Information), where $c_{\text{h}}$ is the depth of the spring laser modulation [33].

In Fig. 3b (left panel), we probe the coherent response of the lower-frequency resonator $2$ in the presence of a coupling modulation with depth $c_{\text{h}}=0.5$, imprinted on the red-detuned ($u_{0}=-1/\sqrt{3}$) spring laser. Here, resonator $2$ is driven around its resonance frequency $\Omega_{2}$ by an additional weak modulation of the spring laser at depth $c_{\text{d}}$. For small $c_{\text{d}}<0.02$, two Lorentzian peaks are observed split by frequency $2J$, signifying hybridisation of the two mechanical resonators by the effective optical coupling. For stronger drives, the larger amplitude induces nonlinearity both in the individual optical springs $k_{\text{opt},j}$, as well as the cross-resonator spring $k_{\leftrightarrow}$. The former increases the effective mechanical resonance frequencies at higher amplitudes, illustrated in Fig. 3b, top, by the overall upshift in peak frequencies, while the latter reduces their effective splitting, which is proportional to the (amplitude-dependent) coupling rate.

Conveniently, the cross-Duffing nonlinearity can be isolated by simultaneously applying two spring lasers at opposite detunings $u_{\pm}=\pm 1/\sqrt{3}$ to cancel the cubic nonlinearity of the individual resonators. As shown in Fig. 3c, a symmetric mechanical response is then restored that exhibits no hysteresis, even at high drive power. By again modulating only one of the spring lasers at $\Delta\Omega_{12}$, the nonlinearity in the cross-resonator spring $k_{\leftrightarrow}$ is, nevertheless, induced (Supplementary Information). As presented in Fig. 3b (right panel), the response of resonator 2 then remains centered around $\Omega_{2}$ with increasing $c_{\text{d}}$. For larger amplitudes, we do, however, see a strong reduction in the frequency splitting between the peaks, indicating a reduction in $k_{\leftrightarrow}$ by the nonlinearity. This is also reflected in the example drive frequency response curves shown in Fig. 3d: For a single modulated laser, these resemble two shark-fin-like curves for the two hybridized normal modes of the pair. With an isolated cross-Duffing nonlinearity, the response spectrum looks more symmetric, characterized by two peaks that tilt towards each other; shifting closer together as the amplitude increases.

Figure 4: Time-evolution of nonlinearly coupled resonators. (a) Ring-down of coupled resonators $1$ (top) and $2$ (bottom) under single-laser driving at detuning $\Delta_{-}=-\kappa/(2\sqrt{3})$. For $t<0$, resonator $1$ is coherently excited by modulating the spring laser at depth $c_{\text{d}}$. At $t=0$, the drive tone is switched off and a coupling tone at the difference frequency $\Delta\Omega_{12}$ is switched on, inducing Rabi oscillations whose period increases and visibility reduces at larger amplitudes. (b) Similar, but with an additional static spring laser at detuning $\Delta_{+}=\kappa/(2\sqrt{3})$. Cancellation of intra-resonator nonlinearity keeps the coupling tone resonant to restore Rabi oscillation visibility, while the retained cross-resonator nonlinearity isolates the amplitude-dependent period lengthening. (c) Instantaneous coupling rates $J_{\text{NL}}(t)$ extracted from the coupled ring-downs in (b). For small drive amplitudes $c_{\text{d}}$, estimated $J_{\text{NL}}(t)$ approach the linear coupling $J/2\pi\approx 7$ kHz (black) throughout. For larger drives, $J_{\text{NL}}$ is initially reduced, recovering toward $J$ as the amplitude dissipates. (d) Plotting the estimated $J_{\text{NL}}$ against the detuning variance induced by both resonators reveals a clear trend independent of initial drive amplitude, with good agreement with Eq. (9) (black) for small to moderate mechanical amplitudes.

Nonlinear coupling also impacts the resonator pair’s temporal dynamics. In the ringdown experiments shown in Fig. 4a for a single spring laser at $u_{0}=-1/\sqrt{3}$, resonator $1$ is resonantly driven until $t=0$, when the coupling tone at $\Delta\Omega_{12}$ is switched on. For small amplitudes, we observe Rabi oscillations in which energy is exchanged between the resonators at the Rabi frequency $J$ during ringdown. For larger driving amplitude $c_{\text{d}}$, the Rabi period increases — indicating a nonlinear reduction in the effective coupling rate $J_{\text{NL}}$ — while visibility reduces. The latter results from unequal intra-resonator nonlinearities, which detunes their effective frequency difference from the coupling tone at $\Delta\Omega_{12}$ causing incomplete Rabi oscillations.

With two spring lasers at $u_{\pm}$ (Fig. 4b), a pure cross-resonator nonlinearity is isolated. The coupling tone then remains resonant, so that full fringe visibility is recovered at high amplitudes. The effective Rabi rate $J_{\text{NL}}$ clearly increases as the resonator amplitudes ring down. Figure 4c quantifies this behavior by extracting the instantaneous value of $J_{\text{NL}}$ from the phase and amplitude of both resonators (Supplementary Information). We see that $J_{\text{NL}}$ starts out at a reduced effective coupling rate at $t=0$ and then gradually increases to $J/2\pi\approx 7$ kHz as the overall amplitude dissipates. Finally, Fig. 4d plots the extracted values of $J_{\text{NL}}$ against the summed squared amplitude of both resonators, expressed as the variance in the optomechanical detuning shift $\Delta_{\text{OM}}^{2}=g_{0,1}^{2}a_{1}^{2}+g_{0,2}^{2}a_{2}^{2}$, where $a_{j}=A_{j}/x_{\text{zpf},j}$ is the normalized amplitude. We see a clear trend independent of the initial amplitude $A_{0}$, which, for small $\Delta_{\text{OM}}^{2}$ agrees well with the theoretically predicted

$$J_{\text{NL}}=J\left[1-\frac{9}{4}\left(\frac{g_{0,1}^{2}}{\kappa^{2}}a_{1}^{2}+\frac{g_{0,2}^{2}}{\kappa^{2}}a_{2}^{2}\right)\right]$$ (9)

that we derive in the Supplementary Information. These results illustrate the capability of optomechanical backaction to precisely control nonlinear coupling rates in multimode mechanical systems.

In conclusion, we have experimentally demonstrated that cavity optomechanical backaction can serve as a powerful and highly tunable source of nanomechanical nonlinearity. In the fast-cavity regime, the nonlinear Lorentzian response of the optical cavity produces a nonlinear optical spring that gives rise to strong Duffing dynamics. Using a sliced photonic crystal nanobeam cavity, we used this mechanism to realize optically induced nonlinearities that exceed the intrinsic mechanical nonlinearity by more than an order of magnitude, yielding sub-nanometer level critical nonlinear amplitudes. We note that the mechanism is very general, and will occur in any cavity optomechanical resonator in the fast-cavity regime with standard dispersive optomechanical coupling. Even the current sliced nanobeam platform could be improved to yield still much stronger nonlinearity, as $g_{0}$ values have been observed that exceed that in the current demonstration by a factor 10, and optical linewidths $\kappa$ that are at least one order of magnitude smaller [22]. As the Duffing coefficient $\beta$ scales with $g_{0}^{4}/\kappa^{3}$, it is expected to be many orders of magnitude larger for optimized parameters. We do note that that increases higher-order nonlinearities by an even greater amount, resulting in an interesting trade-off depending on the specific target application.

A key feature of this nonlinearity is that it is fully controllable through the power and detuning of external laser drives. We showed that the sign and magnitude of both self-Duffing and cross-Duffing interactions can be individually tuned by using multiple lasers. Moreover, specific interactions can be addressed in the synthetic mode dimension through judicious temporal modulation. This points to the possibility of engineering arbitrary nonlinear Hamiltonians in networks of resonators through spectral control of drive fields at frequency scales of the order of optical linewidths and mechanical frequencies. Interactions can be dynamically reconfigured in situ fully optically, without requiring changes to device geometry or material composition. Indeed, while we have focused here on the third-order Duffing nonlinearity, suitable laser drives could generate and optimize different nonlinearity orders related to the spectral shape of the optical force.

Such control could be exploited in a variety of contexts. At the level of individual resonators, optically induced Duffing nonlinearities may enable low-power bifurcation devices, nonlinear sensing schemes, mechanical logic elements, and tunable susceptibilities for noise manipulation and squeezing. In multimode systems, engineered nonlinear couplings provide a promising platform for exploring collective nonlinear dynamics, synchronization, nonlinear transport, and driven-dissipative phases in coupled oscillator networks [28, 24, 34]. In particular, the realization of controllable cross-Duffing interactions could enable nonlinear topological phenomena and soliton formation in optomechanical lattices [7, 9]. In the quantum regime, nonlinear backaction provides a powerful approach to the engineering of nonclassical quantum states and mechanical qubits [3, 39].

More broadly, these results establish cavity optomechanics as a versatile framework for realizing reconfigurable nonlinear dynamics at the nanoscale. Combined with the scalability of nano-optomechanical architectures and the possibility of measurement and control down to the quantum regime, this provides a powerful setting for exploiting controlled nonlinearity for fundamental studies as well as technological applications of nonlinear resonators and networks.

Supplementary Information

I Nonlinear optical backaction

In this section, we present a detailed analysis of optical backaction in fast-cavity optomechanical systems [2]. We start from the equation of motion

$$m\ddot{x}=-kx-m\gamma\dot{x}+F_{\text{opt}},$$ (S1)

for a mechanical resonator subject to radiation pressure force $F_{\text{opt}}$, given in the main text as Eq. (1). Optical backaction arises as the optical force

$$F_{\text{opt}}=\hbar Gn_{\text{c}}=\hbar Gn_{\text{max}}h\left(\frac{2\Delta_{0}}{\kappa}+\frac{2Gx}{\kappa}\right)$$ (S2)

depends instantaneously on the mechanical displacement $x$, through the optomechanical coupling $G$ and the Lorentzian cavity response $h(u)$. This treatment is valid in the fast-cavity limit where the cavity linewidth $\kappa\gg\Omega_{\text{m}}$ is much larger than the mechanical frequency, so that the photon dynamics can be eliminated adiabatically.

To analyse the nonlinear relation between $F_{\text{opt}}$ and $x$, we expand $h(u)$ into its Taylor series around $u_{0}=2\Delta_{0}/\kappa$ and obtain

$$F_{\text{opt}}=\hbar Gn_{\text{max}}\Bigg[h(u_{0})+h^{\prime}(u_{0})\left(\frac{2Gx}{\kappa}\right)$$ (S3)

$$\quad+\frac{h^{\prime\prime}(u_{0})}{2}\left(\frac{2Gx}{\kappa}\right)^{2}+\frac{h^{\prime\prime\prime}(u_{0})}{6}\left(\frac{2Gx}{\kappa}\right)^{3}+\mathcal{O}(x^{4})\Bigg]$$

The constant term $F_{\text{opt}}^{\text{DC}}=\hbar Gn_{\text{max}}h(u_{0})$ represents a static radiation pressure, which displaces the resonator’s equilibrium position $x_{0}$. Conventionally, this is accounted for by redefining $x\to x-x_{0}$.

The linear term

$$F_{\text{opt}}^{(1)}=\frac{2\hbar n_{\text{max}}h^{\prime}(u_{0})G^{2}}{\kappa}x=-k_{\text{opt}}x$$ (S4)

gives rise to the conventional, linear optical spring effect, characterised by spring constant $k_{\text{opt}}$ that supplements the intrinsic mechanical spring constant $k=m\Omega_{\text{m}}^{2}$. It shifts the mechanical frequency by

$$\delta\Omega =\Omega-\Omega_{\text{m}}=\sqrt{\frac{k+k_{\text{opt}}}{m}}-\sqrt{\frac{k}{m}}$$ (S5)

$$=\Omega_{\text{m}}\left(\sqrt{1+k_{\text{opt}}/k}-1\right)$$ (S6)

$$\approx\frac{\Omega_{\text{m}}k_{\text{opt}}}{2k}=\frac{k_{\text{opt}}}{2m\Omega_{\text{m}}},$$ (S7)

where the approximation is valid for small optical spring $k_{\text{opt}}\ll k$. Plugging in $k_{\text{opt}}$ yields

$$\delta\Omega =-\frac{\hbar n_{\text{max}}h^{\prime}(u_{0})G^{2}}{m\Omega_{\text{m}}\kappa}=-\frac{2n_{\text{max}}h^{\prime}(u_{0})g_{0}^{2}}{\kappa},$$ (S8)

where we have substituted the vacuum optomechanical coupling rate $g_{0}=Gx_{\text{zpf}}$ and zero-point fluctuation amplitude $x_{\text{zpf}}=\sqrt{\hbar/(2m\Omega_{\text{m}})}$. This relation for the linear optical spring shift is given in the main text as Eq. (4).

Next, we express the mechanical equation of motion using the normalized displacement $z=x/x_{\text{zpf}}$,

$$\ddot{z}=-\Omega_{\text{m}}^{2}z-\gamma\dot{z}+f_{\text{opt}}(z),$$ (S9)

where $f_{\text{opt}}(z)=F_{\text{opt}}/(mx_{\text{zpf}})$ is the normalized optical force. For reference, we give the linear term in the normalized force as

$$f_{\text{opt}}^{(1)}=\frac{4\Omega_{\text{m}}n_{\text{max}}h^{\prime}(u_{0})g_{0}^{2}}{\kappa}z=-2\Omega_{\text{m}}\delta\Omega\,z$$ (S10)

The quadratic term in Eq. (S3),

$$F_{\text{opt}}^{(2)} =\frac{2\hbar G^{3}n_{\text{max}}h^{\prime\prime}(u_{0})}{\kappa^{2}}x^{2},$$ (S11)

yields a contribution

$$f_{\text{opt}}^{(2)}=\frac{4\Omega_{\text{m}}n_{\text{max}}h^{\prime\prime}(u_{0})g_{0}^{3}}{\kappa^{2}}z^{2}=-\alpha z^{2}$$ (S12)

to the normalized optical force. However, precisely at the spring laser detunings $u_{0}=u_{\pm}=\pm 1/\sqrt{3}$ that we use, where the optical spring shift $\delta\Omega$ is maximally positive or negative, respectively, the second derivative $h^{\prime\prime}(u_{\pm})$ vanishes. Consequently, the quadratic nonlinear coefficient $\alpha$ induced by the spring laser is zero in all of our experiments.

In contrast, the cubic term in Eq. (S3),

$$F_{\text{opt}}^{(3)} =\frac{8\hbar n_{\text{max}}h^{\prime\prime\prime}(u_{0})G^{4}}{6\kappa^{3}}x^{3},$$ (S13)

does not vanish at $u_{0}=u_{\pm}$. When normalized, the cubic optical force term reads

$$f_{\text{opt}}^{(3)} =\frac{8\Omega_{\text{m}}n_{\text{max}}h^{\prime\prime\prime}(u_{0})g_{0}^{4}}{3\kappa^{3}}z^{3}=-\beta z^{3},$$ (S14)

where $\beta$ is the cubic nonlinear coefficient. This coefficient can be conveniently written using the linear spring shift $\delta\Omega$, as

$$\beta=\Omega_{\text{m}}\delta\Omega\frac{4h^{\prime\prime\prime}(u_{0})}{3h^{\prime}(u_{0})}\frac{g_{0}^{2}}{\kappa^{2}}.$$ (S15)

This expression is valid for any spring laser detuning $u_{0}$. At $u_{0}=u_{\pm}$, we find $h^{\prime\prime\prime}(u_{\pm})/h^{\prime}(u_{\pm})=-9/2$, so that we can write

$$\beta=-6\Omega_{\text{m}}\delta\Omega\frac{g_{0}^{2}}{\kappa^{2}},$$ (S16)

given in the main text as Eq. (6).

II Frequency response of a Duffing resonator

In this section, we analyse the frequency response of a Duffing resonator to a strong ‘pump’ force, optionally supplemented by a weak ‘probe’ force. We do so using the method of harmonic balance [21]. This approach, commonly used to evaluate the dynamics of Duffing oscillators [15], relies on expanding the displacement $z(t)$ as a Fourier series and then truncating that expansion to a specified order, while balancing the amplitudes of the remaining harmonics.

We start from the equation of motion

$$\ddot{z}+\gamma\dot{z}+\Omega^{2}z+\beta z^{3}=f(t)$$ (S17)

that describes the evolution of a general Duffing oscillator subject to a two-tone driving force

$$f(t)=f_{1}\cos(\omega_{1}t+\psi_{1})+f_{2}\cos(\omega_{2}t+\psi_{2})$$ (S18)

that includes a strong ‘pump’ tone with amplitude $f_{1}$ and frequency $\omega_{1}$ and a weak ‘probe’ tone with amplitude $f_{2}$ and frequency $\omega_{2}$. For simplicity, we assume that the frequencies $\omega_{1},\omega_{2}$ are incommensurate, such that the phase offsets are not essential and we can set $\phi_{1}=\phi_{2}=0$ without affecting the dynamics.

In response to this composite force, we seek an approximate solution to the resonator evolution of the form

$$z(t)=a_{1}\cos(\omega_{1}t+\phi_{1})+a_{2}\cos(\omega_{2}t+\phi_{2}).$$ (S19)

We insert Eq. (S19) into the Duffing equation (S17) and balance the fundamental harmonics at $\omega_{1}$ and $\omega_{2}$. First, we solve for the pump amplitude $a_{1}$ while neglecting the probe amplitude ($a_{2}=0$). This results in the well-known equation

$$a_{1}^{2}\left[\left(\omega_{1}^{2}-\Omega^{2}-\frac{3}{4}\beta a_{1}^{2}\right)^{2}+(\gamma\omega_{1})^{2}\right]=f_{1}^{2}$$ (S20)

that describes the approximate frequency response of a Duffing oscillator to harmonic driving. This equation underlies the theory curves plotted in Fig. 1d. For driving $f_{1}<f_{1}^{c}$ below the critical driving amplitude, Eq. (S20) admits a single solution for the amplitude $a_{1}$, while for strong driving $f_{1}>f_{1}^{c}$ the amplitude becomes multi-valued and hysteresis emerges.

Next, we solve for the probe amplitude $a_{2}$ in the presence of the pump. After balancing the fundamental harmonics, we find

$$a_{2}^{2}\left[\left(\omega_{2}^{2}-\Omega^{2}-\frac{3}{4}\beta(2a_{1}^{2}+a_{2}^{2})\right)^{2}+(\gamma\omega_{1})^{2}\right]=f_{2}^{2}$$ (S21)

Assuming that $a_{1}\gg a_{2}$ in the term in brackets, we arrive at the approximation

$$a_{2}^{2}=\frac{f_{2}^{2}}{(\mu^{2}-\omega_{2}^{2})^{2}+\gamma^{2}\omega_{2}^{2}}$$ (S22)

for the squared amplitude of motion at the probe frequency. This corresponds exactly to the susceptibility of a linear harmonic oscillator with natural frequency $\mu$. The shifted frequency $\mu$ satisfies the expression

$$\mu^{2}=\Omega^{2}+\frac{3}{2}\beta a_{1}^{2}.$$ (S23)

After plugging in Eq. (S16) for the radiation-pressure-induced $\beta$, we arrive at the expression

$$\delta\Omega_{\text{NL}}=\mu-\Omega\approx\frac{3}{4}\frac{\beta A^{2}}{\Omega_{\text{m}}}=-\frac{9}{2}\delta\Omega\frac{g_{0}^{2}}{\kappa^{2}}A^{2}$$ (S24)

that describes this nonlinear contribution to the spring shift of the probe susceptibility, valid under the assumption $\delta\Omega_{\text{NL}}\ll\Omega$. This equation is given in the main text as Eq. (8). Although Eq. (S24) is derived for a coherent probe tone, it equally predicts the nonlinearly shifted response to weak, incoherent thermal forces, as shown in Fig. 2.

III Effective nonlinear interactions

In this section, we present a detailed analysis of the optical backaction acting between a pair of mechanical modes coupled to a common cavity. We start from the two equations of motion

$$\ddot{z}_{j}=-\Omega_{\text{m},j}^{2}z-\gamma_{j}\dot{z}_{j}+f_{\text{opt},j}$$ (S25)

for the normalized amplitudes $z_{j}=x_{j}/x_{\text{zpf},j}$ of the two resonators $j=1,2$. We write the normalized optical force, arising for now from a single pump laser, as

$$f_{\text{opt},j}=\frac{\hbar Gn_{\text{c}}}{m_{j}x_{\text{zpf},j}}=2\Omega_{\text{m},j}g_{0,j}n_{\text{c}},$$ (S26)

where $n_{\text{c}}$ is total intracavity photon number. This is given by $n_{\text{c}}=n_{\text{max}}h(u_{0}+\delta u)$, where, importantly, the normalized optomechanical shift in the cavity frequency

$$\delta u=\frac{2}{\kappa}\sum_{j}g_{0,j}z_{j}$$ (S27)

now depends on the displacement $z_{j}$ of both resonators.

We first consider driving by a single pump laser tuned to either of the maximal spring shift points $u_{0}=u_{\pm}=\pm 1/\sqrt{3}$. Expanding $n_{\text{c}}$ to third order in $\delta u$ then yields

$$n_{\text{c}} =n_{\text{max}}\left[h(u_{\pm})+h^{\prime}(u_{\pm})\delta u+\frac{1}{6}h^{\prime\prime\prime}(u_{\pm})(\delta u)^{3}\right]$$ (S28)

$$=n_{\text{max}}\left[h(u_{\pm})+h^{\prime}(u_{\pm})\left(\delta u-\frac{3}{4}(\delta u)^{3}\right)\right],$$ (S29)

where we have used $h^{\prime\prime}(u_{\pm})=0$ and $h^{\prime\prime\prime}(u_{\pm})=-(9/2)h^{\prime}(u_{\pm})$. We plug this into the equation of motions to obtain

$$\ddot{z}_{j}= -\Omega_{\text{m},j}^{2}z-\gamma_{j}\dot{z}_{j}+2\Omega_{\text{m},j}g_{0,j}$$ (S30)

$$\times n_{\text{max}}\left[h(u_{\pm})+h^{\prime}(u_{\pm})\left(\delta u-\frac{3}{4}(\delta u)^{3}\right)\right]$$ (S31)

Again, we remove the static radiation pressure contribution by redefining the displacement origin and get

$$\ddot{z}_{j}= -\Omega_{\text{m},j}^{2}z-\gamma_{j}\dot{z}_{j}+2\Omega_{\text{m},j}g_{0,j}$$

$$\times n_{\text{max}}h^{\prime}(u_{\pm})\left(\delta u-\frac{3}{4}(\delta u)^{3}\right).$$ (S32)

We analyse the optical force term linear in $\delta u$,

$$f_{\text{opt},1}^{(1)} =2\Omega_{\text{m},1}\frac{2g_{0,1}n_{\text{max}}h^{\prime}(u_{\pm})}{\kappa}\left(g_{0,1}z_{1}+g_{0,2}z_{2}\right)$$ (S33)

$$=-2\Omega_{\text{m},1}\delta\Omega_{1}\,z_{1}-2\Omega_{\text{m},1}\sqrt{\delta\Omega_{1}\delta\Omega_{2}}\,z_{2},$$ (S34)

that acts on resonator $1$—the equivalent for resonator $2$ is obtained by exchanging indices $1\leftrightarrow 2$. The first term in Eq. (S34) represents the intra-resonator optical spring $k_{\text{opt},1}=2m_{1}\Omega_{\text{m},1}\delta\Omega_{1}$. Meanwhile, the cross-term in the optical force, expressed in physical units as

$$F_{\text{opt},1}^{(1,\leftrightarrow)} =m_{1}x_{\text{zpf},1}f_{\text{opt},1}^{(1)}=-2m_{1}\frac{x_{\text{zpf},1}}{x_{\text{zpf},2}}\Omega_{\text{m},1}\sqrt{\delta\Omega_{1}\delta\Omega_{2}}\,x_{2}$$

$$=-\sqrt{2\Omega_{\text{m},1}m_{1}\delta\Omega_{1}}\sqrt{2\Omega_{\text{m},2}m_{2}\delta\Omega_{2}}\,x_{2}$$

$$=-\sqrt{k_{\text{opt},1}k_{\text{opt},2}}\,x_{2}=-k_{\leftrightarrow}\,x_{2}$$ (S35)

represents an inter-resonator optical spring. This spring dynamically couples the displacements $x_{1},x_{2}$ with spring constant $k_{\leftrightarrow}=\sqrt{k_{\text{opt},1}k_{\text{opt},2}}$.

For a pair of high-quality mechanical resonators separated by a large frequency difference $\Delta\Omega_{12}=\Omega_{1}-\Omega_{2}\gg\gamma_{j}$, as we have in our experiment, inducing a static cross-resonator spring $k_{\leftrightarrow}$ has virtually no effect on their dynamics. We can, however, stimulate frequency conversion between the two modes by modulating $k_{\leftrightarrow}$ at $\Delta\Omega_{12}$. In experiment, this is achieved by modulating the intensity $P_{\text{in}}=P_{\text{in},0}(1+c_{\text{h}}\cos(\Delta\Omega_{12}t))$ of the spring laser.

After absorbing the intra-resonator spring shifts into the mechanical frequencies $\Omega_{j}=\Omega_{\text{m},j}+\delta\Omega_{j}$, neglecting the cubic term for now and removing terms that are not (mutually) resonant, we arrive at

$$\ddot{z}_{1}=-\Omega_{1}^{2}z-\gamma_{1}\dot{z}_{1}-2\Omega_{\text{m},1}c_{\text{h}}\cos(\Delta\Omega_{12}t)\sqrt{\delta\Omega_{1}\delta\Omega_{2}}\,z_{2}$$ (S36)

For convenience, we write the two displacements as

$$z_{j}(t)=\operatorname{Re}[a_{j}(t)e^{-\mathrm{i}\Omega_{j}t}]=\frac{a_{j}(t)e^{-\mathrm{i}\Omega_{j}t}+a_{j}^{*}(t)e^{\mathrm{i}\Omega_{j}t}}{2}$$ (S37)

in terms of the complex amplitudes $a_{j}(t)$. We plug these into the equation of motion (S36), neglect the terms that scale with $\gamma_{j}/\Omega_{j},J/\Omega_{j}\ll 1$ and remove the counter-rotating terms, to obtain the coupled differential equations

$$\begin{aligned} \dot{a}_{1}&=-\frac{\gamma_{1}}{2}a_{1}-\mathrm{i}\frac{c_{\text{h}}\sqrt{\delta\Omega_{1}\delta\Omega_{2}}}{2}a_{2}\\ \dot{a}_{2}&=-\frac{\gamma_{2}}{2}a_{2}-\mathrm{i}\frac{c_{\text{h}}\sqrt{\delta\Omega_{1}\delta\Omega_{2}}}{2}a_{1}\end{aligned}$$ (S38)

for the complex amplitudes. From these, we immediately identify the linear coupling rate

$$J=\frac{c_{\text{h}}\sqrt{\delta\Omega_{1}\delta\Omega_{2}}}{2}.$$ (S39)

Next, we analyse the nonlinearity in the cross-resonator spring $k_{\leftrightarrow}$, induced by the cubic term $\propto(\delta u)^{3}$ in Eq. (S32). Before doing so, we isolate this (cross-Kerr) nonlinearity from the intra-resonator Duffing nonlinearity by driving the cavity with two spring lasers $a$ and $b$ at opposite detunings $u_{0}^{a}=u_{-}=-1/\sqrt{3}$, $u_{0}^{b}=u_{+}=1/\sqrt{3}$. Noting that the difference between $u_{-}$ and $u_{+}$ is larger than any mechanical frequency, we can neglect the optical beating between the two lasers and decompose the total photon number

$$n_{\text{c}}=n_{\text{c}}^{a}+n_{\text{c}}^{b}=n_{\text{max}}^{a}h\left(u_{-}+\delta u\right)+n_{\text{max}}^{b}h\left(u_{+}+\delta u\right)$$ (S40)

into two independent populations $n_{\text{c}}^{a,b}$, one sustained by each laser.

To induce resonant coupling between the mechanical modes, we modulate one of the laser intensities. Here we modulate the positive spring laser intensity $n_{\text{max}}^{a}=n_{\text{max}}(1+f(t))$ while keeping the negative spring laser at fixed intensity $n_{\text{max}}^{b}=n_{\text{max}}$, equal to the average $n_{\text{max}}^{a}$. After expanding Eq. (S40) to third order, we find that in the total photon number

$$n_{\text{c}} =n_{\text{max}}2h(u_{\pm})$$

$$+n_{\text{max}}f(t)\left[h(u_{+})+h^{\prime}(u_{+})\delta u+\frac{1}{6}h^{\prime\prime\prime}(u_{+})(\delta u)^{3}\right]$$

$$=n_{\text{max}}2h(u_{\pm})$$ (S41)

$$+n_{\text{max}}f(t)\left[h(u_{+})+h^{\prime}(u_{+})\left(\delta u-\frac{3}{4}(\delta u)^{3}\right)\right]$$

the static linear and cubic spring shift terms cancel, leaving only the modulated part.

We plug this into the equation of motion, after discarding the static radiation pressure contribution. Moreover, as the modulation $f(t)=c_{\text{h}}\cos(\Delta\Omega_{12}t)$ is not resonant with any of the mechanical frequencies $\Omega_{1},\Omega_{2}$, we also discard the modulated direct radiation pressure contribution. We obtain

$$\ddot{z}_{1}= -\Omega_{\text{m},1}^{2}z_{1}-\gamma_{1}\dot{z}_{1}+2\Omega_{\text{m},1}g_{0,1}$$

$$\times f(t)n_{\text{max}}h^{\prime}(u_{\pm})\delta u\left[1-\frac{3}{4}(\delta u)^{2}\right]$$

$$= -\Omega_{\text{m},1}^{2}z_{1}-\gamma_{1}\dot{z}_{1}-2\Omega_{\text{m},1}f(t)\frac{2g_{0,1}n_{\text{max}}h^{\prime}(u_{\pm})}{\kappa}$$

$$\times\left[g_{0,1}z_{1}+g_{0,2}z_{2}\right]\left(1-\frac{3}{4}(\delta u)^{2}\right)$$

and the equivalent for $z_{2}$ after exchanging indices $1\leftrightarrow 2$. This can be conveniently expressed in terms of the linear spring shifts as

$$\ddot{z}_{1}= -\Omega_{\text{m},1}^{2}z_{1}-\gamma_{1}\dot{z}_{1}-2\Omega_{\text{m},1}f(t)\sqrt{\delta\Omega_{1}}$$

$$\times\Bigg[\sqrt{\delta\Omega_{1}}z_{1}\left(1-\frac{3}{\kappa^{2}}\left(g_{0,1}z_{1}+g_{0,2}z_{2}\right)^{2}\right)$$

$$\quad+\sqrt{\delta\Omega_{2}}z_{2}\left(1-\frac{3}{\kappa^{2}}(g_{0,1}z_{1}+g_{0,2}z_{2})^{2}\right)\Bigg]$$ (S42)

Again, we express the two displacements in terms of the complex amplitudes $a_{j}(t)$, employing the computer algebra software Mathematica to plug these into the coupled nonlinear equation of motion. After neglecting counter-rotating terms and terms scaling with $\gamma_{j}/\Omega_{j},J/\Omega_{j}\ll 1$, we obtain the dynamical equations

$$\dot{a}_{1}= -\frac{\gamma_{1}}{2}a_{1}$$ (S43)

$$-\mathrm{i}Ja_{2}(t)\left[1-\frac{9}{4}\left(\frac{g_{0,1}^{2}}{\kappa^{2}}|a_{1}(t)|^{2}+\frac{g_{0,2}^{2}}{\kappa^{2}}|a_{2}t)|^{2}\right)\right]$$

$$+\mathrm{i}Ja_{1}(t)\frac{9}{2}\left[\frac{g_{0,1}^{2}}{\kappa^{2}}\operatorname{Re}[a_{1}^{*}a_{2}]\right]$$

$$\dot{a}_{2}= -\frac{\gamma_{2}}{2}a_{2}$$

$$-\mathrm{i}Ja_{1}(t)\left[1-\frac{9}{4}\left(\frac{g_{0,1}^{2}}{\kappa^{2}}|a_{1}(t)|^{2}+\frac{g_{0,2}^{2}}{\kappa^{2}}|a_{2}t)|^{2}\right)\right]$$

$$+\mathrm{i}Ja_{2}(t)\frac{9}{2}\left[\frac{g_{0,2}^{2}}{\kappa^{2}}\operatorname{Re}[a_{2}^{*}a_{1}]\right].$$

for the complex amplitudes. Importantly, in each second term we recognize a coupling rate

$$J_{\text{NL}}=J\left[1-\frac{9}{4}\left(\frac{g_{0,1}^{2}}{\kappa^{2}}|a_{1}(t)|^{2}+\frac{g_{0,2}^{2}}{\kappa^{2}}|a_{2}t)|^{2}\right)\right]$$ (S44)

that includes a nonlinear correction depending on the squared amplitudes $|a_{1}(t)|^{2}$ and $|a_{1}(t)|^{2}$.

We initialize the nonlinear ringdown experiment in Fig. 4 with excitation in one mode only. The complex amplitude $a_{1}(t)$ and $a_{2}(t)$ then remain $\pi/2$ out of phase throughout, so that the correlation $\operatorname{Re}[a_{2}^{*}a_{1}]=0$ vanishes. The third term in each equation (S43) is thus identically zero, leaving us only with a nonlinearly shifted coupling rate $J_{\text{NL}}$.

IV Simulating geometric nonlinearity

We estimate the intrinsic geometric nonlinearity from finite element modelling in COMSOL Multiphysics. To do so, we first simulate the linear modeshape of the relevant flexural eigenmode of the device geometry, as shown in Fig. S1. We then perform a series of stationary analyses that include geometric nonlinearity. In this sweep, we impose the linear modeshape scaled by an increasing amplitude factor on the device structure, while for every amplitude, we integrate the total potential energy stored in the strain of the nanobeam [11]. This results in a potential energy $U(x)=\frac{1}{2}k_{1}x^{2}+\frac{1}{4}k_{3}x^{4}$ that is quadratic-quartic in the mode displacement $x$, where $k_{1}$ is the linear spring constant and $k_{3}$ the coefficient for the cubic restoring force.

We fit the potential $U(x)$ obtained from the simulations to extract $k_{1}$ and $k_{3}$. We verify that $k_{1}$, when combined with the effective mass $m_{\text{eff}}$ obtained from the same simulation, produces the correct eigenmode frequency $\Omega_{\text{m}}=\sqrt{k_{1}/m_{\text{eff}}}$. We convert $k_{3}$ into a cubic coefficient using the relation

$$\beta=\frac{k_{3}}{m_{\text{eff}}}x_{\text{zpf}}^{2}$$ (S45)

to produce an estimated value of $\beta\approx 300$ Hz².

V Experimental methods

The experimental set-up, methods, and device employed here have been described in detail in refs. [4, 46, 43]. We briefly summarize the relevant aspects below.

FIG. S1: Nanobeam flexural modes. Simulated modeshapes for the first five flexural mechanical modes. All experiments in this work employed the fourth mode, with the mode coupling experiments additionally employing the third mode.

V.1 Optomechanical system

The nano-optomechanical device under study comprised a suspended, sliced silicon nanobeam hosting a photonic crystal cavity and fabricated in-house using the methods described in [4]. The optical cavity coupled simultaneously to several non-degenerate flexural modes of the beam, whose spatial profiles are shown in Fig. S1. All experiments in this work employed the fourth mode, with the mode coupling experiments additionally employing the third mode.

Optical access was provided by free-space coupling, achieved by focusing a laser at normal incidence onto the device. Direct reflections were suppressed and cavity-emitted light was collected using a cross-polarized detection scheme. The device was mounted inside a room-temperature vacuum chamber (pressure $2\times 10^{-6}$ mbar) and positioned at the laser focus using a piezo-actuated precision stage. All quoted optical and mechanical parameters were characterized following the methods described in [4].

V.2 Excitation and readout

Spring shifts were induced by pumping the cavity with a strong ‘spring’ laser (Toptica CTL 1500/1550). Depending on the experiment, a secondary control laser was either used as a ‘force’ laser, or as an additional, opposing ‘spring’ laser. When used as a weak force laser, it was tuned to the cavity resonance with its intensity modulated to provide radiation pressure driving of the mechanical resonator. When used as an additional spring laser, its power was matched to the primary spring laser but operated at the opposite detuning to cancel the (nonlinear) spring effect. In that case, coherent radiation pressure driving was provided by weakly modulating either of the spring lasers instead. In the multimode coupling experiments, the intensity of a spring laser was modulated at the difference frequency of the resonator pair.

Electronic signals to drive or couple mechanical modes were generated by a high-frequency lock-in amplifier (LIA; Zurich Instruments UHFLI), amplified (Mini-Circuits ZHL-32A+ with 9 dB attenuation) and imprinted on the appropriate laser using fiber coupled electro-optic modulators (Thorlabs LN81S-FC and Covega Mach-10 056)

The instantaneous displacement $x_{j}(t)$ of both mechanical resonators was transduced onto the reflected intensity of light from a tertiary ‘detect’ laser (Toptica CTL 1550, incident power $P_{\text{det}}=1$ mW), far detuned from the cavity resonance ($\Delta_{\text{det}}=-2.5\kappa$,). The reflected detection light was spectrally isolated from control laser light using a tunable band-pass filter (DiCon) before being directed onto a fast, AC-coupled photodiode (New Focus 1811). To extract individual resonator signals, the photodiode output was demodulated by the LIA with a 3-dB-bandwidth of $50$ kHz (third order filter). The resulting complex amplitudes $a_{j}(t)$ were normalized to the detector voltage level corresponding to each resonator’s zero-point motion $x_{\text{zpf}}$, as obtained from thermally driven reference measurements without any modulations.

For time-resolved ring-down measurements, mechanical driving and coupling signals were continuously output by the LIA on separate output channels, which were then routed through individual radio frequency (RF) switches (Mini-Circuits ZYSWA-2-50DR+). An additional signal generator (Siglent SDG series) controlled the switches and triggered LIA data acquisition.

VI Extracting instantaneous oscillation frequencies

The complex resonator amplitude $a(t)=A(t)e^{\mathrm{i}\phi(t)}/x_{\text{zpf}}$ reported by the LIA encodes the amplitude $A(t)$ of resonator motion, as well as its phase relative to the LIA’s electronic local oscillator (LO). The demodulated phase evolves as

$$\phi(t)=\int_{t_{0}}^{t}\left[\Omega(\tau)-\omega_{\text{d}}\right]\,\mathrm{d}\tau$$ (S46)

where $\Omega(t)$ is the instantaneous frequency of the resonator’s motion, $\omega_{\text{d}}$ the LO frequency and $t_{0}$ the time the measurement was started. By taking the derivative of Eq. (S46), the instantaneous frequency is extracted $\Omega(t)=\phi^{\prime}(t)+\omega_{\text{d}}$. This expression was evaluated for the experimental data to obtain the frequency traces shown in Fig. 1f, bottom.

VII Estimating instantaneous coupling rates

In this section, we introduce a numerical method to estimate the instantaneous coupling $J_{\text{NL}}(t)$ from the experimentally obtained complex amplitudes $a_{j}(t)$ of the coupled resonator pair.

The starting point is the coupled system of differential equations

$$\dot{a}_{1} =-\frac{\gamma_{1}}{2}a_{1}-\mathrm{i}J_{\text{NL}}(t)a_{2},$$

$$\dot{a}_{2} =-\frac{\gamma_{2}}{2}a_{2}-\mathrm{i}J_{\text{NL}}(t)a_{1}$$

that describes the evolution of a pair of resonators resonantly coupled by an effective instantaneous rate $J_{\text{NL}}(t)$ (cf. Eq. (S38)). From these equations of motion, we calculate the ‘excess change’ in each resonator’s complex amplitude as

$$\begin{aligned} \dot{a}_{1}+\frac{\gamma_{1}}{2}a_{1}&=-\mathrm{i}J_{\text{NL}}(t)a_{2}\\ \dot{a}_{2}+\frac{\gamma_{1}}{2}a_{2}&=-\mathrm{i}J_{\text{NL}}(t)a_{1}.\end{aligned}$$ (S47)

Here, the left hand sides express the portion of the amplitude change $\dot{a}_{j}$ that is not attributed to the dissipation $\gamma_{1}a_{1}/2$. From these, an estimate for $J_{\text{NL}}(t)$ can be immediately obtained through dividing by $a_{2,1}$. However, these amplitudes will periodically tend to zero, giving imprecise results at those times when applied to real data.

To solve this, we square both and take the sum

$$\left|\dot{a}_{1}+\frac{\Gamma_{1}}{2}a_{1}\right|^{2}+\left|\dot{a}_{2}+\frac{\Gamma_{1}}{2}a_{2}\right|^{2} =J_{\text{NL}}(t)^{2}\left(|a_{1}|^{2}+|a_{2}|^{2}\right)$$ (S48)

The total energy $|a_{1}|^{2}+|a_{2}|^{2}$ does not vanish periodically as the resonators exchange energy, so that $J_{\text{NL}}(t)^{2}$ can be estimated reliably during the entire ringdown.