The walking droplet system has extended the range of classical systems to include several features previously thought to be exclusive to quantum systems. We review the hierarchy of analytic models that have been developed, on the basis of various simplifying assumptions, to describe droplets walking on a vibrating fluid bath. Particular attention is given to detailing their successes and failures in various settings. Finally, we present a theoretical model that may be adopted to explore a more generalized pilot-wave framework capable of further extending the phenomenological range of classical pilot-wave systems beyond that achievable in the laboratory.

The surface of a vertically vibrated bath goes unstable to a pattern of standing waves when the vertical acceleration of the bath exceeds a critical value, known as the Faraday threshold. Below the Faraday threshold, a millimetric droplet may bounce on the surface of the bath indefinitely, provided a thin air layer is maintained between the drop and the bath during impact. In a small parameter regime of vibration frequencies and droplet sizes, this bouncing state may go unstable to a walking state in which the droplet is propelled by its underlying wavefield. We here review the hierarchy of analytic models of increasing complexity that have been developed and proven sufficient for rationalizing various walker behaviors. We further discuss how these models may be adapted in order to explore a more generalized pilot-wave framework capable of achieving quantum-like behavior beyond that accessible with the hydrodynamic pilot-wave system.

## I. INTRODUCTION

The surface of a fluid bath vibrating with vertical acceleration $\gamma (t)=\gamma sin2\pi ft$ goes unstable to a pattern of subharmonic standing waves when the forcing acceleration is increased beyond the Faraday threshold, $\gamma >\gamma F$.^{1,2} Below the Faraday threshold, a droplet may bounce indefinitely on the surface of the fluid bath, provided an air layer is maintained between the drop and the bath during impact.^{3} Coalescence will occur for forcing amplitudes $\gamma $ below a critical bouncing threshold $\gamma B$, while indefinite bouncing will occur for $\gamma >\gamma B$.^{4} The form of the drop’s vertical motion changes as $\gamma $ is increased. Different bouncing states are denoted by $(m,n)$, where $m/f$ is the total period of the bouncing motion, during which the droplet impacts the fluid $n$ times.^{5}

Just above the bouncing threshold, the droplet bounces with the same period as the bath’s oscillation, in a $(1,1)$ state. As the amplitude of the vibration is increased progressively, a bifurcation occurs, giving rise first to a (2,2) bouncer. As the amplitude is increased yet further, and for a certain range of drop sizes, the drop may settle into a fully period-doubled (2,1) bouncing state, wherein the drop impacts the bath precisely once every two driving periods. The vertical motion of the drop is then in resonance with the subharmonic Faraday wave field, leading to a significant increase in the amplitude of the surface waves generated. The transitions between the various bouncing modes were first reported by Protière *et al.*^{6} and Eddi *et al.*^{7} Their observations were then rationalized theoretically and extended to a wider range of system parameters by Moláček and Bush^{8} and Wind-Willassen *et al.*^{9}

As the vibration amplitude approaches the walking threshold, $\gamma W$, the bouncing state is destabilized by the underlying wavefield, causing the droplet to be propelled horizontally in a walking state characterized by rectilinear motion at a constant speed (Fig. 1).^{10} The drop’s free walking speed is determined by a balance between the propulsive wave force that is proportional to the gradient of the wavefield, $\u2207h$, beneath the drop, and a drag force that is linear in the drop speed. This force balance was first discussed by Protière *et al.*,^{6} who wrote a trajectory equation that allowed them to rationalize a bouncing to walking transition. Their trajectory equation was later revisited and further developed by Moláček and Bush,^{11} who gave explicit forms for both the drag term and the wave force in terms of the fluid parameters.

The walker system admits a range of complex dynamical phenomena, many of which were reviewed by Bush.^{12} The dynamics of a single droplet bouncing in the (2,1) mode may be understood using relatively simple models, since it is only necessary to characterize the wave field in the immediate vicinity of the drop. However, when multiple walkers interact, it is necessary to have a more complete wave model that captures both the local wave properties and its far field form. The interactions of multiple walkers (Fig. 2) may induce substantial variations in the vertical dynamics of the droplets, an effect that has been incorporated into several recent studies of droplet pair interactions.^{13–15} Although the wave form generated by drop impact is well approximated by a linear wave,^{16} the coupling between the vertical bouncing dynamics and the wave field is a complex nonlinear problem, which has only recently been addressed systematically.^{17,24,25}

In Sec. II, we review the various models that have been developed for the drop’s vertical and horizontal motion. In Sec. III, we discuss the hierarchy of wave models that have been developed to account for various dynamical phenomena observed for both single droplets and multiple interacting droplets. In Sec. IV, we review recent modeling attempts to describe the interaction of walkers and submerged barriers, which arise in a number of hydrodynamic quantum analogs. Finally, in Sec. V, we present a generalized pilot-wave framework that may be adopted in exploring this new class of dynamical systems.

## II. DROP DYNAMICS

### A. Vertical dynamics

Beyond the bouncing threshold, $\gamma >\gamma B$, the drop bounces indefinitely on the surface of the bath. The bouncing behavior was characterized experimentally by Protière *et al.*,^{6} who presented a regime diagram of the various bouncing modes observed for 50 cSt oil, followed by Eddi *et al.*^{7} who considered 20 cSt oil. In the bouncing regime, the bath and drop never come into contact and are always separated by a thin air layer.^{4} In order to understand the vertical motion of the droplet, a hierarchy of models have been developed. A bouncing droplet typically has a diameter of approximately 1 mm; thus, curvature pressure dominates hydrostatic pressure in its contribution to the reaction force applied by the bath during impact. In order to remove the complications associated with the influence of the bath inertia, Gilet and Bush^{5,20} investigated droplets bouncing on a soap film. They found that as the film distorts during impact, it reacts to the drop like a spring, exerting a force that increases linearly with intrusion depth for small deformations of the film, the spring constant being linearly proportional to the surface tension $\sigma $. We note that the low energy impact of small drops on superhydrophobic surfaces may also be described in terms of a linear spring owing to the form of the increase of surface energy with drop distortion.^{21,22}

