Damped-driven system of bouncing droplets leading to deterministic diffusive behavior

Damped-driven systems in physical and biological sciences often lead to very complex phenomena, which nevertheless can be modeled as mathematically tractable deterministic dynamical systems [1, 2, 3, 4, 5, 6, 7, 8]. It has also been shown that the interplay between damping and driving can induce diffusive behavior such as in the driven Josephson junction [9, 10]. Interestingly, the models for the driven Josephson junction [9, 10] are deterministic, whereas Einstein derived a theory for Brownian motion as a stochastic process [11]. Other examples of diffusive behavior induced through chaotic dynamical systems are treated in reviews by Geisel [12] and by Artuso and Cvitanović [13]. While all chaotic systems are seemingly random [14], deterministically diffusive systems experience Brownian-like behavior [12] where the mean squared displacement of the diffusive variable increases monotonically with time and admits a Gaussian distribution (during normal diffusion) for the variable from the initial point.

Since the seminal work of Couder and coworkers [15, 16], bouncing droplets have demonstrated fundamental physical phenomena at the macro-scale as detailed in the reviews by Bush and coworkers [17, 18, 19]. In this system we drop a silicon oil droplet onto a silicon oil bath sitting on top of a vertically vibrating shaker. It has been known since the 1970s that such a droplet can bounce on the bath without coalescing [20]. If the force on the shaker is increased the droplet starts moving horizontally (and is thus dubbed a walking droplet) [15, 16]. In more recent years these bouncing droplets have shown quite complex self-organizing behavior such as forming crystal-like structures [21, 22, 23]. If the forcing is increased furthermore still, Faraday waves form on the surface of the bath [24, 25, 26, 27, 28], upon which the droplet can bounce around chaotically.

While the detailed hydrodynamic models can often be quite complex, reduced dynamical systems models provide a framework to analytically study the solution space of droplet trajectories thereby facilitating rigorous analysis of long-time statistics of the trajectories. Stroboscopic models [29, 30, 31] were some of the first model reductions allowing rigorous dynamical systems analysis of long-time statistics. These models assume a closed form formula for the eigenmodes that contribute to the wavefield, and a continuous horizontal forcing on the droplet based on this wavefield. Furthermore reductions came from discrete dynamical models [32, 33, 34, 35, 36] and discrete time differential equations [37]. Gilet’s model of a walker in a linear confined geometry [32] considered discrete impacts with a simplified wavefield due to the confined eigenmodes. The model also only kept track of the droplet position at each impact and not the motion between impacts. The simplicity of the model allowed for detailed dynamical systems analysis of the bifurcations and chaotic strange attractors [33, 38, 39], and provided insight for the potential route to chaos. Gilet later developed a discrete dynamical model of walkers on a circular corral [34], which produced long-time statistics similar to that of experiments [40], but was more computationally efficient than the stroboscopic model. In [35] Rahman developed a standard map-like model for walkers on an annulus, which used the experiments of Filoux et al. [41, 42] as inspiration. Through this 1-dimensional model he showed that the walker becomes chaotic through a cascade of period doubling bifurcations, which had previously been observed and hypothesized in experiments [43]. A more detailed overview of the dynamical systems models and techniques used to study walking droplets can be found in the review of Rahman and Blackmore [44].

Recently, researchers have observed potentially diffusive behavior in systems of walking droplets. Tambasco et al. experimentally studied the dynamics of a droplet in 2-dimensional free space above the Faraday wave threshold, and calculated an effective diffusivity for a droplet during erratic walking [28]. Later, Hubert et al. developed an experiment-inspired stroboscope-like model for walkers below the Faraday wave threshold and studied how the effective diffusivity changes with the memory of the system [45, 46]. In more theoretical studies, Gilet observed diffusive behavior in discrete dynamical models below the Faraday wave threshold [32, 34] and Durey observed such behavior in Lorenz-type continuous dynamical systems [47]. Moreover, Durey et al. observed seemingly random walks in a stroboscopic model for speed oscillations [48] and Valani et al. observe diffusive behavior in the timeseries of their stroboscopic model for a hypothetical classical wave-particle entity (akin to walkers) [49]. Diffusive behavior has also been studied in several other works on walkers [50, 51, 52, 53, 54]. Heretofore, diffusive behavior had not been systematically studied for a bouncer. In this manuscript, we study experiments with bouncing droplets in an annulus above the Faraday wave threshold, and construct a hydrodynamic-kinematic model incorporating the evolution of the wavefield and both the horizontal and vertical dynamics of the droplet.

