Two-parameter dynamics of an autonomous mechanical governor system with time delay

A deep understanding of the dynamical behavior in the parameter-state space plays a vital role in both the optimal design and motion control of mechanical systems. By combining the GPU parallel computing technique with two determinate indicators, namely the Lyapunov exponents and Poincaré section, this paper presents a detailed study on the two-parameter dynamics of a mechanical governor system with different time delays. By identifying different responses in the two-parameter plane, the effect of time delay on the complexity of the evolutionary process is fully revealed. The path-following calculation scheme and time domain collocation method are used to explore the detailed bifurcation mechanisms. An interesting phenomenon that the number of intersection points of some periodic responses on the specified Poincaré section differs from the actual period characteristics is found in classifying the dynamic behavior. For example, the commonly exhibited period-one orbit may have two or more intersection points on the Poincaré section rather than one point. The variations of the basins of attraction are also discussed in the plane of initial history conditions to demonstrate the multistability phenomena and chaotic transitions.


Introduction
Time delay effect is widely encountered in many dynamical systems including biological systems [1][2][3], social and economic systems [4,5], neural networks [6,7], chemical reactions [8,9], and mechanical systems [10,11], due to the inherent dependency of present states on past states. Compared with the systems without time delay, the systems involving time delays can exhibit more complicated dynamical phenomena, including the multiple stability switches [12,13], high co-dimension bifurcations [14][15][16], coexisting attractors [17,18], and chaotic transitions [19,20]. In addition to the understanding of delay-induced bifurcations and transitions, the beneficial applications of the time delay effect have also attracted great interest from the research community. For example, the time delay effect has been considered in the design of real-world controllers, due to the unavoidable time consumption in the measurement, signal processing, computing, and transmission processes. The inherent time delays could be used to achieve the desired control objectives, such as for the traffic arrangement [21,22], signal enhancement [23], broadband energy harvesting [24,25], and absorber and isolator design [26,27].
The mechanical centrifugal governor has played a role of milestone in the field of automatic feedback control, since it was invented to maintain the working speed of steam engines. Owing to the simple yet effective characteristics, the centrifugal governor devices have been widely used to prevent potential damage in many rotating machines, such as in the internal combustion engines, fueled turbines, modern striking clocks, and so on. Induced by the inherent nonlinearities, many dynamical phenomena have been observed and reported in the literature, including periodic oscillations [28,29], chaotic motions [30,31], hyperchaotic behavior [32], multistability [15,33], and fractal and mode-locking structures [34,35]. To avoid the damages caused by the undesired system response, various control schemes have been proposed to stabilize the system. For example, the feedback controller [36,37], adaptive controller [37], and time-delayed controller [38] were used to convert the chaotic motions of the governor systems into specified motions. A synchronizationbased adaptive sliding mode controller [39] and a robust adaptive controller [40] were designed to realize the finite-time stabilization of mechanical governor systems with uncertainties. The design of anti-bifurcation controllers [41,42] was presented for a centrifugal governor system to generate desirable oscillations which might be applied in fault diagnosis, monitoring, mixing, and low-energy navigation.
Although remarkable achievements have been made to the understanding of dynamical behavior in the literature, few research works have been devoted to the dynamics of such governor systems in the presence of inherent time delay effect. Actually, time delay can be originated from the gear transmission, spring deformation [37], steam transportation, or chemical combustion [15]. The reaction delay caused by the mass-spring structure was considered in the governor system, and a variety of dynamical responses was obtained by numerical simulations [37]. The reaction delay induced by the working principle of engines was also considered, and the related time delay effect on the system response was revealed in Ref. [15]. It was found that the time delay can shrink the stability region of equilibrium and even induce novel bifurcations, e.g., resonant and non-resonant double Hopf bifurcations [15]. However, the global dynamics and transitions are not fully understood and yet to be studied in the two-parameter region, where the chosen equilibrium loses its stability. To develop an in-depth insight into the open questions unaddressed in Ref. [15], this paper aims to investigate the global dynamics of the autonomous governor system with and without time delay. That is to say, this paper will present detailed studies on the two-parameter classification of system response and the related multistability phenomena, thereby developing a deep understanding of the delay-induced global dynamical behavior of the present governor system.
Notably, the existing studies on the nonlinear dynamics were mainly focused on the local bifurcation analysis during past decades, which only provides a limited view of the system performance. With the development of computer technology, it has been possible for researchers to study a wide range of bifurcation behaviors of certain systems in the parameter and state spaces through extensive numerical simulations. It is worth mentioning that an efficient accelerated tool, i.e., the CUDA-based GPU parallel computing technique [43], has been rapidly developed for the computation of multi-parameter dynamics of the systems and the efficiency of such a numerical tool has been proved in many works [44][45][46]. However, these research works are mainly focused on the dynamical systems described by the ordinary differential equations, and few research works are associated with the multi-parameter dynamics computation for the dynamical systems governed by the delayed differential equations (i.e., functional differential equations). In this paper, the CUDA-based GPU parallel computing technique is adopted to study the classification of two-parameter dynamics in a time-delayed governor system. Therefore, this work can be regarded as an important attempt in the computation of the multi-parameter dynamics of the delayed dynamical system.
Given that the existing numerical methods are primarily for the ordinary differential dynamical systems, the treatment of the delayed terms and approximating the delayed differential dynamical system by the ordinary differential dynamical system then become a noteworthy challenge in the computation of global dynamics of the delayed systems. Currently, the efficient approximation schemes for the delayed dynamical systems include the Galerkin projections [47] and the discretization approaches [48]. Based on these approximation schemes, the delayed dynamical system can be approximated by a high-dimensional ordinary differential system. Significant time consumption will be involved in precisely obtaining the portfolios of global dynamics of such a high-dimensional system in the parameter and state spaces. In addition to the high computation time cost on the global dynamics, the delay-induced winding phenomenon in the identification of bifurcation will be another challenge. As reported in Ref. [49], the selfintersection phenomenon of phase trajectory can be delay-induced in the autonomous and delayed dynamical system, which leads to difficulties in distinguishing periodic response and identifying the bifurcation mechanism.
For the present autonomous and delayed governor system, the winding phenomena of periodic responses are observed for certain values of time delays, and thus, the bifurcation identification becomes a difficult task. To address this problem, a type of semianalytical method, i.e., the time domain collocation method (TDCM) [50][51][52], is adopted to determine the bifurcation types of periodic response for its ease of implementation without the redundant process of formula derivation. The effectiveness of this method has been fully proved in the numerical bifurcation analysis of delayed dynamical systems and the corresponding software packages, e.g., DDE-Biftool [53] and COCO [54,55]. Therefore, the TDCM is chosen to obtain the approximate analytical solution of periodic response and identify the corresponding bifurcations in this study.
The main contributions of this work can be briefly summarized in three aspects. Firstly, although many research works have focused on in the global computations of dynamical systems governed by ordinary differential equations, few works are involved with the multi-parameter simulation of time-delayed dynamics. By combining some indicators with the existing numerical techniques, this paper characterizes the two-parameter dynamics distribution of a type of timedelayed dynamical system. The presented numerical schemes can also be applied to other time-delayed dynamical systems. Secondly, for the mechanical governor system, the study of global dynamics of the system in both the parameter and state spaces by considering the inherent time delay can provide a better understanding of the system performance and avoid some undesired behaviors, e.g., sudden jump phenomena and chaotic oscillations. Thirdly, a novel phenomenon, i.e., the winding phenomenon, is observed in the bifurcation identification of adjacent periodic responses. This phenomenon might cause failure in determining the periodic property of periodic response by using the common Poincaré mapping method.
The rest of this paper is organized as follows. The governing equations of the autonomous mechanical governor system with time delay are briefly introduced in Sect. 2. The methodologies used for computation and bifurcation identification of two-parameter dynamics of the mechanical governor system are presented in Sect. 3, including the numerical integration strategy of time-delayed system, the determinate indicators of dynamical response, GPU parallel computing technique and bifurcation identification methods of periodic response. Based on the abovementioned methods, the classification of the global dynamical behavior of the mechanical governor system with and without time delay is derived and discussed in Sect. 4. In addition, the detailed analyses on the observed dynamical phenomena are also presented in Sect. 4. Finally, the conclusions are drawn in Sect. 5.

