1. Introduction

Cavity optomechanics have been extensively researched in recent years, towards the explorations of the radiation-pressure interaction between mechanical mode and optical field [1]. It provides new possibilities for observing many interesting quantum [2–11] and classical [12–21] nonlinearity phenomena. Due to the optomechanical nonlinearity, bistable, self-induced oscillations, and other classical nonlinear phenomena can be studied in this system. Since the nonlinear is very weak, the occurrence of these phenomena requires a relatively large driving strength. Most interesting, when the system is periodic, by increasing the driving, period-doubling bifurcation occurs. Continuously increasing the driving strength, the system states enter the chaotic region, which has been demonstrated by theoretical [22–31] experimental [32–35] studies. Chaos provides a powerful and reliable tool to send secret messages and realize optical sensing [36–40], because its future behaviors are sensitivity to initial conditions, and is complicated and unpredictable [41]. However, chaos has very high nonlinearity [33]. Due to the weak nonlinear optomechanical coupling, a large drive power is often required to exhibit chaotic motion. Therefore, the generation of chaos under low-power has recently attracted enormous attention in optomechanics [26,42,43].

To enhance the nonlinear of optomechanical system in low-power driving, many schemes have been proposed by different methods. Utilizing the optical parametric amplifier can enhance the intrinsic single-photon optomechanical coupling, then leading to the occurrence of the quantum phenomena such as photon blacked and phonon sideband [44]. In the study of the classical realm, introducing multiple mechanical modes can generate multiple optomechanically induced absorption [45]. Such multimode optomechanical system can be achieved by many promising platforms, such as circuit QED, toroidal optical whispering gallery mode cavity, and optomechanical phononic crystals [46–54]. Therefore, many works based on the multimode optomechancial system are reported in recent years, with finding many interesting physical phenomena such as dynamical phase transitions, mechanical entanglement, coherent state transfer, synchronization and double optomechanically induced transparency [55–79]. Furthermore, the cavity optomechanical system with multiple mechanical modes may have potential applications in amplifying optomechanical frequency conversion [80], storing optical signals [81], structuring quantum networks [82]. The physical phenomena in multimode optomechanical system are surprisingly rich and have not yet been fully explored. Therefore, we take our eye on such system with studying the nonlinear dynamics.

To date, there are few studies that have investigated the influence of multiple mechanical modes on the driving power thresholds of these nonlinear dynamics, especially of the chaotic motion. In this work, by the theoretical calculation and numerical simulations, we find that introducing multiple mechanical oscillators to optomechanical system can reduce the threshold of driving power for generating various nonlinear behaviors, such as bistability, period-doubling bifurcation and chaos. To get this result, firstly we analyze the bistable regime of such multimode optomechanical system in detail. The analytical result and its numerical verification convey that the cooperation of $N$ identical mechanical modes can reduce the bistable driving thresholds to $1/N$. The study result on bistability inspire us to investigate whether the driving power threshold of generating chaos can be decreased by increasing the number of mechanical modes $N$. Our numerical results reveal that the reducing degree of threshold for observing the period-doubling bifurcation and chaos is the same as bistabiliy. In addition, basing on the study of bistable state, we also find that changing the part signs of the optomechanical couplings does not change the nonlinear threshold. Furthermore, we verify again our results by a set of experimentally feasible parameters. This work may have potential applications in the field of modern information physics and optics.

This work is organized as follows: In Sec. 2., we introduce the system Hamiltonian and give the motion equations of the system. In Sec. 3., we discuss the influence of the number of the mechanical modes on the driving power thresholds of generating different nonlinear behaviors. In Sec. 4., we finally summarize our results.

2. Multimode optomechanical model and its motion equations

We consider a multimode optomechanical system where $N$ mechanical modes are individually coupled to cavity mode, as shown in Fig. 1. In a rotating frame with frequency $\omega _{l}$ of the driving field, the Hamiltonian of the system reads

 

Fig. 1. Schematic of the optomechanical system consisting of $N$ mechanical oscillators coupled to a common field.

Download Full Size | PDF

