Fundamental competition of smooth and non-smooth bifurcations and their ghosts in vibro-impact pairs

A combined analysis of smooth and non-smooth bifurcations captures the interplay of different qualitative transitions in a canonical model of an impact pair, a forced capsule in which a ball moves freely between impacts on either end of the capsule. The analysis, generic for the impact pair context, is also relevant for applications. It is applied to a model of an inclined vibro-impact energy harvester device, where the energy is generated via impacts of the ball with a dielectric polymer on the capsule ends. While sequences of bifurcations have been studied extensively in single- degree-of-freedom impacting models, there are limited results for two-degree-of-freedom impacting systems such as the impact pair. Using an analytical characterization of impacting solutions and their stability based on the maps between impacts, we obtain sequences of period doubling and fold bifurcations together with grazing bifurcations, a particular focus here. Grazing occurs when a sequence of impacts on either end of the capsule are augmented by a zero-velocity impact, a transition that is fundamentally different from the smooth bifurcations that are instead characterized by eigenvalues of the local behavior. The combined analyses allow identification of bifurcations also on unstable or unphysical solutions branches, which we term ghost bifurcations. While these ghost bifurcations are not observed experimentally or via simple numerical integration of the model, nevertheless they can influence the birth or death of complex behaviors and additional grazing transitions, as confirmed by comparisons with the numerical results. The competition between the different bifurcations and their ghosts influences the parameter ranges for favorable energy output; thus, the analyses of bifurcation sequences yield important design information.


Introduction
Vibro-impact (VI) systems present a remarkable class of nonlinear systems where impacts drive the interaction between the system components. A classical VI system is comprised of a forced mass and motionless rigid barrier(s) or a pair of moving impacting masses, either of which can be forced. Some examples of VI systems are shown in Fig. 1. The simplest VI models include both that of a ball bouncing on a harmoni-cally moving surface, studied in [1,2] and a pendulum impacting a barrier [3]. A complete dynamical model includes equations governing behavior both between impacts and at the impacts. Thus the dynamics fall into the general class of non-smooth dynamics [4][5][6]. Even when a linear system governs the dynamics between impacts, the VI system exhibits nonlinear behavior that necessarily follows from these non-smooth interactions.
Modeling the impact dynamics can take on different levels of complexity. Using a restitution coefficient r , defined as the ratio between the mass velocity immediately after and before impact, implies that it happens instantaneously or on a negligibly short time scale as compared to the overall system. Typically 0 ≤ r ≤ 1, depending on the material, shape and contact surfaces of the impacting bodies. Systems with instantaneous impacts are strongly nonlinear, given the discontinuous velocities of the impacting bodies before and after any non-elastic impact, since the velocity is reduced for r < 1 and changes direction for r > 0.
The specific type of non-smooth nonlinearity, due to the instantaneous impact approximation in the VI systems, can facilitate analysis in some cases. Indeed, when the energy losses due to impacts (r < 1) dominate other damping mechanisms, the latter are neglected, yielding a conservative system between impacts that can be studied analytically. Complementing the system response with the impact condition(s) for its velocity allows one to derive an explicit analytical solution for its periodic motion, as in [7], and to study the dynamical dependence on system parameters, including the classical nonlinear analyses of smooth bifurcations observed in VI systems. Such analytical results cannot capture the non-smooth behavior explicitly, so that complementary numerical simulations or semianalytical methods are generally required. One type of non-smooth bifurcation, first reported in [8], is a grazing bifurcation characterized by a zero-velocity impact. Later, in [9,10] gave a global bifurcation analysis, indicating the importance of the grazing phenomenon on the global dynamics of VI systems. In the context of a spring-mass system colliding against a motionless barrier [11], A. Nordmark introduced the term grazing impact associated with the flow tangential to the impacting surface. It is treated as a bifurcation, separating non-impacting and vibro-impacting motions. The near-grazing two-dimensional map was also derived in [11], demonstrating the continuity of the bifurcating Fig. 1 Examples for different VI systems. Sketch of a a twodegree-of-freedom (TDOF) impact pair subject to base excitation y 2 and a force y 1 applied to the mass [12], b a single-degreeof-freedom (SDOF) system [13], here k b indicates the spring constant for the barrier, so that as k b → ∞, we recover instantaneous impact. c TDOF impact oscillators with motion limiting constraints for both masses [14], d two oscillators coupled only by impact [15] branch with the square root singularity of the corresponding Jacobian at the grazing point.
These earlier analyses of grazing and related maps generated deep interest in obtaining a better understanding of non-smooth dynamics in related SDOFs.
In [16] the authors compared bifurcations in instantaneous and soft impact models, hypothesizing the grazing as a limiting case of either a fold (FB) or a period doubling (PD) bifurcation. Three bifurcation scenarios in the Nordmark map were discovered in [17], depending on the values of the bifurcation parameter. In [18] the author compared a grazing transition of a stable periodic orbit with FBs, also considering coalescing stable and unstable orbits of the impacting solution. Transitions via grazing to either chaotic attractors or periodic behavior have been studied in terms of generic bifurcation parameters in [19] and in terms of increased excitation amplitude in [18]. Banerjee et al. [20] discussed the grazing for unstable trajectories in an experimental SDOF system, also calling the phenomenon "invisible" grazing, since this grazing phenomenon was not possible to observe in the physical system. Other grazing-related features studied in SDOF VI systems include a dramatic increase in eigenvalues of unstable solutions near grazing, in [18], vanishing singularities of the Jacobian matrix in [21], sufficient conditions for periodic orbits that branch off of the grazing bifurcation [22], and comparisons of Nordmark and Poincaré maps to study grazing born from high-amplitude chaotic attractors [13]. In addition to the above-mentioned studies there are others, including [23][24][25][26] SDOF systems with linear dynamics between one-sided impacts, linear SDOFs with twosided impacts [27][28][29][30][31][32][33][34] and oblique impacts [35], as well as multi-DOF systems with impacts [14,15,[36][37][38][39][40]. Relevant information and references for the dynamics and bifurcations in non-smooth systems and related analytical methods can be found in books [4][5][6]41,42] and in the review paper [43].
In this paper we focus on the analysis of grazing and its interplay with smooth bifurcations, such as period doubling and saddle node or fold bifurcations (FB), in the two-sided vibro-impact pair. The term impact pair is used to clearly distinguish a category of VI systems where the barriers are moving with respect to the main impacting mass. In [12], and also in [16], the authors discussed the different sequences of bifurcation scenarios, including PD, FB, and grazing, in a two-degree-of-freedom (TDOF) impact pair, with an additional forcing on the mass. In [37] the authors considered a vibro-impact pair with a barrier and showed the cases when the PD bifurcation preceded grazing bifurcation in a vibro-impact-slider pair. The authors of [44,45] focused on gear VI dynamics, whose motion resembles that of a system with two moving barriers.
They demonstrated that grazing-grazing and PDgrazing trajectories can appear simultaneously in codimension 2-bifurcation scenarios, creating or destroying the gaps of periodic motion that separate chaotic regions in the bifurcation structure [44], and that chaos may die by grazing or by boundary crisis [45]. Combined analytical and numerical studies in [41] illustrate the potential for complex behaviors and transition, with various sequences of smooth bifurcations, coexisting periodic behaviors, and grazing bifurcations for the horizontal impact pair. Existence and stability conditions for different types of periodic solutions for the impact pair are given in [46] and [47], while sticking and grazing conditions for individual trajectories are discussed in [48]. For complex bifurcation sequences in SDOF systems, [49] find three scenarios of co-dimension 2 bifurcations (amplitude and frequency of the excitation) in cases where grazing may lead to a stable impacting solution and [50] considers degenerate bifurcations, where FB and/or PD bifurcations coincide with the grazing point. The variety of different options observed and reported in the abovementioned publications has been a source of inspiration for this work, where we consider the phenomenon of grazing in a wide range of transitions.
A goal of this paper is to develop analytical results for the impact pair, in particular related to the interplay of smooth bifurcations such as PD and FB and discontinuity-induced grazing bifurcations. The sequence of different types of bifurcations, as well as the settings in which both types of bifurcations occur, has clear implications for the types of solutions that can and cannot be observed in different parameter contexts. In contrast to SDOF systems for which there is a wealth of analyses, as described above, that consider the appearance of both grazing and FB or PD bifurcations, TDOF systems-and in particular the impact pair-have not undergone a systematic study. In the majority of the SDOF studies mentioned above, the grazing corresponds to the transition between nonimpacting and impacting scenarios only. In contrast, in the context of the impact pair, the grazing trajectories separate different types of impacting trajectories and consequently can appear in different sequences that include smooth bifurcations.
The model we consider in this paper takes the form of the "canonical" impact pair, where there is simple unimodal forcing on the motion of the barriers without a separate forcing on the mass. Our analysis provides a framework for studying the general structure that underpins the sequences and "competition" between smooth and non-smooth bifurcations in the impact pair. A formulation in terms of relative displacements facilitates handling the impact as a fixed switching manifold, a valuable viewpoint commonly used throughout studies of non-smooth dynamics. This flexibility allows a thorough analysis of these sequences without the tendency toward more complex behavior that is specific to other impact pair type models with additional inputs and forcing, as in [12,44,45]. At the same time, our framework provides a foundation for studying these other scenarios.
While we focus much on the general question of nonsmooth and smooth bifurcations in an impact pair, we are simultaneously motivated by the relevance of our analysis in the energy harvesting context. We apply the results to describe the dynamics and energy harvesting (EH) capabilities of a VI-EH device whose dynamics follow that of an impact pair. The device is comprised of a capsule, both ends of which are covered with flexible membranes made out of dielectric elastomer (DE) material [51], in which a ball moves freely without any friction. The capsule is excited externally, and the ball is excited by impacting one of the membranes, at the same time deforming the membranes and changing their capacitance. Thus, this concept falls directly into a category of a vibro-impact pair, since both the ball and the membranes, serving as barriers to the ball, move with respect to each other. The collection of work utilizing VI dynamics for energy harvesting before 2017 can be found in [52]. As the research in this area continues, new ideas and concepts have been emerging [53][54][55][56][57].
The paper is structured as follows: In Sect. 2 we provide a framework for studying the relationship between smooth and non-smooth grazing bifurcations, including the context of the VI-EH model and a notation for tracking different behaviors on both stable and unstable or unphysical branches. In Sect. 3.2 we give the analytical results for both smooth and grazing bifurcations. In Sect. 4 we compare the analytical results to numerical simulations, demonstrating how the analytical approach provides insight into the sequences of smooth and non-smooth bifurcations and into additional complex behaviors.