Moláček and Bush^{8} investigated the use of a linear spring model to characterize the interactions between a drop and a vibrating fluid bath. Consider a droplet of mass $m$ bouncing on the surface of a bath vibrating with vertical acceleration $\gamma (t)=\gamma sin\omega t$. Suppose further that there is an underlying wavefield, whose height directly beneath the drop is given by a function $h(t)$. The height of the fluid surface beneath the drop is $a(t)=\u2212\gamma \omega 2sin\omega t+h(t)$, which is a superposition of the shaking of the bath and the underlying wavefield. One then transforms into a reference frame in which the surface of the bath beneath the drop is stationary so that the vertical position of the drop is denoted by $Z(t)$. The origin $Z=0$ is defined as the point when the drop’s center of mass is one radius $R$ above the equilibrium surface of the bath. Moláček and Bush^{8} then used a linear spring model of the form:

where $g\u2217(t)=g+d2adt2$ is an effective gravity that incorporates the influence of the vibrating frame of reference and the underlying wavefield and $H(x)$ is the Heaviside function. $D$, the effective damping constant, and $C$, the spring constant, are the fitting parameters, whose values were inferred by measuring the contact time and coefficient of restitution of drop impacts in the parameter regime of interest. One expects the linear spring model to break down when inertial forces become comparable to the effects of surface tension, i.e., when $We=\rho RV2/\sigma \u226b1$, where $V$ is the impact speed of the droplet. Terwagne *et al.*^{23} described the droplet-bath system in terms of two masses coupled by a damped spring and captured transitions between bouncing modes reminiscent of those observed experimentally.

The linear spring model gives an adequate description of the various bouncing transitions that may occur as the forcing amplitude of the bath is varied. However, Moláček and Bush^{8} found that unrealistic values of $C$ and $D$ were required in order to achieve an acceptable fit with experiments over the entire parameter regime considered. This motivated the authors to adapt their quasi-static model of drop impact on a rigid substrate^{22} to derive a logarithmic spring model that was found to describe more satisfactorily the effect of bath and drop deformations on the vertical dynamics. The logarithmic spring has been adopted by Milewski *et al.*^{24} in their relatively comprehensive numerical model of the walker system, as was used by Galeano-Rios *et al.*^{60} to rationalize the self-propulsion of unequal droplet pairs moving along their line of centers in the so-called *ratcheting* mode.^{7}

The most recent model for the vertical dynamics was developed by Galeano-Rios *et al.*^{25} The authors model the impact of a droplet on the surface of a fluid bath by coupling the free-surface Navier-Stokes equations to the motion of a hydrophobic sphere, whose contact area with the bath is updated dynamically. Their approach allows them to eliminate fitting parameters from their model for the vertical dynamics and so obtain a fully predictive model for the transitions between the wide variety of bouncing states reported by Wind-Willassen *et al.*^{9}

In order to couple the vertical dynamics to the wave form, it is necessary to have a model for the waves, $h(t)$, beneath the droplet. For non-resonant bouncers, this wavefield has a sufficiently small amplitude that its influence on the bouncing dynamics may be negligible.^{8,11} However, for bouncers and walkers in the resonant (2,1)-mode, the amplitude of this wavefield is significantly larger, resulting in a complex coupling between the vertical dynamics and the wave form. The effects of this coupling were found to be negligible in early studies of single walkers,^{18,26,27} but more recent studies^{17} have shown a significant effect of variations in the vertical dynamics for interacting droplets. The manner in which such variations may be treated will be detailed in Sec. III.

### B. Horizontal dynamics

Consider the trajectory of a drop of mass, $m$, bouncing periodically at exactly twice the driving period of the bath. Provided the drop is sufficiently small relative to the wavelength of its Faraday wave field, the drop may be approximated as a point particle following a trajectory $xp(t)=(xp(t),yp(t))$ on the surface of the bath. The drop’s horizontal motion is driven by a propulsive wave force and resisted by drop inertia and drag forces induced during impact and flight. This horizontal force balance was first modeled theoretically by Protière *et al.*,^{6} who considered a one-dimensional trajectory equation of the form:

This equation describes the strobed dynamics, meaning that the forces have been averaged over a period of the drop’s vertical motion. The dominant drag was assumed to arise from the shear drag generated by the thin air layer created during impact and so to scale as $fV\u223c(\mu Gs/hf)(\tau /TF)$, where $\mu G$ is the air viscosity, $s$ is the approximate area of contact, and $hf$ is the approximate air layer thickness. $Fb\u223cm\gamma (A/\lambda F)(\tau /TF)$ was taken to be the effective force on the drop due to the inclined wave field, where $A$ is the wave amplitude and $\tau $ is the contact time between the drop and the bath. Protière *et al.* estimated the local slope of the surface in terms of the ratio of the propagation speed of the drop, $x\u02d9p$, and the phase speed of the waves, $VF$, in which case the force exerted by the wave field on the drop is proportional to $sin2\pi x\u02d9pVF\phi $. By varying $\gamma $, and thus $Fb$, this model captures a key feature of the system, namely, a bouncing to walking transition. However, the wave force in this model does not depend on the history of the drop’s horizontal motion, preventing it from capturing more complex dynamical states reliant on the memory of the system on the droplet’s past trajectory.