The dynamics of the system can be described by the Langevin motion equations: ${\partial }\hat {O}/{\partial }t=i[H,\hat {O}]-\Gamma \hat {O}+\sqrt {2\Gamma }\hat {O}_{\rm {in}}$, where $\Gamma$ represents the cavity dissipation $\kappa /2$ or mechanical damping $\gamma _{k}/2$, respectively; $\hat {O}=\hat {a}$, $\hat {b}_{k}$ with the input quantum noise $\hat {O}_{\rm {in}}=\hat {a}_{\rm {in}}$, $\hat {b}_{k,\rm {in}}$ being $\langle \hat {O}_{\rm {in}}\rangle =0$. Focusing on the mean response of the system, one can safely ignore the quantum correlations of photon-phonon in the semiclassical approximation [25,83] when the driving field is strong in the concerned weak-coupling regime. The Langevin motion equations can be solved under the semiclassical approximations [24,26] with each operator being split into a classical mean value and a quantum fluctuation, i.e., $\hat {a}=\alpha +\delta \hat {a}$ and $\hat {b}_{k}=\beta _{k}+\delta \hat {b}_{k}$. Here, we mainly concern the classical mean parts $\alpha$ and $\beta _{k}$, and the semiclassical Heisenberg-Langevin equations can be described as

For simplify, Eqs. (2) are represented as an evolution of a column real vector $\vec {o}=\left (\alpha _r, \alpha _i,\beta _{k r},\beta _{k i}\right )^{\rm T}$ with $\alpha =\alpha _r+i\alpha _i$ and $\beta _k=\beta _{k r}+i\beta _{k i}$. To fully explore the classical nonlinear behavior, we introduce the evolution of the perturbation $\vec {\delta }_{o}=\left (\delta _{\alpha _{r}}, \delta _{\alpha _{i}},\delta _{\beta _{kr}}, \delta _{\beta _{ki}}\right )^{\rm T}$ of $\vec {o}$ with infinitesimally different initial condition of Eqs. (2). By linearizing Eqs. (2), the dynamic evolution of the nearby trajectories can be expressed by the divergence equation $\dot {\vec {\delta }}_{o}=U\vec {\delta }_{o}$ [33], which characterizes the divergence of the nearby trajectories, where the coefficient matrix $U$ is

3. Analytical analysis and numerical simulation

When the system reaches to steady state, Eqs. (2) can be expressed as

In Fig. 2(a), we present the intensity $I$ as functions of $P$ for $\Delta _{0}/\kappa =-1.2$ in the optomechanical systems with $N=1$ ($N=3$). The pink and red dashed lines correspond to the real part eigenvalues of $U$ being non-negative, corresponding to the unstable solutions. In the regime of $P_{-}<P<P_{+}$, we divide the three solutions of Eq. (5) into upper ($I>I_{+}$), middle ($I_{-}<I<I_{+}$), and lower ($I<I_{-}$) branches. Here the black dots of ($P_{+}, I_{-}$) and ($P_{-}, I_{+}$) are given by Eq. (6), one can see that the numerical values are consistent with the analytical results, and $P_{\pm }$ are inversely proportional to $M$. It clearly shows that the lower and middle branches are the same as the general case. Different from the general case, the upper branch becomes unstable with increasing $P$ in our work, as shown in Fig. 2(a). Therefore, one can see a critical value $P_c$ to distinguish the stability and instability of the upper branch. To give a quantitative relationship between $M$ and $P_{c}$, it needs to give an analytical expression of $P_{c}$. Since $P_{c}$ can only be obtained through $U$, it is difficult for larger $N$. Without loss of generality, we consider all the mechanical modes are the same, i.e., $\beta _1=\beta _2=\cdots =\beta _k=\cdots =\beta _N=\beta _0$ with $g_{k}=g_{0}$, $\omega _{k}=\omega _{0}$ and $\gamma _{k}=\gamma _{0}$. And then all the mechanical modes can be regarded as a collective mechanical mode $\beta =N\beta _0$. In this case, the system evolution can be described by the collective mechanical mode $\beta$ and one optical mode $\alpha$. The perturbation vector can be reduced as $\vec {\delta }_{o}=\left (\delta _{\alpha _{r}}, \delta _{\alpha _{i}},\delta _{\beta _{r}}, \delta _{\beta _{i}}\right )^{\rm T}$, and its divergence equation satisfies $\delta \dot {o}_{\rm tot} = {U}_{\rm tot}\delta o_{\rm tot}$ with a 4$\times$4 coefficient matrix

 

