Crises and chaotic scattering in hydrodynamic pilot-wave experiments

When a fluid bath vibrates vertically with an amplitude above a critical value, the surface undergoes an instability leading to the spontaneous formation of the so-called Faraday waves.^1,2 Nearly two decades ago, Yves Couder and co-workers^3,4 showed that when a bath of silicone oil is vertically vibrating with an amplitude slightly below this point of instability, a millimeter-sized droplet of the same fluid can bounce vertically and “walk” (move horizontally) on the surface without coalescing due to the presence of a thin air layer between the droplet and the bath surface. In this regime, the impulse that the droplet experiences at each impact is determined by the surface topography, which itself is shaped by the waves generated at previous bounces. Although these waves do not persist and decay exponentially, they do so at a rate that vanishes as the vibration amplitude approaches the Faraday threshold, thus retaining a “memory”⁵ of the droplet’s trajectory. Soon after the initial demonstrations of bouncing and walking, Couder and co-workers realized macroscopic analogs to quantum mechanical phenomena such as tunneling and orbital quantization on these setups.^6,7 Their results were extended in subsequent studies^8–12 and the setups were named “hydrodynamic pilot-wave systems”¹³ due to their reminiscence of the de Broglie–Bohm interpretation of quantum mechanics. Today, the exploration of analogies between pilot-wave hydrodynamics and quantum phenomena continues to be an active area of research (see the recent review by Bush and Oza).¹⁴ 

Over the past two decades, pilot-wave hydrodynamics have also received considerable attention for their dynamical properties.¹⁵ Hydrodynamic pilot-wave systems can be viewed from two different theoretical perspectives that are compatible with one another. If one approximates the droplet as a point particle and takes its position and momentum as state variables, then the corresponding dynamical system is non-Markovian since the time-evolution of the system not only depends on its present state but also on its history which determines the surface topography.^5,16 Alternatively, if one takes the droplet and fluid bath as a whole and includes the bath’s configuration in the description of the system’s state,^17,18 then the knowledge of the present state becomes sufficient for predicting its future, thus rendering the system Markovian. In this latter approach, the hydrodynamic pilot-wave setup constitutes an infinite-dimensional dynamical system since the surface field approximated as a continuum introduces infinitely many degrees of freedom. In the following, we adopt this infinite-dimensional Markovian view in order to explain how complex dynamics can arise in a hydrodynamic pilot-wave system.

The view of complex hydrodynamic phenomena as those arising in infinite-dimensional dynamical systems can be traced to Hopf’s early mathematical work in turbulence.¹⁹ Over the past 30 years, this approach enjoyed a resurgent interest mainly due to the advancements in computing hardware and numerical algorithms that made searching for simple time-invariant solutions such as equilibria and periodic orbits in fully resolved Navier–Stokes simulations feasible.^20,21 In addition to numerical simulations, the influence of invariant solutions on weakly turbulent flows was also demonstrated in experiments.^22,23 One of the frequently stated goals of this research field is building low-dimensional models of turbulence using (unstable) periodic orbits as building blocks.²⁴ Although several papers^23,25–28 illustrated the resemblance of turbulent flows and periodic orbits in different fluid flows, building quantitatively accurate models based on periodic orbits has only been possible in highly restricted configurations.²⁹ In this paper, we explore an alternative approach to this problem, namely, one that takes the chaotic repellers as building blocks as opposed to periodic orbits, and demonstrate its success through our analysis of hydrodynamic pilot-wave experiments.

Finite-lifetime chaotic motion can be observed in a variety of settings, such as chemical reactions,³⁰ advection of suspended particles,³¹ and dynamics of transitionally turbulent flows.³² Among these, a well-studied phenomenon is chaotic scattering, that is, the scattering of a particle from the neighborhood of a nonattracting chaotic set in the phase space of a system.^33–36 Although initially studied in open billiard systems,^33–35 chaotic scattering found applications in various other systems with chaotic repellers (for recent examples, see Refs. 37–39). In the present paper, we demonstrate that the chaotic dynamics of a hydrodynamic pilot-wave system can be decomposed into chaotic repellers formed around distinct periodic solutions of the system. Consequently, we show that the observed dynamics can be understood as consecutive scatterings from these repellers.

