Downloaded from https://www.cambridge.org/core. MIT Libraries, on 22 May 2019 at 19:44:15, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2019.293

J. Fluid Mech. (2019), vol. 871, pp. 212–243.
doi:10.1017/jfm.2019.293

c Cambridge University Press 2019212

Bouncing phase variations in pilot-wavehydrodynamics and the stability of droplet pairs
Miles M. P. Couchman1 , Sam E. Turton1 and John W. M. Bush1, †1 Department of Mathematics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA(Received 23 August 2018; revised 16 February 2019; accepted 5 April 2019)

We present the results of an integrated experimental and theoretical investigationof the vertical motion of millimetric droplets bouncing on a vibrating fluid bath.
We characterize experimentally the dependence of the phase of impact and contact
force between a drop and the bath on the drop’s size and the bath’s vibrational
acceleration. This characterization guides the development of a new theoretical modelfor the coupling between a drop’s vertical and horizontal motion. Our model allows
us to relax the assumption of constant impact phase made in models based on thetime-averaged trajectory equation of Moláček and Bush (J. Fluid Mech., vol. 727,
2013b, pp. 612–647) and obtain a robust horizontal trajectory equation for a bouncing
drop that accounts for modulations in the drop’s vertical dynamics as may arise when
it interacts with boundaries or other drops. We demonstrate that such modulations
have a critical influence on the stability and dynamics of interacting droplet pairs.
As the bath’s vibrational acceleration is increased progressively, initially stationary
pairs destabilize into a variety of dynamical states including rectilinear oscillations,circular orbits and side-by-side promenading motion. The theoretical predictionsof our variable-impact-phase model rationalize our observations and underscore thecritical importance of accounting for variability in the vertical motion when modelling
droplet–droplet interactions.
Key words: drops, Faraday waves

1. Introduction

The hydrodynamic pilot-wave system discovered by Couder et al. (2005b) extends
the phenomenological range of classical systems to include behaviour previously
thought to be exclusive to microscopic, quantum systems (Bush 2015b; Bush et al.
2018). The system is characterized by drops, levitated on a vibrating bath, moving
in resonance with waves generated by their bouncing at the Faraday frequency.
It represents a macroscopic realization of de Broglie’s double-solution pilot-wave
theory (de Broglie 1956), wherein quantum particles move in resonance with a field
generated by the particle vibrating at the Compton frequency (Bush 2015a). The
dynamics at the Compton scale, specifically the interaction between the particle andwave, was not resolved by de Broglie. Neither is the Compton-scale interaction† Email address for correspondence: bush@math.mit.edu

1.1

Chaos

Coalescence

1.0