The remainder of the manuscript is organized as follows: Section 3 starts us off with a first principles model of the hydrodynamic-kinematic interactions of the droplet and the bath. We discuss the numerics in Sec. 4 and the theoretical results for the hydrodynamic-kinematic model in Sec. 5. We conclude with some final discussions in Sec. 6.

The schematic for the experimental set-up of Pucci [55], conducted at the MIT Applied Math Laboratory and privately communicated to the author by Pucci himself, are presented in Fig. 1(a,b). A circular container is filled with silicone oil with density $\rho=950$ kg/m³, viscosity $\nu=20.9$ cSt and surface tension $\sigma=20.6$ mN/m. The liquid bath has overall diameter $15.8$ cm and is divided in three regions with different liquid depths: an inner and an outer region (in white in Fig. 1(a)) in which only a thin liquid layer is present, and an annular channel (in black in Fig. 1(a)). The radius of the channel’s centerline is $R=6.63$ cm, the channel’s width is $w=7.5$ mm and its height $H=4.6$ mm. The container is initially filled with silicone oil up to a height $\gtrsim H$. Then some liquid is removed from the channel until the liquid level in the channel reaches the height $h=4.3$ mm. With this procedure, we create a tiny liquid meniscus with height $h_{m}\approx 0.3$ mm that will impede the drop to escape the channel, while the rest of the bath remains covered by a thin layer of oil (Fig. 1(b)). Liquid heights are measured with a dipping tip controlled by a micrometer with measurement error $\Delta h=0.03$ mm.

The bath is vertically driven by an electromagnetic shaker with acceleration $\Gamma(t)=\gamma\cos(2\pi f_{0}t)$, where the vibrational frequency is $f_{0}=80$ Hz. The vibration system ensures constant acceleration amplitude to within $\pm 0.002$ g and spatially uniform vertical vibration to within 0.1$\%$ [56]. The vibration is monitored by two accelerometers placed symmetrically with respect to the center of the container. The experimental uncertainty on the dimensionless acceleration $\gamma/\gamma_{F}$ is less than 0.2%, and primarily due to small variations of $\gamma_{F}$ during the course of the experiments, as result from slow drifts in bath temperature and the dependence of fluid viscosity on temperature. For $\gamma\geq\gamma_{F}\approx 4.2$ g, a quasi-1D standing surface wave with wavelength $\lambda_{F}=4.75$ mm, wave vector $\vec{k}_{F}$ tangential to the channel centerline and period $T_{F}=2/f_{0}$ appears in the annular channel. This is a Faraday wave, that is, a manifestation of the Faraday instability [24], and is consistent with the dispersion relation of gravity-capillary waves [57].

A droplet of the same silicone oil with diameter $0.56$ mm is generated via a droplet generator [58] and placed in the channel while $\gamma<\gamma_{F}$. The container is sealed with a transparent acrylic lid that provides isolation from ambient air currents, a necessary precaution for repeatable experiments [59]. The forcing acceleration is then increased above $\gamma_{F}$ and the droplet starts moving along the channel. The droplet remains at all time in the vicinity of the channel centerline, thus providing a good approximation for one-dimensional motion (Fig. 1(a)). The droplet motion is recorded from above with a CCD camera at 20 frames per second and tracked with an in-house algorithm.