The horizontal dynamics of a bouncing drop were also considered by Shirokoff.^{28} His trajectory equation incorporated a linear drag and a wave forcing proportional to the gradient of the wave field. He assumed that the impacts could be treated as instantaneous and derived the wave form from the linearized capillary-gravity wave equation. However, he neglects the effects of the accelerating frame of reference in the modeling of the wave field and so could not deduce the dependence of the wave amplitude on the forcing acceleration of the bath, $\gamma $. As a result, this model contained several fitting parameters whose values were not prescribed by experiments.

The trajectory equation was revisited by Moláček and Bush.^{11} They found that the dominant contribution to the drag comes from the momentum transfer between the drop and the bath during impact and from the Stokes drag imparted by the air during flight, while the lateral force component is proportional to the gradient of the wavefield, $\u2207h(xp,t)$, at the drop location. This leads to a trajectory equation of the form:

where $m$ is the mass of the drop. $D(t)=C\rho R/\sigma F(t)+6\pi R\mu a$ is the time-dependent drag, whose first term models the influence of momentum transfer between the droplet and the bath during impact, and whose second term is an approximation to the air drag induced during flight. $F(t)=mZ\xa8+mg\u2217(t)$ is the normal force exerted on the drop during contact, which may be calculated from one of the spring models discussed in Sec. II A. Moláček and Bush^{11} pointed out that in order for the drop to bounce indefinitely, it is necessary for $\u222b0TFF(\tau )d\tau =mgTF$ so that the average reaction force of the bath balances the weight of the drop. After determining the full temporal dependence of the drag and wave force in (3), Moláček and Bush restricted attention to resonant $(2,1)$ bouncers, and averaged the full trajectory equation over the bouncing period, to obtain the trajectory equation:

where $D\xaf=1TF\u222b0TFD(t)dt=Cmg\rho R/\sigma +6\pi R\mu a$. $C$ is a fitting parameter with tight bounds, $0.17<C<0.33$, deduced experimentally by Moláček and Bush^{11} through consideration of drop impacts in the parameter regime of interest. Provided $F(t)$ is known, the vertical bouncing dynamics may be eliminated from consideration, and this time-averaged trajectory equation describes the mean horizontal dynamics, which may be viewed experimentally by strobing the system at the Faraday frequency. It is thus referred to as the *strobed* trajectory equation.

## III. WAVE MODELS

The wave forms generated by impacts on a fluid bath have been examined experimentally in various studies. Eddi *et al.*^{19} showed experimentally that a sphere impact on a vibrating bath produces an outward propagating capillary wave that triggers a pattern of standing waves in its wake. This observation motivated the theoretical description of the wave form as a sum of standing waves produced at each impact location. They also studied the wave generated by a walking droplet and observed that the wave form has a clear Doppler effect due to the interference of waves generated at successive impacts. Damiano *et al.*^{16} applied the free-surface Schlieren technique of Moisy *et al.*^{29} in an experimental study of the wave forms excited by bouncing drops. By placing a pattern of random dots at the bottom of the bath, one then measures the refraction of the dots due to the topography of the liquid surface in order to reconstruct the form of the surface waves using a digital image correlation algorithm. Damiano *et al.*^{16} used this technique to examine the spatial decay of the waves generated by a bouncing drop and measured the motion of the radially propagating capillary waves generated at impact.

The models for both the vertical (1) and horizontal dynamics (4) require an accurate description of the underlying wavefield, $h(x,t)$. Consequently, in order to rationalize phenomena such as the transition between bouncing and walking, a hierarchy of wave models has been developed. In all analytic models of the wavefield, it is assumed that the drop is in a $(2,1)$ mode so that it is bouncing in resonance with the subharmonic Faraday wave.

Guided by their experimental observations, Eddi *et al.*^{19} proposed a phenomenological model for the wavefield. The wave after $N$ bounces at impact locations $xn$ was written as a linear superposition of circular waves with the Faraday wavenumber $kF$:

The waves decay exponentially in space over a length-scale $\delta $ that was determined experimentally. The waves decay in time over a time scale $\tau =MeTF$, where $Me$ is a dimensionless memory parameter, which prescribes the proximity to the Faraday threshold, $Me\u221d1\u2212\gamma /\gamma F\u22121$. As shown in the Appendix, this temporal damping rate may be obtained asymptotically by considering the damped Mathieu equation close to the onset of instability. The wave amplitude $\xc5$ is effectively a free parameter, since this wave model does not account for the forcing from the drop. Notably, the wave form (5) is singular at the origin and thus cannot be used to rationalize the bouncing to walking transition. Nonetheless, it captured the essential features of the walker wave form in the far field and introduced the important concept of the *memory* in determining the influence of the droplet’s past on its present. This waveform has been used in conjunction with a discrete iterative model, developed by Emmanuel Fort, in which the horizontal momentum transfer between the drop and the bath is calculated at each impact from the droplet’s speed and the slope of the wave surface. This iterative model has been used to capture a number of walker behaviors, including single-particle diffraction,^{30} orbits in a rotating frame^{31} and harmonic potential,^{32} and the so-called promenading mode, wherein two droplets walk side by side while oscillating laterally along their line of centers.^{33} The details of this model are provided in Appendix A of Labousse *et al.*^{34} and Appendix A of Borgeshi *et al.*^{33}

Moláček and Bush^{11} developed from first principles an improved expression for the wave form, which coupled the vertical dynamics to the wave amplitude beneath the drop. The wave form was expressed in terms of a linear combination of the waves generated at successive impacts, but with a spatial dependence proportional to a zeroth order Bessel function of the first kind:

Although this expression does not contain the exponential spatial damping factor of Eq. (5) important in the far field, it is regular at the origin and so may be used in conjunction with the strobed trajectory equation (4) to rationalize the bouncing to walking transition. The wave amplitude $A~$ was given in terms of known fluid parameters. The memory parameter $Me$ is proportional to the ratio of the viscous decay time of the waves in the absence of vibrational forcing, $Td$, and the Faraday period, $TF$, and is also inversely proportional to the proximity to the Faraday threshold, specifically

We note that the divergence of the memory parameter as $\gamma \u2192\gamma F$ does not indicate the onset of singular behavior but rather signals the limitations of linear wave models in describing this high memory limit.

In Eq. (6), $S$ is one of the two *phase* parameters:

Here again, $F(t)=mZ\xa8+mg\u2217(t)$ is the normal force exerted by the drop on the bath. $S$ and $C$ arise in the wave theory as the first Fourier coefficients in a Fourier series expansion of $F(t)$ divided by the weight of the drop. To gain physical insight into the origin of these parameters, notice that in the limit of an instantaneous impact, $S$ and $C$ are proportional to the sine and cosine of the impact phase relative to the bath acceleration. In general, $S$ determines the effectiveness of the drop in generating waves, while $C$ appears in the strobed trajectory equation (4) when one performs the average $F(t)\u2207h(xp,t)\xaf$ and so determines the effectiveness of the waves in propelling the droplet horizontally. We note that in the strobed trajectory equation (4), the phase parameters appear only as a product, $sin\Phi =2SC$. Moláček and Bush^{11} gave approximate bounds for this product as $0.25<sin\Phi <0.65$, noting that $sin\Phi $ will in general depend on the value of the memory and the fluid parameters. In their study, they used $sin\Phi \u22480.5$ as a fitting parameter chosen to best match the free walking speed for silicon oil with kinematic viscosity $\nu e=20$ cSt.

The discrete nature of the wavefield in Eq. (6) makes study of the trajectory equation analytically intractable for any situation more complex than the single droplet bouncing-to-walking transition. To this end, the trajectory equation developed by Moláček and Bush was transformed by Oza *et al.*^{18} into an integro-differential form. Oza *et al.* demonstrate that, provided the time scale of bouncing $\u223cTF$ is much faster than the time scale of horizontal motion, the sum may be replaced by an integral. The wave form may then be written as an integral over a continuous history of past impacts:

Note that the $1t\u2212s$ factor appearing in the discrete model (6) is replaced in (9) by the constant $TF\u22121/2$ for the sake of analytic expedience, an approximation that was found to have a negligible effect on the model predictions. By inserting the wave form (9) into the trajectory equation (4), one obtains the *stroboscopic* model, so called because the forces have all been averaged over a period of the drop’s vertical motion and the equations contain no notion of the bouncing time scale. In the stroboscopic model, the phase parameters enter as a product $sin\Phi =2SC$, which was taken by Oza *et al.*^{18} to be a constant of approximately 0.3. The lower value used in this study may be attributed to the replacement of the $1t\u2212s$ factor with a constant, $TF\u22121/2$. This approximation alters the amplitude of the waves, thus leading to different inferred values of the phase parameters, as is discussed in detail by Couchman *et al.*^{17} With the simplification of a constant $sin\Phi $, the stroboscopic model was able to accurately predict the bouncing to walking transition and the dependence of the free walking speed on memory.^{18} It also provided theoretical rationale for the stability of circular orbits observed when the bath is subjected to rotation^{26,31,35} or a simple harmonic potential.^{27,32,36–38} Finally, the stroboscopic model was adopted in the first theoretical study of the transitions to chaos in orbital pilot-wave dynamics.^{39}

The success of the stroboscopic model has been limited in rationalizing the stability of configurations with multiple droplets. Oza *et al.*^{13} have recently studied the interactions of two identical droplets above the walking threshold $\gamma W$, where the droplets may settle into orbits with quantized radii, as observed by Couder *et al.*^{10} and Protière *et al.*^{6,40,42} Two additional refinements to the stroboscopic model were required in order to achieve a match between experimental results and theoretical predictions.

The wavefield introduced by Moláček and Bush was found to have insufficient spatial decay to account for the long range interactions of the drops. Consequently, a more accurate spatiotemporal damping was added to the wave kernel of Oza *et al.*,^{18} leading to the following integral expression for the waves:

The spatiotemporal damping term was alluded to, but not fully developed, in Appendix A of Moláček and Bush^{11} [see Eq. (A47)]. We outline a derivation of its form here in the Appendix. Damiano *et al.*^{16} found that when this additional spatial damping is included, the theoretically predicted wave form provides an excellent description of that measured using the free-surface synthetic Schlieren technique of Moisy *et al.*^{29} It should be noted that a wave form with a similar spatiotemporal damping is derived by Tadrist *et al.*^{43} from a fluid mechanical model obtained using an alternative approximation to the Navier-Stokes equations.

Although the assumption that $sin\Phi $ is constant has proven to be sufficient for describing the motion of a single drop, it has limitations when describing the interactions of multiple drops, where modulations of the vertical dynamics are known to arise. The effect of variations in vertical dynamics was seen to be important in both orbiting pairs^{13} and the promenading mode.^{14,33} Oza *et al.*^{13} proposed a linear dependence of $sin\Phi $ on the effective forcing $\Gamma eff=\gamma /\gamma F\u2212phc$, where $p$ is a fitting parameter and $hc$ is the local wave amplitude produced by both drops. The inclusion of spatiotemporal damping and phase modulations was essential in capturing the stability of bound orbits of different radii. Arbelaiz *et al.*^{14} were also obliged to consider variations in $sin\Phi $, deducing a quadratic dependence on an alternative effective forcing $\Gamma eff=\gamma /\gamma F\u2212p~h~$, where $p~$ is another fitting parameter and $h~$ is the component of the wave beneath one drop due to the wave produced by its partner. Consideration of both vertical dynamics and spatial damping was required by Couchman *et al.*^{15} in their study of the instabilities and oscillations of bound droplet pairs, but here the dependence on the vertical dynamics was treated systematically using the procedure described by Couchman *et al.*^{17}

