## 3 Perturbation Analysis and Modulation Equations  
In line with the multiple scales procedure outlined in [45; 46], we expand the unknown displacement coordinates \(q_{k}\) and velocity coordinates \(z_{k}\) (\(k=v,w\)) as follows  
\[q_{k}(t;\varepsilon)=\varepsilon q_{k1}(T_{0},T_{1},T_{2})+\varepsilon^{2}q_{ k2}(T_{0},T_{1},T_{2})+\varepsilon^{3}q_{k3}(T_{0},T_{1},T_{2})+..., \tag{7}\]  
where \(T_{i}\)=\(\varepsilon^{i}t\), \(\partial/\partial t\)=\(\sum_{i=0}^{\infty}\varepsilon^{i}D_{i}\), and \(D_{i}\)=\(\partial/\partial T_{i}\), (\(i=0,1,2\)). By inserting Equation (7) and the time-derivative expansion into the governing equations and equating coefficients of equal powers of \(\varepsilon\), we obtain the following six first-order differential equations of order \(-\varepsilon^{1}\), \(-\varepsilon^{2}\), and \(-\varepsilon^{3}\)  
Order \(\varepsilon^{1}\)  
\[ \begin{split} D_{0}^{2}q_{v1}+\omega_{v}^{2}q_{v1}=0,\\ D_{0}^{2}q_{w1}+\omega_{w1}^{2}q_{w1}=0,\end{split}  \tag{8}\]  
Order \(\varepsilon^{2}\)  
\[ \begin{split} D_{0}^{2}q_{v2}+\omega_{v}^{2}q_{v2}& =-2D_{0}D_{1}q_{v1}+\Lambda_{v1}q_{v1}^{2}+\Lambda_{v2}q_{w1}^{2},\\ D_{0}^{2}q_{w2}+\omega_{w}^{2}q_{w2}&=-2D_{0}D_{1}q _{w1}+\Lambda_{w1}q_{v1}q_{w1},\end{split}  \tag{9}\]  
Order \(\varepsilon^{3}\)  
\[ \begin{split} D_{0}^{2}q_{v3}+\omega_{v}^{2}q_{v3}=& -\left(D_{1}^{2}+2D_{0}D_{2}+2\mu_{v}D_{0}\right)q_{v1}-2D_{0}D_{1}q _{v2}\\ &+2\Lambda_{v1}q_{v1}q_{v2}+2\Lambda_{v2}q_{w1}q_{w2}+\Gamma_{v1 }q_{v1}^{3}+\Gamma_{v2}q_{v1}q_{w1}^{2}+\delta_{kv}f_{v}\cos(\Omega_{v}t),\\ D_{0}^{2}q_{w3}+\omega_{w}^{2}q_{w3}=&-\left(D_{1} ^{2}+2D_{0}D_{2}+2\mu_{w}D_{0}\right)q_{w1}-2D_{0}D_{1}q_{w2}\\ &+\Lambda_{w1}(q_{v1}q_{w2}+q_{v2}q_{w1})+\Gamma_{w1}q_{v1}^{2}q _{w1}+\Gamma_{w2}q_{w1}^{3}+\delta_{kv}f_{w}\cos(\Omega_{w}t),\end{split}  \tag{10}\]  
where \(\delta_{kv}\) and \(\delta_{kv}\) (\(k=v,w\)) are the Kronecker delta functions. One begins by obtaining solutions to the order \(-\varepsilon^{1}\) equations, which can be derived as follows  
\[ \begin{split} q_{v1}&=A_{v}(T_{1},T_{2})\mathrm{e}^{ \mathrm{i}\omega_{v}T_{0}}+c.c.,\\ q_{w1}&=A_{w}(T_{1},T_{2})\mathrm{e}^{\mathrm{i} \omega_{v}T_{0}}+c.c.,\end{split}  \tag{11}\]  
where \(c.c.\) represents the conjugate for the preceding complex terms, and \(A_{k}(k=v,w)\) is the complex valued amplitude of the \(k-\)th mode. Next, we substitute Equation (11) into the order \(-\varepsilon^{2}\) equations and follow the same procedure outlined in [45] to obtain the corresponding second-order solutions\[q_{v2} =-\frac{\Lambda_{v1}}{3\omega_{v}^{2}}A_{v}^{2}\mathrm{e}^{2\mathrm{i }\omega_{v}T_{0}}+\frac{\Lambda_{v1}}{\omega_{v}^{2}}A_{v}\bar{A}_{v}-\frac{ \Lambda_{v2}}{3\omega_{v}^{2}}A_{w}^{2}\mathrm{e}^{2\mathrm{i}\omega_{w}T_{0}}+ \frac{\Lambda_{v2}}{\omega_{w}^{2}}A_{w}\bar{A}_{w}+c.c., \tag{12}\] \[q_{w2} =\frac{\Lambda_{w1}}{\omega_{w}^{2}-(\omega_{v}+\omega_{w})^{2}}A_ {v}A_{w}\mathrm{e}^{\mathrm{i}(\omega_{v}+\omega_{w})t}+\frac{\Lambda_{w1}}{ \omega_{w}^{2}-(\omega_{v}-\omega_{w})^{2}}A_{v}\bar{A}_{w}\mathrm{e}^{ \mathrm{i}(\omega_{v}-\omega_{w})t}+c.c..\]  
To quantify the proximity of external and internal resonances, we introduce two detuning parameters denoted by \(\sigma_{1}\) and \(\sigma_{2}\). \[\Omega_{k} =\omega_{k}+\varepsilon^{2}\sigma_{1},\ \ (k=v,w) \tag{13}\] \[\omega_{w} =\omega_{v}+\varepsilon^{2}\sigma_{2},\]  
where \(\omega_{v}\) and \(\omega_{w}\) are the in-plane and out-of-plane mode frequencies, respectively. \(\Omega_{k}\) is the excitation frequency which is applied on the in-plane or out-of-plane modes directly. \(\sigma_{1}\) and \(\sigma_{2}\) are the external and internal detuning parameters, respectively. By substituting the previous second-order solutions into the order \(-\varepsilon^{3}\) equations and introducing solvability conditions, we obtain the modulation equations governing \(A_{v}\) and \(A_{w}\).