Modal synchronization of coupled bistable Van der Pol oscillators
1 Modal synchronization of coupled bistable Van der Pol oscillators.
I.B. Shiroky and O.V. Gendelman* Faculty of Mechanical Engineering, Technion – Israel Institute of Technology, Haifa, 3200003, Israel * - contacting author, [email protected]
The paper revisits recently revealed regimes of the "nonconventional synchronization" in systems of coupled bi-stable Van der Pol oscillators. These regimes are characterized by periodic (or quasiperiodic) almost complete energy exchanges between the coupled oscillators. In the paper it is demonstrated that such responses correspond to synchronization of the modulation amplitudes between symmetric and antisymmetric modes of the system, with persistent phase drift. This observation substantially simplifies the treatment, reduces the dynamics to simple phase cylinder and allows to reveal a rich set of global and local bifurcations of limit cycles and tori in the system. Among other findings, one encounters an unexpected regime of nonstationary asymmetric synchronization. Keywords: Limit Cycle; Coupled Oscillators; Synchronization; Modal Space; Global Bifurcation Introduction
Van der Pol (VDP) oscillator entered scientific literature in 1920 [1], and became a paradigmatic model for systems with the limit cycle (LC) as a single attractor. In the case of biophysics, the model has been used to study heart beats [2], neurobiology of the lamprey [3], Parkinsonian tremor [4] , eye chemistry [5, 6], EEG dynamics [7], vocal fold oscillator during phonation [8, 9], action potentials of neurons [10] and others. Multiple additional applications of the model include as different subjects as radio electronics, seismology, nonlinear optics and many others. Modified versions of the VDP oscillator can include higher orders of the nonlinear damping. In particular, it is easy to formulate a bi-stable variation of the VDP oscillator with two attractors: a stable fixed point and a stable limit cycle separated by an unstable limit cycle [11]. Discrete chains of bi-stable van der Pol oscillators were studied in [12, 11] where basic regimes of localization and front propagation were demonstrated. . Phenomenon of synchronization in systems of coupled oscillators is well-known in scientific literature [13, 3, 14, 15, 16]. In most cases, one observes a synchronization of stationary responses of the oscillators. As for the coupled VDP oscillators, studies [3, 14] and many others revealed the synchronization regimes when the excitation amplitude of each oscillator remained constant. The case of strongly coupled VDP oscillators is addressed in [17] and compared to the weakly coupled case. This type of synchronization is classified in resent works as being close to nonlinear normal modes [15, 18]. Recently, lot of attention has been paid to nonstationary regimes with strong energy exchange in coupled oscillatory systems [19, 20, 21]. In particular, recently it was revealed [15, 16, 22] that the system of linearly coupled bistable VDP oscillators exhibits quite unusual pattern of nonstationary "non-conventional" synchronization. In this regime the response of the oscillators is strongly modulated; one can say that they exchange excitation. This regime cannot be described in terms of classical synchronization that primarily deals with the stationary responses. The analytical approach used in [15] was based on special symmetry revealed by the system of slow-flow equations for one special combination of the system parameters. In this case, the slow-flow system was reduced to the flow on the phase plane, and the regime of the "non-conventional" synchronization has been clearly identified as the limit cycle of the slow flow. Further numerical studies demonstrated that this regime robustly appeared for finite region in the parametric space, and not only for the special combination of parameters that endowed the slow-flow equations with the special symmetry. Paper [16] explored the effect of detuning in a similar model. In [22] additional effects of the damping and nonlinearity were addressed. The authors showed that for same sets of parameters two regimes can co-exist – the stationary synchronous dynamics and the "non-conventional" synchronization.
Current paper offers an alternative approach to the same model of two weakly coupled bi-stable VDP oscillators. It will be demonstrated that that these special regimes of the "non-conventional" synchronization correspond to synchronization of the modulation amplitudes between modal coordinates of the system, with persistent phase drift. We are going to demonstrate that this observation drastically simplifies the analysis and allows exhaustive parametric study of the system without reliance on any nontrivial and parameter-dependent symmetry. Then, it is possible to reveal rich bifurcation patterns, as well as the regions in the space of parameters where each of the qualitatively different synchronization regimes is stable. The structure of the paper is as follows. In section 2 the model is described, and the regime of "non-conventional" synchronization is presented in physical variables, i.e. in terms of displacements of the individual oscillators. In section 3 the system is analyzed in the modal space and different response regimes and their bifurcations are described. In section 4 the results are verified numerically versus the original model. Section 5 is devoted to conclusions and discussion. Analysis of the model in physical coordinates
We consider a model of two linearly coupled identical bi-stable VDP oscillators similar to models previously discussed in [11, 12, 15, 16, 22]: ( ) ( ) ( ) ( ) y y c y y y y yy y c y y y y y + + − + − + = + + − + − + = (1) Here c is the coupling between the oscillators, mass, linear frequency and the coefficient of quadratic term in the nonlinear damping are set to unity without loss of generality. Both the coupling and the nonlinear damping are assumed as small, and is the appropriate parameter. First, we numerically demonstrate the regime of strong modulation of the fast oscillations with almost complete energy exchange between the two oscillators (Figure 1). The system was deliberately simulated for a large value of = in order to be able to identify visually the fast oscillations in the plot. 3 Figure 1 – The regime of "non-conventional" synchronization, y - blue, y - red (time rescaled by 20), Parameters: c = = = . The following slow-flow analysis is similar to the one used in [15, 16, 22] and many other works, and is presented here for the sake of completeness. In order to separate the slow evolution of the envelope from the fast oscillations through common analytical techniques [19, 20, 21], complex variables , are introduced: * * *1 1 1 1 1 11 1 1 1* * *2 2 2 2 2 22 2 2 2 , ,2 2 2, ,2 2 2 y i y y iy i y y i − + + = − = = − − + + = − = = − (2) Substitution of (2) into (1) yields: ( )( ) i ci i ci + − − − − − − + + + + = + − − − + − − + + + + = (3) By substitution of , it it e e = = into (3) and by dividing the equation by i t e we achieve: ( ) ( ) ( ) ( )( ) ( ) * 2 * 21 1 2 1 2 2 2 42 2 4 22 2 * 2 4 4 2 2 * 2 * 4* 2 1 1 1 1 1 1 1 1 1 11 1 * 2 * 22 1 2 1 2 222 2 ** 2 2 2 22 2 it itit it it it it itit it ititit i c e ee e e e e eei c e ee ee − − − − −− − − −− − − − + + − − − + − ++ + + + = + − − + +− −++ + ( ) ( ) it it it it it e e e e − − − + − + + = (4) The fast oscillations are averaged out from (4) [23, 12], and one obtains the following set of equations for the slow flow: ( )( ) i ci c − − + − + = − − + − + = (5) Let us define a new time scale such that: d dt dt d = → = = = (6) Substitution of (6) into (5) leads to elimination of from the equations: ( )( ) icic − − + − + = − − + − + = (7) Equations (7) are then simplified through standard polar transformations: , , i i R e R e = = = − (8) Here , , ,
R R are real functions. Substitution of (8) into (5) and rearrangement leads to the following set of equations: ( ) sin 12 2 4 8sin 12 2 4 8cos2 R R RcR R R R RcR R R Rc R R = − − − + = − − + −= − (9) The inherent problem in representation (9) is the singularity of when R = or R = . A typical behavior of , , R R in the regime of non-conventional synchronization is presented in Figure 2. One can observe that the phase difference undergoes periodic jumps from / 2 to / 2 − when R or R approach 0, just due to this singularity. Figure 2 – Slow flow of the strongly modulated response (model (9)), R - blue, R - red, - yellow, Parameters: c = = . Analysis in modal coordinates
As an alternative approach to analysis of System (1), we suggest a transformation to modal coordinates. u y yv y y = += − (10) Here ( ), ( ) u t v t represent the symmetric and the antisymmetric modes of linear counterpart in System (1), respectively. The coefficients in (10) are somewhat uncommonly set to unity and the transformation is not area-preserving, but this particularity is immaterial for the following analysis, since we look for the limit cycles and fixed points of the slow flow. In Figure 3 the response regime presented in Figure 1 is re-plotted in the modal coordinates in accordance with (10). One observes that the transition into modal space leads to a major simplification with respect to the physical variables. Namely, the regime corresponds to synchronization of the modulation amplitudes for both modes. In addition, it is easy to see that the relative phase between the modal coordinates is not constant and exhibits the drift. We note in passing that for the synchronization regimes with constant amplitude, the slow-flow modal coordinates will exhibit steady-state patterns with constant amplitudes and constant relative phase.
Figure 3 – Strongly modulated response in modal space, u - blue, v - red (time rescaled by 20), Parameters: c = = = . To explore the regime of the synchronization between the modal modulation amplitudes, we perform the averaging procedure for System (1) in the modal space and introduce new complex variables , : * * *1 1 1 1 1 11* * *2 2 2 2 2 22 , ,2 2 2, ,2 2 2 u i u u iv i v v i − + + = − = = − − + + = − = = − (11) By taking , it it e e = = , substitution of (10), (11) into (1) and averaging out the fast oscillation, the following set of equations for the slow flow of the modal coordinates is obtained: ( )( )( )( ) ic + + − − − ++ + + + + + =− + + − − − ++ + + + + + = (12) Then, we use transformation , , i i N e N e = = = − and obtain: 7 ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) z zz z z z z z z z z z z zz zz z z z z z z z z z z zc z z z z + − + + + + + + + =+ − + + + + + + + = + + + − + = (13) Here , z N z N = = . The structure of system (13) is symmetrical with respect to z and z . Therefore, the amplitude synchronization ansatz z z = always satisfies this system of equations. These equations of the slow flow in modal space do not involve any singularities, contrary to the slow flow in physical coordinates (9). In Figure 4 we re-plot the regime with complete energy exchange in terms of the slow flow in the modal space (13). One observes the synchronization between , z z and the smooth behavior of . Figure 4 – Slow flow of the regime with energy exchange in modal space, , z z - blue, - yellow, Parameters: c = = . Interestingly, fixed points of System (13) can be analyzed in exhaustive way despite high powers of the involved polynomials. From equations (13) one obtains the following algebraic conditions for the fixed points: ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) z zz z z z z z z z z z zz zz z z z z z z z z z zc z z − + + + + + + + =− + + + + + + + =+ + − ( ) z z + = (14) Cancelling , z z from the first and the second equations respectively, and subtracting the remaining equations, one obtains: ( ) ( ) ( ) z z z z − + − + = (15) It is easy to see that the options ( ) cos 2 1 = − or z z + = are inconsistent with the third equation of System (14). Then, for the fixed points only three combinations of solutions are possible:
1) , 2) 0, 3) 0 z z z z = = = . For the modes z = or z = we observed only trivial dynamics of localization on symmetric or antisymmetric mode; these regimes will stay beyond our scope. Moreover, only the dynamics of amplitude synchronization z z z = will be considered. The stable solutions that will be obtained under this assumption need to be validated in the original modal system (1) for global stability. Dynamics of the amplitude synchronization is therefore described by the following reduced model: ( ) ( ) ( ) ( ) ( ) ( ) z zz zc z z + − + + + = + + − = (16) This system defines the phase flow on cylinder z and has two control parameters ( , c ). The map of regimes at the plane of these control parameters is presented in Figure 4. 9 Figure 5 – Bifurcations and solution zones of system (16) for the case z z = In Figure 5 various regimes and their boundaries are described. First, we address non-trivial fixed points of System (16), described by the algebraic equations: ( ) ( ) ( ) ( ) ( ) z zzc z z − + + + = + − = (17) These fixed points correspond to common synchronized limit cycle oscillations (LCO) with constant amplitudes in the original system, since the phase difference between the modes is also stationary. Simple algebraic manipulations yield:
8. Weakly modulated
SN of limit cycles
2) 1 Stable fixed pt.
3) 1 Stable fixed pt.
1) 2 Stable fixed pts.
9. Trivial Solution 10. Trivial Solution 6) Stable limit cycle
7) Stable limit cycle
5) 1 Stable fixed pt.
Detail A
Detail A z zcz z z zcz z − − ++ = − − = − (18) Let us examine the following polynomial obtained from (18): ( ) ( ) ( ) ( ) ( )
P z z z z z z z c z z = − − + − − + − (19)
By taking ( ) ( )
0, 0 dP zP z dz = = one obtains four different saddle-node (SN) limits, at which two fixed points collide and disappear. The SN limits are plotted in Figure 5. The next issue to address is the stability of the fixed points. For this sake, we rewrite system (13) in the following way: ( ) ( ) ( ) ( ) ( ) ( )( ) ( ) ( ) ( ) ( ) ( )( ) ( ) ( ) ( ) z zz z z z z z z z z z z z f z zz zz z z z z z z z z z z z f z zc z z z z f z z = − + + + − + + + + = − + + + − + + + + = − − + − + (20) The Jacobian for the explored case of z z z = = at the fixed points is given by: z z z z z z z z f f fz zf f fJ J z zf f fz z = = = = = = = = (21)
The eigenvalues for the Hopf bifurcation are purely imaginary i = . Therefore, the following two equations can be written at this limit: Re 0Im 0
J i IJ i I − =− = (22)
We solve numerically the set of 4 equations (22), (18) and get the boundary of the Hopf bifurcation illustrated in Figure 5. Regimes of amplitude synchronization with phase drift correspond to limit cycles of System (16). Due to the cylindrical structure of the phase space, one can encounter a LCO of the first type (not winding the phase cylinder) and of the second type (winding the cylinder once) [24, 25]. These LCO are created through bifurcational mechanisms depicted in Figure 5. Namely, the LCO of the first type appears 11 through Hopf bifurcation and disappears through the homoclinic global bifurcation. LCO of the second type appears due to the saddle-node bifurcation of limit cycles. The dashed line in Figure 5 at = and
1/ 4 c corresponds to the unique case of the special symmetry explored in [15], when the LCO can be described at the phase plane even in physical coordinates. On this line the value of z is constant and equal to 12 while is time dependent. The solution for this case is given in (23). Interestingly, despite the additional symmetry, no bifurcation is encountered at this line besides the terminal point c=1/4. However, we observe in Figure 5 that the SN bifurcation of the LCO of the second type occurs very close to this line. This proximity requires further analysis and understanding. ( )tan 16 1 16 1 1412, tan 4 t t c cz c − − − − + = = − (23) One can see that the system exhibits a rich set of responses and bifurcations. The parametric c − space in Figure 5 is split into 10 different regions with the following characteristics: - Nontrivial stable fixed points – regions 1-5 (shaded with green). In these regions one or two fixed points are stable alongside the always stable trivial solution. -
Limit cycle of the second type – regions 6-7 (shaded with red). In these regions a limit cycle is stable alongside with the trivial solution. -
Limit cycle of the first kind – region 8 (shaded with gray). This narrow region is bounded by lines of the Hopf bifurcation and of the global homoclinic bifurcation. This limit cycle is also stable alongside with the trivial solution. -
Trivial solution – regions 9-10 (shaded with blue). This list represents the possible regimes of amplitude synchronization in the modal space of the system. Successive transitions from region 5 to region 9 are illustrated in Figure 6. In the current representative example, parameter c = is kept constant and parameter decreases. In region 5 one encounters a stable focus and a saddle. Transition into region 8 occurs through a Hopf bifurcation, at which the stable focus becomes unstable and a limit cycle of the first kind, marked with green, is formed. Then, a global homoclinic bifurcation is reached and the limit cycle merges with the stable and unstable manifolds of the saddle and disappears. Further decrease of leads into region 9 where one observes an unstable focus and a saddle. All initial conditions on the cylinder lead to the lone stable attractor of the amplitude death. 2 Figure 6 – Transition from region 5 into region 9 through Hopf and homoclinic bifurcations. Parameters: Region 5: c = = , Region 8: c = = , Homoclinic: c = = , Region 9: c = = Then, we describe the creation of the limit cycle of the second kind. At the upper limit, the limit cycle is created through a SN bifurcation of limit cycles in transition between regions 9 and 10 to regions 6 and 7 respectively. The transition from region 10 to 7 is illustrated in Figure 7.
Region 5
Region 8
Homoclinic
Region 9
Region 10
Region 7 Figure 7 – Transition from region 10 into region 7 through a saddle-node of limit cycles bifurcation. Parameters: Region 10: c = = , Region 7: c = = . On the other limits, from the left and from below, the limit cycle of the second kind (marked with green) is created through SN bifurcation – the saddle and the node disappear and form the cycle. This case is demonstrated in Figure 8. Figure 8 – Transition from region 1 into region 6 through saddle-node bifurcations on the limit cycle. Parameters: Region 1: c = = , Region 2: c = = , Region 4: c = = , Region 6: c = = . Numerical verification of the results
So far, the qualitative analysis of system (16) seems complete. In order to reach system (16), two major simplifications were conducted: averaging and the reduction to the special case of amplitude synchronization in modal space z z z = = . The next step is to validate these assumptions, by examining the stability of the unique attractors of the system in the original system ((1) and (10)) and in the full averaged system (13). For this goal we perform numerical tests in the interesting section c = where most of attractors can be found within minor variations of parameter . The first attractor to examine is the limit cycle of the second kind response that is expected according to the analysis, for instance, at = . This case is shown for the original coordinates , y y (according to Equation (1)) in Figure 9a and in modal space Region 2
Region 4
Region 6
Region 1
SN4
SN3
SN4
SN3
15 (10) for coordinate u alongside the envelope obtained from the full averaged model (13) in Figure 9b. The averaged model follows exactly the envelope of the response. We again observe that extreme energy exchanges between the physical oscillators correspond to only relatively weak modulation in the modal space. The second attractor, expected for a representative value = , is a limit cycle of the first kind, contained in a narrow region between the Hopf and homoclinic bifurcations in Figure 5. As expected, we find such a response (Figure 10) for both the modal system and the full averaged system. Yet again a perfect match is obtained between the averaged and the original modal system. Generally, the provided numerical verifications imply that the analysis of the simplified system (16) as outlined in Figure 5 faithfully represents the features of the original system. Figure 9 – Limit cycle of the second kind of the slow flow, (a) Original coordinates, y - blue, y - red (time recalled by 1000), (b) Modal coordinates, full response – blue (time recalled by 1000), averaged envelope (13) – orange, Parameters: c = = = . Figure 10 – Limit cycle of the first kind of the slow flow, (a) Original coordinates (time recalled by 1000), (b) Modal coordinates, full response – blue (time recalled by 1000), averaged envelope (13) – orange, Parameters: c = = = . (a) (b) (a) (b)
6 The results in Figure 9 for the LC of the second kind are similar to large extent to the situation described in Section 2. The results for the LC of the first kind in Figure 10 are somewhat surprising. One can observe that weakly modulated amplitude-synchronized LCO in the modal space (Figure 10b) correspond to asymmetric torus in the physical space. Equations (1) are symmetric, and therefore substitution y y will yield different solution in the physical space with the same slow-flow envelope in the modal space. Mathematically it is possible since the transformation from 4-dimensional state space of the original system to 3-dimensional slow-flow equations in injective, but obviously not bijective. Concluding remarks . The results presented above reveal the underlying mechanism of the "non-conventional" synchronization in a system of weakly coupled bi-stable van der Pol oscillators. This regime, in fact is a LCO in the slow-flow space with synchronization of the modulation amplitudes between the symmetric and antisymmetric modes of the system. Averaging in the modal space substantially simplifies the treatment and removes the singularities from the slow-flow equations. Moreover, these slow-flow equations in the modal space possess sufficient internal symmetry to reduce the dynamics to a simple phase cylinder for any set of the control parameters. Therefore, it is possible to investigate the structure of the parametric space and to explore the set of response scenarios and their local and global bifurcations. Previously observed “nonconventional synchronization” with almost complete energy exchange between the oscillators correspond to the LCO of the second kind on the phase cylinder. In addition, one can observe the regime of asymmetric nonstationary synchronization in the physical space depicted in Figure 10, that corresponds to the LC of the first kind for the slow flow in the modal space. To the best of the author’s knowledge, this latter regime has not been observed in previous studies on the subject. One can expect that the observed regimes will appear generically in similar systems of coupled oscillators with self-excitation. As for the considered system, the unresolved issues include the role of special analytically derivable limit cycle for
1/ 9 = . Both global bifurcations (homoclinic and SN of the LCO) occur not on this line, but in its close vicinity. Unusual sensitivity of responses on miniscule variation of the control parameters allows one to expect that further simplifications in this and similar systems are possible. This issue will be the subject of further investigations. Acknowledgment
The authors are very grateful to Israel Science Foundation (grant 1697/17) for financial support.
Conflict of interests
The authors declare no conflict of interests. References [1]
B. van der Pol, "A theory of the amplitude of free and forced triode vibrations,"
Radio Review, vol. 1, p. 701–710, 1920. [2] B. van der Pol and J. van der Mark, "The heartbeat considered as a relaxation oscillation, and an electrical model of the heart,"
Phil. Mag., vol. 6, p. 763–775, 1928. [3] R. H. Rand and P. J. Holmes, "Bifurcation of periodic motions in two weakly coupled van der Pol oscillators,"
Int J. Non-Linear Mechanics, vol. 15, pp. 387-399, 1980. [4] A. Beuter, R. Edwards and M. Titcombe, "Data analysis and mathematical modeling of human tremor," in
Nonlinear Dynamics in Physiology and Medicine , Springer, 2003. [5] K. Rompala, R. Rand and H. Howland, "Dynamics of three coupled van der Pol oscillators with application to circadian rhythms,"
Commun. Nonlinear Sci. Numer. Simul., vol. 12, pp. 794-803, 2007. [6] E. Camacho, R. Rand and H. Howland, "Dynamics of two van der Pol oscillators coupled via a bath,"
International Journal of Solids and Structures, vol. 41, p. 2133–2143, 2004. [7] K. Kirschfeld, "The physical basis of alpha waves in the electroencephalogram and the origin of the ‘Berger effect”,"
Biol. Cybern., vol. 92, p. 177–185, 2005. [8] R. Garrel, R. Scherer, R. Nicollas, A. Giovanni and M. Ouaknine, "Using the relaxation oscillations principle for simple phonation modeling,"
J. Voice, vol. 22, p. 385–398, 2008. [9] J. C. Lucero, L. L. Koenig, K. G. Lourenco, N. Ruty and X. Pelorson, "A lumped mucosal wave model of the vocal folds revisited: recent extensions and oscillation hysteresis,"
J. Acoust. Soc. Am., vol. 129, p. 1568–1579, 2011. [10] R. FitzHugh, "Impulses and physiological states in theoretical models of nerve membranes,"
Biophysics J, vol. 1, p. 445–466, 1961. [11] Y. D. Defontaines and Y. Pomeau, "Chain of coupled bistable oscillators: a model,"
Physica D, vol. 46, p. 201–216, 1990. [12] I. B. Shiroky, A. Papangelo, N. Hoffmann and O. V. Gendelman, "Nucleation and propagation of excitation fronts in self-excited systems,"
Phys. D, vol. 401, p. 132176, 2020. [13] A. Pikovsky, M. Rosenblum and J. Kurths, Synchronization: A universal concept in nonlinear sciences, Cambridge University Press, 2001. [14] T. Chakraborty and R. H. Rand, "The transition from phase locking to drift in a system of two weakly coupled van der Pol oscillators,"
J. Nonlinear Mech., vol. 23, pp. 369-376, 1988. 8 [15] L. I. Manevitch, A. Kovaleva and V. N. Pilipchuk, "Non-conventional synchronization of weakly coupled active oscillators,"
EPL, vol. 101, p. 50002, 2013. [16] M. Kovaleva, V. Pilipchuk and L. Manevitch, "Nonconventional synchronization and energy localization in weakly coupled autogenerators,"
Phys. Rev. E, vol. 94, p. 032223, 2016. [17] D. W. Storti and R. H. Rand, "Dynamics of two strongly coupled Van der Pal oscillators,"
Int. J. Non-linear Mech., vol. 17, pp. 143-152, 1982. [18] A. F. Vakakis, L. I. Manevitch, Y. Mikhlin, V. N. Pilipchuk and A. A. Zevin, Normal Modes and Localization in Nonlinear Systems, New York: Wiley, 1996. [19] L. I. Manevitch, "New approach to beating phenomenon in coupled nonlinear oscillatory chains,"
Archive of Applied Mechanics, vol. 77, p. 301–312, 2007. [20] L. I. Manevitch and V. V. Smirnov, "Limiting phase trajectories and the origin of energy localization in nonlinear oscillatory chains,"
Phys. Rev. E, vol. 82, p. 036602, 2010. [21] L. I. Manevitch and O. V. Gendelman, Tractable Models of Solid Mechanics: Formulation, New York: Springer, 2011. [22] M. Kovaleva, L. Manevitch and V. Pilipchuk, "Non-conventional phase attractors and repellers in weakly coupled autogenerators with hard excitation,"
EPL, vol. 120, p. 30007, 2017. [23] O. V. Gendelman, L. I. Manevitch, A. F. Vakakis and R. M'Closkey, "Energy Pumping in Nonlinear Mechanical Oscillators: Part I—Dynamics of the Underlying Hamiltonian Systems,"
J. Appl. Mech, vol. 68, no. 1, pp. 34-41, 2000. [24] A. A. Andronov, A. A. Vitt and S. E. Khaikin, Teoriya kolebanii (The Theory of Oscillations), Moscow: Nauka, 1981. [25] L. A. Cherkas and A. A. Grin, "Limit cycle function of the second kind for autonomous systems on the cylinder,"