Bichromatic synchronized chaos in coupled optomechanical nanoresonators
BBichromatic synchronized chaos in coupled optomechanical nanoresonators
Guilhem Madiot, Franck Correia, Sylvain Barbay, and R´emy Braive
1, 2 Centre de Nanosciences et de Nanotechnologies,CNRS, Universit´e Paris-Saclay, Palaiseau, France Universit´e de Paris, F-75006 Paris, France
Synchronization and chaos are two well known and ubiquitous phenomena in nature. Interestingly,under specific conditions, coupled chaotic systems can display synchronization in some of theirobservables. Here, we experimentally investigate bichromatic synchronization on the route to chaosof two non-identical mechanically coupled optomechanical nanocavities. Electromechanical near-resonant excitation of one of the resonators evidences hysteretic behaviors of the coupled mechanicalmodes which can, under amplitude modulation, reach the chaotic regime. The observations, allowingto measure directly the full phase space of the system, are accurately modeled by coupled periodicallyforced Duffing resonators thanks to a complete calibration of the experimental parameters. Thisshows that, besides chaos transfer from the mechanical to the optical frequency domain, spatialchaos transfer between the two nonidentical subsystems occurs. Upon simultaneous excitationsof the coupled membranes modes, we also demonstrate bichromatic chaos synchronization betweenquadratures at the two distinct carrier frequencies of the normal modes. Their respective quadratureamplitudes are consistently synchronized thanks to the modal orthogonality breaking induced by thenonlinearity. Meanwhile, their phases show complex dynamics with imperfect synchronization in thechaotic regime. Our generic model agrees again quantitatively with the observed synchronizationdynamics. These results set the ground for the experimental study of yet unexplored collectivedynamics of e.g synchronization in arrays of strongly coupled, nanoscale nonlinear oscillators forapplications ranging from precise measurements to multispectral chaotic encryption and random bitgeneration, and to analog computing, to mention a few.
I. INTRODUCTION
Synchronization is a common phenomenon where anoscillating physical quantity tends to develop a phasepreference when coupled to a drive or to another os-cillating system [1]. Since Huygens well-known prelim-inary studies of synchronization between two mechani-cally coupled pendulum clocks [2], such phenomenon isubiquitous in many fields of sciences e.g. in physics [3–6], chemistry [7], biology [8], ecology [9], economy [10]and even in sociology [11]. Although it may look like aunique phenomenon, synchronization features may occurwithin different dynamical regime in which the systemis set. Among these regimes, chaotic dynamics, inter-disciplinary by essence, has attracted a lot of attention.Even though synchronization of chaotic system may seemcounterintuitive [12], the investigation and implementa-tion of such a phenomenon is particularly meaningful andhas applications from large scale systems as in meteoro-logic [13, 14] and atmospheric general circulation [15, 16]models, to small scales with atomic clocks [17]. Beyondthis, chaos in arrays of synchronized systems is of greatinterest for generating random numbers [18] potentiallyused for robust and encrypted communications in op-tics [19–22] and astronomy [23]. It could also be rele-vant in artificial neuron networks model [24, 25], ultra-precision measurements [26], fundamental tests on thephysical conditions for classical dynamics [27].Signatures of synchronization can be imprinted on thedynamics of both quadrature of the investigated signal,i.e. the amplitude and the phase. However, this lat-ter exhibits a much richer dynamic which may be evi- denced in time by periodic, erratically distributed jumpsand many others (as described in Kuramoto model [7]).Imperfect phase synchronization [28] refers to the casewhere intermittent phase slips occur while the amplitudeis still synchronized. This phenomenon is barely stud-ied experimentally although it constitutes an universalparadigm anchored in a wide range of disciplines such asphysics [29, 30], electronics [31], and even neurosciences[32] to name a few. A complete understanding of theseconcepts can be reached with a classical model of cou-pled driven resonators. Therefore, the dynamics result-ing from this basic model can be applied and is also rel-evant to many other fundamental studies such as spatio-temporal phase transition physics [33–35], patterns for-mations [36] or synchronized neutrinos oscillations [37]encountered in particle physics. In this frame, mechan-ical and optomechanical systems constitute a platformof choice thanks to their experimental adaptability tohost most of fundamental concepts of nonlinear dynam-ics. Many nonlinear dynamical effects have been demon-strated in single or coupled Nano ElectroMechanical Sys-tems (NEMS) [38–44], including synchronization [45–47]and chaos [48]. In line with these studies, chaos with op-tomechanical systems has been extensively studied the-oretically [49–51]. Only recently and using a single res-onator, experimental demonstrations of chaos in optome-chanical cavity have been reported [52, 53] with a strongoptical driving. Here, thanks to a forcing amplitude mod-ulation technique, we emphasize mechanical chaos usingtwo mechanically coupled optomechanical cavities. Be-yond chaos transfer from mechanics to optics, we evi-dence energy transfer between two non-identical mechan- a r X i v : . [ n li n . C D ] M a y ical resonators through the spring coupling. Despite theirintrinsic natural frequencies mismatch, chaos is simul-taneously generated at two distinct carrier frequenciesopening potential new avenues to synchronized multi-frequency data encryption [19–22] and random numbergeneration [18]. Interestingly, the quadrature amplitudesof the two chaotic signals at two different tones are syn-chronized. In parallel the quadrature phases evidence dif-ferent regimes from phase synchronization to imperfectphase synchronization through phase desynchronization.The unequivocal description of these regimes is enabledby the measurement method which give a direct accessto the dynamical variables, without making any use ofreconstructed signals. These phenomena can be under-stood and fully supported by numerical simulations basedon a classical model using Duffing resonators, which be-yond nanomechanics, can be used in many fields includ-ing superconducting Josephson amplifier [54], ionizationplasma [55] and complex spatiotemporal behaviors suchas chimera sates [33].Our study relies on the mechanical coupling betweentwo non-identical optomechanical systems: electro-capacitively driven Fabry-P´erot cavities which are de-scribed in Sec. II. After a preliminary study of theirlinear mechanical properties in Sec. III A allowing toextract useful modeling parameters, a strong resonantdriving field is applied and evidences hysteretic behav-ior in Sec. III B. As such, we experimentally demonstratein Sec. IV how low-frequency amplitude modulation ex-erted on driven coupled Duffing resonators can boost thenonlinearity in the material so strongly that it leads toa period-doubling cascade and then to chaos. Experi-mental bifurcation diagrams are reconstructed from themeasured time traces when either eigenmode is driven.In each configuration, the largest Lyapunov exponent iscomputed attesting the chaotic behaviors. When simul-taneously driven, the eigenmode orthogonality breakingdue to the Duffing nonlinearity is shown in Sec. V to allowthem to couple. We then investigate the synchroniza-tion of the eigenmode amplitude and phase responses.While the amplitude responses are found to be robustlysynchronized, we evidence several phase synchronizationregimes including complete synchronization, desynchro-nization and imperfect synchronization. The experimen-tal results are in good quantitative agreement with theproposed model. In Sec. VI we conclude. II. EXPERIMENTAL SYSTEM ANDMEASUREMENT SCHEME
The experimental system (see fig. 1a) consists of twocoupled mechanical micro-resonators made of a common260 nm thick InP layer suspended over a 380 nm airgap. Each membrane is a 10 × μ m rectangle piercedwith a square lattice of cylindrical holes. This array bothpermits to reduce the mechanical masses, increasing themechanical resonator natural frequencies and allows an lock-in amplifier AB <10 -5 mbar memb. A memb. B
550 600 650 700 7500,11 n o r m . r e f l e c t a n c e wavelength [nm] n m a)b) c)PdrefV tot spectrosupercont. laserHe-Nelaser x20 Fig. 1: a) excitation and detection setup including anoptical microscope image of the integrated system. Twosuspended membranes are attached by four mesas andbridged by a coupling nanobeam. Gold stripes arevisible underneath the membranes for electro-actuation.b) SEM micrograph of the coupling beam. c)Reflectance spectra for the two Fabry-P´erot cavitiesformed by membranes A (black) and B (blue) and theunderneath substrate.enhancement of the out-of-plane reflectivity [56]. Thetwo membranes are mechanically coupled through a 1 μ m wide and 1.5 μ m long bridge (see figs. 1a and 1b). Apair of gold interdigitated electrodes (IDEs) is positionedon the substrate below each resonator and allow for in-dependent actuation of the mechanical resonators. Thefabrication process is described in [57]. All measurementsare done at room temperature and the chip is placed ina vacuum chamber pumped below 10 − mbar.The system composed of the membrane plus the IDEsconstitutes a low-finesse optical cavity which we canprobe for measuring the mechanical displacement of thesuspended membranes. The reflectance spectra shown infig. 1c are measured by means of a supercontinuum laserreflected on the centers of the at-rest membranes. Thebeam is focused down to 5 μ m on the membrane centerwith a × λ = 633 nm which is usedfor the optical readout the membrane displacements in-duced by the capacitive field generated by the IDEs whensubmitted to a voltage V tot . The electro-capacitive forceexerted on one membrane reads [40]: F = − dCdx V tot (1)where C ( x ) is the position-dependent capacitance of themembrane/IDEs system. The driving voltage is V tot = V dc + V ac cos( ω d t ) where V dc is a static voltage and V ac isthe amplitude of the AC driving at frequency ω d . Eachmembrane can be independently probed. A photodetec-tor converts the reflected optical field into an electricalsignal sent to a lock-in amplifier to access the demodu-lation amplitude ηR A,B and phase θ A,B where A or Bhere refer to the probed resonator and η converts a givenmeasured voltage to the corresponding mechanical dis-placement R A,B . The calibration of η can be obtainedthanks to a Michelson interferometer and we estimate theconversion constant η ≈ . − (see appendix A).The RF input is demodulated with a passband filter cen-tered at ω d enabling a higher signal-to-noise ratio and thedemodulator internal frequency is locked to the appliedharmonic excitation of the IDEs. III. SINGLY DRIVEN COUPLED RESONATORSA. Linear Regime
To build the model describing the system, we first in-troduce the potential energy associated to uncoupled har-monic resonators: U lin ( x A , x B ) = 12 ( ω A, x A + ω B, x B ) (2)with ω A, (resp. ω B, ) the natural frequency of themode resonator A (resp. B). The mechanical beamjoining the two membranes contributes to the dynamicsthrough the coupling spring constant G. This results in acoupling potential U coup ( x A , x B ) = Gx A x B + G ( x A + x B )with the resonators self-coupled frequencies ω A and ω B such that ω A = ω A, + G and ω B = ω B, + G . We con-sider mass normalized physical quantities in our model.We introduce linear mechanical damping terms Γ A andΓ B for each resonator and describe the problem with thefollowing coupled master equations [58]: ¨ x A + Γ A ˙ x A + ∂U tot ∂x A =0¨ x B + Γ B ˙ x B + ∂U tot ∂x B = F B cos( ω d t ) (3)where the total potential energy is U tot = U lin + U coup .Note that the forcing term F B is the only non-zero right-hand side term because we only excite the membraneB. The stationary solutions of eq. (3) can be derived byassuming an oscillatory solution x A,B = r A,B cos( ω d t + θ A,B ) where r A,B is the mechanical displacement of eachmembrane.In this first experiment we will characterize the me-chanical response of the system in the linear regime.Membrane B is driven V dc = 0 .
5V and V ac = 1V to thecorresponding set of IDEs. The He-Ne laser is focusedon membrane B so that we record its response ampli-tude R B while the driving frequency ω d is swept. The driven mechanical system exhibits a large variety of nor-mal modes ranging from 1 to 10 MHz. We focus our at-tention on the lowest frequency modes corresponding tothe coupled fundamental modes of the membranes. Themodes centered respectively at ω − = 2 π × .
161 MHz and ω + = 2 π × .
369 MHz are identified as the symmetrical( − ) and antisymmetrical (+) normal modes.By fitting the theoretical response of resonator B(fig. 2b) we obtain the self-coupled frequencies ω A =2 π × .
187 MHz, ω B = 2 π × .
345 MHz, and the damp-ings Γ A = 2 π × . B = 2 π × . G/ω B ≈ π ×
130 kHz. Themode coupling is further attested on by the presence ofa destructive interference dip around 2.21 MHz which istypical of a Fano resonance [59–61] between the nearlyidentical resonators. In the linear regime and for a lowdriving amplitude, the Fano dip minimum is below thedetection noise floor but its presence on the spectrum isnevertheless clearly visible.The constitutive resonator frequency difference leadsthe normal mode at ω − (resp. ω + ) to be dominated bythe motion of membrane A (resp. B). This conclusion isconfirmed by measuring ω B − ω A with a dielectric tuningtechnique described in appendix B. Altogether it impliesthat ω − ≈ ω A and ω + ≈ ω B . Finite Elements Method(FEM) simulations confirm the normal modes respectivedisplacement fields and the energy imbalance due to anatural frequency mismatch (fig. 2a). The amplitude ofthe normal mode (+) is almost 4 times higher than am-plitude of mode ( − ) due to this imbalance. In the con-text of identical resonators, the strong coupling regimeis established when the criterion G/ω A > Γ A is satisfied[62]. This criterion applies in our experiment, however,since the resonator frequency mismatch is about twiceas large as the minimum normal mode splitting we arerather in an intermediate regime between the strong andweak coupling cases.When expanded, the expression for the electrocapaci-tive force eq. (1) includes a static component ∝ ( V dc + V ac /
2) that displaces the resonator by a negligible off-set plus an off-resonant term at frequency 2 ω d which isignored in our model. A measurement of the displace-ment amplitude at demodulation frequency 2 ω d indeedreveals an amplitude response less than 3% of the one ofthe driven mode amplitude at ω d . Therefore, the driv-ing force amplitude can be related to the experimentalparameters with F B = | m eff dCdx V dc V ac | where m eff =186pg is the effective mass computed at fundamental eigen-frequency of normal mode (+) by FEM. The resonantamplitude R B ( ω d = ω + ) is mapped over the parameterspace { V dc , V ac } as shown in fig. 2c. The stationary so-lutions of eq. (3) indicate that the resonators responsesare both linear with the strength F B . This linear depen-dence is independently checked (red lines in fig. 2c) with data lin. sol. data unstable sol. stable sol.2,1 2,2 2,3 2,4 2,510 -1 -V int V d c = V V ac =3V nonlinear regime B, w + [nm BW] linear regime -2 -1 R B [ n m B W ] noise floor 01234 V dc =2VR B V a c [ V ] -4 -3 -2 -1 0 1 2 3 4 V ac =3V R B V dc [V] R B [ n m B W ] w d / p [MHz] a) c) b) d) 𝜔 − 𝜔 + maxmin0 A BA B Fig. 2: a) FEM simulation of the out-of-plane displacement field for the symmetrical ( ω − ) and the antisymmetrical( ω + ) eigenmode taking into account the frequency mismatch. b) Spectral response of membrane B when driven with V ac = 1V and V dc = 0 . V ac , V dc ) parameter space for ω d = ω + . Horizontal and vertical slices are represented for respectively V ac = 3V and V dc = 2V with a linear fit (red line) in the non saturated regime. The iso- V eff dc V ac curves (black dashed) separate thelinear and Duffing-Duffing regimes. The internal stress (yellow line) parabolically shifts with V ac . d) Spectralresponse of membrane B at V ac = 3V and V dc = 2V. Data are recorded with upward and downward ω d sweeps (bluesymbols) and are fitted with the Duffing-Duffing model with stable (red line) and unstable (red dotted) solutions.All displacement measurement are performed with a domulation BW = 100Hz.varying V dc (by line) or V ac (by column) both allowingthe electro-capacitive force to be consistently calibrated: dC/dx ≈ . μ N/V .The model takes into account a small offset in the ef-fective static voltage V eff dc = V dc + V int where V int cor-responds to the internal stress of the membrane. Weobserve a small shift of the internal stress with increas-ing V ac due to the dielectric tuning induced by the staticcomponent V ac [63]. B. Nonlinear Regime
In fig. 2c, a domain of saturation of the mechanicalresponse settles when the product V eff dc V ac is greater than1 . . The corresponding iso- V eff dc V ac curves delimit athreshold between the linear and the nonlinear regimes.They are fitted (dashed black lines) using the data pointsat the frontier between the two domains.The system response in the nonlinear regime is shownon fig. 2d for V dc = 2V and V ac =3V and displays twohisteretic regions around ω − and ω + that are evidencedby sweeping ω d forward and backward. This saturationarises from intrinsic mechanical nonlinearities [64] whichcan be modeled thanks to a Duffing oscillator model [57, 65].We follow the same approach as in Sec. III A and in-troduce anharmonicity to the uncoupled harmonic res-onators potential energy through the nonlinearity β : U NL = U lin + 14 β [ x A + x B ] (4)The nonlinear dynamics is still governed by eq. (3) butconsidering the new total potential energy U tot = U coup + U NL . The stationary solutions for the resonators am-plitudes r A and r B and phases θ A and θ B can be de-scribed by the following set of equations (see details inappendix C): ˙ r A = − γ A r A + g r B sin( θ A − θ B )˙ r B = − γ B r B − g r A sin( θ A − θ B )+ f B θ B ) r A ˙ θ A = − r A (cid:20) δ − ∆ ω ) + 34 ˜ βr A (cid:21) + g r B cos( θ A − θ B ) r B ˙ θ B = − r B (cid:20) δ + 34 ˜ βr B (cid:21) + g r A cos( θ A − θ B ) + f B θ B ) (5)where we define the normalized forcing strength f B = F B /ω d , dampings γ i = Γ i /ω d , coupling g = G/ω d ,frequency mismatch ∆ ω = ( ω B − ω A ) /ω d , detuning δ = ( ω B − ω d ) /ω d , Duffing nonlinearity ˜ β = β/ω d andrescaled time ω − d . The variables r A and r B are left inunits of nanometer to compare with the experimental re-sults.By assuming the permanent regime ˙ r A = ˙ r B = ˙ θ A =˙ θ B = 0, in eq. (4), the frequency domain response ampli-tudes r B is numerically solved using the experimental pa-rameters extracted from fig. 2b and previously discussed.We adjust the solution on the experimental data by fit-ting with the remaining free parameter β and extract β = (2 π ) × . × − MHz .nm − . The resulting curveshown in fig. 2d is composed of a stable solution (red line)and an unstable solution (red dashed line). It has beenchecked that another model including only one anhar-monic resonator coupled with a harmonic resonator doesnot permit to describe the double bistability we experi-mentally observe. We conclude from the good agreementof our model with the experiment that a linear springcoupling satisfactorily describes the mechanical interac-tion between the membranes. IV. CHAOTIC DYNAMICS UNDERAMPLITUDE MODULATION
An additional low-frequency modulation signal isadded to the total voltage applied to the membrane Bset of IDEs. It now writes: V tot = V dc + V ac cos( ω d t ) + V p cos( ω p t ) with V p and ω p (cid:28) ω d the modulation ampli-tude and frequency respectively.We reduce the set of experimental variables by locking V dc = 2V and V ac = 3V while V p is used as the controlparameter to explore the dynamical changes of the sys-tem. The modulation frequency is set to ω p = 2 π × − ) eigenmode bistability curve at ω d = 2 π × V p is swept from 0 to 3V. Higher valuesare not reached in order to preserve the mechanical sys-tem from failure. For each value of V p , we record thesignal quadratures X A,B = R A,B cos( θ A,B ) and Y A,B = R A,B sin( θ A,B ) in real-time using thus accessing simulta-neously phase and amplitude components. The samplingrate is 500 kHz and each trace has a length of 100ms, en-suring that several hundreds of modulation periods arerecorded. In figs. 3b and 3e, we plot Y B ( t ) as a functionof X B ( t ). It actually corresponds to a 2D projection ofthe whole dynamical phase space. The Poincar´e sectionmade of the local maxima of Y B ( t ) as a function of V p is shown in (fig. 3a). Using Y B ( t ) rather than X B ( t ) isan arbitrary choice motivated by the higher amplitudeof the phase portraits along the Y axis. Additionally,each time trace is used to compute the largest Lyapunovexponent (LLE) shown below the diagram. We use theTISEAN package [66] routine implementing the Rosen-stein algorithm [67]. This calculation essentially relies onthe delay embedding reconstruction of the phase space inwhich initially close trajectories are compared over time.For a low value of the modulation voltage injected intothe normal mode ( − ), the Poincar´e section in fig. 3a re-sults in a closed single loop. In this limit-cycle oscillationregime the membranes oscillation envelopes are modu-lated at ω p . As the amplitude is modulated stronger,we observe two consecutive period-doubling bifurcationsat V p ≈ .
75V and V p ≈ .
5V prior to a window ofchaotic dynamics for a modulation amplitude higher than2 . ω d = 2 π × .
379 MHz.We construct the bifurcation diagrams still reading themotion of membrane B (fig. 3d). The phase portraits as-sociated to this case are shown in fig. 3e. The bifurcationdiagrams of eigenmode (+) also display a period-doublingroute to chaos structure [68] although the chaotic regimenow occurs around V p ≈ V p = 2 .
6V with a period-4 motion.Both experimental diagrams share a common dynamicsbut the bifurcation points significantly differ whether theeigenmode ( − ) or (+) is driven. This quantitative differ-ences between the eigenmodes dynamics result from theimbalanced energy injection in the normal modes sinceonly membrane B is driven. Identical measurements areperformed by reading the membrane A and are shownin appendix D. Both membranes basically settle in thesame dynamical regime under a given excitation.The bifurcation diagram for the ( − ) mechanical mode(resp. for the (+) mode) is numerically reproducedin fig. 3c (resp. in fig. 3f) using the Duffing-Duffingmodel developed previously at the driving frequencyFig. 3: Experimental and numerical bifurcation diagrams under single driving and by reading the displacement ofmembrane B. Measurement and simulations are performed by driving either the symmetrical (left column ω d = 2 π × .
164 MHz) or anti-symmetrical resonance (right column ω d = 2 π × .
379 MHz) with V dc = 2V, V ac =3Vand ω p = 2 π × V p and reading membraneB with the associated largest Lyapunov exponent (LLE). Note the broken axis. b)/e) Phase portraits at differentdynamical regimes. c)/f) Numerical simulation of bifurcation diagrams built from the maxima of quadrature w B with the driving frequencies ω d = 2 π × . ω d = 2 π × . ω d = 2 π × . ω d = 2 π × . f B and injecting the experimentally determinedparameters with˜ f B ( t ) = f B (cid:104) V p V dc cos( ω p t ) (cid:105) (6)This results in a modulated force amplitude F B ( t ) ∝ V ac ( V dc + V p cos( ω p t )). The additional terms arising fromthe expansion of eq. (1) are either static ( V p ) or slowlyoscillating (at ω p and 2 ω p ) and do not significantly par-ticipate to the dynamics in our system. Indeed, the forceassociated to these terms leads to a displacement smallerthan the resonant motion by an amount given by themechanical total quality factor Q m ≈ w B which corresponds to the observable Y B . The route to chaos by period doubling cascade is well captured by our model. The quantitative compar-ison with the experimental results yield a satisfactoryagreement, which is remarkable given the experimentalparameter uncertainties and the simplicity of the model.The period-doubling cascade bifurcation positions arecorroborated with a continuation method. Thanks tothe model, it is possible to track the origin of the chaoticdynamics in the force modulation and not in the cou-pling: indeed uncoupled membranes also display chaosunder modulation. Additional experimental and numer-ical analysis further show that the modulation frequencyplays a role in the appearance of the chaotic regimeand a large window of frequencies around the dampingtimescale Γ − A,B leads to chaos. However, our results showthat the transfer of the chaotic dynamics from one mem-brane to the other is possible in spite of the detuningbetween the membrane resonant frequencies. Further-more, the chaotic dynamics is also imprinted in the laserfield intensity thanks to the integrated Fabry-P´erot cav-ity. This mechanical-to-optical chaos transfer is an inter-esting concept to exploit in chaos based technologies.
V. NONLINEARLY COUPLED NORMALMODES
The previous experiment shows that although the nor-mal modes can be driven to chaotic regime using ampli-tude modulation, they have their own bifurcation points.Additionally, the eigenmodes ( − ) and (+) are no longerexpected to be orthogonal as soon as the Duffing regime isreached. This orthogonality breaking [69, 70] enables theeigenmodes to couple. In order to evidence this coupling,we take advantage of the chaotic dynamics generated byamplitude modulation. We inject energy in both nor-mal modes by using two resonant forces ( V − ac and V + ac )and modulate their evolution with a common modula-tion. The situation is therefore analog to driven coupledchaotic systems in which we will investigate on the syn-chronization properties and where the modulating signalplays the role of the drive. In this experiment, we caninterestingly excite the system and measure all the dy-namical variables exclusively through membrane B. Theresults that one would obtain by probing membrane Awould be identical but with lower amplitude.We adapt the previously derived equations and obtaina new system of governing equations: ¨ x A + Γ A ˙ x A + ω A x A + βx A − Gx B =0¨ x B + Γ B ˙ x B + ω B x B + βx B − Gx A = F − B cos( ω − d t )+ F + B cos( ω + d t )where F − B (resp. F + B ) is the force amplitude exerted oneigenmode ( − ) with driving frequency ω − d (resp. eigen-mode (+) at ω + d ). The stationary solutions are de-tailed in appendix E. These solutions highlight a non-linear coupling between the normal modes of the form˜ βr ± r ∓ . The non-autonomous equations resulting fromthe amplitude-modulation are further used for the nu-merical simulations by implementing the time-dependentstrengths previously presented in eq. (6).The new total voltage applied on membrane B set ofIDEs writes: V tot = V dc + V − ac cos( ω − d t ) + V + ac cos( ω + d t ) + V p cos( ω p t )with V dc = 2 V , V − ac = 3 . V , ω − d = 2 π × .
177 MHZ, V + ac = 0 . V , ω + d = 2 π × .
410 MHz and ω p = 2 π × V + ac < V − ac in order to compensate the re-sponse amplitudes imbalance that results from the mem-branes frequency mismatch. We place the laser spot onmembrane B and use two independent demodulators tosimultaneously access the signal amplitude and phase at ω − d ( R − B and θ − B ) and at ω + d ( R + B and θ + B ). By sweepingthe bifurcation parameter V p , a new diagram is built fromthe local maxima of signal quadratures Y − B = R − B sin( θ − B )and Y + B = R + B sin( θ + B ) that we record with demodula-tion bandwidth of 40 kHz. This allows a reduction ofthe crosstalk between the channels of about -5dB. Wenote that the diagram branches are broader than in the Fig. 4: a) experimental bifurcation diagrams andassociated calculated largest Lyapunov exponent (LLEin ms − ) built from the normal modes responsequadrature Y − B and Y + B in units of nm ×√ BW withswept parameter V p . The modes are driven at V − ac = 3 .
5V and V + ac = 0 .
5V at frequencies ω − d = 2 π × .
177 MHz and ω + d = 2 π × .
410 MHz withamplitude modulation at ω p = 2 π × V p = 1 .
5V (b), V p = 2 . V p = 2 . − ) and (+) settle, more importantlythe bifurcation points are the same. After a limit-cycleregion, both display identical period-doubling route tochaos structure confirmed by the LLE computed for eachdiagram (see fig. 4a).We plot the phase portraits showing the eigenmodesnormalized amplitudes relative dynamics in figs. 4b to 4d(top) for three singular dynamical regimes. The normal-ization is meant to get rid of the unbalanced amplitudesstill present despite the earlier discussed strengths adap-tation. We show ( R ± B − (cid:104) R ± B (cid:105) ) /σ ± with (cid:104) R ± B (cid:105) and σ ± re-spectively the mean value and the standard deviation of R ± B ( t ) calculated over the entire time trace. The dashedblack lines correspond to the synchronization regimewhere both normalized amplitudes are equal. Below theperiod-doubling bifurcation, for V p < V PD = 2 . V p > V PD the responses are now driven ina high-order synchronization regime. Nevertheless theamplitudes are clearly correlated to each other. This iseven more manifest in the chaotic regime where the am-plitudes are still correlated despite their asynchrone be-haviour with the drive. This regime corresponds to thechaotic synchronization of the nonlinearly coupled eigen-modes.We now focus on the phase responses correlationsshown in figs. 4a and 4c (bottom). In each case, we fit thedata with a unit-slope line (black-dashed) correspondingto the synchronization regime ddt ( θ + B − θ − B ) = 0. Theseplots show a tendency of synchrone evolutions of θ − B and θ + B under force modulation for V p < V PD . Contrary tothe amplitudes, the synchronization of the phases is notmaintained for V p > V PD as the trajectory does not onlylie on a unit-slope line.By studying the real-time dynamics of the phase dif-ference (see fig. 5a), we find that 2 π phase slips occurwhen V p > V PD while the resonators are phase synchro-nized (PS) for V p < V PD (green trace). When high-ordersynchronization is established between each mode andthe drive [1], phase slips resulting from phase desynchro-nization (PDS) can come up even if the amplitudes staycorrelated. This process leads to phase slips occurringregularly (red trace) – in this situation, the phases pe-riodically execute one more (or one less) cycle regard-ing the drive – or chaotically (blue trace). In the lattercase, the resonator phases stay synchronized over severalmodulation periods and this regime is interrupted by oc-casional phase slips. This corresponds to the imperfectphase synchronization (IPS) scenario [71].The different synchronization regimes can be describedthrough a statistical study of the durations between twosuccessive phase slips. For a given time trace, we list allthe PS durations τ and calculate both the mean value (cid:104) τ (cid:105) as well as the standard deviation σ τ . In fig. 5b, weplot the scaled mean PS duration (cid:104) τ (cid:105) /σ τ as a function of V p . No value can be estimated below the bifurcation tochaotic regime at V p = 2 . π/ω p are ignoredbecause it can not be qualified as synchronization.The PDS regime is identified by the regularity of the - 2 p p - 4 p
0- 2 4 p
0- 4 8 p p = 1 . 5 V q +B- q -B t i m e [ 2 p / w p ]V p = 2 . 4 2 5 VV p = 2 . 7 5 V P D S I P Sc )b ) < t > / st V p [ V ]P S 0 1 2 3 4 5 6 7 8
5V (green), phasedesynchronization (PDS) at V p = 2 . V p = 2 . (cid:104) τ (cid:105) /σ τ isplotted as a function of V p . c) Experimental probabilitydistributions of the PS durations within the PDSregime (red) or the IPS regime (blue). Distributiongiven by numerical simulation in the chaotic regime(orange dots). The exponential distribution (blackdashed) is shown for comparison.phase slips which implies that the standard deviation ofthe PS durations is near zero and leads to a peak inthe scaled mean PS duration that can be seen around V p = 2 .
425 V. The traces corresponding to this situ-ation (included in the red stripe) are used to built anhistogram in fig. 5c (left) showing the distribution of thePS durations probabilities with the associated 95% con-fidence interval. In order to compare the statistics of τ between the different traces, i.e. for different V p , we nor-malize all the durations found in a given trace by themean duration value for this trace. The probability dis-tribution is concentrated around 1, meaning that all thephase slips have almost equal duration. The correspond-ing mean duration is (cid:104) τ (cid:105) = 3 . τ / (cid:104) τ (cid:105) = 1 .
2) because this PDS occurs while thesystems sets in a period-4 motion dynamics.In the chaotic regime (blue stripe in fig. 5b) we findthat the scaled mean PS duration remains constant andslightly over 1 which tends to indicate an exponential de-cay of the probability distribution of this quantity shownin fig. 5c (right). Note that, due to the normalization bythe mean PS duration, the exponential distribution hasa unique representation (black dashed curve) in this his-togram. In this regime the mean PS duration is (cid:104) τ (cid:105) = 26modulation periods. The probability indeed decays ex-ponentially but we find that the probability around themean PS duration ( τ / (cid:104) τ (cid:105) = 1) is significantly higherthan predicted with this distribution. Additionally theobserved long PS durations occurrences are more un-likely. We conclude that the phase slips constitute anon-Poissonian process due to the deterministic chaoticdynamics and do not result from noise. The experimen-tal histogram is reproduced using a numerical simulationrealized from the new non-autonomous system of equa-tions. We recover a similar bifurcation diagram with arobust scaling of the modulated force as shown in ap-pendix E. We observe 2 π slips of the phase differencewhen the system dynamics is chaotic. From these simu-lations, we reproduce an histogram of the PS durations infig. 5c integrated over a range of modulation amplitudeshowing a chaotic domain. We find a good agreementwith the experimental distribution. Further numericaladjustments of the two driving frequencies in a restrictedrange evidence the possibility to achieve perfect phasesynchronization in the chaotic regime. This could be ofa major interest for applications based on chaos such assynchronized random number at two distinct carrier fre-quencies VI. CONCLUSION
We have analyzed the chaotic and synchronization dy-namics of distinct mechanically coupled optomechanicalresonators under slowly modulated near resonant drives.By setting the driving frequency at a bistability curveturning point, we show how a modulation of the driv-ing force amplitude leads to a period-doubling cascaderoute to chaos. Both resonators display a chaotic dy-namics, though the resonators are not synchronized. Theexperimentally-built bifurcation diagrams, with a directmeasurement of all the dynamical variables, are numeri-cally reproduced using a calibrated model of coupled non-identical Duffing oscillators. Because of the resonatornonlinear behavior modeled by a Duffing nonlinearity,we expect the normal modes not to be orthogonal andtherefore to exchange energy. This nonlinear couplingleads to amplitude synchronization and imperfect phasesynchronization as demonstrated by driving both normal modes and recording all the dynamical variables simulta-neously. Within the chaotic regime, the amplitudes arelocked to each other showing strong correlations. Thenormal modes phases dynamics are also investigated andphase synchronization, phase desynchronization and im-perfect phase synchronization regimes are observed. Inthis last case, we perform a statistical study of the syn-chronization durations and the resulting non-exponentialdistribution is confirmed by our theoretical descriptionattesting the deterministic nature of the dynamics. Here,we have deeply investigated the behavior of mechanicallycoupled optomechanical cavities. Strongly backed bysimulations, we show experimental evidence of amplitudesynchronization and intermittent phase synchronizationof bichromatic chaotic signals opening the path towardcomplex dynamics study in arrays and novel applicationssuch as multispectral synchronized random number gen-eration.
ACKNOWLEDGMENTS
This work is supported by the French RENATECHnetwork, the European Union’s Horizon 2020 researchinnovation program under grant agreement No 732894(FET Proactive HOT), the Agence Nationale de laRecherche as part of the “Investissements d’Avenir” pro-gram (Labex NanoSaclay, ANR-10-LABX-0035) with theflagship project CONDOR and the JCJC project ADOR(ANR-19-CE24-0011-01). We would like to acknowledgeMarcel Clerc for fruitful discussions.
Appendix A: DISPLACEMENT CALIBRATION
A given mechanical displacement of the membrane R A,B converts to a measured voltage ηR A,B with theconversion constant η which can be extracted thanks toa Michelson interferometer. For that purpose, an opened-loop local oscillator is used, whose optical path differencewith the sample arm is set to (cid:96) = λ/
4. Then the reso-nant oscillations amplitude of membrane B induced bya given driving strength, which translates to a voltagevariation δV is compared to the interferometric signalvariation resulting from a small calibrated displacement δ(cid:96) of the path difference. Doing this measurement forseveral excitation stage allows a confident estimation ofthe displacement calibration. We estimate the transduc-tion constant to be η ≈ . − . Appendix B: MECHANICAL MODESIDENTIFICATION
After successfully fitting membrane B response with amodel of coupled harmonic oscillators (fig. 2b), we canconclude that the observed eigenmodes frequency differ-ence is essentially caused by the natural frequency mis-0 w B V d c [ V ] w - / 2 p ( e x p ) w + / 2 p ( f i t ) w + / 2 p ( e x p ) w - / 2 p ( f i t ) u n c o u p l e d eigenfrequencies [MHz] V d c [ V ] w A b ) dw– [kHz] a ) Fig. 6: a) Measurement of the eigenfrequencies understatic voltage V dc applied on membrane B (symbols).The fit consider a parabolic shift of the self-coupledfrequency ω B (black dashed) that we use to solve theeigenmodes ( − ) and (+) (resp. blue and green lines)resulting from the previously measured coupling G. Anavoided crossing is predicted at V dc = 74V when thenatural frequencies are equal (at arrow position). b)identical experimental data presented in terms offrequencies displacement δω ± = ω ± ( V dc ) − ω ± ( V dc = 0).match ω B − ω A ≈ π ×
158 kHz. This implies that theeigenmodes ( − ) and (+) are respectively dominated bythe motion of resonators A and B. In order the verifythis conclusion, we apply a static voltage to the mem-brane B up to 25V and observe how the eigenfrequenciesare affected. It is well known that a static voltage actson a micromechanical resonator as an additional resid-ual stress [63]. This effect can be used to tune the fre-quency of a oscillator and demonstrate strong-couplingbetween several resonators through the observation of anavoided crossing in the mechanical spectrum. We obtainthe eigenfrequencies position by sweeping the drive fre-quency and measuring the response spectrum. A driveamplitude such that V ac < .
1V is set to ensure a linearresponse and avoid a possible confusion with any effectof the Duffing nonlinearity. The frequency displacementswe observe (fig. 6a) are not large enough to observe anavoided crossing. We compare the frequencies displace-ments δω ± = ω ± ( V dc ) − ω ± ( V dc = 0) in fig. 6b and itappears clearly that the eigenmode (+) is mostly affectedby the applied static voltage. In order to fit the data, weuse the Jacobian [58] of the linear system: J = − ω A G − Γ A G − ω B − Γ B The eigenvalues imaginary parts of J correspond tothe eigenfrequencies ω − and ω + while their real partscorrespond to the eigenmodes damping rates. We in-put the V dc -dependent self-coupled frequency ω B ( V dc ) = ω B (0) + αV dc . We assume the coupling G/ω B ≈ π ×
130 kHz and use the unchanged self-coupled frequency ω B (0)and the coefficient α as the fitting parameters. We findwith this second method the frequency mismatch to be ω B (0) − ω A (0) ≈ π ×
168 kHz which confirms our firstresult. The resulting eigenmodes are shown with coloredlines and the self-coupled frequencies with black dashedlines in fig. 6a. It appears that ω A ≈ ω − and ω B ≈ ω + in the range of V dc experimentally checked. The avoidedcrossing is predicted when the natural frequencies areequal, i.e. around V dc = 74V, which is out of reach in ourexperiment. Note that the assumption on the couplingvalue is only necessary to represent the level repulsionone would obtain with this voltage but does not affectour conclusion that the eigenmodes ( − ) and (+) are re-spectively dominated by the motions of membrane A andB. Appendix C: DUFFING-DUFFING MODEL
This appendix describes the derivations of the Duffing-Duffing model shown in eq. (5). We form the masterequations eq. (3) with the nonlinear potential energy: (cid:40) ¨ x A + Γ A ˙ x A + ω A x A + βx A − Gx B =0¨ x B + Γ B ˙ x B + ω B x B + βx B − Gx A = F B cos( ω d t )We start from the ansatz x A = v A cos( ω d t ) + w A sin( ω d t ) and x B = v B cos( ω d t ) + w B sin( ω d t ) wherethe quadratures relate to the response amplitude andphase with r A = v A + w A and θ A = atan2( w A , v A ) (andsimilarly for r B and θ B ). Before injecting the ansatz inthe master equations, we find useful to preliminary calcu-late x A in order to reveal the off-resonant terms oscillat-ing at 3 ω d that we can neglect in the following. We alsoexpand the expression for the derivatives ˙ x A and ¨ x A . Weneglect the quadratures second derivatives by assuming¨ v A , ¨ w A (cid:28) ω d v A , ω d w A .We now inject these preliminary results in the masterequations. We use the normalized quantities ω d t → τ ,( ω B − ω d ) /ω d → δ , ( ω B − ω A ) /ω d → ∆ ω , Γ A,B /ω d → γ A,B , G/ω d → g , F B /ω d → f B and β/ω d → ˜ β . Notethat ω B − ω d ≈ ω d δ and ω A − ω d ≈ ω d ( δ − ∆ ω ). Itwrites:1Fig. 7: Experimental bifurcation diagrams under single driving and by reading the displacement of membrane A.The measurement are performed by driving either the symmetrical (a) ω d = 2 π × .
164 MHz) or anti-symmetricalresonance (b) ω d = 2 π × .
379 MHz) with V dc = 2V, V ac =3V and ω p = 2 π × v A = 12 w A (cid:20) δ − ∆ ω ) + 34 ˜ β ( v A + w A ) (cid:21) − γ A v A − gw B (S.1.1)˙ w A = − v A (cid:20) δ − ∆ ω ) + 34 ˜ β ( v A + w A ) (cid:21) − γ A w A + 12 gv B (S.1.2)˙ v B = 12 w B (cid:20) δ + 34 ˜ β ( v B + w B ) (cid:21) − γ B v B − gw A (S.2.1)˙ w B = − v B (cid:20) δ + 34 ˜ β ( v B + w B ) (cid:21) − γ B w B + 12 gv A + 12 ˜ f B (S.2.2)This system of equations describes the evolution ofour system in terms of quadratures v A , w A , v B and w B . It presents an interest for the numerical simulationssince these quantities are homogeneous. However we canobtain the final set of equations (eq. (5)) expressed interms of the amplitudes and phases by performing the re-verse transformations v A = r A cos( θ A ), w A = r A sin( θ A ), v B = r B cos( θ B ) and w B = r B sin( θ B ) and taking: S. . × cos( θ A ) + S. . × sin( θ A ) S. . × cos( θ B ) + S. . × sin( θ B ) − S. . × sin( θ A ) + S. . × cos( θ A ) − S. . × sin( θ B ) + S. . × cos( θ B )The numerical simulations are performed by integrat-ing the ODEs with the physical quantities found in theexperiments. For a given resonance, the driving fre-quency is adjusted near the jump-up frequency of thebistability for an optimal scaling of the bifurcation struc-ture. The diagrams presented in figs. 3c and 3f are foundfor ω d = 2 π × . ω d = 2 π × . f B = m eff ω d dCdx V dc V ac with V dc = 2 V and V ac = 3V. To ac-count for the dephasing between the experimental driv-ing excitation and the force, we present the simulatedbifurcation diagram in fig. 3f with an offset of 80 ◦ deg in θ B . Appendix D: BIFURCATION DIAGRAMS BUILTFROM THE READING OF MEMBRANE A
The bifurcation diagrams presented in figs. 3a and 3drespectively show the dynamics of the normal modes ( − )and (+) with parameter V p . They are based on the re-sponse of membrane B. Here,we perform the same mea-surement by placing the laser on membrane A to readits displacement. In fig. 7a we show the bifurcation dia-gram made from the response of membrane A when themode ( − ) is submitted to a modulated force by apply-ing V dc = 2V, V ac = 3 V at frequency ω d = 2 π × ω p = 7 kHz. The associated LLE is shownbelow. Similarly in fig. 7b we show the measured dia-gram corresponding to the case where the normal mode(+) is driven with the same parameters but at resonantfrequency ω d = 2 π × − ) and (+)are almost orthogonal and that measuring their responsethrough the membrane A or B gives the same result be-side the amplitude imbalance.Although the bifurcation diagrams are very similarwhether A or B is read, a small shift in the bifurcationpoints positions can be observed and even a regime of pe-riodic oscillations is present around V p =2.3V in fig. 7bthat is not present in fig. 3d. This is a consequence of thephotothermal shift induced by the laser on the eigenfre-quency dominated by the probed membrane [39]. Thisshift is lower than 3 kHz but leads to a significant modifi-cation of the bifurcation diagram. When driving a givennormal mode, we expect the membrane responses to be2perfectly correlated. The normal modes result from thestrong coupling interaction between the membranes andthe fact that they both are identically affected by thedynamics of a normal mode should not be understood assynchronization. This can not be confirmed without asimultaneous lecture of both membranes although it wascorroborated by our numerical simulations. Appendix E: TWO-DRIVES MODEL ANDORTHOGONALITY BREAKING
Driving both normal modes at the same time allowsto evidence synchronization phenomena. We model thesystem with the same master equation but add a secondresonant excitation: ¨ x A + Γ A ˙ x A + ω A (1 + ˜ βx A ) x A − Gx B =0¨ x B + Γ B ˙ x B + ω B (1 + ˜ βx B ) x B − Gx A = F − B cos( ω − d t )+ F + B cos( ω + d t )Since the system is now expected to respond both at ω − d and ω + d , we modify the ansatz: x A = v − A cos( ω − d t ) + w − A sin( ω − d t )+ v + A cos( ω + d t ) + w + A ( ω + d t )The rest of the calculations is essentially the same,except for the development of the cubic terms x A and x B where the nonlinear coupling between r − A,B and r + A,B comes from. We neglect all off-resonant terms includingthe ones oscillating at 2 ω ± d − ω ∓ d . Following the exactsame procedure as in appendix C, we derive a system of8 coupled nonlinear ODEs for the normal modes ( − ) and(+) quadratures accessed either through the membraneA ( v − A , w − A , v + A , w + A ) or B ( v − B , w − B , v + B , w + B ). It writes:˙ v ± A = 12 w ± A (cid:20) ε ± ( δ ± − ∆ ω ) + 34 ε ± ˜ β ( r ± A + 2 r ∓ A ) (cid:21) − ε ± γ A v ± A − ε ± gw ± B ˙ w ± A = − v ± A (cid:20) ε ± ( δ ± − ∆ ω ) + 34 ε ± ˜ β ( r ± A + 2 r ∓ A ) (cid:21) − ε ± γ A w ± A + 12 ε ± gv ± B ˙ v ± B = 12 w ± B (cid:20) δ ± + 34 ε ± ˜ β ( r ± B + 2 r ∓ B ) (cid:21) − ε ± γ B v ± B − ε ± gw ± A ˙ w ± B = − v ± B (cid:20) δ ± + 34 ε ± ˜ β ( r ± B + 2 r ∓ B ) (cid:21) − ε ± γ B w ± B + 12 ε ± gv ± A + ε ± f ± B where we use the amplitudes r ± A = v ± A + w ± A and r ± B = v ± B + w ± B for the compactness of these expres-sions. The parameter ε ± = ω + d /ω ± d aims to regularize the normalization of the parameters such that the timeis arbitrarily chosen to be rescaled with ω + d t . Thus allthe physical are normalized consequently to this choice: δ ± = ( ω B − ω ± d ) /ω + d , ∆ ω = ( ω B − ω A ) /ω + d , γ A,B =Γ A,B /ω + d , g = G/ω + d , ˜ β = β/ω + d and f ± B = f ± B /ω + d .These two systems of 4 equations contains new termsthat allow the normal modes response amplitudes r + A,B and r − A,B to couple. By expressing these equations interms of amplitudes ( r ± A,B ) and phases ( θ ± A,B ), the nor-mal mode coupling takes the form ˜ βr ± r ∓ .Fig. 8: Simulated bifurcation diagrams under twodriving forces with parameter V p . The top diagram(resp. the bottom diagram) uses the Poincar´e sectionmade from the maxima of w − B (resp. the maxima of w + B ). We use the driving frequencies ω − d = 2 π × . ω + d = 2 π × . . This simulation should be compared with the experimentaldiagrams presented in fig. 4a. We integrate these equations using the parameters V dc = 2 V, V − ac = 3 . V + ac = 0 . ω − d = 2 π × . ω + d = 2 π × . ω p = 2 π × V p , a bifurcation diagram is reproducedwith the amplitude response of both normal modes us-ing the maxima of w − B and w + B . The phases θ − B and θ + B are both shifted by 100 ◦ deg thus compensating an ex-perimental dephasing. We recover the period-doublingstructure followed by chaotic intermittency that is ex-perimentally obtained (see fig. 4a). The period doublingoccurs for V p = 1 .
840 V which is slightly below the valuefound experimentally V PD = 2 .
140 V. This can be ex-plained by a limitation in the several experimental cal-ibrations we have performed as well as a drift of someof the mechanical parameters with time. The amplitude3of the local maxima in each simulated diagram is verysimilar to the corresponding experimental diagrams. Fi-nally, the dynamical regimes in which the normal modessettle are perfectly correlated, as observed in the exper-iments. We analyze a thousand chaotic traces between V p = 2 .
742 V and V p = 3 V and extract more than 22000phase slip occurrences on which the statistical study ofthe phase synchronization durations presented in fig. 4cis performed. [1] A. Pikovsky, M. G. Rosenblum, and J. Kurths,Synchronization, A Universal Concept in Nonlinear Sciences(Cambridge University Press, Cambridge, 2001).[2] C. Huygens, Letters to de sluse, Soci´et´e Hollandaise DesSciences, Martinus Nijho, 1895 , letters; no. 1333 of 24February 1665, no. 1335 of 26 February 1665, no. 1345 of6 March 1665 (1665).[3] M. G. Rosenblum, A. S. Pikovsky, and J. Kurths, Phasesynchronization in driven and coupled chaotic oscillators,IEEE Transactions on Circuits and Systems I: Funda-mental Theory and Applications , 874 (1997).[4] A. N. Pisarchik and R. Jaimes-Re´ategui, Intermittent lagsynchronization in a driven system of coupled oscillators,Pramana , 503 (2005).[5] M. Zhang, G. S. Wiederhecker, S. Manipatruni,A. Barnard, P. McEuen, and M. Lipson, Synchronizationof micromechanical oscillators using light, Phys. Rev.Lett. , 233906 (2012).[6] M. Zhang, S. Shah, J. Cardenas, and M. Lipson, Syn-chronization and phase noise reduction in micromechan-ical oscillator arrays coupled through light, Phys. Rev.Lett. , 163902 (2015).[7] Y. Kuramoto and T. Yamada, Pattern formation in oscil-latory chemical reactions, Progress of Theoretical Physics , 724 (1976).[8] L. Glass, Synchronization and rhythmic processes inphysiology, Nature , 277 (2001).[9] B. Blasius, A. Huppert, and L. Stone, Complex dynamicsand phase synchronization in spatially extended ecologi-cal systems, Nature , 354 (1999).[10] C. Volos, I. Kyprianidis, and I. Stouboulos, Synchroniza-tion phenomena in coupled nonlinear systems applied ineconomic cycles, WSEAS Transactions on Systems ,681 (2012).[11] M. Moussaid, S. Garnier, G. Theraulaz, andD. Helbing, Collective information processingand pattern formation in swarms, flocks, andcrowds, Topics in Cognitive Science , 469 (2009),https://onlinelibrary.wiley.com/doi/pdf/10.1111/j.1756-8765.2009.01028.x.[12] L. Pecora and T. Caroll, Synchronization in chaoticsystems, Physical Review Letters , 10.1103/Phys-RevLett.64.821 (1990).[13] G. S. Duane and J. J. Tribbia, Synchronized chaos ingeophysical fluid dynamics, Phys. Rev. Lett. , 4298(2001).[14] G. S. Duane and J. J. Tribbia, Weak atlantic-pacific teleconnections as synchronized chaos,Journal of the Atmospheric Sciences ,2149 (2004), https://doi.org/10.1175/1520-0469(2004)061¡2149:WATASC¿2.0.CO;2.[15] P. H. Hiemstra, N. Fujiwara, F. M. Selten, and J. Kurths,Complete synchronization of chaotic atmospheric modelsby connecting only a subset of state space, Nonlinear Processes in Geophysics , 611 (2012).[16] F. Lunkeit, Synchronization experiments with an atmo-spheric global circulation model, Chaos: An Interdis-ciplinary Journal of Nonlinear Science , 47 (2001),https://aip.scitation.org/doi/pdf/10.1063/1.1338127.[17] A. Patra, B. L. Altshuler, and E. A. Yuzbashyan, Chaoticsynchronization between atomic clocks, Phys. Rev. A , 023418 (2019).[18] M. Sciamanna and K. Shore, Physics and applications oflaser diode chaos, Nature Photonics , 151 (2015).[19] K. M. Cuomo, A. V. Oppenheim, and S. H. Strogatz,Synchronization of lorenz-based chaotic circuits with ap-plications to communications, IEEE Transactions on Cir-cuits and Systems II: Analog and Digital Signal Process-ing , 626 (1993).[20] A. Argyris, D. Syvridis, L. Larger, V. Annovazzi-Lodi,P. Colet, I. Fischer, J. Garc´ıa-Ojalvo, C. R. Mirasso,L. Pesquera, and K. A. Shore, Chaos-based communica-tions at high bit rates using commercial fibre-optic links,Nature , 343 (2005).[21] V. Annovazzi-Lodi, S. Donati, and A. Scire, Synchroniza-tion of chaotic injected-laser systems and its applicationto optical cryptography, IEEE Journal of Quantum Elec-tronics , 953 (1996).[22] C. R. Mirasso, P. Colet, and P. Garcia-Fernandez, Syn-chronization of chaotic semiconductor lasers: applicationto encoded communications, IEEE Photonics TechnologyLetters , 299 (1996).[23] C. M. Marinho, E. E. Macau, and T. Yoneyama, Chaosover chaos: A new approach for satellite communication,Acta Astronautica , 230 (2005), infinite PossibilitiesGlobal Realities, Selected Proceedings of the 55th Inter-national Astronautical Federation Congress, Vancouver,Canada, 4-8 October 2004.[24] C. C. Hsu, D. Gobovic, M. E. Zaghloul, and H. H. Szu,Chaotic neuron models and their vlsi circuit implemen-tations, IEEE Transactions on Neural Networks , 1339(1996).[25] V. Milanovic and M. E. Zaghloul, Synchronization ofchaotic neural networks for secure communications, in1996 IEEE International Symposium on Circuits and Systems. Circuits and Systems Connecting the World. ISCAS 96,Vol. 3 (1996) pp. 28–31 vol.3.[26] L. J. Fiderer and D. Braun, Quantum metrology withquantum-chaotic sensors, Nature Communications ,1351 (2018).[27] L. Bakemeier, A. Alvermann, and H. Fehske, Route tochaos in optomechanics, Phys. Rev. Lett. , 013601(2015).[28] E.-H. Park, M. A. Zaks, and J. Kurths, Phase synchro-nization in the forced lorenz system, Phys. Rev. E ,6627 (1999).[29] R. Lifshitz and M. Cross, Response of parametricallydriven nonlinear coupled oscillators with application tomicromechanical and nanomechanical resonator arrays, Physical Review B , 10.1103/PhysRevB.67.134302(2003).[30] J. A. Blackburn, G. L. Baker, and H. J. T. Smith, In-termittent synchronization of resistively coupled chaoticjosephson junctions, Phys. Rev. B , 5931 (2000).[31] C. Wu and L. Chua, Synchronization in an array of lin-early coupled dynamical systems, IEEE Transactions onCircuits and Systems I: Fundamental Theory and Appli-cations , 10.1109/81.404047 (1995).[32] J.-W. Shuai and D. M. Durand, Phase synchronizationin two coupled chaotic neurons, Physics Letters A ,289 (1999).[33] M. G. Clerc, S. Coulibaly, M. A. Ferr´e, and R. G.Rojas, Chimera states in a duffing oscillators chaincoupled to nearest neighbors, Chaos: An Interdisci-plinary Journal of Nonlinear Science , 083126 (2018),https://doi.org/10.1063/1.5025038.[34] E. A. Martens, S. Thutupalli, A. Fourri`ere,and O. Hallatschek, Chimera states in mechan-ical oscillator networks, Proceedings of the Na-tional Academy of Sciences ,012902 (2015).[37] E. Akhmedov and A. Mirizzi, Another look at synchro-nized neutrino oscillations, Nuclear Physics B , 382(2016).[38] L. Midolo, A. Schliesser, and A. Fiore, Nano-opto-electro-mechanical systems, Nature Nanotechnology ,11 (2018).[39] N. Gao, W. Luo, W. Yan, D. Zhang, and D. Liu, Contin-uously tuning the resonant characteristics of microcan-tilevers by a laser induced photothermal effect, Journalof Physics D: Applied Physics , 385402 (2019).[40] Q. P. Unterreithmeier, E. M. Weig, and J. P. Kotthaus,Universal transduction scheme for nanomechanical sys-tems based on dielectric forces, Nature , 1001 (2009).[41] A. Eichler, M. del ´Alamo Ruiz, J. A. Plaza, and A. Bach-told, Strong coupling between mechanical modes in ananotube resonator, Phys. Rev. Lett. , 025503 (2012).[42] K. Gajo, S. Sch¨uz, and E. M. Weig, Strong 4-mode coupling of nanomechanical string resonators,Applied Physics Letters , 133109 (2017),https://doi.org/10.1063/1.4995230.[43] H. Okamoto, A. Gourgout, C.-Y. Chang, K. Onomitsu,I. Mahboob, E. Y. Chang, and H. Yamaguchi, Coherentphonon manipulation in coupled mechanical resonators,Nature Physics , 480 (2013).[44] A. Chowdhury, S. Barbay, M. G. Clerc, I. Robert-Philip,and R. Braive, Phase stochastic resonance in a forcednanoelectromechanical membrane, Phys. Rev. Lett. ,234101 (2017).[45] S.-B. Shim, M. Imboden, and P. Mohanty,Synchronized oscillation in coupled nanome-chanical oscillators, Science , 95 (2007),https://science.sciencemag.org/content/316/5821/95.full.pdf.[46] D. Pu, X. Wei, L. Xu, Z. Jiang, and R. Huan, Synchro-nization of electrically coupled micromechanical oscilla- tors with a frequency ratio of 3:1, Applied Physics Letters , 013503 (2018), https://doi.org/10.1063/1.5000786.[47] M. H. Matheny, M. Grau, L. G. Villanueva, R. B. Kara-balin, M. C. Cross, and M. L. Roukes, Phase synchro-nization of two anharmonic nanomechanical oscillators,Phys. Rev. Lett. , 014101 (2014).[48] R. B. Karabalin, M. C. Cross, and M. L. Roukes, Nonlin-ear dynamics and chaos in two coupled nanomechanicalresonators, Phys. Rev. B , 165309 (2009).[49] J. Ma, C. You, L.-G. Si, H. Xiong, J. Li, X. Yang, andY. Wu, Formation and manipulation of optomechanicalchaos via a bichromatic driving, Phys. Rev. A , 043839(2014).[50] J. Larson and M. Horsdal, Photonic josephson effect,phase transitions, and chaos in optomechanical systems,Phys. Rev. A , 021804 (2011).[51] L. Jin, Y. Guo, X. Ji, and L. Li, Reconfigurable chaosin electro-optomechanical system with negative duffingresonators, Scientific Reports , 4822 (2017).[52] D. Navarro-Urrios, N. E. Capuj, M. F. Colombano, P. D.Garc´ıa, M. Sledzinska, F. Alzina, A. Griol, A. Mart´ınez,and C. M. Sotomayor-Torres, Nonlinear dynamics andchaos in an optomechanical beam, Nature Communica-tions , 14965 (2017).[53] J. Wu, S.-W. Huang, Y. Huang, H. Zhou, J. Yang, J.-M. Liu, M. Yu, G. Lo, D.-L. Kwong, S. Duan, andC. Wei Wong, Mesoscopic chaos mediated by drudeelectron-hole plasma in silicon optomechanical oscilla-tors, Nature Communications , 15570 (2017).[54] A. Roy and M. Devoret, Introduction to parametric am-plification of quantum signals with josephson circuits,Comptes Rendus Physique , 740 (2016), quantum mi-crowaves / Micro-ondes quantiques.[55] H. C. S. Hsuan, R. C. Ajmera, and K. E. Lonngren, Thenonlinear effect of altering the zeroth order density dis-tribution of a plasma, Applied Physics Letters , 277(1967), https://doi.org/10.1063/1.1755133.[56] T. Antoni, A. G. Kuhn, T. Briant, P.-F. Cohadon,A. Heidmann, R. Braive, A. Beveratos, I. Abram, L. L.Gratiet, I. Sagnes, and I. Robert-Philip, Deformable two-dimensional photonic crystal slab for cavity optomechan-ics, Opt. Lett. , 3434 (2011).[57] A. Chowdhury, I. Yeo, V. Tsvirkun, F. Raineri, G. Beau-doin, I. Sagnes, R. Raj, I. Robert-Philip, and R. Braive,Superharmonic resonances in a two-dimensional non-linear photonic-crystal nano-electro-mechanical oscil-lator, Applied Physics Letters , 163102 (2016),https://doi.org/10.1063/1.4947064.[58] D. H. Zanette, Energy exchange between coupled me-chanical oscillators: linear regimes, Journal of PhysicsCommunications , 095015 (2018).[59] Y. S. Joe, A. M. Satanin, and C. S. Kim, Classical anal-ogy of fano resonances, Physica Scripta , 259 (2006).[60] M. F. Limonov, M. V. Rybin, A. N. Poddubny, and Y. S.Kivshar, Fano resonances in photonics, Nature Photonics , 543 (2017).[61] S. Stassi, A. Chiad`o, G. Calafiore, G. Palmara,S. Cabrini, and C. Ricciardi, Experimental evidence offano resonances in nanomechanical resonators, ScientificReports , 1065 (2017).[62] S. Zanotto, Weak coupling, strong coupling, criticalcoupling and fano resonances: A unifying vision, inFano Resonances in Optics and Microwaves: Physics and Applications,edited by E. Kamenetskii, A. Sadreev, and A. Mirosh- nichenko (Springer International Publishing, Cham,2018) pp. 551–570.[63] J. Rieger, T. Faust, M. J. Seitner, J. P. Kotthaus, andE. M. Weig, Frequency and q factor control of nanome-chanical resonators, Applied Physics Letters , 103110(2012), https://doi.org/10.1063/1.4751351.[64] J. F. Rhoads, S. W. Shaw, and K. L. Turner, Nonlineardynamics and its applications in micro- and nanores-onators, Journal of Dynamic Systems, Measurement,and Control , 10.1115/1.4001333 (2010), 034001,https://asmedigitalcollection.asme.org/dynamicsystems/article-pdf/132/3/034001/5653925/034001 1.pdf.[65] A. Chowdhury, M. G. Clerc, S. Barbay, I. Robert-Philip,and R. Braive, Weak signal enhancement by non-linearresonance control in a forced nano-electromechanical res-onator (2019), arXiv:1910.04686 [physics.app-ph].[66] R. Hegger, H. Kantz, and T. Schreiber, Practi-cal implementation of nonlinear time series meth-ods: The tisean package, Chaos: An Interdisci-plinary Journal of Nonlinear Science , 413 (1999),https://doi.org/10.1063/1.166424.[67] M. T. Rosenstein, J. J. Collins, and C. J. D. Luca, A prac- tical method for calculating largest lyapunov exponentsfrom small data sets, Physica D: Nonlinear Phenomena , 117 (1993).[68] C. Lee, T. Yoon, and S. Shin, Period doubling and chaosin a directly modulated laser diode, Applied Physics Let-ters , 10.1063/1.95810 (1985).[69] D. Cadeddu, F. R. Braakman, G. T¨ut¨unc¨uoglu,F. Matteini, D. R¨uffer, A. Fontcuberta i Morral,and M. Poggio, Time-resolved nonlinear coupling be-tween orthogonal flexural modes of a pristine gaasnanowire, Nano Letters , 926 (2016), pMID: 26785132,https://doi.org/10.1021/acs.nanolett.5b03822.[70] L. Mercier de L´epinay, B. Pigeau, B. Besga, and O. Ar-cizet, Eigenmode orthogonality breaking and anomalousdynamics in multimode nano-optomechanical systemsunder non-reciprocal coupling, Nature Communications , 1401 (2018).[71] S. Boccaletti, J. Kurths, G. Osipov, D. Valladares, andC. Zhou, The synchronization of chaotic systems, PhysicsReports366