A drawback of the stroboscopic model is that it has been derived from purely a standing wave model. This assumption was partially justified by Moláček and Bush^{11} who showed that waves with wavenumber very close to the Faraday wavenumber $kF$ have the smallest temporal decay factor (Fig. 2). However, there are other subdominant wave modes with different wave numbers, the inclusion of which allows the full wave field to be a written as a sum of out-of-phase standing waves. It is possible to include these additional waves within the framework of the stroboscopic model, as will be discussed in a future work. Since the radial motion of the zeros of the wave field does not affect the dynamics of single drops, a standing wave model was sufficient to rationalize the bouncing to walking transition. However, the radially propagating wave front generated by each impact has been seen by Galeano-Rios *et al.*^{60} to play a critical role in the reversal in the direction of ratcheting droplet pairs^{7} as the memory is increased. We thus expect that replacing the monochromatic Bessel function by a sum of Bessel functions with different wave numbers may alter the stability of droplets interacting through their common wave field.

Milewski *et al.*^{24} introduced a more complete Faraday pilot-wave model by coupling the logarithmic spring model of Moláček and Bush to a quasi-potential, weakly viscous wave model. Describing the fluid velocity as a potential flow $u=\u2207\varphi $ and linearizing the incompressible Navier-Stokes with small viscosity, one arrives at the following quasi-potential, weakly viscous system for the wave height $h(x,t)$:

The drop is coupled to the wave evolution via the pressure term $PD=PD(x\u2212xp(t),t)$, which vanishes when the drop is not in contact with the bath, and is otherwise given by a ratio of the normal force exerted by the drop and the contact area between the drop and the bath. This model is able to capture many more subtle features of the walker system, such as modulations in vertical dynamics and the traveling wave fronts reported by Eddi *et al.*,^{19} but neglected in the earlier standing wave models. Their model also does not require the drop to be in a $(2,1)$ mode, thus providing a theoretical framework for examining drops in more complex bouncing modes. Furthermore, through its relatively complete treatment of the wave field, this was the first theoretical model to capture the Doppler effect in the wave field reported by Eddi *et al.*^{19} The wave forms of a bouncing droplet predicted by this model and the stroboscopic model with spatiotemporal damping were found to agree favorably, as demonstrated in Fig. 9 of Milewski *et al.*^{24} Comparisons between the bouncer wave profiles predicted by Milewski *et al.*,^{24} the stroboscopic wave with and without spatial damping^{11,18} and experiments are presented in Figs. 4–6 of Damiano *et al.*^{16}

Simulations of the model of Milewski *et al.*^{24} are numerically expensive. However, Durey and Milewski^{38} derived a reduced model based on Eqs. (11)–(14), by assuming that the drop is in a resonant $(2,1)$ bouncing mode, that the impact is instantaneous, and that the drop contacts the bath at a single point. It is then possible to derive a discrete-time system for the drop’s horizontal position. This numerical model is highly efficient and allowed Durey and Milewski to explore the long-time statistics of a drop moving chaotically in a harmonic potential. Furthermore, it allows for a full stability analysis of different bound states, including promenading and orbiting droplet pairs, and droplet trains.

In the limit of weak acceleration, Bush *et al.*^{44} showed that it is possible to re-express the wave force in terms of a wave-induced added mass and a nonlinear drag. Specifically, taking the weak acceleration limit of the stroboscopic model without spatiotemporal damping or phase variations, the trajectory equation (4) may be reduced to a nonlinear ordinary differential equation, known as the *boost* model

where the form of the hydrodynamic boost factor $\gamma B$ and restoring force $Dw(|x\u02d9|)x\u02d9$ were derived explicitly and $F$ is an externally-imposed, classical force. This equation describes the dynamics of a particle with speed-dependent mass subject to a nonlinear drag force, which acts to drive the particle motion toward its free walking speed. The nonlinear, speed-dependent drag force was first derived by Labousse and Perrard,^{45} who showed that the evolution of walkers in a harmonic potential at low memory can be well described by a two-dimensional Rayleigh oscillator. The boost formulation of the walker dynamics provided insight into the orbital motion of walkers in both a rotating frame and in the presence of a confining spring force:^{26,31,35} the wave-induced added mass accounts for the anomalously large orbital radii observed.^{44} Moreover, when the walker moves at its free walking speed, the nonlinear drag vanishes. The driving of the bath then precisely balances the viscous dissipation, and the boost equation (15) describes the inviscid dynamics of a particle whose mass depends on its speed.

## IV. INTERACTIONS WITH BOUNDARIES

The wave models presented in Sec. III are strictly valid only for drops moving in free space. However, many interesting phenomena arise when walkers interact with abrupt changes in bottom topography. Examples include walker interactions with submerged pillars (see Harris *et al.*,^{61}), walker interactions with a single and double slits,^{30,46,47} unpredictable tunneling of walkers across shallow regions,^{48,49} and the wavelike statistics of walkers in circular corrals.^{50,51}

