Introduction

Synchronization, a fascinating natural phenomenon observed across various systems, manifests in nature through coordinated behaviors such as the flocking of birds or the synchronous flashing of fireflies1, and extends its influence into social systems2, chemical reactions3, biological processes4,5, and so on. The concept of spontaneous synchronization was first proposed in the mid-17th century by Christian Huygens, who observed an intriguing phenomenon where two pendulum clocks, hanging from the same wooden beam, finally swung in perfect unison despite the different natural frequencies and initial positions6. Since then, particularly since the 20th century, the study of synchronization has flourished and expanded significantly, encompassing a broad range of disciplines and applications.

In recent years, cavity optomechanics, which contains nonlinear couplings between mechanical oscillators and optical modes, are widely employed to explore the radiation pressure-mediated interactions between optical and mechanical freedoms7. Due to its unique nonlinear characteristics and the versatility offered by a wide range of adjustable parameters, cavity optomechanics provide an ideal platform for the study of a variety of interesting quantum phenomena such as quantum entanglement8,9, quantum sensing10,11, quantum squeezing12,13, and quantum information processing14,15. Furthermore, studies on quantum entanglement and coherence in various hybrid optomechanical systems16,17,18 have significantly enriched our understanding of quantum interactions. Additionally, research on tunable optical responses19,20 has been pivotal in advancing the capabilities of modern quantum technologies. Benefiting from the enhanced ability to control micro-scale systems21, interest in synchronization has gradually expanded into the quantum ___domain. Unlike classical systems where synchronization concepts are well-established, quantum systems present unique challenges due to the constraint imposed by Heisenberg’s uncertainty principle. To address this challenges, A. Mari and his colleague introduced the notions and metrics of quantum synchronization and quantum phase synchronization22, significantly enriching the field. Recently, quantum synchronization has been explored in various systems, including van der Pol oscillators23,24,25,26,27, atomic ensembles28,29, Bose-Einstein condensates30, and cavity optomechanical systems31,32. Alongside, theoretical exploration in quantum synchronization continues to advance, with the introduction of novel metrics such as quantum \(\phi\) synchronization33,34, providing a comprehensive framework for evaluating synchronization phenomena. Moreover, various methods have been proposed to enhance quantum synchronization35,36,37. These developments present great potential for numerous applications, such as enhancing security and stability in fiber-optic quantum clock synchronization38, facilitating quantum communication39,40, and improving the transmission of quantum information across networks41.

While the study of quantum synchronization is flourishing42,43,44,45, research on quantum synchronization within complex systems remains relatively sparse46,47,48, and investigations into quantum synchronization in many-body optomechanical systems are particularly scarce. To bridge this gap, we propose two methods aimed at enhancing quantum synchronization among three oscillators. Our numerical calculations indicate that, in the absence of periodic modulation, quantum synchronization is nearly negligible. However, appropriate periodic modulation can achieve quantum synchronization levels of up to 0.8 across all subsystems in both open and closed configurations.

Fig. 1
figure 1

Schematic illustration of an optomechanical system comprising three mechanical oscillators with frequencies \(\omega _{m1}\), \(\omega _{m2}\) and \(\omega _{m3}\). Each oscillator is placed in a separate optical cavity, where the cavity modes \(a_1\), \(a_2\) and \(a_3\) are driven by laser beams with common amplitude E. These driving fields induce displacements \(q_1\), \(q_2\) and \(q_3\) in the oscillators from their equilibrium positions caused by the radiation pressure of the driving fields. The subsystems are connected via optical fibers with coupling constant \(\lambda _{1,2,3}\), enabling the system to switch between open and closed configurations.

Full size image

These results may provide a foundational basis for further studies of quantum synchronization in intricate optomechanical networks.

Model and methods

To investigate quantum synchronization in many-body system, we propose a simple prototype consisting of three distinct Fabry-Perot cavities. As depicted in Fig. 1, each cavity contains a tiny mechanical oscillator with unique natural frequency \(\omega _{m1,m2,m3}\). These oscillators interact with their respective cavity modes via radiation pressure. The cavities are interconnected through optical fibers, with coupling strengths \(\lambda _{1}\), \(\lambda _{2}\), and \(\lambda _{3}\). Notably, the connectivity between the first and the last cavities, \(\lambda _{3}\), is adjustable, enabling a switch between an open chain configuration and a closed ring configuration. Each cavity is uniformly driven by laser beams at a common frequency \(\omega\) and strength E, while the driving fields or the internal cavity modes can be periodically modulated via the acousto-optical effect or piezo-electric effect49, respectively. Then the total Hamiltonian after rotating can be expressed as22,35:

