Nonlinear dynamics of weakly dissipative optomechanical systems

Optomechanical systems attract a lot of attention because they provide a novel platform for quantum measurements, transduction, hybrid systems, and fundamental studies of quantum physics. Their classical nonlinear dynamics is surprisingly rich and so far remains underexplored. Works devoted to this subject have typically focussed on dissipation constants which are substantially larger than those encountered in current experiments, such that the nonlinear dynamics of weakly dissipative optomechanical systems is almost uncharted waters. In this work, we fill this gap and investigate the regular and chaotic dynamics in this important regime. To analyze the dynamical attractors, we have extended the"Generalized Alignment Index"method to dissipative systems. We show that, even when chaotic motion is absent, the dynamics in the weakly dissipative regime is extremely sensitive to initial conditions. We argue that reducing dissipation allows chaotic dynamics to appear at a substantially smaller driving strength and enables various routes to chaos. We identify three generic features in weakly dissipative classical optomechanical nonlinear dynamics: the Neimark-Sacker bifurcation between limit cycles and limit tori (leading to a comb of sidebands in the spectrum), the quasiperiodic route to chaos, and the existence of transient chaos.


Introduction
Cavity optomechanics [1] aims to explore and exploit the interaction between radiation fields and mechanical vibrations, with important applications ranging from sensitive measurements to quantum communication. The foundations for this research field were established already at the end of the 60s, when the classical effects of radiation on the motion of a test mass were studied in the context of precision measurements [2,3]. For an extended review we refer the reader to Ref. [1]. In the past few years, a range of impressive achievements has been observed, which includes topological transport in optomechanical arrays [4,5], the engineering of nonreciprocal interactions [6][7][8][9][10][11], the generation of single phonon states using optical control [12], the generation of mechanical squeezed states [13], measurement-based quantum control of mechanical motion [14], conversion of quantum information to mechanical motion [15], conversion between light in the microwave and optical range [16], single photon frequency shifters [17], force measurements using cold-atom optomechanics [18], and the use of unconventional mechanical modes, like high frequency bulk modes of crystals [19], multilayer graphene [20], and the modes of superfluid helium [21].
Classical nonlinear optomechanics is relevant in the case of highly populated optical and mechanical modes. Though it attracted slightly less attention during the initial evolution of modern cavity optomechanics, a number of significant theoretical studies have been devoted to understanding the structure of the phase space, including limit cycles and multistability [22][23][24][25][26], and chaotic dynamics [27,28]. Experimental studies have been relatively rare, but important phenomena have already been observed, including limit cycles [29,30], period doubling and chaos [31][32][33][34][35][36], the predicted multistable attractor diagram [37,38] which is characteristic for optomechanical systems, as well as further aspects [39,40]. More recent studies have exploited the coupling of several OM limit cycle oscillators to explore OM synchronization dynamics. OM synchronization was first predicted theoretically in [41], then observed experimentally for few-mode systems [42][43][44][45], and analyzed in subsequent theoretical studies of large-scale lattice dynamics [46][47][48][49].
Many theoretical works on nonlinear classical OM dynamics have considered mainly systems operating outside the so-called resolved sideband regime. This means that the optical dissipation is assumed to be of the same order or larger than the mechanical frequency. At the same time, the mechanical quality factor is often assumed relatively small, of the order of O(10 3 ). For instance, the authors of Ref. [27] have shown that limit cycles in such strongly dissipative OM systems undergo a period doubling cascade and become chaotic attractors.
On the other hand, most state-of-the-art experiments reach the resolved sideband regime and deal with substantially larger mechanical quality factors, ranging from 10 4 to 10 9 (cf. Figs. 11 and 10 in Ref. [1]). These experiments raise a natural question: do such weakly dissipative systems show something qualitatively new in their classical dynamics? The straightforward guess is: yes, because nonlinear phenomena are expected to be enhanced with decreasing dissipation. For instance, the Hopf bifurcation [23], at which an equilibrium point of the dynamics becomes unstable and a limit cycle emerges, has a clear dependence on the dissipation constants. The smaller the dissipation constants, the weaker the laser pumping needed to observe the Hopf bifurcation. Bistability, which is another nonlinear phenomenon, follows the same rule. Of course, the possible types of attractors are also very sensitive to the dissipation strength. One could take one step further and ask whether the chaotic OM dynamics is enhanced as well and acquires new features in the resolved sideband regime.
In this work we investigate the nonlinear dynamics of weakly dissipative OM systems. Weakly dissipative in this context is the same as sideband resolved, meaning that the optical dissipation is much smaller than the mechanical frequency. Firstly, we are interested in performing a classification of attractors: whether they are chaotic or regular, what is their dimensionality, etc. We show that the weakly dissipative regime is much more complex and nontrivial than the strongly dissipative one. In particular, the OM dynamics becomes very sensitive to the initial conditions in the resolved sideband regime, which represents the first substantial difference between the strongly and weakly dissipative cases.
This sensitivity to initial conditions (as well as the long relaxation times) makes the study of weakly dissipative OM systems computationally very challenging. To overcome this problem, we suggest a new approach to classify the attractors and to detect dynamical chaos. It is based on the GALI (Generalized ALignment Index) method [50][51][52] and has several advantages. Besides being significantly faster than commonly used methods based on the calculation of the maximal Lyapunov exponent, the modified GALI method provides an efficient tool to learn the dimensionality of the attractors. This has allowed us to explore the OM attractors in a large range of parameters and to reveal important phenomena which are well-known in nonlinear science but have been overlooked so far in optomechanics. They include transient classical chaos, quasiperiodic orbits, and routes to chaos beyond the period doubling.
The rest of this paper is organized as follows: In Sect. 2, we introduce the equations of motion of an OM system and discuss the basic differences between the strongly and weakly dissipative regimes. Sect. 3 is devoted to the GALI method and its extension to the analysis of dissipative nonlinear dynamics. We use this method and our numerical simulations to present a diagram that illustrates various regular and chaotic weakly dissipative dynamical regimes in Sect. 4. In particular, we identify two generic features that will become important in the exploration of nonlinear optomechanics: a Neimark-Sacker bifurcation between limit cycles and limit tori (leading to a comb of sidebands in the spectrum) and the existence of transient chaos. In Sect. 5, we discuss the experimental relevance of our results. Finally, Sect. 6 contains our conclusions.
2 Classical dynamics of a weakly dissipative optomechanical system