Revisiting of the time-delayed mechanical governor model
A type of mechanical governor system with reaction time delay generated from the working principle of engines is illustrated in Fig. 1. When the rotating machine is subjected to an unexpected external disturbance, the angular velocity X 2 of the flywheel will change accordingly. To maintain the nearconstant speed operation of the flywheel, an automatic feedback controller is designed to regulate the centrifugal motion of the balls with mass m. According to the Lagrangian method, the governing equation of the balls can be described as follows 2 ml 2 d 2 W dt 02 À mX 2 1 r a l cos W where W represents the angular displacement of the balls. The notations and physical meanings of the parameters involved in Eq. (1) are described in Table 1. Furthermore, by applying the momentum theorem, the motion governing equation of the flywheel can also be derived as: where X 2 is the angular velocity of the flywheel. The notation W a represents Wðt 0 À aÞ, and a is a positive parameter denoting the delayed value of the angular displacement of the balls. The physical meanings of the parameters appeared in Eq. (2) are also given in Table 1. A more detailed modeling process can be found in Ref. [15]. For the convenience of analysis, rescaling operations of time is performed in Eqs. (1) and (2) with t ¼ X n t 0 ; s ¼ X n a; X n ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ð2kl þ mgÞ=ml p . By letting x ¼ W; y ¼ dW=dt, and z ¼ X 2 , the governing equations of the mechanical governor system can be expressed as Here, the six intermediate parameters d; e; p; b; q; Q are given by  Equivalent torque of the load a Propagation delay Fig. 1 The physical model of the mechanical governor system d ¼ For the time-delayed governor system (3), one equilibrium can be easily determined as: By adopting the coordinate transformations where a 0 ¼ ðpz Ã2 þ eÞ=2; a 1 ¼ pz Ã ; a 2 ¼ dz Ã2 ; a 3 ¼ 2dz Ã , and the equilibrium is translated into the origin, i.e., Oð0; 0; 0Þ.
For the mechanical governor system (5), the linear stability, local bifurcation behavior, and the time delay effect have been well investigated in a two-parameter plane, i.e., the ðb; QÞ plane in Ref. [15]. As an extension of Ref. [15], this paper will investigate the global behavior of the time-delayed system and examine the influence of time delay in the same twoparameter plane. More specifically, the system parameters are identical to those set in Ref. [15], i.e., p ¼ 0:04; d ¼ 0:09; e ¼ 0:4; q ¼ 2. The ðb; QÞ parameters are chosen as the primary bifurcation parameters with the ranges of b 2 ½0; 2; Q 2 ð0; 2. The time delay is selected as an auxiliary bifurcation parameter and a set of discrete values, i.e., s ¼ 0; 1; 2; 3; 4:0637988; 5 is considered for studying time delay effects. It should be mentioned that the 1 : 4 resonant double Hopf bifurcation of equilibrium arises with s ¼ 4:0637988, and the detailed analysis of this bifurcation can be found in Ref. [15].

