## 3 Numeric Solution of the Milling Equation  
### Dynamic Model to the Milling Equation  
To validate the proposed EMHPM methods, the numerical solution of the delay differential equation analyzed by [PERSON] et al., in [PERSON] et al. (2018) was calculated, which describes the dynamic model of the milling process in one degree of freedom (DOF):  
\[\ddot{\mathbf{x}}(t)+2\zeta\omega_{n}\dot{\mathbf{x}}(t)+\omega_{n}^{2} \mathbf{x}(t)=-\frac{a_{p}h_{xx}(t)}{m_{m}}(\mathbf{x}(t)-\mathbf{x}(t-\tau)) \tag{23}\]  
where \(\zeta\) is the modal damping ratio, \(\omega_{n}\) is the natural frequency of the workpiece, \(a_{p}\) is the axial depth of cut, \(m_{m}\) is the modal mass, \(\tau\) represents the time delay corresponding to the hitting period between each tooth of the tool and \(h_{xx}(t)\) is the specific cutting force in the \(x\)-direction due to flexibility in \(x\)-direction, which was calculated depending on the position of the tool  
\[h_{xx}(t)=\sum_{iz=1}^{z_{n}}g(\phi_{iz}(t))sin\phi_{iz}(t)(K_{tc}\cos\phi_{iz }(t)+K_{nc}\sin\phi_{iz}(t)) \tag{24}\]  
\(z_{n}\) is the number of edges of the tool, \(K_{tc}\) and \(K_{nc}\) are the average specific cut coefficients in the tangential and normal direction, respectively, and \(\phi_{iz}(t)\) is the angular position of each left edge described by  
\[\phi_{iz}=(2\pi n/60)t+2\pi iz/z_{n} \tag{25}\]  
where \(n\) is the spindle speed in revolution per minute (rpm). The function \(g(\phi_{iz}(t))\) is a window function, which has the value of one when the current edge \(iz\) is cutting material, otherwise it takes the value zero.  
In up-milling \(\phi_{st}=0\) and \(\phi_{cx}=\cos^{-1}(1-2a_{d})\), conversely in down-milling \(\phi_{st}=\cos^{-1}(2a_{d}-1)\) and \(\phi_{cx}=\pi\), \(a_{d}\) is the radial immersion ratio of the cut and \(\phi_{st}\) and \(\phi_{cx}\) are the angular positions where each edge enters and leaves the workpiece.  
The second- and third-order EMHPM is applied to obtain the solution of Equation (23) and it is compared with the solution given by the first-order EMHPM [PERSON] et al. (2018); for a regular tool the matrix \(\mathbf{A}\) and \(\mathbf{B}\) are represented as:  
\[\mathbf{A}_{t}=\left[\begin{array}{cc}0&1\\ -\omega_{n}^{2}-\frac{a_{p}h_{xx}(t)}{m_{m}}&-2\zeta\omega_{n}\end{array} \right],\ \mathbf{B}_{t}=\left[\begin{array}{cc}0&0\\ \frac{a_{p}h_{xx}(t)}{m_{m}}&0\end{array}\right] \tag{26}\]  
\(\mathbf{A}_{t}\) and \(\mathbf{B}_{t}\) correspond to the periodic matrix evaluated at time \(t\). For demonstration purposes, time-domain simulations were computed for a full-immersion down-milling operation. We used the parameters employed by [PERSON] et al., in In [PERSON] et al. (2018) where the stability lobes were also calculated.  
The modal parameters \(f_{n}=922\) Hz, \(\omega_{n}=5793\) rad/s, \(\zeta=0.011\) and \(m_{m}=0.03993\) kg corresponds to a single degree of freedom. The tangential and normal cutting coefficients are \(K_{tc}=6\times 10^{8}\) N/m\({}^{2}\) and \(K_{nc}=2\times 10^{8}\) N/m\({}^{2}\) respectively for an end-mill with \(z_{n}=2\). The time-domain solution was computed using the EMHPM considering \(N=76\) discrete intervals and \(m=7\). Two sets of cutting conditions were chosen for a fixed spindle speed value of \(n=12\),000 rpm where the axial depth of cut of \(a_{p}=1.5\) mm corresponds to a stable cutting operation while that for an unstable operation \(a_{p}=3\) mm was chosen. In Figure 2 we plot the second- and third-order EMHPM solutions and compare it with the first-order EMHPM and the dde23 routine in Matlab, which is used to integrate DDE.