2.1. Methods and Experimental Setup

The experimental setup is illustrated in Fig. 1. It consists of a transparent acrylic container of 9.5 cm in diameter and 20 cm in height; it is permanently fed with water at constant rate of 5.27 mL/s. Inside the recipient there is a plastic tubing with inverted U shape working as siphon. It has 6 mm as inner diameter. One end of the tubing goes through the bottom of the container, and it is through it that the water is discharged. The other end of the tubing is inside the container and when the water level reaches the top of the inverted U, the container empties with an average rate of 20.4 mL/s. The inner end of the siphon is enclosed in a small glass cup ensuring the air income when the tubing empties. Along the container there is a millimeter graduated ruler. At the lowest level the water reaches 33.5 mm; at the highest level 101.5 mm. The discharged volume is 469 mL; we will call it “work volume” (WV). Note that a direct geometric calculation gives us a bigger volume; we need to consider the volume inside the tubing.

Figure 1. Experimental Device. The main element of the experimental device is a glass provided with a siphon inside it. The glass is fed with water at a constant rate of 5.27 mL/seconds, and as the water level increases the siphon is simultaneously filled. When the water reaches the maximum level of the siphon, the glass is rapidly emptied through it. This system is perturbed with water pulses two seconds long injected by means of a pump controlled by an Arduino board. The pulses can be either isolated or periodic. With label P1 the perturbation pump is indicated and with P2 the feeding one. 

There is another tubing inside the container that corresponds to the perturbation hose; it is 9 mm in diameter and it is connected to a Comet Elegant submergible pump controlled by an Arduino board which is controlled by a PC. With this system, we apply water pulses whose volume can be adjusted by the voltage applied to the pump. Pulses can be isolated or periodically applied. Volume perturbations were given as fractions of the work volume, for example a perturbation of 0.1 corresponds to a volume of 46.9 mL.

Water level evolution was recorded by a Sony Digital video camera. To improve the resolution of the images, the water was stained with orange or blue dye. The video acquisition was done at 31 frames per second. Water levels were measured directly on the video recordings and also were analyzed with the free software ImageJ 1.48v by using the luminosity of the water mass, this measurement is proportional to the water height within the container.

2.2. Data Fitting and Simulations

The Tantalus device just described before was used to obtain the natural oscillation of the system. Ten oscillations were recorded to measure the water level evolution directly frame by frame. Data corresponding to three oscillations are shown in panel A of Fig. 2. We define the highest level of the water oscillations as the starting moment of the cycle, and name it marker event ^6,7. By using the ImageJ program, we calculated the mean value for the period resulting 126.7 seconds with standard deviation of 1.24 seconds.

Figure 2. Natural and perturbed oscillations. Panel A. Display of an experimental oscillation. A measure was taken for every 30th of a second. It is considered a complete oscillation going through the points a, b and c. The natural period of the oscillation, T0, shown in this recording has a value of 126.3 seconds. Panel B. Graph of the adjustments to the normalized oscillation. The downward section was adjusted with a parabola, while the upward section with a straight line that reflects the constant rate at which the glass is filled. Panel C. Model of instantaneous perturbation when the water level is rising. The cyan trace shows what would have been the water level evolution if no perturbation was applied. The blue trace is the evolution of the water level resulting from the perturbation. The square pulse in red shows the moment in which the water was injected. It can be noticed that the effect of the perturbation is to reduce the time of the oscillation. According to the moment when the pulse is applied the oscillation can be lengthened, shortened or not present any effect. Panel D. Experimental recording of the effect of an isolated perturbation. T0 corresponds to the natural period of oscillation, T1 to the perturbed period of oscillation which in this case is smaller. Tc is the time elapsed between the “marker event” (corresponding to the beginning of the oscillation) and the water injection. 

In order to fit a function to the water level evolution the experimental data were normalized dividing them by the natural oscillation period. Also the water level variations were normalized with respect to the maximal amplitude. The water discharge, indicated as “a-b” segment in panel A of Fig. 2 has a parabolic shape ²⁰ whose mathematical relationship is given by:

Where N is water level, t is time; (t _v , N _v ) are the coordinates of the parabola vertex; 4p = 5.11, t _v = 0.521 and N _v = -0.39. The fitting function had a correlation R ₂ = 0.996.