Fig. 2. Bistability in the optomechanical systems with one mechanical mode ($N=1$) and three mechanical modes ($N=3$). (a) The dependence of intensity $I$ on the driving strength $P$ for $\Delta _{0}/\kappa =-1.2$, and the solid (dashed) line represents the stable (unstable) state. On the curve of $N=1$, the black dots are $(P_{-},I_{+})=(1.14, 5.38\times 10^{11})$, $(P_{c},I_{c})=(1.20, 6\times 10^{11})$ and $(P_{+},I_{-})=(1.48, 2.61\times 10^{11})$. The inset of (a) shows $P_{-}$, $P_{+}-P_{-}$, and $P_{c}-P_{-}$ versus $N$. The phase diagrams for $N=1$ and $N=3$ versus the detuning $\Delta _0$ and the driving strength $P$ are plotted in (b) and its inset, respectively. Here, the blue ① and green ④ areas represent one solution, and the solution in region ① (④) being stable (unstable) corresponds to $P<P_-$ ($P>P_+$); both the yellow ② and dark yellow ③ areas indicate three solutions, and the region ② (③) having two (one) stable solutions corresponds to $P_{-}<P<P_{c}$ ($P_c<P<P_+$). The parameters are $\omega _{0} = 1$, $\kappa = 1$, $\gamma _{0} = 1\times 10^{-4}$, $g_0 = 1\times 10^{-6}$, $\omega _{k} = \omega _{0}$, and $\gamma _{k}=\gamma _{0}$.

Download Full Size | PDF

Note that in stable regime [blue solid lines in Fig. 2(a)] the systemic steady state is a fixed point, at which the state never change with time. In the unstable regime shown by the red dashed lines, many nonlinear dynamical behaviors such as period-doubling and chaotic motion, may occur. However, the unstable state in the middle branch (pink dashed lines) is a state that the system can’t be touched finally. From the inset of Fig. 2(a), one can obtain that the critical driving powers $P_{\pm }$ as well as $P_c$ between the stable and unstable regimes can be cut down by increasing $N$. Moreover, Fig. 2(b) shows that one also can obtain the same phase diagram in lower driving power with increasing the number of mechanical modes. The explanation can be got from the analytical expressions of $P_{\pm }$ and $P_c$, which are inversely related to $M$ ($M \propto N$). Therefore, the threshold value of driving power may be reduced by increasing the number of mechanical modes ($N$) to push the system into bistable, period-doubling and chaotic regime. In the following, we will discuss the impact of $N$ on the driving threshold in detail, and further show that increasing the number of mechanical modes can significantly decrease the driving threshold.

As we know, in the bistable regime what the steady state the system can arrive at tightly depends on the initial condition. To fully display the nonlinear dynamics, we plot the bifurcation diagrams as the functions of the driving power $P$ with the initial conditions satisfying upper branch and lower branch in Figs. 3(a) and 3(b), respectively. In the regime of $P_{-}<P<P_{+}$, for satisfying the upper (lower) branch the initial conditions are chosen near the upper (lower) solution of Eqs. (4). For $P<P_{-}$ and $P>P_{+}$, the initial conditions can be chosen arbitrarily. In addition, to conveniently distinguish the nonlinear type, the maximal Lyapunov exponents (MLE) calculated by the jacobian matrix algorithm [25] are correspondingly presented in Figs. 3(c) and 3(d). A positive MLE is an indication to refer to a chaotic system. In this work, the bifurcation diagrams are completed by the peak values $x_{k, {\rm ext}}$ of the mechanical coordinate $x_k=(\beta _k+\beta _k^*)/2$. At $P=1.14$, we see that $x_{1,\textrm{ext}}$ suddenly jumps in Fig. 3(a) and the MLE abruptly changes in Fig. 3(b). This is because that before $P=1.14$ the steady state is in the lower branch and when $P>1.14$ the steady state escape to the upper branch, as seen in Fig. 2(a). The dynamical evolution shows that the steady state of the system is a fixed value with the driving power being $P_1=1.08$ or $P_2=1.17$ [see the inset of Fig. 3(a)]. When the driving power is increased beyond $P_c=1.20$, the upper branch goes into unstability and $x_{\rm 1}$ become oscillation on time, as shown in the inset by $P_3=1.30$. Note that in Fig. 3(c), the MLE being zero at $P=P_c$ does not correspond to period-doubling bifurcation, but heralds a new state ($x_{\rm 1}$ being oscillation) the system will enter into from the fixed point ($x_{\rm 1}$ being static). However, MLE becoming zero at $P=1.418$ means the period-doubling bifurcation taking place, which can be clearly seen from Fig. 3(a). This is further demonstrated by the corresponding inset, which shows that the oscillator period is doubled comparing the cases of $P_3=1.30$ and $P_4=1.45$. When the initial condition satisfies the lower branch, the steady state does not change until $P=1.45$ as shown in Fig. 3(b), which shows that $x_{\rm 1}$ is static with $P=1.45$ and becomes oscillation with $P=1.49$. Meanwhile, at $P=1.49$ the sudden change of the system causes the MLE increasing abruptly, seen in Fig. 3(d).

 