Context and framework
In this section we provide the context of the VI-EH model in which we develop our analysis of the impact pair. We also give a framework and notation for studying different families of bifurcations in the general impact pair system. The approach is then demonstrated using this framework in the context of the VI-EH model throughout the paper.

The VI-EH model
We demonstrate the dynamical analyses in the context of an inclined VI-EH system. Figure 2a shows a schematic of the device, consisting of a cylinder of mass M moving under the influence of harmonic excitation F(ωτ + ϕ) with period 2π ω and a ball of mass m freely moving inside the cylinder under the influence of gravity g = 9.8 m/s 2 . The motion of this system is described by two equations for the acceleration of the cylinder and the ball between impacts for X the position of the center of the cylinder, and x the position of the ball. We denote τ = τ j as the jth impact with either of the DE membranes located on the top ∂ T or bottom ∂ B of the cylinder. Then the ball velocityẋ(τ j ) ≡ẋ j and displacement x(τ j ) ≡ x j at the jth impact follow the impact conditioṅ where r is the coefficient of restitution and superscripts − / + denote the ball velocities just before/after each impact. This impact condition reflects the fact that the velocity of the cylinder does not change upon impact, based on our assumption that m is negligible compared to M. Periodic behaviors with different sequences of impacts are shown in Fig. 2b, alternating ∂ T and ∂ B (upper) and periodic behavior with three impacts per forcing period, two on ∂ B and one on ∂ T (lower). Non-dimensionalization of (1) and (2) with (3) allows us to streamline the analysis, by reducing the number of parameters to some key dimensionless quantities. Furthermore, we introduce a relative position (a) 40  Z (t) in terms of the difference of the non-dimensional displacements of the cylinder X * and of the ball x * as in [58,59], In (5) A is an appropriate norm of the strength of the forcing F, so that in the dimensionless setting the forc-ing has the unit norm. Then the impact condition (3) at the dimensionless time t = t j of the j th impact is given bẏ where d is dimensionless capsule length. Throughout our study we fix the parameters M = 124.5 g, ω = 5π Hz, allowing us to focus on variability of the nondimensional capsule length d as it varies with dimensional capsule length s or forcing strength A.
The formulation in terms of the relative variable Z facilitates tracking the dynamics between impacts, which take place on the fixed switching manifolds Z = ±d/2. It is straightforward to find Z between impacts, by integrating (5) over t ∈ (t j , t j+1 ) with initial conditionsŻ + j and Z + j following the impact at t j , yielding equations for the relative velocityŻ (t) and displacement Z (t) at times between two impacts. Replac-ingŻ + j withŻ − j via (6), using continuity of the displacement Z + j = Z − j , and dropping the superscripts yieldṡ where These equations form the basis of our analysis, as we evaluate them at t = t j+1 to get maps for the relative velocity and displacement between two impacts, as described in detail in Sect. 3.1. While our focus is on the analysis of the general bifurcation structure, it is valuable to track the influence of this structure on the energy output for the VI-EH device with DE membranes on both ends of the capsule. Typically DE membranes are made out of acrylic or silicon-based polymers, which are highly stretchable (300-700%), light, and non-conductive, with the dielectric permittivity coefficient being 2-12 times higher than that of air. Both sides of the membranes are covered with compliant electrodes, i.e., electrodes, which can stretch the same amount as the membranes. Thus, each membrane becomes a variable capacitance capacitor. Having applied bias voltage at the point of the maximum membrane deformation, corresponding to the maximum value of the capacitance, the excessive voltage can be harvested from the membrane as it returns to its undeformed state, corresponding to its minimum capacitance value. Thus, the amount of harvested energy is proportional to the difference in the capacitance between the deformed and undeformed states, C = 0 A m / h, where the A m and h are the areas and the thickness of the membrane at each state. Note that the features of the proposed concept align with Fig. 2a and are thus consistent with the relative motion formulation in terms of Z .
To calculate the electrical power generated by the device, we use two quantities: the average output voltage per impact U I and average voltage over a time interval U T . These depend nonlinearly on the relative dimensional velocityŻ at impact. This nonlinear relationship arises from the deformation of the DE membrane upon impact and involves geometric deflection quantities and mechanical properties of the material. These average quantities are calculated by the formulas where U k is the net output voltage at the kth impact representing the difference between the generated voltage U imp k and constant input voltage U in = 2000 V applied to the membranes, N is the number of impacts over a time [t 0 , t f ] and t f − t 0 = ω π (τ f − τ 0 ) is the corresponding dimensionless time interval. Details on the energy harvesting parameters are given in Appendix A.

