Dynamics, emergent statistics, and the mean-pilot-wave potential of walking droplets

A millimetric droplet may bounce on the surface of a vertically vibrating bath of the same fluid; the thin air layer separating the droplet from the bath during impact prevents coalescence.^1,2 Each impact excites a field of temporally decaying Faraday waves, whose longevity depends on the reduced acceleration $Γ=Aω02/g$⁠, where $A$ is the shaking amplitude, $ω0/(2π)$ is the frequency, and $g$ is the gravitational acceleration. As $Γ$ increases, the bouncing may destabilize to horizontal “walking” across the bath, whereby the droplet is propelled at each impact by the slope of its associated Faraday wave field³ [see Fig. 1(a)]. The decay time of the Faraday waves increases with $Γ$ for $Γ<ΓF$⁠, where the Faraday threshold $ΓF$ is the critical vibrational acceleration at which Faraday waves arise in the absence of a droplet. This decay time results in a “path-memory” of previous impacts, where the memory timescale is inversely proportional to the proximity of the Faraday threshold $ΓF$⁠.⁴ The resulting dynamics are similar in many respects to the pilot-wave dynamics envisaged by de Broglie as a physical framework for understanding quantum mechanics.⁵ 

The pilot-wave dynamics of this hydrodynamic system gives rise to quantumlike features in a number of settings, and so has prompted the investigation of several hydrodynamic quantum analogs.^6–10 The Faraday wavelength $λF$ plays a fundamental role in all of the hydrodynamic quantum analogs, imposing a lengthscale on the interaction between droplets, yielding a discrete set of quantized states for orbiting pairs,^3,11–15 promenading pairs,^15–17 and multi-droplet strings.¹⁸ When a walker is confined to a corral, a wavelike statistical pattern emerges.^19,20 A recent study has shown that the statistical wave form is similar to the time-averaged pilot-wave,²¹ but a quantitative relationship between the two was not found. Deducing such a relationship represents one of the key contributions of our study.

Further quantum analogies arise when the droplet is subject to either a central or a Coriolis force, where the latter is realized experimentally in a rotating bath. In both cases, the Faraday wavelength imposes a radial quantization of circular orbits at high wave memory,^22–24 whose stability have been analyzed theoretically.^15,25,26 As the circular orbits destabilize, a new family of stable exotic orbits emerges, revealing a range of extremely rich dynamics.^24,27 In particular, the orbits obtained under a central force exhibit a double quantization in their mean radius and angular momentum, yielding a remarkable analogy to quantum mechanics. The radial quantization may be rationalized in terms of the energy minimization of the mean Faraday wave field, whose form is determined by the orbital symmetry of each eigenstate.²⁴ In the chaotic regime arising at high vibrational forcing, a complicated switching process arises between the system’s underlying orbital states.²⁸ Statistical techniques have demonstrated that the double quantization is still present in the droplet’s chaotic dynamics.^15,29

The tendency of the walker system to self-organize into quantized dynamical states was demonstrated by Perrard et al.^24,28 and Labousse et al.³⁰ The conceptual value of decomposing the instantaneous pilot-wave field into its mean and fluctuating components was further stressed by Labousse.³¹ The merit of this decomposition in connecting the dynamics and statistics of pilot-wave systems is demonstrated here through consideration of a relatively simple geometry.

The complex structure of the exotic (non-circular) orbits has to date prohibited a comprehensive theoretical investigation of their dynamics in the periodic and chaotic regimes. Such a study is likely to shed new light on the quantumlike behavior and the role of the mean wave field in the long path-memory limit. To develop the techniques required for such an analysis, we focus this work on the dynamics in a harmonic potential where the droplet motion is restricted to a line and accompanied by a two-dimensional wave field [Fig. 1(a)]. This system exhibits extremely rich dynamics and analogies to quantum mechanics, whilst remaining simple enough to form the basis of a theoretical investigation that provides a foundational mathematical framework for future studies of more geometrically complex systems.

