Analysis of the non-periodic oscillations of a self-excited friction-damped system with closely spaced modes

It is widely known that dry friction damping can bound the self-excited vibrations induced by negative damping. The vibrations typically take the form of (periodic) limit cycle oscillations. However, when the intensity of the self-excitation reaches a condition of maximum friction damping, the limit cycle loses stability via a fold bifurcation. The behavior may become even more complicated in the presence of any internal resonance conditions. In this work, we consider a two-degree-of-freedom system with an elastic dry friction element (Jenkins element) having closely spaced natural frequencies. The symmetric in-phase motion is subjected to self-excitation by negative (viscous) damping, while the symmetric out-of-phase motion is positively damped. In a previous work, we showed that the limit cycle loses stability via a secondary Hopf bifurcation, giving rise to quasi-periodic oscillations. A further increase in the self-excitation intensity may lead to chaos and finally divergence, long before reaching the fold bifurcation point of the limit cycle. In this work, we use the method of complexification-averaging to obtain the slow flow in the neighborhood of the limit cycle. This way, we show that chaos is reached via a cascade of period-doubling bifurcations on invariant tori. Using perturbation calculus, we establish analytical conditions for the emergence of the secondary Hopf bifurcation and approximate analytically its location. In particular, we show that non-periodic oscillations are the typical case for prominent nonlinearity, mild coupling (controlling the proximity of the modes), and sufficiently light damping. The range of validity of the analytical results presented herein is thoroughly assessed numerically. To the authors’ knowledge, this is the first work that shows how the challenging Jenkins element can be treated formally within a consistent perturbation approach in order to derive closed-form analytical results for limit cycles and their bifurcations.

bifurcation and approximate analytically its location. In particular, we show that non-periodic oscillations are the typical case for prominent nonlinearity, mild coupling (controlling the proximity of the modes), and sufficiently light damping. The range of validity of the analytical results presented herein is thoroughly assessed numerically. To the authors' knowledge, this is the first work that shows how the challenging Jenkins element can be treated formally within a consistent perturbation approach in order to derive closed-form analytical results for limit cycles and their bifurcations.

Introduction
Mechanical systems usually consist of several substructures assembled by joints in which nonlinear contact interactions inevitably occur [1]. The dissipative dry frictional interactions often provide the dominant contribution to the overall energy dissipation in mechanical systems [2]. Besides vibration mitigation, friction damping has acquired high technical relevance because of its stabilizing effect on otherwise unstable self-excited vibrations. These are the consequence of destabilizing negative damping contributions which can be encountered in various applications: flutter vibrations of airfoils [3], chatter in machining processes [4,5] like in rolling mills [6], in control tasks [7][8][9], in robotics [10], and in brake squeal [11]. The driving example for the current work is cascade flutter vibrations of turbomachinery blades [12,13]. Additional frictional interfaces are here often deliberately introduced to mitigate vibrations of the blades [2]. In this application, the self-excitation mechanism is caused by the unstable aerodynamical interference of blades within a blade row. In a recent study, we numerically analyzed the stability of self-excited frictiondamped vibration states in the presence of internal resonance conditions based on a minimal model [14]. The main contribution of the present work is the analytical derivation of existence conditions of quasi-periodic oscillations for closely spaced modes. In the following paragraphs, we will briefly recap the most important results related to this work, before discussing the purpose and outline of this article.