Our experiments consist of a single droplet bouncing on a bath with variable topography that introduces a radially confining force, which is known to enable chaotic dynamics in hydrodynamic pilot-wave systems.^17,40,41 In this setup, we slowly change the control parameter to uncover the series of bifurcations that lead to the formation of the system’s chaotic attractor. In particular, we observe crisis bifurcations,⁴² i.e., discontinuous changes of the system’s attractor upon small changes of the control parameter. In the final chaotic regime, we study lifetime statistics of dynamics in different parts of the attractor and show that they have exponential tails as expected for chaotic repellers.⁴³ Our results suggest that the final crisis bifurcation that leads to the formation of the system’s attractor is mediated by a chaotic repeller, in contrast to well-studied^36,42,44 crisis phenomena that take place when a chaotic attractor of a continuous-(discrete-)time system collides with a periodic orbit (fixed point) of saddle type. We, thus, argue that our findings open new theoretical questions for high-dimensional chaotic systems.

This article is structured as follows. In Sec. II, we discuss the experimental setup, a symmetry-reducing transformation of experimental data, and the construction of Poincaré sections. In Sec. III, we examine how droplet dynamics change as the memory $Me$ of the vibrating fluid bath is varied. We also rationalize our findings using tools from chaos theory. Finally, we discuss the significance of our findings in the broader context of hydrodynamic quantum analogs as well as fluid turbulence in Sec. IV.

Our setup consists of a computer-controlled electromagnetic shaker on which a bath containing silicone oil is mounted and a camera above records the dynamics of the bouncing droplet. As shown in Fig. 1, the bath is in the shape of a circular corral formed by concentric cylinders with a deep inner section surrounded by a shallow damping (overflow) region, a configuration that was shown to yield an effective radially confining force.^8,9,45 All experimental runs are performed by setting the bath’s vibration frequency to $f0=75Hz$ and using a single droplet with a diameter $D=0.85±0.05mm$⁠. The bath acceleration $γ$ is varied to adjust the memory,

$Me$ is a dimensionless control parameter that is proportional to the damping time of the surface waves⁵ and tends to infinity as the bath acceleration $γ$ approaches the Faraday instability threshold at $γF$⁠. In our experiments, we adjust $γ$ according to (1) such that Me is varied in approximately equal steps and study the changes in the system’s behavior. In presenting our results, we convert length and time to dimensionless quantities by measuring them in units of the Faraday wavelength $λF=5.3mm$ and time $tF=2/f0$⁠, respectively. The droplet trajectories are reconstructed via image processing for detecting the droplet center in each frame recorded by the camera [Fig. 1(b)], and the droplet’s instantaneous velocity is estimated by computing the time-derivative of the cubic splines that fit these trajectories. Further experimental details can be found in the  Appendix.

Figure 2(top) shows examples of reconstructed droplet trajectories in our experiments at different Me. While visualizations in the laboratory frame are illustrative, they contain redundant information due to symmetries of the corral, i.e., azimuthal reflection and rotations. As a result of the reflection symmetry, for each clockwise rotating orbit, the system also exhibits a dynamically equivalent counterclockwise-rotating one. In the case of the circular orbit in Fig. 2(a) for instance, the clockwise and counterclockwise orbits can be distinguished using the sign of their angular momentum $L=m(x×v)$⁠, where $m$ is the mass of the silicone oil droplet, and $x(t)=[x(t),y(t)]$ and $v(t)=[vx(t),vy(t)]$ are its instantaneous two-dimensional position and velocity, respectively. These two orbits trace approximately the same trajectory on the $xy$-plane. The continuous rotation symmetry poses a different challenge; any non-circular trajectory [e.g., lemniscate in Fig. 2(a)] can be rotated by an arbitrary angle about the origin to generate a dynamically equivalent trajectory at a different orientation. In order to eliminate this degeneracy, we apply Budanur and Fleury’s¹⁸ continuous symmetry reduction method to our experimental data. The basic idea of this method is to set the polar angle in the velocity plane of the pilot-wave system’s state space to a fixed value so that each rotation-equivalent trajectory is mapped to a single representative. Formally, this transformation is applied to all dynamical degrees of freedom, including those encoding the state of the bath. In the present experimental case, we neglect the bath’s degrees of freedom which are not essential for our discussion to follow and transform the measured coordinate $x(t)$ and velocity $v(t)$ as

where $θ~(t)=arg⁡(vx(t)+ivy(t))$ is the instantaneous polar angle in the velocity plane and $R(θ)$ is the $2×2$ matrix,