Fig. 3. Bifurcation diagram of mechanical oscillator $x_{1,\rm {ext}}$ with the initial conditions satisfying upper branch (a) and lower branch (b) at $\Delta _{0}/\kappa =-1.2$, and the corresponding MLEs are presented in (c) and (d), respectively. The insets of (a) and (b) give the time evolution of mechanical mode for different $P$ marked by the red dots $(P_1, P_2, P_3, P_4, P_5)=(1.08, 1.17, 1.30, 1.45, 1.49)$. Here, $N=1$ and the other parameters are the same as Fig. 2.

Download Full Size | PDF

Further increasing the driving power, the system becomes chaos characterized by MLE above zero near $P=1.52$ as shown both in Figs. 3(a) and 3(b), which also shows whether the system being chaos is irrelevant with the initial condition satisfying the upper branch or lower branch. In order to intuitively show the relationship between the different nonlinear behaviors and the driving power $P$ as well as the number of mechanical modes $N$, we present the evolution of the mechanical oscillator in temporal domain, phase space, and frequency domain in Fig. 4 with $N=1$ and $N=3$. Here, we set the initial condition satisfying the upper branch. For $N=1$, when $P=1.44$ the steady state of the system is periodic as shown in Fig. 4(a), and its phase space only has two limited circles as well as several discrete peaks in frequency domain. Increasing the driving power to $P=1.5$, Fig. 4(d) shows that the oscillation period is doubled, and correspondingly the number of the limited circles in phase space and peaks in frequency domain are also doubled. When $P=1.53$, one can see that the mechanical mode $x_1$ becomes irregular on time, an almost solid swath in phase space, and continuous in frequency domain, i.e., the system is chaotic at $P=1.53$. Interestingly, increasing the number of mechanical modes, one can see the same nonlinear behaviors just in a lower driving. Such as, for $N=3$ the periodic, period-doubling and chaotic motions are observed at $P=0.38$, $0.5$ and $0.51$, respectively. Compared to $N=1$, these values of driving power are reduced by $1/3$.

 

Fig. 4. The evolution of $N x_{1}$ on time in (a), (d) and (g); and of the trajectories in phase space in (b), (e) and (h); and its power spectra in (c), (f) and (i). The driving power in (a)-(c), (d)-(f) and (g)-(i) are $P=1.44/N$, $1.5/N$ and $1.53/N$, respectively. All plots are given by $N=1$ and $N=3$, and the other parameters are the same as Fig. 3.

Download Full Size | PDF

Our discussion about the stability conditions shows that some critical driving threshold values such as $P_\pm$ and $P_c$ are inverse to the number of the mechanical modes $N$, so are the nonlinear behaviors like period-doubling and chaos concluded by Fig. 4. Therefore, we can obtain chaos in a lower driving power by increasing the number of the mechanical modes in optomechanical system. As an example, when $N=10$ we plot $x_{1,\textrm{ext}}$ as the functions of driving power $P$ in Fig. 5(a). At $P=0.142$, Fig. 5(a) shows that the period-doubling bifurcation takes place, which is further evidenced by Fig. 5(b) with MLE being zero. However, the bifurcation taking place needs a larger driving power $P=1.418$ when $N=1$ (see Fig. 3). In addition, Fig. 5(a) indicates that the system enters chaos near $P=0.15$ with the MLE beyond zero begin occurring [see Fig. 5(b)], yet the value is $P=1.49$ for $N=1$ (see Fig. 3). Considering a set of feasible experimental parameters [48–55,85] $\lambda =1550$ nm (the wavelength of driving laser), $\omega _k/2\pi = 51.8$ MHz, $m_k=20$ ng, $G_0/2\pi =-12$ GHz/nm, $g_k= G_0\sqrt {\hbar /{2 m_k \omega _k}}$, $\kappa /2\pi = 51.8$ MHz, $\gamma _k/2\pi = 41$ KHz, and $\Delta = -1.2 \kappa$, for $N=1$ chaotic motion occurs at $P= 18.11$ mW, which is cut down to $0.905$ mW with increasing the number of mechanical modes to $N=20$. This finding may provide theoretical guidance for getting chaos in lower power driving.

 