Preliminaries
Sinha and Griffin [15] were the first to show that friction damping can be used to bound otherwise diverging self-excited vibrations caused by negative (viscous) damping. In case of negative damping of a single mode and well-separated eigenfrequencies, the amplitudedependent friction damping leads to periodic limit cycle oscillations (LCO). If the negative damping exceeds the point of maximum friction damping, a fold bifurcation occurs such that the LCO becomes unstable giving rise to divergent behavior. The dynamics may become even more complicated when several modes are subjected to negative damping. Based on a simplified model of an aerodynamically unstable rotor stage, the stability limit of LCOs was determined by generalized energy relations based on an averaging technique [16]. However, the bifurcation type associated with the stability limit was not addressed in these works. It is well known that friction nonlinearity can trigger modal interactions in the presence of internal resonance conditions (e.g., 1:1 or 1:3 resonances). In [14], we studied numerically a friction-damped two-degree-of-freedom (2DOF) system subjected to negative viscous damping. The configuration of the investigated 2DOF system is motivated by the analysis of rotationally symmetric structures, e.g., bladed disks, where so-called traveling wave modes occur. The 2DOF system's modes are representatives of a forward and backward traveling mode. Neglecting frictional effects, the system remains symmetric, with both modes closely spaced based on the assumption of light aeroelastic coupling. Under the influence of friction damping on one of the oscillators, symmetric behavior will only occur in the case of sticking contact. In the case of sliding friction, energy dissipation affects both modes, but also leads to asymmetrical behavior. This asymmetry causes strong modal coupling, so that a strong energy exchange between both modes is expected close to internal resonance conditions. In Fig. 1, the characteristic behavior for the case of closely spaced modes (near 1:1 internal resonance) is illustrated in a schematic bifurcation diagram with the viscous damping being considered as control parameter. The trivial case of positive damping leads to exponentially decaying oscillations that approach the static equilibrium position (A). Going from positive to negative damping induces a Hopf bifurcation, causing dynamic instability of the static equilibrium. The vibrations grow until dissipative sliding friction compensates the self-excitation and an LCO is formed from the dynamical balance (B). A secondary Hopf bifurcation causes the loss of periodicity of the LCO beyond a critical level of negative damping, as indicated by the splitting of maximum and minimum steady-state vibration amplitudes (C). Then, the oscillations take the form of limit torus oscillations (LTOs) at first, much in line with the well-known case of coupled van-der-Pol oscillators. More severe negative damping may introduce chaos. Finally, the non-periodic attractor collapses yielding divergent behavior well before the fold bifurcation of the LCO (D) is reached. These findings imply that the common practice to compute only the LCOs (using, e.g., Harmonic Balance and numerical contin-uation), without properly analyzing their stability, may lead to incomplete modeling of the nonlinear dynamics, and, hence, to potentially dangerous design decisions in practical applications.

Purpose of the present work
In this work, we seek to better understand the emergence of the non-periodic oscillations encountered for the 2DOF system discussed in [14] for the case of closely spaced modes. In particular, we would like to establish conditions for the secondary Hopf bifurcation and analyze how the quasi-periodic oscillations further bifurcate into chaos. In the analytical development, the key challenge is the treatment of the Jenkins element 1 . To the authors' knowledge, this is the first work that shows how a consistent perturbation approach can be formally applied to analyze this nonlinearity in order to derive closed-form analytical expressions for the emergence of LCOs and their bifurcations.

Outline of the article
The article is structured as follows: In Sect. 2 we describe the investigated 2DOF system. In Sect. 3, we use averaging to derive an approximation for the slow flow dynamics in the neighborhood of the LCO. In Sect. 3.3, we analyze the transition from limit torus oscillations to chaotic oscillations. In Sect. 4, we focus on the secondary Hopf bifurcation and use perturbation calculus to establish the conditions leading to limit torus oscillations. We end this work with some concluding remarks, summarizing the main contributions and providing some related open problems that can be further investigated.