A framework for dynamics of the impact pair
Both TDOF and SDOF systems experience discontinuityinduced bifurcations, such as grazing, and traditional bifurcations such as PD and FB. But the TDOF impact pair, as in the VI-EH model, has fundamentally different sequences of periodic solutions and bifurcations as compared with SDOF systems that are based on a mass impacting a fixed barrier. In SDOF systems, it is common to study bifurcations from a smooth periodic solution to a non-smooth periodic solution via grazing with the barrier. However, for the TDOF impact pair, the simplest periodic solution by definition must involve impact with at least one side of the capsule, similar to the bouncing ball. Furthermore, the sequence of bifurcations in the TDOF impact pair includes a wider range of possibilities, including grazing bifurcations generated via an additional impact per forcing period with either end of the cylinder, as well as smooth bifurcations such as PD appearing without additional impacts.
In order to track the different types of qualitative behavior in the TDOF impact pair, we employ a notation that captures different types of impacting behavior. For example, increasing the forcing amplitude A in the impact pair yields a basic sequence of bifurcations from a low amplitude state, starting with an impact on only one end per forcing period T , to a single impact on each end or to combinations of multiple impacts on each end that may or may not have period T . To differentiate across the range of impacting solutions, we use the notation n:m/ pT to categorize different periodic motions for an impact pair with a T -periodic external excitation, as in [59]. The integers n and m indicate the number of impacts against the bottom and top, respectively, of the capsule per period T , simplified to n:m for the cases where p = 1. This notation also captures sequences of period doubling sequences, such as transitions from 1:1 to 1:1/ pT for p = 2, 4, . . ., as well as non-smooth bifurcations via grazing. Grazing occurs for the critical parameter value where the trajectory has a zero relative velocity at impact, so that beyond this critical value the system dynamics includes an additional impact (with nonzero impact velocity). The notation reflects this as we observe transitions from n:m behavior to (n + 1):m or n:(m + 1) behavior in the impact pair. For example, as an impact pair, the VI-EH system exhibits favorable energy output following a grazing transition from 1:0 to 1:1 behavior, that is, where the system moves from impact on only one end of the capsule to alternating impacts on both ends. This transition is achieved, e.g., through increased forcing amplitude, with further increases driving transitions from 1:1 to either 2:1 or 1:2 behavior via what we term first-order grazing on the bottom or top of the capsule. Alternatively, it may undergo a period doubling transition to 1:1/2T prior to any grazing bifurcation. For completeness, we also include the notation n:m/C, indicating complex, aperiodic behavior, where over long time periods the typical number of impacts against the bottom (top) of the capsule is n(m). Our  Table 1 computational results capture the n:m/C behaviors, but we do not study them analytically in this paper. Figure 3 provides motivating examples of different sequences of bifurcations that commonly appear for the impact pair (5)-(6) for both smooth and grazing bifurcations. Here and throughout the paper we take F = A cos(ωτ + φ) in (4)- (5). The markers on the bifurcation diagrams give values for the impact velocitiesŻ j obtained numerically for large t, that is, for the attracting or stable behavior. These numerical results are obtained by a continuation-type method for d decreasing, choosing a value of d, computing over a sufficiently long time to reach the attracting behavior, from which we obtain the sequence of values ofŻ j at that value of d. This attracting behavior provides the initial condition for computing the attractor at the next value of decreasing d, usually chosen close to the previous value of d. From (6), d varies linearly with dimensional capsule length s and is inversely proportional to forcing amplitude A. Then Fig. 3 shows the results for d decreasing in the continuation-like method, with either decreasing s or increasing A. The notation from Table 1 is used to highlight different types of transitions. Some select analytical results, specifically stable and unstable 1:1 and 2:1 solutions and the PD and G ν:μ n:m bifurcations, are shown in Fig. 3 to provide a more complete picture of the motivating examples. These analyses leading to these results are discussed in detail in Sect. 3.
We note that the detailed structure of the bifurcation diagram follows from the choice to graph relative impact velocity vs. a bifurcation parameter. This choice fits naturally with the expressions forŻ j and Z j in (7)-(8) which form the basis of the maps between impacts defined in Sect. 3.1. One could also choose to use a stroboscopic (Poincaré) map, as in [60] where it is shown to yield a similar bifurcation structure for the VI-EH model as shown in this paper. One difference is that there may be additional discontinuities in the 1:1 branches in the stroboscopic map, resulting from the shift in the asymmetry of the 1:1 solution. That is, the phase of the impacts relative to the periodicity of the stroboscopic map may result in additional discontinuities in the bifurcation diagram besides those associated with grazing or other bifurcations. Figure 3a shows grazing bifurcations for decreasing s, with transitions from 1:1 to 2:1 and 2:1 to 3:1 solutions, denoted by G 2:1 1:1 and G 3:1 2:1 , respectively. Panel (b) shows a period doubling bifurcation from 1:1 to 1:1/2T , labeled PD, with the transition labeled PD G from 1:1/pT to 2:1 behavior following additional period doubling, and finally G 3:1 2:1 . Panel (c) shows PD transitions from 1:1 solutions to 1:1/2T and from 2:1 to 2:1/2T for smaller d, followed by PD G from 2:1/ pT to 3:1.
In between there is a sequence of period doublings of 1:1/ pT behavior, then a window of n:m/C which terminates in 2:1 periodic behavior. Panel (d) is qualitatively similar to Panel (c) for d < .4, but there are two separate regions of stability for the 1:1 solution. For d > .8, the 1:1 solution is stable, followed by a sequence of PD to reach 1:1/ pT and eventually chaotic behavior as d decreases, until a second stability region for the 1:1 solution is reached for d < . 4.
Furthermore, the results in Panels (a)-(d) illustrate several types of transitions to 2:1 behavior, where the low velocity branch of the stable 2:1 solution does not reachŻ j = 0. In Panel (a) at the critical point G 2:1 1:1 , the 1:1 solution meets an unstable 2:1 branch (not shown in Fig. 3), which in turn meets the stable 2:1 solution at a fold bifurcation (FB). A similar transition was observed for a SDOF system in [50]. Figure 3 illustrates that the TDOF impact pair exhibits a variety of transitions in which the FB structure plays a role. In addition to the 1:1 to 2:1 transition in Panel (a), it is present in the transition to 2:1 behavior via PD G in Panel (b), and following n:m/C in Panel (c) via a transient grazing TG defined in Table 1. For simplicity we do not label FB and TG in Fig. 3; the mechanisms and significance of these features in the context of grazing transitions are discussed in detail in Sect. 4, and labeled in the figures there. Table 1 gives the notation to identify and distinguish first-order grazing transitions G ν:μ n:m , period doubling PD, and grazing from period doubled PD G , as shown in Fig. 3. These are all transitions from a stable solution to a different stable state, motivating the analytical study of critical transitions as demonstrated in Sect. 3. In Table 1 we introduce the term first-order grazing, since in the impact pair setting, different types of grazing typically yield transitions that add one impact per forcing period. We also introduce the term higherorder grazing generically for a collection of different types of transitions that can add multiple impacts per forcing period, usually limited to the larger n and/or m settings. While we do not study higher-order grazing analytically, we point out some occurrences observed numerically in Sect. 4.
Beyond G ν:μ n:m and PD, the analysis allows identification of critical parameter values on unstable or unphysical branches. Specifically, we identify firstorder ghost grazing,G ν:μ n:m that is, grazing of unstable solutions, typically those that have lost stability via PD. Furthermore, we identify ghost period doubling PD on unphysical branches, occurring on the dotted branches following grazing, e.g., on the dotted line for the unphysical 1:1 solution shown in Panel (a) for d below G 2:1 1:1 . We refer to these as unphysical solutions since they violate the containment condition, which is | Z (t) | < d/2. Similar terminology was used in [32,50] in a SDOF system, and we discuss it further in Sect. 3.2. These ghost bifurcations are defined in Table 1, and their significance in the dynamical behavior and energy output of the VI-EH device implications are discussed in Sect. 4. One delicate difference to note: The value G ν:μ n:m may correspond to a bifurcation to an unstable ν:μ solution, which generically coexists with a stable ν:μ solution, so that the transition observed numerically is from a stable n:m to the stable ν:μ solution. This was mentioned above in the context of FB for the critical point G 2:1 1:1 in Fig. 3a. This G ν:μ n:m transition contrasts with ghost grazing, which corresponds to grazing of an unstable n:m solution, as discussed in detail in Sect. 4.1.
Grazing on stable n:m branch, witḣ Z j = 0 at impact on ∂ T or ∂ B; generically indicates transition to ν:m periodic solution with one additional impact per period T Grazing from period doubled state PD G Transition via grazing from a n:m/ pT periodic solution to a Tperiodic solution that has an additional impact per period. For example, PD G is used for transitions from 1:1/ pT to 2:1, commonly observed for p > 1 In the following we develop a systematic approach based on the (semi)-analytical results for studying transitions between n:m/ pT solutions for an impact pair. The analytical approach provides a new structure for studying sequences of different solutions on stable and unstable branches, filling a gap in the study of TDOF impact pair-type systems. The new notation provides a framework for tracking the bifurcation sequences of grazing, period doubling, and fold bifurcations, even for the unstable or unphysical solutions, since these bifurcation ghosts contribute to the routes to complex dynamics, as discussed in Sect. 4.2. The added value of our general semi-analytical study of the dynamics is its applicability to explain and predict energy output in the VI-EH system. For example, transitions from n:m/ pT to (n + 1):m solutions frequently appear via grazing bifurcations, with changes in average energy output resulting from the additional impact per forcing period. With the VI-EH system as a motivation for our analysis, we focus on smaller values of n and m, for which that system exhibits larger parameter ranges of higher energy outputs.