In the classical harmonic oscillator $mx″(t)+κx(t)=0$ with spring constant $κ$⁠, a particle of mass $m$ enters into simple harmonic motion with fixed frequency $ω=κ/m$⁠. The energy of the particle varies continuously with the initial conditions, and the motion is entirely deterministic. Conversely, in quantum mechanics, the particle energy $E$ is quantized with $En=ℏω(n+1/2)$ (where $ℏ$ is the reduced Planck’s constant and $n∈N$⁠), where for each energy level there is an associated probability distribution for the particle’s position. In what follows, we will demonstrate that the dynamics of the hydrodynamic pilot-wave system vary from classical to quantumlike, depending on the relative magnitudes of the wave and central forces. At low wave amplitude, the balance of wave and drag forces yields a stable limit cycle, whose oscillation amplitude and period vary continuously with the spring constant. When the waves dominate, the surviving pilot-wave from previous crossings of the bath causes significant variations in the droplet velocity, yielding quantized droplet range and wavelike statistics for the droplet position. The quantization length is $λF/2$⁠, and our study reveals that the Faraday wavelength also plays a pivotal role in the chaotic dynamics emerging near the Faraday threshold.

We herein apply the model of Durey and Milewski¹⁵ to elucidate the emergent quantizations, wavelike statistics, and the role of the mean wave field in the system’s periodic and chaotic dynamics. In Sec. II B, we prove that the droplet’s stationary probability distribution is related to the mean pilot-wave field via a convolution with the wave field of a bouncing droplet. In Sec. III, we extend the methods of earlier work to analyse the amplitude and stability of periodic oscillations, where we see the onset of quantization and wavelike statistics. In the limit of $Γ→ΓF$⁠, periodic wave-trapped solutions arise in which the droplet’s oscillatory motion persists even in the absence of an external force (⁠ $κ=0$⁠) and the mean wave field acts as an effective potential (Sec. III C). In Sec. IV A, we demonstrate that this system exhibits the Ruelle-Takens-Newhouse route to chaos.^32,33 At extremely high memory (as considered in Sec. IV B), the wave field dominates the droplet dynamics, yielding a short-timescale droplet motion similar to a random walk, and a long-timescale behavior influenced by an effective potential prescribed by the mean wave field. By detrending the long-timescale behavior induced by slow variations in the Faraday wave field, we see the emergence of pronounced wavelike statistics whose peaks are determined by the random walk dynamics.

The dynamics of this system are depicted by the schematic diagram in Fig. 1(b). We assume that the droplet and bath are in periodic subharmonic resonance (as observed in experiments over a broad parameter regime³⁴), and we model the impacts as instantaneous and localized at a point. This approximation is reasonable for describing short impacts with a small droplet, which we model as a rigid sphere. A full derivation of this model can be found in Ref. 15.

The semi-infinite fluid bath is governed by linear quasi-potential flow, which includes weak dissipative effects at high Reynolds number.^35–37 The harmonic velocity potential $ϕ$ and wave perturbation $η$ couple with the prescribed impact pressure $PD(x,t)=f(t)δ[x−X(t)]$⁠, where $x=(x,y)$ is the position on the fluid surface and $X(t)$ is the horizontal droplet position. For instantaneous impacts with subharmonic vertical motion, we require $f(t)=mg∑n=0∞δ(t/T−n),$ where $T=4π/ω0$ is the Faraday period and $m$ is the droplet mass.³⁸ The vibrating frame of reference introduces the effective gravity $gΓ(t)=g[1−Γcos⁡(ω0t+β)]$⁠, where $β$ denotes the droplet’s impact phase.

Following the model of Moláček and Bush,³⁸ the horizontal droplet position is governed by

with parameters given in Table I. During flight $(f=0)$⁠, inertia is balanced by the horizontal central force and Stokes’ drag with coefficient $νp=6πR0μair$⁠. During impact $(f>0)$⁠, the reaction force imparts a (linearized) kick to the droplet, which is countered by skidding friction characterized by the dimensionless drag coefficient $c>0$⁠,³⁸ whose value is discussed below.

Henceforth, we describe the dynamics in terms of dimensionless variables, where we scale lengths with the Faraday wavelength $λF=0.51cm$⁠, time $t$ with the subharmonic bouncing period $T=4π/ω0$⁠, force $f$ with $f0=mg$⁠, and pressure $PD$ with $P0=f0/λF2$⁠. This yields the following dimensionless parameters:

Typical parameter values from Table I give reciprocal Reynolds number $ϵ≈0.019$⁠, Bond number $B≈0.102$⁠, $G≈1.201$⁠, $M≈0.0017$⁠, $R≈0.075$⁠, $ν~p≈0.01$⁠, and vibration number $Ω≡4πR3/B=0.8$⁠.³⁹ The dimensionless potential strength $κ~≥0$ is a free parameter of both the model and experiments, with $10−3≲κ~≲10−1$⁠. The dynamics are largely insensitive to changes in the skidding friction $c$ and impact phase $β$⁠; thus, we fix $c=0.33$ and $β=π$⁠.¹⁵ 

To reduce the fluid system from partial to ordinary differential equations, we spectrally decompose $ϕ$ and $η$ in the horizontal plane. The simple “Dirichlet-to-Neumann” map for $ϕ$ under this decomposition allows us to eliminate $ϕ$ in favour of $η$⁠, which we express as

with orthogonal basis functions $Φm(x;k)≡Jm(kr)eimθ$⁠, where $x=(r,θ)$ in polar coordinates and $i$ is the imaginary unit. As $η$ is real and $Jm(z)=(−1)mJ−m(z)$ for all $m∈Z$⁠, the complex coefficients $am$ satisfy the reality condition $a−m=(−1)mam∗$ for all $m$⁠, where $∗$ denotes the complex conjugate. This basis decomposition yields a system of inhomogeneous damped Mathieu equations for the wave amplitudes $am$⁠, where the inhomogeneity arises from the instantaneous forcing at impact. Assuming $X(t)$ and $η(⋅,t)$ are continuous across impacts, we obtain nonlinear jumps in $X′$ and $ηt$ at impact times $t=tn≡n$⁠, which appear in (6) and (7) below.

During flight (⁠ $t≠tn$⁠), the wave and droplet dynamics decouple and evolve according to

where $γk≡2ϵk2$⁠, $ωk2≡Gk+Bk3$⁠, and $ωg,k2≡Gk$⁠. As (3) is of Floquet form with period $1/2$⁠, we evolve the system between impacts (⁠ $tn−1+↦tn−$⁠) using principal fundamental matrices $Mk(Γ)∈R2×2$ for the wave amplitudes [Eq. (3)] and $F(κ~)∈R4×4$ for the droplet dynamics (4). This yields the 3-stage map from $tn−1+↦tn+$⁠:

1. Compute the next droplet position and the pre-impact velocity

   $X(tn)X′(tn−)=F(κ~)X(tn−1)X′(tn−1+).$

    

   (5)
2. Update the wave amplitudes $∀k>0$ and $∀m∈Z$ 

   $am(tn;k)am′(tn+;k)=Mk(Γ)am(tn−1;k)am′(tn−1+;k)−pkΦm∗[X(tn);k]01.$

    

   (6)
3. Apply the droplet jump conditions

   $X′(tn+)=(1−F)X′(tn−)−H∇η[X(tn),tn],$

    

   (7)

where $F=1−exp⁡(−cGR/B)$⁠, $H=(F/c)B/R$⁠, and $pk=kMG/(2π)$⁠. Equations (6) and (7) couple through $∇η$ given by (2).

The wave “memory” $Me$ is defined as the timescale over which the Faraday waves decay, which is a proxy for the number of past impacts that influence the current dynamics.⁴ This appears naturally from the eigenvalues of $Mk(Γ)$⁠, which we write as $exp⁡(−s1)$ and $exp⁡(−s2)$ for $si=si(k,Γ)∈C$⁠, where $0≤Re(s1)≤Re(s2)$⁠. The dominant exponent $s1(k,Γ)$ is real and positive in a neighborhood of $(kF,ΓF)$⁠, with $s1(kF,ΓF)=0$⁠. For $Γ<ΓF$⁠, we thus define

as $Γ→ΓF$⁠, where $Td(Γ)∼0.6$⁠.^15,36 While this parameter diverges as $Γ→ΓF$⁠, we note that the description of the wave field in terms of linear Faraday waves also breaks down in this limit, where nonlinear effects are expected to become significant.