Model description
We investigate a 2DOF oscillator (see Fig. 2a), consisting of two unity masses coupled via linear springs in parallel to viscous damper elements. The nonlinear dry friction force acting on the first oscillator is modeled by where x 1 and x 2 indicate the generalized coordinates and˙ indicates derivative with respect to time t. The nonlinear friction force f t [x 1 ] is defined by the hysteretic differential law where the dependency on time history is indicated by the operator [ · ]. An advantage of this formulation is the straightforward implementability within numerical time-step integration methods by substituting the dif-ferential by a finite time-step. Figure 2c illustrates the corresponding hysteresis cycle under prescribed periodic motion of the coordinate x 1 , including the initial loading curve. The hysteresis cycle divides into sticking and sliding phase of the friction element, transitioning as soon as the spring force exceeds the friction limit force of 1. In the sliding phase, the magnitude of the nonlinear force remains constant. A similar dynamical system was obtained by nondimensionalization in [14]. The non-dimensionalization takes into account time normalization with respect to the natural frequency of the in-phase mode in sticking condition, the normalization of the nonlinear force by the friction limit force, and the normalization of the generalized coordinates by the fraction of nonlinear spring stiffness and friction limit force. This reduces the set of parameters to the coupling spring stiffness κ, the modal damping ξ 1,2 , and the friction intensity γ acting on the first oscillator. Note that the Jenkins element itself remains independent of γ , but it appears as a factor in front of the nonlinear force, preserving the symmetry of the system in sticking condition. Thus, the linear modes of vibration correspond to symmetric in-phase and out-of-phase motions. The corresponding in-phase natural frequency results to ω st = 1 (due to normalization). The natural frequency of the outof-phase mode is given by √ 1 + 2κ [14]. Thus, for κ 1 (i.e., for weak coupling between oscillators) the linear natural frequencies are closely spaced. The (weak) damper coefficients are designed such that ξ 1 is associated with the symmetric in-phase and ξ 2 with the symmetric out-of-phase motion. Throughout this paper, we set ξ 1 < 0 and ξ 2 > 0. For sufficiently large vibrations, the force in the friction element saturates, which leads to nonlinear stick-slip behavior and generally asymmetrical behavior.

Analysis of the slow dynamics in the neighborhood of the limit cycle
In the following, we gain deeper understanding of the dynamics by analyzing the slow dynamics of the system near the LCO. To this end, we use the complexificationaveraging (CX-A) method [17]. Of particular interest is the bifurcation point of the first LCO stability loss, i.e., the secondary Hopf bifurcation point, giving rise to LTOs. Furthermore, the transition process from the LTOs to chaos will be analyzed.