$$\begin{aligned} \begin{aligned} H_{open}=&\sum _{j=1,2,3} \left\{ -\Delta _{j} \left[ 1+\eta _C\cos \left( \Omega _Ct \right) \right] a_j^{\dagger }a_{j}+\frac{\omega _{mj} }{2} \left( p_{j}^{2}+q_{j}^{2}\right) -ga_{j}^{\dagger }a_{j}q_{j}+iE\left[ 1+\eta _{D}\cos \left( \Omega _{D}t \right) \right] \left( a_{j}^{\dagger } -a_{j} \right) \right\} \\ &+\lambda _{1}\left( a _{1}^{\dagger }a_{2}+a _{2}^{\dagger }a_{1} \right) +\lambda _{2}\left( a _{2}^{\dagger }a_{3}+a _{3}^{\dagger }a_{2} \right) , \end{aligned} \end{aligned}$$
(1)

for the open configuration, and:

$$\begin{aligned} \begin{aligned} H_{closed}=&\sum _{j=1,2,3} \left\{ -\Delta _{j} \left[ 1+\eta _C\cos \left( \Omega _Ct \right) \right] a_j^{\dagger }a_j+\frac{\omega _{mj} }{2} \left( p_{j}^{2}+q_{j}^{2}\right) -ga_{j}^{\dagger }a_{j}q_{j}+ iE\left[ 1+\eta _{D}\cos \left( \Omega _{D} t \right) \right] \left( a_{j}^{\dagger } -a_{j} \right) \right\} \\ &+\lambda _{1}\left( a _{1}^{\dagger }a_{2}+a _{2}^{\dagger }a_{1} \right) +\lambda _{2}\left( a _{2}^{\dagger }a_{3}+a _{3}^{\dagger }a_{2} \right) +\lambda _{3}\left( a _{3}^{\dagger }a_{1}+a _{1}^{\dagger }a_{3} \right) , \end{aligned} \end{aligned}$$
(2)

for the closed configuration. For convenience, here we set \(\hbar\) = 1. \(q_{j}\) and \(p_{j}\) are the dimensionless position and momentum operators of the oscillator, \(a_{j}^{\dagger }\) and \(a_{j}\) are the creation and annihilation operators for the cavity modes at frequency \(\omega _{cj}\), the detuning of the \(j-\)th cavity can be defined as \(\Delta _{j} =\omega -\omega _{cj}\), and g is the optomechanical coupling coefficient, which arises due to the radiation pressure. Additionally, we assume that all driving fields (cavity modes) are uniformly modulated at the same frequency \(\Omega _{D}(\Omega _{C})\) and amplitude \(\eta _{D}(\eta _{C})\).

To investigate the dynamical evolution of the system, we take the dissipation into account and derive the Langevin equations from the specified Hamiltonian as follows50,51,52:

$$\begin{aligned} \begin{aligned} \dot{q}_{j}=&\;\omega _{mj}p_{j},\\ \dot{p}_{j}=&-\omega _{mj}q_{j}-\gamma _{m}p_{j}+ga_{j}^{\dagger }a_{j}+\xi _{j},\\ \dot{a}_{1} =&-\left\{ \kappa -i\Delta _{1}\left[ 1+\eta _C\cos \left( \Omega _Ct \right) \right] \right\} a_{1}+iga_{1}q_{1}+E\left[ 1+\eta _D\cos \left( \Omega _Dt \right) \right] \\&-i\lambda _{1}a_{2}-i\lambda _{3}a_{3}+\sqrt{2\kappa }a_{1}^{in},\\ \dot{a}_{2}=&-\left\{ \kappa -i\Delta _{2}\left[ 1+\eta _C\cos \left( \Omega _Ct \right) \right] \right\} a_{2}+iga_{2}q_{2}+E\left[ 1+\eta _D\cos \left( \Omega _Dt \right) \right] \\&-i\lambda _{1}a_{1}-i\lambda _{2}a_{3}+\sqrt{2\kappa }a_{2}^{in},\\ \dot{a}_{3} =&-\left\{ \kappa -i\Delta _{3}\left[ 1+\eta _C\cos \left( \Omega _Ct \right) \right] \right\} a_{3}+iga_{3}q_{3}+E\left[ 1+\eta _D\cos \left( \Omega _Dt \right) \right] \\&-i\lambda _{2}a_{2}-i\lambda _{3}a_{1}+\sqrt{2\kappa }a_{3}^{in}. \end{aligned} \end{aligned}$$
(3)