The fitting for the filling course of the Tantalus, “b-c” segment in panel A of Fig. 2 was done by means of the straight line:

Where m = 1.325 and b = -0.325. Panel B in Fig. 2 displays the normalized and fitted oscillations. The theoretical analysis reported in the following pages is based on these two data fittings. In panel C of Fig. 2 the fitting is used to represent the effect of a volume perturbation of 0.5 times the work volume (WV), and it is applied when the container is filling. The cyan trace corresponds to the unperturbed oscillation and, the blue trace that has been superposed corresponds to the perturbed oscillation. The thin red trace at the lower part indicates, with a square pulse, the moment of the perturbation. We are considering its effect as instantaneously moving upward the water level 0.5 WV. After this change the oscillator continues its evolution starting from the newly reached level. In this case the perturbation effect was to reduce the oscillation period. We will say that in this case the oscillator moved forward. Depending on the precise moment of the perturbation application, the oscillator advances or delays. Moreover there is even a moment with no effect.

To quantitatively measure the effect of isolated perturbations and to predict the periodical perturbations effect we calculated the PTC ^6,7. The central idea for this calculation is that the perturbation instantaneously moves the oscillator from one moment in the natural oscillation to another one. The moment just before the perturbation occurs is called old phase φ, and the moment just after the perturbation which the system moves to is called new phase φ'. In the case of the Tantalus oscillator the new phases can be analytically calculated from the old phases. Let’s denote by g the function or map that relates both quantities:

To predict the effect of periodical perturbations with a fixed period τ, we build a new function which allows to know the phase after perturbation i + 1, knowing the phase after the i.th perturbation with the following relation

Once we know the initial phase, iteration of this expression lets us know the system evolution after n perturbations. Also it lets us infer if there exist coupling rhythms between perturbations and the modified behavior of the oscillator.

2.3. Analytical Method to Find New Phases

To analytically find the dependence of the new phases in terms of old phases we have to consider that subject to perturbations the oscillator may exhibit three operation regions 1) from the beginning of the emptying to the beginning of the filling up. This corresponds to the “a-b” interval in panel A of Fig. 2, this set of phases will be called R_I; 2) from the beginning of filling up to the moment the perturbation initiates the emptying phase. This region will be called R_II; this region depends on the perturbation intensity; 3) from the end of the last described region up to the beginning of the spontaneously emptying, point “c” in panel A of Fig. 2, this set of phases will be called R_III. Given that to do the analytical calculations we only consider the times corresponding to one oscillation. We will take those times described by fittings (1) and (2) as phases.

As it was mentioned in previous section, we consider that both perturbations and its effect on the water level are instantaneous; in the same way the change from old phase to new ones occurs.

The water level N _o , corresponding to the old phase φ is related as follows:

Where the specific form of N _o function depends on the region that we are considering. The water level N _n , corresponding to new phase φ' is given by the relationship:

Where the specific form of N _n function depends on the region that we are considering. If we add δN volume, the water level change becomes:

We will use these expressions to get the relationship between the new phases and old phases. We have to remark that the functional relationship for N _o and N _n depends on the operation region we are considering. For R_I region the water level is decreasing before the perturbation, and after the perturbation it will continue decreasing. In these cases we will use (1). For R_II the water level will be increasing before the perturbation and will continue increasing after the perturbation. Here we will use (2). However, for R_III the water level will be going up before the perturbation, but after the perturbation it will start to going down. In this region we are going to use (1) and (2).

Rearranging the terms we have for each operation region:

2.4. Experimental Method to Find New Phases

There is a huge amount of publications indicating how to calculate new phases from experimental data ²¹. In panel D of Fig. 2 there is a scheme which has been obtained from an experimental recording to illustrate how new phases are calculated. The unperturbed period of oscillation is denoted byT0, the perturbed oscillation period is T1. The time interval between the marker event previous to the perturbation and the perturbation time itself is T _c . The time interval between the perturbation time and the next marker event is called cophase. The old phase is T _c /T0:

While new phase is cophase subtracted from T0 and normalized, resulting in

3.1 Theoretical Phase Transition Curves (PTC)

