## 1 Introduction  
Many mechanical and structural systems under active vibration control can be described as linear time-delayed systems (TDSs) due to sensing and actuation delays. The distribution of the eigenvalues (also known as poles or characteristic roots) of such a system on the complex plane determines the stability and dynamic behaviour of the systems. However, the characteristic equation of a TDS contains an infinite number of roots (i.e., eigenvalues) due to its transcendental nature, it is impossible and unnecessary to determine all its roots. For a TDS of retarded type (which makes up the majority of vibration suppression problems of closed-loop control systems), the overall dynamics is mainly dominated by those eigenvalues that lie closely to the imaginary axis, namely, the dominant eigenvalues. There are some ways to numerically determine a part of eigenvalues within a specified region in the complex plane for a retarded or neutral TDS with constant delays [1]. In addition to traditional methods for computing the zeros of analytic functions, e.g., the well-known Newton's method and recently proposed QPmR algorithm [2], two other groups of methods are usually found in the literature. In addition, the Lambert W function can also be used, e.g., [3]. The first group of methods are based on an approximation of the solution map using linear multi-step methods [4; 5] or Runge-Kutta methods [6]. The second group of methods are based on a discretization of the partial differential equation representation of a TDS, using, e.g., Runge-Kutta methods [7], spectral methods [8], or pseudospectral methods [9]. In both groups of methods, the required eigenvalues of a TDS are determined via the eigenvalues of a matrix resulting from the approximation. It is well known that the aforementioned methods are fast and numerically reliable, and they all depend upon the usage of mathematical models, e.g., delay differential equations, of the system. Thus, it is a prerequisite to require the exact knowledge of system matrices such as, for example, mass, damping and stiffness matrices of the system discussed in this paper, which undoubtedly involves errors in relation to practical systems and quite often are not available. For mechanical and structural systems whose equations of motion are naturally formulated in the second-order setting, the second-order time delayed systems are generally analysed in the first-order state-space framework after a transformation. Another interesting and useful scheme still in development is to use nonparametric models, e.g., measured receptances, which was originally developed to design linear vibration control of the systems without delay and required knowledge of the system matrices [10; 11; 12]. There were some works on full/partial eigenvalue assignment of second-order time delayed systems using the receptance method [13; 14; 15; 16; 17]. Additionally, a simple stability criterion for second-order systems with time-varying delay based on the receptance approach was presented in [18]. When computing the eigenvalues of these multidimensional systems with constant time delay and conducting a posteriori stability analysis, however, these literatures [13; 14; 15; 16; 17] also adopted the aforementioned methods that require the availability of the system matrices. Recently, some authors of this paper proposed a reduced characteristic function derived from the original characteristic function of the second-order time delayed systems. This characteristic function is comprised of \"measured\" open-loop receptances without requiring the knowledge of system matrices. A Newton-type method was adopted to compute the dominant eigenvalues of the system via this characteristic function [19]. This paper addresses the computation of dominant eigenvalues for second-order systems with multiple constant time-delays using a contour integral method.