A millimetric droplet may bounce and self-propel on the surface of a vertically vibrating bath, where its horizontal “walking” motion is induced by repeated impacts with its accompanying Faraday wave field. For ergodic long-time dynamics, we derive the relationship between the droplet’s stationary statistical distribution and its mean wave field in a very general setting. We then focus on the case of a droplet subjected to a harmonic potential with its motion confined to a line. By analyzing the system’s periodic states, we reveal a number of dynamical regimes, including those characterized by stationary bouncing droplets trapped by the harmonic potential, periodic quantized oscillations, chaotic motion and wavelike statistics, and periodic wave-trapped droplet motion that may persist even in the absence of a central force. We demonstrate that as the vibrational forcing is increased progressively, the periodic oscillations become chaotic via the Ruelle-Takens-Newhouse route. We rationalize the role of the local pilot-wave structure on the resulting droplet motion, which is akin to a random walk. We characterize the emergence of wavelike statistics influenced by the effective potential that is induced by the mean Faraday wave field.

A droplet may walk on the surface of a vertically vibrating fluid bath, propelled by the waves generated from all previous impacts. This hydrodynamic pilot-wave system exhibits many features that were previously thought to be exclusive to the quantum realm, such as tunneling, emergent statistics, and quantized droplet dynamics. We herein derive the relationship between the droplet’s statistical distribution and the accompanying mean pilot-wave in a very general setting. When the droplet is subject to a central force with its motion confined to a line, we rationalize a number of regimes, including periodic quantized oscillations, chaotic motion, and the emergence of wavelike statistics. In particular, we demonstrate that the mean-pilot-wave potential has a controlling influence on the droplet’s dynamics at high vibrational forcing, where the resultant droplet motion is similar to a random walk.

## I. INTRODUCTION

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 $\Gamma =A\omega 02/g$, where $A$ is the shaking amplitude, $\omega 0/(2\pi )$ is the frequency, and $g$ is the gravitational acceleration. As $\Gamma $ 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^{3} [see Fig. 1(a)]. The decay time of the Faraday waves increases with $\Gamma $ for $\Gamma <\Gamma F$, where the Faraday threshold $\Gamma 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 $\Gamma F$.^{4} 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.^{5}

Variable | Value | Description |

$\sigma $ | $2.06\xd7 10 \u2212 2 kg s \u2212 2 $ | Surface tension |

$\rho $ | $949kg m \u2212 3 $ | Fluid density |

$\nu $ | $2\xd7 10 \u2212 5 m 2 s \u2212 1 $ | Kinematic viscosity (fluid) |

$ \mu a i r $ | $1.8\xd7 10 \u2212 5 kg m \u2212 1 s \u2212 1 $ | Dynamic viscosity (air) |

$g$ | $9.8m s \u2212 2 $ | Gravity |

$ \omega 0 $ | $80\xd72\pi s \u2212 1 $ | Vibration frequency ( $\xd72\pi $) |

$ R 0 $ | $3.8\xd7 10 \u2212 4 m$ | Droplet radius |

$m$ | $2.2\xd7 10 \u2212 7 kg$ | Droplet mass |

$ \nu p $ | $1.3\xd7 10 \u2212 7 kg s \u2212 1 $ | Stokes’ drag coefficient |

Variable | Value | Description |

$\sigma $ | $2.06\xd7 10 \u2212 2 kg s \u2212 2 $ | Surface tension |

$\rho $ | $949kg m \u2212 3 $ | Fluid density |

$\nu $ | $2\xd7 10 \u2212 5 m 2 s \u2212 1 $ | Kinematic viscosity (fluid) |

$ \mu a i r $ | $1.8\xd7 10 \u2212 5 kg m \u2212 1 s \u2212 1 $ | Dynamic viscosity (air) |

$g$ | $9.8m s \u2212 2 $ | Gravity |

$ \omega 0 $ | $80\xd72\pi s \u2212 1 $ | Vibration frequency ( $\xd72\pi $) |

$ R 0 $ | $3.8\xd7 10 \u2212 4 m$ | Droplet radius |

$m$ | $2.2\xd7 10 \u2212 7 kg$ | Droplet mass |