whose action rotates a two-dimensional vector in counterclockwise direction by $θ$⁠.

By fixing the phase of the measurements on the velocity plane to $0$⁠, transformation (2) maps all rotation-equivalent measurements to one with $v~x=‖v‖$ and $v~y=0$⁠, hence performing a symmetry reduction. Note that the transformation (2) is defined as long as the speed $‖v(t)‖$ of the droplet does not vanish, which is the case for the dynamical regime of interest (see Fig. 8 in the  Appendix). We note that the analogous transformation that maps coordinates $x$ to $(‖x‖,0)$ cannot be applied since this transformation is singular at the origin $x=y=0$⁠, which is approached by the lemniscate-shaped trajectories [Figs. 2(a) and 2(b)].

Figure 2 (bottom row) shows the symmetry-reduced trajectories corresponding to the panels above. As seen in Fig. 2(a) (orange), symmetry reduction maps the circular trajectory to approximately a point corresponding to the one at which the droplet velocity is in positive $x$ direction.

In the case of the lemniscate orbit, the trajectory in Fig. 2(a) and all of its rotation copies are mapped to the symmetry-reduced lemniscate shown in the bottom panel. Another feature of our symmetry reduction can be understood by noting the apparent symmetry of the symmetry-reduced lemniscate in Fig. 2(a) (bottom) with respect to the $y~=0$ line. This follows directly from the fact that the azimuthal reflection symmetry of our pilot-wave system is equivalent to the transformation $y~→−y~$ after symmetry reduction. Consequently, the lower and upper halves of the $x~y~$-plane correspond to droplet motions with positive and negative angular momenta, respectively. This fact is also observed in the circular [Fig. 2(a), bottom] and oval [Fig. 2(c), bottom] trajectories, where the counterclockwise-rotating (⁠$L>0$⁠) trajectories are confined to the lower half of the $x~y~$-plane after symmetry reduction.

In order to identify the sequence of bifurcations, we define the Poincaré section⁴⁶ as the position-velocity pairs $(x~P,v~P)$ on the symmetry-reduced trajectories that satisfy the half-plane condition,

which corresponds to taking the intersections of the trajectories $(x~(t),v~(t))$ with the $y~$ axis in the positive-$x~$ direction. We chose this section since it is intersected by all of the symmetry-reduced trajectories (Fig. 2, bottom row) that we observe. In our results to follow, we approximated these intersections as linear interpolations between the symmetry-reduced data points at times $ti$ and $ti+1$ with $x~(ti)<0$ and $x~(ti+1)>0$⁠.

Because we only measure droplet trajectories, our Poincaré section (4) in the space of symmetry-reduced trajectories $(x~,y~,v~x,0)$ corresponds to a two-dimensional plane $(0,y~,v~x,0)$⁠. Consequently, the Poincaré section (4) retains only a subset of the system’s state due to our neglect of the bath’s degrees of freedom. Formally, an equivalent Poincaré section could have been defined for the infinite-dimensional state measurements containing the bath’s degrees of freedom as done for the numerical data by Budanur and Fleury.¹⁸ Section (4) should be understood as a projection of that Poincaré section onto a two-dimensional plane, and as we shall demonstrate next, the retained information is sufficient for our analysis.

In the following, we present the results of a parametric study wherein we vary Me in small increments to follow the changes in the symmetry-reduced dynamics.

Figure 3 shows the orbit diagram which we obtained by recording the symmetry-reduced droplet trajectories’ intersections with the Poincaré section (4) starting from the lemniscate trajectory at $Me=20$ and increasing $Me$ in small steps of $ΔMe≈0.36$ up to $Me=32.8$ and reversing the direction of parameter sweep afterward. Top (bottom) panels correspond to the experimental run where Me was increased (decreased).