Extraction of the slow dynamics using CX-A
In the following, we briefly explain the main ideas behind the CX-A technique, following [17], while generally there exist several books discussing averaging techniques. CX-A is based on a coordinate transform from a fixed to a rotating coordinate system, which can be expressed by a new complex-valued variable , yielding with the oscillation frequency ω, the complex conjugatē and i = √ −1. This transformation reduces the order of the differential equation system from 2nd to 1st, while retaining the number of equations. At this point, dynamics dominated by a single harmonic are assumed, allowing to introduce the complex-valued ansatz Note, that the integral over frequency ωdt degenerates to ωt only in case of time invariance. Based on this, the dynamics with respect to x j is assumed to be decomposable into a fast oscillation with the instantaneous (fast) oscillation frequency ω and a slow modulation of amplitude a j and phase θ j . Under the premise of well-separable time scales, averaging over the fast dynamics and splitting into real and imaginary parts, leads to a set of explicit first-order differential equations with respect to the slow dynamics. Applied to Eqs. (1) and (2), the following set of equations is obtaineḋ in terms of the unknowns a 1 , a 2 , θ 1 , θ 2 and their associated time derivatives. Without loss of generality, the amplitudes are defined as positive a 1 > 0, a 2 > 0. The terms k eq (a 1 ) and d eq (a 1 ) are the equivalent stiffness and damping of the Jenkins element as defined later. Without loss of generality, one can reduce the problem to an autonomous set of 3 equations governing the slow amplitude modulations a 1 and a 2 and the slow phase difference (lag) which is defined as θ = θ 1 − θ 2 . This leads to what we will refer to as slow flow equations, Averaging the nonlinear force of the Jenkins element, f t [x 1 ], yields the expressions for the amplitudedependent equivalent stiffness k eq and damping d eq [18] as follows: Fig. 3 Behavior of k eq and d eq with respect to τ * ; point of divergence ξ 1,div determined by numerical simulation of Eqs. (1) and The parameter τ * can be interpreted as the stick-slip transition time. In the sticking case it holds that, a 1 ≤ 1. In this case, the stick-slip transition does obviously not occur, which corresponds to τ * = π (due to symmetry of hysteresis cycle). Consequently, the Jenkins element correctly degenerates to a conservative spring with k eq = 1 and d eq = 0. For increasing vibration levels beyond the sticking limit, i.e., a 1 > 1, the stick-slip transition occurs earlier along one period of vibration, i.e., τ * decreases. The maximum equivalent damping is reached for τ * = π/2, where k eq = 0.5. As the amplitude goes to infinity, the friction effect becomes negligible such that both k eq and d eq asymptotically approach zero. The nonlinear dependence of the equivalent damping and stiffness with respect to the parameter τ * is depicted in Fig. 3. The gray shaded area, limited by τ * (a 1 (ξ 1,div )), indicates divergent curve ranges, i.e., τ * related to unbounded amplitude growth. The indicated limit is based on numerical results for the given parameter values. To simulate the slow flow given by Eqs. (10)-(12), the oscillation frequency ω needs to be specified. We are interested in the dynamics near the LCO. Indeed, to compute the LCO with frequency ω, we impose conditions for periodicity, namely,ȧ 1 =ȧ 2 =θ 1 =θ 2 = 0. With this, Eqs. (6)-(9) degenerate to an algebraic system of four equations for the four unknowns a 1 , a 2 , θ and ω, which are solved numerically. The resulting frequency of the LCO is illustrated in Fig. 4, where it is depicted as function of the corresponding mechanical energy of the system averaged over one period of vibration when it oscillates on the LCO, obtained by varying the negative damping ξ 1 < 0. The averaged mechanical energy, i.e., the averaged sum of potential and kinetic energy, is formally given by where τ = ωt. This frequency-energy relation holds strictly only at the LCO, i.e., for the time-periodic solution, which represents a special stationary solution of the mechanical system with negative damping. Perturbing the initial conditions corresponding to the LCO, but staying sufficiently close in the neighborhood of the LCO, we postulate that this frequency-energy relation still provides a reasonable approximation denoted symbolically in the form, Of course, the above described numerical solution method only generates a discrete sequence of points.
To simulate the slow flow equations (10)-(12), however, a continuous expression (i.e., a function) of the oscillation frequency in terms of the states a 1 , a 2 , θ is needed, which we obtained using interpolation. To this end, we applied a piecewise cubic Hermite interpolation between about 400 discrete solution points in the interval E mech ∈ (0, 500). Note, that this interval contains also energy levels beyond the fold bifurcation (i.e., points along the upper half of the LCO curve), cf.  1) and (2)]. We note at this point that the ansatz used for the previous complexification-averaging analysis is based on the assumption that there exists a single fast frequency ω, and that there is a slow-fast separation of the non-stationary dynamics, i.e., that there is no mixing of scales. Whenever either of these conditions is not satisfied we expect that our analytical predictions will diverge compared to the exact solutions derived by direct numerical simulations of the initial system of governing equations. In Fig. 5, we compare the bifurcation diagrams, with respect to the negative damping  (12)) against the initial system (Eqs. (1) and (2)) with respect to the amplitude-damping diagram: The discrete points are obtained by direct numerical integration initialized on the LCO considering a small perturbation; ξ 1 < 0, parameter ξ 1 . The vibration level is depicted in terms of the maximum, a 1,max , and minimum, a 1,min , vibration amplitudes at steady state; we note that the derived steady state is not necessarily periodic, nor even stationary. In spite of the complexity of the dynamics of the system, both the secondary Hopf bifurcation point, ξ 1,Hopf , and the point of divergence, ξ 1,div , are captured by the analysis with remarkable accuracy. The minor deviations are attributed to the aforementioned analytical simplifications: Truncation to a single harmonic, separation of time scales, and presumed invariance of frequency-energy dependence. Figure 6 depicts Poincaré maps associated with a Poincaré section S (associated with the fast time scale).
To this end, given that the initial system [Eqs. (1) and (2)] is four-dimensional, we introduce the following with respect to the initial, fast variables x = [x 1 , x 2 ,ẋ 1 , x 2 ] T . Studying the intersections of the orbits of the dynamical system with the Poincaré section S we obtain a three-dimensional Poincaré map. For simplicity, we depicted just two-dimensional projections of this three-dimensional map, namely theẋ 1 -x 1 -plane. It can be seen that the LTO transitions toward chaotic oscillations for increasing magnitude of the negative damping. Analytical predictions obtained from the slow flow approximation are depicted in Fig. 7. Here, the evolution of the variables x on the fast time scale was recovered using Eqs. (4) and (5). Again, the slow flow approximation is in good agreement with the exact numerical simulation, however, a slight shift of the dynamics appears to be introduced below ξ 1 < −0.1. We conclude that the proposed slow flow approximation, including the assumptions on the frequencyenergy dependence, is sufficiently accurate to further characterize the dynamics on the slow time scale.
with respect to the slow variables a = [a 1 , a 2 ,ȧ 1 ,ȧ 2 ] T . Here, denotes mean-free variables. Note that the amplitudes a 1 and a 2 typically do not have zero crossings. Hence, the steady-state mean value is subtracted so that one can use the common formulation of the Poincaré section in Eq. (18) (involving zero crossings).
The slow time scale provides a clearer picture of the dynamics: The (quasi-periodic) limit torus oscillations are seen as closed curves in the phase projections within Fig. 8. The Poincaré map consists of a single point for ξ 1 = −0.055. Apparently, a cascade of period-doubling bifurcations takes place until the chaotic regime is reached. This behavior is consistent with the results shown in Fig. 5 for the original system, but occurs slightly shifted for ξ 1 < −0.1 2 . In Table 1, we listed estimations for the first successive perioddoubling bifurcation points and their ratios δ i for the slow flow system. Apparently, the period doublings do not approach the Feigenbaum ratio (δ F ≈ 4.6692) [19], at least for the slow flow system.
In [14], the chaotic nature of the original system's dynamics past the secondary Hopf bifurcation was confirmed by numerically estimated Lyapunov exponents following [19]. This method is particularly suitable for smooth systems, while the results for non-smooth systems (e.g., due to frictional forces) should generally be taken with caution. The estimation of Lyapunov exponents for this type of systems is still under research, a summary of different approaches can be found for example in [20].

