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

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 ﬁeld. For ergodic long-time dynamics, we derive the relationship between the droplet’s stationary statistical distribution and its mean wave ﬁeld in a very general setting. We then focus on the case of a droplet subjected to a harmonic potential with its motion conﬁned 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 inﬂuenced by the effective potential that is induced by the mean Faraday wave ﬁeld. © 2018 All otherwise noted, 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 conﬁned 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 inﬂuence on the droplet’s dynamics at high vibrational forcing, where the resultant droplet motion is similar to a random walk.

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 pilotwave 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 = Aω 2 0 /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 3 [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 . 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 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][7][8][9][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][12][13][14][15] promenading pairs, [15][16][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 pilotwave, 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][23][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 pathmemory 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 E n = ω(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 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 pilotwave 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.

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 quasipotential flow, which includes weak dissipative effects at high Reynolds number. [35][36][37] The harmonic velocity potential φ and wave perturbation η couple with the prescribed impact  15 for / F = 0.97. The droplet is located at the origin and walks to the right along the line y = 0, where x = (x, y). (b) Two-dimensional schematic diagram of the fluid system with free surface η 0 ≡ η| y=0 . The forces acting on the droplet are denoted by red arrows, including the central force κX(t) that acts towards the origin O. The system parameters considered in our simulations and analysis are given in Table I. 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. 38 The vibrating frame of reference introduces the effective gravity g (t) = g 1 − cos(ω 0 t + β) , where β 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 ν p = 6π R 0 μ 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 λ F = 0.51 cm, time t with the subharmonic bouncing period T = 4π/ω 0 , force f with f 0 = mg, and pressure P D with P 0 = f 0 /λ 2 F . This yields the following dimensionless parameters: Typical parameter values from Table I 39 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 β = π . 15 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 with orthogonal basis functions m (x; k) ≡ J m (kr)e imθ , where x = (r, θ) in polar coordinates and i is the imaginary unit. As η is real and J m (z) = (−1) m J −m (z) for all m ∈ Z, the complex coefficients a m satisfy the reality condition a −m = (−1) m a * m for all m, where * denotes the complex conjugate. This basis decomposition yields a system of inhomogeneous damped Mathieu equations for the wave amplitudes a m , 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 = t n ≡ n, which appear in (6) and (7) below.
The wave "memory" M e 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 M k ( ), which we write as exp(−s 1 ) and exp(−s 2 ) for s i = s i (k, ) ∈ C, where 0 ≤ Re(s 1 ) ≤ Re(s 2 ). The dominant exponent s 1 (k, ) is real and positive in a neighborhood of (k F , F ), with s 1 (k F , F ) = 0. For < F , we thus define as → F , where T d ( ) ∼ 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. (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 (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 μ(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; 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η = η 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). Theorem 1. Assuming there exists a stationary probability distribution μ(x) for the droplet position and that the system dynamics are ergodic, then the mean wave fieldη(x) [as defined by Eq. (9)] satisfies where η B (x) is the radially symmetric wave field of a bouncer centred at the origin. Proof. We define a n m (k) ≡ a m (t n ; k), a m (t + n ; k) T and rewrite (6) as where e j is the jth basis vector. We then defineā m (k) ≡ lim N→∞ 1 N N n=1 a n m (k) to be the wave amplitudes corresponding to the mean wave fieldη. By taking the mean of (11) over N impacts and considering the limit N → ∞, the ergodic theorem allows for the replacement of time averages in the last term with spatial averages, givinḡ where we have used Graf's addition theorem 41 to write The result (10) follows since the wave field of a bouncer centred at the origin for given is η B (x) = ∞ 0 ka B (k; )J 0 (k|x|) dk. 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 η 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 pilotwave 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(t n+Q ) = X(t n ) 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.
where η B (x) is the radially symmetric wave field of a bouncer centred at the origin and Using the definition of a n m (k) from the proof of Theorem 1, we take the sum Q n=1 of both sides of Eq. (11), giving By the assumed periodicity, we observe that a 0 The conclusion of the proof is identical to that of Theorem 1.
Henceforth, we consider the case where the droplet motion is confined to a line. For , Theorem 1 and Corollary 1 both simplify tō η 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.

III. PERIODIC SOLUTIONS
We seek periodic solutions to the nonlinear discrete-time map (5)-(7) with motion restricted to the x-axis, so X(t) ≡ X (t), 0 and a m ∈ R for all m. For notational convenience, we denote 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 X n = X n+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 θ = X 0 ,κ, G 1 , . . . , G N−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 (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 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Ē 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Ē B at the given memory, whereĒ →Ē 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. 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 → 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 → P c ( ) ∈ (0, ∞); 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κ, period P, and wave energyĒ. 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Ē, 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. 15 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, t n ) and droplet position X n 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.

C. Wave-trapped solutions
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 r 0 satisfies the quantization J 0 (k F r 0 ) = 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 M e [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, t n ) 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 M e 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 M e 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 C X (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 by ]. An application of Corollary 1 thus yieldsη 0 ( 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 M e → ∞ (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 M e → ∞ [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.

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).  0). The impact wave field η 0 (x, t n ) (blue curves) and droplet position X n (dots) are shown over two oscillation periods for P = 42. The black curve is the corresponding mean wave fieldη 0 (x). The wave memory is (a) M e ∼ 47.7 (corresponding to / F = 0.985) and (b) M e ∼ 4.78 × 10 3 . We note that as M e increases, the instantaneous wave field η 0 (x, t n ) approaches its meanη 0 (x) at all times.

A. Transition to chaos
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 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 κ = 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 N 0 + 2 p impacts, where the first N 0 impacts are discarded to remove transient effects. We take the Fourier transform of the droplet position X n for the final 2 p 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 f 1 = 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 f 2 [see Fig. 10(b)]. The resulting instability is saturated by nonlinear effects, leading to the quasi-periodic stable wobbling motion with incommensurate frequencies f 1 and f 2 ≈ f 2 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 f 1 and f 2 .
For / F 0.98050, a third bifurcation (B3) yields the incommensurate frequency f 3 , 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 f 1 and f 2 (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.

B. The high memory limit
We now consider the high memory regime (M e 10 3 ) in which there is a qualitative change in the dynamics.  In the phase-plane plots (left column), the blue dots denote the prior 5000 impacts, and the red lines the final P impacts in the nonperiodic cases. The walker velocity V + n is normalized by the free walking speed V W at the corresponding value of / F . 15 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 V e (x) that influences the chaotic motion of the droplet (Sec. IV C). The additional notation used throughout this section is summarized in Table II.
To gain further understanding of the pilot-wave dynamics in this regime, it is useful to recast the iterative map (5)-(7) as a trajectory equation for the droplet position X n and the mean droplet velocity during flight U n ≡ X n+1 − X n . By computing the droplet fundamental matrix F(κ) analytically, the droplet's evolution may be expressed as where 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| 3 in Fig. 12), the full pilot-wave potential at each impact V p (x, t n ) oscillates in x. However, as |x| → ∞ the instantaneous wave field decays and we observe that V p (x, t n ) ≈ 1 2 Kx 2 for all time. In (15), K(κ) > 0 determines the strength of the time-averaged harmonic potential over one impact period and F(κ) > 0 prescribes the magnitude of the wave force (whose dependence FIG. 11. Route to chaos. (a) Fundamental frequencies f 1 , f 2 , and f 3 (dots) are introduced with successive bifurcations B2-B4 (gray) as the memory is progressively increased. The bifurcation B1 from stationary bouncing to orbiting occurs at / F ≈ 0.81 and is not shown in this figure. The periodic motion has period P, where P = 1/f 1 for stable dynamics (before bifurcation B2). After B2, the periodic orbit is unstable with linear instability frequency (2) (t)| between two trajectories X (1) and X (2) (whose initial position differ by a dimensionless distance of 10 −10 ) is shown for / F = 0.980593 (gray) and / F = 0.980594 (black). These values of correspond, respectively, to 3-frequency quasi-periodic motion and chaotic dynamics. onκ is weak). 51 The system (13) and (14) and the full pilotwave potential V p will be referred to throughout Secs. IV B and IV C.

The random walk dynamics
In Fig. 12, we plot the evolution of the full pilot-wave potential V p (x, t n ) and the corresponding droplet position X n 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 V p (x, t n ) at each impact, from which two important observations emerge. First, the minima far from the droplet (typically 1 Faraday wavelength away) remain at a roughly constant position over time, indicating FIG. 12. Pilot-wave dynamics at high memory (M e ∼ 1.17 × 10 4 withκ = 0.01). Droplet trajectory X n (red dots) and the full pilot-wave potential V p (x, t n ) (blue curves), which is the sum of the harmonic potential and the wave field, as defined in Eq. (15). The black squares denote the spatial minima of V p (x, t n ). 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 V p (x, t n ) on its path.
To analyze these dynamics, we define the set of turning times T ⊂ N to be the times at which the droplet changes direction. That is to say, if τ i ∈ T , then X (τ i ) is a local extremum and T i ≡ X (τ i ) is defined to be a turning point. In the droplet trajectory time-series data in Fig. 13(a), the turning times τ i and positions T i 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 D i = |T i+1 − T i |. Indeed, by plotting the distribution ρ D of distances D i [see Fig. 13(b)], it emerges that the distance between turning points is quantized, where ρ 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 TABLE II. Additional notation used throughout Sec. IV B.

Variable
Description Turning positions (raw data) Trend curve on the slow timescale X R n Residual droplet impact positions Probability distribution of X n (raw data) approximated by J 0 (k F x). From the observations in Fig. 12, it becomes clear that it is the minima of J 0 (k F x) 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 (x n ) n≥0 between turning points, where the jump distances d n ≡ |x n − x n−1 | are exactly restricted to d n ∈ D. 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 α ∈ R as a position visited by the jump process (whose possible values are determined in the following analysis), and without loss of generality, we set x 0 = α and consider x 1 > x 0 . Using the assumed structure of D = 0.6 + N, we define N n ∈ N such that d n = 0.6 + N n . As the droplet changes direction at each turning point, we observe that after an even number of jumps so |x 2n − x 0 | ∈ N for all n ∈ N. 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 x n ∈ M(α) for all n ≥ 0, where α 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 α ∈ [− 1 2 , 1 2 ). For symmetric statistics about the origin, we require that α be such that M(α) is also symmetric, which yields α ∈ {−0.3, 0.2}. For consistency with the jump distances d n ∈ D, the droplet may only leave each mesh point in a fixed direction (as depicted in Fig. 14), namely, to the right for x n ∈ α + Z and to the left for x n ∈ α + 0.6 + Z.

Detrending the long-time statistics
From the analysis in Sec. IV B 1, we expect the turning points T i (and the peaks of the corresponding probability distribution ρ T ) to be determined by the meshes  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 FIG. 14. Schematic diagram for a subset of mesh points for M 0 (top row) and M 1 (bottom row), where both meshes are periodic with period 1 and M 0 is a translation of M 1 by 1/2. The jump distance must lie in D, where turning points necessitate a change in direction after each jump. This evolution is equivalent to leaving each point in the direction of the arrow and changing color at each jump (blue/yellow). The relationship between the random walk dynamics and the derived effective potential V e (x) is evident in Fig. 18. Our study reveals an additional complication; specifically, the finite width about the peaks in ρ 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 M e , the timescale at which the global wave field structure changes. This drift obscures the structure of the underlying statistics induced from the shorttime dynamics; for example, there is only a weak structure apparent in the distribution of turning points ρ 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

Results
We now explore the statistical distributions following the detrending of the slow timescale dynamics. By defining R(t) ≡ X (t) − C(t), we have impact residuals X R n ≡ R(t n ) for all n ≥ 0 and turning point residuals T R i ≡ R(τ i ) for all τ i ∈ T , with respective residual probability distributions ρ XR and ρ 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η 0 (x) and an associated effective potential V e (x) to be defined in Sec. IV C.  Fig. 14). The droplet crosses the origin between crossing quadrants (CQ), in which either T R i > 0 and T R i+1 < 0, or T R i < 0 and T R i+1 > 0.
Examples of the corresponding residual distributions are given in Figs. 15(b) and 15(c), where the residual statistics are symmetric relative to mesh M 0 . The modes of ρ XR correspond to the modes of ρ 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 {±1.2, ±2.2, . . .} 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κ and M e and present the results in Fig. 16. When M e is fixed, the quantization is sharper when the waves dominate the harmonic potential [ Fig. 16(a)], but as κ 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 M 0 (see Fig. 14). However, it is relatively rare for the droplet to cross the centre of the bath (corresponding to T R i T R i+1 < 0), a feature that we rationalize in Sec. IV C.
When the wave memory M e is reduced (withκ fixed), the random walk-like dynamics shift to the  Fig. 17(b)], their presence is obscured in ρ XR [ Fig. 16(b)]. As M e is further decreased, the mesh points that counter the harmonic potential ({±0.7, ±1.7, . . .}) 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 (M 0 or M 1 ) is selected by the relative strength of the central and wave forces, where the mesh M 0 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 λ F /2.

C. The mean-pilot-wave potential
Based on the ideas of Theorem 1, we start by considering an effective potential V e (x) using the stationary residual distribution ρ XR (x) and the applied harmonic potential Remarkably, the direction associated with each mesh point (as given by Fig. 14) corresponds precisely to the gradients of V e , 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 timedependent probability distribution ρ(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 −D(U n )U n (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 |U n | 1 (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 W t over an infinitesimal timestep dt. This yields Langevin evolution equations for the position-velocity process where σ 0 > 0 prescribes the magnitude of the stochastic forcing. Here, we have defined the stochastic potential where ρ(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 η B * ρ (x, t) is for the time-dependent probability distribution and not for the stationary probability distribution. However, an initial condition ρ(x, 0) = δ(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 ρ s (x) were to exist [where ρ(x, t) → ρ s (x) as t → ∞], 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 ρ(x, t) = R p(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 η 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 σ 0 (u), will be the subject of future work. While the case without self-propulsion D(u) = D 0 is not appropriate for modeling the dynamics of walking droplets, we note that the stationary distribution ρ s (x) to Eq. (20) satisfies Kramer's equation for a given potential V, with implicit solution where E = σ 2 0 /(2D 0 ) and Z 0 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][16][17] and wobbling orbits. [13][14][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 ρ s (x) as a (nonlinear) Vlasov-Fokker-Planck equation, we can solve directly for ρ 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.  I (x, y) and Id is the identity operator. From the Neumann series ∞ n=0 L n = (Id − L) −1 , we recognize u B (x, y) as being proportional to the timeperiodic 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 ρ X (x) [ Fig. 20(a)]. Thus, for histogram bin centres ξ j with heights H(ξ j ), we have H(ξ j ) ≈ ρ X (ξ j ). For N X equally spaced points x i in the interval [−3, 3], we compute the mean wave fieldη C 0 (x i ) using the convolution (10), with the midpoint quadrature rulē η C 0 (x i ) ≈ δH j η B (x i − ξ j )H(ξ j ). To compare the simulated mean wave fieldη S 0 , we compute the mean squared error as δH → 0 [ 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(θ ) = 0, where θ = X 0 ,κ, G 1 , . . . , G N−1 and G (of dimension N + 1) is given below. The function G is dependent on several other functions of θ, which are computed at each step of the following algorithm. Hence, computation of the Jacobian ∂G/∂θ T requires an application of the chain rule, where the derivative ∂/∂θ T of each newly defined function is also computed. For an initial guess θ: 1. Useκ and (5) to uniquely find V + 0 (θ ) such that X 1 − X 0 = 0 (X 0 is an extremum). 2. Use the droplet iteration maps (5) and (7) with gradients G n and the initial conditions [X 0 , V + 0 (θ)] T to compute positions X n (θ ) and velocities V − n (θ) for n = 1, . . . , N. 3. For the wave field η 0 to satisfy the reflection conditions (12c) and (12d) with impacts X n (θ), use (6) to find the initial wave amplitudes a m (t 0 ; k, θ ) and a m (t + 0 ; k, θ ), which solve 4. Use (6) to recover the wave field η 0 (x, t n ; θ) and gradients g n (θ) = ∂ x η 0 (X n , t n ; θ). 5. Using g N (θ ) and V − N (θ) with (7), compute V + N (θ ). 6. For consistency with gradients and droplet reflection conditions (12a) and (12b), compute the output 7. If ||G(θ)|| ∞ < TOL, stop. Otherwise, update θ 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 (5)-(7) about the periodic state at times t + n for n = 1, . . . , N. By expressing all the perturbed variables at time t + n−1 as a single column vector, we construct (sparse) transition matrices T n ( ) to map the perturbed variables from t + n−1 → t + n for n = 1, . . . , N. The N-step stability matrix T is the product T = RT N · · · T 1 , 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 ⊂ T , 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 M e 1, we consider a linear combination of K basis functions b j (t), where KM e is the simulation duration. The trend function C(t) is thus given by the linear combination C(t) = K j=1 β j b j (t), where β j are the unknown coefficients. The penalty for over-fitting is chosen to minimize variation in the basis function coefficients β j , where the required smoothness is determined by the parameter θ > 0. (β j+1 − 2β j + β j−1 ) 2 .
Although methods exist to find the "optimal" value of θ for a given dataset, 60 it is sufficient for our purposes to simply fix θ = 500 for all datasets considered, where the residual statistics vary only weakly for 100 θ 1000.