Equations of motion
The classical dynamics of an OM system with one optical mode and one mechanical mode (sometimes referred to as the optical cavity and the mechanical oscillator, respectively) is described by the following equations of motion [1]: Here b = (q + ip)/ √ 2, with q and p being the dimensionless position and momentum of the mechanical oscillator, and a is the suitably normalized complex amplitude of the electric field inside the cavity (|a| 2 and |b| 2 are the photon and phonon number, respectively). The mechanical (optical) mode has frequency Ω m (ω c ) and dissipation constant γ (κ). The optical mode is pumped by an external laser with frequency ω L and amplitude E; ∆ = ω L − ω c denotes the detuning between the laser frequency and the cavity frequency; g 0 is the bare optomechanical coupling constant. We note that E here is normalized such that E 2 /κ is the rate of photons impinging on the cavity. The typical representation of an OM system is shown in Fig. 1. As usual, we work in a reference frame which rotates at the laser frequency 1 . Eqs. (1,2) assume that quantum fluctuations can be neglected, i.e. the dynamics is governed by highly populated optical and mechanical states. These coupled equations have been employed to describe countless experiments to high precision, both in the linearized regime but also in the fully nonlinear regime of interest here.
Further numerical study requires to rewrite Eqs. (1,2) in a dimensionless form. This can be done by defining rescaled variables α = aΩ m /2E and β = g 0 b/Ω m , from which we obtain the Figure 1: Representation of a generic OM system: an optical cavity with a movable mirror driven by an external laser.
following equations: where τ = Ω m t and P = 8g 2 0 E 2 /Ω 4 m . Note that there are fewer parameters in the rescaled Eqs. (3,4) than in the original Eqs. (1,2). This means the qualitative features of the dynamics will only depend on four dimensionless combinations of the original physical parameters: dimensionless power P , normalized detuning ∆/Ω m , normalized cavity decay κ/Ω m , and mechanical dissipation γ/Ω m . For a more extended discussion of the essential dimensionless parameters affecting classical or quantum OM dynamics, we refer the reader to Refs. [23,48].
The parameter P is a dimensionless measure of the laser input power, which also includes the strength of the optomechanical interaction. P can be related to the standard measure of coupling strength vs. dissipation, the so-called OM cooperativity C = 4g 2 0 n c /γκ. Here n c is the mean number of photons stored in the optical cavity. For our purposes the cooperativity is still slightly inconvenient, since n c depends on the detuning (at fixed drive power). For that reason, we rather introduce the maximum cooperativityC = 4g 2 0 n 0 /γκ, where n 0 = 4E 2 /κ 2 is the number of photons in the resonantly pumped optical cavity in the absence of the optomechanical interaction. P is then proportional to the maximum cooperativity as follows: This relation will be useful for comparison with experimental parameters.

Fixed points
Let us start our study of the dynamics with the analysis of the fixed points of the system. Fixed points are points in the phase space which are invariant under time evolution: if we take a fixed point as initial condition of the system, the system stays on the fixed point forever. The analysis of trajectories whose initial conditions are arbitrarily close to the fixed point allows one to classify the fixed point as stable, unstable, or hyperbolic. If any such trajectory is attracted to (repelled from) the fixed point, the fixed point is stable (unstable). If some trajectories are attracted to the fixed point, while other trajectories are repelled from it, the fixed point is hyperbolic. Stable fixed points are the simplest attractors of a dynamical system.
Although the fixed points of the OM systems have been known for a long time [1,53], it is important to understand them in more detail, because this will provide the context for the discussions of the dynamical attractors. The fixed point equations are obtained by setting the time derivatives in Eqs. (3,4) to zero and solving the resulting set of nonlinear equations: Here Q = (β + β * )/ √ 2 is the rescaled position of the mechanical oscillator. After inserting Eq. (6) into (7), we obtain a third order polynomial equation for Q with real coefficients. Since Q is also real, the system has at least one fixed point; the maximum number is obviously three [53]. Figs. 2(a,b) show the fixed point diagram for an OM system with dissipation constants κ = Ω m and γ = 10 −3 Ω m , and for an OM system in the sideband-resolved regime, κ = 10 −1 Ω m and γ = 10 −4 Ω m , respectively. Below, we will refer to these two representative cases as the "strongly dissipative" and the "weakly dissipative" OM systems, respectively. For sufficiently small P there is, as one would expect, just one stable fixed point. As the parameter P is increased, the fixed points can follow two possible scenarios with different bifurcation phenomena. A bifurcation is a qualitative change of the dynamics which occurs as a system parameter is varied [54]. For fixed points, this typically means creation or annihilation of fixed points, or change of the type of a fixed point (whether the fixed point is stable, unstable or hyperbolic). The first scenario is shown in Fig. 2 (c): an (inverse) saddle-node bifurcation 2 takes place at some value of P and a pair of stable-unstable fixed points is created. Increasing P leads to a Hopf bifurcation 3 at which the stable fixed point becomes unstable. Further increase of P results in a saddle-node bifurcation at which a pair of stable-unstable fixed points is annihilated. In some cases, the Hopf bifurcation may occur after the saddle-node bifurcation. This scenario occurs only at ∆ < 0 ("red detuning"). The second scenario is shown in Fig. 2(d): the Hopf bifurcation again occurs at some value of P and makes the stable fixed point unstable. Further increase of P does not change the nature and the number of the fixed points. Even though the above described bifurcations can be observed in both strongly and weakly dissipative OM systems, Figs. 2(a) and 2(b) clearly show the essential difference between them. When the dissipation is weaker, the bifurcations may occur at much smaller values of P and the stability diagram becomes more complex. Since these bifurcations are genuine nonlinear phenomena and P is the strength of the nonlinear interaction, Figs. 2(a,b) provide us with a first indication that nonlinear effects are more pronounced and even qualitatively altered in the weakly dissipative case.

