## 4 Stability Analysis of Multivariable Milling Tools  
The EMHPM can be generalized for stability analysis of DDEs having multiple delays. A multivariable tool contains some of the following characteristics: uneven pitch between teeth, and/or at least one helix angle with a different value from the others. This analysis was developed by [PERSON] et al., in [34] by using the first-order EMHPM, where the methodology for the characterization of the cutting coefficients for a multivariable tool was discussed, and the dynamic behavior was studied from the productivity point of view. Since the angular spacing at the beginning of the edge is different between teeth (pitch) and the different values of helix angles of the edges between adjacent teeth, the angular spacing between teeth at a specific height changes continuously, which produces an infinite number of delays. A common approach to deal with the DDE with an infinite number of delay is to discretize the tool by cutting disks in the axial direction with a thickness \(\Delta a_{disk}\) to induce a DDE with a finite number of delays. A single disk still has the same number of flutes (discrete flutes) and considering that the maximum delay in the process is the period of rotation of the tool or the spindle rotation period \(\tau_{T}\), then, it can be discretized in \(N-1\) intervals.  
The angular position between two adjacent teeth in each cutting disk changes according to the axial position of the referred disk and is related to the expression \(\psi=k_{\beta}a_{p}\), where \(k_{\beta}=2\tan\beta/2D\). Here \(D\) is the diameter of the tool and \(\psi\) represents the cutting-edge offset angle due to the helix angle. A certain interval can be associated with a discrete time delay of each tooth _iz_ and disk \(l\) using the following formulation  
\[N_{iz,l}=round\bigg{(}(N-1)\frac{\delta\phi_{iz,l}}{2\pi}\bigg{)} \tag{27}\]  
where \(\delta\phi_{iz,l}\) is the angular pitch between consecutive teeth for each disk, the round function converts the argument to the nearest integer. In Equation (27) \(N_{iz}\), is a table (matrix) of dimension \(iz\times l\). Since this procedure could generate several delayed terms and some of them with the same value of discrete time delay due to the discretization scheme, it is required to collect all the different (non-repeated) discrete time delays \(d_{n}\) from \(N_{iz,l}\).  
Thus, without loss of generality, the DDE with multiple delays can be written as  
\[\dot{\mathbf{x}}(t)=\mathbf{A}(t)\mathbf{x}(t)+\sum_{d=\min(d_{n})}^{\max(d_{ n})}\mathbf{B}^{d}\mathbf{x}(t-\tau) \tag{28}\]  
\begin{table}
\begin{tabular}{c c c c c} \hline \hline \multirow{2}{*}{**Method**} & \multicolumn{2}{c}{**\(N\) = 60**} & \multicolumn{2}{c}{**\(N\) = 100**} \\  & **Time (s)** & **Error** & **Time (s)** & **Error** \\ \hline First-order EMHPM & 52.39 & 1.6276 & 186.35 & 0.3203 \\ Second-order EMHPM & 54.85 & 0.3507 & 194.65 & 0.0568 \\ Third-order EMHPM & 61.48 & 0.3458 & 195.86 & 0.0563 \\ FDM & 61.26 & 1.4097 & 191.49 & 0.3998 \\ SDM & 278.03 & - & 570.04 & - \\ \hline \hline \end{tabular}
\end{table}
Table 1: Comparison of convergence for different methods for down-milling operation with \(a_{d}=1\)where \(\mathbf{x}\) is the vector of states, \(\mathbf{A}(t+\tau_{T})=\mathbf{A}(\tau_{T})\), \(\mathbf{B}^{d}(t+\tau_{T})=\mathbf{B}^{d}(\tau_{T})\) and \(\tau_{T}\) is the period of rotation of the spindle. Following the EMHPM procedure, Equation (28) can be written equivalently by intervals as:  
\[\dot{\mathbf{x}}_{i}(T)-\mathbf{A}_{t}\mathbf{x}_{i}(T)\approx\sum_{d=\min(d_ {n})}^{\max(d_{n})}\mathbf{B}_{t}^{d}\mathbf{x}_{i}^{\tau_{d}}(T) \tag{29}\]  
being \(\mathbf{x}_{i}(T)\) the solution by intervals of order m for Equation (28) in the \(i-th\) interval that satisfies the initial condition \(\mathbf{x}_{i}(0)=\mathbf{x}_{i-1}\), the matrices \(\mathbf{A_{t}}\) and \(\mathbf{B}_{t}^{d}\) represent the values of the matrices \(\mathbf{A}(t)\) and \(\mathbf{B}^{d}(t)\) evaluated at time \(t\) respectively.