Analytical results for bifurcations
We outline the analytical approaches for identifying smooth bifurcations for n:m periodic solutions for the impact pair, with explicit equations for 1:1 and 2:1 solutions as in [59]. We also provide the analytical approach for identifying first-order grazing transitions. Then these results are combined to capture sequences of the different PD and grazing bifurcations, including their corresponding ghost bifurcations. All results are illustrated in the context of the VI-EH model (7) and (8). We note that similar analyses were applied to a slowly varying subsystem in a vibro-impact nonlinear energy sink, based on an impact pair (see [40] and references therein.)

Maps for periodic solutions
In [58,59] analytical descriptions of n:1 periodic solutions were obtained, based on compositions of maps from impact to impact. There are four basic nonlinear maps P , = 1, 2, 3, 4 for the corresponding transitions between impacts, The mathematical expressions for these maps take different forms depending on whether Z j and Z j+1 are located on either ∂ B or ∂ T . Specifically, we take t = t j+1 in (7)- (8) to get the equations forŻ − j and Z − j , and dropping − as above, we obtain the maps for P , where A n:m periodic solution is then composed of the sequence of n + m maps P n−1 • P 3 that are completed during one forcing period of length T . The sequence involves a transition from ∂ B and ∂ T (P 2 ) and vice versa (P 3 ), as well as consecutive impacts on either ∂ B (P 1 ) and ∂ T (P 4 ). A diagram of such sequences is shown in Fig. 4.
Given the inclination angle β > 0, it is not difficult to imagine that (4)-(6) regularly generate n:1 solutions for n ≥ 1 of the specific form P n−1 1 • P 2 • P 3 , which we focus on here. The maps form the basis for deriving an analytical form of the periodic motion in terms of several quantities: the impact velocitiesŻ k+ , = 0 . . . , n, the time intervals between impacts t k+ −1 = t k+ − t k+ −1 for the n + 1 maps, and the phase difference at the initial impact ϕ k = mod(π t k + ϕ, 2π). The approaches used in [58,59] determine the n + 2-tuple which define the n:1 periodic solution. The n + 2 equations for (13) are obtained through the steps i) summing the equations forŻ j for all P j ; ii) summing the equations for the Z j for all P j ; iii) n − 1 pairwise combinations of the equation forŻ k+ j with that of the impact position Z k+ j+1 , applied successively for j = 1, . . . , n − 1; and iv)Ż k from the equation for position Z k+1 in P 1 . These n + 2 equations take the formŻ k = h (ϕ k , t k , . . . , t k+n−1 ), and the remainingŻ j for j = k + 1, . . . , k + n and t k+n are determined, respectively, from combining (13) with (11) for the P j 's, and with the periodicity conditions n+1 =1 t k+ −1 = T andŻ k+n+1 =Ż k . Note that for this formulation, the analytical solution for the n + 2tuple yields the initial impact velocityŻ k and phase shift ϕ k for the first map in the sequence P n−1 1 • P 2 • P 3 , which is P 1 for n > 1, and P 2 for n = 1.
We focus primarily on transitions from either 1:1 or 2:1 periodic solutions and recall the analytical results that define these two types of periodic solutions for the forcing F = A cos(ωt + φ). The three equations for the triple K * 3 ≡ (Ż k , ϕ k , t k ) that determine the 1:1 periodic solution are [58], The four equations for the quadruple K 4 ≡ (Ż k , ϕ k , t k , t k+1 ) that determine the 2:1 periodic solution are [59], The linear stability analysis of the n:1 periodic solution defined by the n + 2-tuple (13) is based on the linear techniques of [8,41,61]. The analysis provides a linear equation for the small perturbation δH k to the fixed point H * k = (t k ,Ż k ) of the n:1 periodic solution, δH k+(n+1) = J n+1 (H * k+i,i=1,...n )δH k . The matrix J n+1 is based on the product of Jacobians from the composition of maps P j evaluated at the appropriate fixed values H * that define the n:1 solution.
The details for computing J n+1 are given in [58,59] for the 1:1 and 2:1 periodic solutions. The eigenvalues λ i , i = 1, 2 of J n+1 then give the linear (in)stability of the n:1 solution, as well as indicating the nature of the (loss of) stability. For example, |λ i | = 1 yields a PD bifurcation (specifically, λ i = −1), resulting in a transition from 1:1 to 1:1/2T periodic solutions. Figure 6 shows the values d = d PD for transitions to 1:1/2T and 2:1/2T behavior. Depending on the parameter values, there can be additional period doublings to 1:1/ pT and 2:1/ pT solutions for p ≥ 2. We do not study these analytically, but they are observed in the numerical simulations of (5)-(6) shown below.