PTC is the relation between old phases and new phases ^6,7. Figure 3 displays a PTC obtained using the analytical expressions described above with varied intensity: 0.2 WV in panel A; 0.4 WV in panel C; 0.6 WV in panel E. Intensities are indicated in the right lower corner in each panel. In the first panel, R1, R2 and R3 indicate the three different branches in the map depending on the operation region. We can observe that branch R1 has a discontinuity which is induced by the module 1 operator; its position depends on the perturbation intensity. In fact, the discontinuity phase can be calculated, since corresponds to the moment when the descending water level returns to the maximal level. There is another discontinuity in this map, and it is present for all the Tantalus PTC.

Figure 3. Experimental and Theoretical Phase Transition Curves. For all panels the horizontal axis indicates the old phase and the vertical axis the new phase. The panels on the left side correspond to the theoretical results and those on the right side to the experimental results. Each row corresponds to the same perturbation intensity, indicated on the lower right corner of each panel. The diagonal line is the graph of the identity function. In the theoretical panels the dotted line indicates the position of the map discontinuity which is insensitive to the perturbation intensity. Nevertheless, note that the magnitude of the discontinuity grows with the perturbation intensity. The maps are considered to be formed by three branches: R1, R2 and R3. In the experimental panels the different colors indicate experiments done on different days. 

In the three left panels we have marked with a dotted thin line the moment when this discontinuity always occurs. In panel A the magnitude of this discontinuity is 0.22, up to two decimal ciphers; in panel C the magnitude is 0.42; and in panel E is 0.62. This means that the magnitude grows with perturbation intensity. This result can be analytically verified, obtaining the difference of new phase values, between the last point in branch R1 and the first of branch R2. This discontinuity, which we call “essential discontinuity”, occurs because before it perturbations delay the oscillator and, after it, perturbations advance the oscillator. This discontinuity will be the one inducing the BCB.

3.2. Experimental Phase Transition Curves

Following the procedure explained in Sec. 2.4 we did a set of experiments to measure the PTC for the same perturbation intensities as those theoretical results illustrated in Fig. 3A, C, E. Three runs were made on different days for each intensity. The perturbations were applied in five seconds intervals from the moment the container starts to fill and up to the moment it is empty, that is, during an entire cycle. Between consecutive perturbations two unperturbed oscillations were allowed. All this procedure was recorded by a digital camera, and to measure phases we took as a marker event the moment when the container starts emptying. Panels B, D and F in Fig. 3 show the corresponding results for 0.2, 0.4 and 0.6 WV intensities. Different colors represent different experiments done in different days. To facilitate comparison theoretical PTC are superposed. The coincidence between them is remarkable.

3.3. Periodical Rhythms: Theoretical prediction

To explore the dynamical evolution of the system we can choose any intensity, initial condition and a perturbation period using the function (8). In most of the cases we found that orbits converge to in fixed points or periodic orbits, that is, after certain number of perturbations the phases where the system moves around are the same. This repeated phases can be reached in one or several oscillations. The combination of the number of perturbations with the number of oscillations will be called coupling rhythm, and will be identified with N : M notation; where N = number of perturbations and M = number of oscillations. When two rhythms have the same N, we would say that they have the same periodicity.

Figure 4 shows some orbits obtained with 0.5 WV perturbation intensity. Panel A displays the orbit the system follows when 0.7T0 perturbation period and an initial condition of 0.32T0 were used. It can be seen that after a transient course of five iterations, the orbit converges to a point where the identity line intersects the map. Staying in this point means that before each perturbation the system is always at the same phase. This is a rhythm with 1-periodicity. This type of evolution occurs because the intersection of identity line with the map has an absolute value slope lower than one, in consequence, all orbits passing around that vicinity are attracted. This fact allows us to analytically calculate those 1-periodicity cases.

Figure 4. Resulting orbits from iterating the Phase Transition Curve for perturbation intensity of 0.5 WV. Panel A. Resulting orbit from iterating with a perturbation period of 0.7T0 and initial condition 0.32T0. After some iterations the orbit move to the crossing of the map and the identity function. A rhythm with 1-periodicity results given that each iteration falls in the same phase. Panel B. Orbit with 2-periodicity, resulting from iterating with perturbation period of 0.35T0, the orbits visit two phases of the map. In this panel, as in the following ones, the initial condition is omitted and only the asymptotic part of the orbit is illustrated. Panels C and D. Orbits with 3-periodicity and 4-periodicity resulting from iterating with perturbation periods of 0.24T0 and 0.18T0. 