We model the experiment in Sec. 2, a walking droplet on an annulus bouncing on a vibrating fluid bath, as a solid particle interacting inelastically with the waves created by the bath, and subsequently traveling through the air as a projectile. Since the Reynolds number ($Re$) for the droplet during flight will be in the intermediate regime ($1<Re<10^{3}$), there are some nontrivial considerations for the drag on the droplet, which Molacek and Bush [60, 61] analyzed extensively. Through their analysis it is evident that the correction terms to Stoke’s drag is sufficiently small, and does not alter the qualitative behavior of the projectile. Next, we assume the walker comes in contact with the wave at isolated impacts $n$ at a position $(x_{n},y_{n})$; that is, there is a single impact in a temporal neighborhood of the impact $t\in(t_{n}-\varepsilon,t_{n}+\varepsilon)$ for some $\varepsilon>0$. For the intermediate dynamics (the motion between impacts), let $\xi$ be the circumferential direction along the annulus and $\eta$ be the height, then our model becomes

	$\displaystyle\ddot{\xi}$	$\displaystyle=-\nu\dot{\xi};\qquad\xi(0)=x_{n},\quad\dot{\xi}(0)=v_{n}^{x}$		(1a)
	$\displaystyle\ddot{\eta}$	$\displaystyle=-g-\nu\dot{\eta};\qquad\eta(0)=y_{n},\quad\dot{\eta}(0)=v_{n}^{y},$		(1b)

where $\nu=6\pi\mu_{a}/(\rho R_{d}^{2})$ is the coefficient due to Stoke’s drag for a sphere of radius $R_{d}$ with a material density of $\rho$ through air with dynamic viscosity $\mu_{a}$. furthermore, $v_{n}^{x}$ and $v_{n}^{y}$ are the velocities, in the $x$ and $y$ directions respectively, immediately after the $n^{\text{th}}$ impact.

Since the intermediate dynamics only considers the time between impacts, we use the time $t$ for the time between impacts $n$ and $n+1$, and $T=\sum_{i=0}^{n}t_{i}+\Delta t;\qquad\Delta t\in(t_{n},t_{n+1})$ as the elapsed time for the evolution of the wavefield. In the experiments we observe that the wavefield, illustrated in Fig. 2b, is of the form

$$\begin{split}\Psi(T,x)=&A\cos(\pi fT)\cos\left(\frac{2\pi}{\lambda_{f}}x\right)\\ &+\frac{\gamma}{(2\pi f)^{2}}\cos\left(2\pi f[T-\varphi]\right).\end{split}$$

(2)

where $\gamma$ is the maximum bath acceleration, $f$ is the bath frequency, which implies a period of $T_{f}=2/f$, $\varphi$ is the bath vibration phase, and $\lambda_{f}$ is the wavelength.

Next, we nondimensionalize (1) using $\tau=ft/2=t/T_{f}$, $\hat{\xi}=\xi/\lambda_{f}$, and $\hat{\eta}=\eta/Y$ where the choice of $Y$ is $Y=g/(2\pi f)^{2}$ is equivalent to previous studies of McBennett and Harris[62] and Halev and Harris [63]. Then (1) becomes

	$\displaystyle\ddot{\hat{\xi}}$	$\displaystyle=-\nu^{*}\dot{\hat{\xi}};\qquad\hat{\xi}(0)=\hat{x}_{n},\quad\dot{\hat{\xi}}(0)=\hat{v}_{n}^{y}x$		(3)
	$\displaystyle\ddot{\hat{\eta}}$	$\displaystyle=-G-\nu^{*}\dot{\hat{\eta}};\qquad\hat{\eta}(0)=\hat{y}_{n},\quad\dot{\hat{\eta}}(0)=\hat{v}_{n}^{y},$		(4)

with

$$\Psi=A^{*}\cos(2\pi\hat{T})\cos(2\pi\hat{\xi})+\gamma^{*}\cos(2\pi[2\hat{T}-\varphi^{*}]),$$

(5)

where $\nu^{*}=2\nu/f$, $\varphi^{*}=f\varphi$, $G=16\pi^{2}$, $A^{*}=A(2\pi f)^{2}/g$, and $\gamma^{*}=\gamma/g$ are the dimensionless drag coefficient, phase, gravity, wave amplitude, and bath forcing amplitude, respectively. We observe that in our system $\nu^{*}\approx 0.1$, and as we will justify in the sequel, the size of $\nu^{*}$ will exacerbate pathological numerical artifacts that are avoided by considering an undamped model. We write the dimensionless system with the hats and asterisks removed as

	$\displaystyle\ddot{\xi}$	$\displaystyle=0;\qquad\xi(0)=x_{n},\quad\dot{\xi}(0)=v_{n}^{x}$		(6a)
	$\displaystyle\ddot{\eta}$	$\displaystyle=-G;\qquad\eta(0)=y_{n},\quad\dot{\eta}(0)=v_{n}^{y},$		(6b)
	$\displaystyle\Psi$	$\displaystyle=A\cos(2\pi T)\cos(2\pi\xi)+\gamma\cos(2\pi[2T-\varphi]).$		(6c)