Couder and Fort^{30} were the first to numerically investigate the effect of submerged boundaries on the trajectory of a walking droplet in the context of single and double slits. They extended the trajectory equation of Protière *et al.*^{6} to include boundary effects, by modeling the walls as secondary wave sources. The phase of the secondary wave is selected so that the total wave amplitude vanishes on the boundaries, effectively imposing a Dirichlet boundary condition. In the context of walker motion in corrals, Gilet^{52,53} and Blanchette^{54} also treated abrupt variations in depth as effective boundaries. In both studies, the waves are decomposed into a sum of eigenfunctions of the Laplacian on the finite domain of interest. Gilet^{52,53} imposed Neumann boundary conditions, while Blanchette^{54} imposed a Dirichlet condition.

Dubertrand *et al.*^{55} considered the behavior of a walker passing through a slit with a theoretical model similar to Gilet,^{52,53} using Neumann boundary conditions to model effective boundaries, but constructed the wave out of the Green’s function for the Helmholtz equation obeying a far field radiation condition. The models of Gilet,^{52,53} Blanchette^{54}, and Dubertrand *et al.*^{55} were each able to capture certain behaviors of walking droplets interacting with submerged boundaries. However, the use of a Dirichlet or Neumann boundary condition to model the effect of the abrupt change in bottom topography is an assumption that cannot be derived from the fluid mechanical boundary conditions. It is thus not entirely surprising that these models could not capture the full range of behaviors observed in closed geometries.

Nachbin *et al.*^{49} studied tunneling of walkers across submerged boundaries. By restricting their attention to motion in a one-dimensional channel, the authors could use conformal mapping techniques to transform the variable topography onto a flat bottom. Doing so allowed them to reproduce the experimentally observed decay of tunneling probability with barrier width.^{48} Their method is highly appealing for 1D geometries, but its extension to two-dimensional geometries is not straightforward.

Faria^{56} has adapted the quasi-potential model of Milewski *et al.*^{24} to study interactions with submerged barriers. Rather than treating changes in topography as boundaries, he treats them as regions where the local wave speed changes. He thus approximates the vertical derivative of the velocity potential at the surface, $\varphi s(x,0,t)$, in Eq. (14) as

where $b(x)$ is chosen so as to match the local, depth-dependent dispersion relation in all regions of the bath. The model of Milewski *et al.* then reduces to a system of two partial differential equations for the wave form $h(x,t)$ and the velocity potential $\varphi (x,0,t)$ evaluated at the bath surface:

The vertical dynamics are approximated by instantaneous, point impacts. This reduced model was able to capture the non-specular reflection and refraction of walkers approaching discrete steps^{56,57} and circular pillars (Harris *et al.*,^{61}), in addition to the interaction of walkers with slits.^{46,56}

## V. GENERALIZED PILOT-WAVE FRAMEWORK

We proceed by discussing a generalized pilot-wave framework that may be adapted in exploring pilot-wave dynamics in parameter regimes not necessarily accessible with the walker system. Two critical features of the walker system responsible for the quantum-like behavior of the walker system are a monochromatic wavefield and resonance between the drop and the bath. We thus retain these key features while discarding others.

Non-dimensionalizing lengths by $kF$, and time by $TFMe$, the stroboscopic trajectory equation of Oza *et al.*^{18} including spatiotemporal damping and phase variations may be written as

with

When the trajectory equation is written in the form of (19)–(22), we can interpret $\kappa $ as an effective mass of the particle prescribing its inertial response to the wave forcing, whose strength is prescribed by $\beta $. The central feature of the wave form in (20) is a radially symmetric, monochromatic wave kernel, generated at each point along the walker’s trajectory and decaying in time over a time scale $TFMe$. Since the waves in the walker system are in two dimensions, the basic spatial dependence is a monochromatic Bessel function, $J0(r)$. Of course, the trajectory equation (19) may be augmented to include additional classical forces such as the Coriolis force,^{26,39,58} Coulomb force,^{39} or linear spring force.^{32,37,38}

Bush^{12} discussed a relatively simple generalized pilot-wave framework that follows from (19) and (20) if one sets the phase parameters to be constant and the damping kernel to unity. In this limit, for the walker system, Oza *et al.*^{18} showed that $\kappa $ and $\beta $ depend explicitly on the forcing acceleration of the bath $\Gamma $ and hence also the memory $Me$:

In this limit, there is a two-dimensional parameter space $(\beta ,\kappa )$, exploration of which has already yielded some new phenomenology (see Valani *et al.*,^{62}). The droplet system is constrained to a relatively small regime, where $\kappa $ and $\beta $ are $O(1)$, but numerical investigation of the trajectory equation across a larger parameter regime may yield other behavior that is not accessible in the walker system. In particular, Oza *et al.*^{63} consider hydrodynamic spin states, which consist of a droplet trapped in a circular orbit by its own wave field. While such states may appear as transients in the walking drop system,^{27} Oza *et al*.^{63} have shown that in the small $\kappa $ limit of this generalized pilot-wave framework, they are stable.

Spatio-temporal damping in the context of the walker system refers to the addition of a Gaussian profile of time dependent width to the wave kernel:

This specific form of spatiotemporal damping is special to the walker system, since it is derived through an approximation of the effect of viscous damping on the far field of the wave, with the damping factor $\alpha $ given in terms of fluid parameters in Table I. It is an inessential feature of the pilot-wave framework when considering the behavior of isolated walkers, since its inclusion does not substantially alter the local wavefield. However, it has proved essential for rationalizing the interactions of bound droplet pairs.^{13–15} While still an approximation to the true wavefield, the wave form with spatiotemporal damping included matches well with experimental data,^{16} suggesting that it is sufficient for most walker systems.

Variable | Description |

$R$, $m$ | Droplet radius, mass |

$\rho $, $\sigma $, $\nu e$ | Fluid density, surface tension, effective kinematic viscosity |

$g$ | Gravitational acceleration |

$f$, $\omega =2\pi f$ | Frequency and angular frequency of bath oscillations |