Attractors
The Liouville's theorem guarantees that the time evolution of Hamiltonian systems preserves volumes in the phase space. In contrast to Hamiltonian systems, dissipative systems are defined as systems in which volumes shrink over time in some region of the phase space [54]. For these systems, generically speaking, the shrinking volumes collapse, in the long time limit, to the so-called attractors. An attractor has the following properties [55]: 2 In a saddle-node bifurcation a pair of stable-unstable fixed points approach each other as a parameter η is varied (for simplicity and without loss of generality, let us suppose that we are increasing η). At η = η * the two fixed points merge and form one single stable fixed point; if η > η * the fixed points cease to exist. If η is decreased one comes across the inverse saddle-node bifurcation, in which a pair of stable-unstable fixed points is created. 3 In the Hopf bifurcation a stable fixed point becomes unstable and a periodic orbit appears as a parameter η is varied. The periodic orbit can be unstable or stable. In the latter case it is called a limit cycle. The Hopf bifurcation is also known as a Poincaré-Andronov-Hopf bifurcation.   (i) It is a subset of the phase space which is invariant under the dynamics.
(ii) There must exist another (noninvariant) subset of the phase space which defines the initial conditions for the trajectories asymptotically approaching (being "attracted" by) the attractor at t → ∞. The second subset is called the basin of attraction. (iii) An attractor cannot be decomposed in two or more disjoint attractors.
The attractors of a dissipative system typically provide important information about its dynamics. In particular, we expect them to illustrate the differences between the strongly and weakly dissipative nonlinear dynamics of OM systems. As said before, a stable fixed point is the simplest kind of attractor. The Hopf bifurcation leads to the emergence of stable limit cycles, which in turn can undergo transitions to other attractors, including chaotic ones. In the strongly dissipative regime the limit cycles of a OM system undergo the well known "perioddoubling cascade" 4 at P ∼ 1, becoming chaotic attractors. This phenomenon was described theoretically in Ref. [27] and observed in early pioneering experiments [31]. In the weakly dissipative regime, however, where the fixed point analysis suggests stronger nonlinear effects, neither the attractors nor the associated routes to chaos have been studied. Below we focus on this regime.

Basins of attraction and hypersensitivity to the initial state
A nonlinear dissipative system has generally more than one attractor and its long time dynamics depends on initial conditions which can belong to one or another basin of attraction. Some attractors can be very challenging to reach both in numerical simulations and real experiments because their basin of attraction is rather small and their detection would require a nontrivial fine tuning of the initial conditions. We will address the properties of those OM attractors which are easily accessible and, therefore, relevant for experiments. Throughout this section, we focus on the weakly dissipative case.
We have simulated Eqs. (3,4) for different initial conditions of the mechanical oscillator 5 , assuming that the laser is turned on abruptly at t = 0 (thus α(0) = 0). Fig. 3 shows the observed attractors and their basins of attraction. While the strongly dissipative OM dynamics usually reveals just one attractor, the phase space of the weakly dissipative OM systems is much richer. One can observe not only several co-existing attractors, i.e. multistability, but also very complex and entangled basins of attraction, see Fig. 3(b). Fig. 3(c) shows a zoom of a small part of the basin of attraction from Fig. 3(b) (the area within the white square) with a higher resolution. One can see that, even on this scale, the basin of attraction is very complex. This confirms that the weakly dissipative system possesses hypersensitivity to the initial conditions. Thermal mechanical fluctuations σ β = |β| 2 = (g 0 /Ω m ) √ n th would be on the order of 10 −3 for realistic parameters with 100 thermal phonons and g 0 /Ω m ∼ 10 −4 . Panel (c): Zoom of the area within the white square shown in Panel (b). The zoomed picture displays the same degree of complexity as in Panel (b) and illustrate regions where the system is extremely sensitive even to minor changes in the initial conditions. The OM system operates in the weakly dissipative regime (κ = 0.1Ω m and γ = 10 −4 Ω m ).
In a real experiment, the mechanical oscillator's initial state is given by a thermal distribution at a given temperature T . In the classical regime studied here, one can use the Boltzmann (normal) distribution with zero mean and variance σ 2 m = |b| 2 = k B T / Ω m . Note that, though the equations of motion (3,4) contain only the parameter P , we will need also the OM coupling g 0 to calculate the standard deviation of the dimensionless variable β: For a typical value g 0 = 10 −4 Ω m and a thermal phonon number of 100, this amounts to σ β ∼ 10 −3 . As one can see in the simulations, this standard deviation covers a range of different attractors. The hypersensitivity to the initial conditions hampers a comprehensive study of the attractors of the weakly dissipative OM system. In addition to analyzing dynamics for different values of P and ∆, one would need also to consider many different initial conditions. This can be computationally very expensive, especially in the presence of chaotic attractors. The most common way to detect dynamical chaos is to calculate the Lyapunov exponent (LE) of a given trajectory. However, the convergence of the numerical methods available for calculating the LEs is usually slow. This calls for the development of alternative approaches. In the next section, we discuss such an alternative which is faster, reliably detects the chaotic attractors, tolerance set to 10 −9 and absolute tolerance set to 10 −13 . and moreover allows one to determine the dimensionality of the regular attractors.
3 The GALI method 3