To implement the model given by Eqs. d6(5)–(7), we make an appropriate discretisation of the wavenumbers $k$ and truncate the Bessel modes $m$⁠, as detailed in Ref. 15. The diagonal structure in $k$ and $m$ allows for simulations at typically 1000 impacts per second⁴⁰ (this is 25 times faster than the experimental timescale). By using the methods developed in Sec. III, the discrete-time formulation d6(5)–(7) also allows for efficient computation of the system’s periodic states with linear stability analysis.

Previous investigations into the long-time dynamics of this hydrodynamic pilot-wave system have focused primarily on the statistical distribution of the droplet position $μ(x)$⁠, rather than considering the mean pilot-wave $η¯(x)$ at impact, as defined by

A recent study of walker motion in corrals pointed out that the two take a similar form;²¹ however, a quantitative relation between the two was not deduced. We proceed by proving that (in an unbounded domain) these two quantities are in fact related via the convolution $η¯=ηB∗μ$⁠. Here, $ηB(x)$ is the axisymmetric wave field of a stationary bouncing droplet at impact, which is the wave field generated by infinitely many periodic impacts at $x=0$ for a given $Γ<ΓF$ (see Fig. 2). We note that although we focus on the dynamics in the walking regime (in which the bouncing state is unstable), the associated bouncer wave field $ηB$ still plays a pivotal role in the long-time statistics.

The physical intuition behind this convolution result is as follows. For stationary dynamics, each point $x$ in the domain (within the support of $μ$⁠) is visited infinitely many times. If all points were visited equally, each would thus contribute equally to the mean wave field, in the amount of $ηB(x)$⁠. Since they are not visited equally, the contribution of each point $ηB(x)$ must be weighted by $μ(x)$⁠. Our result not only combines three key quantities of this pilot-wave system but is also valid for periodic motion and ergodic dynamics, such as in the chaotic regime. As will be seen, this convolution result is particularly useful for elucidating the dynamics at high memory, including periodic wave-trapped behavior (Sec. III C) and chaotic dynamics near the Faraday threshold (Sec. IV C).

We have proved similar convolution results for the models of Fort et al.²² and Oza et al.,⁴² where different modeling assumptions were made on the wave field dynamics and the droplet-wave coupling. In fact, we generalize the convolution relationship (10) to a wider pilot-wave framework in  Appendix A, which includes the pilot-wave dynamics in a confined geometry. In this more general case, the integral kernel $ηB$ is replaced by a function that no longer exhibits translational invariance.

The result in Theorem 1 rests on the assumptions that a stationary distribution exists and that the pilot-wave dynamics are ergodic. It has been observed experimentally that when the droplet’s motion is confined (by a harmonic potential²⁴ or the boundary walls of a corral^19,21), a stationary distribution may emerge. The ergodicity assumption is more delicate. It has been observed in the one-dimensional tunneling pilot-wave model of Nachbin et al.⁸ that several chaotic trajectories with different initial conditions had the same statistical properties as a single longer run, suggesting that the process is indeed ergodic in that particular configuration. We note, however, that when multiple stable states exist (such as in the case of hysteresis), the long-time behavior may depend on the initialisation of the pilot-wave system, rendering the ergodicity assumption invalid.

To overcome this difficulty, we prove an analogous result to Theorem 1 valid when the pilot-wave dynamics are periodic for all time, namely, $X(tn+Q)=X(tn)$ for all $n$ and some finite $Q∈N$⁠. This corollary does not require any assumptions about the existence or uniqueness of a stationary distribution, nor does it require the ergodic hypothesis.

Henceforth, we consider the case where the droplet motion is confined to a line. For $μ(x)=ρX(x)δ(y)$ [where $x=(x,y)$], Theorem 1 and Corollary 1 both simplify to $η¯0(x)=(ρX∗ηB)(x)$⁠, where $η¯0≡η¯|y=0$ is the mean wave field along the $x$-axis. We demonstrate in  Appendix B that in the case where the period $Q→∞$⁠, the result to Corollary 1 remains robust, even when the probability distribution $ρX(x)$ is approximated by a histogram.