Analysis methodologies
In this section, the numerical schemes and theoretical methods used for the computation and bifurcation identification of the global dynamics of the timedelayed dynamics are briefly introduced. To obtain the distribution diagram of two-parameter dynamics of the time-delayed system, a numerical integration algorithm, classification indicators of dynamical response, and the accelerating technique will be described in Sects. 3.1-3.3, respectively. To analyze the transition mechanism of the adjacent motions, a semi-analytical method, i.e., the time domain collocation method and the calculation of transition matrix of periodic response are subsequently introduced in Sects.3.4 and 3.5.

Numerical integration algorithm
For the sake of description, rewriting the governing Eq. (5) into a compact form yields where X ¼ ðu; v; wÞ T and X s ¼ ðuðt À sÞ; vðt À sÞ; wðt À sÞÞ T represent the present states and the past states, respectively. The newly introduced symbol l ¼ lðb; Q; sÞ denotes the chosen bifurcation parameters, and Fðl; X; X s Þ denotes the right-hand functions of Eq. (5). In addition, uðtÞ is the continuous initial function of the delayed differential Eq. (5). For the dynamical system (6) with time delay and strong nonlinearity, the exact analytical solution is almost unavailable and the direct numerical integration becomes an indispensable technique to solve the complex functional differential equations. In this study, the fourth-order Runge-Kutta method with fixed step and the linear interpolation technique are combined to solve the delayed differential Eq. (6) with the initial value problem.
To begin with, the procedure of dimension reduction of the delayed governor system (6) is implemented via the discretization operation [48]. By dividing the initial interval ½Às; 0 into N s segments with ðN s þ 1Þ time nodes, the discretization step denoted by Dt can then be determined by where N s is a large positive integer. Based on the discretization operation, the continuous initial function uðtÞ can be approximated as an initial vector Y 0 : The component Let h be the integral step size and satisfy the condition h Dt; the singlestep iteration map U can then be conducted as: By selecting the proper numerical integration algorithm, all components of the new iteration vector Y h can be expressed as: where and the component XðÀs þ 0:5hÞ is obtained by using the linear interpolation, namely In this paper, the discretization step Dt is chosen as a constant value in the numerical simulations by error and trial, namely Dt % 0:01 holds for the consideration of the accuracy of numerical results. It should be pointed out that the simulation distortion will occur as the discretization step increases.

Largest Lyapunov exponent and Poincaré section
Lyapunov exponent is an effective indicator for the quantitative classification of different system responses. For an autonomous dynamical system, the system is in the state of equilibrium if the largest Lyapunov exponent (LLE) is negative. Chaotic oscillation appears in the system response if there exists a positive LLE. Otherwise, the system exhibits periodic or quasi-periodic motions, if the LLE equals zero [56]. When the periodic and quasi-periodic motions need to be further distinguished, the second largest Lyapunov exponent or the Poincaré section will be examined. The LLE has received the most attention for its high efficiency in distinguishing different dynamical responses. In this paper, the estimation of LLE is based on the method proposed in Ref. [57]. According to the estimation scheme [57], as shown in Fig. 2, the largest Lyapunov exponent is calculated by where h is the integral step size and n 1 is the total number of iterations. The symbol d 0 represents the initial distance between two nearby phase trajectories, and Y 1 . d i represents the evolutionary distance between the two trajectories at the i th iteration. In this scheme, renormalization is required in 1 as the new disturbed initial vector. In addition to the Lyapunov exponents, the Poincaré section is another widely used tool to distinguish different types of system responses. Generally, the system exhibits chaotic oscillation when a large number of disordered discrete points are found on the Poincaré section. The system demonstrates quasiperiodic oscillation when a closed curve is observed on the Poincaré section. Periodic oscillation occurs in the system response when limited, and countable points are observed on the section. However, it might be inaccurate to determine the period property of periodic response by only counting the points on the Poincaré section for the autonomous system with time delay, since the self-intersection phenomenon of phase trajectory [49] may mislead the classification. To overcome this confusion, the other methods will also be applied and described in the subsequent analysis. In this paper, the Poincaré section is selected as the section of zero angular velocity, which is given by: It should be mentioned that the point of zero angular velocity on Poincaré section is detected from the binary process by varying integral step size h. Moreover, the Poincaré section P is mainly used for the classification of different periodic responses in two-parameter diagram and the presentation of singleparameter bifurcation diagram.