Panels B, C and D illustrate rhythms with higher periodicity, though transitory paths have not been represented, only the final orbits. Panel B presents a 2-periodicity rhythm obtained when perturbation period is 0.35TO. Panel C shows a 3-periodicity orbit occurring for perturbations applied 0.24T0 and, D a 4-periodicity orbit obtained for 0.18T0. It can be observed that between identity line and RI branch in the map there is a channel, and by adequately choosing the perturbation period it can be made as thin as desired. Reduction in the width of this channel will induce very long orbits.

3.4 Periodical Rhythms. Experimental results

To corroborate the theoretical predictions established in the last section, we perturbed the Tantalus oscillator with an intensity equal to 0.5 WV, which matches the case shown in Fig. 4. This means that volume pulses of 234.5 mL are injected in 2 seconds. In all the experiments the initial condition or first perturbation pulse phase was selected among the several points in the studied orbits. Then for perturbation period 0.7TO the first pulse was applied in phase 0.68T0. Panel A of Fig. 5 shows that the result in this case is a coupling with 1:1 rhythm. Panel B in the same figure shows the result for 0.33T0 perturbation period and an initial condition of 0.27T0, the obtained rhythm is 2:1. In panel C the coupling rhythm is 3:1, a result from perturbing with a 0.24T0 period and an initial condition of 0.17T0. Finally, in panel D a coupling rhythm 4:1 produced by perturbation pulses with 0.17T0 period and initial condition 0.81T0 is shown. Once again the coincidence between theory and experiments is remarkable. Not only are the predicted rhythms obtained; but they can also be reached by starting the perturbation protocol in the predicted phases. It has to be noted that each experiment was repeated at least three times.

Figure 5. Resulting couplings from applying periodic perturbations with intensity 0.5 Experimental recordings resulting from applying periodic perturbations with periods and intensities as in Fig. 4. Points of the theoretically calculated orbits were taken as initial conditions for each experiment. In each panel the horizontal axis indicates time in seconds. In the vertical axis indicates the water height in millimeters. In the upper part of each panel the perturbation period and resulting rhythm is indicated, where the first digit indicates the number of perturbations and the second one the number of oscillations in which the patterns is repeated. 

3.5. Bifurcation Diagrams for the Perturbation Period

By repeating the exercise that we have described in Sec. 3.3, but now for intensities going from 0.1 to 1.05 in 0.01 step size; for 1000 periods between 0 and 1 with 0.001 steps size and for 10 initial conditions between 0 and 1 with 0.1 steps size, we completed 950 000 studied cases. We looked for periodicities between 1 and 16 making 2560 iterations in each case. Most of the cases: 871 211 had periodicity between 1 and 16. The remaining cases could have longer periodicities. In Fig. 6, panels A and B, we show the periodicities pattern obtained for perturbation intensities 0.2 WV and 0.4 WV. It can be seen that the periodicity patterns have a general horseshoe shape, built by rungs that constitutes a staircase. At the right and left ends, outside the horseshoe, there are 1-rhythm regions; these rhythms are consequence of the intersection between the identity line and the mapping, as has been described in Sec. 3.3. As perturbation intensity grows these regions grow, which reflects that R3 grows with the perturbation intensity, as can be seen in the left panels of Fig. 3.

Figure 6. One-parametric Bifurcation Diagrams: increasing, additive and with bistabilities zones. In panels A and B are shown the bifurcation diagrams obtained when perturbing with periods between 0 and 1, in steps of 0.0001. In the first panel the perturbations are done with an intensity of 0.2 WV and in the second one of 0.4 WV. Panel C. Detail of panel A, for periods between 0.25T0 and 0.45T0. It can be noted that between the regions with 2-periodicity and 3-periodicity there exists a gap filled with periodicities that sum the periodicities of the edges, this way we have a region with periodicity 3+2=5, another one with 3+5=8, etc. Panel D. Detail of panel B, for periods between 0.17T0 and 0.42T0. It can be noticed that region of 2-periodicity and the one of 3-periodicity overlap, which predicts the existence of bistabilities. 