Analytical results for grazing bifurcations
In order to compare the prevalence of transitions from 1:1 to 1:1/ pT solutions via PD vs. transitions from 1:1 to 2:1 solutions via grazing, we give analytical conditions for determining first-order grazing bifurcations. The linear stability analysis [58,59] based on the eigenvalues for the linearized map composition does not apply in this case, since it assumes the composition of smooth trajectories. In contrast, grazing is a discontinuity-induced bifurcation, corresponding to the critical parameter value at which a smooth trajectory is tangent to the switching surface. Then the location of the grazing bifurcation is obtained by using the maps P j above, combined with conditions for intersection of a trajectory with a switching surface.
As discussed above, analytical studies of grazing in SDOF impacting systems usually consider the first instance of grazing with a single barrier. Then the focus is on the transition from a non-impacting smooth periodic behavior to one that has a zero-velocity impact on the barrier. In contrast, for the TDOF impact pair system, our focus is on the transition between different energy producing n:m/ pT behaviors, for n, m ≥ 1. Therefore we study G ν:μ n:m andG ν:μ n:m for n, m ≥ 1 as described in Table 1 above. We restrict our analysis to first-order grazing as defined there, where one additional impact per period T follows from the grazing.
In particular, we seek the value d G corresponding to grazing of n:m periodic solutions. We apply a condition based on (7), which gives an analytical expression for the behavior ofŻ between impacts on the 1:1 trajectories. In this paper we focus on the condition for the specific type of grazing that occurs when there is a zero-velocity contact on the same surface as that of the previous impact. Then d G is obtained from the generic condition, This definition of d G states that, following an impact on ∂ B (∂ T ) at t j , the ball returns to the same end of the cylinder at t G with zero impact velocityŻ (t G ) and then continues to the next impact t j+1 with nonzero impact velocity. The condition (22) follows from evaluating (7) at t = t G and combining with the requirement thatŻ (t G ) = 0. The conditions in (23) ensure that the remainder of the trajectory starts and ends with nonzero impact velocities on either end of the capsule.
An example of a grazing flow map decomposition is shown in Fig. 5a. The complete orbit for the 1:1 solution is given by the composition P 2 • P 3 . The dashed line shows the map P 2 for d > d G , describing the motion of the ball from point D on ∂ B to point A on ∂ T (no grazing). For d = d G , the solid line shows P 2 that connects points D and A and includes point G wherė This example corresponds to the first-order grazing transition G 2:1 1:1 , for example, as observed in Fig. 3a. Calculating d G from (22)- (23) for this case uses the triple K 3 = (Ż k , ϕ k , t k ) from (14)- (16), which defines the 1:1 solution composed of P 2 • P 3 . Then d G is determined by using this triple in (22). That is, in (22),Ż j =Ż k , ϕ j = ϕ k is used in the definition of F 1 , t j+1 = t j + t k , and the "+" sign is chosen since grazing occurs on the P 2 map, with Z j = Z (t G ) = d/2 on ∂ B. In the calculation we take t k = 0 without loss of generality, since the calculation is based on the 1:1 periodic solution. Then we integrate (7) with these values to obtain d = d G from (22). Values d G for G 2:1 1:1 andG 2:1 1:1 are shown for different values of r and β in Fig. 6a, b.
We also obtain analytical results for first-order (ghost) grazing transitions from 2:1 solutions, of which there are two types: G 3:1 2:1 (G 3:1 2:1 ) as shown in Figs. 3a, b and 7 and G 2:2 2:1 (G 2:2 2:1 ), as shown in Fig. 8. To compute d G from (22), we use the values from the quadruple K 4 = (Ż k , ϕ k , t k , t k+1 ), obtained by solving (17)- (21), which defines the 2:1 solution given by the composition P 1 • P 2 • P 3 . Specifically, for G 3:1 2:1 , grazing occurs when the P 2 map reaches ∂ B with vanishinġ Z j . To compute d G for G 3:1 2:1 , in (22) we setŻ j =Ż k+1 , obtained from K 4 and the appropriate map to getŻ k+1 . The phase ϕ k from K 4 is used in the definition of F 1 , t j = t k+1 and t j+1 = t k+1 + t k+2 , with t k = 0, based on the T -periodic 2:1 solution. The "+" is chosen for Z (t j ) = Z (t G ) = d/2 on ∂ B. For G 2:2 2:1 , grazing occurs when the P 3 map reaches ∂ T withŻ j = 0 . To compute d G for G 2:2 2:1 , in (22) we setŻ j =Ż k+2 obtained from K 4 and the appropriate maps to finḋ Z k+2 . The phase ϕ k from K 4 is used in the definition of F 1 , t j = t k+1 + t k+2 and t j+1 = T , with t k = 0. The "-" sign is chosen for Z (t j ) = Z (t G ) = −d/2 on ∂ T . The analytical results for d G corresponding to G 3:1 2:1 , G 3:1 2:1 , G 2:2 2:1 , andG 2:2 2:1 are shown for different values of r and β in Fig. 6c, d. While in theory G 1:3 1:2 and G 2:2 1:2 are also possible, in general we do not observe these for the parameter ranges of this study, where β > 0 and where the restitution coefficient is the same for both ends of the capsule. The grazing transitions G 1:3 1:2 and G 2:2 1:2 are observed, for example, when the ratio of the restitution coefficients of ∂ T to ∂ B is less than unity [62].
The analytical calculations for smooth bifurcations and for grazing allow us to follow unstable or unphysical n:m solutions, determining G or PD transitions on these branches. Figure 5b illustrates these types of solutions, showing stable 1:1 behavior, an unstable 1:1 solution following a PD bifurcation, and an unphysical 1:1 solution violating containment, following G 2:1 1:1 . Even though these critical points are on unstable or unphysical solutions, they have implications for the impact pair dynamics and energy output of the VI-EH device.