For the stable lemniscate solution at $Me≈20$⁠, we obtain two (localized) intersections on the Poincaré section at $y~P/λF≈±0.025$⁠. Upon increasing the Me, this solution loses stability and a chaotic lemniscate attractor [Fig. 2(b)] takes its place. This transition to chaos results in the expansion of the $y~P$ interval spanned by the intersections in Fig. 3 (top panel), starting at $Me≈22.5$⁠. Upon further increase in Me, the chaotic lemniscate attractor loses stability at $Me≈26.6$⁠, after which we observe circular trajectories similar to the orange curve in Fig. 2(a). On the orbit diagram, the corresponding intersections with the Poincaré section are localized at $y~P/λF≈−0.08$⁠. When we increase Me further, we observe modulations to the circular trajectories and the formation of a new oval chaotic attractor [Fig. 2(c)]. This is once again indicated by the gradual expansion of the markers in the orbit diagram at $Me≈29.5$⁠. Finally, we observe a sudden expansion of the attractor at $Me≈31.65$⁠, which results in droplet trajectories intermittently switching between oval motions and lemniscates as shown in Fig. 2(d). Qualitative features of this chaotic motion remain unchanged for the higher Me values shown in Fig. 3. Decreasing Me from this point (bottom panel of Fig. 3), however, yields a different scenario. Although at $Me≈31.65$⁠, the dynamics fall back onto the oval chaos [Fig. 2(c)], we do not recover the lemniscates upon further decrease of Me . Rather, we follow the branch of stable circular trajectories down to $Me=20.0$ as shown in Fig. 3. In other words, for the parameter interval $Me∈[20.0,26.6]$⁠, we have a multistable system with distinct branches of stable dynamics that can be observed depending on how the system is initiated.

Sudden changes in the orbit diagram (Fig. 3) indicated by the dashed red lines can be understood as crisis bifurcations⁴² at which the attractor of a nonlinear system undergoes a discontinuous change upon a small variation of the control parameter. The first of these bifurcations takes place when the lemniscate chaos [Fig. 2(b)] loses stability at $Me≈26.6$⁠. This so-called boundary crisis⁴² generically takes place when a chaotic set intersects with its basin boundary. Usually, this boundary is the stable manifold of another solution, such as a periodic orbit. Naturally, in an experimental study, we cannot probe unstable solutions. Nevertheless, we note that multistability of the system for $Me∈[20.0,26.6]$ is in agreement with this scenario. The second crisis bifurcation takes place at $Me=31.65$⁠, when the chaotic attractor with oval-shaped trajectories [Fig. 2(c)] suddenly expands into the upper-half of the $y~P$ axis in Fig. 3. Recalling that the azimuthal reflection symmetry is represented by the sign change of the $y~$ coordinate, we conclude that this crisis bifurcation is a symmetry-restoring⁴⁷ one at which the previously disconnected state-space regions that are reflection copies of one another merge to form the final chaotic attractor [Fig. 2(d)] of the pilot-wave system. The appearance of lemniscate-shaped trajectories [Fig. 2(d)] following the symmetry-restoring bifurcation suggests that it is the result of a merger of oval-shaped chaotic trajectories with lemniscate ones, which lost stability via boundary crisis at $Me≈26.6$⁠. This scenario, illustrated as a state-space cartoon in Fig. 4, is markedly different from previously studied^36,42,44 crisis bifurcations that follow a chaotic attractor’s collision with an unstable periodic orbit.

In our bifurcation study, we observed seemingly chaotic trajectories in the neighborhoods of both lemniscate and oval orbits. As our final result and further evidence for the bifurcation scenario depicted in Fig. 4, we demonstrate that the system’s dynamics following the symmetry-increasing bifurcation at $Me≈31.65$ can be understood as chaotic scatterings from the state-space regions corresponding to the lemniscate- and oval-shaped orbits. To begin, we track the time series of the angular momentum $L$ that distinguishes clockwise and counterclockwise-rotating trajectories. Figure 5(a) shows the angular momentum time series (blue) corresponding to a chaotic trajectory segment similar to the one shown in Fig. 2(d). Intervals, such as $t∈[500TF,1000TF]$⁠, during which the angular momentum oscillates between positive and negative values of $L$ correspond to lemniscate-shaped orbits. Conversely, the episodes during which the angular momentum remains either positive or negative correspond to the oval-shaped trajectories. This indeed suggests that the lemniscate orbits should separate the chaotic ovals with opposite senses of rotation, just prior to the crisis bifurcation, and, hence, constitute a chaotic basin boundary.