In general, we have “additive periodicities”, that is, consecutive periodicities do not have contiguous domains of rhythms ²². This characteristic is illustrated in panel C of Fig. 5. Here we show that between 2-periodicity and 3-periodicity there is a gap. This gap is partially occupied by a 2 + 3 = 5-periodicity region. Between 2-periodicity and 5-periodicity zones there is a 2 +5 = 7-periodicity zone; and so on and so forth. Panel D in the same figure shows another type of phenomena commonly found in Periodicity Diagrams or One-parametric Bifurcation Diagrams. In the region between perturbation periods 0.2T0 and 0.4T0 and perturbations intensity 0.4 WV, there is a region in which when changing the initial condition, we get different coupling rhythms, that is, there is bistability and this property is very common in the Tantalus Oscillator diagrams.

3.6 Experimental verification of Bistability and Additive Rhythms

In panel A of Fig. 7 we show a 2-rhythm coupling obtained by using a perturbation period of 0.4T0, we have two perturbations for each Tantalus oscillation. Repeating the procedure for perturbation period 0.25T0 we get a 3-rhythm coupling, three perturbations for one Tantalus oscillation, this is shown in panel B in the same figure. Finally, when we chose an intermediate perturbation period, 0.3T0, which according to the results displayed in panel C of Fig. 6 should correspond to a 5-periodicity pattern, we see that the experimental result is five perturbations coupled to two Tantalus oscillations.

Figure 7. Experimental recordings: additive and bistable rhythms. In the first three panels the addition phenomenon is illustrated, obtained with a perturbation intensity of 0.2 WV. In the last two panels the bistability phenomenon is illustrated, obtained with a perturbation intensity of 0.4 WV. In all the panels the blue trace shows the evolution of the water level, the red trace indicates the perturbation moments. Panel A. Rhythm 2:1 obtained with a perturbation period of 0.4T0; we have two perturbations per oscillation. Panel B. Rhythm 3:1 obtained with a perturbation period of 0.25T0; we have three perturbations per oscillation. Panel C. Applying perturbations with period in between the previous ones: 0.3, a rhythm 5:2 is obtained, which is the combination of the elemental blocks for 2:1 and 3:1. We have 5 perturbations in two oscillations. In panels D and E perturbations with the same intensity and period were applied: 0.4WV and 0.28T0. In panel D the initial condition was 0.25T0 and the resulting rhythm was 2:1. In panel E the initial condition was 0.22T0 and the obtained rhythm 3:1. 

In relation to the experimental verification for the existence of predicted bistabilities (panel D of Fig. 6), we have explored for perturbation intensity 0.4 WV and different perturbation periods overlapping their respective regions. In panels D and E of Fig. 7 we show the results obtained after applying repetitive pulses with a perturbation period of 0.28T0 and initial condition of 0.25T0 (panel D) and the same parameters values but with initial condition 0.22T0 (panel E). It can be observed that we get two different couplings: 2:1 and 3:1.

3.7. Border Collision Bifurcations

Bifurcation Diagrams described in the previous paragraph in which the system “jumps” from one periodicity to another pose the question: what kind of bifurcations occur between one periodicity and the other? In most cases we are studying the bifurcations type is Border Collision and it happens because the orbits collide with the essential discontinuity described at the beginning of this results section. In Fig. 8 we illustrate how this kind of bifurcation occurs when the perturbation intensity is 0.2 WV, perturbation period goes from 0.182T0 to 0.154T0 and the initial condition is 0.1. In panel A the 4-period orbit when perturbation period is 0.182T0 is shown.

Figure 8. Border Collision Bifurcations. Transition between 4-periodicity and 5-periodicity during the period reduction. In Panel A the obtained orbit for a perturbation period of 0.182T0 is shown; it has 4-periodicity. The third element of the orbit is identified with a blue arrow and the letter “a”. Panel B. The reduction the perturbation period to 0.156T0 approaches this third element to the map discontinuity, now marked labeled with letter “b”. Panel C. A very small reduction makes the third element of the map collide with the map discontinuity and reorganizes the orbit with now five elements. Panel D. In the horizontal axis the perturbation period is indicated, in the vertical axis the phases the iteration visits are indicated. Phases with 4-periodicity are in black and those of 5-periodicity in red. The horizontal blue line indicates the discontinuity position. Observing the graph from right to left it can be seen that the 4 periodicity starts in the fourth phase and ends in the third one. At that moment the 5 periodicity starts in its fourth phase to end in its third one. We are counting the phases from top to bottom. 