.1 Indicators of dynamical chaos
An important task of any study of nonlinear dynamics is to distinguish regular and chaotic parts of the phase space in the most efficient way. A standard procedure for detecting chaotic trajectories is based on calculations of the maximal Lyapunov exponent (mLE). Let us consider the following general dynamical equations: One can start from a given trajectory x(t) and focus on small deviations w(t) from that trajectory. The linearized dynamics of w(t) is described by where Clearly, the mLE reflects the sensitivity of the trajectory x(t) to perturbations. A chaotic trajectory has positive mLE while regular trajectories have nonpositive mLE, making λ 1 a good indicator of chaotic dynamics. A numerical approximation for λ 1 can be obtained by calculating Λ(t) in Eq. (10) for a sufficiently large t, at which Λ(t) converges. This approach, however, has the drawback that the convergence of Λ(t) can be rather slow, and a long computation time is needed to learn whether λ 1 is positive or not. Many chaos indicators have been suggested to work around this problem; see Ref. [58]. We have used two of them: the SALI (Smaller ALignment Index) [50] and the GALI (Generalized ALignment Index) [52], which are especially well-suited for our goals. Before we discuss the SALI and the GALI, we have to define all LEs. Firstly, let us replace the n-dimensional vector w(t) in Eq. (9) by a n × n time-dependent matrix W (t), whose initial condition is W (0) = 1. The i-th column of W (t) describes the propagation of a perturbation acting in the i-th direction of the phase space at t = 0 (i.e. a perturbation proportional to the vector with components v j = δ j,i , where δ i,j is the Kronecker delta). Using the singular value decomposition, one can show that there is a set of n nonnegative real numbers {σ 1 , . . . , σ n }, and two sets of n orthonormal vectors, { v 1 , . . . , v n } and { u 1 , . . . , u n }, which satisfy the following equation [58]: This means that a perturbation in the direction of v i at t = 0 is mapped to a perturbation in the direction of u i multiplied by σ i at time t. The definition of the LEs reads where {σ j } are sorted in decreasing order. Eq. (12) gives all LEs of the dynamical system, and not only λ 1 . Figure 4: Panel (a): Evolution of two deviation vectors along a chaotic trajectory. Even if w 1 (0) ⊥ w 2 (0), the chaotic dynamics ensures that w 1 (t) w 2 (t) at t → ∞ provided that λ 1 > λ 2 . Panel (b): Evolution of the normalized deviation vectorŵ(t) in the case λ 1 = λ 2 . When t → ∞,ŵ(t) approaches the plane defined by u 1 and u 2 .
Let us return to the n-dimensional vector w(t), which satisfies Eq. (9). Using Eq. (11), we can rewrite w(t) for t → ∞ in the following way: where a, b denotes the inner product between a and b. Since t is very large, the term proportional to e λ 1 t dominates the time dependence of w(t) (provided that λ 2 < λ 1 ), such that Eqs. (10) and (12) are consistent. Now, we are in a position to introduce the SALI and the GALI. These indicators of chaos have been initially suggested for Hamiltonian systems, whose evolution preserves areas in the phase space. This means that the LEs are either zero, or appear in pairs with the same absolute value and opposite signs. The SALI and the GALI are constructed in a similar way, but the SALI is simpler; therefore, we start with the SALI: Consider two orthogonal initial conditions for Eq. (9), w 1 (0) ⊥ w 2 (0). Their evolution yields vectors w 1,2 (t) which become parallel to u 1 , and consequently to each other, at t → ∞; see Fig. 4(a). This holds true if λ 1 > λ 2 regardless of the initial condition. The SALI method uses this property to distinguish the chaotic dynamics from the regular one. Let us define whereŵ 1,2 (t) = w 1,2 /| w 1,2 | are unit vectors, and ( w 1 (0), w 2 (0)) = 0. The above discussion suggests that, if the dynamics is chaotic, the SALI tends to zero as t tends to infinity. In fact, the SALI decays exponentially to zero at the rate λ 1 − λ 2 [51]. If the dynamics is regular, all LEs are zero, and there is no reason for the alignment of vectors w 1 (t) and w 2 (t). The SALI does not decay to zero in this case. Thus, the SALI is a good chaos indicator for Hamiltonian systems provided that λ 1 = λ 2 . In the opposite case, where λ 1 = λ 2 , Eq. (13) suggests that w(t) tends to c 1 u 1 + c 2 u 2 , with c 1,2 depending on w(0). Therefore, w 1 (t) and w 2 (t) do not become parallel at t → ∞ but rather approach the plane defined by u 1 and u 2 ; see Fig. 4(b). The SALI does not decay to zero and a more advanced chaos indicator is needed. To construct it, we calculate the time evolution of a third deviation vector, w 3 (t), satisfying ( w 3 (0), w 1,2 (0)) = 0. We then compute the volume of the parallelepiped defined by the vectorsŵ 1,2,3 (t). It is given by the so-called GALI 3 : Hereŵ i = w i /| w i | is again the unit vector, and a ∧ b is the exterior product between the vectors a and b. One can show that GALI 3 ∝ exp(−2λ 1 t + λ 2 t + λ 3 t) [52], and it decays to zero exponentially quickly unless λ 1 = λ 2 = λ 3 . It can be shown that the GALI 3 decays to zero also on some regular orbits. However, such a non-chaotic decay is much slower as it follows a power law. This allows one to distinguish the chaotic and regular motion [52]. If the first (k − 1) LEs are equal to each other and positive, the chaotic and regular motion are distinguished by the GALI k [52]: (16) It is clear that GALI k (t) does not decay exponentially if and only if λ 1 = λ 2 = . . . = λ k . This applies to regular orbits where λ 1...k = 0. If the trajectory is chaotic, there exists a k which is smaller than the phase space dimension such that GALI k decays exponentially. One can show that SALI ∝ GALI 2 [52]. Therefore, we will refer only to the GALI in what follows.