If the chaotic attractor can be decomposed into distinct chaotic repellers, then we would expect the droplet to spend exponentially distributed times within each of these sets.⁴³ In order to test this hypothesis, we first computed the moving average $⟨L⟩Tw$ of the angular momentum, shown in Fig. 5(a), with a window length $Tw∈[104,113]TF$ (see the  Appendix for details). During the lemniscate episodes, the window-averaged angular momentum remains near zero, whereas for the ovals, it takes a near constant value $⟨L⟩Tw≈±0.05mλF2/TF$⁠. Following this observation, we choose the thresholds $Ll=0.025mλF2/TF$ and $Lo=0.045mλF2/TF$ and assume that the episodes with $|⟨L⟩Tw|<Ll$ [Fig. 5(a), transparent red] and $|⟨L⟩Tw|>Lo$ [Fig. 5(a), transparent cyan] correspond lemniscate and oval trajectories, respectively. Under this assumption, we estimated the distribution of the lifetime $τ$ for the lemniscate trajectories at four Me values as shown in Fig. 5(b). The survival function $S(τ)$ is the probability of the droplet to remain in lemniscate-like motion for a time $τ$⁠, and for all Me values, it appears to have an exponential tail for $τ>100TF$⁠. Interestingly, this distribution shows very little variation upon changing Me; thus, we fit an exponential to the mean of $S(τ)$ in Fig. 5(b) for $τ≥125TF$⁠. The lifetimes of the oval trajectories shown in Fig. 5(c), on the other hand, become progressively shorter as we increase Me. Similar to the lifetime distributions of the lemniscates, the ovals also show exponential tails for high $τ$⁠.

A lifetime distribution with an exponential tail is the hallmark⁴³ of transient chaos.^48–50 Intuitively, one can understand the necessity of an exponential-tailed lifetime distribution by arguing that a chaotic system is memoryless⁵¹ for time scales much longer than the Lyapunov time $μL−1$⁠, where $μL$ is the leading Lyapunov exponent. Since any noise in the system is typically amplified as $exp⁡μLt$⁠, the system’s memory of its present state will be completely lost after a time $t≫μL−1$⁠. Therefore, the exponential tails of the lifetime distributions in Figs. 5(b) and 5(c) lead us to the conclusion that the lemniscate- and oval-shaped motions of the droplet correspond to distinct chaotic repellers, which are visited transiently by the dynamics. The overall dynamics can, thus, be viewed as consecutive scatterings between these strange repellers. In Fig. 6(a), we illustrate one such scattering event where the droplet transitions from the oval-repeller with positive angular momentum (⁠$⟨L⟩Tw>Lo$⁠) to the negative side after spending some time on the lemniscate repeller. In order to illustrate these different regions in Fig. 6, we also plotted samples from a very long trajectory with different colors corresponding to distinct state-space regions (blue: lemniscate, orange: oval with $⟨L⟩Tw>Lo$⁠, green: oval with $⟨L⟩Tw<−Lo$⁠).

Prior to the symmetry-increasing crisis bifurcation at $Me≈31.65$⁠, the oval-shaped chaotic trajectories constitute an attractor of our pilot-wave system. The symmetry-increasing bifurcation is, therefore, one at which the chaotic ovals lose stability and become transient. The post-crisis reduction of the mean lifetimes of oval-shaped orbits [Fig. 5(c)] is reminiscent of the well-understood behavior in dissipative two-dimensional maps.^42,44 In those systems, attractors generically lose stability through homoclinic or heteroclinic tangencies of periodic orbits. Subsequent mean lifetimes of the chaotic transients drop off following a power-law as a function of the distance from the critical parameter value at which the crisis bifurcation takes place. In the present case, our measurements near the critical Me are not dense enough to quantitatively test a power-law behavior. Moreover, the chaotic attractor in the present case loses stability through its merger with the chaotic lemniscates, for which, to the best of our knowledge, there is no obvious reason to expect a power-law scaling. Nevertheless, it is clear that the lifetimes of ovals become progressively shorter as we increase the Me further from its critical value. We would like to note that this behavior can be exploited to tune the probabilities of observing the droplets in oval and lemniscate states.

In this paper, we presented a parametric study of a hydrodynamic pilot-wave system, in which we identified two crisis bifurcations. The first of these occurs when chaotic lemniscate trajectories lose stability giving way to stable circular-like trajectories. The second symmetry-increasing one occurs when a chaotic attractor corresponding to oval trajectories merges with the lemniscate repeller, yielding the final attractor of the system. We have also demonstrated that the proximity of the control parameter Me to its critical bifurcation value determines the lifetime distribution of the oval-shaped motions, providing a useful tool for adjusting the global statistics of similar pilot-wave systems. We note that global bifurcations due to the merger of distinct chaotic sets were previously observed in the discrete-^52–54 and continuous-time¹⁸ models of the pilot-wave hydrodynamics.