Analytical investigation of the secondary Hopf bifurcation
We now focus on the transition from limit cycle oscillations to limit torus oscillations. We seek to establish analytically the conditions for the emergence of non-periodic oscillations in the mechanical system of Fig. 2, and to obtain a closed-form approximation for Table 1 Estimation of the first period-doubling bifurcations i with the associated self-excitation level ξ 1,i and their ratio δ i for the slow flow approximation to zero. It turned out to be useful to transform the slow flow by introducing a new set of variables: Instead of a 1 , a 2 , we introduce the new variables ρ = a 2 /a 1 and τ * (stick-slip transition time) which is uniquely related to a 1 as defined in Eq. (14). With this variable transformation, the governing algebraic equations determining the fixed points of the slow flow can be cast into the form: We use perturbation calculus to approximate the fixed points of this set of transformed algebraic equa-tions. To this end, we introduce the small parameter ε which will act as the formal perturbation parameter of our analysis. To introduce the small parameter we assume that the coupling and damping terms are small, scaled according to the following relations: Recall that the condition of weak coupling is needed to have closely spaced natural frequencies in the mechanical system of Fig. 2: The ratio of the natural frequencies in the linear case when only sticking contact occurs equals √ 1 + 2κ. A typical value for κ is 0.1 leading to about 10% difference in the natural frequencies. Moreover, the secondary Hopf bifurcation occurs already for very light damping, with typical values for ξ 1 of − 0.035 and for ξ 2 of 0.011. Thus, it holds that ξ 1,2 κ 1, as defined in Eq. (23). It should be emphasized that nonlinearity is not assumed as weak, i.e., γ is considered to be a system parameter of order O(ε 0 ). Yet we impose the condition γ ≤ 1 which is necessary to have positive stiffness values for all springs in the model (Fig. 2). We seek an approximation to the fixed point in terms ofτ * ,ρ, θ andω, whereˆ denotes fixed point quantities. In accordance with regular perturbation theory (see, e.g., [17]), these quantities are expanded in regular power series in ε, of the form As will be shown later, reasonable results are achieved already for a second-order approximation (up to and including the ε 2 term). The power series and Eq. (23) are substituted into Eqs. (19)-(22). The equations contain several nonlinear terms, namely, sin( θ ), cos( θ ), ωd eq (τ * ), k eq (τ * ), and 1/ρ. In accordance with perturbation calculus, each nonlinear term is expanded into a Taylor series of sufficient order, and the power series of the unknowns (Eq. (24)) is substi-tuted. The unknown coefficients of the power series are then determined successively by matching (balancing) the coefficients at each order separately, starting with the leading-order O(ε 0 ) terms. With this, the following second-order approximation of the fixed point of the slow flow equations is obtained: As expected, non-trivial fixed points only exist for ξ 1 ≤ 0.