The GALI method for dissipative systems
We have already mentioned that the GALI has been developed as an indicator of chaos for Hamiltonian systems, and its archetypal treatment generally does not work in the presence of dissipation and attractors.
Before extending the GALI to dissipative dynamics, let us first comment on the relation between attractors and LEs. The "attraction" of nearby orbits by the attractor comes from the fact that some LEs are negative (when the system is near the attractor). If the attractor is regular, all LEs are non-positive, and the number of zero-valued LEs is equal to the dimension of the attractor, see Chapter 10 of Ref. [59]. If all LEs are negative, the attractor is a fixed point. An attractor with only one zero-valued LE is a 1D curve in phase space, that is commonly called a limit cycle. An attractor which has p zero-valued LEs is a p-dimensional torus in phase space, that is dubbed a limit torus. The most complex attractors have positive and negative LEs, such that "attraction" co-exists with chaotic divergence of the trajectories. Those are called chaotic or strange attractors.
Consider now the GALI 2 in a dissipative system. On a limit cycle, w(t) approaches v 1 (regardless of the initial condition), and the GALI 2 decays exponentially to zero at a rate −λ 2 . On the other hand w(t) approaches v 1 on a chaotic attractor with λ 1 > λ 2 (again regardless of the initial conditions) and GALI 2 decays exponentially to zero at a rate λ 1 − λ 2 . We can conclude that the GALI k decays to zero both on the limit cycle and on the chaotic attractor for all possible values of k. Therefore, the GALI method cannot distinguish between the limit cycle and the chaotic attractor. We argue that the GALI is nevertheless useful for the study of dissipative systems because it is able to distinguish dynamics in the vicinity of the attractor from transient dynamics. Let us use Eqs. (13,16) to analyse the behaviour of the deviation vector modulus, | w(t)|, and of the GALI k on the different kinds of attractors: • Fixed point: all LEs are negative. Consequently, | w(t)| decays to zero exponentially quickly. The GALI k do not necessarily decay to zero since some of the LEs may have the same value.
• Limit cycle: λ 1 = 0, while all other LEs are negative. Consequently, | w(t)| does not decay to zero. The GALI k , on the other hand, decay to zero exponentially quickly for k ≥ 2.
• p-dimensional limit torus: λ 1 = . . . = λ p = 0, while all other LEs are negative. Consequently, | w(t)| and the GALI k for k ≤ p do not decay to zero. The GALI k for k > p, on the other hand, decays to zero exponentially quickly.
• Chaotic attractor : there are generically N 1 positive LEs and N 2 negative ones, where N 1,2 > 0. Consequently, | w(t)| grows and the GALI k>N 1 decay exponentially quickly. The behaviour of the GALI 2≤k≤N 1 (whether or not they decay to zero) depends on the degeneracy of the positive LEs.
Hence, when the trajectory is in the vicinity of an attractor, either | w(t)| or some GALI k must decay exponentially. Note that the inverse statement does not hold true: the fast decay of either | w(t)| or some GALI k cannot prove that the trajectory is in the vicinity of the attractor. The transient dynamics is more difficult for the analysis since one cannot make any general statement about the behaviour of | w(t)| or the GALI k when the trajectory is not close to any attractor. In principle, there is a possibility that | w(t)| or the GALI k could decay to very small values during the transient dynamics. On the other hand there is no generic reason for such a behaviour and it seems unlikely that many different deviation vectors would behave in such a way. Therefore, we will assume that whenever either | w(t)| or the GALI k decays to zero, the trajectory is in the vicinity of an attractor.
Once we know that the trajectory is in the vicinity of the attractor, knowing the properties of | w(t)| suffices to distinguish the chaotic attractors from the regular ones. If | w(t)| grows exponentially the attractor is chaotic; if there is no exponential growth of | w(t)| the attractor is regular. In the latter case, the GALI k provides the information about the dimensionality of the attractor. The ability of the GALI to detect the transient dynamics is especially important for a blue detuned OM system, since deterministic (non-chaotic) amplification of the mechanical motion represents the default behaviour in this regime and the growth of | w(t)| could be easily misinterpreted as a signature of chaos.
Armed with this novel understanding, we have successfully applied the GALI method to the dynamics of weakly dissipative OM systems. This will be the focus of the next section.
4 Applying the GALI method to OM systems 4