The idea of decomposing high-dimensional chaos into subunits of chaotic repellers by identifying crisis bifurcations has recently been explored in the numerical studies of transitionally turbulent shear flows.^55–57 Similarly, chaotic basin boundaries were also investigated in the numerical studies of laminar-turbulence transition in shear flows.^58–61 The instability of these basin boundaries, however, renders them inaccessible to direct observations in the laboratory. To the best of our knowledge, our study is the first experimental demonstration of transitions mediated by a chaotic basin boundary and the high-dimensional attractor formation through the merger of distinct chaotic sets. A future research direction that is motivated by these results is a theoretical treatment of high-dimensional attractors as those formed by chaotic repellers around distinct periodic orbits. In particular, it might be possible to adapt the periodic orbit expansions for chaotic repellers⁵⁰ to predict lifetime distributions such as those shown in Fig. 5.

One of the most striking features of hydrodynamic pilot-wave systems is the wavelike patterns that are reminiscent of the quantum wave functions which emerge in long-time statistics of the droplet position at high Me.^9,10,62 Already in Ref. 9, Harris et al. argued that the emergent statistics can be understood as a droplet’s transient visit of unstable (quasi-)periodic orbits. Our results are consistent with their insights and suggest new methods to understand the nature of deterministic dynamics underlying the emergent wavelike statistics. Specifically, lifetime measurements, such as those in Fig. 5, can be used to characterize individual neighborhoods in these experiments. Moreover, we would like to note that tuning the lifetime distribution of a particular droplet motion by varying a control parameter’s proximity to a bifurcation value might be relevant for quantum analogies since the adjustment of the lifetime distribution of some state directly influences the probability of observing the droplet in that state. In our case, we achieved this by varying Me, which as we moved farther away from the symmetry-increasing crisis bifurcation, resulted in shorter and shorter oval lifetimes, hence a lower probability of observing it.

This work was partially funded by the Institute of Science and Technology Austria Interdisciplinary Project Committee Grant “Pilot-Wave Hydrodynamics: Chaos and Quantum Analogies.”

The authors have no conflicts to disclose.

G.C. and B.S. contributed equally to this work.

George Choueiri: Data curation (equal); Methodology (equal); Software (equal); Validation (equal); Writing – review & editing (equal). Balachandra Suri: Data curation (equal); Software (equal); Visualization (equal); Writing – original draft (supporting). Jack Merrin: Methodology (equal); Resources (equal); Writing – review & editing (equal). Maksym Serbyn: Funding acquisition (supporting); Supervision (supporting); Writing – review & editing (lead). Björn Hof: Funding acquisition (equal); Resources (lead); Supervision (supporting); Writing – review & editing (equal). Nazmi Burak Budanur: Conceptualization (lead); Data curation (equal); Funding acquisition (equal); Investigation (equal); Software (equal); Supervision (equal); Visualization (equal); Writing – original draft (lead).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

The vibrating bath in the experiment is made from black anodized aluminum. The circular corral located at its center has a diameter $19.95±0.05mm$ and depth $6±0.05mm$⁠, as shown in Fig. 1. The corral is filled to a height $6.9±0.1mm$ with silicone oil (polydimethysiloxane), which has kinematic viscosity $ν=21.5×10−6m2$/s, density $ρ=953kg/m3$⁠, and surface tension $σ=20.8×10−3$ N/m at room temperature (⁠$T=21.6±0.1°$C). For this fluid layer depth, a thin $0.9±0.1$ mm overflow layer—that serves to dampen surface waves—is formed outside the corral. A transparent plastic lid placed on top of the aluminum bath shields the corral (and the droplet) from stray air currents.