Solving (6) yields

	$\displaystyle\xi$	$\displaystyle=v_{n}^{x}t+x_{n},$		(7a)
	$\displaystyle\eta$	$\displaystyle=-\frac{1}{2}Gt^{2}+v_{n}^{y}t+y_{n},$		(7b)

which explicitly determines the path of flight of the particle (Fig. 2a red dashed line); thereby mitigating the computational expense of solving a system of Ordinary Differential Equations (ODEs).

Now, given a particular point of intersection between the path of the projectile and wavefield, we can model the inelastic collision between the droplet and the bath. As an illustrative aid we employ Fig. 3. Suppose the droplet is already in flight with a velocity of $v_{\text{in}}$, just before the $n^{\text{th}}$ impact, approaching the point of impact at a slope of $d\eta/d\xi=(d\eta/dt)/(d\xi/dt)=\dot{\eta}/\dot{\xi}$ yielding an approach angle (with respect to the abscissa) of

$$\theta_{\text{in}}=\tan^{-1}\left(\dot{\eta}/\dot{\xi}\right)\in[-\pi/2,\pi/2].$$

(8)

furthermore, assume that at the point of impact the wavefield has a tangent line with a slope of $d\Psi/dx$, and hence an angle

$$\phi=\tan^{-1}\left(\frac{d\Psi}{dx}\right)\in[-\pi/2,\pi/2],$$

(9)

with respect to the abscissa. The angle between the tangent lines of the path and wavefield is

$$\alpha_{\text{in}}=(\pi+\theta_{\text{in}})-\phi.$$

(10)

This gives us the components of the approach velocity $v_{\text{in}}$ that are tangential ($v_{\text{in}}\cos(\alpha_{\text{in}})$) and normal ($v_{\text{in}}\sin(\alpha_{\text{in}})$) to the wavefield. In Sec. 4, we use the velocity components to compute the tangential and normal components of the coefficient of restitution, $C_{R}^{T}$ and $C_{R}^{N}$. Then the tangential, which has the opposite sign to that of the approach, and normal components of the launch velocity are

	$\displaystyle v_{n}^{T}$	$\displaystyle=-C_{R}^{T}v_{\text{in}}\cos(\alpha_{\text{in}}),\quad\text{and}$		(11a)
	$\displaystyle v_{n}^{N}$	$\displaystyle=C_{R}^{N}v_{\text{in}}\sin(\alpha_{\text{in}}),$		(11b)

yielding a launch angle of

$$\alpha_{\text{out}}=\tan^{-1}\left(\frac{v_{n}^{N}}{v_{n}^{T}}\right)\in[0,\pi],\quad\text{and}\quad\theta_{n}=\alpha_{\text{out}}+\phi,$$

(12)

with respect to the gradient of the wavefield for $\alpha_{\text{out}}$ and the abscissa for $\theta_{n}$. Finally, we arrive at the launch velocity components,

	$\displaystyle v_{n}^{x}$	$\displaystyle=\sqrt{\left(v_{n}^{T}\right)^{2}+\left(v_{n}^{N}\right)^{2}}\cos\theta_{n},$		(13a)
	$\displaystyle v_{n}^{y}$	$\displaystyle=\sqrt{\left(v_{n}^{T}\right)^{2}+\left(v_{n}^{N}\right)^{2}}\sin\theta_{n},$		(13b)

required for (6).

Since the ODEs in (6) have explicit solutions, our numerical methods are dedicated to tracking the evolution of the wavefield with respect to the path of the projectile, resolving the effects of each impact, and calculating various statistics to compare with experiments. The event-driven nature of the phenomenon [64, 65, 66], and inherent discrepancies in sampling between experimental data and numerics, have necessitated novel computational remedies to oft-overlooked simulation problems. We approach the computational issues from a physics-inspired point of view with the aim of preserving as much of the kinematic realism as possible.