Comparisons of grazing versus period doubling bifurcations
Based on the results obtained from the analytical approaches described above, we compare the sequences of the grazing and PD bifurcations. In Fig. 6a, b, for different values of r , we compare values d G at which there is a first-order grazing transition G 2:1 1:1 with the values d PD for the PD bifurcation from 1:1 to 1:1/2T . In Fig. 6c, d we compare values d G for which we have a first-order grazing transition from a 2:1 solution, either G 3:1 2:1 or G 2:2 2:1 , vs. d PD for the PD bifurcation from 2:1 to 2:1/2T . The analytical results are shown in the d vs. r plane for different values of β, illustrating the relative location of the grazing and period doubling bifurcations, and their ghost counterparts. For example, if d G > d PD for a given value of r in Fig. 6a, b, then as we decrease d, the stable 1:1 solution undergoes a G 2:1 1:1 at d G , while the PD bifurcation of 1:1 to 1:1/2T is not observed, thus leading to a ghost PD. If d G < d PD in Fig. 6a, b, then PD from stable 1:1 to 1:1/2T behavior is observed, while the analytical grazing result occurs on the unstable portion of the 1:1 branch, corresponding to a firstorder grazing ghostG 2:1 1:1 . Similar conclusions follow from Fig. 6c, d for the 2:1 solution.
Combinations of the Panels (a)-(b) with (c)-(d) in Fig. 6 show the analytical results for the bifurcations observed in Fig. 3, which we recall here.     Figure 6b shows the largest d in the sequence of d PD from 1:1 to 1:1/2T , where we see that this range is substantially larger for β = 7π/18, π/2 than for smaller β. While not shown in Fig. 6a or b, there are ranges of d for which there are separate sequences of 1:1 solutions, followed by 1:1/ pT and in some cases 1:1/C, for larger values of r and other values of β, as discussed further in Fig. 10 where the analytical results are compared with numerical transitions to 2:1. The complex behavior observed in Fig. 3d is consistent with the multiple 1:1 periodic solutions and larger ranges of bi-stability between different types of periodic solutions that generically appear for larger values of r [62]. We do not pursue a detailed analysis of these here, given that larger values of r are less common for the VI-EH device of interest. However, for decreasing s there is no sequence of different 1:1/ pT bifurcations, as observed for increasing A.
The relative location of the curves in Fig. 6 indicates that for some values of parameters we cannot observe first-order grazing at all, for decreasing d, while for other parameters, grazing dominates the transitions. For larger β and r and decreasing s shown in Fig. 6a, c, PD bifurcations are generally observed, rather than first-order grazing. Then the grazing transitions shown correspond to ghost grazing. For smaller values of β, the grazing bifurcations of periodic 2:1 solutions tend to beG 2:2 2:1 or G 2:2 2:1 , which follows from the inclination angle closer to the symmetric case of β = 0. The prominence of PD andG for larger values of β has been observed in other studies [62], also showing G and PD for certain combinations of asymmetric r (different r on ∂ B and ∂ T ) not shown here. In addition to the influence of β, results in [59] demonstrated that changes in s for the varying A case can shift the PD curves as shown in Fig. 6b, d. We highlight a number of subtle points that are discussed in detail in the next section. The grazing bifurcation G ν:μ n:m may correspond to an unstable ν:μ branch which coexists with a stable ν:μ branch via a FB, as shown in Figs. 7 and 8 for ν = n +1 and μ = m. In that case the low velocity impactŻ k may be bounded away from zero. These transitions may also be associated with bi-stability of the n:m and ν:μ branches. In that case, hysteresis may be observed numerically, i.e., the observed transitions may be shifted relative to the fixed analytical results for G ν:μ n:m and PD shown in Fig. 6. That is, using a continuation-type method for decreasing d, the n:m branch is followed until the G ν:μ n:m transition is reached, while for increasing d, the ν:μ branch may be followed until it loses stability, e.g., at an FB above that of the G ν:μ n:m . This hysteresis is not shown in Figs. 3, 7, and 8, since the numerical continuation-type method is focused on decreasing d. A final observation from Fig. 6 is that the analytical results for G ν:μ n:m rely on a geometric analysis, in contrast to the analytical results for smooth bifurcations, based on the eigenvalues. This difference is typical of non-smooth bifurcations, which in general cannot be obtained via the techniques for smooth bifurcations.   Table 1. Here we discuss the sequence of these critical points and transitions, as is consistent with Fig. 6. Below (Sects. 4.1 and 4.2) we compare and contrast the multiple routes to grazing and the implications of these critical points.
As shown in Fig. 7a the PD transition from 1:1 to 1:1/2T occurs for larger d than forG 2:1 1:1 ghost grazing. Decreasing d belowG 2:1 1:1 in Fig. 7a, the system follows the stable 2:1 solution, until it undergoes an additional grazing at G 3:1 2:1 . Thus, even though the eigenvalues of the 2:1 solution |λ i | < 1 between G 3:1 2:1 and PD (see Fig. 9a), that solution is unphysical as it violates containment following grazing, analogous to Fig. 5. Figure 7a indicates that the resulting 3:1 solution is stable at G 3:1 2:1 , with the fourth branch of this solution emerging fromŻ k = 0 at this point. The stability of the 3:1 solution is also illustrated in the phase plane and time series in Panels (k) and (l).
In contrast to Fig. 7a, in which s decreases, there is a different sequence of PD and G bifurcations (and other transitions) in Fig. 8a, in which A increases for decreasing d with slightly smaller r . Figures 8b, c illustrate G 2:1 1:1 , which, as d decreases, is reached before PD. Then the ghost PD of the 1:1 solution occurs for smaller d and is not observed in the dynamics. At the grazing point G 2:1 1:1 there is an unstable 2:1 solution (magenta dashed line) that bifurcates subcritically from the 1:1 solution, and there is a FB at which a stable 2:1 solution meets the unstable one. Corresponding to this subcritical bifurcation, the value of d for FB is greater than that for G 2:1 1:1 . For d below G 2:1 1:1 the 1:1 solution is unphysical (red trajectory in Panel (d)) and the solution moves to the stable 2:1 solution (blue trajectory in Panel (d), time series in (e)). Figure 9b shows the eigenvalues for the two 2:1 solutions that exist for d between G 2:1 1:1 and FB. The unstable branch of 2:1 solutions has an eigenvalue |λ i | > 1, while for the second branch both eigenvalues |λ i | < 1 for values of d between FB and PD. Following PD there is a loss of stability from 2:1 to 2:1/2T behavior, as shown in Panels (f)-(g). Consistent with Fig. 6 for smaller β, we observeG 2:2 2:1 in Fig. 8h-j, in contrast with G 3:1 2:1 (orG 3:1 2:1 ) solutions determined for larger β. Since PD for 2:1 precedes the grazing for decreasing d, then the ghost grazingG 2:2 2:1 (red trajectory in Panel (i), with time series in (h)) is obtained theoretically, while the dynamics attract to a more complex alternating 2:1 and 2:2 solution with period 2T (blue trajectory in Panel (i), time series in Panel (j)).

Multiple routes to first-order transitions
Comparison of Figs. 7a and 8a illustrates four different mechanisms of first-order transitions that increase the number of impacts per forcing period by one, e.g., . These all involve grazing or ghost grazing bifurcations, and we list the different signature for each that is apparent in the bifurcation diagram. Some of these types of non-smooth bifurcations have been observed in other studies, and we highlight these connections as well.
1. The stable n:m solution experiences first-order grazing at G n+1:m n:m or G n:m+1 n:m , at which an additional impact withŻ k = 0 occurs. If this solution is stable, then an additional low velocity bifurcation branch, initiated atŻ k = 0, is also observed in the numerics, e.g., as in G 3:1 2:1 in Fig. 7a, k, and l. 2. In contrast to the previous item 1., it is also possible that an unstable (n + 1):m (or n:(m + 1)) solution follows from the first-order grazing of a stable n:m solution, yielding G n+1:m n:m or G n:m+1 n:m . This type of grazing bifurcation has been called "invisible" grazing in the experimental context in [20], since the grazing was predicted analytically but not observable, presumably because of the instability. That case commonly yields a scenario where both a subcritical unstable bifurcating branch and a stable branch meet at a  Fig. 7, where the FB, corresponding to the birth of a stable 2:1, quenches the 1:1/C behavior. A similar transition to a 2:1 solution following 1:1/C appears in Fig. 3c. Similar transitions were observed in [44], explaining some of the gaps in chaotic behavior for a model of gear dynamics. This model is similar to that of the impact pair, but it includes additional forcing frequencies that facilitate larger windows of chaos.