.1 Details of the implementation
In the previous Section, we have explained that the GALI method is a powerful tool for the analysis of the attractors of weakly dissipative OM systems because it allows one to detect the chaotic attractors very efficiently and to distinguish the regular attractors of different dimensionality. To study the nonlinear OM dynamics, we have solved the equations of motion (3,4) and analyzed the evolution of three deviation vectors w 1,2,3 (t) whose initial conditions are orthogonal. After this, we have calculated the GALI 2,3 (t). Three different pairs chosen from the three deviation vectors can generate three GALI 2 . We have calculated the GALI (w 1 ,w 2 ) 2 (t) based on w 1,2 (t) and the GALI (w 1 ,w 3 ) 2 (t) based on w 1,3 (t). We have used the average norm, the average GALI 2 , and the GALI 3 (t) for classification of the attractors. Specifically, we have assumed that any of these quantities has effectively "decayed to zero" when it becomes smaller than a given cutoff . We have chosen = 10 −6 . Our operational rules are: • If w(t) < , the attractor is a fixed point.
Note that, since the OM phase space is 4-dimensional, we could, in principle, come across limit tori with higher dimensionality. Their detection would require using the fourth vector w 4 (t) and constructing GALI 4 (t) because neither w(t) nor GALI 2 (t) nor GALI 3 (t) would drop below . We will show, however, that this is not the case for our choice of the parameters and of the initial conditions and, thus, the selected indicators suffice for our purposes. Fig. 5(a) shows a diagram as a function of the detuning ∆ and the drive power P which confirms the existence of various attractors in the phase space of a weakly dissipative OM system. We have already discussed that OM systems possess multistability: several attractors of different dimension can co-exists at given values of ∆ and P . Therefore, each pixel of the diagram has been obtained by solving the equations of motion for ten different initial conditions. Its color corresponds to the most "complex" attractor observed in these ten simulations. The attractors, sorted by increasing "complexity", are: fixed points, limit cycles, limit tori, transiently chaotic attractors, and chaotic attractors.

Attractors of weakly dissipative OM systems
To prove that our implementation of the GALI method yields reliable results, we show in Fig. 5(c) a similar diagram which has been obtained by calculating the mLE. One can observe qualitative similarity of the results generated by the two different methods, which confirms the validity of the diagram 5(a). On the other hand, this comparison also shows that the mLE method yields less detailed information and is unable to distinguish between the limit cycles and the limit tori.

Limit tori, Neimark-Sacker bifurcation, and transition to chaos
The diagram 5(a) displays the presence of four OM attractors with different dimensions: fixed points, limit cycles, limit tori and chaotic attractors. While there is a number of works addressing OM limit cycles (see e.g. Refs. [22][23][24][25]29,30,48]), and some works devoted to chaos in OM systems (see e.g. Refs. [27, 28,31-36, 60, 61]), studies of the OM limit tori, or quasiperiodic orbits, are scarce. We are aware of only one paper, Ref. [62], which reports the theoretical prediction of quasiperiodic OM orbits for parameters close to our choice. Quasiperiodic orbits were not observed in the strongly dissipative OM system, cf. Ref. [27].
We have detected the limit tori mostly in the range 0.6Ω m ≤ ∆ ≤ Ω m , which corresponds to the blue-detuned regime. The limit tori can be also found in the red detuned region, but they are rather rare there. Fig. 6(a) shows how a quasiperiodic orbit appears and disappears when the detuning is changed adiabatically 7 . We have plotted the spectrum of the position of the mechanical oscillator when the OM system is close to some attractor. There is only one peak in the spectrum at ∆ ∼ 0.7Ω m which means that the attractor is a limit cycle. A qualitative change occurs at ∆ 0.77Ω m and several peaks appear at larger ∆.  Figure 5: Attractor diagrams for the weakly dissipative OM system which have been generated by using the GALI method (Panels (a) and (b)) and the mLE method (c). Each pixel has been obtained after solving the equations of motion for ten different initial conditions. Its color corresponds to the the most "complex" attractor which we have detected for given ∆ and P . Note that the mLE method does not distinguish between the limit cycles and the limit tori. Panel (b) shows a zoom of the area within the white box in Panel (a). The OM system operates in the weakly dissipative regime (κ = 0.1Ω m and γ = 10 −4 Ω m ).
is quasiperiodic in this range and the attractor is now a limit torus 8 . All secondary peaks disappear at ∆ 0.89Ω m , and again only one peak is visible 9 ; the attractor becomes a limit cycle at ∆ > 0.89Ω m . These two transitions between a limit cycle and a limit torus agree with the diagram 5(a) and are known in the literature as the Neimark-Sacker bifurcation [63,64]. Fig. 6(b) shows time evolution of the mechanical degree of freedom. This trajectory is in the vicinity of a limit torus. The beating created by the sidebands is clearly visible. A similar phenomenon has been observed in Refs. [43,65].
Our remarkable finding is that the critical value of P , at which chaos appears, becomes considerably smaller when the dissipation is weak; compare the value P c ≈ 0.1 from Fig. 5(a) with P c ≈ 1.4 reported in Ref. [27] for the strongly dissipative case. We have discovered another qualitative difference between the strongly and weakly dissipative chaotic OM dynamics: chaos is observed mostly in the red detuned regime (∆ < 0) in the former case, while in the latter case it is observed mostly in the blue detuned regime (∆ > 0). In order to support the claim that the differences between the strongly and the weakly dissipative regimes depend on the sideband parameter κ/Ω m only, we have obtained the attractor diagram also for κ = 10 −1 Ω m and γ = 10 −3 Ω m . This diagram, which is not shown here, displays the same qualitative features The mechanical spectrum for an OM system close to an attractor. We have initially simulated the dynamics for P = 0.075Ω m and ∆ = 0.74Ω m until the system reached an attractor. Afterwards, the detuning ∆ has been slowly increased while P was kept fixed. The spectrum of the position of the mechanical oscillator has been computed during this process. Colors denote the absolute value of the spectrum, |S(ω)|. Only one peak is observed at ∆ < 0.77Ω m , i.e. the attractor is a limit cycle. Several other peaks appear at ∆ 0.77Ω m , i.e. the attractor becomes a limit torus. These secondary peaks disappear at ∆ 0.89Ω m ; for ∆ > 0.89Ω m , the attractor is again a limit cycle. Panel observed in Fig. 5(a), differing markedly from the results reported in [27]. This shows that the sideband parameter κ/Ω m is the only relevant parameter in our classification of strongly and weakly dissipative regimes. Though we have detected chaos for both positive and negative detuning, the chaotic region in the blue detuned part of the diagram looks more "dense" because pixels denoting chaotic dynamics agglomerate and are not isolated. We expect that small changes of the parameters inside the agglomerates cannot destroy chaotic dynamics. The chaotic region in the red detuned regime is very sparse and even subtle changes of the parameters are likely to convert the chaotic attractor to a regular one. Such a "sparse chaotic region" is shown in Fig. 5(b) which displays a zoomed part of Fig. 5(a) (the area within the white box in the red detuned region). One can see that the sparse chaotic region consists of very thin chaotic layers.
Close proximity of chaotic and quasiperiodic regions in the diagram Fig. 5(a) at ∆ > 0 provides a hint that OM systems can reach dynamical chaos via a route involving quasiperiodic orbits. To test this guess, we have investigated how an OM attractor behaves when the detuning is changed adiabatically such that the system starts in a quasiperiodic region of the parameters space and ends in a chaotic region. The results are shown in Fig. 7. The time evolution of | w|, the GALI 2 , and the GALI 3 are given in the upper Panels, while the lower Panels present the spectrum of the motion of the mechanical oscillator. At ∆ = 0.74Ω m , in Panels (a,e), | w| and GALI 2 oscillate around some nonzero values, while GALI 3 decays to zero exponentially quickly. Simultaneously, the mechanical spectrum has only two independent frequencies. Therefore, the attractor is a 2-dimensional torus. When ∆ is decreased (down to ∆ = 0.735Ω m , Panels (b,f), and further to ∆ = 0.732Ω m , Panels (c,g)), the behaviour of all three indicators remains qualitatively the same though more and more additional peaks (marking more frequencies) become visible and pronounced in the spectrum. The dynamical picture becomes qualitatively different at the smallest chosen detuning (∆ = 0.727Ω m , Panels (d,h)): | w| increases while the GALI 2,3 decay to zero exponentially. It means that the attractor is chaotic. This conclusion is confirmed by the dense nature of the mechanical spectrum. The transition to chaos depicted in Fig. 7 is called the quasiperiodic route to chaos [66]. It is characterized by the appearance of new frequencies when the control parameter (∆ in our study) is changed. The new frequencies must be commensurate with the basic two frequencies, see Fig. 7(e,f,g). If the new frequencies were incommensurate we would come across a higher dimensional torus and the GALI 3 would not vanish. One can notice a similarity between the quasiperiodic route to chaos and the period doubling cascade. Indeed, the dense chaotic spectrum is reached via an increasing number of new frequencies which are commensurate.

