Physics Letters A 379 (2015) 2830–2838 Contents lists available at ScienceDirect Physics Letters A www.elsevier.com/locate/pla Crises in a dissipative bouncing ball model André L.P. Livorati a,b,c,∗, Iberê L. Caldas c, Carl P. Dettmann b, Edson D. Leonel a a Departamento de Física, UNESP, Universidade Estadual Paulista, Av. 24A, 1515, Bela Vista, 13506-900, Rio Claro, SP, Brazil b School of Mathematics, University of Bristol, Bristol, BS8 1TW, United Kingdom c Instituto de Física, IFUSP, Universidade de São Paulo, USP, Rua do Matão, Tr.R 187, Cidade Universitária, 05314-970, São Paulo, SP, Brazil a r t i c l e i n f o a b s t r a c t Article history: Received 16 July 2015 Received in revised form 9 September 2015 Accepted 10 September 2015 Available online 14 September 2015 Communicated by C.R. Doering The dynamics of a bouncing ball model under the influence of dissipation is investigated by using a two-dimensional nonlinear mapping. When high dissipation is considered, the dynamics evolves to different attractors. The evolution of the basins of the attracting fixed points is characterized, as we vary the control parameters. Crises between the attractors and their boundaries are observed. We found that the multiple attractors are intertwined, and when the boundary crisis between their stable and unstable manifolds occurs, it creates a successive mechanism of destruction for all attractors originated by the sinks. Also, a physical impact crisis is described, an important mechanism in the reduction of the number of attractors. © 2015 Elsevier B.V. All rights reserved. 1. Introduction Modeling of dynamical systems is one of the most embracing areas of interest among physicists and mathematicians in general [1]. Very popular among these models are low-dimensional sys- tems [2,3], whose complex dynamics leading to a rich variety of nonlinear phenomena [3–6], including bifurcations in non-smooth dynamical systems [7]. Here we study the problem of a bouncing ball model, where a free particle is suffering collisions with a vibrating wall under the presence of a constant gravitational field. Holmes [8,9] and Pustylnikov [10,11] were among the first to study the bouncing ball dynamics. This model has been used in many physical and engineering applications. For instance, it describes a similar accel- eration phenomenon that cosmic rays experience to acquire high energies, known as Fermi acceleration [12] (considered as the first attempt of prototype for the bouncing ball dynamics); the dynamic stability in human performance, where a human tries to stabilize a ball on a vibrating tennis racket [13]; and the subharmonic vibra- tion waves in a nanometer-sized mechanical contact system [14]. One can also find studies in granular materials [15–18], experi- mental devices concerning normal coefficient of restitution [19, 20], mechanical vibrations [21–23], anomalous transport and diffu- sion [24,25], thermodynamics [26], chaos control [27–29], besides the well known connection with the standard mapping [2], which leads to other several applications. * Corresponding author. E-mail address: livorati@usp.br (A.L.P. Livorati). http://dx.doi.org/10.1016/j.physleta.2015.09.016 0375-9601/© 2015 Elsevier B.V. All rights reserved. Although the bouncing ball problem has been studied for many years [8–11,30,31], concerning different aspects and applications, the implications of the nonlinear perturbation requires an exten- sive and complex analysis where some chaotic properties are not yet fully understood. In this paper we consider a high dissipative bouncing ball model where a coefficient of restitution plays the role of dissipation, and the perturbation parameter is physically in- terpreted as a ratio between the moving plate acceleration and the gravitational field. For some combinations of parameters, plenty of attractors can coexist [32–34]. We found that these attractors in the phase space are intertwined, and varying the value of the con- trol parameter of perturbation, we characterize a boundary crisis [6,35–37] between the stable and unstable manifold of the same saddle point. Such a crisis leads to successive destruction of these intertwined attractors and is a mechanism that allows the lowest energy attractor, which is related to the vibrating wall, to continue to exist, giving it the status of a robust attractor. In addition, we describe a physical impact crisis, between the real vibrating plate and the border of an attractor. This crisis, as yet unclassified, re- duces the number of attractors dramatically at a single parameter value. The organization of the paper is given as follows. In Section 2 we describe the dynamical system under study and its chaotic properties. Section 3.1 is devoted to the numerical analysis of the average velocity, in Section 3.2 we study the basin of attraction of the fixed points and set up the impact physical crisis, and in Sec- tion 3.3 we discuss the relation between the manifolds boundary crisis and the attractors; finally in Section 4 we draw some final remarks and conclusions. http://dx.doi.org/10.1016/j.physleta.2015.09.016 http://www.ScienceDirect.com/ http://www.elsevier.com/locate/pla mailto:livorati@usp.br http://dx.doi.org/10.1016/j.physleta.2015.09.016 http://crossmark.crossref.org/dialog/?doi=10.1016/j.physleta.2015.09.016&domain=pdf A.L.P. Livorati et al. / Physics Letters A 379 (2015) 2830–2838 2831 2. The model, the mapping and chaotic properties In this section we describe the model under study, the bounc- ing model, which consists of a particle, under the influence of a constant gravitational field, that suffers inelastic collisions with a heavy oscillating wall. Dissipation is introduced via a restitution coefficient γ ∈ [0, 1], where γ = 1 recovers the conservative case, where Fermi Acceleration (FA) is inherent [25,38]. The introduction of dissipation can be considered as a suppression mechanism for this unlimited energy growth [39,40]. The system is oriented along the vertical axis, where the upward direction is said to be positive, the wall equilibrium position is set at y = 0, and the dynamics is basically described by a non-linear mapping for the variables ve- locity of the particle v and time t immediately after a nth collision of the particle with the vibrating wall. There are two distinct versions of the dynamics description: (i) complete one, which consists in considering the complete move- ment of the time-dependent wall, and (ii) simplified, where the wall is assumed to be fixed, but exchanges momentum and en- ergy with the particle upon collision. Both approaches produce a very similar dynamic considering conservative [25] and dissipa- tive cases [26,39–41]. In the complete version, the vibrating wall obeys the equation yw(tn) = ε cos wtn , where ε and w are re- spectively, the amplitude and the frequency of oscillation of the vibrating wall. In the simplified version, the vibrating wall is said to be fixed at y = 0, but when the particle collides with it, they ex- change momentum and energy as if the wall were vibrating. Thus, the simplified approach keeps the nonlinearity of the model and significantly speeds up the numerical simulations, as well allows easier analytical calculations. In this paper and from this point be- yond, we only deal with the complete version of the mapping. Considering the flight time, which is the time that the particle spends to go up, stop with zero velocity, starts falling and collides again with the vibrating wall, we define some dimensionless and more convenient variables as: Vn = vn w/g , ε = εw2/g , where Vn is the “new dimensionless velocity”, g is the gravitational field and ε can be understood as a ratio between accelerations of the vi- brating wall and the gravitational field. For instance, one can set some real values for the dimensional variables, as g = 10 m/s2, ε = 0.001 m, w = 2π f , where f = 100 Hz, and obtain the dimen- sionless variable ε ≈ 0.1591. Some real devices concerning impact experiments with granular material can be found in Refs. [19,20]. Also, measuring the time in terms of the number of oscillations of the vibrating wall φn = wtn , we obtain the mapping T : { Vn+1 = −γ (V ∗ n − φc) − (1 + γ )ε sin(φn+1), φn+1 = [φn + �Tn] mod (2π), (1) where the expressions for V ∗ n and �Tn depend on the kind of the considered collision. For the case of multiple collisions inside the collision zone [−ε, +ε], the expressions are V ∗ n = Vn and �Tn = φc where φc is obtained from the condition that matches the same position for the particle and the vibrating wall, expressed as G(φc) = ε cos(φn + φc) − ε cos(φn) − Vnφc + 1 2 φ2 c , (2) where this transcendental equation must be solved numerically for G(φc) = 0, with φc ∈ (0, 2π ]. If the particle leaves the collision zone case after a colli- sion, goes up, reach null velocity, and falls for an another col- lision, we have indirect collisions and the expressions are V ∗ n = − √ V 2 n + 2ε(cos(φn) − 1) and �Tn = φu + φd + φc with φu = Vn denoting the time spent by the particle in the upward direction up to reaching the null velocity, φd = √ V 2 n + 2ε(cos(φn) − 1) cor- responds to the time that the particle spends from the place where it had zero velocity up to the entrance of the collision zone at ε . Fig. 1. Comparison between phase space for conservative and dissipative dynamics. In (a) ε = 0.6 and γ = 1.0, and in (b) ε = 0.6 and γ = 0.9. In (b) the thick black regions are the sinks and the bottom attractor. Also, all the spread dots are the transient. Finally the term φc has to be obtained numerically from the equa- tion F (φc) = ε cos(φn + φu + φd + φc) − ε − V ∗ n φc + 1 2 φ2 c , (3) where F (φc) represents a transcendental equation that must be solved numerically in order to find the exact “time” of collision, as F (φc) = 0, with φc ∈ [0, 2π ]. The obtainment of the numerical root φc is done considering at first G(φc) = 0. If we did not find any root for G(φc), we start to evaluate F (φc) = 0. The root seeking process is made by solving the transcendental equations via bisection method, with a preci- sion of 10−14. Taking the determinant of the Jacobian matrix of both kinds of collisions (see Ref. [40] for details), and after a straightforward al- gebra, it is easy to show that the mapping (1) shrinks the phase space measure since the determinant of the Jacobian matrix is given by Det J = γ 2 [ Vn + ε sin(φn) Vn+1 + ε sin(φn+1) ] . (4) Here, if γ = 1 we recover the non-dissipative version of the map- ping, in fact, as velocity and phase are not canonical pairs in the complete version, the determinant of J is not the unity, but rather it leads to the following measure to be preserved, dμ = (V +ε sin φ)dV dφ. Indeed, the extended phase space for the whole version of the model considers four variables namely: (1) yw de- noting the position of the vibrating wall; (2) V p corresponding to the velocity of the particle; (3) E p which is the mechanical en- ergy (kinetic+gravitational) of the particle and (4) the time t . The canonical pairs however are: position and velocity (yw , V p) and energy and time (E p, t). Another useful property for the dynamics evolution, as func- tion of the control parameters, is the analysis of the fixed points and their stability. For the bouncing ball model the period-1 fixed points can be obtained by doing Vn+1 = Vn = V ∗ and φn+1 = φn = φ∗ + 2mπ in Eq. (1). For both kinds of collisions, successive and indirect, the fixed points are V ∗ = mπ ;m = 1,2, . . . , φ∗ = arcsin ( V ∗(γ − 1) ) . (5) (1 + γ )ε 2832 A.L.P. Livorati et al. / Physics Letters A 379 (2015) 2830–2838 Fig. 2. (Color online.) Behavior of the average velocity curves for an extensive range of ε and γ . In (a) and (c), we have a high dissipation and low ε , and the dynamics are basically controlled by the attracting fixed points. In (b) and (d), we have small dissipation and high ε . Here the dynamics is ruled by the attractor localized in the bottom, embedded with the vibrating wall. Also, in (d) a rearrange in the horizontal axis is made, in order to make all the curves start to grow together. This rearrange is convenient, when we make use of the scaling techniques. Their stability are given by Det( J − λI) = 0, evaluated over the fixed points, where λ are the eigenvalues of the Jacobian matrix and I is the identity matrix. The eigenvalues can be found solving the expression λ1,2=Tr J±√ Tr J 2−4 Det J 2 [2]. Fig. 1 shows a comparison between the phase space for the conservative and dissipative versions. As the dynamics evolves, the introduction of dissipation destroys the invariant curves in the sta- bility islands, and the stable fixed points become sinks [1,2,6]. Depending of the control parameter, there may have a plenty of these attractors, where the orbits converge to they. In Fig. 1 one can see that after the dissipation was introduced the first three stability islands, denoted by white regions among the chaotic sea in Fig. 1(a), became attracting fixed points. Also, there is the pres- ence of an attractor on the bottom of Fig. 1(b), near where the vibrating wall is located. For a better understanding and visualiza- tion, we are going to use in all figures the phase representation between −π and +π . 3. Results and discussions In this section, we describe the results obtained by the numer- ical simulations. First we draw some average velocity curves for a combination of the two control parameters, ε and γ . Then, the basins of attraction of some attracting fixed points are obtained. We investigate how these basins of attraction behave, as we range the control parameters. We constructed some bifurcation diagrams for the fixed points and analyzed how their stability vary. By draw- ing the unstable and stable manifolds we notice that the attractors are intertwined in the phase space, and there is a boundary crises, a crossing of their manifolds, which creates a successive mecha- nism of destruction for all attractors originated by the sinks. Also, we made a study over the bifurcation process and the chaotic properties of it. 3.1. Average velocities Let us set the equation for the average velocity, which depends on both ε and γ . The statistics was made considering two steps, an average taken along the orbit, evolved until a finite high number of collisions n and an average taken along the ensemble of initial conditions. So, we may define V i(n, ε,γ ) = 1 n n∑ j=1 V j , (6) and hence V = 1 M M∑ i=1 V i(n, ε,γ ) , (7) where M represents an ensemble of initial conditions. The index j and i represent respectively the collision number, and the number of initial conditions. Fig. 2 shows the evolution of the average velocity for an exten- sive range of the control parameters ε and γ as function of the number of collisions. Here, we consider two regimes: (i) high dis- sipation and small ε , and (ii) small dissipation and high ε . One can identify very distinct behavior between both regimes. Consid- ering regime (i), as shown in Figs. 2(a, c), we see that the V curves start with high initial velocity near V 0, given as an initial condi- tion, and experience an exponential decay, [41]. After a transient they bend towards different regimes of saturation in low energy levels. Basically, the orbits are attracted by the several fixed points that coexist in the phase space, when the perturbation parame- ter ε is still small enough, therefore marking their convergence to different plateaus. On the other hand, in regime (ii), as present in Fig. 2(b, d), the V curves start with low initial velocity, near V 0 given as an initial condition, and they experience a growth for short time according A.L.P. Livorati et al. / Physics Letters A 379 (2015) 2830–2838 2833 Fig. 3. (Color online.) Basins of attraction for the fixed points considering ε = 1.0 and γ = 0.95. The colors represent the fixed points were the initial conditions are being attracted to. Black → π , red → 2π , green, → 3π , blue → 4π , yellow → 5π and brown → 6π , and the white regions denotes initial conditions that converged to the bottom attractor. to a power law with exponent β ≈ 0.5, until they bend towards a saturation regime, in high energy levels. For this case, the sat- uration plateaus of high energy are connected with the chaotic attractor, and there are no attracting fixed points for such high energy regime. Also, if we rescale the horizontal axis by nε2 all V curves, start to grow together. Fig. 2(b, d) contain previously re- sults observed in [40], for regime (ii). We recommend for a more complete investigation on the regime (ii), concerning the chaotic attractor behavior and a fully scaling analysis on the V curves, in Refs. [39,40] for a numerical point of view; and in Refs. [26,41] for an analytical interpretation. Also, one could check Refs. [42–44] for formal and analytical points of view for statistical properties for similar dynamical systems. However, in this paper, we are interested in the range of high dissipation and small ε (regime (i)), shown in Figs. 2(a, c), where the attracting fixed points still exist and play an important role in the dynamics. Indeed, this high dissipation analysis can be very useful in a further experimental analysis of the bouncing ball model, once they are easily obtained, instead of the tiny dissi- pations as used in Refs. [26,39–41], that would be impossible to obtain in a laboratory. Also, one could think about experience with granular materials interacting with vibrating plates [16–20], as a direct application of the bouncing ball model and its phenomena. 3.2. Basins of attraction and bifurcations Depending on the combination of the control parameters ε and γ , several attracting fixed points can coexist in the phase space, where some of them can be more influential to the dynamics than other. In order to understand how these attracting fixed points be- have as related to initial conditions, Fig. 3 shows the basins for the periodic attractors for ε = 1.0 and γ = 0.95, where a grid of 1000 × 1000 initial conditions were equally split and set in the axis of velocity V ∈ [−ε, 7π ] and phase φ ∈ [−π, +π ]. The differ- ent colors represent the basins of attraction for each attractor of period-1 located in V ∗ = mπ , according the fixed points obtained in Eq. (5). Here, each initial condition were evaluated up to a tran- sient of 105 collisions, and then we evaluate another more 105 collisions, and marked its final velocity in the phase space. Since, we have a positive restitution coefficient, one may think that the particle would “glue” on the vibrating wall for long times. This indeed can happen, if it lands deep enough in the absorbing region of the phase space, the particle will perform a large num- ber of smaller and smaller bounces, that could follow progressive geometric conversions [31]. Such peculiar behavior is known as locking [31]. The white regions in Fig. 3 denotes the initial condi- tions that converged to the locking region attractor, i.e. the bottom attractor embedded with the vibrating wall. However, if the par- ticle has a positive relative velocity, it will not be glued to the wall. And, depending of the control parameters ε and γ , this ve- locity can acquire multiplicity as function of π , as set by Eq. (5), in different regions of the phase space, giving birth to period-1 attrac- tors for higher velocities. In Fig. 3, each color denotes a different period-1 fixed point basin of attraction. As m is increased, it seems that the basins are getting smaller, which could be a possible indi- cation of less influence in the dynamics. Also, for the higher values of m, the boundaries of the basins, that are limited by the unstable manifolds of each fixed point [1,6,9], behave in a very complicated stretching and mixing way, folding themselves embedded like a horseshoe [1,6,9]. The beautiful mazy behavior of the basins of attraction shown in Fig. 3, can drastically change, when ε is increased. As one can see in Fig. 4, where the evolution of the basin of attractions for γ = 0.8 and some values of ε is shown, and a very dynamic scenario between the attractors are set. Each item of Fig. 4 was constructed considering a grid of initial conditions of 500 × 500 equally distributed along the velocity V ∈ [−ε, 3π ] and phase φ ∈ [−π, +π ]. The color scale represents the final velocity of each initial condition after a transient of 105 iterations. The color scale was kept fixed as the final velocity ranges. For low velocities and orbits that were attracted to the bottom attractor, we have the darker colors as black, dark blue (black). The orbits that were at- tracted by the first two sinks, V ∗ = π and V ∗ = 2π , the color scale ranges from blue (dark gray) to yellow (light gray), and for the high velocity attracting orbits (said above V ∗ = 2π ), we have the color range scale between yellow (light gray) and red (dark gray). In Fig. 4, as ε is increased, the boundaries of the basins of attraction start to grow following the stretching and mixing be- havior, just like the Smale horseshoes, [1,6,9], as Figs. 4(a, b, c) dis- plays, where in particular for Fig. 4(a), the sink located in V ∗ = π did not exist yet, the same applies for Fig. 4(b), where the sink located in V ∗ = π also did not exist. One can check that by the period-1 one fixed points expressions, specially for φ∗ , given in Eq. (5). As we increase ε , the sinks located in V ∗ = mπ start to notice the consequences of being lied intertwined with the bottom attractor. Considering the first sink V ∗ = π , it suffers a tangent bifurcation in εt ≈ 0.964 (this bifurcation will be explained in the manifolds section), which creates three new attracting zones inside its own basin of attraction, as shows Fig. 4(d). Raising the value of ε a little bit, the sink in V ∗ = π , bifurcates, creating, together with the other sinks, a plenty of attracting regions in the phase space, as shown in Fig. 4(e). The crisis happens near εc ≈ 1.22635, shown in Fig. 4(f). After a tangent bifurcation in εt , three branches evolves as ε increases. The two upper branches collide with each other, gen- erating a boundary crisis that destroys both branches. At the same parameter εc , the lower branch, suffers a physical collision with the vibrating boundary. This is a different crisis, once the attractor is colliding with a physical structure, instead of a another attractor or manifold. This non-categorized crisis can be better visualized in Fig. 5. Finally, in Figs. 4(g, h), the basin of attraction of the fixed 2834 A.L.P. Livorati et al. / Physics Letters A 379 (2015) 2830–2838 Fig. 4. (Color online.) Characterization of the dynamic scenario between the attractors, by the evolution of the basins of attraction for γ = 0.8 for a grid of 500 × 500 initial conditions. The perturbation parameter follows as: (a) ε = 0.32005, (b) ε = 0.60095, (c) ε = 0.88715, (d) ε ≈ εt = 0.964, (e) ε = 1.08590, (f) ε ≈ εc = 1.22635, (g) ε ≈ εd = 1.25815 and (h) ε = 1.50195. The color scale denotes the final velocity of each pair of initial condition after a transient of 105 iterations. One can see that the basins of attraction suffer huge transformations, as successively bifurcation process as ε is increased, until they are destroyed when ε acquires critical values. Fig. 5. (Color online.) Bifurcation diagram as function of ε for γ = 0.8. It is shown a scenario of the evolution of the attracting fixed points, where attractors originated by the sinks V ∗ = mπ , coexist with period-3 orbits originated by tangent bifurca- tions, the bottom attractor (given by the dashed line) and others attractors. This figure illustrates the plenty of attractors and the high active dynamic structure that orbits may experience. Also, three critical values of ε are set concerning each char- acterized crisis. point V ∗ = π (said together with the branches of the tangent bi- furcation) is totally destroyed, and the other basins of attraction for the other upper sinks started their successive destruction pro- cess. In the end, when we have a large enough value of ε , only the bottom attractor remains in the system. For a better understanding and visualization of the crises, sinks bifurcations and the evolution of the basins of attraction, we con- structed a bifurcation diagram for some fixed points, as shown in Fig. 5, for a fixed dissipation γ = 0.8. The diagram illustrates the plenty of attractors and the high active dynamic structure that or- bits may experience during the time evolution. This diagram was constructed in two different ways: (i) we have an initial condition very near the sinks, and let it follows the attractor, as ε increases; and (ii) we kept the same initial condition for every value of ε . In both cases, the range of ε was split in 5000 equal parts. Consider- ing the evolution in (i), the dynamics follows a regular bifurcation diagram, where the period-1 sinks located in V ∗ = mπ , stay sta- ble until they suffer successive bifurcations as ε increases, until it finally disappear. In the same manner, tangent bifurcations hap- pen for V ∗ = π , V ∗ = 2π and so on. These new branches for each tangent bifurcation, have the same fate of the attractors originated by the sinks. These fixed points of the tangent bifurcation were obtained considering case (ii), where the evolution of the same initial condition, gives rise to different attractors in the bifurcation diagram, for example, the one located near V ∗ ∼= 1.7, basically suf- fers the same bifurcation process as the main ones and then it is destroyed. One can notice in Fig. 5, that there are three critical values: (i) εt ≈ 0.964, which is the value of ε where the tangent bifur- cation occurs for the first sink V ∗ = π ; (ii) εc ≈ 1.22635, which is the critical parameter where the crisis occurs for the three branches originated in the tangent bifurcation. One can see in Fig. 5, that in ε = εc , there are two simultaneous crises. The up- per two branches collide, and destroy each other, while the lower branch physically collides with the vibrating wall, denoted here by the dashed line, characterizing an yet unclassified physical crisis, between the real structure (vibrating wall), and the border of an attractor. Finally, (iii) εd , is the value where the attractor origi- nated by the sink V ∗ = π is destroyed. Beyond εd , a successive destruction mechanism takes place in the dynamics, destroying all the other attractors, except by the bottom one, that starts to rule the dynamics. The values of these critical ε parameters may vary as we range the dissipation. Let us now address to the chaotic attractor born from the bi- furcation process of the fixed point V ∗ = π . As one can see from Fig. 5, when the bifurcation process occurs, we have the creation A.L.P. Livorati et al. / Physics Letters A 379 (2015) 2830–2838 2835 Fig. 6. In (a) we show an amplification of the successive bifurcation process that the first sink suffers as ε is increased, in (b) we made a bigger zoom yet in this process. In (c), (d), (e) and (f), we display the shape of the local chaotic attractors in the phase space, for the lower branch of the first sink. The values of ε used were: in (c) ε = 1.2520, (d) ε = 1.253631, (e) ε = 1.253632, and in (f) ε = 1.2550. One can see that the attractors pass through a bifurcation process, and after a collision between the branches, they merged into a bigger chaotic attractor. In all items we considered γ = 0.8. Fig. 7. Successive zoom-in windows of the chaotic attractor for ε = 1.2549. It seems that the chaotic attractor has a fractal-like shape very similar to one found in the Henon–Heiles map [31]. of two main symmetric branches. According to Fig. 5, each branch, will begin its own bifurcation process, which further will collide with itself in a boundary crisis, generating its own local chaotic at- tractor. Fig. 6 shows how these processes occur for the first sink, presenting the shape of the attractors in the phase space, also we compare it with amplifications of the bifurcation diagram. Once, both branches are symmetric, let us do this analysis considering just one of them. One can see in Figs. 6(a, b) the behavior of bifurcation process of the branches as ε is increased. Comparing them with Fig. 6(c), for ε = 1.2520, we can see that the attractor is under a period-8 bifurcation. Increasing ε , and after a compar- ison between Figs. 6(a, b) and Figs. 6(d, e); the two branches of the attractor merge together into one local chaotic attractor. Here we can characterize another crisis between the attractors, where two attractors become one, known as fusion crisis [1,6]. Finally, in Fig. 6(f) we have the behavior of the local chaotic attractor in its final stage for ε = εd ≈ 1.2550, where it seems to take the shape of the old basin of attraction, originated by the tangent bifurcation, and then it vanishes. Here, another crisis can be characterized. The attractor (darker region), collides with the transient basin of the old attractor already destroyed (originated by the tangent bifurca- tion). This ghost behavior of the transient basin [6], allows to the attractor to interact with the region that suffered the physical colli- sion with the vibrating wall. This interaction, destroys the attractor in the same manner that destroy the previous one in ε ≈ εc . One can ask about the nature of this local chaotic attractor. Indeed, it seems to have fractal shape, as one can see in the suc- cessive amplification of Figs. 7(a, b, c, d). The shape of the chaotic attractor is basically the same, no matter how much further in the zoom windows. It is interesting, the fact that the kind of crises we are seeing here, are basically the same one found in the Henon– Heiles attractor [6,45,46], for which chaotic attractor collides with the stable manifold and is destroyed. Also, the same fractal-like shape of the chaotic attractor is found. It would be interesting to investigate later, if any other chaotic common property related to both attractors can be found. 2836 A.L.P. Livorati et al. / Physics Letters A 379 (2015) 2830–2838 Fig. 8. (Color online.) Boundary crises and the embedded behavior of the stable and unstable manifolds for γ = 0.8. For the V ∗ = π saddle point, unstable manifold is drawn in black, and stable in red (dark gray); and for the V ∗ = 2π saddle, unstable manifold is drawn in green (light gray), and stable in blue (gray). In (a) ε = 1.0 and in (b) ε = 1.07. One can see that in (a) there is no crossing between the manifolds of either saddle fixed points, and in (b) the boundary crisis occurred for both saddle points. 3.3. Manifolds Let us address now to the crises related with the manifolds. From the literature, we expect that a saddle fixed point, in the plane V vs. t has stable and unstable manifolds [29,35–37]. The unstable manifolds are formed by a family of trajectories that turn away from the saddle fixed point. One of them evolves to the chaotic bottom attractor, or visit the region of the chaotic bot- tom attractor after the event of crisis; while the other one evolves towards an attracting fixed point. These unstable manifolds are obtained from the iteration of the map T with appropriate ini- tial conditions. Similarly, the construction of stable manifold re- quires that the inverse of the mapping, say T −1, must be obtained. Here the operator follows T −1(Vn+1, φn+1) = (Vn, φn). Basically, one must replace every pair (Vn, φn) to (Vn+1, φn+1) and vice- versa in Eq. (1), and rearrange the terms as function of Vn and φn . Also, the transcendental equations should be solved as second degree ordinary functions. So, after some algebra we have T −1 : ⎧⎨ ⎩ φn = φn+1 − Vn − ν, h(Vn) = V 2 n + 2ε[cos(φn+1 − Vn − ν) − cos(φn+1)] − ν2, (8) where ν = [Vn+1 + (1 +γ )ε sin(φn+1)]/γ . Here, the function h(Vn) must be solved numerically, once it depends on both Vn and Vn+1. So, we made use of the Newton’s method to find the root, where Vn = Vn+1 −h(Vn)/h′(Vn), and h′(Vn) = 2Vn +2ε sin(φn+1 − Vn − ν). The procedure for obtaining the stable manifolds is the same as that one used for the unstable manifolds, however, instead of it- erating the map T we must iterate its inverse T −1. We set a tiny circle of radius δ = 10−4, and split 104 initial conditions around the respective saddle point and iterate the normal and reverse dy- namics. After that, we make a zoom-in the region of the saddle point, and made a linear fit in both branches of the unstable and stable manifolds, in order to find the eigenvectors. After that, we just evolve the normal and reverse dynamics again, but consid- ering the distribution of initial conditions along the linear fit of the manifolds branches. Just for notice the saddles are localized in V ∗ = mπ and φ∗ → φ∗ − π . For closed domain dynamical systems, such as the Fermi–Ulam model and other time-dependent billiards [47–50], the unstable manifolds generate the border of the basin of attraction of the chaotic attractor and the stable manifolds draw the boundaries of the basin of the attracting sink. A boundary crisis happens when the stable manifold touches the unstable manifold of the same sad- dle fixed points due to a modification of the control parameter. In this case, the chaotic attractor is destroyed, and only the sink re- mains in the system. This collision implies in a sudden destruction of the chaotic attractor and also of its basin of attraction [35,36]. This destruction can be very useful as a mechanism of controlling chaos in this dissipative version of the model since, after the crisis event, the particle is captured by an attracting fixed point (sink). In the literature, one can still find others examples of crises [6,35, 36], as crises between attractors, where a fusion between two or more attractors into a bigger one is observed, or even inner crises [6,35,36], where a chaotic attractor increases its size by colliding with a stable periodic orbit inside its own basin. However, the kind of crisis that happens in the unbounded bouncing ball model is a bit more complicated, once we have plenty of attractors, and their manifolds found themselves embedded. Stable and unstable manifolds for the first two saddle points for a dissipation value of γ = 0.8 are displayed in Figs. 8(a, b). For V ∗ = π saddle point, unstable manifold is drawn in black, and stable in red (dark gray); and for V ∗ = 2π saddle, unstable mani- fold is drawn in green (light gray), and stable in blue (gray). Both branches of the stable manifold behave as follows: the upward branch evolves to the attracting fixed point while the downward branch evolves to the chaotic attractor. In both Figs. 8(a, b), we see that the stable manifold for the saddle in V ∗ = 2π , are embed- ded with the stable manifold of the saddle V ∗ = π . Both stable manifolds are also intertwined with the attractor in the bottom. We believe that this peculiar behavior happens for all the saddles located above in the phase space for all values of V ∗ = mπ , creat- ing a whole chain of iteration between the attracting fixed points and the attractor in the bottom. Also, in both Figs. 8(a, b) the two branches of the unstable manifold generate the basin boundaries for both attracting fixed points. Here, there is no iteration between the unstable manifolds of different sinks. The unstable manifolds go up and up in the velocity axis, in a stretching and mixing way, drawing the limits of the attracting boundaries of each fixed point. One can imagine how they would behave just looking at Fig. 3, where the basin of attraction until V ∗ = 6π is drawn. It would be interesting to compare the relations between the stretching and mixing of the manifolds with the Smale horseshoe mapping. Also, in Fig. 8(a) shows the embedded manifolds ε = 1.0. One can notice here, that there is no crossings between the bound- aries of the stable and unstable manifold of both saddle points yet, but we can visualize the tangent bifurcation occurring in the sink located in V ∗ = 2π , where three branches rise and start to grow from the place where the sink should be located. One can also compare with Fig. 5, and confirm that this is the exactly control parameter of the tangent bifurcation. Once we raise the value of the control parameter ε , the respective unstable and sta- ble manifold from both saddle points cross each other, as show Fig. 8(b), for ε = 1.07. One could imagine that the crossing be- havior of the manifolds would destroy the attractor in the bottom, and let only the sink in the phase space. This indeed happens, but once the manifolds for all saddles found themselves embedded, this destruction seems not to cause greater effects in the dynamics. Because even that for one saddle the attractor is destroyed, there will be the manifold from upper saddle points, embedding them- selves with the region where the attractor was, giving an “extra life time” for the attractor, until the next boundary crisis from the upper saddles, where the same process will occur in a successive way. We think, that the boundary crisis between the manifolds, are serving as a trigger for the unbounded growth of the chaotic attractor in the bottom. One could imagine for higher values of ε , A.L.P. Livorati et al. / Physics Letters A 379 (2015) 2830–2838 2837 Fig. 9. (Color online.) Basin of attraction overlapped with the stable and unstable manifolds from the first sink, and the vibrating wall itself. One can notice that for ε ≈ εt and γ = 0.8, the basin of attraction of the first sink is tangent to the growing branch of the stable manifold embedded with the chaotic bottom attractor, gener- ating a tangent bifurcation inside the attracting basin of the first sink. Yet, one can look at and see that from the place where the sink should be located, are rising three new branches, as a result from the tangent behavior between the manifolds. The comparison with the basin of attraction overlapped, shows better this behavior. that the bottom chaotic can be very influential in the dynamics. If, we take ε = 10, the bottom chaotic attractor would occupy the region where the first three sinks of period-1 should be located, and the crossing between the manifolds that we are seeing here in Fig. 8 for low values of ε , would be happening for very high sad- dle points in the phase space. We dare to call this mechanism, as a regenerating process, or “robust attractor”. In time, this ‘robust’ term is not yet fully understood in the literature and still needs more deeper investigations. Still, we display the tangent bifurcation for the first sink in Fig. 9, where for ε ≈ εt , we overlap the basin of attraction with the stable and unstable manifold from the first saddle and the vi- brating wall itself. We can see that the basin of the first sink is tangent the growing branch of the bottom attractor. At the same time, we can see three branches rising from the place where the sink should be located inside its basin of attraction. We also stress that if noise were added to the dynamics [29,51–53], the scenario of the crises might would change, once the perturbation would be different. Of course, that noise should be consider in real experi- mental devices, so as a possible future work, we would be consider the introduction of noise in the system. 4. Final remarks and conclusions The dynamics of a bouncing ball model is investigated for a high dissipation regime. Some chaotic properties were set up, and average properties of the velocity were obtained. Depend- ing of the combination of control parameters many attractors can coexist, leading to a very complex dynamics in the phase space. Basins of attraction for some attractors and their evolution with the increase of the perturbation parameter were character- ized. Increasing the perturbation parameter leads us to characterize an unusual class of boundary crises, where the attractor collides physically with the vibrating wall. In these crises we observed a reduction in the number of attractors in the phase space. These phenomena could be extended to other vibrating and impact sys- tems with high dissipation regimes, as in human stability perfor- mance [13] and microscopic vibration systems [14], in order to reduce the number of attractors (stable vibration modes) in such systems. Also, we found the attractors are in an intertwined form. This was confirmed by the drawn of stable and unstable manifolds. Here, when a crisis between manifolds occurs, it creates a suc- cessive destruction mechanism for all the attractors originated by the sinks, giving an “extra life time” to the bottom chaotic attrac- tor, turning it into a robust attractor. As a future work, it would be interesting to investigate shrimp-like structures in the param- eter space [54–56]. Also, we could consider the introduction of noise in the dynamics, to make a link with real experimental de- vices. Acknowledgements A.L.P.L. acknowledges FAPESP (2014/25316-3), CNPq and CAPES – Programa Ciências sem Fronteiras – CsF (0287-13-0) for fi- nancial support. I.L.C. thanks FAPESP (2011/19296-1) and E.D.L. thanks FAPESP (2012/23688-5), CNPq and CAPES, Brazilian agen- cies. A.L.P.L. also thanks the University of Bristol for the kindly hospitality during his stay in UK. This research was supported by resources supplied by the Center for Scientific Computing (NCC/GridUNESP) of the São Paulo State University (UNESP). References [1] R.C. Hilborn, Chaos and Nonlinear Dynamics: An Introduction for Scientists and Engineers, Oxford University Press, New York, 1994. [2] A.J. Lichtenberg, M.A. Lieberman, Regular and Chaotic Dynamics, Appl. Math. Sci., vol. 38, Springer Verlag, New York, 1992. [3] R.M. May, Nature 261 (1976) 459. [4] G.M. Zaslasvsky, Physics of Chaos in Hamiltonian Systems, Imperial College Press, New York, 2007. [5] G.M. Zaslasvsky, Hamiltonian Chaos and Fractional Dynamics, Oxford University Press, New York, 2008. [6] K.T. Alligood, T.D. Sauer, J.A. Yorke, Chaos: An Introduction to Dynamical Sys- tems, Springer Verlag, New York, 1996. [7] R.I. Leine, H. Nijmeijer, Dynamics and Bifurcations of Non-Smooth Mechanical Systems, Lect. Notes Appl. Comput. Mech., vol. 18, Springer Science & Business, 2013. [8] P.J. Holmes, J. Sound Vib. 84 (1982) 173. [9] J. Guckenheimer, P.J. Holmes, Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields, Appl. Math. Sci., vol. 42, Springer Verlag, New York, 1983. [10] L.D. Pustilnikov, Theor. Math. Phys. 57 (1983) 1035. [11] L.D. Pustilnikov, Sov. Math. Dokl. 35 (1987) 88. [12] E. Fermi, Phys. Rev. 75 (1949) 1169. [13] D. Sternad, M. Duarte, H. Katsumata, S. Schaal, Phys. Rev. E 63 (2000) 011902. [14] N.A. Burnham, A.J. Kulik, G. Gremaud, G.A.D. Briggs, Phys. Rev. Lett. 74 (1995) 5092. [15] P. Dainese, P.St.J. Russell, N. Joly, J.C. Knight, G.S. Wiederhecker, H.L. Fragnito, V. Laude, A. Khelif, Nat. Phys. 2 (2006) 388. [16] M. Scheel, R. Seemann, M. Brinkmann, M. Di Michiel, A. Sheppard, B. Breiden- bach, S. Herminghaus, Nat. Mater. 7 (2008) 189. [17] M.K. Müller, S. Ludinga, T. Pöschel, Chem. Phys. 375 (2010) 600. [18] F. Spahn, U. Schwarz, J. Kurths, Phys. Rev. Lett. 78 (1997) 1596. [19] P. Müller, M. Heckel, A. Sack, T. Pöschel, Phys. Rev. Lett. 110 (2013) 254301. [20] F. Pacheco-Vazquez, F. Ludewig, S. Dorbolo, Phys. Rev. Lett. 113 (2014) 118001. [21] A.C.J. Luo, R.P.S. Han, Nonlinear Dyn. 10 (1996) 1. [22] J.J. Barroso, M.V. Carneiro, E.E.N. Macau, Phys. Rev. E 79 (2009) 026206. [23] A. Okniński, B. Radziszewski, Int. J. Non-Linear Mech. 65 (2014) 226. [24] L. Mátyás, R. Klanges, Physica D 187 (2004) 165. [25] A.L.P. Livorati, T. Kroetz, C.P. Dettmann, I.L. Caldas, E.D. Leonel, Phys. Rev. E 86 (2012) 036203. [26] E.D. Leonel, A.L.P. Livorati, Commun. Nonlinear Sci. Numer. Simul. 20 (2015) 159. [27] T.L. Vincent, A.I. Mees, Int. J. Bifurc. Chaos Appl. Sci. Eng. 10 (2000) 579. [28] T.L. Vincent, Nonlinear Dyn. Syst. Theory 1 (2001) 205. [29] S.K. Joseph, I.P. Mariño, M.A.F. Sanjuán, Commun. Nonlinear Sci. Numer. Simul. 17 (2012) 3279. http://refhub.elsevier.com/S0375-9601(15)00803-8/bib72656631s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib72656631s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib72656632s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib72656632s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib72656633s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib72656634s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib72656634s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib72656635s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib72656635s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib72656636s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib72656636s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib61646431s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib61646431s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib61646431s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib72656637s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib72656638s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib72656638s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib72656638s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib72656639s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib7265663130s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib7265663131s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib7265663132s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib7265663133s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib7265663133s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib7265663134s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib7265663134s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib7265663135s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib7265663135s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib7265663136s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib7265663137s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib7265663138s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib7265663139s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib7265663230s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib7265663231s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib7265663232s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib7265663233s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib7265663234s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib7265663234s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib7265663235s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib7265663235s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib7265663236s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib7265663237s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib7265663238s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib7265663238s1 2838 A.L.P. Livorati et al. / Physics Letters A 379 (2015) 2830–2838 [30] R.M. Everson, Physica D 19 (1986) 355. [31] J.M. Luck, A. Mehta, Phys. Rev. E 48 (1993) 3988. [32] U. Feudel, C. Grebogi, B.R. Hunt, J.A. Yorke, Phys. Rev. E 54 (1996) 71. [33] U. Feudel, C. Grebogi, L. Poon, J.A. Yorke, Chaos Solitons Fractals 9 (1998) 171. [34] L.C. Martins, J.A.C. Gallas, Int. J. Bifurc. Chaos Appl. Sci. Eng. 18 (2008) 1705. [35] C. Grebogi, E. Ott, J.A. Yorke, Phys. Rev. Lett. 48 (1982) 1507. [36] C. Grebogi, E. Ott, J.A. Yorke, Physica D 7 (1983) 181. [37] E.J. Kostelich, J.A. Yorke, Z. You, Physica D 93 (1996) 210. [38] T. Kroetz, A.L.P. Livorati, E.D. Leonel, I.L. Caldas, Phys. Rev. E 91 (2015) 012905. [39] E.D. Leonel, A.L.P. Livorati, Physica A 387 (2008) 1155. [40] A.L.P. Livorati, D.G. Ladeira, E.D. Leonel, Phys. Rev. E 78 (2008) 056205. [41] E.D. Leonel, A.L.P. Livorati, A.M. Cespedes, Physica A 404 (2014) 279. [42] D.F.M. Oliveira, M. Robnik, Phys. Rev. E 83 (2011) 026202. [43] D.F.M. Oliveira, M. Robnik, E.D. Leonel, Phys. Lett. A 376 (2012) 723. [44] D.F.M. Oliveira, E.D. Leonel, Physica A 413 (2014) 498. [45] M. Hénon, Commun. Math. Phys. 50 (1976) 69. [46] V. Franceschini, L. Russo, J. Stat. Phys. 25 (1981) 757. [47] E.D. Leonel, P.V.E. McClintock, J. Phys. A 38 (2005) L425. [48] E.D. Leonel, R.E. de Carvalho, Phys. Lett. A 364 (2007) 475. [49] D.F.M. Oliveira, E.D. Leonel, Physica A 389 (2010) 1009. [50] D.F.M. Oliveira, E.D. Leonel, Phys. Lett. A 374 (2010) 3016. [51] N.B. Tufillaro, T.M. Mello, Y.M. Choi, A.M. Albano, J. Phys. (Paris) 47 (1986) 1477. [52] K. Wiesenfeld, N.B. Tufillaro, Physica D 26 (1987) 321. [53] E.G. Altmann, A. Endler, Phys. Rev. Lett. 105 (2010) 244102. [54] J.A.C. Gallas, Phys. Rev. Lett. 70 (1993) 2714. [55] C. Bonatto, J.C. Garreau, J.A.C. Gallas, Phys. Rev. Lett. 95 (2005) 143905. [56] E.S. Medeiros, S.L.T. de Souza, R.O. Medrano-T, I.L. Caldas, Chaos Solitons Frac- tals 44 (2011) 982. http://refhub.elsevier.com/S0375-9601(15)00803-8/bib7265663239s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib7265663330s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib7265663331s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib7265663332s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib7265663333s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib7265663334s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib7265663335s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib7265663336s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib6D61s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib7265663337s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib7265663338s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib7265663339s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib666F726D616C31s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib666F726D616C32s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib666F726D616C33s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib7265663430s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib7265663431s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib7265663432s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib7265663433s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib7265663434s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib7265663435s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib7265663436s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib7265663437s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib7265663438s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib7265663439s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib7265663530s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib7265663531s1 http://refhub.elsevier.com/S0375-9601(15)00803-8/bib7265663531s1 Crises in a dissipative bouncing ball model 1 Introduction 2 The model, the mapping and chaotic properties 3 Results and discussions 3.1 Average velocities 3.2 Basins of attraction and bifurcations 3.3 Manifolds 4 Final remarks and conclusions Acknowledgements References