GPU parallel computing
Using the above-mentioned numerical computation schemes, the single-parameter bifurcation diagram of the delayed mechanical governor system can be easily obtained. However, it is often a partial understanding of the delayed dynamics by only varying a single parameter. To explore the variation of the dynamical behavior of delayed dynamics with simultaneous changes of multiple system parameters, the GPU parallel computing technique is considered here to obtain the two-parameter diagram.
By dividing the two-parameter ðb; QÞ plane into L Â L grid points uniformly, L Â L groups of independent cases can be obtained. Thus, the global dynamics can be further determined by solving the L Â L deterministic governing equations. Benefiting from the NVIDIA company introduced CUDA-capable GPU parallel computing technique, this task can be implemented. In the present GPU parallel computing program, the total integration time n 1 h is set as 10,000, where h is the integral step size and n 1 is the number of iterations. The first 3000 integration time is neglected for the elimination of transient effect, and the subsequent 7000 integration time is applied to calculate the LLE. After that, a limited integration time is taken to calculate the points on Poincaré section and symbolized by Ni, where i is the number of intersection points between the phase trajectory and Poincaré section. By running the designed parallel computing program, the two indicators are recorded for each grid point and the category of system response can be accordingly classified in the two-parameter plane. More specifically, the equilibrium point and chaos are firstly distinguished by the value of LLE. Namely, the equilibrium point is separated from other system responses by condition LLE\ À e, and chaos is detected by LLE [ e, where e is a user-defined small positive constant with e ¼ 0:001 in this paper for the consideration of numerical error. Furthermore, the simple periodic responses with LLE j j e and i 6 are characterized. For the convenience of distinction of all the presented colors, the rest system responses, e.g., the long periodic responses with i [ 6 and the quasiperiodic responses, are not specifically distinguished and represented by the same color.

Time domain collocation method
In this study, the time domain collocation method (TDCM) [50] is chosen as a theoretical tool to identify the bifurcation types of periodic response. The principle of TDCM will be briefly described in this subsection. Since the mechanical governor system (6) is autonomous and the frequency is unknown, the rescaling operation of time in Eq. (6) is conducted firstly for the convenience of the use of TDCM. Let t 1 ¼ X t; d ¼ Xs, then the time rescaling governor system becomes: where X 0 ¼ dX=dt 1 and X d ¼ Xðt 1 À dÞ. The period of system (14) has been transformed into 2p. Approximating the solution Xðt 1 Þ of Eq. (14) by the Fourier series expansion yields: where k ; a ðwÞ k Þ T ; k ¼ 0; 1; . . .; 2N are the 3 Á ð2N þ 1Þ undetermined coefficients and N is the harmonic number in the truncated Fourier series solution. Correspondingly, the derivative and the delayed state of the approximate solution can be expressed as: Substituting Eqs. (15)-(17) into the rescaled governing Eq. (14) leads to the expression of the residual function: By dividing the period interval ½0; 2p into ðM À 1Þ segments, M equidistant points t 1j ; j ¼ 1; 2; . . .; M can be obtained. Thereafter, a series of residual functions for all the equidistant points are derived as: To obtain the above 3 Á ð2N þ 1Þþ1 unknowns, i.e., 3 Á ð2N þ 1Þ harmonic coefficients A k ; k ¼ 0; 1; . . .; 2N and frequency X, the number of algebraic Eq. (18) should be larger than the number of the unknowns. That is, the number of collocation points M should satisfy the condition of 3 Á M [ 3 Á ð2N þ 1Þ þ 1. For the consideration of the balance between the simulation accuracy and time consumption, the value of M is determined by the specific case. The detailed description of the overdetermined algebraic Eq. (18) can be found in Ref. [51]. After that, by numerically solving Eq. (18), the harmonic coefficient vector A k and the frequency X can be confirmed.

Transition matrix for stability analysis
The state transition matrix is needed to analyze the stability and bifurcation of the detected periodic response. Based on the obtained harmonic solution of system (14) and regarding it as the initial condition, the state transition matrix G i at a single integral step can be theoretically determined as Accordingly, the stability and bifurcation of the periodic response can be determined by checking the eigenvalues of G, termed as the Floquet multipliers. From Ref. [58], it is known that there are n Floquet multipliers for the n-dimensional dynamical system and one of them called the trivial Floquet multiplier which is always equal to the constant 1. The stability and bifurcation of periodic orbit can be determined by the distribution of the rest ðn À 1Þ nontrivial Floquet multipliers [58]. In the present work, the two-parameter ðb; QÞ plane is uniformly divided into 500 Â 500 grid points, by considering a balance between the computation cost and resolution of two-parameter diagram. By simulating 500 Â 500 deterministic governing differential Eq. (5), the two indicators, i.e., the LLE and the number of intersection points on Poincaré section, are calculated and recorded. After that, the motion classification of the mechanical governor system can be obtained in terms of these two indicators, and the results are shown in Fig. 3a-f. In Fig. 3a-f, ten colors are used and the corresponding colored bars are presented in Fig. 3g. Here, symbol D indicates the diverging response, symbol C denotes the chaotic motions, symbol E indicates the equilibrium state, N1 $ N6 correspond to the simple periodic responses, and NQ represents the long periodic motions with N7; N8;. . . types and quasi-periodic motions. In Ref. [15], the stability boundaries of equilibrium were derived by studying the characteristic equation. And six cases were investigated including the cases of non-delay and five discrete delays, which are the same as the six cases in Fig. 3a-f. By comparing the theoretical results of Ref. [15] with the numerical results obtained in this study, the shapes and variation tendencies of the stability boundaries of equilibrium with respect to the time delay match well in these two studies. Additionally, as reported in Ref. [15], the equilibrium would lose its stability mainly through Hopf bifurcation. The bifurcating limit cycle (indicated by symbol N1 and pink region) is also identified near the stability boundaries of the equilibrium state, as shown in Fig. 3a-f. Moreover, the fine structures lying on the left bottom of the two-parameter plane and the steady-state response originating from a certain co-dimension two bifurcation of equilibrium, which were reported in Ref. [15], are also detected and shown in Fig. 3a-f. These features indicate the effectiveness of the numerical computation strategies designed in this paper.
For the nonlinear dynamics not yet uncovered in Ref. [15], the numerical results presented in Fig. 3a-f can now provide a clear interpretation and bring new insights. For example, the N2 $ N6 types of periodic responses and chaotic motions revealed here were not reported in Ref. [15]. In the parameter region of chaotic motion in Fig. 3a-c, the existence region of such motion shrinks when the time delay varies from 0 to 2. Specifically, no chaos region is found in the twoparameter diagram in Fig. 3c with s ¼ 2 and the stable region of equilibrium marked by the gray shadow is divided into two parts in this case. To clarify this ''break'' phenomenon, the effect of time delay on the evolution of the stability boundary of the equilibrium is further investigated. It is found that the observed ''break'' phenomenon is caused by the range limitation of parameter b, rather than by any novel bifurcation of the equilibrium. To better illustrate such phenomenon, the stability boundary of equilibrium with s ¼ 1:7 is plotted in Fig. 4. As indicated in Fig. 4, the left part of the vertical black solid line is the region of interest under the consideration of nonlinear dynamics analysis and the right part is not considered in this study. Besides, although no chaotic response occurs in the system with s ¼ 2 in Fig. 3c, the existence region of chaotic motion would emerge again with the increase of time delay, as can be seen from Fig. 3d-f.
In addition to chaotic motion, different periodic responses are also revealed in the two-parameter region, where the equilibrium loses its stability. According to the distributions of periodic response, the transitions of periodic responses seem to become more and more complicated with the increase of time delay, especially for the cases of Fig. 3e-f. It should be noted that the boundaries of different colors shown in Fig. 3a-f do not indicate the actual bifurcation curves.
To determine the bifurcations shown in Fig. 3a-f, the single-parameter bifurcation analysis would be required for a better understanding of the obtained two-parameter results.

Continuation analyses for the cases of s ¼ 0 and s ¼ 1
For the case of s ¼ 0 in Fig. 3a, the single-parameter bifurcation diagrams are obtained by using the direct numerical integration (DNI) and path-following strategy. The corresponding results in terms of the chosen Poincaré section determined by Eq. (13) are marked by the blue dots and shown in Fig. 5a-b. It can be seen that there exist two types of jump phenomena in these diagrams, i.e., the jump between two branches of stable N1 type of periodic response and the jump between N3 type of periodic response and chaos. To conduct the continuation analysis and reveal the transitions of stabilities, the time domain collocation method (TDCM) and the transition matrix mentioned in Sects. 3.4 and 3.5 are applied to identify the bifurcations. Firstly, for the jump between two branches of stable N1 type of periodic response in Fig. 5a-b, the unstable solutions are tracked and marked by the red dotted lines. It is found that such jump phenomena are mainly caused by the saddlenode (SN) bifurcations of the N1 type of periodic response. Secondly, for the jump between the N3 type of periodic response and chaotic motion, the continuation analysis is also implemented to investigate the generation mechanism of the N3 type of periodic solution. For the case shown in Fig. 5a, the appearance and disappearance of the N3 type of periodic response are both induced by the SN bifurcation. However, some new coexisting attractors can be discovered and displayed in the inserted windows in Fig. 5b. These newly detected attractors might be often neglected if the continuation analysis of the N3 type of periodic response is not fully considered. Upon increasing the bifurcation parameter b of Fig. 5b, the stable N3 type of periodic response shall lose its stability through the subcritical period-doubling (PD) bifurcation on the vertical red dotted line, as shown in the inserted windows. The unstable N3 type of periodic solution marked by the red dots appears accordingly. Thereafter, the unstable N3 type of periodic solution will regain its stability by a supercritical PD bifurcation.
What should be noted is that a chaotic attractor will appear near such a supercritical PD bifurcation. And the newly stable N3 type of periodic response will finally disappear via a SN bifurcation.
For the case of s ¼ 1 of Fig. 3b, the singleparameter bifurcation diagrams and the continuation results are presented in Fig. 6a-b. From Fig. 6b, it can be seen that a typical PD bifurcation route from periodic solution to chaotic solution is observed. Therefore, the focus will be put on the jump phenomena in Fig. 6a. Different from the jump mechanism of the N1 type of periodic response revealed in Fig. 5a, the case in Fig. 6a is more complicated. Two branches Fig. 6 Bifurcation analyses for the black solid lines in Fig. 3b. a, b Single-parameter bifurcation diagrams with b ¼ 0:09 and Q ¼ 1:6, respectively. c, d The top two nontrivial Floquet multipliers along with the unstable branches 1 and 2 in Fig. 6a, respectively, where k 1 j j! k 2 j j holds of unstable periodic solutions, i.e., branches 1 and 2, are traced. To better illustrate the evolutionary process in the case, variation tendencies of the largest two Floquet multipliers with the change of parameter Q are presented in Fig. 6c-d. Here, the trivial Floquet multiplier 1 (which always exists) is not displayed in the diagrams. Apart from this trivial multiplier 1, the rest multipliers are ranked and the top two multipliers, i.e., k 1 ; k 2 with k 1 j j! k 2 j j, are displayed. The ranked multipliers corresponding to branch 1 are shown in Fig. 6c, while those of branch 2 are shown in Fig. 6d. It can be seen that a limit point exists at the intersection point between these two branches, since a nontrivial Floquet multiplier þ1 is calculated and the largest Floquet multiplier is larger than 1, i.e., k 1 [ 1; k 2 % þ1. Besides, one pair of purely imaginary eigenvalues are calculated in Fig. 6d, which indicates a Neimark-Sacker (NS) bifurcation of branch 2. As for the jump between the N3 type of periodic motion and chaos in Fig. 6a, similar scenarios in Fig. 5b are found and the newly bifurcating system responses are shown in the inserted windows.
Based on the above analyses, it is found that rich phenomena of multistability can be observed in the system, such as the coexistence of periodic and chaotic responses, and the coexistence of varieties of periodic responses. The details on these coexisting phenomena will be discussed later.

Winding phenomena and bifurcations with s ! 2
Increasing the value of time delay to s ¼ 2, the corresponding single-parameter bifurcation diagrams are plotted in Fig. 7a-b. It is easy to notice that a supercritical PD bifurcation will arise in these cases.
And the system will exhibit simple behaviors, e.g., the equilibrium state and simple periodic oscillations. When the time delay is further increased to s ¼ 3, the single-parameter bifurcation diagrams shown in Fig. 8a-b become more complicated and the transition mechanisms become hard to be identified. From Fig. 8a-b, it seems that a supercritical PD bifurcation occurs in these cases, while some extra response dots would exist along the supercritical PD bifurcation tree. To verify the evolutionary process in these cases, a set of phase diagrams are inserted in Fig. 8a-b. From the inserted diagrams, it can be seen that the winding phenomena of phase trajectories appear, i.e., the period property of periodic responses cannot be clearly determined by the number of intersection points on the Poincaré section. To address this problem, the TDCM will be applied to obtain the analytical solutions of the periodic orbits and further clarify the bifurcation types.
Taking region I (black dotted circle) in Fig. 8a as an example, N ¼ 12 is chosen in the Fourier series expansion, i.e., Eq. (15), and the number of collocation points M is set as M ¼ 34 in the case. Accordingly, the analytical solutions of the periodic responses Fig. 7 Single-parameter bifurcation diagrams along the black solid lines in Fig. 3c. a, b the case of s ¼ 2 with b ¼ 0:1 and Q ¼ 1:6, respectively corresponding to region I can be derived and marked as red circles in Fig. 8c. It can be seen that the theoretical results are in good agreements with the numerical results. The relevant top two Floquet multipliers are shown in Fig. 8d. From Fig. 8d, it is known that a PD bifurcation can occur on the red dotted line and no bifurcations are detected on the black dotted lines. To further explore the transition mechanism of just before and after the black dotted lines, the relationships between the bifurcation parameter and the harmonic components of analytical solutions, i.e., the harmonic orders and amplitudes, are subsequently studied.
To begin with, re-expressing the harmonic amplitudes in Eq. (15) as where i ¼ 1; 2; 3 corresponds to the state variables u; v; w, respectively. The symbol k represents the kth harmonic order. For the periodic responses shown in Fig. 8c, the distribution of points on Poincaré section varies in the sequence of N1 ! N2 ! N3 ! N2. The variations of all thirteen harmonic components in the analytical solutions, i.e., k ¼ 0; . . .; 12, with respect to the bifurcation parameter Q are demonstrated in Fig. 8e-f, where the harmonic components for the state u are displayed in Fig. 8e and harmonic components for the state v are shown in Fig. 8f. It is noted that only the first few harmonic components are marked by c i;k in Fig. 8e-f, for their significance in characterizing the periodic orbits. By comparing the harmonic components before and after the critical lines represented by the back dotted lines, it is found that the change in the number of intersection points on Poincaré section is caused by the varying harmonic amplitudes rather than the bifurcation. Besides, the fact that the existence of even harmonic components in the analytical solutions can induce the asymmetry of periodic orbits [59] is also witnessed in this analysis. Same as the bifurcation identification procedure for the complex periodic response in Fig. 8a, the analytical solutions and the Floquet multipliers are also calculated for the case of Fig. 8b, and a PD bifurcation is determined in this case. Furthermore, the bifurcation types of periodic orbits along the chosen bifurcation lines for the cases of s ¼ 4:0637988 and s ¼ 5 in Fig. 3e-f are clarified accordingly.
For the case of s ¼ 4:0637988, the corresponding single-parameter bifurcation diagrams are shown in Fig. 9a-b. By comparing the evolutionary process of Fig. 9a-b with Fig. 8a-b, the asymmetry feature of the periodic orbits becomes more and more complicated as time delay increases. To identify the bifurcation mechanisms in Fig. 9a-b, the Floquet multipliers are calculated and the associated bifurcation types are indicated in the two diagrams. Specifically, the identifications of bifurcations by checking the Floquet multipliers in regions II and III of Fig. 9a are presented in Fig. 9c-e, respectively. In region II, the N3 type of periodic response loses its stability by the subcritical NS bifurcation, as one pair of conjugate complex Floquet multipliers with modulus 1 are detected in Fig. 9c. Meanwhile, the coexistence of attractors can be observed in this region, as shown in the inserted phase diagram. For the bifurcations in region III, the N1 and N3 types of periodic responses will both disappear by the SN bifurcation, since the approximated Floquet multiplier þ1 can be predicted in Fig. 9d, e. In region III, the coexistence of these independent periodic solutions is observed and the coexisting attractors are shown in the inserted phase portrait of Fig. 9a. For the bifurcation diagram shown in Fig. 9b, the supercritical PD bifurcation route arises. However, it is mentioned that the evolutionary process of periodic responses shown in Fig. 9b is more complicated than the case in Fig. 8b.
For the case of s ¼ 5, the resultant single-parameter bifurcation diagrams are displayed in Fig. 10a-b. As mentioned previously, the evolutionary process is much more complicated than that in Fig. 9a-b. That is, b Fig. 8 Bifurcation analyses for the black solid lines in Fig. 3d.
a, b Single-parameter bifurcation diagrams with b ¼ 0:06 and Q ¼ 1:56, respectively. The inserted windows represent the phase diagrams, where the X-axis denotes the state variable v and the Y-axis denotes the state variable u. The red circles in the phase diagrams indicate the mapping point in Poincaré section. c Enlargement view of region I of Fig. 8a, where the blue dots are derived from the direct numerical integration and the red circles are derived from the time domain collocation method. d The top two nontrivial Floquet multipliers along with the red circles in Fig. 8c. e, f Variation of the harmonic components in the analytical solutions of the red periodic responses in Fig. 8c with respect to parameter Q. The symbol c i;k ; i ¼ 1; 2 represents the harmonic amplitude for the state variables u and v, and k denotes the kth harmonic order the relationship between the distribution of intersection points on the chosen Poincaré sections and periodic characteristic becomes more indistinguishable. To distinguish different responses and clarify bifurcation behaviors in these cases, the Floquet multipliers are also calculated and the corresponding bifurcation types are inserted in Fig. 10a-b. The PD bifurcation is the main bifurcation type in these two cases. Specifically, a supercritical NS bifurcation tree is found in Fig. 10a, which coexists with a PD bifurcation tree. Accordingly, the multistabilities, i.e., periodic response coexisting with quasi-periodic response, can be observed. Moreover, the coexistence pattern of two independent periodic responses is detected in region IV of Fig. 10a and will be discussed in the following subsection.  Fig. 9 Bifurcation analyses for the black solid lines in Fig. 3e. a, b Single-parameter bifurcation diagrams with b ¼ 0:15 and Q ¼ 1:3, respectively, c The top two nontrivial Floquet multipliers for the N3 type of periodic response in region II of Fig. 9a, d, e The top two nontrivial Floquet multipliers for the N1 and N3 type of periodic responses in region III of Fig. 9a, respectively 4.2.3 Co-dimension two bifurcation of equilibrium for s ¼ 5 A co-dimension two Hopf bifurcation point of equilibrium could be found in the neighborhood (in terms of parameters) of region IV of Fig. 10a and was also reported in Ref. [15]. Based on the results of Ref. [15], four types of responses can arise in region IV, i.e., equilibrium E 0 , two periodic solutions E 1 and E 2 , and one unstable quasi-periodic solution E 3 . The bifurcations and transitions between the four responses are schematically presented in Fig. 10c. Using the DNI and path-following methods, the steady response of the system in region IV is traced and indicated by blue dots in Fig. 10d. It can be seen that the transitions of Fig. 10d are in accordance with the transitions marked by the black solid lines in Fig. 10c. As for the jump phenomenon arising between E 1 and E 2 in Fig. 10c, it is also verified by numerical results in Fig. 10d.
To verify the transition mechanisms in terms of the periodic solutions E 1 and E 2 , the analytical solutions are determined as the red circles shown in Fig. 10d. The upper branch of analytical solutions marked by intensive red circles is derived by tracing the periodic solution E 1 . The lower branch of analytical solutions marked by sparse red circles is derived by tracing the periodic solution E 2 . For the upper branch E 1 , one pair of purely imaginary Floquet multipliers k 1;2 % À0:560 AE 0:858 i can be calculated on the line l 3 and the periodic solution E 1 becomes unstable via the subcritical NS bifurcation. Furthermore, the unstable periodic solution E 1 will disappear upon passing the line l 2 and one limit nontrivial Floquet multiplier þ1 arises. For the lower branch E 2 , the SN bifurcation occurs on the line l 1 and then E 2 disappears.

Multistable phenomena and chaos routes
This subsection discusses the qualitative or quantitative changes of the system responses in the state plane. As mentioned previously, some coexisting attractors have been observed in the mechanical governor system with and without time delay. For the nondelay cases shown in Fig. 5a-b, two main bifurcation trees are observed under different initial perturbations. One tree originates from the supercritical PD bifurcation of the N1 type of periodic response, and the other is caused by the supercritical PD bifurcation of the N3 type of periodic response. Accordingly, rich coexistence phenomena can be induced in the governor system without time delay, such as the coexistence of different periodic attractors, the coexistence of different periodic attractors and chaos. In addition to the above-mentioned coexistence patterns, the coexistence of equilibrium and periodic response can also be observed in this case.
To study the evolutionary process of the basins of attraction, the coexistence of the N3 type of periodic response and chaos arising in Fig. 5a is taken as an illustrative example. Four groups of parameters are selected to compute the basins of attraction and the results are presented in Fig. 11a-d, accordingly. It is mentioned that the basin of attraction is also obtained via the GPU parallel computing technique. The state space is considered as follows: And P 1 is divided into 300 Â 300 grid points in the parallel computing program for the consideration of balance between computational accuracy and time consumption. In Fig. 11a-d, the red domain indicates the basin of attraction of the N3 type of periodic response, while the dark green domain stands for the basin of attraction of chaos. And the rest rufous domain corresponds to the divergent response. By analyzing the basins of attraction shown in Fig. 11a-d, it can be seen that the phenomenon of self-similarity structure arises. The basin of attraction of the N3 type of periodic response gradually shrinks and the basin of attraction of chaos will expand when the bifurcation parameter Q increases. Eventually, the basin of attraction of the N3 type of periodic motion disappears through a SN bifurcation, as indicated in Fig. 5a. For the mechanical governor system with time delay, some coexistence phenomena are newly found. When the time delay is small as the case shown in Fig. 6a, a new coexistence of stable N1 type and N3 type of periodic responses emerges. As the time delay increases, new dynamic phenomena are observed, e.g., the coexistence of two branches of NS bifurcation tree in region II of Fig. 9a and the co-dimension two bifurcation of equilibrium in Fig. 10a. Accordingly, two coexistence patterns, i.e., the coexistence of the N1 type and N3 type of periodic attractors, the coexistence of periodic and quasi-periodic attractors, occur in the cases of Fig. 9a and Fig. 10a, respectively. To illustrate the coexistence of the N1 type and N3 type of periodic attractors for the case of Fig. 9a, the basin of attraction is considered to be computed in the following state space: . . .; N s . In this case, the state plane ðu 0 ; w 0 Þ is displayed and also divided into 300 Â 300 grid points. Moreover, four groups of parameters are chosen to study the evolution of the basin of attraction.
The corresponding results are displayed in Fig. 12, where the red domain is the basin of attraction of the N3 type of periodic response, and the green domain represents the basin of attraction of the N1 type of periodic response. The rest rufous domain corresponds to the diverging motion under the chosen initial condition. Same as the case in Fig. 11, the self-similar basin structure is also observed in Fig. 12. Due to the range limitation of the selected state plane, the selfsimilar basin structures located on the right-hand sides of Fig. 12a-d are partially visible and would not be further studied in this paper. It is easy to find that the evolution tendency of the basin of attraction of the two different periodic solutions is identical with the bifurcation diagram in Fig. 9a.
In the chaotic transitions exhibited in the above results, three main routes leading to chaos can be found in the system response. The first route is that the chaotic motion is generated from a PD bifurcation and then disappears by the way of boundary crisis [60], e.g., Fig. 5a, b and Fig. 6a. In these cases, the PD bifurcation is involved with the supercritical type. The second route is that the chaotic motion is born and destroyed both via a PD bifurcation, e.g., Figs. 5b, 6b and 8a. The third route is that the chaotic motion arises from a NS bifurcation, such as Fig. 9a. In Fig. 9a, the subcritical NS bifurcation of the N1 type of periodic response is involved and such a route has ever been reported in Ref. [60].

Conclusion
The motion patterns in terms of ðb; QÞ two-parameter domains of a mechanical governor system have been studied and analyzed with the aid of Lyapunov exponents and the Poincaré section. The correlations between reaction delay s and governor system dynamics were discussed. It was found that many rich nonlinear dynamic behaviors could occur in this system, including simple periodic motion, complex periodic motion, quasi-periodic motion, and chaos. As time delay s varies from zero to a small value, the parameter domains corresponding to the set of complex periodic motions, quasi-periodic motions and chaos would shrink gradually and the nonlinear behaviors of the system are mainly characterized by the simple periodic motion, such as motion patterns shown in Fig. 3c. However, as the time delay s is further increased, rich nonlinear responses arise again and the transition characteristics become more and more complicated. Indeed, such complicated transitions are caused by the winding characteristics of certain periodic motions rather than bifurcations. In addition to the transition of adjacent motions, some multistability behaviors were also discussed and the basin analyses of attractors were presented in the state plane with certain cases. The typical self-similarity characteristic was observed in the basin structure of attraction.
For these rotating machines, e.g., internal combustion engines, fueled turbines, modern striking clocks, the occurrence of unstable rotating speed can greatly affect their work efficiency, operation reliability, and even service life. Therefore, reducing or eliminating undesired effects on mechanical devices is a vital goal of dynamic research. Actually, an exploration of motion patterns and bifurcation laws of the system in the multi-parameter plane can reveal the correlations between occurrence domains of motion patterns and system parameters. In this way, desirable motion can be achieved by optimizing some parameters in the matching range, according to the results of presented bifurcation analyses. In addition, the system behavior would remain uncertain due to the appearance of multistability. To ensure the better operation performance and reliability of the mechanical system, the multi-parameter, multi-state plane, and the time delay effect should be fully considered to achieve the desirable state with the optional ways.
In summary, the two-parameter dynamics of the delayed governor system have been numerically investigated to uncover the new dynamic behavior which was not explored in the previous work. It is anticipated that these research results would provide a meaningful reference to the optimal design and motion control of the mechanical systems and the exploration of local and global dynamics of the delayed systems.