## 3 Mathematical Investigations
### Stability Analysis via Linearizing the above System  
The equilibrium solution's stability was evaluated by analyzing the eigenvalues of the Jacobian matrix presented in Equations (27)-(30). Stability is assessed by focusing on the real components of these eigenvalues: if all eigenvalues exhibit a negative real part, the equilibrium solution is deemed asymptotically stable. Conversely, if any eigenvalue possesses a positive real component, the equilibrium solution is regarded as unstable. To provide further clarity on stability, the analysis initially involves examining small deviations from steady-state solutions to determine the stability criteria. This examination enables the understanding of the system's response to small disturbances around the equilibrium \(a_{10},a_{20},\phi_{10}\) and \(\phi_{20}\). Thus, we assume that  
\[\left.\begin{array}{ll}a_{1}=a_{11}+a_{10},&a_{2}=a_{21}+a_{20},\\ \phi_{1}=\phi_{11}+\phi_{10},&\phi_{2}=\phi_{21}+\phi_{20},\\ \dot{a}_{1}=\dot{a}_{11},&\dot{a}_{2}=\dot{a}_{21},\dot{\phi}_{1}=\dot{\phi}_{ 11},\\ \dot{\phi}_{2}=\dot{\phi}_{21}.\end{array}\right\} \tag{43}\]  
where \(a_{10},a_{20},\phi_{10}\) and \(\phi_{20}\) satisfy (35) and (42). Also, \(a_{11},a_{21},\phi_{11}\) and \(\phi_{21}\) are perturbations assumed to be small compared to \(a_{10},a_{20},\phi_{10}\) and \(\phi_{20}\). Substituting (43) into (27)-(30), expanding for small \(a_{11},a_{21},\phi_{11}\) and \(\phi_{21}\), and keeping linear terms in \(a_{11},a_{21},\phi_{11}\) and \(\phi_{21}\), we obtain  
\[\dot{a}_{11}=r_{11}a_{11}+r_{12}\phi_{11}+r_{13}a_{21}+r_{14}\phi_{21}, \tag{44}\]  
\[\dot{\phi}_{11}=r_{21}a_{11}+r_{22}\phi_{11}+r_{23}a_{21}+r_{24}\phi_{21}, \tag{45}\]  
\[\dot{a}_{21}=r_{31}a_{11}+r_{32}\phi_{11}+r_{33}a_{21}+r_{34}\phi_{21}, \tag{46}\]  
\[\dot{\phi}_{21}=r_{41}a_{11}+r_{42}\phi_{11}+r_{43}a_{21}+r_{44}\phi_{21}. \tag{47}\]  
where \(r_{ij},i=1,2,3,4\) and \(j=1,2,3,4\) are given in the Appendix A.  
Equations (43)-(46) can be presented in the following matrix form:  
\[\left[\dot{a}_{11}\;\dot{\phi}_{11}\dot{a}_{21}\;\dot{\phi}_{21}\right]^{T}=[ J][a_{11}\;\phi_{11}a_{21}\;\phi_{21}]^{T}, \tag{48}\]  
\[[J]=\left[\begin{array}{c}r_{11}\;r_{12}\;r_{13}\;r_{14}\\ r_{21}\;r_{22}\;r_{23}\;r_{24}\\ r_{31}\;r_{32}\;r_{33}\;r_{34}\\ r_{41}\;r_{42}\;r_{43}\;r_{44}\end{array}\right]. \tag{49}\]  
where \([J]\) is the Jacobian matrix.  
Therefore, the stability of the steady-state solutions is determined by examining the eigenvalues of the Jacobian matrix. This analysis yields the following eigenvalue equation:  
\[\left|\begin{array}{cccc}r_{11}-\lambda&r_{12}&r_{13}&r_{14}\\ r_{21}&r_{22}-\lambda&r_{23}&r_{24}\\ r_{31}&r_{32}&r_{33}-\lambda&r_{34}\\ r_{41}&r_{42}&r_{43}&r_{44}-\lambda\end{array}\right|=0, \tag{50}\]  
where the above polynomial has the following roots:  
\[\lambda^{4}+\Gamma_{1}\lambda^{3}+\Gamma_{2}\lambda^{2}+\Gamma_{3}\lambda+ \Gamma_{4}=0, \tag{51}\]  
where \(\Gamma_{i};\;(i=1,\ldots,4)\) are the coefficients of Equation (50) defined in the appendix. The Routh-Hurwitz criterion must be satisfied for the abovementioned system's solution to be stable:  
\[\Gamma_{1}>0,\;\Gamma_{1}\Gamma_{2}-\Gamma_{3}>0,\Gamma_{3}(\Gamma_{1}\Gamma _{2}-\Gamma_{3})-\Gamma_{1}^{2}\Gamma_{4}>0,\Gamma_{4}>0. \tag{52}\]