$\mu a$ | Dynamic viscosity of air |

$\gamma $, $\Gamma =\gamma /g$ | Amplitude of bath acceleration, dimensionless amplitude |

$kF$, $TF=4\pi /\omega $ | Faraday wavenumber, Faraday period |

$D\xaf=Cmg\rho R/\sigma +6\pi R\mu a$ | Drag coefficient |

$S,C,sin\Phi =2SC$ | Phase parameters (see Sec. III) |

$\omega D=\sigma \rho R3$ | Droplet’s natural oscillation frequency |

$Td=\nu ekF2$ | Decay time |

$Bo=\rho gR2/\sigma $ | Bond number |

$Ohe=\nu e\rho (\sigma \rho R)\u22121/2$ | Ohensorge number |

$\Omega =2\pi f\rho R3/\sigma $ | Vibration number |

$Me=TdTF(1\u2212\gamma /\gamma F)$ | Memory parameter |

$A=2\pi \nu eTF\pi mgkF3R2\sigma (3kF2R2+Bo)$ | Wave amplitude in stroboscopic model^{18} |

$\kappa =mD\xafTFMe$ | Non-dimensional mass |

$\beta =mgAkF2TFMe2D\xaf$ | Non-dimensional memory force coefficient |

$\alpha =\u03f522\nu e(1+2\u03f52)$, | Spatial damping factor |

$\u03f5=Ohe\Omega RkF3kF2R2+Bo$ | Viscosity induced wave-number correction^{13} |

Variable | Description |

$R$, $m$ | Droplet radius, mass |

$\rho $, $\sigma $, $\nu e$ | Fluid density, surface tension, effective kinematic viscosity |

$g$ | Gravitational acceleration |

$f$, $\omega =2\pi f$ | Frequency and angular frequency of bath oscillations |

$\mu a$ | Dynamic viscosity of air |

$\gamma $, $\Gamma =\gamma /g$ | Amplitude of bath acceleration, dimensionless amplitude |

$kF$, $TF=4\pi /\omega $ | Faraday wavenumber, Faraday period |

$D\xaf=Cmg\rho R/\sigma +6\pi R\mu a$ | Drag coefficient |

$S,C,sin\Phi =2SC$ | Phase parameters (see Sec. III) |

$\omega D=\sigma \rho R3$ | Droplet’s natural oscillation frequency |

$Td=\nu ekF2$ | Decay time |

$Bo=\rho gR2/\sigma $ | Bond number |

$Ohe=\nu e\rho (\sigma \rho R)\u22121/2$ | Ohensorge number |

$\Omega =2\pi f\rho R3/\sigma $ | Vibration number |

$Me=TdTF(1\u2212\gamma /\gamma F)$ | Memory parameter |

$A=2\pi \nu eTF\pi mgkF3R2\sigma (3kF2R2+Bo)$ | Wave amplitude in stroboscopic model^{18} |

$\kappa =mD\xafTFMe$ | Non-dimensional mass |

$\beta =mgAkF2TFMe2D\xaf$ | Non-dimensional memory force coefficient |

$\alpha =\u03f522\nu e(1+2\u03f52)$, | Spatial damping factor |

$\u03f5=Ohe\Omega RkF3kF2R2+Bo$ | Viscosity induced wave-number correction^{13} |

The phase parameters are more subtle but provide a means to couple the oscillations of a particle to those of its spatially extended wavefield. $S$ describes how effective the particle oscillations are in generating waves at each impact, while $C$ describes how effectively the waves impart horizontal momentum to the particle. For the droplet system, Couchman *et al.*^{17} find that, in the single droplet walking regime, $S$ is approximately constant, while $C$ may be approximated as a linear function of $\Gamma $ and the local wave amplitude $h~=h(xp(t),t)$. For an isolated walker, the product $1/2sin\Phi =SC$ is bounded between $0.1$ and $0.4$, which accounts for the successes of the stroboscopic model when $sin\Phi \u22480.5$ was adopted.^{11} One thus expects that including phase modulations will have a negligible effect on single droplet dynamics. However, both parameters $S$ and $C$ may become time-dependent when multiple walkers interact. Including the effect of phase modulations has been essential in several recent studies of droplet pair interactions.^{13–15}

The progressive theoretical modeling of this subtle hydrodynamic system suggests a number of new directions to explore with a generalized pilot-wave framework. For a droplet bouncing on a fluid bath, the wave-field can only depend on two spatial dimensions, giving rise to a wave force in the direction of its gradient. One could then first extend the trajectory equation to three dimensions by replacing the Bessel function with the spherically symmetric eigenfunctions of the Laplacian in higher dimensions. The form of spatiotemporal damping in the wave kernel of Eq. (20) arises due to the action of viscous dissipation on the waves. Within the context of a generalized pilot-wave framework, the form of the spatial damping kernel may be prescribed by different physical effects and thus take on different functional forms, a direction which has yet to be fully explored.

Another direction in which the generalized pilot-wave framework (19)–(22) could be extended would be the consideration of more complex couplings between the particle and the wave. Fort and Couder^{59} considered a deviation from the hydrodynamic walker system, by considering *inertial* walkers, wherein the particle generates a wave at each impact that continues to translate with the velocity of the particle at the time of its generation. These inertial walkers were found to admit circular orbital trajectories, whose radii satisfy the Bohr-Sommerfeld quantization condition. One expects that explorations of a generalized pilot-wave framework may yield a new array of dynamical phenomena not necessarily achievable with the walker system.

## VI. DISCUSSION AND CONCLUSIONS