Implications for bi-stability, transient behavior, and complex dynamics
The analytical approach provides knowledge about the sequence of smooth bifurcations, e.g., PD and FB, and non-smooth grazing bifurcations. While this analysis does not provide exact results for the more complex bifurcations, it nevertheless provides insight into the birth and death of these more complex dynamics. We highlight this insight in a few specific contexts. Ghost grazing: As mentioned above in the context of different routes to grazing, in instances where there is PD of the stable solution andG of the unstable solution, typically PD G or TG leads to a transition from n:m to (n + 1):m or n:(m + 1) for values of d abovẽ G. The basis for this conclusion is illustrated in the phase planes where we observe that a trajectory for the n:m/ pT solution is naturally closer to grazing than the n:m trajectory. This follows from the "loops" in the phase plane trajectories as shown in Figs. 5 and 7d. For example, in the map P 2 from ∂ B → ∂ T , these loops correspond to intervals of increased Z as the bottom of the capsule approaches the ball. Since PD leads to multiple loops in the phase plane trajectory, one of these loops on the n:m/2T solution must take values of Z (t) closer to d/2 (∂ B) than the values of Z (t) on the P 2 map for the n:m solution (compare Fig. 7b, d). Then, the resulting loops from n:m/ pT or from n:m/C yield PD G or TG, respectively, at a larger value of d than for ghost grazingG ν:μ n:m of n:m solutions. Thus being able to locate ghost grazing analytically provides a lower bound on the critical parameter value at which the n:m-type solution transitions to one with an additional impact per period. Figure 10 illustrates this for different parameter combinations as d decreases, e.g., PD G or TG transitions to 2:1 behavior occur at or beforẽ G 2:1 1:1 , and transitions to 3:1 or 2:2 periodic solutions precedeG 3:1 2:1 andG 2:2 2:1 , respectively. Fold bifurcations and bi-stability: The analysis allows the calculation of both unstable and stable branches of n:m solutions, and with the corresponding stability analysis we can identify any FB points. The location of FB indicates the potential for either intermediate regions of TG or bi-stability between different periodic solutions.
Regions of bi-stability occur for grazing transitions to unstable (subcritical) solutions, as shown in Fig. 8a. A typical signature in a computationally or experimentally observed transition from n:m or n:m/ pT to (n + 1):m is that the additional impact has nonzero impact velocity at the transition. This feature suggests an underlying grazing bifurcation to an unstable subcritical (n + 1):m, which in turn meets the stable (n + 1):m solution at FB. These transitions, in the context of bi-stability, may occur from G or from PD G as shown in Fig. 3a, b. An example of such a grazing point is given by G 2:1 1:1 in Fig. 8a. While the location of G 2:1 1:1 is a fixed value determined analytically, the transition between 1:1 and 2:1 solutions observed numerically may occur at a different value due to the bi-stability. For example, one could use the numerical continuationtype method to follow the stable 2:1 branch for increasing d starting below G 2:1 1:1 . Then it is possible to follow this branch to FB rather than reaching the 1:1 solution at G 2:1 1:1 , as was observed for decreasing d. This hysteresis in the numerical results illustrates how the analytical results are necessary to provide the full picture, independent of increasing or decreasing d.
In contrast to the bi-stable case, TG tends to occur in the absence of bi-stability, i.e., when there is no overlap of stable n:m/ pT and (n + 1):m/ pT behavior. Then grazing may follow the sequence of PD's leading to n:m/C, so that the grazing occurs intermittently on a complex trajectory, without a sustained first-order grazing transition. The contrast can be seen by comparing transitions from 1:1/ pT to 2:1 behavior in Fig. 3b and from 1:1/C to 2:1 in Fig. 3c, as well as transitions from 1:1/C to 2:1 in Fig. 7 and from 1:1 to 2:1 in Fig. 8. Larger ranges of bi-stability are prevalent for certain combinations of larger r and β, with other factors such as increasing A also playing a role. These aspects have been observed in other studies [58], [62], but we do not explore them here.
Ghost period doubling: For smaller r , we observe ranges of β in Fig. 6 where PD of the n:m solutions occur well below G ν:μ n:m , possibly with PD following G ν:μ n:m yielding ν:μ/ pT . In this setting PD of n:m does not influence the ν:μ behavior, since the n:m and ν:μ are different periodic solutions with different dynamics. This separation follows from the dynamical route to grazing, where the bifurcation is induced purely by the discontinuity of impact rather than by eigenvalues or other smooth parametric variation.
The mechanism of a PD, as an unphysical solution, suggests changes in the device design to avoid (or facilitate) G occurring before the PD. This is particularly relevant for parameter ranges where the PD and G ν:μ n:m are in close proximity to each other, near the intersection of different colored curves as shown in Fig. 6. Increased r or β near these intersections could yield a PD rather than PD, which may or may not be desirable. In cases where PD and G ν:μ n:m (nearly) coincide, there may be complex transients or chaotic behavior. For example, in Fig. 10c, where the PD of the 2:1 solu-tion and G 2:2 2:1 occur in close proximity, there are complex higher-order grazing transitions for small windows of r near these bifurcations. Figure 6 shows that the proximity of PD and G depends on how d changes with its underlying dimensional quantities (s and A), indicating multiple scenarios where small variation in these parameters can advance PD ahead of the G ν:μ n:m . Figure 10 combines the analytical results from Fig. 6 together with numerically obtained transitions. As described above, we see that PD G and TG transitions are bounded below by the analytically determined ghost grazingG. These transitions do not necessarily depend on r continuously, as seen for the PD G /TG transitions to 2:1 solutions in Panels (b) and (c). Recall that for larger β and r with decreasing A, there are two separate ranges of d with stable 1:1 solutions, as discussed for Fig. 3d. In these cases there is a larger range of d over which 1:1, 1:1/ pT , or 1:1/C solutions are sustained, so that the PD G or TG transition to 2:1 behavior is shifted to smaller values of d as discussed in the contexts of Figs. 3 and 6. Similarly, in Fig. 10c there is a range of d below the PD of the 1:1 curve for which 1:1/ pT and 2:1 solutions are bi-stable. In the figures we give the PD G /TG curve for decreasing d, indicating quenching of 1:1/ pT or 1:1/C in this setting. Furthermore in Panel (c) for larger r , fluctuations in the PD G transition to 2:2 correspond to some narrow windows of bi-stable 3:1 transitions, apparent with slight variation in initial conditions but not shown in these figures. As discussed above, here we focus on the case of decreasing d, but in the case of bi-stability, the numerically obtained transitions may shift for increasing d due to hysteresis. Other influences of complexity appear in Fig. 10a, where larger r and β contribute to fluctuations in transitions to 3:1 from 2:1/ pT via PD G or from 2:1/C following TG, indicating dependence on initial conditions and numerical sensitivities in the simulations. Additional irregularities appear in Panel (c), near the transition from G 3:1 2:1 to PD of 2:1 and PD G to 2:2 values (near r = .2). These follow from the coincidence of PD of the 2:1 solution and G 2:2 2:1 , discussed above in the context of PD. Indeed, a closer examination of the solution near this transition indicates a narrow region that takes values very close to that of 3:1 behavior, but in fact is of the form of alternating 3:1 and 2:1 behavior over a 2T period.