Here, \(\kappa\) represents the decay rate of the cavity, while \(\gamma _m\) denotes the mechanical damping rate of the m-th oscillator. The input noise operator \(a_{j}^{in}\) for the j-th cavity mode and the random noise operator \(\xi _{j}\) for the j-th mechanical oscillator both exhibit zero mean values. Under the Markovian approximation, they respectively satisfy the following correlation formulas53:

$$\begin{aligned} \begin{aligned} \begin{array}{c} \left\langle a_{j}^{in\dagger }\left( t \right) a_{j^{'} }^{in}\left( t^{'} \right) + a_{j^{'} }^{in}\left( t^{'} \right) a_{j}^{in\dagger }\left( t \right) \right\rangle =\delta _{jj^{'}}\delta \left( t-t^{'} \right) , \\ \left\langle \xi _{j}\left( t \right) \xi _{j^{'} }\left( t^{'} \right) + \xi _{j^{'} }\left( t^{'} \right) \xi _{j}\left( t \right) \right\rangle =\gamma _{m}\left( 2\bar{n}_{bath}+1 \right) \delta _{jj^{'}}\delta \left( t-t^{'} \right) , \end{array} \end{aligned} \end{aligned}$$
(4)

where \(\bar{n}_{bath}=\left[ \textrm{exp}\left( \hbar \omega _{mj}/k_{b}T \right) -1\right] ^{-1}\) is the mean phonon number of the mechanical bath at temperature T54,55,56.

To characterize the quantum dynamical features of the system, we employ the mean-field approximation to linearize the equations50,57. Consequently, we express the operators in Eq. (3) as the sum of their mean values O and the quantum fluctuation terms \(\delta o\) around these means. This approach enables us to derive the mean-value equations:

$$\begin{aligned} \begin{aligned} \dot{Q}_{j}=&\;\omega _{mj}P_{j},\\ \dot{P}_{j}=&-\omega _{mj}Q _{j}-\gamma _{m} P_{j}+g\left| A_{j}\right| ^{2},\\ \dot{A}_{1}=&-\left\{ \kappa -i\Delta _{1}\left[ 1+\eta _{C}\cos \left( \Omega _{C}t\right) \right] \right\} A_{1}+igA_{1}Q_{1}+E\left[ 1+\eta _{D}\cos \left( \Omega _{D}t\right) \right] \\ &-i\lambda _{1}A_{2}-i\lambda _{3}A_{3},\\ \dot{A}_{2}=&-\left\{ \kappa -i\Delta _{2}\left[ 1+\eta _{C}\cos \left( \Omega _{C}t\right) \right] \right\} A_{2}+igA_{2}Q_{2}+E\left[ 1+\eta _{D}\cos \left( \Omega _{D}t\right) \right] \\ &-i\lambda _{1}A_{1}-i\lambda _{2}A_{3},\\ \dot{A}_{3}=&-\left\{ \kappa -i\Delta _{3}\left[ 1+\eta _{C}\cos \left( \Omega _{C}t\right) \right] \right\} A_{3}+igA_{3}Q_{3}+E\left[ 1+\eta _{D}\cos \left( \Omega _{D}t\right) \right] \\ &-i\lambda _{2}A_{2}-i\lambda _{3}A_{1}, \end{aligned} \end{aligned}$$
(5)

and the linear equations for the quantum fluctuation terms:

$$\begin{aligned} \begin{aligned} \dot{\delta q_{j}} =&\;\omega _{mj}\delta p_{j},\\ \dot{\delta p_{j}} =&-\omega _{mj}\delta q_{j}-\gamma _{mj}\delta p_{j}+g\left( A_{j}\delta a_{j}^{\dagger } +A_{j}^{*}\delta a_{j}\right) +\xi _{j},\\ \dot{\delta a_{1}}=&-\left\{ \kappa -i\Delta _{1}\left[ 1+\eta _{C}\cos \left( \Omega _{C}t\right) \right] \right\} \delta a_{1}+ig\left( A_{1}\delta q_{1}+Q_{1}\delta a_{1}\right) \\ &-i\lambda _{1}\delta a_{2}-i\lambda _{3}\delta a_{3}+\sqrt{2\kappa }a_{1}^{in},\\ \dot{\delta a_{2}}=&-\left\{ \kappa -i\Delta _{2}\left[ 1+\eta _{C}\cos \left( \Omega _{C}t\right) \right] \right\} \delta a_{2}+ig\left( A_{2}\delta q_{2}+Q_{2}\delta a_{2}\right) \\ &-i\lambda _{1}\delta a_{1}-i\lambda _{2}\delta a_{3}+\sqrt{2\kappa }a_{2}^{in},\\ \dot{\delta a_{3}}=&-\left\{ \kappa -i\Delta _{3}\left[ 1+\eta _{C}\cos \left( \Omega _{C}t\right) \right] \right\} \delta a_{3}+ig\left( A_{3}\delta q_{3}+Q_{3}\delta a_{3}\right) \\ &-i\lambda _{2}\delta a_{2}-i\lambda _{3}\delta a_{1}+\sqrt{2\kappa }a_{3}^{in}. \end{aligned} \end{aligned}$$
(6)

Further, we introduce the quadrature components \(\delta x_{j}=\frac{1}{\sqrt{2}}\left( \delta a_{j}^{\dagger }+\delta a_{j}\right)\) and \(\delta y_{j}=\frac{i}{\sqrt{2}}\left( \delta a_{j}^{\dagger }-\delta a_{j}\right)\) for the fluctuation operators, and the corresponding input noise terms for the cavity modes are represented as \(x_{j}^{in} =\frac{1}{\sqrt{2}}\left( a_{j}^{in\dagger }+a_{j}^{in}\right)\) and \(y_{j}^{in} =\frac{i}{\sqrt{2}}\left( a_{j}^{in\dagger }-a_{j}^{in}\right)\), then we can rewrite Eq. (6) as:

$$\begin{aligned} \dot{u}=\;Mu+n, \end{aligned}$$
(7)

where u is a column vector combining all quadrature components:

$$\begin{aligned} u = (\delta q_{1}, \delta p_{1}, \delta x_{1}, \delta y_{1}, \delta q_{2}, \delta p_{2}, \delta x_{2}, \delta y_{2}, \delta q_{3}, \delta p_{3}, \delta x_{3}, \delta y_{3})^T, \end{aligned}$$
(8)

and n denotes the combined noise vector:

$$\begin{aligned} n = (0, \xi _1, \sqrt{2\kappa }x_1^{in}, \sqrt{2\kappa }y_1^{in}, 0, \xi _2, \sqrt{2\kappa }x_2^{in}, \sqrt{2\kappa }y_2^{in}, 0, \xi _3, \sqrt{2\kappa }x_3^{in}, \sqrt{2\kappa }y_3^{in})^T. \end{aligned}$$
(9)

M is a \(12\times 12\) time-dependent coefficient matrix, defined differently for closed and open configurations as:

$$\begin{aligned} M=\begin{pmatrix} M_{1}& M_{01}& M_{03}\\ M_{01}& M_{2}& M_{02}\\ M_{03}& M_{02}& M_{3} \end{pmatrix}, \end{aligned}$$
(10)

applicable when \(\lambda _{3} \ne 0\), corresponds to the closed configuration, and:

$$\begin{aligned} M=\begin{pmatrix} M_{1}& M_{01}& M_{\textrm{zero}}\\ M_{01}& M_{2}& M_{02}\\ M_{\textrm{zero}}& M_{02}& M_{3} \end{pmatrix}, \end{aligned}$$
(11)

applicable when \(\lambda _3 = 0\), relating to the open configuration. The corresponding matrix elements are expressed as:

$$\begin{aligned} M_{j\left( j=1,2,3\right) } =\begin{pmatrix} 0& \omega _{mj}& 0& 0\\ -\omega _{mj}& -\gamma _{m}& \sqrt{2}g\textrm{Re}\left( A_{j}\right) & \sqrt{2}g\textrm{Im}\left( A_{j}\right) \\ -\sqrt{2}g\textrm{Im}\left( A_{j}\right) & 0& -\kappa & -F_{j}\\ \sqrt{2}g\textrm{Re}\left( A_{j}\right) & 0& -F_{j}& -\kappa \end{pmatrix}, \end{aligned}$$
(12)

with

$$\begin{aligned} M_{01,02,03}=\begin{pmatrix} 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& \lambda _{1,2,3}\\ 0& 0& -\lambda _{1,2,3}& 0 \end{pmatrix}, \end{aligned}$$
(13)