1(2, 1)
2),1
2((2, 1)1

(2 2)
(2,

0.6

(4,
(4
4 3)

(4

(2, 1)2

0.5
0.4

Chaos

0.3
0.2

os

(1, 1)

R = 0.40 mm
R = 0.36 mm
R = 0.32 mm2)ØBounce

0.7

Walk

Cha

0.8

,2

)

0.9

213

(4,

Downloaded from https://www.cambridge.org/core. MIT Libraries, on 22 May 2019 at 19:44:15, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2019.293

Bouncing phase variations in pilot-wave hydrodynamics(1, 1)10(1, 1)

0.5

1.0

Chaos

2

1.5

2.0

2.5

3.0

3.5

4.0

©/g

F IGURE 1. A regime diagram describing the various bouncing modes of a single drop ofradius R bouncing on a bath vibrating at a frequency of f = 80 Hz with fluid density ρ =949 kg m−3 , surface tension σ =p20.6 × 10−3 N m−1 and kinematic viscosity ν = 20 cSt.
The vibration number, Ω = 2πf ρR3 /σ , represents the non-dimensional drop size andγ /g denotes the dimensionless driving acceleration of the bath, where g is the gravitational
acceleration. For our experiments, the Faraday threshold was γF ≈ 4.2g. Shaded regions
denote theoretical bouncing modes obtained using the model of Moláček & Bush (2013a)
and markers denote experimentally measured thresholds between them (Wind-Willassen
et al. 2013). The drop sizes considered in our study are marked for reference by whitelines.

between a charge and its own electromagnetic field adequately resolved by the
Lorentz–Abraham–Dirac equation, the solutions of which break down on the Compton
time scale (Hammond 2010). An attractive feature of the walker system is that thefast dynamics responsible for both wave generation and particle propulsion may
be resolved both experimentally and theoretically. We thus proceed by resolving
experimentally and modelling theoretically the fast, vertical dynamics of bouncing
droplets with hopes of providing insights into the walking-droplet system in particularand pilot-wave systems in general.
The surface of a fluid bath vibrating vertically with an acceleration γ sin(ωt)
remains flat until the Faraday threshold, γF , above which it destabilizes into apattern of standing, subharmonic Faraday waves with period TF = 4π/ω (Benjamin &Ursell 1954; Miles & Henderson 1990). Millimetric droplets may bounce indefinitelyon the surface of the bath for driving accelerations above the bouncing threshold,γB < γF , exciting spatially extended temporally decaying waves at each impact (Walker
1978; Couder et al. 2005a; Damiano et al. 2016). A regime diagram characterizing
the dependence of the bouncing motion on a drop’s size and the bath’s driving
acceleration is shown in figure 1. The notation (i, j)k denotes that a drop’s verticalmotion has a period of i times that of the bath vibration, TF /2, and impacts thebath j times during this period; k = {1, 2} distinguishes between relatively small- andlarge-amplitude bouncers with the same periodicity (i, j) (Moláček & Bush 2013a). As
γ is increased progressively beyond γB , the bouncing amplitude of a drop increases
and the drop’s vertical motion eventually becomes synchronized with the Faraday
waves produced at each bounce, leading to resonant bouncing in the (2, 1) mode andDownloaded from https://www.cambridge.org/core. MIT Libraries, on 22 May 2019 at 19:44:15, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2019.293

214

M. M. P. Couchman, S. E. Turton and J. W. M. Busha substantial increase in the wave amplitude (Protière, Boudaoud & Couder 2006;
Eddi et al. 2011b). As γ is further increased beyond the walking threshold, γW , thedrop destabilizes into a ‘walker’, a dynamic state in which it moves steadily acrossthe surface of the bath, propelled by its own wavefield (Couder et al. 2005b; Protière
et al. 2006). The local wavefield, the slope of which prescribes the force acting on
the walking drop at each impact, depends not only on the drop’s present position
but on its past trajectory (Eddi et al. 2011b). The closer the driving acceleration isto the Faraday threshold, the longer the surface waves persist after each impact andthe more strongly the drop’s dynamics is affected by its past. This non-Markovian
feature of the walking-droplet system is generally referred to as ‘path memory’.
As recently reviewed by Turton, Couchman & Bush (2018), a hierarchy oftheoretical models of various degrees of sophistication has been developed to describe
walking droplets in a variety of settings (Protière et al. 2006; Eddi et al. 2011b;
Moláček & Bush 2013a,b; Oza, Rosales & Bush 2013; Bush, Oza & Moláček 2014;
Milewski et al. 2015; Dubertrand et al. 2016; Durey & Milewski 2017; Faria 2017;
Galeano-Rios, Milewski & Vanden-Broeck 2017; Nachbin, Milewski & Bush 2017).
The most sophisticated of such models (Milewski et al. 2015; Galeano-Rios et al.
2017) involve full numerical simulations of the coupled vertical and horizontal motionof a droplet and so are computationally expensive, with simulations taking 102 –104
times longer than real time. Moláček & Bush (2013b) developed coupled equations
describing a walking droplet’s vertical and horizontal motion. By time averaging theseequations over a bouncing period, they obtained a ‘stroboscopic’ horizontal trajectory
equation, simulation of which requires a computational time comparable to real time.
Oza et al. (2013) adapted the time-averaged model of Moláček & Bush (2013b)
into an integral model, by treating the walker as a continuous, rather than discrete,source of waves. In the stroboscopic model of Oza et al. (2013), all information
about the vertical dynamics is contained in a single fitting parameter −1 6 sin Φ 6 1,
as describes the phase shift between the periodic vertical motion of the bath and thedrop at impact.
While the time-averaged, stroboscopic models worked relatively well for capturing
the horizontal dynamics of single drops in various scenarios (Moláček & Bush2013b; Oza et al. 2013, 2014a,b; Labousse et al. 2016), significantly different valuesof sin Φ (in the range 0.16–0.5) were used in each of these studies. Further, owing
to the assumption of a constant impact phase, the stroboscopic approximation was
unable to capture accurately the observed stability of more complicated horizontal
trajectories such as the wobbling orbital motion of single drops in a rotating frame(Oza et al. 2014a; Labousse et al. 2016). Oza et al. (2017) and Arbelaiz, Oza & Bush(2018) demonstrated that for orbiting (Protière et al. 2006; Protière, Bohn & Couder2008) and promenading (Borghesi et al. 2014) droplet pairs, respectively, variationsin the vertical dynamics are apparent and significantly influence the pair’s horizontal
motion, features that cannot be captured with stroboscopic models. Galeano-Rios
et al. (2018) have likewise highlighted how subtle variations in a drop’s bouncing
phase can change both the speed and the direction of ratcheting droplet pairs (Eddiet al. 2008). Tadrist et al. (2018a) have demonstrated that such variations strongly
influence the outcome of droplet–droplet scattering events, and may be the source ofchaos in a number of hydrodynamic quantum analogues. Perrard & Labousse (2018)suggest that variations in the vertical dynamics may also serve as a mechanism forswitching between unstable periodic orbits in orbital pilot-wave dynamics.
In the current study, we develop a reduced model for the coupling between a drop’s
vertical and horizontal motion by using the linear spring model developed by Moláček

Downloaded from https://www.cambridge.org/core. MIT Libraries, on 22 May 2019 at 19:44:15, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2019.293

Bouncing phase variations in pilot-wave hydrodynamics

215& Bush (2013a) to characterize a drop’s interaction with the bath. This allows us
to extend the stroboscopic model of Oza et al. (2013) to account for variations in a
drop’s vertical dynamics, without having to resort to a full numerical simulation of the
problem. Our model is tested against the results of an experimental investigation ofidentical droplet pairs, and is found adequate in rationalizing the observed behaviourof the pair as the bath’s driving acceleration, γ , is increased progressively.
In § 2, we briefly review the horizontal trajectory equations of Moláček & Bush(2013b) and Oza et al. (2013) as they form the starting point for our theoreticaldevelopments. In § 3, we use high-speed imaging to characterize the dependence ofa single drop’s impact phase and contact force with the bath on the drop size andthe bath’s driving acceleration. Subsequently, we characterize experimentally how
neighbouring droplet pairs destabilize into horizontal motion as the bath’s driving
acceleration is increased progressively. In § 4, we develop a theoretical model forthe impact phase that is consistent with our experimental measurements, and develop
a horizontal trajectory equation that accounts for variations in the drop’s vertical
dynamics. Using this trajectory equation, we perform a linear stability analysis for
identical droplet pairs and show that our variable-phase model provides rationale
for the observed instabilities. Finally, in § 5, we discuss applications of our model
to systems containing multiple droplets, such as droplet rings and lattices, where
variations in the vertical dynamics are also expected to be significant.
2. Prior trajectory equations

We begin with a brief review of the model of Moláček & Bush (2013b) forthe horizontal motion of a bouncing droplet. Definitions of relevant variables andparameters are provided in table 1. The vertical position of the vibrating bath isdefined to be B (t) = −(γ /ω2 )sin(ωt). The wavefield generated by a single bounceof a drop at (x0 , t0 ) is described as a standing wave of the form (Moláček & Bush2013b)r ωt  J (k |x − x |)2νe mgTF
kF3 R20 F
0√H0 (x, x0 , t, t0 ) =e−(t−t0 )/(TF Me ) . (2.1)
S
cos
2 2
π σ 3kF R + Bo2
t − t0
The amplitude of the wave generated at each bounce depends on the drop’s verticalmotion through the phase parameter 0
Z t+TFωt0
dt0
FN (t )sin
2
S = t Z t+TF
.
(2.2)
0
0
FN (t ) dttHere, FN (t) is the contact force acting between the drop and the bath and is zero when
the drop is in free flight. One contribution of the current study is the first experimental
measurements of FN (t) for bouncing and walking droplets.
The phase parameter S encodes how efficiently the drop generates waves at each
impact. In order for the drop to bounce indefinitely, the average force imparted by
the bath on the drop must balance the drop’s
R t+Tweight. For a drop in a (2,1) bouncing
mode, this balance may be expressed as t F FN (t0 ) dt0 = mgTF (Moláček & Bush2013b). Evidently, it is the phase of the contact force relative to the bath’s oscillationthat governs the amplitude of the wave generated at each bounce. Moláček & Bush

Downloaded from https://www.cambridge.org/core. MIT Libraries, on 22 May 2019 at 19:44:15, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2019.293

216

M. M. P. Couchman, S. E. Turton and J. W. M. Bush

SymbolDefinition

ρ, σ , µ, ν = µ/ρ, νeµagf , ω = 2πf
TF = 2/f , λF , kF = 2π/λF
γ , γB , γW , γF
Γ = γ /gp
R, ωD = σ /ρR3 , Ω = ω/ωD
m = 4ρπR3 /3
Bo = ρgR2 /σ
xp = (xp , yp ), zpdt
H, h = H/cos(ωt/2)
hp = h(xp , t)
B
A=B+H
Z =z−AFN (t)
Λ1 = 0.48, Λ2 = 0.59
Td = 1/(νe kF2 )Td
Me =TF (1 − γ /γF )r
νe T F
mgkF3 R2A=
2π σ (3kF2 R2 + Bo)r
ρR
C = 0.17, D = Cmg+ 6πRµaσm
TF Me DmgkF2 TF Me R
β=D
S , C , sin(Φ)/2 = SC2ωρνe kFα=, = 222νe (1 + 2 )3kF σ + ρgr
2αξ=
kF TF Meκ=

Fluid density, surface tension, dynamic viscosity,
kinematic viscosity, effective kinematic viscosity(Moláček & Bush 2013b)
Air viscosity
Gravitational acceleration
Driving frequency of bath, angular frequency
Faraday period, wavelength, wavenumber
Peak driving acceleration of bath, bouncing threshold,walking threshold, Faraday threshold
Non-dimensional peak driving acceleration of bath
Drop radius, characteristic drop oscillation frequency,
vibration number
Droplet mass
Bond number
Horizontal, vertical position of drop
Inter-drop distance of a droplet pairTime
Wave height, strobed wave height
Local wave amplitude
Vertical position of bath
Vertical position of fluid surface
Vertical position of drop in reference frame where fluid
surface beneath drop is stationary
Contact force acting between drop and bath
Damping, spring coefficient for linear spring model(Moláček & Bush 2013a)
Wave decay time (Moláček & Bush 2013b)
Memory parameter
Wave-amplitude coefficientContact drag coefficient (Moláček & Bush 2013b),
averaged horizontal drag coefficient
Non-dimensional mass
Non-dimensional wave-force coefficient
Impact-phase parameters
Spatio-temporal damping coefficient (Turton et al. 2018)
Non-dimensional spatial damping coefficient

TABLE 1. Definitions of relevant variables and parameters.

Downloaded from https://www.cambridge.org/core. MIT Libraries, on 22 May 2019 at 19:44:15, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2019.293

Bouncing phase variations in pilot-wave hydrodynamics

217(2013b) assumed that the waves generated at each impact are dominated by modes
with a wavenumber close to the Faraday wavenumber, kF . Consequently, the phaseparameter S , the first coefficient in a Fourier series expansion of FN (t), captures thedominant response of the bath to the impacting droplet. A recent review of the wave
model of Moláček & Bush (2013b) and a further discussion of S is provided in Turton
et al. (2018).
We note that the wavefield described by (2.1) has been found to have insufficientspatial decay for |x − x0 | & λF , prompting a correction to the wave model inthe form of a spatio-temporal damping factor (Damiano et al. 2016; Oza et al.
2017; Turton et al. 2018). When modelling the dynamics of single drops, we
neglect this correction as the drop is primarily influenced by waves produced
within a distance λF of the drop’s current position. For instance, consider a drop
walking at a speed u = 10 mm s−1 at γ /γF = 0.85. Once the drop has moved a
distance λF from a previous impact, the wave produced at that impact will havedecayed to exp[−λF /(uTF Me )] ≈ 2 % of its original amplitude. Hence, the far-field
behaviour of the pilot wave has a negligible effect on the motion of a single walker
executing rectilinear motion. However, it is significant when modelling droplet–droplet
interactions, as is evident in our study of droplet pairs (see § 4.2).
Moláček & Bush (2013b) derived the following trajectory equation for thehorizontal motion of a droplet:!
rρRmẍp + C
FN (t) + 6πRµa ẋp = −FN (t)∇ H(xp , t),(2.3)
σ
where the net wavefield,

H(x, t) =X

Hn (x, xn , t, tn ),

(2.4)

n

with Hn described by (2.1), is assumed to be the linear superposition of the waves
generated by the drop at each previous bounce. When in contact with the bath, thedrop receives a horizontal impulse proportional to the gradient of the wavefield at
the point of impact and weighted by FN (t). The drop’s horizontal motion is resisted
by air drag and an additional drag force imparted by the bath during impact that is
proportional to FN (t). In order to solve equation (2.3), a model for a drop’s verticalmotion is required to obtain FN (t). In the inelastic ball model of Gilet, Vandewalle &Dorbolo (2009), FN (t) is not specified; however, the form of FN (t) could be deduced
from the linear spring model of Terwagne et al. (2013). Moláček & Bush (2013a)
have developed both linear and logarithmic spring models for the drop’s verticalmotion from which FN (t) can be inferred directly. It is computationally expensive
to simulate simultaneously the horizontal and vertical dynamics of a drop (Milewski
et al. 2015; Galeano-Rios et al. 2017); thus, the majority of subsequent models (Ozaet al. 2013; Bush et al. 2014; Dubertrand et al. 2016; Durey & Milewski 2017;
Faria 2017; Nachbin et al. 2017) have not accounted for variability in the vertical
dynamics.
To make theoretical headway, Moláček & Bush (2013b) focused on drops in
the resonant (2, 1) bouncing mode, where the time scale of a drop’s verticalmotion, TF = 0.025 s, is much smaller than that of the drop’s horizontal motion,
λF /|ẋp | ≈ 0.5 s. Time-averaging equation (2.3) over a (2, 1) bouncing period, TF ,yields a ‘stroboscopic’ horizontal trajectory equation (Moláček & Bush 2013b),mẍp + Dẋp = −mgC ∇h(xp , t),(2.5)

Downloaded from https://www.cambridge.org/core. MIT Libraries, on 22 May 2019 at 19:44:15, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2019.293

218

M. M. P. Couchman, S. E. Turton and J. W. M. Bush

where h = H/cos(ωt/2) represents the wavefield strobed at the drop’s bouncing period.
As the wavefield, H, on the bath surface (2.1) oscillates in time as cos(ωt/2), thegradient of the wavefield experienced by the droplet will depend on when in the
wavefield’s cycle the drop impacts. Since the impact is not instantaneous, the average
horizontal impulse exerted by the wave on the drop is modelled as −mgC ∇h(xp , t),where the phase parameter 0
Z t+TFωt
dt0
FN (t0 )cos
2
(2.6)
C= t
Z t+TF00
FN (t ) dtt

represents the average phase of impact with respect to the oscillating wavefield,
weighted by the normal force FN (t). As a first-order approximation, Moláček & Bush(2013b) assumed that the phase parameters S and C were constants and replaced the
product SC in (2.5) by the fitting parameter sin(Φ)/2, which then contains all theinformation about the bouncing phase.
Finally, for analytic simplicity, Oza et al. (2013) approximated the drop to be acontinuous rather than discrete source of waves, resulting in the trajectory equationmẍp + Dẋp = −mg sin Φ∇h(xp , t),
2

(2.7)

where
s
h(x, t) =mgkF3 R22νe
πTF σ (3kF2 R2 + Bo)

Z t

J0 (kF |x − xp (s)|)e−(t−s)/(TF Me ) ds.

(2.8)

−∞

We note that when Oza et al. (2013) approximated the sum in the wave modelof√Moláček & Bush (2013b)
(2.4)
by
an
integral(2.8),
thedecay
factor1/t
was√approximated by 1/ TF for simplicity, as it was assumed to be sub-dominant to the
exponential temporal decay. However, as shown in appendix A, this approximation
causes the amplitude of the wave built up after many bounces to be approximately
twice that predicted by Moláček & Bush (2013b). Thus, Moláček & Bush (2013b)
and Oza et al. (2013) were obliged to use different values of sin Φ, 0.5 and 0.3
respectively, in order to fit the same experimental data detailing the dependence ofa drop’s walking speed on γ . In this study, we wish to adopt the continuous modelof Oza et al. (2013), as it is more analytically tractable than the discrete modelof Moláček & Bush (2013b). However, as we wish to eliminate sin Φ as a fitting
parameter, we require that the wave amplitude of our model be consistent with thatof Moláček & Bush (2013b), which has been found to be in good agreement withexperimentally measured wavefields (Damiano et al. 2016). Thus, we proceed by
dividing the wave amplitude of Oza et al. (2013) in (2.8) by a factor of 2, so that it
approximately matches that of Moláček & Bush (2013b).
Using the dimensionless variables x̄ = kF x, h̄ = h/R and t̄ = t/(TF Me ), the trajectory
equation (2.7) of Oza et al. (2013) becomesκ x̄¨ p + x̄˙ p = −β C ∇ h̄(x̄p , t̄),
where
AMe
h̄(x̄, t̄) =
R

Z t̄
−∞

S J0 (|x̄ − x̄p (s̄)|)e−(t̄−s̄) ds̄.

(2.9)

(2.10)

Downloaded from https://www.cambridge.org/core. MIT Libraries, on 22 May 2019 at 19:44:15, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2019.293

Bouncing phase variations in pilot-wave hydrodynamics
(a)

(b)

Overhead camera
Enclosure

219

(c)

High-speed
cameraLight
Fluid bathbzT
Vertical acceleration©sin (ø t)

1 mm

F IGURE 2. (a) A schematic diagram of the experimental set-up. A detailed descriptionof the shaker used to vibrate the bath is given in Harris & Bush (2015). A high-speed
camera and an overhead camera are used to record the vertical and horizontal motion of adrop, respectively. A Plexiglas enclosure surrounds the bath to eliminate the influence of
air currents on the motion of the droplet. (b) A snapshot of a bouncing droplet captured
with the high-speed camera. The dark region at the top of the image is the far edge of thebath and appears out of focus due to the small depth of field of the magnifying lens. (c)Thresholding the image in (b) allows us to track the vertical positions of both the drop’s
top, zT , and the bath’s edge, B, and so infer the phase difference between the bouncing
motion and the bath vibration.

Note that the wave-amplitude coefficient, A (see table 1), has been corrected to beconsistent with that of Moláček & Bush (2013b). We proceed by determining thedependence of the phase parameters S and C on the driving acceleration of the bath,γ , the local wave amplitude at the location of a drop’s impact, hp = h(xp , t) and thedrop radius, R, both experimentally and theoretically. This will allow the horizontal
trajectory equation (2.9) to capture the subtle coupling between the drop’s vertical andhorizontal motion.
3. Experiments

Experiments were performed using the set-up described in Harris & Bush (2015)(see figure 2a). A bath of silicon oil (density ρ = 949 kg m−3 , surface tension
σ = 20.6 × 10−3 N m−1 , kinematic viscosity ν = 20 cSt) was sinusoidally shaken
in the vertical direction with a frequency of f = 80 Hz. Droplets of a desired size,
composed of the same silicon oil, were created using a piezoelectric droplet generator(Harris, Liu & Bush 2015). In § 3.1, we characterize the dependence of the phaseparameters S and C on the drop radius, R, and the driving acceleration of the bath,γ , for single bouncers and walkers. In § 3.2, we consider bound pairs of identical
drops and characterize both their vertical motion and destabilization into horizontal
motion as γ is increased progressively.
3.1. Single drops
We measure the vertical dynamics of drops of three different radii, R = 0.32 mm,
R = 0.36 mm and R = 0.40 mm, across the range of driving accelerations fromγ = 2.5g to γ = 4.2g, in increments of 0.1g. This parameter space was chosen tofocus on drops in the (2, 1) bouncing and walking regimes (see figure 1) in which

Bath height, b(a)(2, 2)

Borderline (2, 1)1(2, 1)1(2, 1)2(4, 2)

R = 0.40 mm, ©/©F = 0.58

R = 0.32 mm, ©/©F = 0.58

R = 0.32 mm, ©/©F = 0.67

R = 0.32 mm, ©/©F = 0.81

R = 0.32 mm, ©/©F = 0.95

0 0.5 1.0 1.5 2.0 0

0.5 1.0 1.5 2.0 0

0.5 1.0 1.5 2.0 0

0.5 1.0 1.5 2.0 0

0.5 1.0 1.5 2.00

0.5 1.0 1.5 2.0 0

0.5 1.0 1.5 2.0 0

0.5 1.0 1.5 2.0 0

0.5 1.0 1.5 2.00-g
0.5 1.0 1.5 2.0 0

0.5 1.0 1.5 2.0

Drop height, zT(b)

0.5 1.0 1.5 2.0 0(c)..

Vertical accel., z T

Downloaded from https://www.cambridge.org/core. MIT Libraries, on 22 May 2019 at 19:44:15, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2019.293

M. M. P. Couchman, S. E. Turton and J. W. M. Bush

220

tI
0
-g

0

0
-g
0.5 1.0 1.5 2.0 0

t/TF

tR

0
0
-g
-g
0.5 1.0 1.5 2.0 0 0.5 1.0 1.5 2.0 0

t/TF

t/TF

t/TF

t/TF

F IGURE 3. Examples of the time evolution of (a) the bath’s vertical position, (b) thedrop’s vertical position and (c) the drop’s vertical acceleration for the five bouncing modes
observed in our experiments. The Faraday period is TF = 1/40 s. The red data pointsin (a) and (b) are obtained by tracking the vertical position of the bath’s edge and the
top of the drop respectively. For the sake of comparison, pure sinusoids with the bath’s
oscillation frequency of 80 Hz are plotted with solid lines in (a). The solid lines in (b)are the result of applying the smoothing Savitzky–Golay filter to the experimental data,
which are then differentiated twice to obtain the drop’s vertical acceleration reported in(c). The red dashed lines in (c) indicate the free-fall acceleration, −g. Each green shaded
region indicates a calculated period of contact between the drop and the bath, starting at
the time of drop impact, tI , and ending at the time of drop release, tR . The vertical scalesin (a) and (b) are arbitrary.

the resonance between drop and wave assumed in the derivations of the trajectory
equations of Moláček & Bush (2013b) and Oza et al. (2013) is achieved.
In order to compute S and C for a given R and γ , we used the following procedure.
A high-speed camera (Phantom v5.2, macro lens, 2× teleconverter, tube extension)
was used to record simultaneously the motion of the drop and the bath at a frame
rate of 3000 fps, corresponding to 75 frames per Faraday period, TF . A snapshot fromone of these videos is shown in figure 2(b). The dark stripe at the top of the frame
corresponds to the far edge of the bath and is out of focus due to the small depthof field of the magnifying lens. As shown in figure 2(c), this same frame can bethresholded to yield an outline of the drop’s edge, the top of which we denote by zT ,and a horizontal line near the top of the frame, representing the vertical position of thebath’s edge, B . By thresholding each frame in the video, we obtain time series B (t)and zT (t), examples of which are shown in figures 3(a,b). If a droplet was walking,it was recorded as it moved through the high-speed camera’s field of view while theoverhead camera simultaneously recorded its horizontal motion.

Downloaded from https://www.cambridge.org/core. MIT Libraries, on 22 May 2019 at 19:44:15, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2019.293

Bouncing phase variations in pilot-wave hydrodynamics

221

The contact force, FN (t), can be computed directly from zT (t). A Savitzky–Golay
filter is used to compute the drop’s vertical acceleration, z̈T (figure 3c). The
acceleration of the top of the drop provides a good approximation to the accelerationof the drop’s centre of mass provided the drop deformation during impact is small.
In our experiments, the maximum deformation observed is ≈ 5 % of the undeformed
drop diameter. The regions of constant acceleration in the time series of z̈T occur
when the drop is in free flight, with acceleration −g, while the regions of increased
acceleration, shaded in green in figure 3(c), occur when the drop is in contact withthe bath; FN (t) is then obtained by subtracting the free-fall acceleration, −g, from z̈Tand is necessarily zero when the drop is not in contact with the bath. We note thatwhen a drop is walking in a straight line at a constant speed, the vertical position
recorded by the high-speed camera will be of the form zT (t) + at + b, where zT (t)would be the vertical position recorded if the drop was bouncing in place and thelinear function reflects the drop’s translation across the camera’s field of view. Since
the linear function at + b makes no contribution to the second derivative, the drop’s
horizontal motion does not affect the computed vertical acceleration. When a drop isnear the (2, 2) to (2, 1)1 transition, it is difficult to characterize the bouncing mode
(Galeano-Rios et al. 2017). We classify a drop as being in a ‘borderline’ (2, 1)1
mode when two distinct impacts per Faraday period are discernible by analysing
FN (t) in figure 3(c), but the droplet does not rebound upwards between these impacts,
as determined by analysing zT (t) in figure 3(b).
Our direct measurements of the vertical motion of both bouncing and walking
droplets will be valuable in benchmarking and guiding the development of theoreticalmodels. Our measurements of the contact times and contact forces across a rangeof bouncing modes (figure 3) are very similar to those predicted theoretically by
the droplet impact models of Milewski et al. (2015) and Galeano-Rios et al. (2017).
For example, comparing figure 3c to figures 8 and 11 in Galeano-Rios et al. (2017)indicates that the model of Galeano-Rios et al. (2017) accurately captures the (2, 2)
to (2, 1)1 to (2, 1)2 progression in bouncing modes as γ is increased. For a drop inthe (2, 2) mode, as γ is increased progressively, the contact force during the first
impact decreases relative to that during the second impact, and the time of flightbetween these two impacts also decreases. At the critical γ value at which these two
initially separate impacts merge, the drop becomes a (2, 1)1 bouncer. As γ is further
increased, the time of the drop’s release from the bath, tR , remains approximately
constant, while the time of impact, tI , moves progressively closer to tR until the drop
eventually enters a (2, 1)2 mode.
Once FN (t) is computed, the trapezoidal rule is used to numerically calculate S(2.2) and C (2.6). When performing the integrals, the time t = 0 is defined such that
B (t) = −(γ /ω2 )sin(ωt). For each R and γ , approximately 20 bounces were recorded
and the phase parameters were computed for each bounce and then averaged. The
standard deviation of the phase parameters before averaging is of the order of 1 %.
To compute S and C for drops in a (4, 2) mode, the integrals in (2.2) and (2.6)
were evaluated from t to t + 2TF . We emphasize that S and C represent average
impact phases weighted by the contact force FN (t). We have observed that within the(2, 1)2 mode, the times of drop impact, tI , and release, tR , from the bath may remain
approximately constant as γ is increased. Thus, one might infer that the drop’s impact
phase is independent of γ . However, even when tI and tR remain constant, the profile
FN (t) for tI 6 t 6 tR shifts as γ increases, resulting in changes to both S and C . Hence,
knowledge of the time dependence of the contact force is critical in order to model
effectively the coupling between the drop’s horizontal and vertical motion.

(a) 1.0(b) 1.0

0.8

0.8

s

0.6

c

(c)

0.6

0.4

0.4

0.2

0.2

0

u (mm s-1)

Downloaded from https://www.cambridge.org/core. MIT Libraries, on 22 May 2019 at 19:44:15, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2019.293

M. M. P. Couchman, S. E. Turton and J. W. M. Bush

222

0.6

0.7

0.8

©/©F

0.9

0

1.0

0.6

0.7

0.8

©/©F

0.9

1.0

20
15
10
5
0
0.70

0.75

0.80

0.85

0.90

0.95

1.00

©/©F

F IGURE 4. The dependence of (a) the phase parameter S , (b) the phase parameter C
and (c) the walking speed, u, on the non-dimensional driving acceleration of the bath,
γ /γF , for single bouncing drops of size R = 0.32 mm (red), R = 0.36 mm (green) and
R = 0.40 mm (blue). The unfilled and filled markers indicate experimental data anddenote whether the droplet is bouncing or walking, respectively. The shape of the marker
denotes the experimentally observed bouncing mode: (2, 2) mode (∗), borderline (2, 1)1
mode (?), (2, 1)1 mode (@), (2, 1)2 mode (A), (4, 2) mode (E). Characteristic error bars
are shown. In (a) and (b) the dotted and solid curves are the corresponding theoreticalpredictions (see § 4.1) and denote whether the droplet is predicted to bounce or walk,
respectively. For values of γ where (∗), (?) or (@) bouncing modes were experimentallyobserved, the (2, 1)1 phase functions of the form given in table 2 are used for thetheoretical predictions. For (A) and (E) modes, the (2, 1)2 phase functions are used.
The discontinuities in the theoretical curves occur at the (2, 1)1 to (2, 1)2 transitions.
In (c), the solid and dashed lines represent, respectively, the walking speeds predictedby our variable-impact-phase trajectory equation (see § 4.1) and the constant-impact-phasetrajectory equation of Oza et al. (2013).

The dependences of S , C and the walking speed, u, on R and γ are shown infigure 4. Both S and C are seen to have two distinct trends separated by the (2, 1)1 to(2, 1)2 transition. Before this transition, both S and C increase approximately linearlywith γ , while after the transition S saturates to a constant value of approximately 1
while C decreases monotonically.
3.2. Pairs of identical drops
We proceed by highlighting the critical role of bouncing phase variations on the
bouncing to walking transition of identical droplet pairs. As shown by Eddi et al.
(2009), droplet pairs can exist in a discrete set of bound states corresponding to223(2, 1) mode(n = 1)(3/2) (2) (5/2) (3) (7/2) (4) (9/2) (5) (11/2) (6)

1.0

©W

0.8
(2, 2) → (2, 1)©/©F

Downloaded from https://www.cambridge.org/core. MIT Libraries, on 22 May 2019 at 19:44:15, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2019.293

Bouncing phase variations in pilot-wave hydrodynamics

0.6(1, 1) → (2, 2)

0.4
0.2

©B
(n = 1) (2)

0
0

1

(3)

(4)
2

(5)

(6)
(7)
(1, 1) mode
3
4

(8)

(9)
5

(10)

(11)
6

7

d/¬F

F IGURE 5. The dependence of the inter-drop separation distance, d, on the dimensionlessdriving acceleration, γ /γF , for stable, bound droplet pairs of radius R = 0.32 mm. Circles
indicate data obtained by increasing the driving acceleration and triangles indicate dataobtained by lowering it again. The top-most data points indicate the instability thresholds,γ∗ , at which the stationary pairs go unstable to horizontal motion, in a manner detailed
in figure 7. A characteristic error bar is shown for an n = 1 data point. The bouncing
threshold γB , (1, 1) → (2, 2) transition, (2, 2) → (2, 1) transition, and walking threshold,γW , for a single drop of radius R = 0.32 mm are shown for reference. The vertical ticks
along the horizontal lines denoting γB and γW indicate the stable separation distances, dn ,predicted by our theoretical model (see § 4.3).

each drop bouncing in a minimum of its neighbour’s wavefield. The quantized set of
possible inter-drop separation distances, dn , is characterized by the binding numbern = {1, 3/2, 2, 5/2, . . . , }, where integers and half-integers denote pairs bouncing
in- and out-of-phase, respectively (Eddi et al. 2009). We consider the behaviour ofbouncing pairs of two different radii, R = 0.32 mm and R = 0.40 mm, with various
binding numbers n, as γ is increased progressively. These sizes were chosen to ensure
that, close to their instability thresholds, the pairs had different vertical dynamics,
specifically, the larger and smaller pairs are in (2, 1)1 and (2, 1)2 bouncing modes,
respectively.
For each pair size, γ is set slightly above the bouncing threshold for a single drop,γB , and two identical drops are created that bounce in a (1, 1) mode. As (1, 1) pairs
bounce with the same period as the vibrating bath they can only bounce in phase. It
was possible to initialize (1, 1) pairs in six binding numbers n ∈ [1, 6] (see figure 5).
Our findings extend the work of Eddi et al. (2009) who reported three possible n ∈[1, 3], for (1, 1) pairs. Pairs starting with a different, intermediate d, spontaneously
adjust their spacing to reach the closest stable dn . For each (1, 1) pair in n ∈ [1, 6],
we increased γ incrementally, causing the inter-drop distances to evolve as reported infigure 5, as the pairs transitioned from (1, 1) to (2, 1) bouncing modes. We note thatas γ is increased, the bouncing phase can flip spontaneously. For example, in-phasen = 2 and n = 4 (1, 1) pairs end up as out-of-phase n = 3/2 and n = 5/2 (2, 1) pairs,
respectively. The (2, 1) pairs ultimately go unstable to horizontal motion at γ = γ∗ ,indicated by the top-most data points in figure 5. For pairs of size R = 0.32 mm, theinstability threshold, γ∗ , occurs above the walking threshold of a single drop, γW , asshown in detail in figure 7.

Downloaded from https://www.cambridge.org/core. MIT Libraries, on 22 May 2019 at 19:44:15, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2019.293

224

M. M. P. Couchman, S. E. Turton and J. W. M. Bush

Due to the increased amplitude and spatial extent of a drop’s wavefield close to γ∗ ,(2, 1) pairs could be initialized with larger binding numbers n that were not accessible
in the vicinity of γB . As shown in figure 5, eleven possible n were found for (2, 1)
pairs, as compared to the six reported by Eddi et al. (2009). For each new n, a (2, 1)
pair was initialized at a driving acceleration slightly below the anticipated instability
threshold, and γ was then increased until γ∗ was found. For each of the eleven (2, 1)
binding numbers n, we also progressively decreased γ until the pairs coalesced withthe bath. This progression revealed hysteresis in the dependence of d on γ as well as
the existence of additional (1, 1) binding numbers up to n = 11. These additional states
were impossible to achieve by initializing the pairs in the vicinity of γB due to theweakness of the droplet–droplet interaction at such large n and low γ . For drops ofradius R = 0.40 mm, the dependence of d on γ and n is very similar to that presented
in figure 5, except that the instability thresholds, γ∗ , occur below the walking thresholdof an individual drop, γW , as shown in figure 7.
In addition to measuring the inter-drop distance d, we used the technique presented
in § 3.1 to measure the phase parameters S and C for drops in period-doubled (2, 1)
pairs with binding numbers n = 1 and n = 2. The results are presented in figure 6.
For drops in the (2, 1)1 mode, both S and C increase with decreasing n. For drops inthe (2, 1)2 mode, S has saturated to 1 and does not depend strongly on n, while C
now decreases with decreasing n. Notably, the bouncing mode transitions are different
when the drops are in pairs. Specifically, the (2, 2) → (2, 1)1 , (2, 1)1 → (2, 1)2 , and(2, 1)2 → (4, 2) transitions occur at lower values of γ than for single drops, and thesevalues decrease monotonically with decreasing n.
In figure 7, we show the dependence of a pair’s instability threshold, γ∗ , on thebinding number, n, for the two drop sizes considered. The larger pair, of radius
R = 0.40 mm, bounces in a (2, 1)1 mode in the vicinity of γ∗ and destabilizes below
the walking threshold of an individual drop, γW . For all n, the pair destabilizes into
oscillations along the line connecting the drops, which we refer to henceforth asco-linear oscillations. An example of the onset of co-linear oscillations is shownin figure 8(a). Conversely, the smaller pair, of radius R = 0.32 mm, bounces in a(2, 1)2 mode in the vicinity of γ∗ and goes unstable above γW . For the smaller pair,the type of instability now depends on n. An n = 1 pair goes unstable to orbital
motion (figure 8b) while pairs in higher n > 1, go unstable to promenading states,in which the drops walk side-by-side (Borghesi et al. 2014; Arbelaiz et al. 2018).
Promenading pairs with n = {3/2, 2, 5/2} move in straight lines (figure 8c), whilepairs with n > 3 exhibit lateral oscillations (figure 8d). We note that for pairs of size
R = 0.32 mm, the type of instability is sensitive to small size variations between
the two drops. Specifically, if the size variation is &5 %, the difference in walking
thresholds between the two drops causes one drop to orbit the other, which initiallyremains stationary. Orbital motion then occurs at most separation distances, not justfor the n = 1 pair. As seen in figure 7, for both pair sizes, γ∗ tends to γW as nincreases, due to the interaction between the drops decreasing with increasing n.
4. Theoretical models for the phase parameters S and C

In this section, we deduce theoretical relations for the dependence of the phaseparameters S and C on the bath’s driving acceleration, γ , the local wave amplitude,
hp = h(xp , t), and the drop radius, R, for drops in a (2, 1) bouncing mode. We begin byconsidering a stationary (2, 1) bouncer. Using the non-dimensional variables z̄ = z/Rand τ = ωD t, the vertical position of the fluid surface beneath the drop, A(t), is theDownloaded from https://www.cambridge.org/core. MIT Libraries, on 22 May 2019 at 19:44:15, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2019.293

Bouncing phase variations in pilot-wave hydrodynamics(a)

R = 0.32 mm

R = 0.32 mm

(b)

1.0

225

0.6
0.8

s

c

0.6

0.4

0.2
0.4
0.6

(c)

0.7

0.8

0

0.9

R = 0.40 mm

(d)

1.0

0.6

0.7

0.8

0.9

R = 0.40 mm
0.6

0.8

s

0.4

0.6

c

0.4

0.2

0.2
0
0.55

0.60

0.65

0.70

©/©F

0.75

0.80

0
0.55

0.60

0.65

0.70

0.75

0.80

©/©F

F IGURE 6. The dependence of the phase parameters S and C on the non-dimensional
driving acceleration, γ /γF , for bouncing pairs of size R = 0.32 mm and R = 0.40 mm
with binding numbers n = 1 (red) and n = 2 (green). The corresponding data for a single
bouncer, from figure 4(a,b), are shown for reference, in black. The shape of each data
point indicates the bouncing mode: (2, 2) mode (∗), borderline (2, 1)1 mode (?), (2, 1)1
mode (@), (2, 1)2 mode (A), (4, 2) mode (E). Characteristic error bars are shown. Data
were collected up until the instability threshold, γ∗ , at which the pair went unstable tohorizontal motion. Theoretical predictions of our variable-phase model (see § 4.3) are
shown for reference by solid curves.

superposition of the harmonic, vertical shaking of the bath, B (t), and the sub-harmonicstanding wavefield on the surface of the bath, H(t), generated by the (2, 1) bouncer:γΩτĀ = B̄ + H̄ = − 2 sin(Ωτ ) + h̄p cos
.
(4.1)
Rω
2
The vibration number, Ω = ω/ωD , is the ratio pof the bath’s driving frequency to thedrop’s internal frequency of oscillation, ωD = σ /ρR3 . We define Z = z − A to be
the height of the base of the drop above the fluid surface (Moláček & Bush 2013b).
To model the vertical motion of the droplet we use the linear spring model ofMoláček & Bush (2013a),
!

!
d2 Z̄
dZ̄
γ
h̄p Rω2Ωτ+ H(−Z̄ ) Λ1
+ Λ2 Z̄ = −Bo 1 + sin(Ωτ ) −
cos
,
dτ 2
dτ
g
4g
2
(4.2)

1.10
Orbit

1.05

©*/©W

Downloaded from https://www.cambridge.org/core. MIT Libraries, on 22 May 2019 at 19:44:15, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2019.293

M. M. P. Couchman, S. E. Turton and J. W. M. Bush

226

Promenade (no osc.)

Promenade (with osc.)

R = 0.32 mm
1.00
Co-linear oscillations

0.95
0.90

R = 0.40 mm
0

1

2

3

4

5

6

d/¬F

F IGURE 7. The instability thresholds, γ∗ , and type of instability for identical droplet
pairs of radius R = 0.40 mm (blue) and R = 0.32 mm (red), as are in (2, 1)1 and(2, 1)2 bouncing modes, respectively. Circles denote experimental data and the squares
connected by dashed lines are the theoretical predictions of the variable-phase model
presented in § 4.3. Characteristic error bars are shown. Note that the experimental data
for R = 0.32 mm correspond to the top most data points in figure 5. The instability
thresholds are normalized by the walking threshold of a single drop, γW , as was measured
to be γW ≈ 0.77γF and γW ≈ 0.88γF for drops of size R = 0.40 mm and R = 0.32 mm,
respectively. The type of instability predicted by the theoretical model is discussed in § 4.3.

where H(x) denotes the Heaviside function and the right-hand side of (4.2) is anexpression for the effective gravity in the frame of reference where the fluid surface
beneath the droplet is stationary. When in contact with the bath, the drop feels botha drag force and a spring force whose magnitudes, prescribed by the parameters Λ1
and Λ2 , respectively, depend on the Weber number characterizing the droplet’s impact(Moláček & Bush 2013a). For simplicity, we fix these parameters to be Λ1 = 0.48
and Λ2 = 0.59, which match the (2, 1)1 to (2, 1)2 transition predicted by the linear
spring model to that observed in our experiments. We note that these values of Λ1and Λ2 are slightly different than those used in the study of Moláček & Bush (2013a)(Λ1 = 0.41 and Λ2 = 0.60), owing to their approximating h̄p = 0 in (4.2), which one
expects to hold only at low memory.
It is worth pointing out that Moláček & Bush (2013a) also developed a nonlinear
logarithmic spring model in an attempt to correct discrepancies between experimentallymeasured bouncing mode transitions and those predicted by the linear spring model.
However, as these discrepancies arose only at high γ , they might also be attributable
to the influence of surface waves on the drop’s vertical dynamics, which was notconsidered by Moláček & Bush (2013a). Thus, it is not entirely clear that thelogarithmic spring model is more accurate than the linear spring model; therefore, inthis study, we adopt the linear spring model, which can be solved analytically.
For Z̄ > 0, when the drop is in free flight, the general solution to (4.2) is


γ
Ωτ
Bo 2
sin(Ωτ ) − h̄p cos+ A1 τ + A2 ,(4.3)
Z̄A (τ ) = − τ +2Rω22
where A1 and A2 are constants of integration. For Z̄ < 0, when the drop is in contact
with the bath, the general solution to (4.2) is227(b) 0.5(a) 3.5

y/¬F

3.4

0

3.2

4

3.1

d
0

3

25

50

75

100

-0.5

0

0.5

x/¬F

t (s)
(d) 10

8

8

6

6

y/¬F

(c) 104

21Drop speed (mm s-1)

d/¬F

3.3

y/¬F

Downloaded from https://www.cambridge.org/core. MIT Libraries, on 22 May 2019 at 19:44:15, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2019.293

Bouncing phase variations in pilot-wave hydrodynamics

0

4
2

2
0
-2

0

2

0
-4

-2

0

x/¬F

2

4

x/¬F

F IGURE 8. Examples of the onset of instability of bound droplet pairs: (a) co-linear
oscillations (R = 0.40 mm, n = 3.5), (b) orbiting (R = 0.32 mm, n = 1), (c) promenading(R = 0.32 mm, n = 2) and (d) promenading with lateral oscillations (R = 0.32 mm, n = 5).
Panel (a) shows the evolution of the inter-drop distance d. Panels (b–d) show trajectoriesof the droplets coloured according to the droplet speed.

Z̄B (τ ) = B1 e

−Λ1 τ /2pcos

4Λ2 − Λ21 τ
+ B2
2

!
−

Bo
Λ2


−

−

2
γ Bo 
 (Λ2 − Ω ) sin(Ωτ ) − Λ1 Ω cos(Ωτ )g 
Λ22 − 2Λ2 Ω 2 + Λ21 Ω 2 + Ω 4

h̄p Rω2
γ




Ωτ
Ωτ
+ 2Λ1 Ωsin
2
2
2 
,
2
2

16Λ2 − 8Λ2 Ω 2 + 4Λ1 Ω 2 + Ω 4! (4Λ − Ω 2 )cos



(4.4)

where B1 and B2 are constants of integration.
To identify (2, 1) bouncing solutions, we seek periodic solutions to (4.2) with period
τ = 4π/Ω where the droplet impacts the bath once during this period. Denoting thetimes at which the droplet leaves and hits the surface of the bath during a bouncing
period by τ1 and τ2 , respectively, we require that Z̄A (τ ) and Z̄B (τ ) satisfy theDownloaded from https://www.cambridge.org/core. MIT Libraries, on 22 May 2019 at 19:44:15, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2019.293

M. M. P. Couchman, S. E. Turton and J. W. M. Bush

228
Ø†1

(a) 1.5Ø†1 + 4π

(b)

-

1.0

z-

Ø†2

Ø†1

1.5

Ø†2

Ø†1 + 4π

-

zA

1.0

zA

0.5

0.5

0

0
-

a

-0.5
0

2π

4π

-

-0.5

zB

6π

8π

a
0

-

zB
2π

Ø†

4π

6π

8π

Ø†

F IGURE 9. One period of the vertical displacement of a drop of radius R = 0.36 mm
predicted by the linear spring model in (a) a (2, 1)1 bouncing mode at γ = 3g with local
wave amplitude h̄p = 0.03, and (b) a (2, 1)2 bouncing mode at γ = 4g with h̄p = 0.1, as
viewed in the laboratory frame. The trajectories of the droplet in free flight, zA , and incontact with the bath, zB , are shown by the red and blue curves, respectively. The black
curves denote the vertical position of the fluid surface, A, beneath the droplet.

following system of equations:Z̄A (τ1 ) = 0,

Z̄A (τ2 ) = 0,Z̄˙A (τ2 ) = Z̄˙B (τ2 ),

4π
Z̄B (τ2 ) = 0, Z̄B τ1 +Ω

4π
Z̄˙A (τ1 ) = Z̄˙B τ1 +
.
Ω






= 0,


(4.5)






These conditions ensure that the drop’s vertical position and speed are continuous as
it enters and leaves the bath, and that the bouncing motion has a period of τ = 4π/Ω.
For given values of γ , hp , and R, we can numerically solve the nonlinear system (4.5)
to obtain τ1 , τ2 , A1 , A2 , B1 and B2 . Figure 9 shows examples of the bouncing solutions
so deduced for droplets in (2, 1)1 and (2, 1)2 modes.
For a given drop radius R, we can then determine the regions in the (γ , hp )-planewhere (2, 1)1 and (2, 1)2 bouncing solutions exist, and compute the corresponding
values of the phase parameters S and C . We thus obtain surfaces for S and C as
functions of γ and hp . Next, we observe that each of these surfaces can be collapsed
onto a curve by expressing S and C as functions of a new variable, a suitable linear
combination of γ and hp . Examples of these collapses for R = 0.36 mm are shown infigure 10(a–d), where the linear combination that best collapses the data is indicatedon each x-axis. We observe that for S and C in the (2, 1)1 mode, the data are bothwell approximated by linear relations, the coefficients extracted from a least squares
fit. For S in the (2, 1)2 mode, we chose to fit the data with an exponential curvethat asymptotes to one at high memories, as this is the maximum value that S cantheoretically attain. Similarly, for C in the (2, 1)2 mode, we fit the data with an
exponential curve that asymptotes to zero at high memories, on the grounds that we
do not expect C to become negative, as would correspond to a reversal in walking
direction. The same analysis is then repeated for drops of sizes R = 0.32 mm and
R = 0.40 mm, and the dependencies of the function coefficients a, b and c on the nondimensional drop size Ω are presented in figure 10(e–h). We see that the dependenciesof the coefficients on Ω are well approximated by linear functions. Based on the data(b) 0.53(c) 1.00(d) 0.4s(2,1)1c(2,1)1s(2,1)2c(2,1)2229(a) 0.8

0.45--

a + b˝ + ch p

0.37
2.40

3.1
-

2.65

2.90

-2
-4
-6
0.6

c
0.7

0.8

0.9

c(2,1)1 coefficients

b
a

3.8

b

0

a
c

-2
-4
-6
0.6

Ø

0.7

0.8

Ø

0.9

3.0

3.5

4.0
-

˝ - 5.84hp

(g)
2

0
2.5

4.4

˝ - 5.73hp

(f)
2

0.2

-

˝ - 7.57hp(e)0-b(˝+chp-2)

0.88 1 - ae
2.6 3.2-

˝ - 4.23hp-

ae-b(˝+chp-2)-

a + b˝ + ch ps(2,1)2 coefficients

2.8

0.94

(h)
5

b

0

a

-5

c

-10
0.6

0.7

0.8

0.9

Ø

c(2,1)2 coefficients

0.6

0.4
2.5

s(2,1)1 coefficients

Downloaded from https://www.cambridge.org/core. MIT Libraries, on 22 May 2019 at 19:44:15, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2019.293

Bouncing phase variations in pilot-wave hydrodynamics

5 b
0
-5

a
c

-10
0.6

0.7

0.8

0.9

Ø

F IGURE 10. The dependence of the phase parameters S and C on Γ = γ /g and h̄p = hp /R
for a drop of radius R = 0.36 mm in (a,b) the (2, 1)1 mode and (c,d) the (2, 1)2 mode.
The red crosses are the values of the phase parameters obtained using the linear springmodel and the black curves indicate approximate functional forms. Panels (e–h) show thedependencies of the coefficients for the functional forms in (a–d) on the non-dimensional
drop radius Ω. Lines of best fit are shown for reference.
PhaseparameterS(2,1)1
C(2,1)1
S(2,1)2
C(2,1)2

Functional
form

Parameter values

a + bΓ + ch̄pa + bΓ + ch̄p
1 − ae−b(Γ +ch̄p −2)ae−b(Γ +ch̄p −2)

a = −3.71Ω + 1.35, b = 1.24Ω − 0.224, c = −13.6Ω + 6.83
a = −1.92Ω + 1.17, b = 0.490Ω − 0.108, c = −7.29Ω + 3.32
a = 1.79Ω, b = −5.60Ω + 7.65, c = −8.00Ω + 0.168
a = −3.55Ω + 4.60, b = −6.06Ω + 6.84, c = −8.57Ω + 0.453

TABLE 2. Approximate functional dependences of the phase parameters S and C on the
non-dimensional driving acceleration of the bath, Γ , the local wave amplitude, h̄p , andthe dimensionless drop radius, Ω, for drops in a (2, 1)1 or (2, 1)2 bouncing mode. These
relations are valid across the range of drop radii 0.6 6 Ω 6 0.9 considered in this study.

in figure 10, we are able to develop functional forms for S and C , which depend on
γ , hp and R, as presented in table 2.
The relations for S and C given in table 2 were developed by considering ahorizontally stationary bouncing droplet. They are also expected to be valid forwalking drops provided that changes in γ and hp occur slowly relative to thebouncing period, TF , in which case the droplet’s vertical dynamics will closely
resemble that of a stable bouncing state. Specifically, for fixed γ , we require thatthe time scale of bouncing, TF , be small relative to that over which the local wave
amplitude changes by a characteristic value 1h, Thp = 1h/[(∂hp /∂x)|ẋp |]. In a typicalexperiment, ∂hp /∂x ≈ 1h/λF and |ẋp | ≈ 10 mm s−1 which yields Thp /TF ≈ 20  1.
The strobed trajectory equation (2.9) coupled with the functions S = S (Γ , h̄p , Ω)and C = C (Γ , h̄p , Ω) presented in table 2, can now be used to compute the horizontal
trajectory of a (2, 1) walker while simultaneously accounting for variations in theDownloaded from https://www.cambridge.org/core. MIT Libraries, on 22 May 2019 at 19:44:15, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2019.293

230

M. M. P. Couchman, S. E. Turton and J. W. M. Bushdrop’s vertical dynamics. Except when modelling simple scenarios, such as single
bouncers or straight-line walkers in unbounded geometries, incorporation of theinfluence of impact phase variations requires that the trajectory equation be solved
numerically. Specifically, one must first discretize equation (2.9) using a numericaltime-stepping scheme. At each time step, the wave amplitude at the position of thedrop, hp , must be computed in order to obtain values for S and C , after which thedrop’s horizontal position can be updated via the trajectory equation (2.9).
4.1. Impact-phase variations for single drops
We proceed by calculating how S and C vary with γ for single bouncers and straightline walkers, in order to compare the predictions of our theoretical models for S andC with our experimental results in figure 4. According to the trajectory equation (2.9),
the following equations must be satisfied for a single drop walking in a straight linewith constant speed u:
AMe S (Γ , h̄p , Ω)
√
h̄p =
,
(4.6a)
R 1 + ū2rp1
AMe β S (Γ , h̄p , Ω)C (Γ , h̄p , Ω).
(4.6b,c)ū = φ − (1 + 1 + 4φ), φ =2RFor a given drop radius R, the walking threshold, γW , is obtained by finding the valueof γ that satisfies both φ = 2 and (4.6a) with ū = 0. To compute the phase parameters
for a bouncer, at a given value of γ , we first solve equation (4.6a) for h̄p , setting ū = 0.
Knowing h̄p then allows for S and C to be computed. In the case of a walker, we firstsolve the system of (4.6a) and (4.6b) in order to obtain ū and h̄p . From h̄p , we canthen compute S and C .
In figure 4, we compare our theoretical predictions for S and C with ourexperimental data. For values of γ where (2, 2) and (4, 2) modes were experimentallyobserved, we used our phase relations for (2, 1)1 and (2, 1)2 modes, respectively, as
an approximation. Our model is seen to capture all of the observed trends in S and C ,however, there is a horizontal shift between our experimental and theoretical data forthe phase parameter S in the (2, 1)1 mode. This shift is a result of slight differencesin the profile of the contact force, FN (t), observed experimentally and predicted
theoretically by the linear spring model. In particular, S is highly sensitive to therelative heights of the two maxima in FN (t) in the (2, 1)1 mode seen in figure 3,
while C is not.
In figure 4(c), we compare our experimental data of walking thresholds and speeds
for single drops with the theoretical predictions of both the constant-impact-phasemodel of Oza et al. (2013) and the variable-phase model developed in this study.
Oza et al. (2013) used a fixed value of sin(Φ)/2 = SC = 0.30, chosen to best matchthe dependence of the walking speed on γ for a drop of radius R = 0.40 mm. As
seen in figure 4(c), this assumption of constant-impact-phase results in significant
discrepancies between the theoretical predictions and experimental data for otherdrop sizes. By capturing the dependence of S and C on R and γ , our model isable to adequately predict the walking thresholds across a range of drop sizes,
R = 0.32, 0.36 and 0.40 mm. We note that for a drop radius of R = 0.36 mm,
the predicted walking speed at the walking threshold is non-zero, as the walking
threshold occurs at the locus of the (2, 1)1 to (2, 1)2 transition, where the phaseparameters change discontinuously. We also note that the model of Oza et al. (2013)

Downloaded from https://www.cambridge.org/core. MIT Libraries, on 22 May 2019 at 19:44:15, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2019.293

Bouncing phase variations in pilot-wave hydrodynamics

231

significantly over-predicts the experimentally measured walking speeds as γ /γFapproaches 1. By capturing the decrease in C with increasing γ in the (2, 1)2 mode,
our variable-phase model is able to better capture the experimentally observed plateau
in walking speeds at high memories. As shown in figure 12 in appendix A, our wave
model underpredicts the local wave amplitude, hp , for fast walkers at high memories,
compared to the model of Moláček & Bush (2013b) which is in agreement withexperimental wavefield measurements. This results in the slight decrease in thetheoretically predicted walking speed with increasing γ that is evident in figure 4(c)for a drop of radius R = 0.40 mm.
4.2. Impact-phase variations for multiple drops
When modelling droplet–droplet interactions, it is necessary to use a wave model that
accurately captures the far-field decay of the waves. Damiano et al. (2016), Tadrist
et al. (2018b) and Turton et al. (2018) all proposed the inclusion of a spatio-temporal
damping factor which modifies the wave kernel in (2.10) to an expression of the form
J0 (kF |x − xp (s)|) exp{−α|x − xp (s)|2 /(t − s)}. However, this factor is derived using along-time asymptotic expansion that introduces an unphysical singularity in the second
spatial derivative at |x − xp (s)| = 0. The presence of this singularity greatly complicatesa stability analysis of the trajectory equation (2.9). In order to analyse the stability oforbiting pairs, Oza et al. (2017) modified this damping factor by setting (t − s) →(t − s + TF ). In their study of promenading pairs, Arbelaiz et al. (2018) neglected thespatio-temporal damping factor altogether. To eliminate these mathematical difficulties,
in appendix B we use a quasi-static approximation to derive a purely spatial damping
factor that is smooth at |x − xp (s)| = 0, while remaining consistent with the spatiotemporal damping factor of Turton et al. (2018) in the far field.
Based on (2.9), the equations governing a system of N interacting drops areκ x̄¨ i + x̄˙ i = −βσi Ci ∇ h̄(x̄i , t̄),
where

N

AMe X
h̄(x̄, t̄) =
σj
R j=1

Z t̄

i = 1, 2, . . . , N,

(4.7)

Sj f (|x̄ − x̄j (s̄)|)e−(t̄−s̄) ds̄

(4.8)

−∞

and the wave kernel, with the spatial damping factor derived in appendix B, is−2

f (r̄) = J0 (r̄)[1 + (ξ K1 (ξ r̄)r̄ − 1)e−r̄ ],(4.9)

where ξ is defined in table 1. As a (2, 1) drop can either impact the bath during
the first or second cycle of the bath’s oscillation, the parameters σi = ±1 are chosen
to describe the relative bouncing phase of the drops: σi σj = 1 and σi σj = −1 indicatethat drops i and j are bouncing in-phase and out-of-phase, respectively. The phaseparameters Si and Ci describe modulations in the bouncing phase of the ith drop
within a given bath cycle.
The trajectory equation (4.7) highlights two main mechanisms through which
droplets interact. Firstly, the gradient of the strobed wavefield, ∇h, beneath a dropis altered due to waves generated by neighbouring drops. Secondly, neighbouringdrops also alter the local wave amplitude beneath a drop, hp , which will change thedrop’s vertical dynamics and hence S and C . This second interaction mechanism is
not captured by a constant phase model but can significantly affect droplet–droplet
interactions as we shall demonstrate in our analysis of the stability of droplet pairs.

Downloaded from https://www.cambridge.org/core. MIT Libraries, on 22 May 2019 at 19:44:15, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2019.293

M. M. P. Couchman, S. E. Turton and J. W. M. Bush232

4.3. Stability of bound droplet pairs
We analyse the stability of a bound pair of identical droplets using the system of(4.7). Pairs will exist in stationary bound states when the right-hand sides of (4.7)
are zero, corresponding to the discrete set of separation distances, dn , that satisfy
f 0 (dn ) = 0. In order for a separation distance dn to be stable, it must correspond toeach drop bouncing in a minimum of its neighbour’s wavefield (Eddi et al. 2009).
Thus, for drop’s bouncing in-phase (σ1 σ2 = 1) or out-of-phase (σ1 σ2 = −1), the stable
separation distances correspond to the minima or maxima of f (r), respectively. In
figure 5, we plot our theoretical predictions for stable values of dn at two values of γ ,the bouncing threshold, γB , and the walking threshold, γW , of an individual drop. The
values of dn are dependent on the vertical bouncing mode of the drops, as (1, 1) and(2, 1) drops excite waves of different wavelengths. Using the deep water dispersion
relation for gravity–capillary waves, ω2 = gk + σ k3 /ρ, we expect the wavelengths of
the wavefields produced by (1, 1) and (2, 1) bouncers to be λ(1,1) = 2.86 mm andλ(2,1) = λF = 4.75 mm, respectively. For the (1, 1) pairs in the vicinity of γB , ourtheoretical predictions for dn are found to be in good agreement with our experimental
data, except at the high binding numbers n = 10 and 11. This discrepancy presumably
arises because the drop–drop interaction is so weak at such low γ and large n thatthe drops had not fully settled into their equilibrium separation distances. For the(2, 1) pairs in the vicinity of γW , theoretical predictions are again found to be ingood agreement with the experimental data, but slightly underpredict the experimental
values of dn at small n. One expects this slight discrepancy to be caused by thetravelling wave fronts generated at each impact, an effect detailed by Galeano-Rios
et al. (2018) that is not captured by the wave model (4.8).
Next, we examine how stationary (2, 1) bound pairs destabilize into horizontal
motion as γ is increased. We assume that the drops start at stable positions(x1 , y1 ) = (0, 0) and (x2 , y2 ) = (dn , 0) and introduce arbitrary horizontal perturbations
x1 = δx1 , y1 = δy1 , x2 = dn + δx2 , y2 = δy2 . As shown in appendix C, substituting
these perturbations into (4.7) and retaining only terms up to O() yields the following
12 × 12 block diagonal, linear system:A3×3

dQ  0
=
0
dt
0

0
B 3×3

0
0

0
0
C 3×3

0


0
0 
Q,
0 

(4.10)

C 3×3

where
δx1 + δx2 , δu1 + δu2 , δX1 + δX2 , δx2 − δx1 , δu2 − δu1 , δX2 − δX1 ,δy1 , δv1 , δY1 , δy2 , δv2 , δY2 ,(4.11)




0
1
0
0
1
0


1
1



A = 
−χ − ψ − κ χ + ψ  , B = −χ − ψ − κ χ − ψ  ,
1
0
−1
1
0
−1


0
1
0
1


χ ,
C = −χ −(4.12)
κ
1
0 −1

Q =

Downloaded from https://www.cambridge.org/core. MIT Libraries, on 22 May 2019 at 19:44:15, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2019.293

Bouncing phase variations in pilot-wave hydrodynamics

233

AMe β Cn Sn
AMe β Cn Sn 00
f (0), ψ =σ1 σ2 f 00 (d̄n ).
(4.13a,b)Rκ
Rκ
The effect of the vertical dynamics on the pair’s horizontal stability is captured by
Sn and Cn , which are the values of the phase parameters for each drop in an initiallyhorizontally static pair with separation distance dn . Following the same procedure asin § 4.1, Sn and Cn can be found by solving the implicit equation (4.8) for the local
wave amplitude, h̄p :
χ=

h̄p =

AMe S (Γ , h̄p , Ω)
(1 + σ1 σ2 f (d̄n )).
R(4.14)

Having obtained h̄p for a given separation distance d̄n , the corresponding values of Sn
and Cn can be obtained using the relations in table 2.
In figure 6, we compare our theoretical predictions for Sn and Cn with ourexperimental data for pairs of radius R = 0.32 mm and 0.40 mm with binding
numbers n = 1 and 2. Both our theoretical model and experimental data highlightimportant features of the phase parameters close to the pair instability thresholds.
The pairs of radius R = 0.40 mm are in a (2, 1)1 bouncing mode and both Sn and
Cn increase with decreasing n. Conversely, the pairs of radius R = 0.32 mm are ina (2, 1)2 bouncing mode and while S has saturated to 1 for all n, C now decreases
with decreasing n. We proceed by demonstrating that this subtle difference in thedependences of S and C on n has a dramatic effect on a pair’s horizontal stability.
The block diagonal matrix in (4.10) indicates that there are three ways in whicha pair of initially stationary drops may destabilize into horizontal motion. The submatrix A governs the stability of the quantity δx1 + δx2 which corresponds to motionof the pair’s centre of mass in the x̂ -direction. All of the eigenvalues of A are realand one is zero, indicating that the pair is invariant to changes in the x̂ position of
its centre of mass. As γ is increased across the instability threshold γ∗ , the two nonzero eigenvalues switch from both being negative, to one being negative and the otherpositive. The resulting instability corresponds to a translation of the pair’s centre ofmass in the x̂-direction, where one drop follows the other while the inter-drop spacing,d, remains fixed. We refer to this instability as a co-linear translation.
The sub-matrix B governs the stability of the quantity δx2 − δx1 , which corresponds
to changes in the inter-drop distance, d. B has a real eigenvalue and a pair ofcomplex conjugate eigenvalues. As γ is increased across the instability threshold γ∗ ,the real eigenvalue remains negative, while the real parts of both complex conjugate
eigenvalues switch from negative to positive. The resulting instability corresponds toco-linear oscillations in the x̂-direction about a fixed centre of mass.
The sub-matrix C governs each drop’s stability to perturbations in the ŷ-direction,
which we will refer to as transverse perturbations. The onset of orbital andpromenading motion are both examples of such transverse instabilities. We notethat C is exactly the matrix obtained by Oza et al. (2013) that governs the bouncing
to walking transition of a single droplet. As the linear system (4.10) predicts thatthe transverse perturbations of each drop in a pair are independent of the motion of
the neighbouring droplet, it might seem that each drop will destabilize into walking
in the ŷ-direction at γW , the walking threshold of an individual droplet. However,this is not the case: although the pairs do not interact through the gradient of eachother’s wavefield, they modify each other’s local wave amplitude, hp , as captured by
the phase parameters S and C . Therefore, for transverse perturbations, the dropletsM. M. P. Couchman, S. E. Turton and J. W. M. Bush

©*/©W

(a)

R = 0.40 mm, (2, 1)1 mode
Variable-phase model

1.10
1.05
1.00
0.95
0.90

0

1

2

©*/©W

(b) 1.20
1.15
1.10
1.05
1.00
0.95
0.90

©*/©W

1.08
1.06
1.04
1.02
1.00
0.98

3

4

5

6

R = 0.40 mm, (2, 1)1 mode
Constant-phase model

0

1

2

(c) 1.10

3

4

5

6

R = 0.32 mm, (2, 1)2 mode
Variable-phase model

0

1

2

3

4

5

6

(d) 1.15©*/©W

Downloaded from https://www.cambridge.org/core. MIT Libraries, on 22 May 2019 at 19:44:15, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2019.293

234

R = 0.32 mm, (2, 1)2 mode
Constant-phase model

1.10
1.05
1.00
0.95

0

1

2

3

4

5

6

d/¬F

F IGURE 11. The theoretical instability thresholds, γ∗ , for co-linear translations of thepair’s centre of mass (green, matrix A), co-linear oscillations (blue, matrix B ) andtransverse motion (red, matrix C ) as a function of the pair’s initial separation distance,d. The corresponding experimental data from figure 7 are indicated by the black circles.
The panels illustrate predicted thresholds for (a) R = 0.40 mm with a variable-phase model,(b) R = 0.40 mm with a constant-phase model, (c) R = 0.32 mm with a variable-phasemodel and (d) R = 0.32 mm with a constant-phase model. The variable- and constant-phasemodels correspond to predictions based on the linear system (4.10) using the phasefunctions S and C listed in table 2, and the constant value of the phase parametersin(Φ)/2 = SC = 0.30 used by Oza et al. (2013), respectively.

interact solely through modulations in their vertical dynamics, an effect that cannot
be captured by a constant-impact-phase model.
In figure 11, we present the theoretical instability thresholds, γ∗ , for co-linear
translations of the pair’s centre of mass (matrix A), co-linear oscillations (matrix
B ) and transverse motions (matrix C ) as a function of the pair’s initial separation
distance, dn . We compare our experimental measurements from figure 7 to thetheoretical predictions based on both our variable-impact-phase functions S and C

Downloaded from https://www.cambridge.org/core. MIT Libraries, on 22 May 2019 at 19:44:15, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2019.293

Bouncing phase variations in pilot-wave hydrodynamics

235

listed in table 2, and the constant impact phase sin(Φ)/2 = SC = 0.30 used by Oza
et al. (2013). For pairs in the (2, 1)2 mode, our variable-phase model predicts thatthe pairs destabilize above the walking threshold of an individual drop, γW , and thatthe most unstable motion is transverse, in agreement with our experiments. We notethat for higher binding numbers n, the instability thresholds for transverse motion
and co-linear oscillations become progressively closer. This is consistent with ourexperimental observations that the promenading states begin to oscillate for n > 2.5.
For pairs in the (2, 1)1 mode, our variable-phase model predicts that the pairs
now destabilize below γW with the most unstable motion corresponding to co-linear
oscillations, again in agreement with our experimental results.
Comparing the predictions of our variable-phase model with those of a constantphase model makes clear that modulations in the vertical dynamics strongly affectthe stability of droplet pairs. For example, a constant-phase model cannot predict thedifference in behaviour between the (2, 1)1 and (2, 1)2 pairs: it predicts that both pairsizes should destabilize into co-linear oscillations below γW . The fact that (2, 1)1 and(2, 1)2 pairs destabilize below and above γW , respectively, is a direct consequence ofthe phase variations induced by the presence of a neighbouring drop. As the binding
number, n, of a pair decreases, so does the local wave amplitude beneath each drop,
hp , as each drop sits in a deeper minimum of its neighbour’s wavefield. Based ontable 2, we see that the phase parameters have opposite dependencies on hp for(2, 1)1 and (2, 1)2 drops. In the (2, 1)1 mode, both S and C increase with decreasing
hp . Therefore, pairs with smaller binding numbers n generate relatively large waves
and receive relatively large horizontal impulses from the bath, prompting an onset of
motion below γW . Conversely, for pairs in a (2, 1)2 mode, S has saturated to 1 while
C now decreases with decreasing hp causing the drops to stabilize each other above
γW .
5. Discussion

The trajectory equation (2.9), coupled with the impact phase functions S (Γ , h̄p , Ω)and C (Γ , h̄p , Ω) presented in table 2, provides a stroboscopic model for a drop’s
horizontal motion that accounts for modulations in the drop’s vertical dynamics. Ourmodel is valid for drops in the resonant (2, 1) bouncing mode, the regime of interestin the study of hydrodynamic quantum analogues. We have highlighted how variations
in a drop’s vertical motion affect both the amplitude of the wave generated at eachbounce, as described by S , as well as the horizontal impulse imparted on the dropby the bath, as described by C , quantities that strongly influence the drop’s horizontal
motion.
Previous models for the horizontal motion of walking droplets can be categorizedaccording to whether or not they account for variations in the vertical dynamics.
Models that neglect such variations (Oza et al. 2013; Bush et al. 2014; Dubertrand
et al. 2016; Durey & Milewski 2017; Faria 2017; Nachbin et al. 2017) are easier
to work with mathematically and have numerical simulation times comparable to,or faster than real time. However, they have limited predictive power in that afitting parameter resulting from the unresolved vertical dynamics must be used tomatch experimental data. Conversely, models that explicitly simulate a drop’s coupled
vertical and horizontal motion (Milewski et al. 2015; Galeano-Rios et al. 2017) haveno such fitting parameter but are computationally intensive, with numerical simulations
102 –104 times slower than real time. The new model developed in this study bridges
the gap between these two categories; specifically, it accounts for variations in a drop’sDownloaded from https://www.cambridge.org/core. MIT Libraries, on 22 May 2019 at 19:44:15, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2019.293

236

M. M. P. Couchman, S. E. Turton and J. W. M. Bushvertical dynamics while maintaining the simplicity of the stroboscopic models. Ourmodel is particularly suited to examining the role played by the vertical dynamics
in droplet chains (Filoux, Hubert & Vandewalle 2015) and lattices (Protière et al.
2005; Lieber et al. 2007; Eddi et al. 2008, 2009; Eddi, Boudaoud & Couder 2011a),systems that would require considerable amounts of computational power to perform
full numerical simulations of each drop’s coupled horizontal and vertical motion.
For single bouncers and steady, rectilinear walkers in a homogeneous bath, thephase parameters are uniquely determined by γ and R. As shown in figure 4, ourmodel was able to adequately predict the experimentally measured dependences of S
and C on γ for single bouncers and walkers of three different sizes across a rangeof driving accelerations spanning the (2, 1) mode. Based on these results, we can
now justify why the majority of droplets begin to walk only once they have enteredthe (2, 1)2 mode, as shown in the regime diagram in figure 1. In order for a drop
to walk, the wave amplitude must be sufficiently large to overcome the drag forces
resisting horizontal motion (Moláček & Bush 2013b). In the (2, 1)2 mode, S isclose to one, meaning that at each bounce the droplet produces a wave of maximum
possible amplitude, compared to a (2, 1)1 drop where S ≈ 0.5 and waves are only
generated with ≈50 % efficiency. Figure 4 also highlights the strong dependence ofthe phase parameters on the drop radius R. This explains why both Moláček & Bush(2013b) and Oza et al. (2013) could only successfully fit experimentally measured
walking thresholds and speeds for one drop size, using a constant-phase model.
By incorporating a variable phase, our model can now reproduce this experimental
data across a range of drop sizes. We note that Moláček & Bush (2013b) used avalue of sin (Φ)/2 = SC = 0.25 to best fit the experimental data for a drop of size
R = 0.40 mm. In figure 4, the product SC is seen to vary from approximately 0.35
to 0.2 for R = 0.40 mm in the walking regime, explaining why Moláček & Bush(2013b) were able to obtain an adequate fit using a value of 0.25. Furthermore,using a constant-phase model, both Moláček & Bush (2013b) and Oza et al. (2013)
overpredicted the walking speeds at high memories. By capturing the decrease in C
with increasing γ , our model now captures the experimentally observed plateau in
walking speeds at high memories, as is evident in figure 4.
For systems containing multiple droplets, the phase parameters can change
independently of γ and R, as shown in figure 6, due to changes in the local wave
amplitude induced by neighbouring drops. By characterizing how bound droplet pairs
destabilize as γ is increased, we find that variations in the pair’s vertical dynamics
strongly influence the nature of the droplet–droplet interaction. For example, figure 7
highlights that the interaction between droplet pairs in a (2, 1)1 bouncing mode
promotes an onset of instability below the walking threshold, γW , of an individual
drop. Conversely, the interaction between droplet pairs in a (2, 1)2 bouncing mode
has the opposite effect, stabilizing the pair above γW . This difference in behaviour
is due to variations in the pair’s impact phase. For (2, 1)1 bound pairs, both S andC increase with decreasing binding number n which augments the strength of thedroplet–droplet interaction and promotes instability. Conversely, for (2, 1)2 boundpairs, S ≈ 1 remains roughly constant while C decreases with decreasing n, resulting
in stabilization. As shown in figure 11, our variable-impact-phase model is able tocapture both the type and thresholds of instability observed in our experiments, whilea constant-impact-phase model is not. This distinction highlights the importance ofaccounting for phase variations when modelling systems containing multiple droplets.
Our study also allows us to better understand the results of the studies of Oza et al.
(2017) and Arbelaiz et al. (2018) on orbiting and promenading droplet pairs. In both

Downloaded from https://www.cambridge.org/core. MIT Libraries, on 22 May 2019 at 19:44:15, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2019.293

Bouncing phase variations in pilot-wave hydrodynamics

237

studies, the walking speed of the pairs was observed to decrease below the free speedof an individual drop with decreasing binding number n, prompting both studies todevelop ad hoc empirical functions for sin(Φ)/2 = SC that decreased with decreasingn in order to fit their experimental data. Our study now rationalizes these trends interms of the droplet pair’s variable vertical dynamics. For (2, 1)2 drops, the local wave
amplitude and hence the phase parameter C decrease with decreasing n, reducing thehorizontal impulse the pair receives from the bath. We note that for drops in a (2, 1)1
mode, the opposite behaviour is expected on the basis of our study: the walking speedof the pair should increase with decreasing n, thus exceeding the walking speed ofa single drop. Accounting for a droplet’s vertical dynamics is also expected to beimportant when modelling a drop’s interaction with a submerged barrier. For example,
Pucci et al. (2016) observed experimentally that a walking droplet will slow down inthe vicinity of a submerged barrier and Damiano et al. (2016) observed experimentallythat the wavefield gets significantly damped in the vicinity of such a barrier. The
results of our study could help link these two observations: for a (2, 1)2 walker, ourmodel predicts that C and hence a drop’s walking speed will decrease when the local
wave amplitude decreases, as one expects to arise in the vicinity of barriers.
Although our derivations of the phase functions S and C were restricted to drops
in a (2, 1) bouncing mode, our analysis of the vertical dynamics can be readilyextended to other modes. As shown in figure 4, our models of S and C for (2, 1)
bouncers already seem to approximate well the experimentally observed impact
phases in the neighbouring (2, 2) and (4, 2) regimes. Moreover, it would be relatively
straightforward to develop phase functions for (1, 1) bouncers, for instance. First,
the approximation hp = 0 could be used in (4.1) as the surface wave amplitude isnegligible compared to the vibrational amplitude of the bath, γ /ω2 , for the smallvalues of γ for which (1, 1) bouncers exist. Then, a similar analysis to that presented
in §4 could be performed, this time seeking periodic solutions to the linear springmodel (4.2) with period TF /2. We also note that Tadrist et al. (2018a) have recently
found that, in certain parameter regimes, two interacting droplets can suddenly
flip their relative bouncing phases. In our model, the relative bouncing phases areprescribed through the parameters σi = ±1 in (4.7), and so cannot change in time.
However, if a mechanism causing such reversals in bouncing phase were to be
identified, it could be incorporated into (4.7).
When modelling multiple droplets, it may also be necessary to use a more complex
wave model than that presented in (4.8). For instance, Galeano-Rios et al. (2018)
found that travelling wave fronts have a significant influence on the horizontal motionof ratcheting pairs (Eddi et al. 2008). We note that the phase functions in table 2 arevalid for any wave model in which the drop produces a standing wave at each impactthat oscillates with half the frequency of the bath. To develop phase functions for a
traveling wave model, the more complicated time dependence of the surface waves
would first have to be incorporated into (4.1), after which the same analysis presented
in this study could be followed.
Finally, a subject of current interest is the droplet–droplet correlations that
may be induced by wave-mediated forces, such as drops communicating between
one-dimensional slots (Nachbin 2018) and in hydrodynamic spin lattices (Sáenz et al.
2018). As these correlations are expected to depend on the drop’s vertical dynamics,
our model is likely to yield valuable insight into the interactions and statisticalcorrelations emerging in these subtle multiple-drop systems.

Downloaded from https://www.cambridge.org/core. MIT Libraries, on 22 May 2019 at 19:44:15, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2019.293

238

M. M. P. Couchman, S. E. Turton and J. W. M. BushAcknowledgements

M.M.P.C. gratefully acknowledges the financial support of the Natural Sciencesand Engineering Research Council (NSERC) through grant no. 502891. J.W.M.B.
gratefully acknowledges the financial support of the National Science Foundation(NSF) through grants DMS-1614043 and CMMI-1727565. We thank Professor R. R.
Rosales for valuable discussions, and particularly his contribution to the theoreticaldevelopments in appendix B.
Appendix A. Wave-amplitude correction

We here compare the local wave amplitudes, hp , predicted by the discrete wave
model of Moláček & Bush (2013b) (2.1) and the continuous wave model of Oza et al.
(2013) (2.8) for single bouncers and walkers. For a droplet walking with speed u ina straight line, the model of Moláček & Bush (2013b) predicts thatr∞X
2νe mgkF3 R2 TF
J0 (kF unTF ) −n/Me
√
hp =
e
,
(A 1)
S
2 2
π σ (3kF R + Bo) n=1
nTF
where the sum must be performed numerically. The model of Oza et al. (2013)
approximates the sum in the model of Moláček & Bush (2013b) by an integral which
results in a different local wave amplitude for a droplet walking with speed u:
s
Z ∞
2νe
mgkF3 R2hp =S
J0 (kF uz)e−z/(TF Me ) dz
(A 2)
πTF σ (3kF2 R2 + Bo)0s
mgkF3 R2
2νe
TF Me
√
S.
(A 3)
=
2 2
πTF σ (3kF R + Bo) 1 + (kF TF Me u)2
Figure 12 shows the dependence of hp on γ as predicted by the models of Moláček& Bush (2013b) (A 1) and Oza et al. (2013) (A 3) for bouncers (u = 0) and walkers
with prescribed speeds of u = 5 and 10 mm s−1 . In all cases, the wave amplitudeof Oza et al. (2013) is seen to overpredict that of Moláček & Bush (2013b). This√
discrepancy is due to the approximation of Oza et al. (2013)√ to replace the 1/ t
decay factor in the model of Moláček & Bush (2013b) by 1/ TF , in order to avoid
the singularity that would otherwise appear in the integral model at t = 0.
For a bouncer, one can calculate the factor needed to correct this discrepancy.
Specifically, the local wave amplitudes predicted by the model of Moláček & Bush(2013b) (A 1) and Oza et al. (2013) (A 3) differ by the factor Li1/2 (e−1/Me )/Me ,where Lin (x) denotes the polylogarithm function of order n. However, for a walker,
a different correction factor would be required that would depend not only on thememory, Me , but also on the walker’s speed, u, and past trajectory. For example, thescaling factor would be different for a drop walking in a circle compared to a drop
walking in a straight line. To more accurately approximate the discrete wave modelof Moláček & Bush (2013b) using an integral would require the development of anon-singular model for the behaviour of the wavefield for times t < TF , immediately
following a drop’s impact. This is beyond the scope of this paper which is focusedon modelling the coupling between a drop’s vertical and horizontal dynamics. Instead,in figure 12 we observe that for both bouncers and walkers of various speeds acrossa range of memories, the integral model of Oza et al. (2013) overpredicts hp by

(a) 0.6(b) 0.6(c) 0.6

0.4

0.4

0.4

hp/R

Downloaded from https://www.cambridge.org/core. MIT Libraries, on 22 May 2019 at 19:44:15, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2019.293

Bouncing phase variations in pilot-wave hydrodynamics

0.2
0
0.5

0.2
-1

u = 0 mm s
0.6

0.7

©/©F

0.8

0.9

0
0.7

239

0.2
-1

u = 5 mm s
0.8

©/©F

0.9

1.0

0
0.7

u = 10 mm s-1
0.8

©/©F

0.9

1.0

F IGURE 12. The non-dimensional local wave amplitude, hp /R, as a function of thenon-dimensional driving acceleration, γ /γF , for a drop of radius R = 0.36 mm that is (a)bouncing, (b) walking at speed u = 5 mm s−1 and (c) walking at speed u = 10 mm s−1 .
The black line corresponds to the discrete model of Moláček & Bush (2013b), the blueline to the continuous model of Oza et al. (2013) and the red line to the model ofOza et al. (2013) with the wave amplitude divided by 2, as used in this study to better
approximate the model of Moláček & Bush (2013b); S = 1 has been used for the sakeof comparison.

approximately a factor of two relative to the predictions of the discrete model ofMoláček & Bush (2013b). Thus, in order to obtain a qualitative agreement between
the wave amplitudes of Oza et al. (2013) and Moláček & Bush (2013b), we proceedby dividing the wave amplitude of Oza et al. (2013) by a factor of two.
Appendix B. Quasi-static spatial damping factor

We use a quasi-static approximation to derive a purely spatial damping factor thatremains consistent with the spatio-temporal damping factor proposed by Turton et al.
(2018) but that is smooth at the origin. Including the spatio-temporal damping factor
of Turton et al. (2018) in the wave model (2.10) yields
Z tAh(x, t) =S f (|x − xp (s)|, t − s)e−(t−s)/(TF Me ) ds,
(B 1)
TF −∞
where

f (r, t) = J0 (kF r)e−αr /t ,
2

(B 2)

and α is defined in table 1.
Consider the wavefield produced by a drop bouncing in place at the position xp = 0:Z ∞AS
2
J0 (kF |x|)e−α|x| /z−z/(TF Me ) dz
(B 3)
hB (|x|) =
TF
0
r
 r

αMe
α
= 2AS
J0 (kF |x|)K1 2|x| |x|.
(B 4)
TF
TF Me
√
We note that the integrand in (B 3) is sharply peaked around z = |x| αTF Me meaningthat the bouncer√ wavefield, hB , at a distance |x| from the drop forms on a time scaleof τF (|x|) = |x| αTF Me . To get a sense of this time scale, τF ≈ 3TF at the intermediate
distance |x| = 2λF for γ /γF = 0.85. In our experiments, a typical horizontal speed ofa droplet in a pair is |ẋ| ≈ 5 mm s−1 , as shown in figure 8. At this speed, a drop will
only move horizontally by approximately one drop radius, R, over the time scale, τF ,

Downloaded from https://www.cambridge.org/core. MIT Libraries, on 22 May 2019 at 19:44:15, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2019.293

240

M. M. P. Couchman, S. E. Turton and J. W. M. Bushtaken for a bouncer wavefield to form. We can thus use a quasi-static approximation
to simplify the spatio-temporal damping factor of Turton et al. (2018). Specifically,
we can assume that at each instant a walking droplet locally generates a bouncer
wavefield, allowing the spatio-temporal damping kernel (B 2) to be approximated bya simpler, purely spatial damping kernel
 r

r
α
α
rJ0 (kF r)K1 2
r .
(B 5)
f (r) = 2
TF Me
TF Me
As with the original spatio-temporal damping factor (B 2), the purely spatial
damping factor (B 5) still has a singularity at r = 0. This singularity arises because
equation (B 2) was derived using a far-field asymptotic expansion invalid for small r,
where we expect f (r) = J0 (kF r) which is not singular at r = 0. Therefore, we wish tosmoothly connect the near field J0 (kF r) and far-field (B 5) profile of the wave, where
the transition between the two profiles should occur at a distance r ∼ λF from thedrop. We thus adopt the following wave kernel

 r
 r



α
α
−1/(kF r)2
f (r) = J0 (kF r) 1 + 2
K1 2
r r−1 e
,
(B 6)
TF Me
TF Me
where the function exp[−1/(kF r)2 ] smoothly transitions from 0 to 1 at a characteristic
distance r = 1/kF .
A comparison between the wavefields predicted by the spatio-temporal damping
kernel of Turton et al. (2018) (B 2) and the purely spatial damping kernel adopted
in this study (B 6) is shown in figure 13 for walkers travelling at various speeds.
For a bouncer, the wavefields predicted using the wave kernels (B 2) and (B 6)
are indistinguishable. For a walker, the agreement between the wavefields remains
excellent for the typical walking speeds in our experiments (u . 5 mm s−1 ). For
larger u, the agreement becomes progressively worse, especially at large x, owing to
the limitations of the quasi-static approximation.
Appendix C. Linear stability analysis for droplet pairs

We here derive the linear system (4.10) that governs the stability of a pair ofidentical stationary bouncers separated by a distance dn . For notational simplicity, we
omit the overbars denoting non-dimensional variables. Substituting the perturbations
x1 = δx1 , y1 = δy1 , x2 = dn + δx2 and y2 = δy2 into the system (4.7), using the factthat f 0 (0) = f 0 (dn ) = 0, and only keeping terms up to O() givesZ
AMe β Cn Sn t  00
f (0)(δx1 (t) − δx1 (s))κδ ẍ1 (t) + δ ẋ1 (t) = −R−∞+ σ1 σ2 f 00 (dn )(δx1 (t) − δx2 (s)) · x̂ e−(t−s) ds,Z
AMe β Cn Sn t  00
f (0)(δx2 (t) − δx2 (s))κδ ẍ2 (t) + δ ẋ2 (t) = −R−∞+ σ1 σ2 f 00 (dn )(δx2 (t) − δx1 (s)) · x̂ e−(t−s) ds,

(C 1)

(C 2)

where x̂ = (1, 0) and Sn and Cn denote the values of the phase parameters for drops
separated by a distance dn in a horizontally stationary pair.

(a)

h/R

0.2

0.1

0.1

0

0

-0.1
-5

0

241

(b)

u = 0 mm s-1

0.2

5

(c)

-0.1
-5

u = 5 mm s-1

0

5

(d)
0.2

0.2
u = 10 mm s-1

h/R

Downloaded from https://www.cambridge.org/core. MIT Libraries, on 22 May 2019 at 19:44:15, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2019.293

Bouncing phase variations in pilot-wave hydrodynamics

0.1

0.1

0

0

-0.1
-5

0

5

-0.1
-5

x/¬F

u = 20 mm s-1

0

5

x/¬F

F IGURE 13. A comparison of the wavefields predicted by (B 1) using the spatio-temporal
damping kernel (B 2) of Turton et al. (2018) (black) and the spatial damping wave kernel(B 6) used in this study (red). The comparison is made for a drop of size R = 0.40 mm
at γ /γF = 0.85 with S = 1 for (a) a bouncer, and walkers of speed (b) u = 5 mm s−1 ,(c) u = 10 mm s−1 , and (d) u = 20 mm s−1 . The droplet is at the position x = 0 in each
plot.

Rt
Following Oza et al. (2013), we then apply the substitution δX(t)= −∞ δx(s)e−(t−s) ds
to (C 1) and (C 2) which yields

AMe β Cn Sn  00
f (0)(δx1 − δX1 ) + σ1 σ2 f 00 (dn )(δx1 − δX2 ) · x̂ , (C 3)R
AMe β Cn Sn  00
κδ ẍ2 + δ ẋ2 = −f (0) (δx2 − δX2 ) + σ1 σ2 f 00 (dn ) (δx2 − δX1 ) · x̂ , (C 4)R
δ Ẋ1 = δx1 − δX1 ,(C 5)
δ Ẋ2 = δx2 − δX2 .
(C 6)
κδ ẍ1 + δ ẋ1 = −Finally, by making the substitution δui = δ ẋi , the linear system (4.10) is obtained.
REFERENCES

A RBELAIZ , J., O ZA , A. U. & B USH , J. W. M. 2018 Promenading pairs of walking droplets: dynamics
and stability. Phys. Rev. Fluids 3 (1), 013604.
B ENJAMIN , T. B. & U RSELL , F. 1954 The stability of the plane free surface of a liquid in vertical
periodic motion. Proc. R. Soc. Lond. A 225, 505–515.

Downloaded from https://www.cambridge.org/core. MIT Libraries, on 22 May 2019 at 19:44:15, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2019.293

242

M. M. P. Couchman, S. E. Turton and J. W. M. BushB ORGHESI , C., M OUKHTAR , J., L ABOUSSE , M., E DDI , A., F ORT, E. & C OUDER , Y. 2014 Interactionof two walkers: wave-mediated energy and force. Phys. Rev. E 90 (6), 063017.
DE B ROGLIE , L. 1956 Une Tentative d’Interprétation Causale et non Linéaire de la Mécanique
Ondulatoire: la théorie de la double solution. Gauthier-Villars.
B USH , J. W. M. 2015a The new wave of pilot-wave theory. Phys. Today 68 (8), 47–53.
B USH , J. W. M. 2015b Pilot-wave hydrodynamics. Annu. Rev. Fluid Mech. 47, 269–292.
B USH , J. W. M., C OUDER , Y., G ILET, T., M ILEWSKI , P. A. & NACHBIN , A. 2018 Introduction to
focus issue on hydrodynamic quantum analogs. Chaos 28 (9), 096001.
B USH , J. W. M., O ZA , A. U. & M OLÁČEK , J. 2014 The wave-induced added mass of walking
droplets. J. Fluid Mech. 755, R7.
C OUDER , Y., F ORT, E., G AUTIER , C.-H. & B OUDAOUD , A. 2005a From bouncing to floating:noncoalescence of drops on a fluid bath. Phys. Rev. Lett. 94 (17), 177801.
C OUDER , Y., P ROTIERE , S., F ORT, E. & B OUDAOUD , A. 2005b Dynamical phenomena: walking andorbiting droplets. Nature 437 (7056), 208.
DAMIANO , A. P., B RUN , P-T., H ARRIS , D. M., G ALEANO -R IOS , C. A. & B USH , J. W. M. 2016
Surface topography measurements of the bouncing droplet experiment. Exp. Fluids 57 (10),163.
D UBERTRAND , R., H UBERT, M., S CHLAGHECK , P., VANDEWALLE , N., BASTIN , T. & M ARTIN , J.
2016 Scattering theory of walking droplets in the presence of obstacles. New J. Phys. 18(11), 113037.
D UREY, M. & M ILEWSKI , P. A. 2017 Faraday wave–droplet dynamics: discrete-time analysis. J. Fluid
Mech. 821, 296–329.
E DDI , A., B OUDAOUD , A. & C OUDER , Y. 2011a Oscillating instability in bouncing droplet crystals.
Europhys. Lett. 94 (2), 20004.
E DDI , A., D ECELLE , A., F ORT, E. & C OUDER , Y. 2009 Archimedean lattices in the bound states ofwave interacting particles. Europhys. Lett, 87 (5), 56002.
E DDI , A., S ULTAN , E., M OUKHTAR , J., F ORT, E., ROSSI , M. & C OUDER , Y. 2011b Information
stored in Faraday waves: the origin of a path memory. J. Fluid Mech. 674, 433–463.
E DDI , A., T ERWAGNE , D., F ORT, E. & C OUDER , Y. 2008 Wave propelled ratchets and drifting rafts.
Europhys. Lett. 82 (4), 44001.
FARIA , L. M. 2017 A model for Faraday pilot waves over variable topography. J. Fluid Mech. 811,
51–66.
F ILOUX , B., H UBERT, M. & VANDEWALLE , N. 2015 Strings of droplets propelled by coherent waves.
Phys. Rev. E 92 (4), 041004.
G ALEANO -R IOS , C. A., C OUCHMAN , M. M. P., C ALDAIROU , P. & B USH , J. W. M. 2018 Ratcheting
droplet pairs. Chaos 28 (9), 096112.
G ALEANO -R IOS , C. A., M ILEWSKI , P. A. & VANDEN -B ROECK , J.-M. 2017 Non-wetting impact ofa sphere onto a bath and its application to bouncing droplets. J. Fluid Mech. 826, 97–127.
G ILET, T., VANDEWALLE , N. & D ORBOLO , S. 2009 Completely inelastic ball. Phys. Rev. E 79 (5),055201.
H AMMOND , R. T. 2010 Relativistic particle motion and radiation reaction in electrodynamics. Electron.
J. Theoret. Phys. 7 (23), 221–258.
H ARRIS , D. M. & B USH , J. W. M. 2015 Generating uniaxial vibration with an electrodynamic
shaker and external air bearing. J. Sound Vib. 334, 255–269.
H ARRIS , D. M., L IU , T. & B USH , J. W. M. 2015 A low-cost, precise piezoelectric droplet-on-demand
generator. Exp. Fluids 56 (4), 83.
L ABOUSSE , M., O ZA , A. U., P ERRARD , S. & B USH , J. W. M. 2016 Pilot-wave dynamics in aharmonic potential: quantization and stability of circular orbits. Phys. Rev. E 93 (3), 033122.
L IEBER , S. I., H ENDERSHOTT, M. C., PATTANAPORKRATANA , A. & M ACLENNAN , J. E. 2007
Self-organization of bouncing oil drops: two-dimensional lattices and spinning clusters. Phys.
Rev. E 75 (5), 056308.
M ILES , J. & H ENDERSON , D. 1990 Parametrically forced surface waves. Annu. Rev. Fluid Mech. 22(1), 143–165.

Downloaded from https://www.cambridge.org/core. MIT Libraries, on 22 May 2019 at 19:44:15, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2019.293

Bouncing phase variations in pilot-wave hydrodynamics

243

M ILEWSKI , P. A., G ALEANO -R IOS , C. A., NACHBIN , A. & B USH , J. W. M. 2015 Faraday pilot-wave
dynamics: modelling and computation. J. Fluid Mech. 778, 361–388.
M OLÁČEK , J. & B USH , J. W. M. 2013a Drops bouncing on a vibrating bath. J. Fluid Mech. 727,
582–611.
M OLÁČEK , J. & B USH , J. W. M. 2013b Drops walking on a vibrating bath: towards a hydrodynamic
pilot-wave theory. J. Fluid Mech. 727, 612–647.
NACHBIN , A. 2018 Walking droplets correlated at a distance. Chaos 28 (9), 096110.
NACHBIN , A., M ILEWSKI , P. A. & B USH , J. W. M. 2017 Tunneling with a hydrodynamic pilot-wave
model. Phys. Rev. Fluids 2 (3), 034801.
O ZA , A. U., H ARRIS , D. M., ROSALES , R. R. & B USH , J. W. M. 2014a Pilot-wave dynamics ina rotating frame: on the emergence of orbital quantization. J. Fluid Mech. 744, 404–429.
O ZA , A. U., ROSALES , R. R. & B USH , J. W. M. 2013 A trajectory equation for walking droplets:
hydrodynamic pilot-wave theory. J. Fluid Mech. 737, 552–570.
O ZA , A. U., S IÉFERT, E., H ARRIS , D. M., M OLÁČEK , J. & B USH , J. W. M. 2017 Orbiting pairsof walking droplets: dynamics and stability. Phys. Rev. Fluids 2, 053601.
O ZA , A. U., W IND -W ILLASSEN , Ø., H ARRIS , D. M., ROSALES , R. R. & B USH , J. W. M. 2014b
Pilot-wave hydrodynamics in a rotating frame: exotic orbits. Phys. Fluids 26 (8), 082101.
P ERRARD , S. & L ABOUSSE , M. 2018 Transition to chaos in wave memory dynamics in a harmonic
well: deterministic and noise-driven behavior. Chaos 28 (9), 096109.
P ROTIÈRE , S., B OHN , S. & C OUDER , Y. 2008 Exotic orbits of two interacting wave sources. Phys.
Rev. E 78 (3), 036204.
P ROTIÈRE , S., B OUDAOUD , A. & C OUDER , Y. 2006 Particle–wave association on a fluid interface.
J. Fluid Mech. 554, 85–108.
P ROTIÈRE , S., C OUDER , Y., F ORT, E. & B OUDAOUD , A. 2005 The self-organization of capillary
wave sources. J. Phys.: Condens. Matter 17 (45), S3529–S3535.
P UCCI , G., S ÁENZ , P. J., FARIA , L. M. & B USH , J. W. M. 2016 Non-specular reflection of walking
droplets. J. Fluid Mech. 804, R3.
S ÁENZ , P. J., P UCCI , G., G OUJON , A., C RISTEA -P LATON , T., D UNKEL , J. & B USH , J. W. M.
2018 Spin lattices of walking droplets. Phys. Rev. Fluids 3 (10), 100508.
TADRIST, L., S AMPARA , N., S CHLAGHECK , P. & G ILET, T. 2018a Interaction of two walkers:
perturbed vertical dynamics as a source of chaos. Chaos 28 (9), 096113.
TADRIST, L., S HIM , J.-B., G ILET, T. & S CHLAGHECK , P. 2018b Faraday instability and subthreshold
Faraday waves: surface waves emitted by walkers. J. Fluid Mech. 848, 906–945.
T ERWAGNE , D., L UDEWIG , F., VANDEWALLE , N. & D ORBOLO , S. 2013 The role of the droplet
deformations in the bouncing droplet dynamics. Phys. Fluids 25 (12), 122101.
T URTON , S. E., C OUCHMAN , M. M. P. & B USH , J. W. M. 2018 A review of the theoreticalmodeling of walking droplets: towards a generalized pilot-wave framework. Chaos 28 (9),096111.
WALKER , J. 1978 The amateur scientist. Drops of liquid can be made to float on liquid. Whatenables them to do so? Sci. Am. 238 (6), 151.
W IND -W ILLASSEN , Ø., M OLÁČEK , J., H ARRIS , D. M. & B USH , J. W. M. 2013 Exotic states ofbouncing and walking droplets. Phys. Fluids 25 (8), 082002.