Implications for the VI-EH device
The PD, G ν:μ n:m , PD G and TG transitions all have implications for energy output in the VI-EH device. In general, we see that any function dependent on the relative impact velocityŻ k is influenced by these transitions, e.g., any type of energy measure that is dependent oṅ Z k , as shown in A. Each of these types of grazing bifurcations yields abrupt changes inŻ k directly, as well as changes in the underlying dynamics. In addition, the stability or complexity properties of these solutions influence the robustness of the behavior. Figures 7a and 8a give the bifurcations ofŻ k in terms of d in two different settings: for decreasing s, and increasing A. We emphasize that this is not just a matter of d increasing with s or decreasing with A, but rather the locations of the bifurcations are different in these two cases; this follows from the fact that increasing A not only reduces d, but also reduces the influence of gravity, as evidenced from the non-dimensional gravity termḡ in (5). Likewise the energy output has different sequences in the two different scenarios, following from the different bifurcations which indicate critical parameters at which the energy output may change abruptly. Figure 11 shows the energy output corresponding to the bifurcations in Figs. 7a and 8a in terms of the average output voltages per impact U I and per time interval U T , both of which exhibit abrupt changes at parameter values corresponding to different grazing transitions. These changes follow directly from the definitions of U I and U T in (9), since via each grazing transition, an additional low velocity impact is added per forcing period. In the transitions from n:m to (n + 1):m or n:(m + 1) behavior, typically the impact velocity for the additional impact decreases with n, m, leading to larger drops in U I , while the U T increases with successively smaller jumps. In contrast, following the PD of the 1:1 solution shown in Fig. 7a, the corresponding U I and U T shown in Fig. 11a change gradually for successive PD solutions, since the branches forŻ k of the 1:1/ pT solutions are changing gradually. Likewise for a sequence of PD leading to complex behavior, as in Fig. 8a for small d < 0.12, there is a gradual drop in U I for complex behavior that includes several impacts with smallŻ k . Then there is a discontinuous decrease in U I , shown numerically in Fig. 11b, following the subsequent higher-order grazing transitions to 3:2, 4:2,  Figure 12 compares the influence of the different types of bifurcations on the energy output, for smaller and larger values of r . Figure 6 indicates that the sequence of bifurcations in terms of d can be completely different for smaller or larger r . These differences influence the sequence of energy output. For decreasing s shown in Fig. 12a, for smaller r Fig. 12b this ordering is different, as follows from the relative positions of G, PD, and PD G transitions shown in Figs. 6 and 10. The G 2:1 1:1 transition for smaller r occurs for smaller d than the PD G transition from 1:1/ pT to 2:1 for larger r , while the PD G transition to 2:2 behavior for smaller r occurs for larger d than for the PD G with larger r . Furthermore, increasing A (smaller d) drives sequences of PD and more complex solutions, particularly for larger r . As shown in Fig. 12b for smaller d, U I increases for complex solutions where larger r expands the range ofŻ k , while for smaller r these complex solutions add more frequentŻ k near 0, so that U I decreases.
While the differences in U I and U T follow naturally in regions where different values of r yield different n:m/ pT solutions for the same value of d, it may seem counterintuitive that the energy output for the different values of r is comparable for the same type of n:m behavior, or in some cases larger for smaller r . This observation, apparent in U T in Fig. 12a and in both U T and U I in Fig. 12b, follows directly from the phase differences φ k of impacts, relative to minima of the capsule trajectory, as shown in Fig. 2b. Impacts that occur away from the minima, where the capsule and ball move in opposite directions, generate largeṙ Z k . In contrast, for impacts that occur closer to these minima, with smaller differences in the capsule and ball velocities, the relative impact velocitiesŻ k are smaller. Typically larger r favors solutions where the impact is near the minima of the capsule trajectory, thus yielding reduced energy output.
We recall the analytical results for PD, G 2:1 1:1 and G 3:1 2:1 from Fig. 6 highlighting their connection to the influence of device features, such as its length s, the angle β and the restitution coefficient r . Since in a real application the excitation amplitude is given, the device and its energy performance should be designed using the three parameters mentioned above. Obviously the 1:1 response is the most beneficial EH option, thus where possible, the device should be designed to sustain this response under a given excitation. Small values of the angle β reduce the asymmetry of the device, thereby expanding the range of a 1:1 response in a shorter device. The value of the restitution coefficient, which can be controlled by the thickness of the DE membranes, can also influence the device performance. This may be especially helpful when the transition to 2:1 motion is inevitable, as demonstrated in Fig. 12.
For larger values of β with s (and thus d) decreasing, the PD transition from 1:1 to 1:1/2T is observed, with transitions to 2:1 behavior via PD G or TG occurring for s above that for ghost grazingG 2:1 1:1 , as shown  Fig. 6a, as β decreases, the value of d PD decreases until it is below d G given by G 2:1 1:1 for a range of r . This trade-off between the variation of β, s and r is then captured by the intersection of the analytically provided curves for d PD and d G , indicating parameter combinations that limit the potential for different types of grazing transitions that yield additional impacts per forcing period.

Conclusions
A combined analysis of smooth and non-smooth bifurcations provides new insight into the interplay of smooth and non-smooth bifurcations. The novel approach is developed in the context of a canonical model of an impact pair, a forced capsule in which a ball moves freely between impacts on either end of the capsule. The results provide a systematic approach to study the interplay of different qualitative transitions in TDOF systems, which have received limited attention as compared with SDOF systems.
The integrated analytical approach provides sequences of bifurcations on stable branches as well as complementary ghost bifurcations on unstable or unphysical branches. The foundation of the results is an analytical characterization of impacting solutions and their stability based on the maps between impacts. In the language of non-smooth dynamics, we use compositions of maps based on the smooth dynamics between intersections with switching manifolds. This foundation provides two different bifurcation analyses, one traditionally based on eigenvalues as is used widely for smooth bifurcations, and a second based on the geometry of the intersection of trajectories with the impacting surface, yielding grazing for zero-velocity impacts. Combining these two approaches, we obtain sequences of period doubling and fold bifurcations together with grazing bifurcations, the type of non-smooth bifurcation naturally occurring in impacting systems. With this analytical basis, we are also able to identify ghost bifurcations for both smooth and non-smooth bifurcations which can occur on unstable or unphysical solutions.
Even though they cannot be observed numerically or experimentally, ghost bifurcations influence the birth or quenching of complex behaviors and additional grazing transitions. For example, the ghost grazing bifurcation provides a lower bound on a variety of transitions in which an additional impact is added per forcing period. Tracking stable and unstable-sometimes termed "invisible"-branches allows the identification of fold bifurcations, which are critical mechanisms in multiple types of grazing transitions. The proximity of ghost period doubling bifurcations to grazing bifurcations can yield different transition scenarios with complex quasi-periodic behavior, where a small variation in the parameters can replace the grazing with period doubling in the observed dynamics. Transient grazing, which can occur following sequences of period doubling that culminate in chaotic windows, can eventually lead to a quenching of the complex behavior when it reaches a branch from a grazing bifurcation, typically at a fold point. All of these influences, indicated by the analytical results, are confirmed by comparisons with numerical results.
While the analysis is developed with broader TDOF systems in mind, we demonstrate that it is also relevant for applications. It is applied to a model of inclined vibro-impact energy harvesting (VI-EH) device, where the energy is generated by the variable capacitance principle via impacts of the ball with a dielectric polymer on the capsule ends. The desire for new sustainable and renewable power technologies as alternatives to conventional batteries has generated strong interest in energy harvesting, usually characterized by energy scavenging from vibrations occurring in manmade machines and structures. As is typical for an energy harvesting system, the VI-EH system considered here includes both a mechanical system to absorb the energy of vibrations, and a transduction mechanism to convert it to electrical energy. The saturation of ideas in the linear theory of energy harvesting, i.e., using linear mechanical and electrical systems, has shifted the focus of scientists to nonlinear systems, exploiting their relatively broad bandwidth that allows efficient energy absorption over a wide range of the excitation frequency. VI systems, as a class of strongly nonlinear system, can provide this preferred wide bandwidth once the mass is engaged in the vibro-impact regime.
The competition between the different bifurcations and their ghosts influences the parameter ranges for favorable energy output of the VI-EH. The analyses of the bifurcation sequences provide their parametric dependencies, indicating the types of transitions that are more likely for given parameter combinations. For small restitution coefficients r and small inclined angle β, first-order grazing bifurcations dominate the lowerorder transitions, i.e., with a small number of impacts per forcing period being increased by one with a small relative impact velocity. Then we observe jumps in the energy output, increasing for output averaged over time, and decreasing for output averaged per impact. In contrast, for larger values of r and β, sequences of period doubling bifurcations are more prevalent, which can result in aperiodic or chaotic dynamics. In that case the change in energy output is more gradual during these sequences, but there are still opportunities for complex grazing-type transitions which can generate additional impacts per period with small relative impact velocities. Then the ghost grazing bifurcation gives a lower bound on these transitions, grazing period doubled solutions PD G and transient grazing TG, at which there are abrupt changes in energy output. These transitions are often affiliated with fold bifurcations, near which there is the potential for bi-stability of behaviors with different energy outputs. By capturing the different sequences of various types of bifurcations, the analysis provides direction on how to avoid undesirable behavior that may occur via grazing or period doubling.
In addition to the valuable insights provided by this analysis, it also indicates several areas of future study that are needed, related to bi-stability, global stability analyses, and aperiodic behavior. Furthermore, given the important role of the restitution coefficient, the results point to the need for more realistic models that take into account interactions of the ball with the membrane.