Fig. 5. (a) Bifurcation diagram of the mechanical oscillator $x_{1,\rm {ext}}$ versus $P$, and (b) the corresponding MLE as the function of $P$ with the initial value satisfying upper branch and $N=10$. The two red dots in each plot represent $P=0.142$ and $0.149$. The other parameters are the same as Fig. 3.

Download Full Size | PDF

Furthermore, one should note that partly changing the sign of $g_k$ does not change the stability conditions due to $M\propto g_k^2$, so $P_\pm$ and $P_c$ are independent on the sign of $g_k$. Naturally, a hypothesis that the sign of $g_k$ has no influence on the power driving threshold of generating chaos comes. To verify this hypothesis, firstly we plot the temporal evolution of the intracavity field intensity ($I$) and the mechanical modes ($x_1$ and $x_2$) in Fig. 6(a) and 6(b), respectively. For $g_{1}=g_{0}$ the evolution of $I$ is periodic, and for $g_{1}=-g_{0}$ the evolution is the same but only a phase different [see Fig. 6(a)]. In addition, when $g_{1}=g_{0}$ the oscillator of $x_1$ is synchronous with $x_2$, but only with a sign different when $g_{1}=-g_{0}$, i.e., $g_1 x_1=g_2 x_2$. To a certain extent, these demonstrate that partially changing the sign of $g_k$ does not affect the property of nonlinear. The underlying reason is that the choice of the positive direction of the mechanical motion has no influence on the physical effects. Coincidentally, the sign of $g_{k}$ just depends on that choice. In order to visually view the influence of partially changing the sign of $g_k$ on nonlinear behaviors, we present the bifurcation diagram of the mechanical oscillator $x_{1,\textrm{ext}}$ as the function of $P$ in Fig. 7(a) for $g_1=-g_0$. For convenient comparison, the MLE versus $P$ is plotted with $g_1=g_0$ in Fig. 7(b). From Fig. 7(a), one can see a sudden transition at $P=0.491$, and in Fig. 7(b) the MLE manifests a sharp increase implying that the $x_{1,\textrm{ext}}$ at $P=0.491$ is also saltatory. Moreover, at $P=0.504$, Fig. 7(a) ($g_1=-g_0$) shows that the system is chaotic, and for $g_1=g_0$ the MLE is above zero [see Fig. 7(b)].

 

Fig. 6. (a) The evolution of intracavity field intensity $I$ with $g_1=g_0$ (solid blue) and $g_1=-g_0$ (solid green), and (b) of the first mechanical mode $x_1$ with $g_1=g_0$ (dashed green) and $g_1=-g_0$ (solid green) and second mechanical mode $x_2$ with $g_1=g_0$ (solid blue) and $g_1=-g_0$ (light red) for $N=3$. The other parameters are the same as Fig. 3.

Download Full Size | PDF

 

Fig. 7. (a) Bifurcation diagram of the amplitude of the mechanical oscillator $x_{1,\rm {ext}}$ versus $P$ for $g_1=-g_0$, and (b) the MLE versus $P$ for $g_1=g_0$. In both (a) and (b), the initial conditions satisfying the lower branch are set, and the two red dots represent $P=0.491$ and $0.504$. Here $N=3$ and the other parameters are the same as Fig.  3.

Download Full Size | PDF

4. Conclusion

In conclusion, a set of nonlinear dynamics are investigated in a multimode optomechanical system with $N$ mechanical oscillators coupled to a common cavity field. We show that the values of critical driving power are scaled by $1/N$ when increasing the mechanical modes to $N$. Basing on this feature, we further demonstrate that the threshold value of driving power is reduced to $1/N$ in such multimode optomechanical system with $N$ mechanical modes for generating chaos. Moreover, we find that the signs of the coupling strength partially changed has no influence on the driving threshold value for generating corresponding nonlinear dynamics. Our work may enrich the nonlinear dynamics of cavity optomechanics and have potential application in modern information science and the technology of optical sensing.