and

$$\begin{aligned} M_{\textrm{zero}}=\begin{pmatrix} 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0 \end{pmatrix}, \end{aligned}$$
(14)

where \(F_{j}=\Delta _{j}\left[ 1+\eta _{C}\cos \left( \Omega _{C}t\right) \right] +g Q_{j}\).

As proposed by A. Mari et al, the degree of quantum complete synchronization between arbitrary two subsystems can be quantified by the general metric22:

$$\begin{aligned} S_{q_{ij}}\left( t\right) =\frac{1}{\left\langle \delta q_{ij,-}\left( t\right) ^{2}+\delta p_{ij,-}\left( t\right) ^{2}\right\rangle }, \end{aligned}$$
(15)

with \(\delta q_{ij,-}(t) = \left[ \delta q_i(t) - \delta q_j(t) \right] / \sqrt{2}\) and \(\delta p_{ij,-}(t) = \left[ \delta p_i(t) - \delta p_j(t) \right] / \sqrt{2}\) being the related quantum errors. Accordingly, \(S_{qij}(t)\rightarrow 1\) indicates the highest quantum synchronization level, while \(S_{qij}(t)\rightarrow 0\) denotes the lowest level.

Furthermore, for effective assessment of quantum synchronization across multiple subsystems, it would be beneficial to adopt the following system-wide quantum complete synchronization metric, which reflects the system-wide synchronization level of the entire system and offers significant advantages for analyzing quantum synchronization in complex networks48:

$$\begin{aligned} S_{q}\left( t\right) =\frac{C_{n}^{2} }{\left\langle \sum _{i \ne j}\left( \delta q_{ij,-}\left( t\right) ^{2}+\delta p_{ij,-}\left( t\right) ^{2}\right) \right\rangle }, \end{aligned}$$
(16)

where \(C_{n}^{2}\) denotes the binomial coefficient associated with the system size n, while \(\delta q_{ij,-}\) and \(\delta p_{ij,-}\) are defined as the error operators between the i-th and j-th subsystems. Analogous to \(S_{qij}(t)\), the range of \(S_{q}(t)\) is also from 0 to 1, with values nearing 1 indicating higher level of synchronization and vice versa. This metric enables effective evaluation of quantum synchronization across the entire many-body system.

Fig. 2
figure 2

Dynamical evolution of quantum synchronization \(S_{qij}(t)\) and \(S_{q}(t)\) under driving field modulation for (a) the open configuration and (b) the closed configuration. In the unmodulated scenario, quantum synchronization between subsystems are depicted by red dashed line (\(S_{q12}^{0}(t)\)), blue dashed line (\(S_{q23}^{0}(t)\)), and green dashed line (\(S_{q13}^{0}(t)\)), along with the magenta dashed line for the system-wide synchronization (\(S_{q}^{0}(t)\)). In the driving-field modulated scenario, they are represented by red solid line (\(S_{q12}(t)\)), blue solid line (\(S_{q23}(t)\)), and green solid line (\(S_{q13}(t)\)), with the system-wide synchronization (\(S_{q}(t)\)) shown by the magenta solid line. Relevant parameters are: \(\Delta _1 = \omega _{m1} = 0.995\), \(\Delta _2 = \omega _{m2} = 1.000\), \(\Delta _3 = \omega _{m3} = 1.005\), \(\lambda _{1,2} = 0.05\), \(\lambda _{3} = 0 (0.05)\) for open (closed) configuration, \(g = 0.005\), \(\gamma _m = 0.005\), \(\kappa = 0.15\), \(E = 100\), \(\Omega _D = 0(4)\), \(\eta _D = 0(0.8)\) in the unmodulated (driving-field modulated) case, \(\Omega _{C}=0\) and \(\eta _{C}=0\).

Full size image

The numerical calculations of \(S_{qij}\left( t\right)\) and \(S_{q}\left( t\right)\) inevitably involves the quadratic terms \(\delta q_{j}(t)^{2}\), \(p_{j}(t)^{2}\), \(\delta q_{i}\delta q_{j}\left( t\right)\) and \(\delta p_{i}\delta p_{j}\left( t\right)\). Consequently, for the triple-cavity system, we introduce the \(12 \times 12\) covariance matrix as:

$$\begin{aligned} V_{ij}\left( t\right) =\frac{1}{2}\left\langle u_{i}\left( t\right) u_{j}\left( t\right) +u_{j}\left( t\right) u_{i}\left( t\right) \right\rangle , \end{aligned}$$
(17)

which allows us to transform Eq. (7) into the dynamical equation32,58,59:

$$\begin{aligned} \dot{V} \left( t\right) =MV\left( t\right) +V\left( t\right) M+N, \end{aligned}$$
(18)

where \(N=\textrm{diag}\left[ 0,\gamma _{m}\left( 2\bar{n}_{bath}+1\right) ,\kappa ,\kappa ,0,\gamma _{m}\left( 2\bar{n}_{bath}+1\right) ,\kappa ,\kappa ,0,\gamma _{m}\left( 2\bar{n}_{bath}+1\right) ,\kappa ,\kappa \right]\) is the diagonal noise matrix satisfying the correlation relation:

$$\begin{aligned} N_{ij}\delta (t-t') = \frac{1}{2} \langle n_i(t) n_j(t') + n_j(t') n_i(t) \rangle . \end{aligned}$$
(19)

Consequently, both two-body quantum synchronization, \(S_{qij}(t)\), and system-wide quantum synchronization, \(S_q(t)\), can be expressed in terms of \(V_{ij}\) as:

$$\begin{aligned} \begin{aligned} S_{q12}\left( t\right) =&\left\{ \frac{1}{2}\left[ V_{1,1} \left( t\right) +V_{5,5}\left( t\right) +V_{2,2}\left( t\right) +V_{6,6}\left( t\right) -V_{1,5}\left( t\right) -V_{5,1}\left( t\right) -V_{2,6}\left( t\right) \right. \left. -V_{6,2}\left( t\right) \right] \right\} ^{-1},\\ S_{q23}\left( t\right) =&\left\{ \frac{1}{2}\left[ V_{5,5} \left( t\right) +V_{9,9}\left( t\right) +V_{6,6}\left( t\right) +V_{10,10}\left( t\right) -V_{5,9}\left( t\right) -V_{9,5}\left( t\right) -V_{6,10}\left( t\right) \right. \left. -V_{10,6}\left( t\right) \right] \right\} ^{-1},\\ S_{q13}\left( t\right) =&\left\{ \frac{1}{2}\left[ V_{1,1} \left( t\right) +V_{9,9}\left( t\right) +V_{2,2}\left( t\right) +V_{10,10}\left( t\right) -V_{1,9}\left( t\right) -V_{9,1}\left( t\right) -V_{2,10}\left( t\right) \right. \left. -V_{10,2}\left( t\right) \right] \right\} ^{-1},\\ S_{q}\left( t\right) =&\left\{ \frac{1}{6}\left[ V_{1,1} \left( t\right) +V_{5,5}\left( t\right) +V_{2,2}\left( t\right) +V_{6,6}\left( t\right) -V_{1,5}\left( t\right) -V_{5,1}\left( t\right) -V_{2,6}\left( t\right) -V_{6,2}\left( t\right) +V_{5,5} \left( t\right) \right. \right. \\ &\left. \left. +V_{9,9}\left( t\right) +V_{6,6}\left( t\right) +V_{10,10}\left( t\right) -V_{5,9}\left( t\right) -V_{6,10}\left( t\right) -V_{9,5}\left( t\right) -V_{10,6}\left( t\right) +V_{1,1} \left( t\right) +V_{9,9}\left( t\right) \right. \right. \\ &\left. \left. +V_{2,2}\left( t\right) +V_{10,10}\left( t\right) -V_{1,9}\left( t\right) -V_{9,1}\left( t\right) -V_{2,10}\left( t\right) -V_{10,2}\left( t\right) \right] \right\} ^{-1}. \end{aligned} \end{aligned}$$
(20)

By solving Eq. (5), Eq. (7) and Eq. (18), we can straightforwardly derive the evolution of the aforementioned two synchronization metrics. To facilitate the assessment of the impact of modulations on quantum synchronization, we further introduce:

Fig. 3
figure 3

Mean values of the system-wide synchronization \(\bar{S}_{q}(t)\) plotted against driving field modulation amplitude \(\eta _D\) with fixed \(\Omega _D = 4\) for (a) the open configuration and (c) the closed configuration. Similarly, \(\bar{S}_{q}(t)\) are shown as a function of driving field modulation frequency \(\Omega _D\) for fixed amplitudes (b) \(\eta _D = 0.82\) in the open configuration and (d) \(\eta _D = 0.84\) in the closed configuration. Other parameters are as specified in the caption of Fig. 2.

Full size image
$$\begin{aligned} \bar{S}_{q}\left( t\right) =\lim _{T \rightarrow \infty } \frac{1}{T} \int _{t=0}^{T} S_{q}(t) \textrm{d}t \end{aligned}$$
(21)

as the corresponding time-averaged quantum synchronization. This metric allows for an intuitive comparison of quantum synchronization levels under varying modulation parameters.

Results and discussion

Driving modulation

In this section, we explore the potential to achieve high-level quantum synchronization within triple-coupled optomechanical systems. Initially, it is crucial to note that in the absence of modulation, quantum synchronization remains notably low in both open and closed configurations, as evidenced by the dashed curves in Fig. 2(a) and Fig. 2(b).To achieve high-level quantum synchronization within an optomechanical network, our first strategy involves introducing simultaneous periodic modulation across all driving fields.As illustrated in Fig. 2, the pairwise quantum synchronization \(S_{qij}\left( t\right)\), which measures the quantum synchronization between arbitrary two subsystems, along with the system-wide quantum synchronization \(S_{q}\left( t\right)\), which quantifies the quantum synchronization of the entire system, achieve a dynamic equilibrium characterized by minor oscillations following a brief period of evolution. Notably, when reaching the steady state, all metrics exhibit significant improvement over their respective unmodulated baselines. In the open-chain scenario, the steady-state pairwise quantum synchronization values oscillate around average values of \(\bar{S}_{q12}\left( t\right) \approx 0.675\), \(\bar{S}_{q23}\left( t\right) \approx 0.645\) and \(\bar{S}_{q13}\left( t\right) \approx 0.685\), while the system-wide quantum synchronization oscillates around \(\bar{S}_{q}\left( t\right) \approx 0.665\). Comparable enhancements are observed in the closed system scenario, with pairwise quantum synchronization values reaching \(\bar{S}_{q12}\left( t\right) \approx 0.724\), \(\bar{S}_{q23}\left( t\right) \approx 0.661\) and \(\bar{S}_{q13}\left( t\right) \approx 0.696\), and the system-wide quantum synchronization averaging \(\bar{S}_{q}\left( t\right) \approx 0.692\). These results, on one hand, demonstrate that greater uniformity in the synchronization of individual subsystem pairs correlates with enhanced overall system synchronization, and on the other hand, confirm the effectiveness of our modulation strategy for driving fields, which significantly improves quantum synchronization in multi-body optomechanical systems.

Fig. 4
figure 4

Dynamical evolution of quantum synchronization \(S_{qij}(t)\) under cavity mode modulation for (a) the open configuration and (b) the closed configuration. In the unmodulated scenario, quantum synchronization between subsystems are depicted by red dashed line (\(S_{q12}^{0}(t)\)), blue dashed line (\(S_{q23}^{0}(t)\)), and green dashed line (\(S_{q13}^{0}(t)\)), along with the magenta dashed line for the system-wide synchronization (\(S_{q}^{0}(t)\)). In the cavity-mode modulated scenario, they are represented by red solid line (\(S_{q12}(t)\)), blue solid line (\(S_{q23}(t)\)), and green solid line (\(S_{q13}(t)\)), with the system-wide synchronization (\(S_{q}(t)\)) shown by the magenta solid line. Relevant parameters are: \(\Omega _D = 0\), and \(\eta _D = 0\), \(\Omega _{C}=0(3)\) and \(\eta _{C}=0 (2)\) for the unmodulated (cavity mode modulation) case, other parameters are as specified in the caption of Fig. 2.

Full size image

To gain deeper insight of how the driving field modulation can improve the quantum synchronization in triple-coupled optomechanical systems, the mean value of the extended measurements, \(\bar{S}_{q}(t)\), as a function of the driving field modulation parameters \(\Omega _{D}\) or \(\eta _{D}\) are presented in Fig. 3 respectively. As illustrated by Fig. 3(a), \(\bar{S}_{q} \left( t\right)\) exhibits stable plateaus within the \(\eta _{D}\) range of 0.5 to 0.7, and then followed by drastic fluctuations and several discrete peaks. This pattern suggest a region of reliable enhancement of synchronization with such modulation scheme, albeit with heightened sensitivity to changes in modulation amplitude beyond this range. Additionally, as shown in Fig. 3(b), we note that satisfactory synchronization can only be achieved near integer values of \(\Omega _{D}\), with a peak at \(\Omega _{D}=4\) and decreasing towards both extremes. In the closed configuration, as shown in Figs. 3(c) and 3(d), broader peaks of synchronization appear across a more extensive range of \(\eta _{D}\). This behavior suggests that the closed system’s response to driving field modulation is not only robust but also exhibits reduced sensitivity to minor variations in \(\eta _{D}\), likely due to its higher symmetry.

Cavity mode modulation

In the following, we explore the effect of periodic cavity mode modulation on achieving satisfactory quantum synchronization in triple-coupled optomechanical systems. As depicted in Fig. 4, the quantum synchronization dynamics, shown by \(S_{qij}(t)\) for subsystems and \(S_{q}(t)\) for the overall system, initially display transient behaviors that stabilize into minor oscillations. Specifically, as shown in Fig. 4(a), the steady-state pairwise quantum synchronization achieve the average values as \(\bar{S}_{q12}(t) \approx 0.779\), \(\bar{S}_{q23}(t) \approx 0.732\) and \(\bar{S}_{q13}(t) \approx 0.765\), with the average system-wide quantum synchronization \(\bar{S}_{q}(t) \approx 0.758\) in the open configuration. Similarly, as presented in Fig. 4(b), the steady-state quantum synchronization achieve the average values as \(\bar{S}_{q12}(t) \approx 0.805\), \(\bar{S}_{q23}(t) \approx 0.736\), and \(\bar{S}_{q13}(t) \approx 0.772\), with the average system-wide synchronization as \(\bar{S}_{q}(t) \approx 0.770\) in the closed configuration. These results suggest that the periodic cavity modulation, in comparison to driving field modulation, is more effective at enhancing quantum synchronization. Therefore, the strategy of periodic modulation on cavity modes proves to be a highly effective method for improving synchronization in multi-coupled optomechanical networks.

Fig. 5
figure 5

Mean values of the system-wide synchronization \(\bar{S}_{q}(t)\) plotted against cavity mode modulation amplitude \(\eta _C\) with fixed \(\Omega _C = 3\) for (a) the open configuration and (c) the closed configuration. Similarly, \(\bar{S}_{q}(t)\) are shown as a function of the cavity mode modulation frequency \(\Omega _C\) for fixed amplitudes (b) \(\eta _C = 2.5\) in the open configuration and (d) \(\eta _C = 2.4\) in the closed configuration. Other parameters are as specified in the caption of Fig. 4.

Full size image

For a more detailed exploration of how the cavity mode modulation scheme improves the quantum synchronization, we further show in Fig. 5 the time-averaged quantum synchronization \(\bar{S}_{q} \left( t\right)\) against the modulation amplitude \(\eta _{C}\) (frequency \(\Omega _{C}\)) with appropriate selected modulation frequency \(\Omega _{C}\) (amplitude \(\eta _{C}\)). Unlike the narrow platforms and discrete peaks observed in the driving field modulation case, appropriate selected cavity modulation can lead to a higher and more continuous region of quantum synchronization depicted in Fig. 5(a). This continuous region indicates a wider range for parameter choice of the cavity modulation. Additionally, we note from Fig. 5(b) that the optimal quantum synchronization occurs when \(\Omega _{C}\) near an integer, and achieves peak of value \(\bar{S}_{q} \left( t\right) \approx 0.800\) at \(\Omega _{C}=3\). Similar results can also be observed in the closed configuration, as illustrated by the bottom panels in Fig. 5. Therefore, these results suggest that the cavity modulation might be a more stable and effective synchronization enhancement scheme than the driving field modulation in multi-coupled optomechanical networks.

Conclusions

In summary, we have proposed two periodic modulation strategies aimed at enhancing quantum synchronization in two configurations of indirectly coupled optomechanical systems, each comprising three oscillators with distinct frequencies. Numerical simulations demonstrate that satisfactory quantum synchronization is achievable only in the presence of periodic modulation on either driving fields or cavity modes, in both open and closed configurations. Furthermore, it is critical to select appropriate modulation parameters to ensure effective quantum synchronization in both configurations, while modulation on cavity modes appears to be a more stable and effective strategy, owing to its potential for higher levels of quantum synchronization and greater robustness against fluctuations in modulation amplitude. Our approaches could be directly promoted to longer optomechanical oscillator chains or larger networks, and the relevant results should potentially advance the study of quantum synchronization in complex networks.