We seek periodic solutions to the nonlinear discrete-time map d6(5)–(7) with motion restricted to the $x$-axis, so $X(t)≡[X(t),0]$ and $am∈R$ for all $m$⁠. For notational convenience, we denote

where $Gn=∂xη0(Xn,tn)$ is the wave field gradient along the $x$-axis at impact $n$⁠.

For any given $(Γ,κ~)$⁠, the frequency of the periodic oscillation is generally incommensurate with the Faraday frequency, which complicates the analysis for our discrete-time system. To resolve this, we exploit continuity of the parameter space to seek a subset of solutions where the oscillation period $P$ satisfies $P=φN$ (⁠ $N∈N$ and $φ∈Q$⁠) for a given $Γ$⁠, and solve for $κ~$ (it should be noted that in this case, there is a relationship between the oscillation period $P∈Q$ and the number of impacts $Q∈N$ such that $Xn=Xn+Q$ for all $n$⁠). Typically, $φ=2$ is sufficient to resolve the solution curve, which corresponds to the droplet crossing the bath once after $N$ impacts. This case yields reflection conditions (for all $m∈Z$ and $k>0$⁠)

For given $Γ$ and $N$⁠, we use a Newton method to compute the periodic states for $(N+1)$ unknowns $θ=(X0,κ~,G1,…,GN−1)$⁠, with the details given in  Appendix C. We exploit continuity of the solution branch by using as an initial guess a converged solution along the same branch. The idea is to use the iterative map d6(5)–(7) to first obtain droplet positions at each impact and then use the reflection conditions (12c) and (12d) to find the unique corresponding wave field. This gives the gradients at each impact, which need to be consistent with the initial guess, and also the final droplet position and velocity, which need to be consistent with the reflection conditions (12a) and (12b). The stability is analyzed through computing the eigenvalues of the linearized $N$-fold iterative map for perturbations about the periodic state, where the periodic solution is defined as asymptotically unstable if an eigenvalue lies outside the unit disc in the complex plane.

We characterize the periodic solutions in terms of the period $P$⁠, amplitude $A$⁠, and the mean energy of the wave field $E¯=P−1∫0PE(t)dt$⁠, where $E(t)$ is the wave field energy at time $t$⁠, as defined in Ref. 15. This is the additional energy of the fluid induced by the past droplet impacts, which has components of gravitational potential energy, surface energy, and the kinetic energy contribution from the potential flow within the bath. The energy $E¯$ also includes the wave field energy during droplet flight, which cannot be captured in models that neglect the oscillatory motion of the wave field between impacts.^22,42 We compare the energy to the mean energy of a bouncing droplet $E¯B$ at the given memory, where $E¯→E¯B$ as $A→0+$⁠. We also neglect the mean energy contribution from the droplet’s horizontal and vertical motions; the former is several magnitudes smaller than the mean wave energy, and the latter is constant in our model due to the imposed periodic vertical motion.¹⁵ 

We first consider the onset of small-amplitude oscillation that arises for a sufficiently weak spring constant. In the limit $A→0$⁠, the degenerate case $P=1$ describes a bouncer at the origin for a given $Γ$⁠, which is stable for $κ~>κ~c(Γ)$⁠. Thus, the bouncing state can persist beyond the free-space (⁠ $κ~=0$⁠) walking threshold; a sufficiently steep harmonic potential may trap the droplet at the origin. For $κ~<κ~c$⁠, the bouncing destabilizes via a supercritical Neimark-Sacker bifurcation, where the period of unstable oscillation $P⋆>0$ is given by the argument of the unstable complex conjugate eigenvalues of the stability matrix. A stable limit cycle forms after an initial transient, whose period $P$ and amplitude $A$ we compute directly. For sufficiently small oscillations (⁠ $A≲0.15$⁠), the period associated with the destabilising mode of the bouncer is well approximated by the limit cycle period, with $|P−P⋆|≲1$⁠, as shown in Fig. 3. In the limit $A→0+$⁠, we have $P→Pc(Γ)∈(0,∞)$⁠; this infinitesimal oscillation amplitude with a finite frequency is analogous to the small radius limit of circular orbits.²⁵ 