Classical transient chaos
Let us finally discuss another nonlinear phenomenon, which is captured by Fig. 5(a) but has not been revealed in previous studies of classical OM chaos. This is the well-known transient chaos. Its name perfectly reflects its main features: a dynamical system can display chaotic motion for a finite time interval after which its dynamics becomes regular. Dissipative transient chaos can be explained by a coexistence of different attractors, e.g. one attractor is chaotic and the other regular. Each attractor has its own basin of attraction. The basins could be separated by only an unstable periodic orbit. One can tune a parameter of the nonlinear system, η, such that the chaotic attractor approaches the unstable periodic orbit. At a critical value η = η c , the chaotic attractor "touches" the unstable periodic orbit. This phenomenon is called the boundary crisis [67,68] and it is one possible mechanism underlying transient chaos. One can imagine that, at η η c , a tiny fraction of the chaotic attractor penetrates the basin of attraction of the regular attractor. If a chaotic trajectory reaches this intersection region, where the chaotic attractor is entangled with the regular basin of attraction, it can be intercepted and "dragged" into the regular attractor. In other words, the chaotic attractor becomes leaky 10 .
Transient chaos also provides a possible route to chaos: changing η in the opposite direction results in the creation of a chaotic attractor at η η c . The time that a trajectory spends on the chaotic attractor before the leakage is typically very sensitive to the initial conditions. Nevertheless, we can define the average escape time τ esc . To this end, we select N 0 points in the leaky attractor, and use them as initial conditions of the nonlinear system. We then compute N (t), the number of trajectories remaining on the chaotic attractor at time t. The escape time can be found from the approximation N (t) N 0 exp(−t/τ esc ). Clearly, numerical approaches cannot distinguish the genuine chaotic trajectories and the transient trajectories with very large t esc 11 . We have used an empirical criterion: (i) trajectories which display chaotic motion during a time interval larger than T cutof f = 10 5 Ω −1 m are labeled "chaotic"; (ii) trajectories whose dynamics remains chaotic only for shorter times and becomes regular afterwards are labeled "transiently chaotic". Interested readers can find more details on transient chaos in the book [69].
Transient chaos in OM systems has been discussed for the first time in Ref. [60]. The authors of this paper argue that transient chaos underlies the breakdown of the quantumclassical correspondence in strongly dissipative OM systems, which display chaotic evolution in the classical regime and regular dynamics in the quantum one. To the best of our knowledge, the purely classical OM chaos has not yet been studied. We explore it in the weakly dissipative case. A representative example of transient chaos in classical OM is shown in Fig. 8. The time evolution of the optical and mechanical variables, Figs. 8(a,b), clearly manifests a crossover from the initially stochastic dynamics to subsequent regular motion. The crossover is obvious also in the behaviour of the modulus of the deviation vector | w|, and the GALI 2,3 , Fig. 8(c). Before the crossover, | w| increases while the GALI 2,3 decays exponentially, confirming that the trajectory is chaotic. | w| stops increasing at some time instant and oscillates around a nonzero value at longer times. This means that the trajectory becomes regular. The decay of the GALI 2,3 is cut at even much shorter times because of the finite numerical precision of the method which has been used to solve the equations of motion.
The chaotic fractions of the phase space are elaborately intertwined with the regions of transient chaos; see Fig. 9. We expect that this is a generic property, though details of the phase space (whether a given pixel belong to the genuine or transient chaos) are certainly sensitive to the cutoff time used in the empirical criterion explained above. Fig. 9(c) shows the same basin of attraction as that drawn in Fig. 9(a), but now with the doubled cutoff T cutof f = 2 × 10 5 Ω m . We note that many pixels, which were classified as chaotic in Fig. 9(a), are now classified as transiently chaotic in Fig. 9(c). Thus, many chaotic trajectories are actually transiently chaotic, but with a large escape time t esc .
The high complexity of the phase space results in hypersensitivity of the dynamics to the initial conditions. We have discussed this phenomenon already in Sect. 2.3.1; see Fig. 3. Fig. 9 suggests that the hypersensitivity is generic in weakly dissipative OM systems which possess multistability (co-existence of different attractors).