We first need to have a method to find the point of impact as the path continues indefinitely as shown in Fig. 2. To do this we need to calculate the point of intersection between the flight path and the wavefield, which would be easy if the wavefield stood still. What makes it even more challenging is the evolution of the wavefield shape in addition to the vertical oscillations. So, the impact involves three actions: the flight of the particle, the vertical motion of the table, and the evolution of the shape of the wavefield. Now, we have the quintessential “not knowing where the edge is until passing it” problem, and off the edge we go by letting the particle go through the bath. As soon as the particle height is less than the bath height, $\eta(t_{i+1})<\Psi(t_{i+1},\xi(t_{i+1}))$, at the horizontal location of the particle, $\xi(t_{i+1})$, the code is paused. Then it is necessary to work backwards to find the intersection; however, since so many aspects of the system are changing around that impact point, going one time step back to $t_{i}$ would cause the particle, $(\xi(t_{i}),\eta(t_{i}))$, to be too far from the point of intersection. One remedy is to use a finer global time step, but this would prohibitively increase the computational expense, and instead we choose to use a bisection scheme in the domain $(t_{i},t_{i+1})$. The scheme produces the time $t_{*}$ such that $\eta(t_{*})-\Psi(t_{*},\xi(t_{*}))<10^{-4}$, which implies $\eta(t_{*})>\Psi(t_{*},\xi(t_{*}))$; a necessary condition to keep the simulation physically realistic. Finally, the launch variables can be calculated as shown in Sec. 3.

The issue of chattering arises when the velocity of the bath is asymptotically equivalent to the vertical velocity of the droplet at impact; that is, $d\Psi/dt\sim d\eta/dt$ at $t=t_{*}$. Then the droplet can numerically go through the bath more than once within one time step, no matter how small we make the time step, which creates a non-isolated intersection; i.e., there is some $t_{**}\neq t_{*}$ such that $t_{i}\leq t_{*}\leq t_{**}\leq t_{i+1}=t_{i}+\varepsilon$ for all $\varepsilon$ larger than some computationally prohibitive time step, and often occurs for all $\varepsilon$ larger than machine accuracy. That is, numerically, the algorithm records repeated intersections when no such intersection is present physically, which causes the algorithm to effectively halt. Fortunately, the physical observations offer guidance. In the laboratory, the impacts are not instantaneous, and therefore a close encounter with the bath may not behave like an impact in our computational setup. Numerically, we let the droplet shadow the bath for a small period of time until the droplet velocity and bath velocity diverge. Furthermore, in the previous paragraph we calculated our impact with respect to the center of the droplet. However, in reality the contact is not with the center, but rather the surface of the droplet, and hence we do the same calculations except with $\Psi+R$ in place of $\Psi$ where $R$ is the radius of the droplet. The actual physics is more nuanced as the radial separation should be calculated in the direction of the flight path and we do not consider deformations, however this would add considerable numerical complexity with minimal additional insight about the longtime statistics. Now, if the newly calculated vertical velocity (just after impact) is less than the bath velocity, $v_{n}^{y}<d\Psi/dt|_{t=t_{*}}$, which implies that there will be a non-isolated numerical intersection, we let the velocity and position shadow that of the bath, $v_{n}^{y}=d\Psi/dt|_{t=t_{*}}$ and $y_{n}=\Psi(t_{*},\xi(t_{*}))+R$.

Next, due to the inelastic nature of the collision between the droplet and the bath, we calculate the coefficients of restitution, $C_{R}^{T}$ and $C_{R}^{N}$, as functions of the Weber number, $We=\rho R_{d}\left(v_{\text{in}}^{N}\right)^{2}/\sigma$, where $\sigma$ is the surface tension of the droplet and $v_{\text{in}}^{N}$ is the droplet velocity in the normal direction with respect to the wave surface. In the experiments of Moláček and Bush [60, 61], the experimental coefficients of restitution were plotted against the Weber number. We fit various polynomials, with respect to the log of the Weber number, to the data points (Fig. 4) and choose a linear fit for $C_{R}^{N}$ and a quartic fit for $C_{R}^{T}$,

	$\displaystyle C_{R}^{N}\approx-0.0779\log_{10}We$	$\displaystyle+0.2198,$		(14a)
	$\displaystyle\begin{split}C_{R}^{T}\approx-0.0044(\log_{10}We)^{4}&-0.0454(\log_{10}We)^{3}\\ &-0.1906(\log_{10}We)^{2}\\ &-0.3948\log_{10}We+0.6209.\end{split}$			(14b)