The arrow and “a” letter are indicating one of the four elements in the orbit. As the perturbation periodicity is decreasing, this element is approaching the discontinuity. In panel B we see its new position indicated with an arrow and “b” letter, when the perturbation period is 0.156T0. At that moment it is almost touching the discontinuity. Another small reduction in the perturbation period causes the orbit to touch the discontinuity and moves from branch R2 to branch R1, reorganizing the orbit path that now has 5-periodicity. In panel D of the same figure, we show in a general fashion these bifurcations. We plot the phase values for each orbit element for perturbation periods from 0.12 to 0.2. In red we show phase values corresponding to 5-rhythm, in black the phase values corresponding to 4-rhythm. A thin horizontal blue line marks the essential discontinuity position; arrows indicate the position of a, b and c points marked in panels A, B and C. If we analyze the pattern of moving phases going from right to left along the blue line, we can see that 4-rhythm begins at the discontinuity point for the fourth phase (if considered from top to bottom), and it finishes to the end of third phase in the same discontinuity line. For the fourth phase, the 5-periodicity rhythm emerges which ends in the same discontinuity for the third phase.

A detailed inspection of the phase diagrams for each studied perturbation intensity and for each initial condition shows all bifurcations occur in this way when we plot the phases between 0 and 1. Then we can say that most of the bifurcations occurring in this system are BCB, without discarding that there may be some other bifurcation types.

3.8 Two-parametric Bifurcation Diagram

We can group the One-parametric Bifurcations Diagrams as those shown in panels A and B of Fig. 6, to build a TBD. In panel A of Fig. 9 this construction is shown. The horizontal axis indicates perturbation periods from 0 to 1 in step size of 0.001; in the vertical axis we represent the perturbation intensity between 0.1 and 1.05 WV in step size of 0.001; in total we have 950 000 points. Periodicities up to fourteen are represented giving the same color to these rhythms with the same periodicity. It can be seen that in general, the diagram has the same appearance as those obtained by Chialvo ¹⁸ and Santillán ¹⁹. There are large zones with 1:0 and 1:1 rhythm at both sides. The region with higher rhythm has a triangular shape with its vertex pointing towards higher intensities. The border zone with 1:0 rhythms is barely curved, while the border with the 1:1 region is a straight-line.

Figure 9. Two-parametric Bifurcation Diagram and Bistabilities Diagram. In both panels the horizontal axis indicates the periodicity and the vertical axis the perturbation intensity. Panel A. Regions with similar periodicity are indicated with the same color. Yellow for periodicity 1. White regions have periodicity greater than 14. Panel B. Regions with different periodicities overlap, this indicates that there are many bistability regions. These regions are indicated in black. 

The pattern inside the triangle is formed by bands, that is, sets of points with the same rhythm. As some of the rhythms share the same periodicity, we assigned all of them the same color. For example, looking at the diagram from left to right we see 5:1, 5:2, 5:3 and 5:4 rhythm bands. To all of them we assigned the red color. At the right side of the figure it appears the color code for periodicities, reported up to 14. Yellow color corresponds to 1-periodicity. On the left side of the figure is the 1:0 coupling, meaning that the system receives periodic perturbation in the same phase but never completes an oscillation. On the right side the coupling is 1:1, implying an oscillation is produced for each volume perturbation applied in the same T0 phase.

Although the diagram changes as the intensity perturbation grows, there are some general characteristics. A notable one is that periodicities accumulate to the left and right of the diagram boundaries. In both cases, as we approach to the borders, longer periodicities are observed, but each time in narrower zones. This property could be observed in the One-parametric Bifurcation Diagrams of panels A and B of Fig. 6. Another general characteristic is that a band rhythms goes to the vertex of the triangle becoming thinner but preserving their ordering. Zooming the visualization of the diagram at different heights, the same band ordering is observed. A novel result respect to Chialvo and Santillán ^18,19 papers, is the existence of bistabilities in our system. In panel B of Fig. 9 they are shown. We can see that besides a wide zone between 2:1 and 3:1 rhythm; in the middle of the diagram there are bistabilities too. They accumulate towards the boundaries. White zones in the diagram are partially ought to bistabilities of very high order not detected in our simulations.