Experimental relevance of our results
There are several experimental works devoted to OM systems that report dissipation constants similar to (or even smaller than) those we have used for our numerical simulations. Weakly dissipative OM resonators can be fabricated in microwave systems [70][71][72], microresonators [73,74] and photonic crystals [75,76], to name just a few platforms. The detuning can usually be changed in a broad range. More important for investigations of the nonlinear effects is the accessible range of the driving strength, which governs the values of the parameter P . P itself is not convenient to describe the experiments, and it is better to consider the maximum cooperativityC. The value P ∼ 0.1 corresponds to a maximum cooperativityC ∼ 10 6 . This value agrees, for example, with the experimental value reported in Ref. [77]. We thus believe that the nonlinear phenomena described in the current paper can be explored in the near future in modern experiments. In particular, phenomena similar to the Neimark-Sacker bifurcation have been already observed experimentally in Refs. [43,65].
We note also that some platforms have dissipation constants substantially smaller than the values chosen for our study [71,72]. We have not considered such a weak dissipation but we think that nontrivial nonlinear phenomena could be found in less dissipative OM samples for substantially smaller values of P .

Conclusions
We have demonstrated that the classical nonlinear dynamics of an optomechanical resonator shows a great variety of nontrivial properties when the dissipation is weak. This regime had not received proper attention in the few previous studies dedicated to nonlinear OM dynamics, though it is of great experimental significance.
The phase space of the simplest OM system is four dimensional and includes two mechanical and two optical variables. High dimensionality and the presence of dissipation bring an extreme level of complexity to any systematic study. This is because analytical methods are basically unavailable while standard numerical approaches converge rather slowly. To overcome these technical difficulties, we have suggested a novel application of the GALI method, which was initially developed for Hamiltonian systems, to study attractors of the dissipative nonlinear OM system. Our approach has several advantages. Firstly, it has proved to be substantially faster than that based on an analysis of the maximal Lyapunov exponent. Even more importantly for our goals, it allows one to easily distinguish attractors of different dimensionality.
We have shown that weak dissipation strongly facilitates various nonlinear OM effects, which can appear at substantially lower laser power as compared to the previously studied strongly dissipative OM dynamics. In particular, weakly dissipative dynamics becomes chaotic at P ≈ 0.1 (see the definition in Sect. 2), one order of magnitude smaller than the typical values of P needed for chaos in the strongly dissipative case.
Our choice of parameters has allowed us to reveal multistability, i.e. the co-existence of different attractors. Their basins of attraction are very complex and entangled. As a result, a tiny variation of the initial conditions can completely change the dynamics on long time scales, since the trajectory is driven to a different attractor. Such a hypersensitivity to the initial conditions occurs even when the dynamics is regular and there are no chaotic attractors. We believe this to be a generic property of weakly dissipative OM systems.
Another generic feature reported in the current paper is the existence of quasiperiodic attractors, or 2-dimensional tori, in the OM phase space. We have investigated the transition from limit cycles to quasiperiodic orbits, which, in turn, can undergo a transition to chaos. The latter transition has some similarities to the well known period doubling cascade and provides a new route to chaos for OM systems. Finally, we have detected transient chaos. To the best of our knowledge, transient chaos has not been observed in previous studies of classical OM dynamics.
In spite of the great power of our numerical approach, we have not been able to obtain completely exhaustive information about weakly dissipative OM dynamics. This is because scanning all possible combinations of the four dimensionless parameters (rescaled power, detuning, mechanical and optical dissipation) and a broader range of the initial conditions is simply not feasible. We have focussed on exploring the phase diagram in terms of power and detuning, while keeping the dissipation values fixed. Thus, any complementary analytical method could be of great importance. We believe that an extension of the method suggested in Refs. [78][79][80] might help to achieve further progress. This method is based on the analysis of hyperbolic trajectories in phase space. It has initially been developed for ac driven dissipation-less dynamics. However, its generalization to the weakly dissipative case seems to be possible and promising.
We have argued that all the nonlinear OM phenomena which we have described are within the reach of state-of-the-art experiments in optomechanics. Moreover, it would be interesting to extend the present analysis to OM arrays, which are known to have a tendency towards complex and chaotic motion [41]. This could lead to exploring the complex interplay of the Anderson localization physics, first predicted in Ref. [81], and nonlinear OM dynamics.