Since the polynomial approximations do not asymptote, and it is shown that $0.1\lesssim We\lesssim 1$ [60, 61] (and also in Fig. 4(c)), we restrict the ranges to $0.2\leq C_{R}^{N}\leq 0.34$ and $0.26\leq C_{R}^{T}\leq 0.8$. Interestingly, the motion of the particle depends more on the tangential coefficient of restitution than on the normal coefficient.

In this section we use the parameter values listed in Table 1 for the simulations unless otherwise stated. We first present qualitative similarities between the simulations and experiments. Then we show statistical results to illustrate the amount of the dynamics that is captured by the model. Finally, we use the model to discuss subtleties in the dynamics that are inaccessible to experiments.

In addition to the fundamental parameters in Table 1, we also must resolve the wave amplitude, $A$, which depends on the amplitude of forcing acceleration, $\gamma$. Since our maximum acceleration is $\gamma=1.047\gamma_{f}$, it is convenient to write the coupling factor [26] as $\Gamma\approx 1/2\sqrt{1.047-1}$, which allows the model to capture the qualitative features observed in experiments.

In both experiments and numerical simulations we identify two regimes of droplet motion, which depend on the dimensionless acceleration $\gamma/\gamma_{F}$ and thus on the wave amplitude $A$ (Fig. 5). At low amplitude, the droplet executes jumps with length $\sim\lambda_{F}$ between sites spaced by multiples of $0.5\lambda_{F}$. The droplet lingers around each site for a typical time $T_{\rm trap}\gg 10T_{F}$ within a region with width $<0.5\lambda_{F}$. Each site is presumably a wave trough within which the droplet is transiently confined in a regime of chaotic bouncing. We call this regime “transient confinement” (Fig. 5(a,b)), which is reminiscent of Lévy flights [67]. At higher amplitude the droplet motion is erratic and each site is visited during a typical time $T_{\rm err}\lesssim 10T_{F}$ (Fig. 5(c,d)), which is reminiscent of deterministic diffusion [12].

The walker appears to prefer velocities near zero and chooses higher speeds with decreasing probability as illustrated in Fig. 6(a, b), where the trajectory from Fig. 5(c, d) is one of the trajectories from Fig. 6(a, b). This is very similar to particles experiencing Brownian motion. Moreover, Fig. 6(c, d) shows that the distribution of the velocities for $\gamma=1.047\gamma_{f}$ appear to be Gaussian, which is indicative of Brownian motion.

In order to quantify the behavior of the system in terms of diffusion we compute the quantity

$$d(t)=\frac{T_{F}}{\lambda_{F}^{2}}\frac{\sigma^{2}(t)}{2t},$$

(15)

where $\sigma^{2}(t)=\sum_{i=1}^{N}|x_{i}-\bar{x}|^{2}/(N-1)$, with $N$ the number of points and $\bar{x}$ the mean position. Experimentally, this is measured by splitting a video with duration $\geq 60$ s in samples with duration $2$ s and computing the variance across the samples. We observe that $d(t)$ increases with time to saturate around a constant mean value (Fig. 7(a)). Using the dimensionless model in (7), $d(t)$ is computed numerically by first applying the model to $300$ equally-spaced initial points within an interval corresponding to one wavelength, $x_{0}\in[0,1)$, for an elapsed time corresponding to $80$ bath oscillations, $T=80$. We then compute the variance $\sigma^{2}$ for each run and divide by $2T$ to obtain $d(t)$ as in (15).

The temporal evolution and saturation value of $d(t)$ in experiments and numerical simulations are similar. We call

$$D=\lim_{t\to\infty}d(t)$$

(16)

the diffusion coefficient, that we compute in both experiments and numerical simulations as the average value of $d(t)$ over the last $20$ time steps.