In Fig. 4, we show the dependence of the oscillation amplitude $A$ on the spring constant $κ~$⁠, period $P$⁠, and wave energy $E¯$⁠. For weak memory (⁠ $Γ/ΓF=0.9$⁠), all oscillations are stable (blue curves) and the amplitude grows monotonically as $κ~$ decreases. The period increases approximately linearly with the amplitude for large oscillations under a weak central force, which dominates the wave force only at the extrema of the periodic motion.

As the memory is increased (⁠ $Γ/ΓF=0.94$⁠), unstable oscillations emerge (red curves), corresponding to forbidden oscillation amplitudes. Strikingly, the unstable oscillations have the largest mean energy $E¯$⁠, and as more oscillations destabilize for $Γ/ΓF=0.96$⁠, the remaining stable oscillations (blue curves) have the lowest mean wave field energy, suggesting an underlying energy minimization principle. A similar energy minimization was also observed for circular orbits in a harmonic potential and at the bifurcation between bouncing and walking.¹⁵ The remaining stable oscillations exhibit quantization of the oscillation amplitude, with a large number of stable plateaus (blue) in the $(κ~,A)$-plane emerging for a fixed memory, as apparent in Fig. 4(c). There are, moreover, several examples of hysteresis [Fig. 4(c)]. The emergent quantization is analogous to that arising in the quantum harmonic oscillator, where the increment between energy levels $δE=ℏω$ is fixed. Similarly, the fluid system exhibits a quantization in the oscillation amplitude $A$ with fixed increment $δA≈1/2$ equal to the radial quantization increment observed for circular orbits.^15,24,26

In Fig. 5, we plot the computed pilot-wave field $η0(x,tn)$ and droplet position $Xn$ at impact over two periods of the oscillatory periodic state. When the central force dominates the wave force [Fig. 5(a)], the droplet motion is approximately sinusoidal. In contrast, at larger wave memory [Fig. 5(b)], the pilot-wave has a strong influence on the droplet’s oscillatory motion, resulting in a pronounced departure from the sinusoidal behavior.

In Fig. 6, we plot the phase space and corresponding probability distribution for simulation of the stable oscillation states with $κ~=0.012$ (corresponding to the black circles in Fig. 4). At the point of maximum range, the droplet reverses, turning over the back of its pilot-wave field, causing a sharp increase in the droplet speed, to approximately twice the free walking speed (see supplementary material). The wave field generated during previous crossings of the bath thus substantially modulates the droplet speed during transit, indicating that the weak-acceleration limit approximation is not valid in this regime.^43,44 As reported for the case of corrals,^19,21 this speed-modulation is responsible for the emergence of wavelike statistics, where the maxima of the stationary probability distribution $ρX(x)$ arise when the droplet speed is lowest. Through its modulation of the droplet speed, the wavelength of the pilot-wave thus prescribes the wavelength of the statistical wave, as is most apparent in Fig. 6(c). We see that for all values of $Γ/ΓF$⁠, the mean wave field $η¯0(x)$ and probability distribution $ρX(x)$ take a similar form on the interval $x∈[−A,A]$⁠, as expected on the basis of our convolution relationship $η¯0=ρX∗ηB$⁠.

For $Γ/ΓF=0.96$⁠, we plot the mean Faraday wave field $η¯(x)$ in Fig. 7. Since $ρX(x)$ is largest near the oscillation extrema, we see corresponding peaks in $η¯$ near the points $(x,y)=(±A,0)$⁠. Furthermore, we typically see $η¯0(x)>0$ for all $x∈[−A,A]$ since the local wave field is generally maximal near the droplet [for example, see the free-walker wave field in Fig. 1(a)]. Moreover, the symmetry about $x=0$ of the statistical distribution ensures symmetry in the mean wave field.

As $Γ/ΓF$ increases, we observe that the plateaus of stable oscillations in the $(κ~,A)$-plane become flatter and wider (see Fig. 4). We thus seek solutions where the periodic motion is sustained even in the absence of a harmonic potential (⁠ $κ~=0$⁠), in which the mean wave field traps the droplet. We note that analogous solutions exist for circular orbits at high wave memory, where the orbital radius $r0$ satisfies the quantization $J0(kFr0)=0$⁠.^{15,25,26,45,46} The periodic wave-trapped solutions of interest here are a version of these “hydrodynamic spin states” for motion confined to a line.