The aluminum bath is mounted on an air-cooled electromagnetic shaker (Data Physics V55) and is leveled perpendicular to gravity. The shaker oscillates vertically at a frequency $f0=75Hz$⁠, when a sinusoidally varying voltage signal of the same frequency is input. Three piezoelectric sensors (PCB 352C65), mounted as shown in Fig. 1, measure the vertical (⁠$γ$⁠) and horizontal (⁠$γx,γy$⁠) accelerations of the vibrating aluminum bath. The setup is aligned such that $|γx|/|γ|,|γy|/|γ|≲0.01$⁠. For a given memory $Me$⁠, the target vertical acceleration $γMe$ is computed using Eq. (1), i.e., $γMe=γF(1−Me−1)$⁠. Here, $γF=4.32g$ is the experimentally measured critical acceleration (in units of $g=9.8m2$/s) for the onset of Faraday instability at $21.7°C$⁠. A feedback loop controls the amplitude of the sinusoidal voltage signal driving the shaker such that the measured vertical acceleration (⁠$γ$⁠) deviates from $γMe$ by less than $±0.1$%, i.e., $|γ−γMe|/γMe≤0.001$⁠.

To generate droplets of a desired size, silicone oil—filled in a syringe to a fixed height—was drained through a 33 gauge needle for a fixed time. The syringe was then touched against the vibrating bath to dislodge the droplet onto the fluid surface. The longer (shorter) the duration of draining, the larger (smaller) is the drop size. For a fixed drain-time, droplets generated using this technique varied in diameter by about $±$0.05 mm. All experimental runs reported in this study were performed with a single silicone oil droplet of diameter $0.85±0.05mm$⁠. Images of the bouncing droplet were recorded at intervals $Δt≈57$ ms, using a CMOS camera (Basler acA2000-165um) mounted above the bath [cf. Fig. 1(b)]. In pixel units, the diameters of the corral and the droplet are $468±1$ and $20±1$⁠, respectively. To track the position of the droplet in real-time, a gaussian filter (with a 6 pixel standard deviation) was applied to each image and the location of the brightest pixel—approximating the coordinates (⁠$xc,yc$⁠) of the droplet center—was measured. The time series $xc(t),yc(t)$ were then interpolated onto a temporal grid with spacing $Δt=5.7$ ms, using a cubic spline interpolation. Instantaneous droplet velocities $vx(t),vy(t)$ were then computed by computing derivatives of the spline interpolation. An example, Fig. 7 shows droplet trajectories reconstructed by overlaying successive images, each approximately 150 ms apart in time.

The horizontal speed $‖v‖$ of the droplet remains fairly constant (⁠$≈0.06$⁠) for the dynamical regimes explored in this article. To demonstrate this, Fig. 8 shows that the probability density function of $‖v‖$ for representative values of $Me∈[20,33.64]$ is normalized such that $∫Pd‖v‖=1$ in each case. Recall that symmetry reduction in the velocity plane requires that the speed of the droplet does not vanish. Indeed, the probability of the droplet speed $‖v‖≤0.01$ is smaller than $10−3$ for $Me≥32.1$⁠. Even though this is an insignificant fraction, such a decrease in droplet speed is possibly due to the droplet (very rarely) approaching the boundary of the inner corral and subsequently bouncing back.

The sequence of bifurcations (Figs. 2 and 3) was identified using an experimental run, where memory $Me∈[20,33]$ was increased (decreased) in steps of $ΔMe≈0.36$ (⁠$ΔMe≈−0.37$⁠) and the droplet position was tracked for a duration of 1800 s at each $Me$⁠. The lifetime distributions (cf. Fig. 5) were estimated from separate runs, each approximately $9×104$ s long. The ambient room temperature (measured using a PT-100 probe placed beside the shaker) over the duration of these experimental runs was fairly constant (⁠$20.6±0.1°C$⁠). Nevertheless, we measured the critical memory $Mec$ corresponding to the crisis bifurcation before and after each experimental run. $Mec$ was reproducible to within $ΔMe≈±0.15$ for different runs, which suggests that the temperature of the silicone oil, its physical properties, and, consequently, $γF$ do not vary significantly across the various experiments runs. Last, although all experimental runs reported in this article correspond to a single droplet, the repeatability of results was validated using (at least five) different experiments performed with different droplets and replacing the silicone oil in the bath each time.

The survival probabilities in Figs. 5(b) and 5(c), as mentioned in the main text, were estimated from window-averaged angular momentum $⟨L⟩Tw$ of the droplet. Here, $Tw$ is the average characteristic time-scale for tracing the lemniscates and ovals at each value of Me. $Tw$ was estimated by computing the Fourier transform of the time series of $L$ for the entire duration (⁠$≈9×104$ s) of each experimental run. Table I lists the values of $Tw$ for lemniscates and ovals at each $Me$⁠.