## 3 Methodology
### HB-AFT Method for Obtaining VC Period-\(M\) Solution  
To improve the accuracy and efficiency of obtaining the VC period-\(M\) solution, one can introduce a dimensionless displacement \(\mathbf{X}=[X,Y]^{\mathrm{T}}=[x/\delta_{0},y/\delta_{0}]^{\mathrm{T}}\), and dimensionless time \(\tau=N_{b}\cdot\Omega\cdot t/\mathrm{M}\); then, Equation (5) can be expressed as:  
\[\mathbf{X}^{\prime\prime}(\tau)=\mathbf{F}(\mathbf{X}(\tau),\mathbf{X}^{\prime}(\tau),\lambda) \tag{12}\]  
Here, let the control parameter \(\lambda=\Omega\) for tracing the amplitude-frequency response curves, and for System (12) the period \(T_{VC}\) of the VC period-\(M\) solution is \(2\pi\). Then, \(\mathbf{X}(\tau)\) and \(\mathbf{F}(\cdot)\) can be expressed as:  
\[\left[\begin{array}{c}\mathbf{X}\\ \mathbf{F}\end{array}\right]=\left[\begin{array}{c}\mathbf{a}_{0}\\ \mathbf{c}_{0}\end{array}\right]+\sum_{k=1}^{K}\left[\begin{array}{c}\mathbf{a}_{k} \\ \mathbf{c}_{k}\end{array}\right]\cos(k\tau)-\left[\begin{array}{c}\mathbf{b}_{k}\\ \mathbf{d}_{k}\end{array}\right]\sin(k\tau)\right] \tag{13}\]  
According to the process of harmonic balance, insert Equation (13) into Equation (12) and obtain the following \(4K+2\) algebraic relationships:  
\[g(\mathbf{P},\mathbf{Q},\lambda)=0 \tag{14}\]  
where the harmonic coefficients of \(\mathbf{X}(\tau)\) and \(\mathbf{F}(\cdot)\) are:  
\[\left[\begin{array}{c}\mathbf{P}\\ \mathbf{Q}\end{array}\right]^{\mathrm{T}}=\left[\begin{array}{cc}\mathbf{a}_{0},\ \mathbf{a}_{1},\ \mathbf{b}_{1},\ \mathbf{a}_{2},\ \mathbf{b}_{2},\ \ldots,\ \mathbf{a}_{K},\ \mathbf{b}_{K}\\ \mathbf{c}_{0},\ \mathbf{c}_{1},\ \mathbf{d}_{1},\ \mathbf{c}_{2},\ \mathbf{d}_{2},\ \ldots,\ \mathbf{c}_{K},\ \mathbf{d}_{K}\end{array}\right]^{ \mathrm{T}} \tag{15}\]  
To obtain \(\mathbf{P}\) from Equation (14), it should first formulate \(\mathbf{Q}\) by \(\mathbf{P}\), because Equation (14) contains \(8K+4\) unknown harmonic coefficients of \(\mathbf{P}\) and \(\mathbf{Q}\), and this process can be completed by the following AFT technology.  
Figure 1: Schematic diagram of (**a**) the ball bearing system, and (**b**) its two degrees-of-freedom spring model.  
The discrete sample of \(\mathbf{X}(\tau)\) is given by the inverse discrete Fourier transform as:  
\[\mathbf{X}(n)=\text{Real}[\sum_{k=0}^{K}P_{k}e^{i(2\pi\text{k}n/N)}] \tag{16}\]  
and the discrete time domain sample of Equation (12) is:  
\[\mathbf{F}(n)=\mathbf{F}(\mathbf{X}(n),\mathbf{X}^{\prime}(n),\lambda) \tag{17}\]  
Then, by discrete Fourier transform, \(\mathbf{Q}\) can be expressed as:  
\[\mathbf{Q}_{k}=\frac{\phi}{N}\!\sum_{n=0}^{N-1}\mathbf{F}(n)\ e^{i(-2\pi\text{k}n/N)} \tag{18}\]  
Herein, \(\mathbf{P}_{k}=\mathbf{a}_{k}+i\ \mathbf{b}_{k}\); \(\mathbf{Q}_{k}=\mathbf{c}_{k}+i\ d_{k}\); \(N\) is the numbers of samplings in the time domain, and \(n=0,\ \dots,\ N-1\); \(\phi=1\) if \(n=0\), or else \(\phi=2\).  
At this moment, \(\mathbf{Q}\) is formulated by \(\mathbf{P}\) as shown in Equation (18). Taking \(\mathbf{P}\) as an unknown variable, one can obtain the harmonic solution of \(\mathbf{X}(\tau)\) by solving the nonlinear algebraic Equation (14). In order to avoid the singular case at the turning point and track \(\mathbf{P}\) automatically, the arc-length continuation is introduced by transforming Equation (14) into following form [22]:  
\[\left\{\begin{array}{l}\frac{dg}{d\mathbf{P}}\frac{d\mathbf{P}}{d\mathbf{S}}+\frac{dg}{ d\lambda}\frac{d\lambda}{d\lambda}=0\\ \frac{dg}{ds}\frac{d\mathbf{P}}{ds}+\left(\frac{d\lambda}{d\lambda}\right)^{2}=1 \end{array}\right. \tag{19}\]  
where \(s\) denotes the arc-length parameter of curve (\(\mathbf{P}(s)\), \(\lambda(s)\)).  
Finally, Equation (19) can be solved by the predictor-corrector procedure [50], where the value of the initial condition \(\left(\mathbf{P}(s_{0}),\lambda(s_{0})\right)=\left(\mathbf{P},\lambda\right)_{0}\) is searched for using the HB-AFT process above.