We have reviewed the various models for the vertical and horizontal dynamics of a droplet bouncing on the surface of a vertically vibrated bath. In particular, we have considered the hierarchy of wave models that have been developed to capture the range of complex dynamical phenomena observed in the walker system. For describing the dynamics of single droplets in the resonant (2,1) bouncing mode in unbounded domains, the relatively simple stroboscopic models derived by Fort^{33,34} and Oza *et al.*^{18} suffice. Specifically, they have provided an adequate description of the bouncing to walking transition and the dynamics of circular orbits in the presence of an externally imposed classical potential. The wave model of Milewski *et al.*^{24} provides a more comprehensive framework in which to study drops in other bouncing modes or drop-drop interactions. The extension of this model to a discrete form by Durey and Milewski^{38} provides an efficient method for performing stability analysis of more complex orbits and interacting droplet pairs. In order to describe walkers interacting with submerged boundaries, simple analytic models for the wave field have been found lacking, but numerical methods^{49,56} have proven to be successful.

To accurately model interactions between walkers using analytic wave models, it is essential to consider both the effects of phase modulations and the spatiotemporal wave damping.^{13–15} The spatial damping in the walker system arises from the specific nature of the viscous damping. Phase modulations, however, arise due to a nonlinear coupling between dynamics on the scale of the particle and the wavefield. While this coupling may be explicitly understood for the walker system,^{17} exploring generalized dependencies of the phase parameters on the local wave field and other system parameters may allow us to rationalize other exotic pilot-wave systems where the physics at the lengthscale of the particle and over the time scale of its vibration has yet to be fully resolved.

## ACKNOWLEDGMENTS

J.W.M.B. gratefully acknowledges the financial support of the National Science Foundation (NSF) through Grant Nos. DMS-1614043 and CMMI-1727565. M.M.P.C. gratefully acknowledges the financial support of the Natural Sciences and Engineering Research Council (NSERC) through Grant No. 502891. We thank Professor R. R. Rosales for valuable discussions.

### APPENDIX: STROBOSCOPIC MODEL WITH SPATIOTEMPORAL DAMPING

We proceed by outlining a derivation of the stroboscopic model with spatial damping from the quasi-potential weakly viscous wave model (Eqs. (11)–(14)) presented by Milewski *et al.*^{24} Since we are interested in solving the problem in free space, we may introduce a Fourier transform so that the velocity potential $\varphi $ and wave height $h$ are written as

where $k=|k|$. This form ensures that $\varphi $ satisfies the Laplace equation and decays as $z\u2192\u2212\u221e$. We can then write the system (11)–(14) as two equations in spectral space for $\Phi $ and $a$:

To evaluate the pressure term, we note that the pressure exerted by the drop on the bath is zero outside of the contact area. We can further simplify the calculation by supposing that the pressure distribution inside the contact area does not depend on position so that

We can then evaluate the integral in Eq. (A3) exactly to yield

By non-dimensionalizing lengths by the droplet radius $x\u2032=x/R$ and time by the droplet’s natural oscillation frequency $t\u2032=\omega Dt$, it is then possible to rewrite the system (11)–(14) as a single forced, damped Mathieu equation for the non-dimensional $a(k,t)$. The non-dimensional equation contains four parameters: an effective Ohnesorge number, $Ohe$, determining the viscous damping of the waves; the Bond number, $Bo$, determining the relative importance of gravity to surface tension; the dimensionless bath acceleration, $\Gamma $; and the *vibration number*, $\Omega $, determining the ratio of the bath oscillation frequency to the natural oscillation frequency of the drop (see Table I). Dropping the primes and taking the limit $w\u21920$, we can write this equation as

The left-hand side is the same as Eq. (2.33) of Milewski *et al.*^{24}

A solution for this equation can formally be written in terms of the Green’s function:

We will now restrict attention to $(2,1)$ bouncers so that $F(\tau )$ is $4\pi /\Omega $ periodic. Moláček and Bush^{11} showed that in the limit of small $Oh$ and $Bo$ and for $k\u2248kF$, $G(t,\tau )$ may be approximated by

where $\beta (k)=1Me+\beta 1(k\u2212kF)2$. Under these assumptions, we may then approximate each wavemode in the neighborhood of $k=kF$ by the following integral:

In the limit of $\Gamma \u2192\Gamma F$, we note that the integrand of (A10) depends on two separate time scales. There is a fast time scale associated with the bath oscillations and the forcing due to the droplet, which scales with the Faraday period $TF$, and a slow time scale associated with the viscous decay of the waves and the horizontal trajectory of the droplet, which both vary over the memory time, $Me$. Assuming the vertical dynamics remain constant, we may approximate (A10) for large memory by averaging over the fast time scale so that

By taking the inverse Fourier transform of (A11) and switching orders of integration, we formally recover an analytic approximation to the wavefield generated by a walking drop:

where

In the limit of $\beta 1s\u2192\u221e$, the asymptotic approximation

was alluded to by Moláček and Bush^{11} and leads to the wave form included in the stroboscopic model with spatiotemporal damping, as presented in Eq. (10). A derivation of a similar result has been presented in Appendix C of Tadrist *et al.*^{43}

Using the approximation for $I(s,x)$ in Eq. (A14) leads to waves whose derivatives are singular at the origin. The bouncer wave field obtained from (A14) was found to give good comparisons with experimental observations as reported by Damiano *et al.*,^{16} provided the $1s$ factor was neglected in the same manner as presented by Oza *et al.*^{18} Spatio-temporal damping of this form was used explicitly by Oza *et al.*^{13} to study orbiting pairs of droplets, but in addition to neglecting the $1s$ factor, it was also necessary to adjust the Gaussian factor to $e\u2212|x|24\beta 1(t+TF)$ to remove singularities from the gradient of the wavefield underneath the drop. Addressing these singularities in a rigorous manner will be the subject of future work.