We proceed by varying the forcing acceleration $\gamma/\gamma_{F}$ in both experiments and simulations and computing $D(\gamma/\gamma_{F})$. The results are presented in Fig. 8. In both experiments and numerical simulations we observe an overall increase of the diffusion coefficient with the forcing acceleration, in both the transient confinement and the erratic motion regimes. As an exception to this general trend, we observe $D$ to decrease at the edge between the two regimes. In experiments, this is due to long-lasting confinement in which the droplet presumably bounces only once over four bath vibration periods. In the model when the particle starts moving, its velocity is nearly constant yielding an appreciable $\sigma(t)$, however when the motion starts to become erratic (at the edge between the two regimes), the droplet velocity has equal probability of being positive and negative, effectively reducing $\sigma(t)$. As $\gamma/\gamma_{F}$ is increased furthermore, the droplet velocity is large enough that once again the deviation from most initial points is appreciable, which contributes to an increasing $\sigma(t)$ even though there will still be some initial points that do not contribute.

We often find chaotic dynamics, induced by simple interactions, giving rise to quite complex phenomena. This is particularly prevalent in deterministically diffusive systems as described in a review by Geisel [12] and by Artuso and Cvitanović [13]. This is also true for a larger class of damped-driven systems [68, 69, 8], which can be studied using an energy gain-loss formulation.

In this work we studied the motion of a droplet bouncing on a standing wave. As the wave amplitude is increased the motion of the droplet is observed to go through several bifurcations until it eventually becomes erratic. At this erratic regime the distribution of droplet velocities tends to a Gaussian, which indicates potentially diffusive behavior. In fact, the effective diffusivity can be computed using (15), and we observed that the diffusivity is not a monotonic function of the bath forcing. We then presented a detailed derivation of the hydrodynamic-kinematic model of walking droplets on an annular cavity with bath forcing above the Faraday wave threshold. Then we went over the numerical methods used to simulate the circumferential droplet motion leading to long-time statistics. The numerical methods require unique solutions to mitigate the chattering effect due to the event-driven nature of the problem. These mitigation strategies were inspired by the physical process being investigated thereby preserving much of the physical agreement with experiments. The theoretical results are studied in more detail in Sec. 5. We showed that there are three dynamical regimes in the model: trapped (minor excitation around a fixed velocity), transient confinement (Levy-like flights [67] where the droplet is trapped for variable amounts of time), and the diffusive state (droplet motion is seemingly random with chaotically changing velocity). Notice that the trapped state is just a special case of the transient confinement state. We show that the qualitative behavior of the model is very close to that of experiments. This modeling framework can also be used for fast computation of long-time statistical behavior when the vertical motion of the particle is necessary, which is something a stroboscopic model may overlook.

Interestingly, our hydrodynamic-kinematic model exhibits properties, such as Brownian-like motion and Levy flights, similar to those of other deterministically induced diffusive systems [12]. In addition to future dynamical systems studies, the present experimental and modeling framework can be used to conduct a more detailed 2-dimensional investigation of droplets above the Faraday wave threshold. While Tambasco et al. [28] also observed diffusive behavior above the Faraday wave threshold, they use a much larger droplet ($20\%$ larger at the lower end). By using a smaller droplet we are able to ignore deformation and drag during the modeling process. Furthermore, we foresee a more detailed investigation of diffusive properties for the motion of the droplet than done here yielding far more nuanced and unexpected dynamics.

Finally, there are also intriguing pedagogical possibilities. Both the experiment and model can be reproduced to be used in a teaching lab as a visual demonstration of a system exhibiting diffusion. The experimental setup can be reproduced using an inexpensive audio speaker and a wave generator on a mobile phone. Using 3-dimensional printing, a dish with an annular cavity to hold the bath can be produced. The dish would then be attached to the speaker (either by glue or a 3-D printed adapter). When the waves are sent through the mobile phone to the speaker the bath will be accelerated vertically creating the environment for walking droplets. Moreover, the full hydrodynamic-kinematic model of Sec. 3 can be easily simulated with the available computationally efficient MATLAB codes. To study a variety of parameters one can modify the easy to follow preamble of the code before running it as usual. Most importantly, the model is accessible to a wide range of individuals, and we hope that it will be particularly useful as a pedagogical example of a damped-driven system for undergraduate students and even high school students.

First and foremost, this work could not have been accomplished without the experimental data from Giuseppe Pucci, to whom A.R. is very grateful. A.R. also thanks John Bush for the use of the facilities in MIT’s Applied Math Lab. In addition, this work has benefited from fruitful discussions with Giuseppe Pucci, Dan Harris, and J. Nathan Kutz. Finally, A.R. appreciates the support of the Amath department at UW.