Calculation of the secondary Hopf bifurcation point
Clearly, the limit torus on the fast time scale corresponds to a periodic orbit on the slow time scale, and the bifurcation of the limit cycle to the limit torus corresponds to a secondary Hopf bifurcation. Consequently, on the slow time scale, this point is associated with (primary) Hopf bifurcation of the fixed point of the slow flow equations given by Eqs. (25)-(28). Hence, to locate this point and establish conditions for its existence, it suffices to analyze the asymptotic stability of the approximated fixed point of the slow flow and prove that the conditions for Hopf bifurcation hold at the point of instability of the fixed point. This is done in the conventional way, based on the eigenvalues of the Jacobian with respect to the slow flow equations (10)- (12). Accordingly, the characteristic equation governing the eigenvalues, λ, can be explicitly computed by imposing the following condition, where the following abbreviations have been introduced, d 1 = 1 2 ξ 1 +ξ 2 , d 2 = 1 2 ξ 1 −ξ 2 , s = sin θ and c = cos θ , together with the additional notations, The solution of the characteristic equation yields the following second-order analytical approximation for the eigenvalues: Recalling that ξ 1 < 0, one can infer that λ (1) remains real-valued, whereas λ (2,3) appear as a complex-valued pair. Thus, the sought bifurcation occurs when λ (2,3) crosses the imaginary axis. This is the necessary condition for the realization of the Hopf bifurcation; it can also be shown (not reported here) that the additional sufficient conditions are also met Apparently, besides the expected transition from limit cycle to limit torus, another transition back to a stable cycle is possible. This will be further discussed in the next subsection. It is interesting to note that the bifurcation valueξ start/end 1,Hopf is a function of the outof-phase damping ξ 2 and the coupling-to-nonlinearityratio κ 2 /γ . This relation is illustrated in Fig. 9a. Solutions forξ start/end 1,Hopf only exist for ν ≥ 3, i.e., if κ 2 γ ≥ 3ξ 2 π . This means that the positive damping of the (outof-phase) higher-frequency mode must be sufficiently small (for given γ and κ) for this bifurcation to occur. We now compare the previous analytical results to the results of exact numerical simulations. The location of the secondary Hopf bifurcation, in terms of the critical damping value ξ 1,Hopf , is illustrated in Fig. 9. The corresponding analytical result is given in Eq. (31). As throughout this work, the numerical reference was obtained by simulation of the initial system [Eqs. (1) and (2)].
Overall the analytical result is in good agreement with the numerical reference. Deviations occur at the ranges where the assumptions underlying the perturbation approach are violated. More specifically, the assumption concerning the hierarchy of the orders of magnitude, ξ 1,2 κ γ (cf. Eq. (23)) must hold to ensure reasonable accuracy. For instance, the analytical results are accurate for ξ 2 < 0.02 (Fig. 9b) and ξ 2 < κ < 0.2 (Fig. 9c). In contrast, if ξ 1 or ξ 2 approach the order of κ, the accuracy is poor (as expected). Similarly, if κ approaches the order of γ (order ε 0 ), the analytical result is no longer valid. As noted earlier, the analytical approximation yields two critical values of ξ 1 , as expressed in Eq. (31). The one with smaller magnitude,ξ start 1,Hopf , corresponds to the point where the limit cycle becomes unstable giving rise to a stable limit torus. The one with larger magnitude,ξ end 1,Hopf , refers to another secondary Hopf bifurcation where the limit cycle regains stability. As seen from these results, for a wide range of parameters,ξ end 1,Hopf lies beyond the point of divergence, ξ 1,div , and, therefore, it is never reached. The point of divergence and the corresponding gray shaded area indicate parameter ranges leading to unbounded amplitude growth based on numerical simulation of Eqs. (1) and (2). For other parameter combinations, however, the secondary Hopf bifurcation can be reached, e.g., for 0.6 < γ < 0.8 in Fig. 9d, or for κ ≈ 0.074 and κ ≈ 0.21 in Fig. 9c. The resulting bifurcation diagram is illustrated in Fig. 10 for three different parameter sets. Again, reasonable agreement between the analytical results and numerical solutions is achieved with respect to the bifurcation point.

