Although these methods are efficient, the implementation of the algorithms tends to be cumbersome because the semi-analytical equations used to obtain the sensitivities change depending on the position of the largest magnitude component of each eigenvector being differentiated. In some cases, these changes result in the solution of ill-conditioned systems of equations [30]. [PERSON] et al. [33] addressed these problems by introducing a general semi-analytical expression to obtain sensitivities that are independent of the characteristics of the eigenvectors. In this approach, only the eigenvalues and eigenvectors of the modes of interest, and the sensitivities for \([M]\) and \([K]\) are needed to calculate the eigenvector sensitivities. This methodology is available for first and second-order derivatives and has been demonstrated in test cases that exhibit distinct and repeated eigenvalues. However, although the semi-analytical expression can be expanded to arbitrary-order, it becomes lengthy, making it difficult to implement within commercial software. In addition, obtaining the sensitivities for \([M]\) and \([K]\) is not trivial and it is prone to error, especially when estimating higher-order sensitivities [41]. The sensitivities for \([M]\) and \([K]\) are calculated either analytically [35; 42], or by using first-order FD approximations [43]. The analytical approach is limited because \([K]\) and \([M]\) are typically estimated numerically. On the other hand, FD presents the same inaccuracy and step size dependency problems highlighted before [44]. Several attempts have been made to develop alternative methodologies to improve the accuracy of the numerical differentiation of the system matrices using FD. For instance, in the works by [PERSON] et al. [45; 46], a machine precision method for numerical differentiation of system matrices was presented. In this method, the first-order forward difference expression was corrected by using a set of scalar and matrix functions. In addition, in [PERSON] et al. [47], correction factors related to the rigid body rotations of the specific finite element formulations were introduced to the FD to calculate machine precision sensitivities. Despite these efforts, the methods are limited to finite element formulations based on Hermite Shape functions and only work for first-order approximations of the sensitivities. As noted in our review, the calculation of higher-order sensitivities in eigenfrequency problems is limited by several technical issues. These include the unavailability of expressions to calculate higher-order sensitivities of the eigenvalues and eigenvectors that are general, computationally efficient, and easy-to-implement; and the absence of a general method to compute high-order derivatives for \([K]\) and \([M]\). These facts have limited the widespread implementation of the current methods into computational practice and commercial software. In this work, these technical challenges are addressed by introducing a novel methodology to obtain arbitrary-order sensitivities for undamped eigenfrequency problems that takes advantage of the hypercomplex automatic differentiation (HYPAD) method. Recently, [PERSON] et al. introduced HYPAD to obtain arbitrary-order sensitivities without truncation errors [48] using multicomplex algebra. HYPAD consists of multiple non-parallel independent imaginary axes used to perform function evaluations. The calculation of high-order and mixed partial sensitivities is obtained by perturbing the parameter of interest along the various imaginary axes, evaluating the function of interest, and extracting the imaginary components of the solution. Furthermore, when the subset of hypercomplex numbers, known as multidual as opposed to multicomplex, are used in HYPAD, no truncation and subtraction cancellation errors are produced in computing the derivatives [49].