In Fig. 8, we plot the wave profile over time for two periods of a periodic wave-trapped solution at high memory $Me$ [as defined in Eq. (8)], which is a more useful measure of the vibrational forcing in the limit $Γ→ΓF−$⁠. Strikingly, we observe that at high memory, the wave at each impact $η0(x,tn)$ differs from the mean wave field $η¯0(x)$ only by a small perturbation. The unstable nature of this periodic state is emphasized by the fact that $η¯0(x)$ decreases rapidly for $|x|≳A$⁠, which is to say that the droplet could escape the potential trap imposed by its mean wave field for sufficiently large perturbations.

From Fig. 9(a), we observe that the amplitude $A$ of the periodic oscillation decreases as the wave memory $Me$ increases, while the oscillation period $P$ attains a minimum value before increasing at high vibrational forcing. We rationalize these dependencies in terms of the effective potential induced by the mean wave field. By applying Corollary 1, we use the convolution result to obtain the mean wave field $η¯0(x)$ over one period of the oscillatory motion, with results shown in Fig. 9(b). As $Me$ increases, $η¯0(x)$ becomes increasingly flat for $|x|≤A$⁠, resulting in a decrease in the propulsive force provided by the mean wave field. This reduces the average droplet speed, and thus the oscillation period $P$ increases. Furthermore, the steepness of the stationary cumulative probability distribution $CX(x)$ at high vibrational forcing for $|x|≈A$ [see Fig. 9(c)] indicates that the droplet spends a significant portion of the oscillation bouncing near its maximum range, which further increases the oscillation period.

To postulate a lower bound of the oscillation amplitude $A$ in the high memory limit, we exploit the fact that the droplet spends significant time near its oscillation extrema [see Fig. 9(c)], and so approximate its probability density function by $ρX(x)≈12[δ(x+A)+δ(x−A)]$⁠. An application of Corollary 1 thus yields $η¯0(x)≈12[ηB(x+A)+ηB(x−A)]$⁠. For oscillatory motion to persist, it is natural to require that the extremum at $x=0$ is a local minimum, corresponding to $η¯0″(0)>0$⁠, or equivalently, $ηB″(A)≳0$⁠. A second natural requirement is for $η¯$ to slope inwards at the point of maximum oscillation amplitude, corresponding to $η¯0′(A)>0$ (by the symmetry of $η¯0$⁠), or equivalently, $ηB′(2A)≳0$⁠. From the computation of $ηB(x)$ in the limit $Me→∞$ (as depicted in Fig. 2), the two conditions on $ηB$ are both satisfied for $0.3≲A≲0.5$⁠; we thus postulate that $A≈0.3$ is a lower bound for the amplitude of periodic wave-trapped solutions as $Me→∞$ [see Fig. 9(a)], a limit prescribed by the length scale of the Faraday waves.

Although the wave-trapped solutions are unstable in the parameter regime explored experimentally, they demonstrate that in the high memory limit, the mean Faraday wave field may trap the droplet in periodic motion. In a sense, the mean wave field $η¯$ then acts as a potential, related by Corollary 1 to the droplet’s statistical distribution through $η¯0=ρX∗ηB$⁠. Hence, the periodic motion of the droplet is in effect driven by its own stationary probability distribution $ρX$⁠. We re-explore this concept in Sec. IV C for the case of chaotic dynamics in the high memory limit.

We now consider the chaotic dynamics arising at sufficiently high memory that the periodic states destabilize via the Ruelle-Takens-Newhouse scenario (Sec. IV A). In the high memory limit, we rationalize the form of the chaotic dynamics and emergent statistics (Sec. IV B) and propose a stochastic reformulation of the pilot-wave dynamics (Sec. IV C).

As $Γ$ is increased, the periodic phase-plane orbits may destabilize into regular wobbling orbits, before transitioning to chaos. The route to chaos for circular orbits in a harmonic potential has been explored experimentally^28,47 and theoretically⁴⁸ using the stroboscopic trajectory equation.⁴² In both cases, the Ruelle-Takens-Newhouse route to chaos^32,33 was observed. According to this scenario, from a fixed point, three bifurcations induce additional incommensurate frequencies into the spectrum, after which it is likely (but not guaranteed) that a strange attractor appears in the phase space.⁴⁹ 