$ \nu p $ | $1.3\xd7 10 \u2212 7 kg s \u2212 1 $ | Stokes’ drag coefficient |

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 $\lambda 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.^{18} 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,^{21} 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.^{24} In the chaotic regime arising at high vibrational forcing, a complicated switching process arises between the system’s underlying orbital states.^{28} 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.*^{30} The conceptual value of decomposing the instantaneous pilot-wave field into its mean and fluctuating components was further stressed by Labousse.^{31} 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\u2033(t)+\kappa x(t)=0$ with spring constant $\kappa $, a particle of mass $m$ enters into simple harmonic motion with fixed frequency $\omega =\kappa /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=\u210f\omega (n+1/2)$ (where $\u210f$ is the reduced Planck’s constant and $n\u2208N$), 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 $\lambda 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^{15} 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 $\Gamma \u2192\Gamma F$, periodic wave-trapped solutions arise in which the droplet’s oscillatory motion persists even in the absence of an external force ( $\kappa =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.

## II. DISCRETE-TIME MODEL

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^{34}), 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 $\varphi $ and wave perturbation $\eta $ couple with the prescribed impact pressure $PD(x,t)=f(t)\delta [x\u2212X(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\u2211n=0\u221e\delta (t/T\u2212n),$ where $T=4\pi /\omega 0$ is the Faraday period and $m$ is the droplet mass.^{38} The vibrating frame of reference introduces the effective gravity $g\Gamma (t)=g[1\u2212\Gamma cos\u2061(\omega 0t+\beta )]$, where $\beta $ denotes the droplet’s impact phase.

Following the model of Moláček and Bush,^{38} 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 $\nu p=6\pi R0\mu 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$,^{38} whose value is discussed below.

### A. Dimensionless variables

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

Typical parameter values from Table I give reciprocal Reynolds number $\u03f5\u22480.019$, Bond number $B\u22480.102$, $G\u22481.201$, $M\u22480.0017$, $R\u22480.075$, $\nu ~p\u22480.01$, and vibration number $\Omega \u22614\pi R3/B=0.8$.^{39} The dimensionless potential strength $\kappa ~\u22650$ is a free parameter of both the model and experiments, with $10\u22123\u2272\kappa ~\u227210\u22121$. The dynamics are largely insensitive to changes in the skidding friction $c$ and impact phase $\beta $; thus, we fix $c=0.33$ and $\beta =\pi $.^{15}

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

with orthogonal basis functions $\Phi m(x;k)\u2261Jm(kr)eim\theta $, where $x=(r,\theta )$ in polar coordinates and $i$ is the imaginary unit. As $\eta $ is real and $Jm(z)=(\u22121)mJ\u2212m(z)$ for all $m\u2208Z$, the complex coefficients $am$ satisfy the reality condition $a\u2212m=(\u22121)mam\u2217$ for all $m$, where $\u2217$ 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 $\eta (\u22c5,t)$ are continuous across impacts, we obtain nonlinear jumps in $X\u2032$ and $\eta t$ at impact times $t=tn\u2261n$, which appear in (6) and (7) below.

During flight ( $t\u2260tn$), the wave and droplet dynamics decouple and evolve according to

where $\gamma k\u22612\u03f5k2$, $\omega k2\u2261Gk+Bk3$, and $\omega g,k2\u2261Gk$. As (3) is of Floquet form with period $1/2$, we evolve the system between impacts ( $tn\u22121+\u21a6tn\u2212$) using principal fundamental matrices $Mk(\Gamma )\u2208R2\xd72$ for the wave amplitudes [Eq. (3)] and $F(\kappa ~)\u2208R4\xd74$ for the droplet dynamics (4). This yields the 3-stage map from $tn\u22121+\u21a6tn+$:

- Compute the next droplet position and the pre-impact velocity$X(tn)X\u2032(tn\u2212)=F(\kappa ~)X(tn\u22121)X\u2032(tn\u22121+).$
- Update the wave amplitudes $\u2200k>0$ and $\u2200m\u2208Z$$am(tn;k)am\u2032(tn+;k)=Mk(\Gamma )am(tn\u22121;k)am\u2032(tn\u22121+;k)\u2212pk\Phi m\u2217[X(tn);k]01.$
- Apply the droplet jump conditions$X\u2032(tn+)=(1\u2212F)X\u2032(tn\u2212)\u2212H\u2207\eta [X(tn),tn],$

where $F=1\u2212exp\u2061(\u2212cGR/B)$, $H=(F/c)B/R$, and $pk=kMG/(2\pi )$. Equations (6) and (7) couple through $\u2207\eta $ 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.^{4} This appears naturally from the eigenvalues of $Mk(\Gamma )$, which we write as $exp\u2061(\u2212s1)$ and $exp\u2061(\u2212s2)$ for $si=si(k,\Gamma )\u2208C$, where $0\u2264Re(s1)\u2264Re(s2)$. The dominant exponent $s1(k,\Gamma )$ is real and positive in a neighborhood of $(kF,\Gamma F)$, with $s1(kF,\Gamma F)=0$. For $\Gamma <\Gamma F$, we thus define

as $\Gamma \u2192\Gamma F$, where $Td(\Gamma )\u223c0.6$.^{15,36} While this parameter diverges as $\Gamma \u2192\Gamma 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^{40} (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.

### B. Long-time statistical behavior

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

A recent study of walker motion in corrals pointed out that the two take a similar form;^{21} 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 $\eta \xaf=\eta B\u2217\mu $. Here, $\eta 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 $\Gamma <\Gamma 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 $\eta 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 $\mu $) 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 $\eta B(x)$. Since they are not visited equally, the contribution of each point $\eta B(x)$ must be weighted by $\mu (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).

*Assuming there exists a stationary probability distribution $\mu (x)$ for the droplet position and that the system dynamics are ergodic $,$ then the mean wave field $\eta \xaf(x)$ $[$as defined by Eq.*

*(9)*$]$ satisfies*where $\eta B(x)$ is the radially symmetric wave field of a bouncer centred at the origin.*

#### Proof

^{41}to write

^{15}

We have proved similar convolution results for the models of Fort *et al.*^{22} and Oza *et al.*,^{42} 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 $\eta 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^{24} 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.*^{8} 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\u2208N$. This corollary does not require any assumptions about the existence or uniqueness of a stationary distribution, nor does it require the ergodic hypothesis.

*If there exists $Q\u2208N$ such that $X(tn+Q)=X(tn)$ for all $n,$ then the mean wave field $\eta \xaf(x)=Q\u22121\u2211n=1Q\eta (x,tn)$ satisfies*

*where $\eta B(x)$ is the radially symmetric wave field of a bouncer centred at the origin and $\mu (x)=Q\u22121\u2211n=1Q\delta [x\u2212X(tn)].$*

#### Proof

Henceforth, we consider the case where the droplet motion is confined to a line. For $\mu (x)=\rho X(x)\delta (y)$ [where $x=(x,y)$], Theorem 1 and Corollary 1 both simplify to $\eta \xaf0(x)=(\rho X\u2217\eta B)(x)$, where $\eta \xaf0\u2261\eta \xaf|y=0$ is the mean wave field along the $x$-axis. We demonstrate in Appendix B that in the case where the period $Q\u2192\u221e$, the result to Corollary 1 remains robust, even when the probability distribution $\rho X(x)$ is approximated by a histogram.

## III. PERIODIC SOLUTIONS

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

where $Gn=\u2202x\eta 0(Xn,tn)$ is the wave field gradient along the $x$-axis at impact $n$.

For any given $(\Gamma ,\kappa ~)$, 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=\phi N$ ( $N\u2208N$ and $\phi \u2208Q$) for a given $\Gamma $, and solve for $\kappa ~$ (it should be noted that in this case, there is a relationship between the oscillation period $P\u2208Q$ and the number of impacts $Q\u2208N$ such that $Xn=Xn+Q$ for all $n$). Typically, $\phi =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\u2208Z$ and $k>0$)

For given $\Gamma $ and $N$, we use a Newton method to compute the periodic states for $(N+1)$ unknowns $\theta =(X0,\kappa ~,G1,\u2026,GN\u22121)$, 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\xaf=P\u22121\u222b0PE(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\xaf$ 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\xafB$ at the given memory, where $E\xaf\u2192E\xafB$ as $A\u21920+$. 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.^{15}

### A. From bouncing to oscillating

We first consider the onset of small-amplitude oscillation that arises for a sufficiently weak spring constant. In the limit $A\u21920$, the degenerate case $P=1$ describes a bouncer at the origin for a given $\Gamma $, which is stable for $\kappa ~>\kappa ~c(\Gamma )$. Thus, the bouncing state can persist beyond the free-space ( $\kappa ~=0$) walking threshold; a sufficiently steep harmonic potential may trap the droplet at the origin. For $\kappa ~<\kappa ~c$, the bouncing destabilizes via a supercritical Neimark-Sacker bifurcation, where the period of unstable oscillation $P\u22c6>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\u22720.15$), the period associated with the destabilising mode of the bouncer is well approximated by the limit cycle period, with $|P\u2212P\u22c6|\u22721$, as shown in Fig. 3. In the limit $A\u21920+$, we have $P\u2192Pc(\Gamma )\u2208(0,\u221e)$; this infinitesimal oscillation amplitude with a finite frequency is analogous to the small radius limit of circular orbits.^{25}

### B. From classical to quantized dynamics

In Fig. 4, we show the dependence of the oscillation amplitude $A$ on the spring constant $\kappa ~$, period $P$, and wave energy $E\xaf$. For weak memory ( $\Gamma /\Gamma F=0.9$), all oscillations are stable (blue curves) and the amplitude grows monotonically as $\kappa ~$ 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 ( $\Gamma /\Gamma F=0.94$), unstable oscillations emerge (red curves), corresponding to forbidden oscillation amplitudes. Strikingly, the unstable oscillations have the largest mean energy $E\xaf$, and as more oscillations destabilize for $\Gamma /\Gamma 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.^{15} The remaining stable oscillations exhibit quantization of the oscillation amplitude, with a large number of stable plateaus (blue) in the $(\kappa ~,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 $\delta E=\u210f\omega $ is fixed. Similarly, the fluid system exhibits a quantization in the oscillation amplitude $A$ with fixed increment $\delta A\u22481/2$ equal to the radial quantization increment observed for circular orbits.^{15,24,26}

In Fig. 5, we plot the computed pilot-wave field $\eta 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 $\kappa ~=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 $\rho 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 $\Gamma /\Gamma F$, the mean wave field $\eta \xaf0(x)$ and probability distribution $\rho X(x)$ take a similar form on the interval $x\u2208[\u2212A,A]$, as expected on the basis of our convolution relationship $\eta \xaf0=\rho X\u2217\eta B$.

For $\Gamma /\Gamma F=0.96$, we plot the mean Faraday wave field $\eta \xaf(x)$ in Fig. 7. Since $\rho X(x)$ is largest near the oscillation extrema, we see corresponding peaks in $\eta \xaf$ near the points $(x,y)=(\xb1A,0)$. Furthermore, we typically see $\eta \xaf0(x)>0$ for all $x\u2208[\u2212A,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.

### C. Wave-trapped solutions

As $\Gamma /\Gamma F$ increases, we observe that the plateaus of stable oscillations in the $(\kappa ~,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 ( $\kappa ~=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 $\Gamma \u2192\Gamma F\u2212$. Strikingly, we observe that at high memory, the wave at each impact $\eta 0(x,tn)$ differs from the mean wave field $\eta \xaf0(x)$ only by a small perturbation. The unstable nature of this periodic state is emphasized by the fact that $\eta \xaf0(x)$ decreases rapidly for $|x|\u2273A$, 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 $\eta \xaf0(x)$ over one period of the oscillatory motion, with results shown in Fig. 9(b). As $Me$ increases, $\eta \xaf0(x)$ becomes increasingly flat for $|x|\u2264A$, 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|\u2248A$ [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 $\rho X(x)\u224812[\delta (x+A)+\delta (x\u2212A)]$. An application of Corollary 1 thus yields $\eta \xaf0(x)\u224812[\eta B(x+A)+\eta B(x\u2212A)]$. For oscillatory motion to persist, it is natural to require that the extremum at $x=0$ is a local minimum, corresponding to $\eta \xaf0\u2033(0)>0$, or equivalently, $\eta B\u2033(A)\u22730$. A second natural requirement is for $\eta \xaf$ to slope inwards at the point of maximum oscillation amplitude, corresponding to $\eta \xaf0\u2032(A)>0$ (by the symmetry of $\eta \xaf0$), or equivalently, $\eta B\u2032(2A)\u22730$. From the computation of $\eta B(x)$ in the limit $Me\u2192\u221e$ (as depicted in Fig. 2), the two conditions on $\eta B$ are both satisfied for $0.3\u2272A\u22720.5$; we thus postulate that $A\u22480.3$ is a lower bound for the amplitude of periodic wave-trapped solutions as $Me\u2192\u221e$ [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 $\eta \xaf$ then acts as a potential, related by Corollary 1 to the droplet’s statistical distribution through $\eta \xaf0=\rho X\u2217\eta B$. Hence, the periodic motion of the droplet is in effect driven by its own stationary probability distribution $\rho X$. We re-explore this concept in Sec. IV C for the case of chaotic dynamics in the high memory limit.

## IV. CHAOTIC DYNAMICS

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).

### A. Transition to chaos

As $\Gamma $ 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^{48} using the stroboscopic trajectory equation.^{42} 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.^{49}

Following the methodology of Tambasco *et al.*,^{48} we fix $\kappa ~=0.03$ and initialize a simulation for a value of $\Gamma $ 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 $\Gamma \u21a6\Gamma +\Delta (\Gamma /\Gamma F)\Gamma F$, where $\Delta (\Gamma /\Gamma 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\u224863$ 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 $\Gamma /\Gamma F\u22480.980447$ (B2), this motion destabilizes through the emergence of complex conjugate unstable eigenvalues with oscillatory frequency $f2\u22c6$ [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\u2248f2\u22c6$ 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,^{48} we do not observe any frequency locking between $f1$ and $f2$.

For $\Gamma /\Gamma F\u22730.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 $\Gamma /\Gamma F\u22730.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\u221210$, yielding trajectories $X(1)(t)$ and $X(2)(t)$. As shown in Fig. 11(b), the difference $\chi \u2261|X(1)\u2212X(2)|$ oscillates in the interval $10\u221211\u2272\chi \u227210\u22126$ just before B4 ( $\Gamma /\Gamma F=0.980593$), but grows to be of order 1 just after B4 ( $\Gamma /\Gamma F=0.980594$), indicating a positive Lyapunov exponent and the onset of chaos.

### B. The high memory limit

We now consider the high memory regime ( $Me\u2273103$) 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.

Variable | Description |

$ \tau i \u2208 T $ | Turning times (raw data) |

$ T i \u2261X( \tau i )$ | Turning positions (raw data) |

$C(t)$ | Trend curve on the slow timescale |

$ X n R $ | Residual droplet impact positions |

$ T i R $ | Residual turning positions |

$ \rho X (x)$ | Probability distribution of $Xn$ (raw data) |

$ \rho T (x)$ | Probability distribution of $Ti$ (raw data) |

$ \rho X R (x)$ | Probability distribution of $XnR$ (residuals) |

$ \rho T R (x)$ | Probability distribution of $TiR$ (residuals) |

Variable | Description |

$ \tau i \u2208 T $ | Turning times (raw data) |

$ T i \u2261X( \tau i )$ | Turning positions (raw data) |

$C(t)$ | Trend curve on the slow timescale |

$ X n R $ | Residual droplet impact positions |

$ T i R $ | Residual turning positions |

$ \rho X (x)$ | Probability distribution of $Xn$ (raw data) |

$ \rho T (x)$ | Probability distribution of $Ti$ (raw data) |

$ \rho X R (x)$ | Probability distribution of $XnR$ (residuals) |

$ \rho T R (x)$ | Probability distribution of $TiR$ (residuals) |

To gain further understanding of the pilot-wave dynamics in this regime, it is useful to recast the iterative map d6(5)–(7) as a trajectory equation for the droplet position $Xn$ and the mean droplet velocity during flight $Un\u2261Xn+1\u2212Xn$. By computing the droplet fundamental matrix $F(\kappa ~)$ analytically, the droplet’s evolution may be expressed as

where $D=1\u2212(1\u2212F)e\u2212\nu ~p$ is a drag coefficient, and

is the time-dependent full pilot-wave potential, which is the sum of the applied harmonic potential and the wave field at each impact. In the vicinity of the origin ( $|x|\u22723$ in Fig. 12), the full pilot-wave potential at each impact $Vp(x,tn)$ oscillates in *x*. However, as $|x|\u2192\u221e$ the instantaneous wave field decays and we observe that $Vp(x,tn)\u224812Kx2$ for all time. In (15), $K(\kappa ~)>0$ determines the strength of the time-averaged harmonic potential over one impact period and $F(\kappa ~)>0$ prescribes the magnitude of the wave force (whose dependence on $\kappa ~$ is weak).^{51} The system (13) and (14) and the full pilot-wave potential $Vp$ will be referred to throughout Secs. IV B and IV C.

#### 1. The random walk dynamics

In Fig. 12, we plot the evolution of the full pilot-wave potential $Vp(x,tn)$ and the corresponding droplet position $Xn$ at successive impacts in the high memory regime. To understand the role of the long-lived Faraday waves in this regime, we plot the spatial minima of $Vp(x,tn)$ at each impact, from which two important observations emerge. First, the minima far from the droplet (typically $\u22731$ Faraday wavelength away) remain at a roughly constant position over time, indicating the potential has an underlying stationary structure induced by the wave field. Second, when the droplet changes direction (at which point it is moving slowly), the local-pilot wave accumulates, increasing the droplet’s potential energy, from which the droplet departs and heads towards one of the neighboring potential minima. Depending on the prior dynamics, the droplet will turn around again at one of the minima of $Vp(x,tn)$ on its path.

To analyze these dynamics, we define the set of turning times $T\u2282N$ to be the times at which the droplet changes direction. That is to say, if $\tau i\u2208T$, then $X(\tau i)$ is a local extremum and $Ti\u2261X(\tau i)$ is defined to be a turning point. In the droplet trajectory time-series data in Fig. 13(a), the turning times $\tau i$ and positions $Ti$ correspond to the red dots. Furthermore, it appears that the droplet changes position only in the vicinity of specific points on the bath and that there is an apparent structure to the distance between turning points $Di=|Ti+1\u2212Ti|$. Indeed, by plotting the distribution $\rho D$ of distances $Di$ [see Fig. 13(b)], it emerges that the distance between turning points is quantized, where $\rho D$ has sharp maxima at points approximated by the set

The emergence of this quantization lies in the combined structure of the global standing wave field and the wave field generated by the droplet at each impact, whose shape is approximated by $J0(kFx)$. From the observations in Fig. 12, it becomes clear that it is the minima of $J0(kFx)$ that plays a role in prescribing the quantized distance between turning points, with values in the set $D$. This correspondence is shown in Fig. 13(b). In what follows, we rationalize these dynamics by considering a jump process, before postulating a stochastic model in Sec. IV C.

We proceed by presenting a simple geometric argument that demonstrates the role of the quantized distance between turning points in the long-time statistics. We consider a Markovian jump process $(xn)n\u22650$ between turning points, where the jump distances $dn\u2261|xn\u2212xn\u22121|$ are *exactly* restricted to $dn\u2208D$. In accordance with our observations of the pilot-wave system in the high memory limit (see Figs. 12 and 13), we require that the set of possible points visited by the jump process forms a communicating class with symmetry preserved about the origin.

We denote $\alpha \u2208R$ as a position visited by the jump process (whose possible values are determined in the following analysis), and without loss of generality, we set $x0=\alpha $ and consider $x1>x0$. Using the assumed structure of $D=0.6+N$, we define $Nn\u2208N$ such that $dn=0.6+Nn$. As the droplet changes direction at each turning point, we observe that after an even number of jumps

so $|x2n\u2212x0|\u2208N$ for all $n\u2208N$. By a similar calculation, we find that after an odd number of jumps

Thus, for all points in the jump process to form a communicating class, we require $xn\u2208M(\alpha )$ for all $n\u22650$, where $\alpha $ parametrizes the mesh

We note that this mesh is periodic with period 1, so without loss of generality, we restrict the displacement of the mesh to $\alpha \u2208[\u221212,12)$. For symmetric statistics about the origin, we require that $\alpha $ be such that $M(\alpha )$ is also symmetric, which yields $\alpha \u2208{\u22120.3,0.2}$. For consistency with the jump distances $dn\u2208D$, the droplet may only leave each mesh point in a fixed direction (as depicted in Fig. 14), namely, to the right for $xn\u2208\alpha +Z$ and to the left for $xn\u2208\alpha +0.6+Z$.

#### 2. Detrending the long-time statistics

From the analysis in Sec. IV B 1, we expect the turning points $Ti$ (and the peaks of the corresponding probability distribution $\rho T$) to be determined by the meshes $M0=M(0.2)$ and $M1=M(\u22120.3)$, where

and $M0$ is a translation of $M1$ by $1/2$. The jump lengths $dn\u2208D$ impose a unique direction to leave each mesh point, yielding qualitatively different dynamics, as highlighted by the schematic diagram in Fig. 14. Due to the finite width of the distribution $\rho D$ about each of its modes [with typical value $\u22720.25$ — see Fig. 13(b)], there is a corresponding finite width in the turning point distribution about each predicted mesh point. Hence, these distributions may overlap for mesh points spaced 0.4 apart but are well separated for mesh points spaced 0.6 apart. In the turning points’ time series, this yields a thicker “band structure” between mesh points spaced 0.4 apart, as seen in Fig. 15(a), where the central mesh points are visited most frequently. By symmetry, we expect the central band $[\u22120.2,0.2]$ to dominate the statistics in the case of $M0$ but the two bands $[\u22120.7,\u22120.3]$ and $[0.3,0.7]$ to be equally dominant for $M1$.

Our study reveals an additional complication; specifically, the finite width about the peaks in $\rho D$ allows for a slow translation in the dominant turning point locations, as is evident in Fig. 15(a). The translation occurs on a slow timescale, comparable to the memory time $Me$, the timescale at which the global wave field structure changes. This drift obscures the structure of the underlying statistics induced from the short-time dynamics; for example, there is only a weak structure apparent in the distribution of turning points $\rho T$ in Fig. 15(b).

To remedy this, we *detrend* the time-series data using statistical methods and then analyse the residuals. This detrending involves finding a smooth best fit $C(t)$ for the time varying drift and re-expressing the variation in the data about $C(t)$. This technique visibly enhances the wavelike nature of the droplet’s statistical distribution [see Figs. 15(b) and 15(c)] and allows us to compare the resultant dynamics to the predicted random walk meshes $M0$ and $M1$. The details of the statistical techniques used are given in Appendix D.

#### 3. Results

We now explore the statistical distributions following the detrending of the slow timescale dynamics. By defining $R(t)\u2261X(t)\u2212C(t)$, we have impact residuals $XnR\u2261R(tn)$ for all $n\u22650$ and turning point residuals $TiR\u2261R(\tau i)$ for all $\tau i\u2208T$, with respective residual probability distributions $\rho XR$ and $\rho TR$. We demonstrate that the distribution modes vary with the relative strength of the central and wave forces, and are intrinsically linked to the mean wave field $\eta \xaf0(x)$ and an associated effective potential $Ve(x)$ to be defined in Sec. IV C.

Examples of the corresponding residual distributions are given in Figs. 15(b) and 15(c), where the residual statistics are symmetric relative to mesh $M0$. The modes of $\rho XR$ correspond to the modes of $\rho TR$ since the droplet is moving slowest at the turning points, so spends most of its time in their vicinity. The harmonic potential dominates the wave field far from the origin, which explains the slight discrepancy between the distribution modes and the mesh points for large $|x|$. We note that the sub-mesh points ${\xb11.2,\xb12.2,\u2026}$ are visited less frequently as these drive the droplet away from the origin (see Fig. 14), countering the harmonic potential.

To explore the extent of the random walk-like dynamics, we vary the parameters $\kappa ~$ and $Me$ and present the results in Fig. 16. When $Me$ is fixed, the quantization is sharper when the waves dominate the harmonic potential [Fig. 16(a)], but as $\kappa ~$ is increased, the peaks become broader and the quantization loses clarity [Fig. 16(c)]. The plot of successive turning points [Fig. 17(a)] confirms that the droplet motion is consistent with the directional arrows predicted by the mesh $M0$ (see Fig. 14). However, it is relatively rare for the droplet to cross the centre of the bath (corresponding to $TiRTi+1R<0$), a feature that we rationalize in Sec. IV C.

When the wave memory $Me$ is reduced (with $\kappa ~$ fixed), the random walk-like dynamics shift to the $M1$ mesh, where the sub-mesh ${\xb10.3,\xb11.3,\u2026}$ dominates the distribution modes [see Figs. 16(b) and 16(d)]. Although the mesh points ${\xb10.7,\xb11.7,\u2026}$ remain apparent in the residual turning point distribution $\rho TR$ [Fig. 17(b)], their presence is obscured in $\rho XR$ [Fig. 16(b)]. As $Me$ is further decreased, the mesh points that counter the harmonic potential ( ${\xb10.7,\xb11.7,\u2026}$) are visited less frequently [Fig. 16(d)]. Indeed, it appears from Fig. 17(b) that the centre of the bath is crossed more frequently in this regime, as the relative strength of the central force is more pronounced at lower wave memory.

These random walk-like dynamics differ substantially from those arising in a bath driven at two incommensurate frequencies^{52} and those in a corral given by the toy model of Gilet.^{20,50} In our case, the domain is unbounded, so the allowable steps between turning points are dominated by the structure of the droplet’s local wave field. The associated random walk mesh ( $M0$ or $M1$) is selected by the relative strength of the central and wave forces, where the mesh $M0$ is dominant in the high memory limit. In contrast, the random walk-like motion observed by Gilet is instead induced by the global wave field given by the corral’s cavity modes, with a fixed random walk step size of $\lambda F/2$.

### C. The mean-pilot-wave potential

Based on the ideas of Theorem 1, we start by considering an effective potential $Ve(x)$ using the stationary residual distribution $\rho XR(x)$ and the applied harmonic potential $12Kx2$, defining $Ve(x)=12Kx2+F(\rho XR\u2217\eta B)(x)$. Remarkably, the direction associated with each mesh point (as given by Fig. 14) corresponds precisely to the gradients of $Ve$, as indicated by the arrows in Fig. 18. This correspondence provides a strong indicator that the chaotic motion of the droplet is driven by an effective potential induced by the slow decay of the pilot-wave field in the high memory limit. With this observation in mind, we sketch a stochastic reformulation of the long-time pilot-wave dynamics in the high memory limit, from which we aim to derive an equation for the time-dependent probability distribution $\rho (x,t)$ for the droplet’s position.

Following a similar idea to that proposed by Labousse *et al.*,^{30} we decompose the pilot-wave dynamics using its contrasting short and long timescale behavior. Specifically, we model the contribution of the wave field to the pilot-wave dynamics in terms of a propulsive nonlinear drag $\u2212D(Un)Un$ (similar to that used in the weak acceleration limit),^{43,44} an approximation for the effect of the long-lived Faraday waves, and a mean-zero normally distributed random noise that accounts for the local fluctuations of the pilot-wave. Using the fact that $|Un|\u226a1$ (i.e., the distance between successive impacts is small relative to the Faraday wavelength), we approximate (13) and (14) by the continuous limit, in which the Gaussian noise is replaced by an increment of the Wiener process $Wt$ over an infinitesimal timestep $dt$. This yields Langevin evolution equations for the position-velocity process $(Xt,Ut)$ at time $t>0$

where $\sigma 0>0$ prescribes the magnitude of the stochastic forcing. Here, we have defined the stochastic potential

where $\rho (x,t)$ is the time-dependent probability distribution for the droplet’s position.

The system (17) and (18) is speculative, and it should be noted that, unlike Theorem 1, the convolution $(\eta B\u2217\rho )(x,t)$ is for the time-dependent probability distribution and not for the stationary probability distribution. However, an initial condition $\rho (x,0)=\delta (x)$ would correspond to prescribing the initial pilot-wave field as that of a bouncer, which is consistent with the numerical simulations of Sec. IV B. Moreover, if a stationary probability distribution $\rho s(x)$ were to exist [where $\rho (x,t)\u2192\rho s(x)$ as $t\u2192\u221e$], then the system (17) and (18) would be consistent with the results of Theorem 1.

The evolution of the time-dependent joint probability distribution $p(x,u,t)$ corresponding to (17) and (18) is governed by the Vlasov-Fokker-Planck equation

where $\rho (x,t)=\u222bRp(x,u,t)du$ is the marginal distribution and $V(x,t)$ is defined in Eq. (19). An interesting aspect of this equation is the nonlinearity and spatial nonlocality in $p(x,u,t)$ arising through $V(x,t)$. Indeed, similar equations have been used in granular flow,^{53} and it has been proved that such equations yield a unique stationary probability under suitable assumptions for the nonlinear drag $D$, the applied potential, and the convolution kernel $\eta B$.^{54,55} Self-propulsive particles in the case of no spatial nonlocality ( $F=0$) have also been studied in a biological context.^{56} The numerical solution to (20), with the possible inclusion of a velocity-dependent multiplicative noise $\sigma 0(u)$, will be the subject of future work.

While the case without self-propulsion $D(u)=D0$ is not appropriate for modeling the dynamics of walking droplets, we note that the stationary distribution $\rho s(x)$ to Eq. (20) satisfies Kramer’s equation for a given potential $V$, with implicit solution

where $E=\sigma 02/(2D0)$ and $Z0$ is a normalisation constant.^{57} In Fig. 19, the numerical solution to (21) with different parameter values (solved using a Newton method) yields wavelike stationary statistics, a feature consistent with not only the pilot-wave dynamics of this system (Fig. 16) but also pilot-wave dynamics under a Coriolis force^{23,27} and motion confined to a corral.^{19,21} This provides a strong indication that the stochastic system (17) and (18) with the corresponding Vlasov-Fokker-Planck equation (20) will still exhibit wavelike statistics when the nonlinear drag $D(u)$ is included.

## V. DISCUSSION

We have studied the dynamics of a droplet walking in a harmonic potential with its motion confined to a line. By performing linear stability analysis of the periodic states, we have captured the changes to the limit cycle dynamics as the wave force begins to dominate the harmonic potential. In particular, we have elucidated the oscillation amplitude quantization that appears at higher wave memory, which is analogous to the energy quantization in the quantum harmonic oscillator. We have also demonstrated that the pilot-wave has the lowest mean energy for stable oscillations, suggesting the significance of an underlying energy minimization principle in rationalising the quantized states.

The methods developed herein for analyzing periodic orbits are readily adaptable for studying the droplet motion in a harmonic potential without restricting the motion to a line, which will be useful for further characterization of the more exotic periodic orbital states observed in the laboratory (e.g., lemniscates and trefoils).^{24,28} We expect some of these orbital states to be related by a (currently unknown) unstable branch in the parameter space, which is likely to connect two local minima of the wave’s mean energy. Additionally, this methodology will allow for further analysis of the periodic motion observed between two droplets (in free-space), such as promenading pairs^{15–17} and wobbling orbits.^{13–15}

We have demonstrated that this system follows the Ruelle-Takens-Newhouse route to chaos, provided that the periodic state destabilizes via a pair of complex-conjugate eigenvalues. Furthermore, each of the new incommensurate frequencies that emerges after each of the first two bifurcations is approximated by the frequency of the corresponding unstable state, as predicted by the linear stability analysis. This result is a useful verification of our stability analysis and allows us to predict the dynamics of the quasi-periodic orbits.

Finally, we have uncovered the relationship between the mean wave field and the droplet statistics (Theorem 1), which represents a powerful diagnostic tool at extremely high wave memory. In this high memory regime, the droplet motion is reminiscent of a random walk, where the distance between successive turning points is prescribed by the minima of the local pilot-wave. By detrending the slow-timescale variations in the droplet’s trajectory, we have highlighted the wavelike nature of the statistics, as becomes more pronounced at higher memory. We expect our approach to reveal the underlying statistical structure in other experimental configurations of this pilot-wave system, such as tunneling^{8,10} and in corrals.^{19,21}

Remarkably, the mean wave field yields an effective potential that has a controlling influence on the droplet dynamics and thus the emergent statistics. This draws further parallels to Bohmian mechanics, in which the statistical and guiding wave fields are identical.^{58} Furthermore, we have proposed a Langevin equation to describe the dynamics in the high memory limit, where the motion is subject to an effective potential. By expressing the stationary probability distribution $\rho s(x)$ as a (nonlinear) Vlasov-Fokker-Planck equation, we can solve directly for $\rho s(x)$. We hope that these developments will lead to a fruitful comparison of the long-time behavior of this pilot-wave system in the chaotic regime to both statistical mechanics and Bohmian mechanics.

We expect the connection between the dynamics and statistics elucidated here to apply in other experimental configurations (such as corrals^{19,21}) or indeed in a more generalized pilot-wave framework.^{59} The generalization of Theorem 1 (as given by Appendix A) will play a key role in elucidating the link between the dynamics and statistics of pilot-wave systems and may provide a tool for better understanding the ingredients required for observing quantumlike behavior on a classical scale.

## ACKNOWLEDGMENTS

M.D. gratefully acknowledges support through a scholarship from the EPSRC Centre for Doctoral Training in Statistical Applied Mathematics at Bath (SAMBa) under the project EP/L015684/1. J.W.M.B. gratefully acknowledges the support of the NSF through Grant Nos. DMS-1614043 and CMMI-1727565. P.A.M. gratefully acknowledges support through the EPSRC project EP/N018176/1.

## SUPPLEMENTARY MATERIAL

See supplementary material for videos of the periodic oscillatory pilot-wave dynamics.

### APPENDIX A: STATISTICS IN A GENERALIZED PILOT-WAVE FRAMEWORK

In a generalized (linear) pilot-wave setting (with instantaneous impacts), the evolution of the fluid variables $u(x,t)$ may be written as $u(x,tn+)=Lu(x,tn\u22121+)+uI[x,X(tn)]$, where $uI(x,y)$ is the (not necessarily symmetric) pilot-wave field generated by a single impact at time $t=0$ centred at $x=y$. The linear operator $L$ maps the fluid variables from $tn\u22121+\u21a6tn\u2212$ with operator norm $||L||op<1$ corresponding to a dissipative system. Assuming a stationary distribution $\mu (y)$ exists and that the dynamics are ergodic, then ergodicity gives $u\xaf(x)=Lu\xaf(x)+\u222bR2uI(x,y)\mu (y)dy$, where $u\xaf(x)\u2261limN\u2192\u221e1N\u2211n=1Nu(x,tn+)$. Thus,

where $uB(x,y)\u2261(Id\u2212L)\u22121uI(x,y)$ and $Id$ is the identity operator. From the Neumann series $\u2211n=0\u221eLn=(Id\u2212L)\u22121$, we recognize $uB(x,y)$ as being proportional to the time-periodic Green’s function for the domain centred at $x=y$, which is analogous to the wave field of a bouncer in a generalized framework.

### APPENDIX B: ROBUSTNESS OF THE CONVOLUTION RESULT

To demonstrate the robustness of Corollary 1, we simulate the droplet motion in a parameter regime that corresponds to the stable periodic motion (see Sec. III) and compute the corresponding histogram $H(x)$ to approximate the droplet’s probability distribution $\rho X(x)$ [Fig. 20(a)]. Thus, for histogram bin centres $\xi j$ with heights $H(\xi j)$, we have $H(\xi j)\u2248\rho X(\xi j)$. For $NX$ equally spaced points $xi$ in the interval $[\u22123,3]$, we compute the mean wave field $\eta \xaf0C(xi)$ using the convolution (10), with the midpoint quadrature rule $\eta \xaf0C(xi)\u2248\delta H\u2211j\eta B(xi\u2212\xi j)H(\xi j)$. To compare the simulated mean wave field $\eta \xaf0S$, we compute the mean squared error $MSE\u2261NX\u22121\u2211i=1NX[\eta \xaf0S(xi)\u2212\eta \xaf0C(xi)]2$, which has size $O(\delta H3)$ as $\delta H\u21920$ [Fig. 20(c)], thus indicating convergence.

### APPENDIX C: ANALYSIS OF THE PERIODIC STATES

Following from Sec. III, we perform a Newton iteration to find the periodic states. Specifically, we solve $G(\theta )=0$, where $\theta =(X0,\kappa ~,G1,\u2026,GN\u22121)$ and $G$ (of dimension $N+1$) is given below. The function $G$ is dependent on several other functions of $\theta $, which are computed at each step of the following algorithm. Hence, computation of the Jacobian $\u2202G/\u2202\theta T$ requires an application of the chain rule, where the derivative $\u2202/\u2202\theta T$ of each newly defined function is also computed. For an initial guess $\theta $:

Use $\kappa ~$ and (5) to uniquely find $V0+(\theta )$ such that $X1\u2212X0=0$ ( $X0$ is an extremum).

Use the droplet iteration maps (5) and (7) with gradients $Gn$ and the initial conditions $[X0,V0+(\theta )]T$ to compute positions $Xn(\theta )$ and velocities $Vn\u2212(\theta )$ for $n=1,\u2026,N$.

- For the wave field $\eta 0$ to satisfy the reflection conditions (12c) and (12d) with impacts $Xn(\theta )$, use (6) to find the initial wave amplitudes $am(t0;k,\theta )$ and $am\u2032(t0+;k,\theta )$, which solvewhere$[(\u22121)mI\u2212MkN]am(t0;k,\theta )am\u2032(t0+;k,\theta )=\u2212Hm(k;\theta )01,$$Hm(k;\theta )=pk\u2211n=1NJm[kXn(\theta )]MkN\u2212n.$
Use (6) to recover the wave field $\eta 0(x,tn;\theta )$ and gradients $gn(\theta )=\u2202x\eta 0(Xn,tn;\theta )$.

Using $gN(\theta )$ and $VN\u2212(\theta )$ with (7), compute $VN+(\theta )$.

If $||G(\theta )||\u221e<TOL$, stop. Otherwise, update $\theta $ with a Newton iteration and return to step 1.

To analyse the stability of the $N$-step periodic states, we extend the method used for the 1-step stability maps explored by Durey and Milewski,^{15} where perturbations are now restricted to the $x$-axis. In brief, we linearize the map d6(5)–(7) about the periodic state at times $tn+$ for $n=1,\u2026,N$. By expressing all the perturbed variables at time $tn\u22121+$ as a single column vector, we construct (sparse) transition matrices $Tn(\Gamma )$ to map the perturbed variables from $tn\u22121+\u21a6tn+$ for $n=1,\u2026,N$. The $N$-step stability matrix $T$ is the product $T=RTN\cdots T1,$ where $R$ is the diagonal reflection matrix about the $x$-axis. The eigenvalues of $T$ are computed numerically, and the periodic state is defined to be asymptotically unstable if at least one eigenvalue lies outside the unit disc in the complex plane.

### APPENDIX D: DETRENDING THE LONG-TIME STATISTICS

To detrend the data, we fit a simple version of a generalized additive model^{60} to one of the aforementioned central bands of turning points. This yields a subset of turning point times $S\u2282T$, which corresponds to the black data points in Fig. 15(a). This detrending technique is a form of regression, in which the trend curve $C(t)$ is expressed as a linear combination of smooth linearly independent basis functions (in this case, B-splines) whose weights are computed to give a least-squares fit of the data. However, to avoid over-fitting of the data [characterized by an excessively “wiggly” function $C(t)$], we introduce a smoothing penalization term.

As the trend changes over a timescale comparable to the memory time $Me\u226b1$, we consider a linear combination of $K$ basis functions $bj(t)$, where $KMe$ is the simulation duration. The trend function $C(t)$ is thus given by the linear combination $C(t)=\u2211j=1K\beta jbj(t),$ where $\beta j$ are the unknown coefficients. The penalty for over-fitting is chosen to minimize variation in the basis function coefficients $\beta j$, where the required smoothness is determined by the parameter $\theta >0$. The coefficients $\beta =(\beta 1,\u2026,\beta K)$ are then defined as the minimizer

with the smoothness penalization term

Although methods exist to find the “optimal” value of $\theta $ for a given dataset,^{60} it is sufficient for our purposes to simply fix $\theta =500$ for all datasets considered, where the residual statistics vary only weakly for $100\u2272\theta \u22721000$.

## References

*Encyclopedia of Complexity and Systems Science*, edited by R. Meyers (Springer, New York, 2009), pp. 3097–3113.