Conclusions
We considered a system of two weakly coupled oscillators, with an elastic dry friction element attached to one of the oscillators, subjected to negative viscous damping of the symmetric in-phase and positive damping of the out-of-phase motion. When increasing the magnitude of the negative damping, the well-understood limit cycle oscillations lose stability and bifurcate into limit torus oscillations. For some parameter ranges, an additional secondary Hopf bifurcation may occur beyond which the limit cycle regains stability. In other cases, the quasi-periodic oscillations further bifurcate into chaotic dynamics. The analysis of the slow flow in the neighborhood of the limit cycle revealed that chaos is reached via a cascade of period-doubling bifurcations. Assuming weak coupling, of O(ε), (leading to closely spaced natural frequencies) and even weaker damping, of O(ε 2 ), we derived an analytical approximation of the limit cycle using perturbation calculus. A stability analysis provided a closed-form expression for the point of the secondary Hopf bifurcation and conditions for its existence. In particular, it was established that non-periodic oscillations occur for sufficiently weak coupling and sufficiently strong nonlinearity. The quality of the analytical results was thoroughly validated by comparing them against exact numerical simulations of the original equations of motion of the initial dynamical system.
The key challenge for deriving the closed-form approximations was the treatment of the dry friction nonlinearity (Jenkins element). To tackle this, a change of variables was carried out, involving the normalized stick-slip transition time. As the Jenkins element is commonly used to model nonlinear hysteresis (including tribology and material behavior), we believe that the proposed approach will be instrumental to analyze systems far beyond the one considered in this work. Here, we anticipate that our results will trigger more efforts within our scientific community to derive analytical approximations in order to gain deep understanding of parameter dependencies, in addition to systematic numerical simulation. This will help in predictive design of oscillating systems with dry friction undergoing instabilities, LCOs, LTOs, and chaotic motions.
Acknowledgements The author Alexander F. Vakakis would like to acknowledge the support of a von Humboldt Award which made possible his visit at the University of Stuttgart.
Funding Open Access funding enabled and organized by Projekt DEAL.

Availability of data and material
The datasets generated during the current study are available from the corresponding author on reasonable request.

Conflict of interest
The authors declare that they have no conflict of interest.

Code availability
The codes generated during the current study are available from the corresponding author on reasonable request.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/.