Following the methodology of Tambasco et al.,⁴⁸ we fix $κ~=0.03$ and initialize a simulation for a value of $Γ$ where the periodic motion is stable, as indicated by the linear stability analysis. The simulation runs for $N0+2p$ impacts, where the first $N0$ impacts are discarded to remove transient effects. We take the Fourier transform of the droplet position $Xn$ for the final $2p$ impacts (typically $p=17$⁠) and locate the frequencies $f$ corresponding to the peaks in the power spectrum $P$⁠. At the end of the simulation, we increment $Γ↦Γ+Δ(Γ/ΓF)ΓF$⁠, where $Δ(Γ/ΓF)$ is chosen adaptively to capture the bifurcations.

The fixed point of this system is a bouncer at the origin, which destabilizes via a Neimark-Sacker bifurcation (bifurcation B1), as discussed in Sec. III A. Beyond this threshold, the frequency spectrum of the resulting stable limit cycle is dominated by $f1=1/P$ and its harmonics, where $P≈63$ is as computed in Sec. III. This is highlighted by the frequency spectrum in Fig. 10(a) with accompanying phase portraits and probability density functions. At $Γ/ΓF≈0.980447$ (B2), this motion destabilizes through the emergence of complex conjugate unstable eigenvalues with oscillatory frequency $f2⋆$ [see Fig. 10(b)]. The resulting instability is saturated by nonlinear effects, leading to the quasi-periodic stable wobbling motion with incommensurate frequencies $f1$ and $f2≈f2⋆$ and their integer combinations [see Fig. 11(a)]. This evolution invokes a qualitative change in the statistics, with several peaks emerging in the droplet position stationary distribution [Fig. 10(b)]. Unlike the route to chaos of circular orbits,⁴⁸ we do not observe any frequency locking between $f1$ and $f2$⁠.

For $Γ/ΓF≳0.98050$⁠, a third bifurcation (B3) yields the incommensurate frequency $f3$⁠, as is typical of the Ruelle-Takens-Newhouse route to chaos^32,33 [see Fig. 10(c)]. While several additional peaks arise in the frequency spectrum following this bifurcation, the dynamics are still dominated by the frequencies $f1$ and $f2$ (and their harmonics), yielding a qualitatively similar probability distribution. For $Γ/ΓF≳0.980594$ (B4), additional peaks emerge in the probability distribution and the phase-portrait appears less regular [Fig. 10(d)]. In particular, the broad-banded frequency spectrum suggests chaotic dynamics, which we verify by considering the Lyapunov exponent. We follow Gilet^20,50 and consider two simulations from the same initial conditions, except for an initial perturbation in the dimensionless droplet position of $10−10$⁠, yielding trajectories $X(1)(t)$ and $X(2)(t)$⁠. As shown in Fig. 11(b), the difference $χ≡|X(1)−X(2)|$ oscillates in the interval $10−11≲χ≲10−6$ just before B4 (⁠ $Γ/ΓF=0.980593$⁠), but grows to be of order 1 just after B4 (⁠ $Γ/ΓF=0.980594$⁠), indicating a positive Lyapunov exponent and the onset of chaos.

We now consider the high memory regime (⁠ $Me≳103$⁠) in which there is a qualitative change in the dynamics. Specifically, the wave field dominates the harmonic potential so that the droplet may change the direction several times before crossing the origin, as indicated in Fig. 12. We find that the mean Faraday wave field plays a crucial role in these chaotic dynamics, giving rise to a jump-like process between a discrete set of points, the locations of which we rationalize in Sec. IV B1. In Secs. IV B2 and IV B3, we see the emergence of wavelike statistics, where the peaks correspond to the discrete turning points of the droplet motion. We then use the relationship between the droplet statistics and the mean wave field (Theorem 1) to postulate an effective potential $Ve(x)$ that influences the chaotic motion of the droplet (Sec. IV C). The additional notation used throughout this section is summarized in Table II.